© The Authors 2023

1. Introduction

Forecasting solar energetic particle (SEP) events is a highly intricate task that presents one of the most complex challenges in space weather modelling. Such events involve a range of associated phenomena with vastly different lengths and timescales, while the accelerated particles themselves can traverse the distance from the Sun to Earth in several minutes. Two different types of SEPs exist, proton-rich and electron-rich SEPs (Reames 1999). The focus of this study is on the proton-rich SEPs, which are generally considered to be accelerated at shock waves driven by coronal mass ejections (CMEs) and are often referred to as gradual SEP events (see e.g. Desai & Giacalone 2016, and references therein). The magnetic connection of an observer to the CME-driven shock wave is essential in shaping the measured particle intensity profiles (e.g. Cane et al. 1988; Reames et al. 1997). Under nominal solar wind conditions, when a CME occurs an observer situated to the east of the CME will establish the most favourable magnetic connection at the event onset, owing to the spiral shape of the interplanetary magnetic field (IMF). However, deviations from the nominal solar wind conditions, such as the presence of preceding interplanetary CMEs, often complicate the precise determination of an observer’s magnetic connection to the shock wave (e.g. Palmerio et al. 2022). Furthermore, recent studies by Ding et al. (2022a) and Wijsen et al. (2023) have demonstrated the significant influence of a non-uniform solar wind on the properties of CME-driven shocks, which subsequently impacts the observed intensity profiles of SEPs. These publications demonstrate the crucial importance of accurately modelling the solar wind in any physics-based SEP forecasting tool.

Current global 3D magnetohydrodynamic (MHD) models of the inner heliosphere are mostly capable of reproducing the large-scale solar wind phenomena, but they typically lack the resolution to reproduce the small-scale structures. Enlil (Odstrcil 2003) was among the first 3D MHD heliospheric models used for operational as well as research purposes. The EUropean Heliospheric FORecasitng Information Asset (EUHFORIA; Pomoell & Poedts 2018) is another state-of-the-art MHD forecasting-targeted model that couples a data-driven coronal model with a physics-based time-dependent 3D heliospheric model.

Coupling full heliospheric models, such as EUHFORIA, with SEP models to obtain more accurate prediction tools has been a challenge. The SEPMOD model (Luhmann et al. 2007, 2010, 2017) propagates energetic protons scatter-free along magnetic field lines extracted from an Enlil simulation. In this model, the SEPs are accelerated by the CME-driven shock wave. The SEP emission function intensity depends on the shock’s velocity jump via the semi-empirical relationship presented in Lario et al. (1998) and its energy spectrum on the compression ratio (Luhmann et al. 2010). This model is currently available to users through the Community Coordinated Modelling Center¹ (CCMC). In a study conducted by Kozarev et al. (2010), the Enlil model was coupled with the Energetic Particle Radiation Environment Module (EPREM; Schwadron et al. 2010) to investigate the SEP transport during the complex Halloween 2003 SEP events. The MHD simulations captured the fast and slow solar wind streams observed, but CME simulations were not included at that time. These series of events were observed by both near-Earth spacecraft and the Ulysses spacecraft, which was positioned at a distance of approximately 5 au from the Sun. EPREM models the transport of SEPs by solving the focussed transport equation (FTE; Skilling 1971; Isenberg 1997; le Roux & Webb 2009, 2012; Zank et al. 2014) using a finite difference approach on a Lagrangian mesh. Similarly, the EUHFORIA model has been coupled with PArticle Radiation Asset Directed at Interplanetary Space Exploration (PARADISE; Wijsen et al. 2019a; Wijsen 2020), which solves the FTE in a stochastic manner by integrating the equivalent set of Itô stochastic differential equations forward in time. The novelty of the PARADISE model is that it solves the full FTE while using the time-dependent EUHFORIA simulations as the background solar wind. In addition, the model can include pitch-angle-dependent cross-field diffusion and the particles’ guiding centre drifts (van den Berg et al. 2021). The PARADISE model has recently been employed to simulate gradual SEP events (Wijsen et al. 2022, 2023). These studies show the significant impact of the high-speed streams and the stream interaction regions on the observed SEP intensities (Wijsen et al. 2022, 2023). In this work, we build on this by using EUHFORIA and PARADISE to investigate the gradual SEP event that occurred on 2013 March 15, which took place within a highly perturbed solar wind environment due to the presence of multiple preceding CMEs.

When simulating CMEs with global MHD codes such as Enlil or EUHFORIA, several different CME models are available that vary in their treatment of the internal magnetic field of the CME. The simplest model, referred to as the cone model (e.g. Xie et al. 2004; Millward et al. 2013), is often employed with both Enlil (e.g. Soni et al. 2023; Odstrcil et al. 2020) and EUHFORIA (e.g. Scolini et al. 2019; Verbeke et al. 2022b). This model considers CMEs as a hydrodynamic cloud characterised by elevated densities, temperatures, and speeds. However, due to its simplicity and the absence of an internal magnetic field, it fails to capture important phenomena such as the expansion of the CME resulting from magnetic pressure. Moreover, the lack of a complex internal magnetic field in the cone model implies that it is permeated by the IMF, with minimal distortion. Consequently, the cone model is not very appropriate for accurately modelling the transport of SEPs downstream of the CMEs’ shock wave. This inaccuracy arises because SEPs, being charged particles, gyrate along the magnetic field lines and are thus strongly influenced by the underlying magnetic field structure.

To some extent, the limitations of the cone model can be mitigated by utilising more sophisticated magnetised CME models, such as FRi3D (Isavnin 2016; Maharana et al. 2022) or the linear force-free spheromak (LFF-spheromak) model (Verbeke et al. 2019; both these models recently became available in EUHFORIA). These advanced CME models offer a closer approximation and more realistic description of the complex CME structures we see in observations (Scolini et al. 2019; Asvestari et al. 2021). The difficulty that arises with these complex CME models is the necessity to incorporate several parameters that describe the CMEs’ internal magnetic field. Namely, accurate estimation of these parameters from remote observations is exceedingly challenging. When the objective is not real-time forecasting, the model parameters of studied CMEs can be further constrained using in situ observations. However, it is important to note that these observations are limited to a few specific points in the heliosphere (i.e. at the spacecraft positions), making it challenging to assess the model’s performance in other regions of the heliosphere.

Solar energetic particle observations offer a valuable opportunity to explore solar wind regions that are not directly observed by spacecraft. Solar energetic particles propagate at a significantly higher speed than the underlying solar wind and closely follow trajectories along magnetic field lines, which is somewhat different to the predominantly radial solar wind flows. As a result, SEP observations can provide insights into extended regions of the heliosphere that are not covered by direct solar wind and IMF spacecraft observations. Furthermore, by assessing the performance of the SEP models, we can gain a more comprehensive understanding of the overall underlying solar wind and CME models.

In this study we combine the EUHFORIA and PARADISE models to simulate the gradual SEP event that took place on 2013 March 15. This event was closely preceded by several CMEs, which resulted in strongly perturbed solar wind conditions. To investigate the impact of these preceding CMEs, we conducted a series of EUHFORIA simulations using different CME models. Subsequently, we employed the PARADISE model to propagate SEPs through these varying solar wind configurations. By comparing the simulation results with the observations, we gain valuable insights into the performance of both EUHFORIA and PARADISE on a larger part of the computational domain (as opposed to the more standard point-like performance analysis). Our findings highlight the significant influence of preceding CMEs on the observed SEP intensities. Moreover, we demonstrate that the choice of the CME models employed in the simulations plays a crucial role in shaping the behaviour of the modelled SEPs. Our findings emphasise the ongoing need for advancements in both CME and SEP modelling techniques in order to enhance our understanding and predictive capabilities of SEP events.

The rest of this article is organised as follows. Section 2 presents the relevant observations and a description of the models and methodology employed for the analysis. In Sect. 3 a detailed description of the remote sensing observations and in situ measurements is provided. Section 4 focusses on the discussion of solar wind modelling. To ensure an accurate EUHFORIA background solar wind, a CME reconstruction is performed using the graduated cylindrical shell (GCS; Thernisien et al. 2009) model, which is elaborated upon in this section. The results of the PARADISE simulation and relevant discussions are presented in Sect. 5. Lastly, Sect. 6 summarises the results and provides our concluding remarks.

2. Data, models, and methods

In this study we used observations at a number of different wavelengths and from different instruments. in situ measurements are taken from spacecraft positioned at the Lagrangian point L1. The 3D MHD model EUHFORIA was used for modelling solar wind and the propagation of CMEs previously reconstructed with the GCS model. To study the particle propagation we used the recently developed model PARADISE.

2.1. Data

We used data from X-Ray Sensor (XRS; Hanser & Sellers 1996) aboard Geostationary Operational Environmental Satellite (GOES) to determine the peak intensity of the flare associated with the event. We also used extreme ultraviolet (EUV) observations from Atmospheric Imaging Assembly (AIA) aboard the Solar Dynamics Observatory (SDO; Lemen et al. 2012). The EUV observations were employed to locate the source region of the flare, the possible existence of the on-disk signatures of the CMEs (waves and coronal dimmings) and the prominence associated with the eruption. In particular, flare source regions and post-eruptive signatures were detected using the AIA 94 Å channel, while coronal hole extension and other post-eruptive CME signatures were analysed using AIA 193 Å and 304 Å channels.

For studying the propagation of the CME and making the reconstructions necessary for modelling the CMEs with EUHFORIA, coronagraph data were extensively analysed. Coronagraph images were taken from COR2-A and COR2-B, aboard the Solar TErrestrial RElations Observatory (STEREO; Howard et al. 2008) and by Large Angle Spectroscopic COronagraph (LASCO) aboard the SOlar and Heliospheric Observatory (SOHO; Brueckner et al. 1995). All three coronagraphs had available images within a 5-day window to the event of interest, with minimal or no data gaps. STEREO-A and B at the time of the event had a separation angle with Earth of 131 and 140°, respectively.

Information about the IMF and solar wind plasma was used, including the magnitude of the IMF and its components (B_x, B_y, and B_z) and the velocity, density and temperature of the solar wind. The in situ measurements of the solar wind characteristics were taken from instruments aboard the Wind spacecraft (Acuña et al. 1995), located in the L1 Lagrangian point. We used them to estimate the arrival of the CME and associated CME-driven shock wave to Earth.

In situ measurements of energetic proton intensities, obtained by two different instruments, were used to study the associated SEP event. We employed the particle observations at 15 energy channels between 2.2 MeV and 58.4 MeV recorded by both detectors of Energetic Relativistic Nuclei and Electron (ERNE; Torsti et al. 1995) instrument, Low-Energy Detector (LED) and High-Energy Detector (HED), all aboard SOHO. We also used data from Low-Energy Magnetic Spectrometer 120 (LEMS120) of the Electron, Proton, and Alpha Monitor (EPAM; Gold et al. 1998) on board the Advanced Composition Explorer (ACE).

2.2. Models

The first model used for this study was the GCS model, which provides the 3D morphology, position, and kinematics of CMEs. It uses observations from at least two coronagraphs at different positions in order to constrain the CME shape, using two conical legs joined by a pseudo-circular front with a circular cross-section. A full description of the model can be found in Thernisien (2011).

The solar wind, CMEs, and a CME-driven shock wave were modelled with EUHFORIA, a 3D MHD heliospheric model focussed on both the research and operational space weather forecasts. EUHFORIA is a data-driven model that consists of a semi-empirical coronal module and an MHD heliospheric module. The coronal module computes the magnetic field topology from the solar surface up to 0.1 au, the so-called inner boundary of EUHFORIA, by combining the potential field source surface (PFSS) model with the Schatten current sheet (SCS) model. The starting point and the main input to the coronal module are synoptic magnetograms. In a next step, the Wang-Sheeley-Arge (WSA; Arge et al. 2003) empirical relation uses the information of the calculated coronal magnetic field to generate the inner-boundary conditions for the heliospheric module of EUHFORIA. This module solves the 3D time-dependent ideal MHD equations in the heliocentric Earth equatorial (HEEQ) coordinate system. In this study, we employed EUHFORIA with an angular resolution of 1° in the longitudinal and latitudinal direction, and a radial resolution of 0.0039 au, with a modelling domain extending up to 4 au.

To model the SEPs, the recently developed particle transport model PARADISE (Wijsen et al. 2019a; Wijsen 2020) is used. PARADISE evolves energetic particle distributions in the solar wind configurations generated by EUHFORIA by solving the FTE. Specifically, the FTE is solved in a stochastic manner, by using a time-forward Monte Carlo simulation. The time-dependent FTE solved by PARADISE includes the effect of the large-scale solar wind structures, modelled by EUHFORIA, on the SEP distributions. Additionally, the FTE takes into account the effects of solar wind turbulence through a set of phase-space diffusion processes (see e.g. van den Berg et al. 2020, for a recent review). In this work, we include a pitch-angle diffusion process and a spatial diffusion process perpendicular to the background magnetic field, the so-called cross-field diffusion. We refer to Wijsen et al. (2019a) and Chapter 3 of Wijsen (2020) for a more detailed description of the PARADISE model.

3. Event description

The focus of this study is on the SEP event starting on 2013 March 15. An enhancement in the particle count, starting at 19:00 UT (top panel of Fig. 1) was observed by ERNE. This initial fast increase was more evident in the lower and intermediate energy channels (from 2.2 to 17.0 MeV and from 17.0 to 32.0 MeV, respectively). After the initial increase, the particle count remained constant for approximately 4 h. A second increase is registered at 02:00 UT on March 16, this time seen in all available channels. The lower energy particles showed a gradual increase until the CME-shock passage time. On the other hand, the higher energy channels (from 32.0 to 68.0 MeV) showed a fast increase, followed by a gradual decay (see the top panel Fig. 1). Observations exhibit smooth and simple time profiles for the higher energy particles. As the energy decreases to intermediate and lower levels, the profile became more complex, with the peak intensity shifting towards later times with respect to the start of the particle event.

Fig. 1. In situ data at L1. Top panel: intensity time profiles for solar energetic protons observed by ERNE aboard SOHO for energy channels ranging between 2.2 MeV and 68.0 MeV (labels correspond to the average energy of the channels). Lower panels, from top to bottom: solar wind speed, density, temperature, and the different components of the magnetic field from the Wind spacecraft. The vertical dotted red line shows the launch time of the Main CME from the Sun. Grey and green shaded areas correspond to the duration of the First and Main CMEs’ effects at Earth, respectively.

The studied particle event was associated with a long-duration M1.1 GOES X-ray flare on 2013 March 15 at 05:46 UT, originating from an active region (AR) with a rather simple configuration of the photospheric magnetic field, NOAA AR 11692 (N11E12). AIA images are unavailable from 06:15 UT to 07:45 UT. After this data gap, the first available images show the existence of irregular, narrow and rather weak coronal dimmings, which can be observed in the 193 Å channel, shown in Fig. 2b.

Fig. 2. SDO/AIA images at different wavelengths. (a) Image at 19.3 nm for 2013 March 12 at 12:15 UT. The coronal hole situated close to the northern polar coronal hole is labelled CH. The coronal hole has a positive polarity, which is opposite to the expected polarity of the nearby northern polar coronal hole. During the year 2013, the polarity inversion was still in progress, so the northern polar coronal hole was very small and hard to detect (Janardhan et al. 2018). (b) Image at 19.3 nm for 2013 March 15 at 07:45 UT. The inverted sigmoid position is highlighted with the yellow box. (c) Image at 9.4 nm for 2013 March 15 at 06:03 UT. The inverted sigmoid position is again highlighted with a yellow square.

The CME associated with the M1.1 flare (hereafter referred to as Main CME) was an asymmetric halo CME first seen in the LASCO C2 coronagraph field of view at 07:12 UT. Its main propagation direction had a slight north-east component, its projected plane-of-the-sky velocity was 1060 km s⁻¹. This Main CME was also observed by both STEREO spacecraft. The halo CME was first observed at 07:39 UT in the STEREO COR2-A field of view, and at 07:09 UT in the COR2-B field of view, above the west solar limb as seen in STEREO-B. All of the above, together with the positioning of the spacecraft indicates that the CME had a strong component of its 3D velocity towards Earth.

All the CMEs observed within a 24-h period after the Main CME were slower and narrower in comparison to the studied CME. As these succeeding CMEs could not influence the Main CME we did not consider them in the study. We also examined CMEs in the time window of 4 days prior to the Main CME in order to understand if the Main CME interacted with some of the earlier and slower CMEs.

The first additional CME we analysed was ejected from the Sun on March 12 at 11:12 UT (as seen in SOHO/LASCO C2 field of view). The associated flare was a GOES C2.0, with the peak intensity registered at 11:06 UT, originating from NOAA AR 11 690 positioned at N24 E04. This event will be referred to as the First CME. The second additionally studied CME was ejected from the Sun on March 13 at 00:36 UT (as seen in SOHO/LASCO C2 field of view). The associated flare was a GOES C1.9, with the peak intensity registered at 00:17 UT. The flare originated from NOAA AR 11 692 positioned at the moment of flaring at N09 E35. This event will be referred to as the Second CME.

In Fig. 1 the in situ observations by Wind spacecraft for the period of interest are shown (panels 2 to 8). Before the arrival of the First CME (marked by the grey-shaded area), the solar wind speed value was about 300 km s⁻¹. An interval of increased density, starting on March 14 at about 12:00 UT, was followed by the velocity increase (at about 22:00 UT) that indicates the arrival of the solar wind originating from the coronal hole. The source of this wind was probably the elongated and narrow positive polarity coronal hole, stretching from about 30 to 70° north (seen in Fig. 2a). The arrival of the First CME is registered at 04:40 UT and lasts for about 6 h.

After the First CME passed and until March 16 the solar wind speed remained very stable with the value of about 450 km s⁻¹. On March 16, a weak velocity enhancement of about 25 km s⁻¹ was registered. The enhancement in velocity lasted around 9 h and was accompanied by fluctuations in the temperature, magnetic field, and density. Although the changes were relatively small, this behaviour of in situ data indicates the presence of a magnetic structure, which could possibly be associated with the Second CME. After the passage of this magnetic structure, the solar wind speed decreases gradually until reaching 400 km s⁻¹.

The Main CME-driven shock arrived at L1 on March 17 at 05:22 UT (marked with the beginning of the green shaded area). The in situ signatures of the Main CME started with the shock arrival, characterised by the sudden simultaneous increase in the plasma parameters (i.e. solar wind velocity, density, temperature, and total magnetic field magnitude). Velocity and density increased up to 750 km s⁻¹ and 28 particles cm⁻³, respectively.

Magnetic field components show very small variations (of the order of few nT) until the arrival of the main CME-driven shock. During the passage of the sheath region, following the shock arrival, all magnetic field components (B_x, B_y, and B_z) showed fluctuations on the level of ±10 nT. The B_x and B_y components were stable, having values of about +7 nT and −7 nT, respectively, over around 10 h. Both of the components showed a rotation starting at about 23:00 UT on March 17 and ending at about 10:00 UT. Similarly to the other components, the B_z component also shows fluctuations during the passage of the sheath region, reaching −20 nT. The sheath region was followed by a smooth rotation of B_z (from −10 nT to +10 nT) starting on March 17 at about 17:00 UT and ending on March 18 at about 05:00 UT. We consider the end time of this interplanetary coronal mass ejection (ICME) was about 60 h after its arrival.

4. The modelling of the CMEs

In this section we discuss the modelling of the Main CME event and two preceding CMEs that potentially influenced its propagation.

4.1. GCS model

A GCS fit (Thernisien 2011) was performed for all three studied events. Calculated parameters are presented in Table 1. Figure 3 shows the GCS fitting results obtained with coronagraph images from all three spacecraft (STEREO-COR2A and COR2B and SOHO/LASCO-C2). The fitting was done at three different times, keeping constant all the parameter values except for the height, which was used for the calculation of the CME propagation speed (starting at 07:36 UT and ending at 08:42 UT).

Fig. 3. White light coronagraph images. Top panels correspond to the GCS fitting (green mesh) for the Main CME as seen by (a) LASCO-C2, (b) STEREO-COR2A, and (c) STEREO-COR2B. The bottom panels correspond to the images of (d) LASCO-C2, (e) STEREO-COR2A, and (f) STEREO-COR2B, used for the fitting procedure.

Quantifying the GCS fitting errors is a difficult task. Geometric parameters coming from the 3D reconstructions can strongly vary, and the observer’s expertise is among the main sources of errors (Thernisien et al. 2009; Verbeke et al. 2022a; Singh et al. 2022). Our GCS results were used in EUHFORIA as input for two different CME models. We performed several EUHFORIA simulations, in order to optimise the GCS parameters needed for accurate modelling of the in situ observations. This iteration process was performed in order to obtain as accurate as possible solar wind and CME modelling inputs for PARADISE (see Sect. 5) since the focus of this work is on the modelling of the SEP event. The CME event on 2013 March 15 (Main CME) was also studied by Pal et al. (2017) with the focus of the study on the low coronal modelling of the event and the estimation of the magnetic helicity of the flux rope. The authors also used the GCS model for constraining the CME parameters, and they obtained somewhat different but comparable results.

4.2. Modelling with EUHFORIA

Accurate modelling of the solar wind characteristics with EUHFORIA is a complex task. In order to do so, we explored different configurations for the different building blocks of the model. This was done in the particular order presented below.

4.2.1. Coronal model optimisation

Magnetogram selection. In order to optimise the background solar wind modelling we used GONG magnetograms and ADAPT maps (provided by Global Oscillation Network Group and Air Force Data Assimilative Photospheric Flux Transport, respectively) within a time window up to 3 days prior to the Main CME. EUHFORIA simulations with the default coronal model parameters and GONG magnetograms do not capture any high speed stream (HSS) arriving on 2013 March 14. On the other hand, when using the GONG-ADAPT magnetogram for 2013 March 12 at 14:00 UT, the model manages to capture signatures of the HSS, comparable to the ones seen in situ.

PFSS and SCS. Asvestari et al. (2019) present a comprehensive study on the influence of different combinations of the height R_ss of the source surface radius (in the PFSS) and R_i (the inner boundary of the SCS model) that are employed in the coronal model of EUHFORIA. The authors concluded that the default pair of values – R_ss = 2.6 R_⊙ and R_i = 2.3 R_⊙ for SCS – does not reproduce correctly the background solar wind and that lowering the R_ss and R_i heights yields better results. We performed a similar analysis for the herein studied events. Simulation results using the default parameters underestimate the background solar wind speed and overestimate the density. In order to increase the background solar wind speed while also getting more open field lines into the heliosphere, we lowered the position of these two surfaces and used 1.810 R_⊙ and 1.805 R_⊙ for the PFSS and SCS, respectively. These yield more accurate solar wind modelling results, as described in Sect. 4.2.3.

Other parameters. For further improvement in the accuracy of solar wind modelling, the minimum fast solar wind speed output at the inner boundary of EUHFORIA’s coronal model was lowered to 300 km s⁻¹ from the default value of 330 km s⁻¹. We also reduced the fast solar wind plasma density input at EUHFORIA’s inner boundary to 250 cm⁻³ (default value is 300 cm⁻³). These additional modifications of EUHFORIA’s input resulted in the better reproduced HSS arrival (see Sect. 4.2.3). It is important to note that all the steps described here are not done when using EUHFORIA in an operational environment.

4.2.2. LFF-spheromak parameters

After optimising the background solar wind modelling, CMEs are inserted into the heliospheric domain of EUHFORIA to have a complete representation of the heliosphere at the time of the SEP event. EUHFORIA has a few different CME models implemented, such as the cone model (Xie et al. 2004), FRi3D (Maharana et al. 2022), and the LFF-spheromak (Verbeke et al. 2019). Each of the models has its strengths and weaknesses. As the FRi3D model is very computationally demanding, the LFF-spheromak and the cone model were selected for this study. The cone model treats ejecta as a hydrodynamic pulse with its geometry characterised by self-similar expansion. This model does not contain any intrinsic magnetic field, and the CME parameters needed as input are the time of insertion at 21.5 R_s, density, temperature, velocity, half-width, latitude, and longitude are presented in Table 1. The last four are calculated using the GCS model. The LFF-spheromak model considers the CME to be a magnetised sphere with a predetermined radius when it is injected into the heliospheric domain and needs nine predetermined parameters: latitude, longitude, the radius of the CME at 21.5 R_s, 3D speed, tilt angle, density, temperature, helicity sign, and toroidal flux. In order to calculate the first five parameters, the GCS fitting was used and the results are presented in Table 1. Density and temperature are prescribed at 10⁻¹⁸ kg m⁻³ and 0.8 MK, respectively. The last two parameters are obtained from the observations.

It is important to note that the definition of the tilt angle in the LFF-spheromak model differs from the one calculated by GCS (see e.g. Asvestari et al. 2022). For this event, the magnetic axis of LFF-spheromak is aligned with the GCS fitting results, but due to the negative helicity in all modelled CMEs, the tilt angle set for the model is multiplied by −1.

When using the GCS model for constraining the CME and obtaining its main characteristics, we are assuming that all parameters of the CME, except the height, remain constant. This is not always true, as deflections or deformations of the CMEs can occur in their transit up to 21.5 R_s. These effects are not included in either the CME fitting model nor EUHFORIA simulations. Images from the different coronagraphs suggest that the main direction of propagation of the First CME is towards the north-east (as seen by both LASCO C2 and C3 coronagraphs). We extracted parameters from the GCS model and include them in the different CME models available in EUHFORIA. Preliminary simulations indicate that this event is not reaching Earth’s position, while the in situ data show the arrival of a CME together with the HSS. This suggests that with EUHFORIA we are not able to reproduce some deflections and changes in the CMEs direction of propagation that happened in reality. To account for it, we use the geometric parameters for this event and aim at fitting just the leg of the CME using the GCS geometry as presented in Thernisien (2011). So the First CME’s GCS fit results presented in Table 1 are adapted to this specific case.

We do not address the magnetic flux calculations in detail as it is not the main focus of the study. The authors in Pal et al. (2017) explored two ways of calculating the magnetic reconnection flux, using AIA 193 Å and GOES Solar X-ray Imager observations. Employing these techniques they estimated an upper and lower boundary of the magnetic flux and the average value of these two values was used in our study. If the reconnected flux is assumed to go into the poloidal magnetic flux, then the toroidal flux is calculated using Eq. (13) from Scolini et al. (2019). For the First and Second CME, we used the Kazachenko et al. (2017) flare ribbon database². Employing the same relation as for the Main CME, we converted the reconnected flux into toroidal flux.

4.2.3. Background solar wind with CMEs

In order to understand the effects of the different CME models on particle propagation, three different EUHFORIA simulations were made, keeping the modelling of the First and Main CME unchanged. The Second CME is slow and has a radius comparable to the Main CME, as shown in Table 1. This makes it a good candidate for generating a blockade for propagation of particles accelerated by the Main CME. The first simulation (Case 1) was aimed at generating a baseline on how the particles would behave when only the First and Main CMEs are included. For the second and third EUHFORIA simulations, we included the Second CME using the cone model (Case 2) and the LFF-spheromak (Case 3), respectively.

Figure 4 shows a snapshot at the same simulation time for all considered cases. Panels a and b correspond to Case 1, panels c and d correspond to Case 2 and panels e and f correspond to Case 3. The figure’s top row shows the solar wind radial velocity, while the figure’s bottom row displays scaled number density in the equatorial plane. The white, red and black arrows mark the approximate position of the First, Second and Main CME. The red circle corresponds to the Earth’s position while the green square and triangle show the position of two virtual spacecraft positioned 15° east and west of Earth, respectively. For all cases, the First CME is already beyond Earth’s position and interacting with the HSS. At the same time (for Cases 2 and 3), the Second CME’s main propagation direction is characterised by a southward and eastward components (see Table 1). This makes it a flank encounter at Earth in both cases. At this point in the simulation, the Main CME is positioned at ∼0.5 au, along the Sun-Earth line. Figure 5 shows modelling results together with the in situ data. The results for Cases 1–3 are presented with the green, blue and red curves, respectively. Figure 5 shows an accurate arrival time for the CME-driven shock of the Main CME, with all of the runs modelling the arrival within a 2-h time window with respect to observations (see the top panel of Fig. 5).

Fig. 4. Snapshot of the EUHFORIA simulations for Case 1 (a and b), Case 2 (c and d), and Case 3 (e and f) for 2013 March 16 at 06:44 UT. The top row shows the radial solar wind speed, and the bottom row shows the scaled number density in the equatorial plane. White, red, and black arrows point at the approximate position of the First, Second, and Main CMEs (respectively). Dashed lines correspond to the magnetic field line connection from the E15, Earth, and W15 positions, backtracked to 0.1 au. The temporal evolution for (a), (c), and (e) are available as an online movie.

Fig. 5. EUHFORIA simulations compared with in situ data. Top to bottom: solar wind speed, density, temperature, total magnetic field, and its three components. In all panels, the observations from Wind spacecraft (black dots) and EUHFORIA simulation results are given for Case 1 (green), Case 2 (blue), and Case 3 (red). The vertical dotted red line shows the launch time of the Main CME event from the Sun. Grey and green shaded areas correspond to the duration of the First and Main CME’s impact at Earth.

All simulations show a slight overestimation of the modelled speed after the shock arrival, with the red line peaking ∼50 km s⁻¹ higher than the blue line (with the green line in the middle). The same behaviour is also observed for other variables, namely density and temperature, which show accurate arrival time but slightly overestimated amplitude.

The HSS is reproduced by the model but not very accurately for all cases, with a late arrival and underestimated speed increase. The modelled density is strongly overestimated. We think that the First CME event piled up material in front of it. As the arrival in the simulations is a glancing blow that coincides with the arrival of the HSS, this enhances even further the density, without a clear jump in speed.

The magnetic field components (B_x, B_y, and B_z) show the arrival of the Main CME in good accordance with the in situ data. Correlation is especially good for the simulations with 3 LFF-spheromaks (Case 3), where the rotation in all three components begins earlier than for the other cases presented with a blue and green line (Cases 1 and 2). The total magnetic field (B_tot) shows a late arrival for the blue and green lines and no overestimation in the amplitude for all the cases. The B_x component is not well reproduced by any simulation setups, as after the sheath passage the data show positive values, and simulations remain to have negative values.

The modelled B_y magnetic field component shows a dip in intensity similar to the observations after the sheath region. A shift to the positive values earlier than in the case of observations can be noticed in all cases. For the B_z component, simulations show a similar intensity at the minimum value, with the red curve reaching the minimum value before the others. After the sheath, all simulations show a rotation, which starts earlier than in the in situ observations but neither of them reaches positive values.

5. The modelling of the SEP event

The three MHD simulation cases described in the previous section were also used as input for PARADISE to model the SEP event. In this section we first describe the used particle emission source, and then how we determine the location of the particle injection for the Main CME, and the particle diffusion conditions used in the model. Next, we describe and discuss the resulting spatial distribution of SEPs and the proton intensity-time profiles observed by three virtual observers located at 1 au and within 30° in longitude. This illustrates how strongly the shape of the intensity profiles can vary in the different modelled solar winds.

5.1. Simulations with PARADISE

5.1.1. The source of particles

To mimic the emission of shock-accelerated protons, we inject particles continuously at the Main CME-driven shock, with an energy spectrum that follows a power law with an exponential rollover:

$$\begin{array}{cc}& Q(E,r)\propto E^{-2}exp(-\frac{E}{E_0(r)}){\left(\frac{r_0}{r}\right)}^2\mathrm{\nabla }·{\mathit{V}}_{\mathrm{SW}},\hfill \end{array}$$(1)

where

$$\begin{array}{c}\hfill E_0(r)=5\mathrm{MeV}\left(\frac{r}{r_0}\right).\end{array}$$(2)

In this equation, E denotes the particle’s kinetic energy, V_SW is the solar wind speed, r is the heliocentric radial distance, and r₀ = 1 au. The power-law index in Eq. (1) is determined by fitting a power law, E^−γ, to the 0.047−6.4 MeV proton differential intensities measured at the in situ shock-passage time (05:22 UT on 2013 March 17). For this fitting, we used the eight energy channels of the EPAM LEMS120 (Gold et al. 1998) on board ACE and the first seven nominal proton channels of the LED of SOHO/ERNE (Torsti et al. 1995). The observed energy spectrum at the shock decreases faster for E > 5 MeV, and hence we took this energy as the exponential turnover energy at 1 au, to simulate the spectrum of the emitted particles by the shock.

The inverse dependence of the roll-over energy E₀ on r aims at mimicking the gradual weakening of the particle acceleration efficiency of the shock wave as it moves away from the Sun. This means that at the inner boundary of EUHFORIA, r = 0.1 au, the energy spectrum of the particle emission source extends as a power law almost up to 50 MeV protons. In this way, we ensure the presence of protons as energetic as those observed at 1 au during the 2013 March 15 event (see the top panel of Fig. 1) in the simulations. In a nominal solar wind configuration, the first magnetic connection for an observer at Earth would be established at the west flank of the shock driven by the main CME. However, by this point, the shock is already several solar radii away from the Sun, and it may not be efficient enough to accelerate protons to a few tenths of MeV. This interpretation is supported by the lack of a prompt increase in the particles in the observed event (see the top panel of Fig. 1) and by the inverse velocity dispersion observed at the onset of the 14−32 MeV proton intensities measured by SOHO/ERNE/HED on March 15. Hence, the constraint of injecting the CME at 21 R_⊙ may be much less critical in our event than for some other well-connected cases.

The modelled particle intensities, described in the following sections, are translated to physical units using three different normalisation factors. We scaled the 2.4 MeV intensity value to the observed value in the 2.2−2.7 MeV channel of SOHO/ERNE/LED at the time of the shock passing by the Earth. The scaling procedure is conducted individually for each case to ensure a fair comparison with the observations of the SEP event corresponding to each specific solar wind scenario. For Case 2, the ratio value used in comparison with the normalisation factor for Case 1 is 0.58, while for Case 3, the ratio value is 2.52.

The dependence of Q on the divergence of the solar wind velocity and the heliocentric radial distance in Eq. (1) has been previously used by Wijsen et al. (2021, 2022, 2023) to successfully reproduce the < 10 Mev intensities of gradual SEP events and an energetic particle event associated with a corrotating interaction region, respectively. We refer the reader to those previous works and to Prinsloo et al. (2019) for a detailed description.

5.1.2. Location of the shock

To locate the shock in the EUHFORIA’s 3D data cube at one selected time instant, we applied the automated shock tracer described in Appendix A of Wijsen et al. (2022). Firstly, this code subtracts the unperturbed background solar wind (i.e. the solar wind without any CMEs) at a selected time instant. As a result, the plasma variables become zero everywhere except in the regions containing CMEs.

In the next step, the shock position is determined by identifying where the divergence of the solar wind velocity vector becomes negative and reaches a specified tolerance level. Unlike the simulations discussed in Wijsen et al. (2022), our simulations involve an additional complication due to the presence of preceding CMEs that perturb the background solar wind and generate their own shock waves. To isolate the shock driven by the primary CME, the shock tracer estimates the maximum radial distance of the shock nose at each time instance and disregards shock waves beyond that distance.

Figure 6 shows a snapshot of the EUHFORIA background solar wind speed, in which the black curve indicates where the tracing algorithm has located the shock of the main CME. This figure illustrates the performance of the algorithm, which succeeds reasonably well at locating the shock of the primary CME. We assume that this shock is the main source of energetic particles observed near Earth and, consequently, shocks generated by other CMEs or by interactions of solar wind streams were not traced.

Fig. 6. EUHFORIA background solar wind snapshot with shock tracer function applied at the time instant 10:14:22 UT. White, red, and black arrows point at the approximate position of the First, Second, and Main CME (respectively).

5.1.3. Diffusion conditions

The turbulence that is omnipresent in the solar wind affects the transport of particles through interplanetary space (Dmitruk et al. 2003; le Roux 2023). In the FTE, the effect of turbulence on particle transport is taken into account through a set of diffusion processes in phase-space. The PARADISE simulations presented here include two of such diffusion processes: Acuña et al. (1995) a diffusion process in the pitch-angle coordinate and Aran et al. (2005) a cross-field diffusion process, that is, a spatial diffusion process perpendicular to the background magnetic field. For the pitch-angle diffusion process, we utilised the results of quasi-linear theory (Jokipii 1966) as specified in Wijsen et al. (2019b). The magnitude of the pitch-angle diffusion coefficient was fixed by assuming a constant parallel mean free path λ_∥ = 0.3(R/R₀)^2 − q au, where q = 5/3 is the spectral index of a Kolmogorov turbulence spectrum, R is the particle rigidity, and R₀ = 97 MV corresponds to the rigidity of a 5 MeV proton. The transport perpendicular to the magnetic field was modelled via the isotropic cross-field diffusion coefficient characterised by a constant perpendicular mean free path λ_⊥.

We simulated two different cases for the perpendicular transport of the particles. In the first case, the ratio between the perpendicular and parallel mean free paths is λ_⊥/λ_∥ ∼ 10⁻³ and in the second case it is λ_⊥/λ_∥ ∼ 10⁻⁴. As λ_⊥ ≪ λ_∥, we ensured that the particles primarily propagate along the IMF lines. The main difference between the two cases is that the cross-field diffusion of particles is more effective in the former than in the latter.

5.2. The effect of the different models for the preceding CME on the spreading of SEPs

A snapshot of the simulated omnidirectional SEP intensities for the energy channel 2.40−3.00 MeV (shown in a colour scale) in the solar equatorial plane is presented in Fig. 7. The snapshot was taken 21.75 h after the launch of the Main CME into the computational domain³. The simulations in the top, middle, and bottom rows of the figure correspond to the EUHFORIA solar winds for Cases 1–3, respectively. The left column of Fig. 7 corresponds to the simulations where λ_⊥/λ_∥ ∼ 10⁻³, while the right column show results for the simulations using λ_⊥/λ_∥ ∼ 10⁻⁴. As for Figs. 4 and 6, the black, white, and red arrows indicate the location of the Main, First, and Second CMEs. The Earth’s position is marked by a circle, while the square and triangle mark the positions of two virtual observers at 1 au, at the same HEEQ latitude as Earth, but at 15° east and west with respect to Earth, respectively.

Fig. 7. PARADISE simulation snapshot for 2013 March 16, 06:44 UT. Panels a and b correspond to Case 1, panels c and d correspond to Case 2, and panels e and f correspond to Case 3 background solar wind. The left column shows the simulations using λ⊥/λ∥ ∼ 10−3, while the right column corresponds to λ⊥/λ∥ ∼ 10−4. E15, Earth, and W15 virtual spacecraft are represented with a green square, a red circle, and a green triangle (respectively). White, red, and black arrows point at the approximate position of the First, Second, and Main CME, respectively. The temporal evolution of each case is available as an online movie.

The two top panels of Figs. 7a and b show the spatial distribution of particle intensity for Case 1, where the First CME was the only preceding CME included in the simulation. It is important to note that the LFF-spheromak model was used to simulate the preceding CME, and the IMF in front of the Main CME’s centre-to-western wing is connected to this First CME. Particles that propagate upstream of the Main CME encounter thus the spheromak of the First CME, which is characterised by a (mostly⁴) closed magnetic topology. A consequence of this closed magnetic topology is that the energetic particles can only access the spheromak because of the cross-field diffusion, which is a much slower process than the parallel transport (recall that λ_⊥ ≪ λ_∥). This effect is less noticeable (although not absent) in panel a of Fig. 7, where the First CME is filled with particles earlier, since in that simulation λ_⊥ is one order of magnitude larger. Upstream of the eastern wing of the Main CME, the absence of the Second CME in the EUHFORIA simulation causes the IMF to have a Parker-like structure, resulting in a smooth gradient in the particle intensities observed in panels a and b.

In Case 2, panels c and d of Fig. 7, the Second CME is also included and simulated with a cone model. In this case, the intensity of the energetic particles in front of the Main CME is enhanced in a more extended region towards the west than in Case 1 (panels a and b). Downstream of the shock of the Second CME, there is a compression region, where the magnetic field intensity shows increased values, acting as a magnetic mirror for the particles. Therefore, the higher intensities can be attributed to particles being trapped in a bigger region between the Second and the Main CMEs, and not only between the First and the Main CME as in Case 1. The increased particle trapping in Case 2 allows a greater number of particles to have more time to diffuse into the LFF-spheromak structure of the First CME, in contrast to Case 1. Consequently, the intensities observed in panel d of the LFF-spheromak for the First CME exhibit similarities to those in panel a, despite the fact that λ_⊥ is an order of magnitude smaller. In addition, the interaction between the First and Second CMEs may also modify the topology of the magnetic field in such a way that more particles get access to the First CME. In addition, we notice that the cone CME model of the Second CME modifies the structure of the IMF in such a way that it prevents a portion of the particle population from propagating through it. This allows more particles to interact during a more extended period of time, which enables the particles to fill the LFF-spheromak structure of the First CME (see panel d), with similar intensities as in panel a, but with a ten times smaller λ_⊥. In the case of the larger λ_⊥ (panel c) the intensities inside the First CME (white arrow) are only about one order of magnitude lower than in the surrounding area due to the larger cross-field diffusion simulated. Moreover, panels c and d also show that the Second CME cone model contributes to spreading the particles to farther longitudes than in Case 1. This is consistent with the study by Lario & Karelitz (2014), where they concluded that the presence of preceding CMEs while an SEP event is developing, on average, have larger peak intensities than the events observed in isolation.

In Case 3, panels e and f of Fig. 7, the Second CME is simulated with an LFF-spheromak CME model. The intensity levels attained at distances beyond the position of the Second CME are significantly lower than in the previous two cases. This is due to the complex interaction of the Second CME’s spheromak with the solar wind’s magnetic field and the First CME’s spheromak. In addition to the sheath region behind its shock, the big extent of this LFF-spheromak has two distinct effects on particles. First, it blocks some particles from crossing it, similar to what happens in Case 1 for the First CME. Additionally, as the IMF interacts and reconnects with field lines with non-zero latitudinal component inside the structure, particles move away from the equatorial plane. As a result, the streaming of particles upstream of the Main CME in the solar equatorial plane is significantly obstructed by the presence of the LFF-spheromaks.

5.3. The effect of the preceding CME models on the shape of the SEP intensity-time profiles

Earlier observational and simulation studies show that the intensity-time profiles of SEP events observed within a narrow range of helio-longitudes may strongly vary depending on the different solar wind streams in which the observers are embedded during the development of the short-duration SEP events (Pacheco et al. 2017) or due to the presence of a high-speed solar wind stream affecting both the acceleration and transport of the particles during the gradual SEP events (Lario et al. 2022; Wijsen et al. 2023; Ding et al. 2022b). In this section we focus on understanding the variations in the modelled proton intensity-time profiles that were recorded by spacecraft separated by only 30° in longitude. These variations were caused by the presence of the preceding CMEs, which caused the solar wind to strongly deviate from nominal conditions, as discussed in the previous sections. We analyse the modelled proton intensity-time profiles for the three virtual spacecraft marked in Figs. 4, 6, and 7.

The three rows of the panels in Fig. 8 present the time profiles obtained from simulations corresponding to EUHFORIA solar wind Case 1 (top row), Case 2 (middle row), and Case 3 (bottom row). Each panel presents time profiles of, in the top proton intensity, and on the bottom the first-order parallel anisotropy⁵. The time profiles are shown for every other simulated differential energy channel for the sake of clarity, and the labels in the inset at the top of the figure indicate their geometric mean value. The three columns correspond to different virtual spacecraft: E15 (left), Earth (centre), and W15 (right).

Fig. 8. PARADISE simulated proton time-intensity profiles (top panels) for energy channels ranging from 2.2 MeV to 51.1 MeV, with their corresponding first-order parallel anisotropy profiles (bottom panels). The top, middle, and bottom rows correspond to Cases 1–3 background solar wind for runs where λ⊥/λ∥ ∼ 10−3 at the different virtual spacecraft: E15 (left), Earth (centre), and W15 (right) columns. Dashed green, blue, and red lines correspond to the time of the Main CME-driven shock at the position of the virtual spacecraft for each of the cases. The vertical dash-dotted black line marks the time of injection of the Main CME at the inner boundary of the EUHFORIA model. The bottom red (blue) line indicates the periods of positive (negative) values of the radial component of the simulated magnetic field at the position of each spacecraft.

The 2013 March 15 SEP event was simultaneously detected by STEREO B, located 140° west and 5° north from Earth and at the distance of 1.04 au from the Sun. In our simulations, only the case with the largest perpendicular mean free path results in a particle enhancement at the position of STEREO-B. Therefore, the profiles displayed in Fig. 8 correspond to the simulations with λ_⊥/λ_∥ ∼ 10⁻³. For completeness, we provide the resulting profiles for the lower cross-field diffusion case in Fig. A.1, where we discuss the differences in the profiles due to the two different values adopted for the perpendicular mean free path.

Figure 8 shows that for each observer (different columns) and after the initial phase of the particle event (approximately after 12:00 UT on March 16), the energy spectrum of the particle intensities is very similar (i.e. up to ∼32 MeV) for the three cases (different rows). For the higher energies, the energy spectrum is softer in Case 3 (bottom row), as high-energy particles escape away from low latitudes more easily than the lower-energy particles because of their larger mean free paths. The similarity of the energy spectrum obeys the use of the same energy dependence in the simulated particle emission source (Eq. (1)). In fact, the energy spectrum of the event-fluence (up to ∼32 MeV, not shown here) is similar for all cases and across the three observers. This similarity indicates that the impact of the different solar wind backgrounds on altering the energy spectrum of particle intensities is relatively minor compared to the influence of the particle emission source. This finding is consistent with some previous results. In particular, Aran et al. (2005, 2008) conducted 1D particle transport simulations, investigating the impact of various particle mean free paths, and obtained results that align with our study. The results of a study by Ding et al. (2022b), in which a 2D particle transport model that accounted for cross-field diffusion alongside parallel transport was utilised, also emphasised the significance of the particle emission source on the observed energy spectrum.

For each case (each row in Fig. 8), the onset time of the particles occurs later for the E15 observer (left column) than for the Earth (middle column) and the W15 (right column) observers. The reason for this later onset is that the E15 observer establishes the magnetic connection with the Main CME shock approximately 20 h after the Main CME is injected into the simulations. This delayed magnetic connection is better inferred from the corresponding intensity-time profiles in Fig. A.1, as in that simulation the effect of the cross-field diffusion of particles is reduced. In Fig. 8 instead, we observe an earlier onset of the particles for the E15 observers due to the more efficient perpendicular transport in these simulations. The magnetic connection for the Earth and W15 observers is established earlier, approximately 2.5 h and 1.5 h after the Main CME injection, respectively. The magnetic connection of the observer with the Main shock not only improves in time from the E15 to the W15 observer but also the observers establish a connection with different portions of the shock front with increasing efficiency in particle emission (from E15 to W15), as expected under unperturbed solar wind conditions (e.g. Heras et al. 1994, 1995; Lario et al. 1998; Fry et al. 2005). The connection to stronger regions of the shock front results in less noisy time profiles at the beginning of the particle events for the Earth (middle panels) and even less for the W15 observers (right panels) when compared to the E15 virtual spacecraft (left panels). Also, the intensity-time profiles show more prominent prompt components for the W15 observers in Cases 1 and 2. Contrary to what we can expect in nominal or slightly perturbed solar wind conditions, the onset times for the W15 observers occur at the same or slightly later times than for the Earth-observers, despite their earlier magnetic connection. The reason for this is that, at the initial stages of the particle event, the First CME passes through the location of the W15 observer, causing a significant perturbation in the magnetic field line crossing this particular observer. As a result, the particles reaching the W15 observer follow a longer trajectory along the IMF lines compared to the particles travelling towards Earth.

For the E15 and Earth observers, the first-order parallel anisotropies show a similar evolution for the three cases, once there are enough statistics in all energy channels. This similar evolution suggests that the transport of the particles is not largely affected by the presence of the second preceding CME, which was launched approximately 34° westwards from the W15 observer (see Table 1). Moreover, we have also found that employing the specific CME model also leads to significantly different particle transport conditions.

Figure 8 clearly shows that for the E15 observer, the inclusion of the Second CME shifts the onset of the simulated particle event earlier in time, regardless of the CME model employed. The magnetic connection for Case 1 is established on March 16 at 06:30 UT, approximately 1.5 h later than for Cases 2 and 3. The earlier onsets seen in Figs. 8 and A.1 are due to particles transported perpendicularly to the IMF line connecting with the E15 observer. Particles reach this IMF line, by perpendicular diffusion, both from distances beyond the spacecraft and from shorter radial distances, which explains the oscillating signs of the anisotropy of the first arriving particles, for each energy channel.

For Case 1, after the onset of the lower energy particles (see the top-left panel of Fig. 8), the proton intensity-time profiles show a plateau that starts with the arrival of the First CME’s shock and lasts approximately half a day, during the passage of the associated ICME. The plateau-shaped region is followed by a sharp rise in intensity at lower energy channels (≲13 MeV). As we move to more energetic channels, the increase becomes less pronounced. This second increase corresponds to the crossing of the HSS, as shown in Fig. 4, which causes a better magnetic connection to the CME-driven shock due to the smaller curvature of the IMF in the fast solar wind. As mentioned in the previous section, the closed magnetic topology of the First CME traps particles between the First and Main CMEs. This particle confinement leads to increased intensity profiles for the lower energy channels (≲50 MeV) and a constant profile in the highest energy channel (≳50 MeV). The effect of the particles’ entrapment is also observed in Cases 2 and 3, but the profiles show a smoother increase due to the observers’ different magnetic connection to the shock (caused by the presence of the Second CME that distorts the IMF). Despite the presence of reflected particles, the anisotropy remains low (< 0.5) but still positive, which results from the continuous emission of particles from the Main shock as it approaches the observer (Heras et al. 1994). We also want to point out that the particle anisotropy reverses sign for all energies (within minutes) at the shock passage by all the simulated observers, indicating a continuous connection with the shock (Sanahuja & Domingo 1987). Following the passage of the shock, the negative anisotropies reach significant negative values (∼ − 0.5). These values indicate that the observer, positioned downstream of the shock, maintains magnetic connectivity with the shock structure. The post-shock portion of the proton intensity-time profiles differs greatly among the three considered cases. In Case 1, the peak intensity occurs after 12:00 UT on March 17, and shows a velocity dispersion, indicating (together with the anisotropy information) that particles come from the distant Main shock. In contrast, in Case 2, the low-energy intensities peak at the shock passage. Hence, the Second CME affects the configuration of the IMF, and the connection downstream of the shock is changed with respect to Case 1. Finally, for Case 3, the peak intensity for E < 50 MeV is reached close after the shock passage, peaking at the time the magnetic field intensity reaches its maximum value in the simulation. The shape of the intensity-time profiles depends thus on the presence or not of intervening large-scale structures in the solar wind, and on how these structures are modelled.

For an Earth-positioned observer (middle column in Fig. 8), the three cases show a similar initial phase, with intensity-time profiles increasing fast for all energies. For Cases 2 and 3 more particles reach the virtual spacecraft than for Case 1, because, as explained for the E15 observers, the presence of the Second CME improves the magnetic connection with the Main CME-driven shock. This translates into higher intensities for the lower-energy channels and a slightly faster increase for the high-energy profiles. After the initial enhancement of the intensities, the shape of the profiles gradually changes with the energy of the particle profiles from a continuous increase in the intensity for the lower-energy channels to decreasing intensities for the higher-energy channels (E > 45 MeV). For all cases, the first-order parallel anisotropies show positive values (≲0.5) up to the shock passage. As for E15, the sustained positive values reflect the continuous emission of particles at the shock wave driven by the Main CME, consistent with the Richardson & Cane (1996) analysis of in situ observations.

The differences in the evolution of the proton intensity- and anisotropy-time profiles among the three cases increase as the Main shock approaches the observers. For Case 1 (top middle panel), at around 15:00 UT on March 16, the intensity-time profiles exhibit an increase across all energies that lasts until the shock passage. This increase is accompanied by a slight rise of the anisotropy, which peaks before decreasing to low values again, just before the shock arrival. These changes in the profiles correspond to the end of the passage of the simulated First CME and the connection with a stronger portion of the Main CME-driven shock front as it approaches the observer. The anisotropies show positive values because particles are injected from the distant Main CME shock, but a fraction of the particles previously injected are reflected at the rear of the First CME and travel back, towards the Earth observer, thereby leading to a reduction in the first-order parallel anisotropy.

For Case 2, the increase seen for Case 1 is not noticeable. Instead, the E < 50 MeV proton profiles show a local maximum value before 20:00 UT on March 16, and then continue with a slight decrease, or with a roughly constant profile for the lower-energy channels, up to the shock arrival. This behaviour is due to the presence of the Second CME, which deforms the IMF lines in such a way that some of the particles that are streaming upstream of the Main CME-driven shock are spread over a larger region to the west of an Earth-positioned observer.

On the other hand, in Case 3, the Second CME’s LFF-spheromak traps particles closer to the main CME-driven shock, making the increase in the intensities as the shock approaches the Earth observer seen up to ∼45 MeV, and narrower as compared to Case 1. For 24 < E < 45 MeV (orange to purple curves), the intensities at the shock are of similar value to Case 2. In contrast to Case 2, however, the intensities show a local peak during the passage of the post-shock’s sheath region, at the time the IMF intensity profile shows a maximum value that is larger than for Cases 1 and 2. Hence, the LFF-spheromak model used for the Second CME modifies the magnetic configuration downstream of the main CME-driven shock. Downstream of the Main shock, the three cases show different intensity-time profiles. For Case 1, intensities remain roughly constant until the observer enters inside the magnetic structure of the Main CME, where intensities drop suddenly to almost zero, similar to an extreme Forbush decrease. Such a decrease is also seen in the other two cases, but particles do fill the region inside the Main CME spheromak structure thanks to the IMF modified by the presence of the Second CME. The decrease in the intensities happens slightly faster for Case 2 than for Case 3. This feature is better seen for the W15 observers. Finally, the intensities of the higher-energy channels decay faster for Case 3, after the peak at the prompt phase, due to particles transported to other latitudes.

Finally, the profiles for the W15 observer (right column in Fig. 8) show the larger variation among the three solar wind cases, especially in the anisotropy profiles. Figure 7 clearly demonstrates that among the three observers, the W15 observer is the most significantly impacted by the First CME. At the time of the Main CME insertion into the computational domain, the W15 observer is situated within the magnetic structure of the First CME. This results in a perturbed magnetic connection with the Main CME shock for all three solar wind cases. In Cases 1 and 2, this leads to a smooth prompt phase in the intensity-time profiles at the onset of the event. Conversely, in Case 3, the presence of the Second spheromak CME creates a more complex ambient IMF structure, causing particles to take a longer time to reach the observer. This is more clearly seen in Fig. A.1, where the onset of the intensities occurs noticeably later in Case 3. The magnetic connection to the shock is slightly better in the presence of the Second cone CME than for Cases 1 and 3, explaining the larger peak intensities at the prompt phase. These peak intensities are about one order of magnitude smaller in Case 3 due to particles escaping to other latitudes.

After the passage of the First CME at the W15 observer in Case 1, particle intensities begin to rise, similar to the observer positioned at Earth. However, particles reaching the W15 observer follow a longer path. The IMF line connecting with the shock initially follows a similar path as Earth’s IMF line up to 1 au, but then the IMF extends to greater radial distances and curves inwards towards the W15 observer, located at 1 au (see, for instance, the top-left panel of Fig. 4). This configuration of the IMF causes the positive sign of the particle anisotropies when the IMF points towards the Sun (blue horizontal line in the top-left panel of Fig. 8). Next, as the First CME moves farther beyond the observer, the IMF changes its shape and orientation, with particles streaming outwards from the Main shock towards the observer; hence, positive anisotropy values are obtained. At this time, the proton intensities start to slightly decrease, indicating a lack of significant acceleration or trapping at the shock arrival. This is due to the magnetic connection shifting from the nose towards the western wing of the Main shock.

The profiles for Case 2 evolve similarly to Case 1 until the arrival of the Second CME on ∼04:00 UT on March 16. Then the anisotropy profile changes to positive values, indicating a direct connection with the Main Shock. The intensities for E ≲ 13 MeV protons show a plateau, resulting from the particles trapped between the Main shock and the Second CME. At higher energies, the trapping was not so efficient, and thus the intensity-time profiles decrease.

On the other hand, Case 3 shows a different time-history for both intensities and anisotropies. Early on March 16, the situation is similar to Case 1, exhibiting positive anisotropies. This is because the IMF initially crosses 1 au east of the observer, followed by a loop-back and subsequent crossing of the W15 observer. Next, everything changes with the arrival of the Second CME at ∼06:00 UT on March 16. The magnetic configuration just after the arrival of the Second CME causes the magnetic connection with the Main shock to largely vary for the different observers. First, the W15 observer is connected with the nose of the shock and, hence, the increase in the intensities and the larger positive anisotropies (from 06:00 to 10:00 UT on March 16). Next, at about 10:15 UT, the magnetic connection is established with the western flank of the Main CME shock (where the ∇ ⋅ V_SW is smaller than in the nose of the shock), which explains the decrease in the particle intensities in all energy channels. After the Second CME has passed the W15 observer, particle intensities stop decreasing and start rising again. This behaviour is due to particle entrapment between the Main CME’s shock and the rear of the magnetic structure of the Second CME. The increase at all energies sharpens with the arrival of the Main shock, with all intensities showing a peak value after the shock passage, in coincidence with the IMF intensity’s peak value.

The presence of the Second CME, and more importantly, the CME model used, largely affects the downstream portion of the intensity-time profiles. In the absence of the Second CME, the particle intensities drop a couple of orders of magnitude for all energies when the Main CME LFF-spheromak magnetic structure passes through the observer. This decrease is softened in the presence of the Second CME. When using the cone CME model (Case 2), particles behave similarly to Case 1 after the shock, but the decrease is less steep as more particles reach inside the Main CME magnetic structure. For Case 3, particle intensities start to decay just after the peak value, due to a fraction of particles diffusing into the magnetic structure of the CME.

5.4. Comparison with the observed SEP event

When we compare PARADISE results for an Earth observer with in situ observations (top panel of Fig. 1), we see some similarities, but overall the initial phase of the intensity-time profiles is not reproduced by any of the synthetic instances. In contrast, from approximately 08:00 UT on March 16 up to the shock arrival, the intensity-time profiles compare well with observations when using Case 2 as the background solar wind. The E ≥ 21.6 MeV proton intensities for the W15 observer (middle right panel of Fig. 8) are very similar to the observed profiles, whereas for E < 21.6 MeV the virtual Earth spacecraft matches better the observations (middle centre panel of Fig. 8). In addition, for the Earth Case 2 observer, the peak intensities attained for each energy channel are within a factor of ∼2 from the actual peak values. This suggests that the energy spectrum chosen for the particle source function is rather adequate. For the downstream portion of the profiles, however, the observer that better reproduces the observations is the E15 Case 3 observer.

Regarding the initial phase of the intensity-time profiles, the plateau-shape profile registered by ERNE between 16:00 UT on March 15 and 02:00 UT on March 16 is fairly well reproduced for E < 21.6 MeV by the Earth Case 2 observer. This feature is even better reproduced for the simulation with the smaller perpendicular mean free path (see Fig. A.1). In fact, this initial plateau-shape of the profiles is best reproduced for the E15 Case 2 observer in Fig. A.1, although too late in time.

The high-energy (E ≳ 30 MeV) intensity increase observed on March 16 at ∼01:00 UT is not reproduced by any of the cases. Nevertheless, a similar increase is obtained in our simulations for the W15 observers in Cases 1 and 2 (top and middle right panels of Fig. 8), although the simulated intensities are lower than the observed values. This second increase appears more clearly for the W15 Case 1 observer in Fig. A.1, but still it is low. It should be noted that the instrumental background for the 51.1 MeV channel is higher than the simulated intensities for Case 1. If the magnetic connection with the Main shock – after the First CME passes by the W15 observer – were established closer to the nose of the shock, higher intensities could have been attained during this second increase, and the earlier bump would have remained masked by the background intensities; thus resembling the observations. The smoother profiles obtained with the larger perpendicular mean free path and the somewhat better simulation of the shape of the intensity-time profiles suggest that a better modelling of the perpendicular spatial diffusion coefficient would correspond with a perpendicular mean free path between the two values used in our simulations. Finally, as already mentioned in the previous section, the onset of the simulated SEP events occurs early than in observations because the magnetic connection with the emitting Main CME shock is established early. In the presence of the Second preceding CME, this connection even improves, as it happens closer to the nose of the shock, contrary to what observations suggest. Thus, having an accurate modelling of the region of the space where particles are transported is key to reproduce SEP events.

On the other hand, we used the GCS model to constrain the CME and obtain the parameters needed by the CME models, as shown in Sect. 4. As discussed in Singh et al. (2022), the longitude parameter is susceptible to errors associated with the different spacecraft white-light images available and used for the fitting. The width of the CME calculated with the GCS model is another source of error that can affect the results. Furthermore, in Verbeke et al. (2022a) it is shown that the human-in-the-loop has a quantifiable influence over the accuracy of the fitting. Modifying different parameters in the fitting and including them in the different CME models, might not yield worse simulation results (when compared with in situ data), but the resulting simulated SEP profiles can be strongly affected, as shown in Sect. 5.3. For example, if we overestimated the width of the Main CME, and it was actually narrower, then magnetic connection with Earth would have been established later in time in the simulations, and so the proton intensity onset. The same is true for Cases 2 and 3, if we would have launched the Second CME more towards the eastern limb of the Sun. More importantly, the CME models used in this study are approximations of the real CME shape and internal magnetic structure. As SEPs stream along the magnetic field lines, it is extremely important to capture the magnetic field topology of the preceding CMEs and the ambient solar wind in which particles will propagate. This on its own is a major challenge in heliospheric modelling. The improved magnetic connectivity that the Second CME introduces to observers in the centre-to-western flank of the Main CME, results in earlier onsets of the simulated events. This is another indication that not only the Main CME and its shock may propagate in a different direction than the GCS reconstruction suggests, but that these errors can also apply to the preceding CMEs. Furthermore, the LFF-spheromak model is likely not an accurate global model of the CME. Instead, this CME model is better at reproducing a rather small part of the magnetic cloud around the centre of the cloud, but not the CME flanks (Scolini et al. 2019). Solar energetic particle modelling is very susceptible to this, as particles propagate through large portions of solar wind on their way to the observer.

6. Summary and conclusions

In this work we have presented a comprehensive study on how preceding CMEs can affect the propagation of energetic particles in the heliosphere, using the EUHFORIA+PARADISE modelling chain. To do so, we selected the gradual SEP event that occurred on 2013 March 15 as a case study. This event is highly complex and involves two preceding CMEs (the First and Second CME) along with a Main CME, which was likely responsible for the observed SEP event (see Sect. 3). As presented in Sect. 4, the GCS model was used to constrain the studied CMEs and obtain the geometric and kinematic parameters necessary for inserting the CME into the heliospheric model. Three different background solar wind conditions were generated by EUHFORIA and used as input for the particle transport model, PARADISE. The First and Main CMEs were inserted into EUHFORIA’s heliospheric domain using the LFF-spheromak model. When the Second CME was considered, two distinct CME models were utilised. Accordingly, we distinguish three different cases:

Case 1: No Second CME.

Case 2: Second CME using the cone model.

Case 3: Second CME using the LFF-spheromak model.

As presented in Sect. 4.2, we find no major difference between the results of different EUHFORIA simulation setups at the location of Earth. Furthermore, EUHFORIA modelling results show a good agreement with in situ observations by Wind spacecraft. These results indicate that the Second CME has no direct impact on Earth, because it does not propagate along the Sun-Earth line (see Table 1).

The three cases of EUHFORIA simulations served as input for PARADISE to model the propagation of energetic protons with an energy range spanning from 2.0 MeV to 58 MeV, as shown in Sect. 5.1. Protons were continuously injected at the shock driven by the Main CME. To determine the position of the Main CME-driven shock in each snapshot, we applied a shock tracing model to the EUHFORIA results.

Two different perpendicular diffusion conditions were used in this work. The ratio of the perpendicular to the parallel mean free paths were assumed to be λ_⊥/λ_∥ ∼ 10⁻³ and λ_⊥/λ_∥ ∼ 10⁻⁴. A smaller perpendicular mean free path, compared to the parallel mean free path, restricts particle diffusion across magnetic field lines, resulting in their predominant propagation direction being along the field lines. Effects of the different cross-field diffusion conditions on the particles’ spatial distribution are presented in Sect. 5.2, and their effects on the intensity- and anisotropy-time profiles are described in Sect. 5.3 and Appendix A.

The main conclusions from this study are:

• Current CME models are capable of reproducing single-spacecraft in situ data. However, their validation at other locations in the heliosphere, crucial for understanding phenomena such as SEP transport, remains challenging due to limited available spacecraft data. This limitation prevents us from assessing the model’s performance in regions where preceding CMEs alter the IMF but do not directly impact any spacecraft. In our study (see Sect. 4.2), we have demonstrated that including preceding CMEs in simulations may not always have a significant impact on the accuracy of the solar wind model when compared to in situ data at the spacecraft location. However, we have also highlighted (see Sects. 5.2 and 5.3) the substantial influence of preceding CMEs on the transport of SEPs. Consequently, the comparison between SEP simulations and observations can be utilised to assess the performance of the solar wind and CME models.

• Depending on the CME model used, the structure of the IMF is modified differently, which in turn influences the transport conditions of the energetic particles. However, the CME models used in this study lack the necessary realism to fully capture the complete CME structure. These models are either too simplistic, such as the cone model, or provide only a local approximation of a flux rope, such as the LFF-spheromak model. Given that CMEs play a crucial role in heliospheric modelling and are the primary drivers of space weather, more realistic models are needed. Although improved flux-rope models, such as those proposed by Maharana et al. (2022), are available, they currently require significantly more computational resources. More efficient, self-similar toroidal CME models (e.g. Linan et al. 2023) are currently being developed and will be included in future versions of EUHFORIA.

• Different CME models require different parameters describing the CME dynamics and geometry. To determine the parameters, simplifying assumptions are needed. As CME models become more sophisticated, the number of parameters required to accurately fit the observations grows. However, this increase in parameters can introduce the possibility of degenerate cases and pose challenges in effectively constraining the parameters. In this work we used the GCS model to constrain the parameters for the cone and LFF-spheromak CME. This model on its own can introduce uncertainties into the CME model, as shown by Singh et al. (2022). Also, the accuracy of this type of fitting is influenced by the subjectivity of the human-in-the-loop, as discussed in Verbeke et al. (2022a). Errors in constraining the kinematic parameters will have effects on the SEP propagation. In this work we show qualitatively how different structures present in the MHD simulations affect particle propagation. The most straightforward features to compare quantitatively are the observed and modelled onset times of the SEP event. In a case when observations show the SEPs with a prompt onset and the model does not manage to reproduce the observations, the most probable reason is that the MHD simulation did not correctly capture the magnetic connection between the CME-driven shock and the observer. In the simple case of quiet solar wind conditions (i.e. in the absence of preceding CMEs and HSSs), this generally implies that modifications of the assumed half-width and/or propagation direction of the injected CME can improve the accuracy of MHD modelling and, accordingly, the SEP modelling results. This iterative process allows one to adjust and account for the errors inherited by the coronagraph fitting models, such as the GCS model, which depend solely on coronagraph images of the CME and are susceptible to projection effects.

• In the simulation, we assumed that the parallel mean free path is three orders of magnitude larger than the perpendicular mean free path, favouring particle propagation along the IMF lines. At first glance, this setup may seem to restrict the entry of SEPs into the closed magnetic topologies of preceding CMEs. However, as demonstrated in Sect. 5.2, the presence of these preceding CMEs has an additional effect of trapping the SEPs in the inner heliosphere. As a result, the parallel transport of particles is constrained, preventing the SEPs from escaping into the outer heliosphere. This confinement provides a longer duration for cross-field diffusion to occur, allowing particles to gradually fill the closed magnetic structures associated with the preceding CMEs.

This work has primarily focussed on determining the effects of preceding CMEs on the propagation of SEPs in interplanetary space by using EUHFORIA. However, it is important to note that many crucial processes involved in SEP generation occur closer to the Sun than the current capabilities of EUHFORIA allow us to model. To capture the complexity of the solar corona, other models such as COCONUT (Perri et al. 2022, 2023; Kuźma et al. 2023) or MAS (Mikić et al. 1999; Linker et al. 1999) have been developed. These coronal models can provide a more realistic representation of the early stages of SEP events by more accurately capturing the intricate physical processes occurring during a CME eruption.

Another significant challenge in heliospheric modelling is the limited availability of in situ data for model validation. While recent spacecraft missions such as Solar Orbiter (Müller et al. 2020), BepiColombo (Benkhoff et al. 2021), and Parker Solar Probe (Fox et al. 2016) aim at filling this observational gap, the number of observational points in the heliosphere remains limited. This scarcity of data hinders our ability to comprehensively assess the accuracy of heliospheric models. In this work we took a novel approach by using SEP modelling and observations as a means to evaluate the performance of heliospheric models on larger portions of the computational domain. As energetic particles travel larger distances from the Sun, they might interact with solar wind structures not directed to the instruments we have in space, which modifies their properties. Since the energetic particles mostly travel along the IMF lines, we can use the SEP observations not only to evaluate the performance of the SEP modelling but also to indirectly evaluate the accuracy of the 3D MHD models and their CME models on a more global scale. Every SEP event is unique. Starting by evaluating the onset time of the simulated event and the measured one is a good first approach, but accounting for different features in the intensity-time or anisotropy-time profiles is strongly case-dependent. In summary, the abrupt variations in intensity levels and directional characteristics of SEPs (as indicated by, for example, the anisotropy) frequently serve as indicators of the underlying solar wind conditions. Consequently, they harbour a wealth of (non-local) information that can and should be harnessed to enhance the accuracy of MHD simulations. The complexity of the particle behaviour and the solar wind configuration makes the retrieval of this valuable information from SEP observations difficult. To a certain extent, SEP models such as PARADISE facilitate this, as they enable a direct comparison between the SEP data and simulations. A simultaneous evaluation of the MHD and SEP simulations is important. These types of methods that include an iterative process, focussed on addressing the inaccuracies of both the MHD and SEP model, can lead to improvements and more globally accurate heliospheric modelling.

Movies

Movie 1 associated with Fig. 4 (Fig4_a) Access here

Movie 2 associated with Fig. 4 (Fig4_c) Access here

Movie 3 associated with Fig. 4 (Fig4_e) Access here

Movie 4 associated with Fig. 7 (Fig7_a_Case1_10_3) Access here

Movie 5 associated with Fig. 7 (Fig7_b_Case1_10_4) Access here

Movie 6 associated with Fig. 7 (Fig7_c_Case2_10_3) Access here

Movie 7 associated with Fig. 7 (Fig7_d_Case2_10_4) Access here

Movie 8 associated with Fig. 7 (Fig7_e_Case3_10_3) Access here

Movie 9 associated with Fig. 7 (Fig7_f_Case3_10_4) Access here

Acknowledgments

EUHFORIA was created as a joint effort between KU Leuven and the University of Helsinki and was developed further by the project EUHFORIA 2.0, a European Union’s Horizon 2020 research and innovation programs under grant agreement No. 870405 and the ESA project “Heliospheric modelling techniques” (Contract No. 4000133080/20/NL/CRS). These results were also obtained in the framework of the projects C14/19/089 (C1 project Internal Funds KU Leuven), G.0B58.23N and G.0025.23N (FWO-Vlaanderen), SIDC Data Exploitation (ESA Prodex-12), and BelSPo project B2/191/P1/SWiM. The ROB team thanks the Belgian Federal Science Policy Office (BELSPO) for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA) under contract numbers 4000134474. N.W. acknowledges support from NASA program NNH17ZDA001N-LWS and from the Research Foundation – Flanders (FWO-Vlaanderen, fellowship no. 1184319N). A.A. acknowledges the support by the Spanish Ministerio de Ciencia e Innovación (MICINN) under grant PID2019-105510GB-C31 and through the “Center of Excellence María de Maeztu 2020-2023” award to the ICCUB (CEX2019-000918-M) Computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Centre), funded by the FWO and the Flemish Government-Department EWI.

References

Appendix A: Effect of the preceding CME models on the shape of the SEP intensity-time profiles when λ_⊥/λ_∥ ∼ 10⁻⁴

In this section we build upon the discussion presented in Sect. 5.3 by examining the intensity and anisotropy-time profiles in the simulations where the perpendicular mean free path is reduced by one order of magnitude. Specifically, we present Fig. A.1, which displays the same information as Fig. 8, but for λ_⊥/λ_∥ ∼ 10⁻⁴ rather than λ_⊥/λ_∥ ∼ 10⁻³. As mentioned in Sect. 5.1.3, the simulated intensities are translated into physical units by scaling the 2.4 MeV intensity value (blue curve) to the observed value in the 2.2–2.7 MeV channel of SOHO/ERNE/LED at the time of the shock crossing by the Earth. As for the simulations with the larger λ_⊥, the scaling procedure is conducted individually for each case to ensure a fair comparison with the observations of the SEP event corresponding to each specific solar wind scenario. For Case 2, the ratio value used in comparison with the normalisation factor for Case 1 is 0.70, while for Case 3, the ratio value is 10.4.

Fig. A.1. PARADISE simulated proton time-intensity profiles (top panels) for energy channels ranging from 2.2 MeV to 51.1 MeV, with their corresponding first-order parallel anisotropy profiles (bottom panels). The top, middle, and bottom rows correspond to Case 1, Case 2, and Case 3 background solar wind for runs where λ⊥/λ∥ ∼ 10−4 at the different virtual spacecraft: E15 (left), Earth (centre), and W15 (right) columns. Dashed green, blue, and red lines correspond to the time of the Main CME-driven shock at the position of the virtual spacecraft for each of the cases. The vertical dash-dotted black line marks the time of injection of the Main CME at the inner boundary of the EUHFORIA model. The bottom red (blue) line indicates the periods of positive (negative) values of the radial component of the simulated magnetic field at the position of each spacecraft.

In general, the obtained characteristics of the intensity-time profiles shown in Figs. 8 and A.1 are similar. However, since particles are less prone to diffuse into closed magnetic structures, as in the latter case, the intensity-time profiles exhibit greater variability when such structures cross the observer’s path. Furthermore, all cases demonstrate higher levels of noise at the onset of the SEP event, especially when the magnetic connection with the main CME-driven shock is poor. This increased noise can be attributed to a reduced influx of particles at the observer’s location, resulting in inferior statistics. In the anisotropy-time profiles, this is more evident, as the features are very similar to those in Fig. 8 but with higher noise levels, making the features less clear.

The reduced cross-field diffusion in simulations with lower perpendicular mean free path values highlights that the incorporation of the Second CME, regardless of the CME model employed, results in earlier onset times for SEPs. As previously discussed in Sect. 5.3, this phenomenon arises due to the improved magnetic connectivity between the observer and the Main CME in the cases when the Second CME is present. Consequently, in Cases 2 and 3, the arrival of the particles occurs earlier than in Case 1.

The effect of the First CME is more evident in Case 1 as it somewhat blocks the arrival of the particles to the three observers, at the beginning of the event. Case 2 shows smoother profiles, which evidences an improved magnetic connection with respect to Case 1. Case 3 shows the more pronounced variation in the profiles due to the more intricate IMF structure resulting from the three LFF-spheromak CMEs. In particular, the W15 observer is the most affected by the two preceding CMEs. This observer shows more variations in the intensity-time profiles than for the E15 and Earth observer because of the more direct impact of the Second CME.

For an E15 observer, the effect of the worse initial magnetic connection to the Main CME-driven shock is more noticeable, as it makes particles arrive even later than in Fig. 8. Due to the lower λ_⊥, particles are not able to diffuse into regions with a better connection to this observer and reach the observer. Regardless of the CME model, the inclusion of the Second CME modifies the IMF so that the E15 observer establishes an improved connection with the Main shock.

As discussed in Sect. 5.3, the W15 observed is located within the magnetic structure of the First CME at the beginning of the particle event. The lower intensities shown during this time when comparing with the corresponding cases in Fig. 8 indicate that most of the particles reach the W15 observer through perpendicular diffusion.

For the Earth observers in Fig. A.1, the smaller perpendicular mean free path translates into lower proton intensity-values in the initial phase of the profiles when compared to Fig. 8, and more pronounced variations in the profiles responding to the presence of the preceding CMEs. For instance, as particles are not able to easily diffuse into the First CME, they bounce back towards the Earth. This can be seen in Case 1 as a second increase at approximately 16:00 UT on March 16, and in Case 2, with rather flat profiles (even for the highest energy channel) up to the shock passage. Finally, in Case 3, particle intensities react to the passage of the Second CME by exhibiting a decrease at high-energies and a plateau at the lower-energies. Later, the intensities show a larger increase towards the shock arrival than in Fig. 8 due to SEPs bouncing between the Main and the Second CME.

All Tables

All Figures

Fig. 1. In situ data at L1. Top panel: intensity time profiles for solar energetic protons observed by ERNE aboard SOHO for energy channels ranging between 2.2 MeV and 68.0 MeV (labels correspond to the average energy of the channels). Lower panels, from top to bottom: solar wind speed, density, temperature, and the different components of the magnetic field from the Wind spacecraft. The vertical dotted red line shows the launch time of the Main CME from the Sun. Grey and green shaded areas correspond to the duration of the First and Main CMEs’ effects at Earth, respectively.

In the text

Fig. 2. SDO/AIA images at different wavelengths. (a) Image at 19.3 nm for 2013 March 12 at 12:15 UT. The coronal hole situated close to the northern polar coronal hole is labelled CH. The coronal hole has a positive polarity, which is opposite to the expected polarity of the nearby northern polar coronal hole. During the year 2013, the polarity inversion was still in progress, so the northern polar coronal hole was very small and hard to detect (Janardhan et al. 2018). (b) Image at 19.3 nm for 2013 March 15 at 07:45 UT. The inverted sigmoid position is highlighted with the yellow box. (c) Image at 9.4 nm for 2013 March 15 at 06:03 UT. The inverted sigmoid position is again highlighted with a yellow square.

In the text

	Fig. 3. White light coronagraph images. Top panels correspond to the GCS fitting (green mesh) for the Main CME as seen by (a) LASCO-C2, (b) STEREO-COR2A, and (c) STEREO-COR2B. The bottom panels correspond to the images of (d) LASCO-C2, (e) STEREO-COR2A, and (f) STEREO-COR2B, used for the fitting procedure.
In the text

Fig. 4. Snapshot of the EUHFORIA simulations for Case 1 (a and b), Case 2 (c and d), and Case 3 (e and f) for 2013 March 16 at 06:44 UT. The top row shows the radial solar wind speed, and the bottom row shows the scaled number density in the equatorial plane. White, red, and black arrows point at the approximate position of the First, Second, and Main CMEs (respectively). Dashed lines correspond to the magnetic field line connection from the E15, Earth, and W15 positions, backtracked to 0.1 au. The temporal evolution for (a), (c), and (e) are available as an online movie.

In the text

Fig. 5. EUHFORIA simulations compared with in situ data. Top to bottom: solar wind speed, density, temperature, total magnetic field, and its three components. In all panels, the observations from Wind spacecraft (black dots) and EUHFORIA simulation results are given for Case 1 (green), Case 2 (blue), and Case 3 (red). The vertical dotted red line shows the launch time of the Main CME event from the Sun. Grey and green shaded areas correspond to the duration of the First and Main CME’s impact at Earth.

In the text

	Fig. 6. EUHFORIA background solar wind snapshot with shock tracer function applied at the time instant 10:14:22 UT. White, red, and black arrows point at the approximate position of the First, Second, and Main CME (respectively).
In the text

Fig. 7. PARADISE simulation snapshot for 2013 March 16, 06:44 UT. Panels a and b correspond to Case 1, panels c and d correspond to Case 2, and panels e and f correspond to Case 3 background solar wind. The left column shows the simulations using λ⊥/λ∥ ∼ 10−3, while the right column corresponds to λ⊥/λ∥ ∼ 10−4. E15, Earth, and W15 virtual spacecraft are represented with a green square, a red circle, and a green triangle (respectively). White, red, and black arrows point at the approximate position of the First, Second, and Main CME, respectively. The temporal evolution of each case is available as an online movie.

In the text

Fig. 8. PARADISE simulated proton time-intensity profiles (top panels) for energy channels ranging from 2.2 MeV to 51.1 MeV, with their corresponding first-order parallel anisotropy profiles (bottom panels). The top, middle, and bottom rows correspond to Cases 1–3 background solar wind for runs where λ⊥/λ∥ ∼ 10−3 at the different virtual spacecraft: E15 (left), Earth (centre), and W15 (right) columns. Dashed green, blue, and red lines correspond to the time of the Main CME-driven shock at the position of the virtual spacecraft for each of the cases. The vertical dash-dotted black line marks the time of injection of the Main CME at the inner boundary of the EUHFORIA model. The bottom red (blue) line indicates the periods of positive (negative) values of the radial component of the simulated magnetic field at the position of each spacecraft.

In the text

Fig. A.1. PARADISE simulated proton time-intensity profiles (top panels) for energy channels ranging from 2.2 MeV to 51.1 MeV, with their corresponding first-order parallel anisotropy profiles (bottom panels). The top, middle, and bottom rows correspond to Case 1, Case 2, and Case 3 background solar wind for runs where λ⊥/λ∥ ∼ 10−4 at the different virtual spacecraft: E15 (left), Earth (centre), and W15 (right) columns. Dashed green, blue, and red lines correspond to the time of the Main CME-driven shock at the position of the virtual spacecraft for each of the cases. The vertical dash-dotted black line marks the time of injection of the Main CME at the inner boundary of the EUHFORIA model. The bottom red (blue) line indicates the periods of positive (negative) values of the radial component of the simulated magnetic field at the position of each spacecraft.

In the text