© J. Pomoell et al., Published by EDP Sciences 2015

1. Introduction

Large solar energetic particle (SEP) events result mainly from the acceleration of particles at shocks driven by coronal mass ejections (CMEs), propagating from the low solar corona out into interplanetary (IP) space (e.g., Cliver et al. 1982; Sanahuja et al. 1983; Reames 1999; Kahler 2001). The accelerated particles travel along the interplanetary magnetic field (IMF) and are finally detected by space-based observatories as an SEP event. The plasma conditions at the IP shock front are known to be fundamental in determining the resulting characteristics of the energetic particle population (e.g., Bell 1978; Lee 1983a, 1983b, 2005; Heras et al. 1995; Zank et al. 2000; Vainio & Laitinen 2007). In order to be able to test current models – a crucial step for constructing tools for predicting fluences and peak fluxes of SEP events – the plasma conditions at the propagating shock front need to be determined (e.g., Cane & Lario 2006).

A fundamental difficulty in such an endeavour is the fact that systematic observations of the plasma and energetic particles at different radial distances away from the Sun are scarce, and will continue to be so for the foreseeable future. One way of circumventing this problem, as adopted here (see also Aran et al. 2007; Luhmann et al. 2010; Rodrìguez-Gasén et al. 2014, and references therein), is by performing magnetohydrodynamic (MHD) simulations of the IP propagation of a CME-driven shock in order to obtain the plasma parameters at the shock front, and feeding the acquired shock parameters into a simulation describing the transport of particles from the shock to the observer.

When modelling SEP events using simulations of IP shock propagation and particle transport models, a key issue is the identification of the shock front and the extraction of the plasma parameters across it. Various methodologies and different plasma variables have been used by several authors aiming at achieving a proper determination of the strength of the shock (e.g., Lario et al. 1998; Manchester et al. 2005; Luhmann et al. 2010; Rodríguez-Gasén et al. 2011; Rouillard et al. 2011). Unfortunately, even in the absence of shocks, the IP medium is never uniform and the location where to evaluate the values of the plasma variables and compute the jump conditions is not an easy task. This is because usually, contrary to what is depicted in theoretical cartoons, there is no clear plateau (or it is very small) indicating the downstream (post-shock) region from which the peak value can be identified. As Luhmann et al. (2010) stated: “(…) In any case shock identification in the MHD codes remains a key challenge for SEP event modelling.”. It is the goal of the current work to address this problem.

The compound “Shock-and-Particle” model (SaP; Lario et al. 1998) combines an MHD simulation of the propagation of a CME-driven shock and a simulation of the transport of particles along IMF lines via the concept of the cobpoint (connecting with the OBsever POINT, Heras et al. 1995). The cobpoint is the point at the shock front that is magnetically connected to the observer. It is determined from the simulation of the shock propagation which utilizes observational data to constrain the free parameters of the model in order to reproduce the IP transit time of the shock as well as the observed plasma jumps at the shock as it crosses the observer’s position. On the other hand, the particle transport model uses SEP data to derive the evolution of the source function of shock-accelerated protons, Q, at the cobpoint. In order to fit the observed particle intensity-time profiles (and the first-order parallel anisotropies) of a gradual SEP event, we assume that the detected particles have been injected at the cobpoint and have propagated along the magnetic flux tube connected to the observer. As the shock propagates and expands, the cobpoint slides eastward along the shock front and, hence, the observer is connected to different regions of the shock front. In our simulation approach we consider a unique magnetic flux tube corotating with the Sun. Thus, the protons successively injected at the cobpoint populate the same magnetic flux tube¹. The MHD simulation of the propagation of the IP shock gives the position of the shock front and the location of the cobpoint in it. The solar wind and magnetic field are assumed to be stable in the upstream region of the shock and the description of the particle transport is limited to this upstream region. Under such conditions, the large-scale upstream IMF can be described by a Parker spiral.

By modelling several events, a semi-empirical relation between Q and the normalized radial velocity jump at the cobpoint, VR, was derived (see e.g., Lario et al. 1998; Aran 2007). Once a functional dependence, Q = Q(VR), is established, it is then possible to invert the procedure and produce synthetic SEP flux profiles for given CME-driven shocks (Aran et al. 2007; Rodríguez-Gasén et al. 2011). In this way the SaP model allowed a space weather related tool to be developed with the capability to provide quantitative predictions of <90 MeV proton flux and cumulative fluence profiles in gradual SEP events, known as SOLPENCO (SOLar Particle Engineering Code, Aran et al. 2006; Aran 2007; Watermann et al. 2009). Recently, in the framework of the European Space Agency’s Solar Energetic Particle Environment Modelling (SEPEM) project², an upgrade of the SaP model was developed (Aran et al. 2011; Jacobs & Poedts 2011). This version incorporates a two-dimensional (2D) MHD model of shock propagation that enables the application of the Q(VR) relation to the prediction of proton intensity-time profiles up to 200 MeV. Based on this model, we developed the SOLPENCO2 tool (Crosby et al. 2015). This tool provides predictions of 5–200 MeV proton peak intensities and event fluences of SEP events for virtual observers located at heliocentric radial distances between 0.2 AU and 1.6 AU. The derived radial dependencies are used as input to the SEPEM statistical modelling tools for interplanetary missions³. Results from the SEPEM project are presented in Crosby et al. (2015).

In this work, we describe the newly developed algorithms by which the necessary plasma characteristics are extracted from the MHD shock simulation. An outstanding feature of the tool is that the shock parameters are obtained in an automated fashion. This is necessary not only in order to ensure objectiveness in the determination of the plasma conditions, but also necessary from a practical viewpoint, as a large number of shock crossings need to be processed for each simulation. Thus, the tool allows the simulation of a virtual armada of spacecraft scanning at a distance the front of the expanding IP shock and observing the plasma features at different locations, as if they were real spaceborne observatories. Moreover, we apply the new algorithms to two SEP events simulated with a new solar wind model, and compare the results with those obtained by Aran et al. (2011). The aim is to investigate the robustness of the Q(VR) relation, relying on the new solar wind model and the methodology used to determine the parameters of the shock.

This article is structured as follows: an overview of the MHD solar wind and the CME-driven shock model is described in Section 2. The process employed by the shock extraction tool is presented in Section 3. In Section 4, we show the application of these models and of the shock extraction tool to the simulation of two large gradual SEP events and we discuss the correlation found between Q and VR after applying the particle transport model to reproduce the observed proton intensity-time profiles. Finally, a summary and conclusions are provided in Section 5.

2. Sun-to-Earth MHD modelling of the solar wind and CME-driven shock

Similar to our previous work (i.e. Heras et al. 1995; Lario et al. 1998; Aran et al. 2007, 2011; Rodríguez-Gasén 2011), we employ 2D magnetohydrodynamic simulations in order to model the evolution of CME-driven shocks from the Sun to the Earth in the ecliptic plane. The equations used are those of ideal MHD augmented with momentum and energy sources describing gravity and a parametrized coronal heating term, and can be written (in SI units) as the set of equations:(1)

(2)

(3)

(4)giving the time dependence of the mass density (ρ), velocity (v), magnetic field (B) and total energy density (E). Furthermore, I is the unit dyad, is the gravitational acceleration, while S is an energy source term (specified in Sect. 2.1), and the total energy density is related to the thermal pressure (P) by(5)where γ is the adiabatic index which is chosen to be constant. Additionally, an ideal gas law is invoked for computing the temperature,(6)where m is the mean molecular mass taken to be where m_p is the proton mass.

2.1. Solar wind model

In order to accurately simulate the evolution of a CME-driven shock propagating from the corona to 1 AU and beyond, a model describing the ambient interplanetary plasma characterized by the solar wind in which the shock expands is needed. In application-oriented modelling where the focus is not on studying the coronal heating problem, the fundamental problem of accelerating the solar wind is typically circumvented by employing a semi-empirical approach in which the required coronal heating is added as parametrized sources into the momentum and/or energy equations.

However, in doing so, care must taken so that the addition of momentum and/or energy does not affect the dynamics of the simulation in an unrealistic way. For instance, utilizing a much lowered polytropic index or a temperature-dependent cooling/heating energy source term has been shown to interfere detrimentally with the dynamics of shocks in particular (Pomoell & Vainio 2012). Therefore, in the present version of the SaP model, we have implemented a new solar wind model in order to obtain more realistic compression and heating of the plasma. Specifically, we adopt a simple coronal energy source term(7)allowing us to choose the adiabatic index to be constant and closer to the value of γ = 5/3 expected for a monoatomic plasma.

The background solar wind model utilized in this work, similar to the solar wind model by Jacobs & Poedts (2011), is time-independent and 1.5-dimensional (i.e. the MHD variables are functions of radial distance only, but vector quantities have two components). The solution is obtained by running a time-dependent computation by solving Eqs. (1)–(4) until a steady-state solution is obtained. In order to do so, the density ρ₀ ≡ ρ(r = r₀) and temperature T₀ ≡ T(r = r₀) in the low corona at r₀ = 1.03 R_⨀ as well as the heating parameters S₀ and L and the adiabatic index γ need to be specified. These are to be chosen so that the resulting steady-state (time-independent) azimuthally-symmetric solar wind solution in the ecliptic plane gives a realistic radial dependence of the plasma while simultaneously producing values representative of the solar wind conditions at 1 AU during the time from the onset of the considered event in the corona to the time of arrival of the shock at the spacecraft.

In principle, determining the free parameters ρ₀, T₀, S₀, L, γ in order to reproduce the observed solar wind values at 1 AU prior to the shock arrival could be done by running a large number of time-consuming MHD relaxation simulations. The process can, however, be made significantly faster by considering a simplified solar wind model that can be solved essentially instantaneously.

The basic realization allowing to simplify the solar wind model is the well-known result (e.g., Weber & Davis 1967) that the density and radial velocity profiles are altered only slightly when solar rotation and magnetic field (using values typical for the Sun) are taken into account. Thus, the assumption of a non-magnetic non-rotating solar wind is a good approximation of the full system. Moreover, in this approximation, the steady-state set of equations can be solved directly without the need of a costly time-dependent relaxation simulation (de Sterck et al. 2009).

With the 1.5D MHD solution obtained, the initial steady-state solar wind solution for the entire ecliptic plane for use in the shock simulation is constructed by simply replicating the solution in the azimuthal direction. Note that, as a consequence, a single solar wind state is present in the entire ecliptic plane and, e.g., no stream interfaces between winds of differing speeds are present. Therefore, the current wind model is suited for time intervals when the solar wind is steady and does not change significantly.

2.2. CME-driven shock model

The shock simulation is performed by superimposing a perturbation on the background solar wind, following a similar approach to that adopted by, e.g., Jacobs et al. (2005, 2007); i.e. in the same way as described in Rodríguez-Gasén (2011) and Aran et al. (2011). The modelled perturbation consists of an enhancement in both solar wind density and radial velocity. The reason for using such a simple model is that we are interested in the propagation of the IP shock rather than in the detailed initiation process of the CME. Moreover, the restriction of considering only dynamics in the ecliptic plane places limitations on the types of CME models that are reasonable to pursue.

The initial disturbance is described by seven parameters, which may be varied in order to obtain an agreement between the estimated transit time as well as the in-situ data as the shock crossed the spacecraft. The profiles of the initial perturbation, both in density (ρ_extra) and radial velocity (v_extra), have the same functional form as Rodríguez-Gasén (2011) and Aran et al. (2011). That is:(8)

(9)where(10)

(11)with

The initial disturbance is located at r_cme = 1.75 R_⨀ and centred around φ = 180°. The disturbance only affects the region and The parameters a_cme and d_cme control the angular extent and thickness of the initial perturbation, respectively, and are fixed to a_cme = 160° and d_cme = 0.5 R_⨀.

The parameter A controls the shape of the exponential pulse profile in the radial direction, and it is set to a fixed value, The parameter ΔΦ defines the distribution of the density and velocity in the angular direction and may vary from event to event. We have opted for this profile of the initial disturbance (i.e. Eq. (8)), instead of utilizing a simple spherical blob, in order to have more control over the angular extent of the propagating shock and to keep the analogy with the previous simulations of shock propagation used in SOLPENCO (Aran et al. 2006; Aran 2007). By adjusting the value of ΔΦ the shape of the shock front can be varied from parabolic-like (ΔΦ = 0.5 a_cme) to circular-like (ΔΦ = 0). Thereby, the model can be used for narrow as well as for wide events allowing, for instance, the simulation of 180°-wide events as were analysed earlier using multi-spacecraft observations (Zurbuchen & Richardson 2006).

Table 1 lists the seven parameters that describe the initial perturbation. Four out of the seven parameters are fixed to a given value, as commented previously and indicated in this table. The parameters indicated by “free” are determined differently for each event considered (see Sect. 4.1).

The simulation of the shock propagation is performed on a 2D polar grid, covering the equatorial plane. In this 2D approach we identify the ecliptic with the equatorial plane. The nominal grid consists of 868 cells in the radial direction and 404 cells in the longitudinal direction, covering 360° with two ghost cells at each boundary. The number of cells in the grid can be, however, easily changed. In the radial direction the grid extends from 1.2 R_⨀ up to 215 R_⨀ (= 1.0 AU). This radial extent is flexible in the sense that it could be enlarged to accommodate the simulation of SEP events seen by STEREO-B or up to the Mars’ orbit, if required. The radial size of cells is non-uniform, varying from Δr = 0.003 R_⨀ at the inner boundary to Δr = 1.09 R_⨀ at 215.3 R_⨀. For 1.2 < r < 4 R_⨀, the cell size increases by 0.50% from one cell to the next; for 4 < r < 18 R_⨀, such an increase is of the order of 1.5% and for r > 18 R_⨀ it is of the order of 0.44%. In the longitudinal direction, the grid is equidistant, with Δφ = 0.9°.

The MHD equations are solved using a conservative second-order finite volume method in which Eqs. (1), (2) and (4) are solved using the methods detailed in Kissmann et al. (2009), while the induction equation is solved by employing a consistent constrained transport scheme (Kissmann & Pomoell 2012) so that the divergence of the magnetic field remains zero to machine precision throughout the computation.

3. Extracting the shock data

In our previous SaP modelling, we developed an automated method to extract the parameters of the MHD shock (Rodríguez-Gasén 2011). In the present work we employ a newly developed alternative algorithm instead. The motivation for developing a second method was to have two independent computations of the shock parameters thereby allowing us to estimate the errors involved in the extraction process itself. A goal was also to find a method using a reduced number of free parameters as compared with the previous one. In Section 4, we employ both methods and discuss their differences.

The task of the extraction tool is essentially to locate the point on the shock front that is magnetically connected to the observer, the cobpoint, and thereafter determine the plasma parameters of the shock at the cobpoint. The flow of the process is as follows:

1. Determine locations of the virtual spacecraft.

2. For each MHD snapshot in time and each virtual spacecraft:

   1. Compute the magnetic field line passing through the virtual spacecraft.

   2. Determine the location of the cobpoint, i.e. the location where the field line computed in the previous step crosses the shock.

   3. Compute the shock normal at the cobpoint.

   4. Determine the locations of the regions upstream and downstream of the shock.

   5. Compute the properties of the shock at the cobpoint, e.g., density compression ratio, speed of shock, the normalized downstream-to-upstream radial velocity ratio, etc.

   6. The methods we employ to perform the individual tasks in step 2 are detailed in Appendix A.

3. For each virtual spacecraft, output the characteristic parameters of the shock as a function of time to a file suitable for use by the particle transport code.

4. Application to two large SEP events

We have applied the new MHD shock-and-particle model to simulate two large SEP events, namely the 2000 April 4 (Apr00) and 2006 December 13 (Dec06) events. These two particular events were chosen in order to compare the results with our previous simulations of the events, thereby allowing us to study the influence of the new codes on the Q(VR) relations on which the SEPEM/SOLPENCO2 tool is based (Aran et al. 2011; Crosby et al. 2015).

4.1. Overviews of the events

The 2000 April 4 event occurred during the maximum phase of Solar Cycle 23. This SEP event is associated with a strong IP shock detected at ACE at 16:00 UT on April 6, driven by a full halo CME observed by SOHO/LASCO at 16:32 UT on April 4. The CME was in temporal association with a C9.7/2F flare observed at 15:12 UT on April 4, located in the AR 9833 region (N16 W66). The transit time of the shock, measured as the time elapsed from the start time of the X-ray flare until the shock passage by the ACE spacecraft, is 48.82 h, with an average speed of 843 km s⁻¹ (Aran et al. 2011). For a detailed analysis of the solar origin and the interplanetary conditions during this SEP event we refer to Rodriguez et al. (2009). The proton differential intensity enhancements extend only up to 80 MeV protons as measured by the SOHO/ERNE experiment (Torsti et al. 1995). On the other hand, the 2006 December 13 SEP event happened close to the minimum phase of Solar Cycle 23 and it was the last event of the cycle to produce a Ground Level Event (GLE; e.g., Storini et al. 2008). This event was associated with an X3.4/4B flare starting at 02:14 UT on December 13 in AR 10,930 (S06 W23) and a Halo CME detected by LASCO C2 at 02:54 UT on the same day with a plane of sky speed of 1774 km s⁻¹. An interplanetary shock associated with the interplanetary counterpart of the halo CME was detected by ACE at 13:52 UT on December 14. Hence it had a transit time of 35.64 h and an average speed of 1155 km s⁻¹. Details on the solar origin of this event and a description of the measured intensities can be, respectively, found in e.g., Liu et al. (2008) and Verkhoglyadova et al. (2010).

4.2. Modelling the solar wind and the IP shock

In order to model the shock propagation we first reproduce the solar wind conditions prior to the shock arrival at 1 AU. Figures 1 and 2 display the solar wind density, radial velocity and temperature as measured by the ACE/SWEPAM instrument (McComas et al. 1998) as well as the magnetic field strength measured by the ACE/MAG instrument (Smith et al. 1998), all plotted in red, for the Apr00 and Dec06 events, respectively. In these plots, time is measured from the shock passage detected by ACE through the simultaneously measured jumps in these plasma variables. Therefore, negative times represent pre-event solar wind conditions.

Fig. 1. 2000 April 4 event. From left to right, top to bottom: the solar wind density and radial velocity, and the magnetic field strength and the solar wind temperature measured by the ACE spacecraft (red curves). Corresponding values from the simulation are plotted in black. Time is counted from the time of the interplanetary shock passage by the ACE spacecraft.

Fig. 2. 2006 December 13 event. From left to right, top to bottom: the solar wind density and radial velocity, and the magnetic field strength and the solar wind temperature measured by the ACE spacecraft (red curves). Corresponding values from the simulation are plotted in black. Time is counted from the time of the interplanetary shock passage by the ACE spacecraft.

As can be seen in these figures, the simulated solar wind (black curves) reproduces the large-scale values of the pre-event solar wind conditions for both events. The time-independent nature of the solar wind model used in the MHD computation is clearly visible until the arrival of the shock at the observer. The values of the initial parameters used to obtain the two simulated solar wind conditions are listed in Table 2, while Figure 3 shows the radial distance profiles of the plasma variables in the case of the Apr00 (left) and Dec06 (right) events.

Fig. 3. The plasma variables as a function of radial distance from the Sun preceding the 2000 April 4 (left) and 2006 December 13 (right) CME-driven shocks.

Note that the Apr00 event developed under slow solar wind conditions (v_sw ~ 380 km s⁻¹) whereas the Dec06 event occurred under fast solar wind conditions (v_sw ~ 595 km s⁻¹). As can be seen in Figure 3, for both events the solar wind was accelerated within 20 R_⨀ from the Sun with a temperature profile already decreasing at ≤4 R_⨀ in contrast with the previous polytropic solar wind model (>8 R_⨀) (Jacobs & Poedts 2011).

With the 1.5D solar wind modelled, we next proceed to simulate the IP shock propagation using the disturbance model described in Section 2.2. When modelling the IP shock associated with an SEP event the values of the free parameters defining the initial shock characteristics (ρ_cme, v_cme and ΔΦ) are obtained following a process of trial and error. The observables used in order to restrict the possible values of these parameters are: the transit time of the shock and jumps in the plasma variables, the radial distance of the spacecraft from the Sun and the longitudinal separation of the spacecraft from the leading edge of the shock or CME. If this latter information is not available, then the longitude of the accompanying X-ray (or H_α in lieu of that) flare is used.

First, we keep the value of ΔΦ fixed, and we produce different test runs by varying ρ_cme and v_cme in order to approximately match the observed transit time of the shock. Under these conditions, the transit times of the disturbances are ordered by increasing values of the quantity (e.g., Aran et al. 2011). In this step we have found it beneficial to first perform 1.5D simulations of the shock propagation in order to limit the parameter space to be explored using 2D computations. Finally, we vary the obtained values of ρ_cme and v_cme and that of ΔΦ, in order to reproduce as closely as possible the observed plasma jumps, but putting an emphasis on capturing the jumps in radial velocity and number density. In this step it is important to know the relative position of the observer with respect to the incoming shock because the initial parameters may vary depending on whether the observer is crossed by a flank of the shock or by the central region of the shock front. Table 3 lists the values of the parameters describing the initial shock conditions for both events.

Figures 1 and 2 illustrate the level of accuracy that we aim at when reproducing the observed plasma jumps and the time of shock arrival. These figures show the resulting simulation profiles (solid black lines) matching the plasma data (red dots) measured aboard the ACE spacecraft. As can be seen, the observed jumps are reproduced for all variables in the case of the Dec06 event and for the solar wind density and radial velocity in the case of the Apr00 event. In the latter case the model fails to reproduce the observed jumps in magnetic field and temperature. This is mainly due to the assumption of a merely hydrodynamic initial shock and to the 2D geometry adopted. Note that this level of accuracy in matching the time of shock passage at the spacecraft (i.e. the observer) and the jumps is needed in order to derive as correctly as possible the values of the plasma and magnetic field jumps across the shock. This is still rather unaffordable due to the amount of computational resources needed for a 3D MHD model. Most of the currently available 3D MHD models applied to real SEP events up to 1 AU use a coarser resolution (e.g., Luhmann et al. 2010; Shen et al. 2011). For instance, Figure 10 of Shen et al. (2011) shows a 3D simulation of the shock associated with the 2000 April 4 event, for which the simulated plasma jumps at the shock are not captured with enough resolution to be used in conjunction with the proton transport model, since we are seeking a relation between the source of particles at the shock and the plasma jumps.

In order to better constrain the values of the three free parameters of the disturbance model and to verify the adopted model, it would be desirable to perform simulations of multi-spacecraft SEP events such as that performed by Rodríguez-Gasén (2011). Nevertheless, we have started by modelling two single-point-measured SEP events, Apr00 and Dec06, in order to control the changes introduced by the new solar wind model described in Section 2 with respect to the simulations of these two SEP events performed with the previous version of the SaP model (Aran et al. 2011; Jacobs & Poedts 2011).

To illustrate the simulations, Figure 4 shows the radial velocity contours at 20.5 h for the Dec06 event. The observer in the middle corresponds to the near-Earth spacecraft that measured this event. The interplanetary magnetic field line shows that at that time the observer was connected to the central region of the propagating shock. The point of intersection is the “cobpoint”. From the simulation of the shock propagation we obtain the position of this cobpoint from about 3.5 R_⨀ up to the observer’s location. This information together with the shock average velocity is stored every 10 min. These outputs of the MHD simulation are used as inputs to the particle transport code (Lario et al. 1998).

Fig. 4. Radial velocity at t = 20.5 h after the CME onset for the 2006 December 13 event. Five observers (dots) with their corresponding magnetic field lines connecting with the shock front are displayed in grey. The observer, corresponding to the near-Earth spacecraft, is the one in the middle. Their corresponding cobpoints are the circles on the shock front. Black lines across them indicate the shock normal whereas magenta lines show the radial direction.

4.3. Evolution of the plasma parameters at the cobpoint

As in the works of Heras et al. (1995) and Lario et al. (1998), along with the position of the cobpoint, the extraction procedure also outputs the values of the following ratios computed along the radial direction through the cobpoint:(12)where “d” and “u” stand for the downstream and upstream points, respectively. Note that the velocity components are measured in a frame of reference at rest, centred at the Sun and not in the shock reference frame. The previous (SEPEM) as well as the current procedures to determine the cobpoint and the plasma characteristics at the shock also allow the extraction of these variables along the shock normal direction as well as the retrieval of the angle between the shock normal and the upstream magnetic field direction, i.e. the θ_Bn parameter. The motivation for using the parameters in Eq. (12) stems from the experience obtained since the first SaP models (i.e. Heras et al. 1995; Lario et al. 1998).

In particular, as a result of the modelling of a number of SEP events, a semi-empirical relation was established between the injection rate of shock-accelerated particles, Q, and the VR ratio, both obtained at the cobpoint position (Lario et al. 1998; Aran 2007; Aran et al. 2007):(13)

This association, termed as the “Q(VR) relation” (Lario et al. 1998), permits the calculation of proton synthetic intensity-time profiles at a given location in the ecliptic plane, once VR(t) is known for the virtual observer of interest. Note that this Q(VR) relation is only applicable while there is a magnetic connection between the observer and the shock front. It is the base of the SOLPENCO and SOLPENCO2 tools. These tools provide predictions of SEP intensity-time profiles for different virtual observers and under various IP shocks and transport conditions (see details in Aran et al. 2006, 2008; Watermann et al. 2009; Aran et al. 2011).

Figure 5 shows the comparison of the position of the cobpoint and the evolution of VR, BR and r, and of the transit speed of the shock for the Apr00 (left) and Dec06 (right) SEP events using two different models. The red curves are the results obtained from the MHD simulations performed with the previous polytropic solar wind model (Aran et al. 2011; Jacobs & Poedts 2011) and applying the SEPEM cobpoint algorithm described in p. 132 of Rodríguez-Gasén (2011). The black curves are obtained with the solar wind model and the cobpoint procedure described in Sections 2 and 3, respectively. Furthermore, blue curves depict results obtained with the new solar wind model but using a slightly modified version of the shock parameter extraction tool employed in SEPEM. In the latter, we use the relative entropy in addition to the relative density and velocity to determine the cobpoint candidate (Horne et al. 2013).

Fig. 5. The cobpoint position and plasma jumps at the cobpoint for the 2000 April 4 (left) and 2006 December 13 (right) SEP events. From top to bottom: The first and second panels show the radial position and the angular position of the cobpoint with respect to the launch direction of the shock. The third, fourth and fifth panels, respectively, show the normalized downstream-to-upstream radial velocity ratio, VR, the jump in magnetic field strength, BR, and the density compression ratio, r, respectively. The sixth panel displays the Mach number, M; the dashed horizontal line marks M = 1. The bottom panel shows the transit speed of the shock up to each cobpoint position. The black curves correspond to the results applying the methodology described in Section 3. Blue curves display the results obtained with the cobpoint procedure used in SEPEM but with the present MHD simulation model, and the red curves show the same quantities as derived from the model developed during the SEPEM project.

The two methods for extracting the shock data (black and blue curves in Fig. 5) yield almost the same values in the case of the Apr00 event. The SEPEM method (blue) presents larger fluctuations when determining the plasma parameters than the method presented in this work (Sect. 3). It is worthwhile to note that the fluctuations occur due to inherent difficulties in extracting the shock characteristics, and are not a signature of oscillatory behaviour in the MHD simulation. In particular, at the end of the event, the fluctuations in the blue curve are more frequent because the cobpoint algorithm cannot detect the downstream point and then the code outputs a zero value. In the case of the Dec06 SEP event, the differences between both algorithms are larger, and particularly from ~7 to ~13 h. These differences result mainly from the different methods used in determining the downstream point (i.e. the start of the downstream region as described in Appendix A). Note that the new methodology uses the normal direction through the cobpoint to determine these parameters whereas the SEPEM method uses the radial direction. For the case of VR, the methods differ up to ΔVR = 0.5 in this interval, while for the remainder of the event we have ΔVR < 0.1. These values serve as a valuable indication of the uncertainty involved in calculating VR.

Comparing the results obtained by using the two different solar wind models (red and black/blue curves in Fig. 5), we observe that the differences are larger in the case of the Apr00 event than for the Dec06 event. In the case of the former simulation of the Apr00 event, in order to match the time of the shock arrival and the plasma jumps, the values of the parameters describing the initial disturbance were larger than in the present simulation (a factor of 3.76 for ρ_cme and a factor of 1.98 for v_cme, with the shock being narrower since ΔΦ = 0.25 a_cme). This means that the energy input was larger than that used in the present simulation. On the other hand, the use of low values of the polytropic index for r < 50 R_⨀ (Aran et al. 2011; Jacobs & Poedts 2011) translates into the kinetic energy of the shock being used to compress the plasma rather than heating it (Pomoell & Vainio 2012). These two facts, the larger kinetic energy in the initial input disturbance and the low polytropic index, yield the large compression ratios, r and BR, and in turn the large values of the Mach number (red curves), shown in Figure 5. The values for these parameters, obtained with the present model (black/blue curves), are more in accordance with theoretical jump predictions for MHD shocks. The same conclusions apply also in the case of the simulation of the IP shock in the Dec06 event. In this case, the peak density of the initial disturbance superimposed on top of the polytropic wind (Aran et al. 2011), ρ_cme, was a factor of 0.53 smaller than in the present simulation (see Table 3) whereas v_cme was only 1.47 times larger (ΔΦ had the same values).

Despite the large differences in the evolution of BR and r, the overall evolution of VR is similar in the two solar wind models. However, the new solar wind yields, at the beginning of the event, smaller values due to the less compressed plasma (Pomoell & Vainio 2012). The large values of VR at the beginning of the event (the two-three first points, r < 6 R_⨀) are mainly due to two facts: (1) the solar wind is still accelerating, and hence v _r (u) is small, and (2) the fact that at these distances the shock has not yet detached from the actual compressed material forming the CME ejecta, a phenomenon found in other simulation studies as well (Manchester et al. 2008; Luhmann et al. 2010; Rodríguez-Gasén et al. 2011).

4.4. The injection rate of escaping protons

The description of SEP events associated with travelling shocks requires consideration of a moving source of particles (e.g., Heras et al. 1994; Cane 1995; Reames et al. 1996). We use the model by Lario et al. (1998) to solve the transport of protons from the cobpoint to the observer. The focused-diffusion transport equation used (Ruffolo 1995) accounts for the effects of particle focusing in the diverging IMF, pitch-angle scattering by the IMF irregularities, solar wind convection and adiabatic deceleration. A varying source term is added to the model to describe the injection of particles at the shock front (Lario et al. 1998).

The transport equation we use is Eq. (1) in Lario et al. (1998). In this equation, F(t, μ, r, p) is the distribution function of the particles in the simulated magnetic flux tube, and it is related to the distribution function in the phase-space, f, by F(t, μ, r, p) = A(r) f(t, μ, r, p), where A(r) is the cross-sectional area of the IMF flux tube (Lario et al. 1998). The local source of particles per unit area in the magnetic flux tube is described in this equation by the function G(t, μ, r, p). The injection rate of shock-accelerated particles in the phase-space (Heras et al. 1992), Q(cm⁻⁶ s⁻³/s) (i.e. the distribution function rate), is given by Q(t, μ, r, p) = G(t, μ, r, p)/A(r), with A(r) estimated at the point where particles are injected by the shock, i.e. at the cobpoint. We identify Q with the “efficiency” of the shock as a particle accelerator, which means that Q accounts for the effectiveness of the shock in accelerating protons and the efficiency of the shock at injecting these protons into IP space. According to the theoretical models of particle acceleration, in oblique and quasi-parallel shocks, a turbulent foreshock develops where the scattering of particles is enhanced (e.g., Vainio & Laitinen 2007, 2008). The transport model assumes an absorbing inner boundary in the radial direction placed behind the cobpoint; that is, the influx of particles is set to zero (Lario et al. 1998). In the case of the presence of a shock + foreshock structure where the scattering processes are more effective, the particles would be scattered back (reflected) towards the upstream region. These particles are considered in Q in our model. That is, Q quantifies the rate at which shock-accelerated protons escape from the shock structure or are reflected by it and evolves in time and space as the shock propagates and the cobpoint scans different regions of the shock front (Lario et al. 1998; Aran 2007; Aran et al. 2007). In our modelling approach, the physical processes of particle acceleration by the shock are not explicitly considered (like in e.g., Heras et al. 1992; Kallenrode 1997; Wang et al. 2012). This approach, where we put the emphasis on deconvoluting the transport processes to derive the “effective” injection of particles (i.e. with the aim to quantify the number of particles (re-)accelerated by and escaping from the shock) without assuming any specific particle acceleration mechanism, can be considered as complementary to other models that put the emphasis on the modelling of the particle acceleration (see comments by Zank et al. 2000). A relevant model pertaining to the latter category that has been applied to analyse actual SEP events from a space weather viewpoint is the PATH (Particle Acceleration and Transport in the Heliosphere) code (Verkhoglyadova et al. 2009, 2010, 2012). Being a realization of the model by Zank et al. (2000), Rice et al. (2003) and Li et al. (2003), this code uses the theory of diffusive shock acceleration to model particle acceleration at shocks. It assumes an analytic approximation of the spectrum of the accelerated particles derived from the properties of the shock obtained through an MHD model. The PATH code has been used to provide predictions of the radial dependence (from 0.5 to 2 AU) of the peak intensities of SEP events generated by shocks (with the observer always connected to the same region of the shock front), but these dependences are restricted to <5 MeV/nuc particles since the inner boundary of the MHD shock modelling is located at 0.1 AU (Verkhoglyadova et al. 2012).

Besides the injection rate of shock-accelerated particles, Q, the other main parameter of the transport model is the particle parallel mean free path, λ_||. We consider processes of diffusion in pitch-angle, and the mean free path to be given by λ_|| = λ_||0 (R/R₀)^2−q (Hasselmann & Wibberenz 1970), where R is the particle magnetic rigidity and q is the power index of magnetic field fluctuations (see details in Lario et al. 1998, and references therein).

When analysing a particular SEP event, the model allows the user to simulate any given number of energy channels, n_p, depending on the available measurements. In the case of the Apr00 event we have reproduced the 0.59–80 MeV proton differential intensity-time profiles measured by ACE/EPAM, 0.59–1.90 MeV (Gold et al. 1998) and SOHO/ERNE, 1.8–80 MeV protons (Torsti et al. 1995), summing up a total number of n_p = 19 energy channels. In the case of the Dec06 event, we have used n_p = 25 channels, reproducing the 0.3–500 MeV proton differential intensities measured from various instruments: 0.58–1.90 MeV from ACE/EPAM, 1.8–13 MeV from SOHO/ERNE, 14–100 MeV STEREOA/IMPACT/HET (von Rosenvinge et al. 2008) and 80–500 MeV GOES/EPS (Sauer 1993).

The methodology followed to perform the fittings of the observed proton intensities is fully described in Lario et al. (1998) and Aran et al. (2007) and can be summarized as follows:

1. We fit the intensity-time profile (and the anisotropy profile if available) of one energy channel with mean energy, E ₀, that we use as reference. E ₀ varies from event to event depending on the available measured differential intensity channels. For example, in the case of the Apr00 event, E ₀ = 9.1 MeV and in the case of the Dec06 event, we chose E ₀ = 14.3 MeV. This step yields the values of G and λ _|| at this energy as well as their evolution.

2. We simultaneously fit the intensity-time profiles (and first-order anisotropies when available) of the other energy channels by assuming that G scales with energy as G(E) = G(E₀)(E/E₀)^−γ and that the mean free path scales with the particle rigidity where we take q = 1.6 (Kunow et al. 1991). This step often requires fine tuning of the adopted value of λ_|| in step 1. The process is iterative, until we find the values of γ that may vary in time (two or three values depending on the event) that fit the whole set of differential intensities. The values of γ might also change for different energy ranges, assuming a double (or even a triple) power-law for some time intervals.

At the end of this process, the model outputs are the intensity and first-order parallel anisotropy time profiles for the whole set of simulated energies as well as the evolution of λ_|| and G. For each energy, we determine the continuous evolution and spectrum of the injection rate of shock-accelerated particles, Q, from applying a polynomial fit to the evolution of G.

A complete description of the particle intensities and modelling of these events can be found in Aran et al. (2011). In the simulations presented here we have used the evolution of the cobpoint position derived from the MHD simulations and extraction procedure presented in Sections 2 and 3 (i.e. the black curves in Fig. 5). In the case of the Apr00 SEP event, the parameters defining the particle transport (i.e. the injection rate, G, and λ_||) are the same as those used in Aran et al. (2011), with the difference that now the magnetic connection with the shock front is established closer to the Sun (at 3.5 R_⨀ instead of 5.8 R_⨀) which yields a larger value of Q at the beginning of the event due to the smaller value of A(r). The top left panel of Figure 6 displays the evolution of the injection rate Q for nine of the modelled energy channels (colour coded). For both events, we do not show the whole set of modelled energies in order to render the graphics clear. Figure 6 shows that the energy spectrum of the source of injected particles for the Apr00 event softens as the shock propagates towards the observer and the cobpoint slides along the eastern flank of the shock. This means that the shock efficiency at accelerating particles diminishes gradually in such a way that the bulk of >30 MeV protons are accelerated and escape from the shock during the first 7 h (i.e. r < 46 R_⨀). Note that for all energies, Q decreases rapidly during the first 2 h (first vertical dotted line) and coincides with the rapid decrease of VR, displayed in the bottom left panel of Figure 6. Later on, Q remains roughly constant from 2 to 11.8 h for <21 MeV protons and after that, for r > 70 R_⨀ (second vertical dotted line in Fig. 6) the particle injection by the shock decreases for all energies while VR is slowly decreasing.

Fig. 6. Top: Evolution of the injection rate of shock-accelerated protons at the cobpoint, Q, for a subset of modelled energy channels: left, for the 2000 April 4 SEP event and right for the 2006 December 13 SEP event. The bottom panels display the evolution of the VR parameter (see text for details).

In the case of the Dec06 event, the connection with the shock front is established 25 min after the start of the associated X-ray flare, when the shock is at 3.4 R_⨀ from the Sun, closer than in the former simulation, 4.7 R_⨀ (Aran et al. 2011). As for the Apr00 event, this yields a larger value of Q at the beginning of the event. The transport parameters are similar to those employed in the former simulation, the main differences being the following: (i) we use a slightly smaller mean free path, λ_||0 = 0.2 AU instead of 0.3 AU, that we assume constant throughout the simulation; and (ii) we assume that the shock is the only source of particles as we do in the case of the Apr00 event; i.e. we do not assume any additional source fixed at the Sun for the first 25 min. In order to model the energetic storm component (ESP) seen at 1 AU for E < 1.8 MeV, we have assumed a region ahead of the shock with enhanced turbulence. Such an ad hoc foreshock is 0.07 AU wide, it is switched on 14 h prior to the shock arrival and it is characterized by a proton mean free path given by λ_||c = λ_||c0 (R/R₀)^α, where λ_||c = 0.01 AU and α = +2 (e.g., Heras et al. 1992).

The top right panel of Figure 6 shows the resulting evolution of Q for 12 of the 25 modelled energy channels. The energy spectrum at the beginning of the event is harder (less steep decline) than in the April00 event. Similarly to April00, the injection rate derived for the Dec06 event decreases quickly during the first 2 h and coincides with the decrease of VR (see the bottom panel). The injection of >70 MeV protons happens during this first period whereas most of the 70 > E > 9 MeV particles are injected by the shock during approximately the six first hours (second vertical dotted line), when the shock has propagated up to 46 R_⨀. For the lower energies, Q remains constant and only increases up to the shock arrival for E < 5 MeV (see discussion in Sect. 4.5).

4.5. Synthetic intensity-time profiles

From Figure 6 one can infer that there does not exist a single linear relation between log Q and VR that applies for the whole duration of the event and for all energies, especially in the case of the Dec06 event. This is confirmed when plotting Q vs. VR, as depicted by cross-symbols in Figure 7, for the two events (left Apr00, right Dec06). However, these plots do suggest that it is possible to define two sets of linear relations that fit the data, as was obtained by Aran et al. (2011). The resulting linear fits are plotted as solid lines in Figure 7. The parameters defining the Q(VR) relation (Eq. (13)) found per each energy are listed in Table B.1 for the April00 event and in Table B.2 for the Dec06 event. The regions of applicability of each Q(VR) relation are marked by blue solid vertical lines in Figure 7 that coincide with the dotted vertical lines in Figure 6.

Fig. 7. Correlation between Q and VR at the cobpoint position for selected modelled energy channels (colour coded crosses): left, for the 2000 April 4 SEP event and right for the 2006 December 13 SEP event. Straight lines show the obtained Q(VR) relations (detailed in Appendix B). Vertical solid lines mark the radial distance boundaries of application of the semi-empirical relations found for each event.

In the case of the Dec06 event, to study the correlation between Q and VR, we have applied a hyperbolic fit to the values of VR (orange curve in the right bottom panel of Fig. 6) in order to deal with the inherent errors involved in determining VR in the interval from 7 to 13 h (see the discussion in Sect. 4.3). During this time interval, the cobpoint moves from the main body of the modelled disturbance to the small shock front developing ahead of its western wing (the case for the lowest observer in Fig. 4) and as a consequence the shock extraction procedures are more sensitive to the threshold values used when determining automatically both the candidate cobpoint and the location of the start of the downstream region.

For the Apr00 event, the linear fits are performed for the values of VR and log Q corresponding to the cobpoint radial positions from 3.5 R_⨀ to 18 R_⨀ (the inner region) and from 70 R_⨀ to 1 AU (the outer region). The correlation coefficients, ξ, of the inner region (column 4 of Table B.1) are larger than 0.96 for all energies. The values of the slopes k (see Eq. (13)) vary between 0.07 and 0.15 (column 3 of Table B.1). They are very similar for all energies and slightly larger than those determined by Aran et al. (2011) for the same event (0.03 < k < 0.08). In the case of the outer region, the correlation coefficients are less good than for the inner region, but ξ > 0.85 for all energies. These fits could be improved by smoothing the values of VR, but we did not perform that as we want to highlight the effect of the very small oscillations in VR on the derived intensity-time profiles (see below). For E ≤ 9 MeV, the values of k increase with the energy of protons, from 1.67 to 6.99, whereas for higher energies, k remains roughly constant. These values are larger than those found in Aran et al. (2011), where k varied from 0.48 to 1.06. Note that for E > 40 MeV and for r > 70 R_⨀ the values of Q are very low, more than three orders of magnitude lower than those at the beginning of the event. Hence, for these energies, the Q(VR) relation found in the outer region, functions only as a means to provide a numerical continuity for the source function. The shock is not accelerating these particles efficiently any more. In the intermediate region, 18 R_⨀ ≤ r ≤70 R_⨀, we have connected the two Q(VR) relations by means of a straight line also shown in Figure 7, with slope,

In order to investigate how the intensity-time profiles obtained by using the assumption that the source of particles is given by the Q(VR) relations compare with observations, we use VR(t) to compute Q(t) at the cobpoint, and use this source as input to the particle transport model. Figure 8 shows the resulting synthetic intensity-time profiles (black lines) and the observed data (red dotted curves) for 10 energy channels. In the case of the two lowest energy channels, we also show the fit to the evolution of the parallel (to the magnetic field) component of the first-order anisotropy (A1/A0), derived from the ACE/EPAM/LEMS30 and LEMS120 sectored data computed as in Lario et al. (2004). It can be seen that the synthetic profiles reproduce the observations despite the assumption made on Q in the intermediate region. The small irregularities shown by the synthetic time profiles at E > 5 MeV are a direct reflection of the oscillations present in the profile of the VR parameter (see bottom left panel in Fig. 6). At lower energies, these irregularities are not seen because the transport processes – that are more effective on low-energy protons – have smoothed out the changes of the values of Q.

Fig. 8. Observed (dotted red curves) and calculated (solid black curves) proton differential intensities for 10 out of 19 energy channels simulated for the 2000 April 4 SEP event. The calculated profiles are obtained by using the Q(VR) fits displayed in the left panel of Figure 7 and listed in Appendix B. The two top left panels also show the modelled and observed first-order anisotropy component parallel to the magnetic field, A1/A0. Arrows mark the starting time of the X-ray flare and vertical lines the time of the shock crossing by the ACE spacecraft.

Similarly as for the Apr00 event, we have performed two piecewise linear fits to the values of log Q and VR obtained for the Dec06 event. The first piecewise fit covers the spatial region from 3.4 R_⨀ to 16 R_⨀. The values of Q₀, k and ξ, defining the Q(VR) relation for each energy, are shown in columns 2–4 of Table B.2, respectively. The second fit starts from 46 R_⨀ and extends up to 1 AU, and the values of the corresponding parameters are listed in columns 5–7. The straight lines in the right panel of Figure 7 depict the two fits, together with the linear connection used in the intermediate region. The first vertical line marks the value of VR at r = 46 R_⨀ whereas the second vertical line marks the value of VR at r = 16 R_⨀. In the first region, k takes slightly different values across energy, ranging from 0.33 to 0.68. These values are larger than in the modelling by Aran et al. (2011), where k varied between 0.07 and 0.24. The linear correlation coefficients are larger than 0.94 for all energy channels with the exception of the 165–500 MeV channel, for which the injection rate has already decreased at r = 16 R_⨀ (see right panel of Fig. 7). In the case of the outer region, we only find a Q(VR) relation for E ≥ 6.4 MeV, for which ξ ≥ 0.96, except for the 6.4–8.1 MeV channel, for which ξ = 0.83. For E < 5.1 MeV, the slope k adopts negative values. This means that Q is increasing while VR is decreasing, indicating that for these energies the evolution of Q might not be determined by VR. Two factors may explain this: (i) as can be seen in Figure 9, the pre-event intensity levels for <10 MeV protons were high (red dotted curves). This is due to the fact that the Dec06 event occurred after a series of events on 5–8 December, originating from the same active region. These previous SEP events may have provided a supply of energetic seed particles that were re-accelerated at the Dec06 IP shock (e.g., Gopalswamy et al. 2004; Desai & Burgess 2008; Lee et al. 2012, and references therein); in this way, the increase of Q cannot be ascribed to the strength of the shock but rather to the extra pool of available particles, (ii) on the other hand, the way we model the foreshock region in conjunction with the assumption of an inner absorbing boundary (the effects of which are more noticeable on the lower-energy population due to the enhanced scattering processes at the foreshock Aran 2007), may account for an artificial increase in Q. We are currently working on a way to improve these two assumptions of the model. In particular, we plan to substitute our ad hoc foreshock by a semi-analytic foreshock modelling approach based on self-consistent simulations of diffusive coronal shock acceleration (Vainio et al. 2014).

Fig. 9. Observed (dotted red curves) and calculated (solid black curves) proton differential intensities for 10 out of 25 energy channels simulated for the 2006 December 13 SEP event. The calculated profiles are obtained by using the Q(VR) fits displayed in the right panel of Figure 7 and listed in Appendix B. The top left panel also shows the modelled (black trace) and observed (red dotted trace) first order anisotropy component parallel to the magnetic field, A1/A0. Arrows mark the starting time of the X-ray flare and vertical lines the time of the shock crossing by the ACE spacecraft.

Figure 9 shows the comparison between observations (red) and the synthetic proton intensity-time profiles (black curves) obtained by using the Q(VR) relations for the 2006 December 13 event. For E > 7.20 MeV, the synthetic profiles reproduce the observed intensities, and hence the adopted Q(VR) relations work for the higher energies. In the case of the two lowest energy channels shown in Figure 9, we have shown the resulting profiles taking them as simple mathematical fits (|ξ| > 0.86). The synthetic profiles reproduce the observed ones, even for the 5.1–6.4 MeV channel for which the correlation coefficient of the fitted Q(VR) relation for the region r > 46 R_⨀ is close to zero, ξ = 0.088 (Table B.2). This is because the main injection of particles occurs when the shock is close to the Sun, with the magnetic flux tube becoming populated with protons injected within the first region, i.e. 3.4 R_⨀ ≤ r ≤ 16 R_⨀. This highlights the fact that, in order to apply physics-based models for space weather predictions of large gradual proton events, it is important that these models allow the description of moving particle sources starting from a few solar radii, i.e. from coronal distances.

As shown in this work and in previous studies (e.g., Lario et al. 1998; Aran 2007; Aran et al. 2007, the values of the parameters defining the Q(VR) relation vary from event to event and with the energy channel. In order to use this relation in a predictive way (i.e. to build a tool useful for space weather applications), one can employ average values (among the modelled events) of the slope k and of Q₀(E) and apply the relation to different shocks and locations of the observers. In this way, Aran et al. (2006) built the SOLPENCO tool that was later on verified against data (Aran 2007; Aran et al. 2008). In another approach, the SaP model and the Q(VR) relation were used to generate 5–200 MeV proton synthetic intensity-time profiles in the SEPEM/SOLPENCO2 tool (Crosby et al. 2015). In this case, the outputs of the model, using average values of k (across a range of energies) of the Q(VR) relation, were compared with six SEP events, and the synthetic peak intensities and fluences were scaled to the observed 1 AU values in order to obtain predictions of these quantities for observers located at other radial distances (from 0.2 AU to 1.6 AU) along the same IMF line passing through the observer at 1 AU. We classified 204 SEP events in the event list of the SEPEM statistical model according to the six types of modelled events (Crosby et al. 2015). In this way SOLPENCO2 provided the SEPEM statistical model (Jiggens et al. 2012) with radial dependencies of peak intensities and fluences that consider the contribution of IP shocks to SEP events. This enables the user of SEPEM to evaluate the SEP radiation environment that a given interplanetary mission would encounter during its orbit (Crosby et al. 2015). These results are publicly available on the SEPEM server website. At present, however, these two tools cannot be used in real-time prediction mode until a solar proxy for the initial speed and density of the initial disturbance that generates the travelling shock is found, and until a statistical validation to estimate the errors in the prediction of the timing and size of the events is performed.

5. Summary and conclusions

We have upgraded the “Shock-and-Particle” model (Aran et al. 2011; Jacobs & Poedts 2011) that forms the basis of the SEPEM/SOLPENCO2 tool. This tool provides predictions of the peak intensity and event fluence radial dependencies for 5–200 MeV protons in gradual SEP events. We have developed a new 2D MHD solar wind model from 1 R_⨀ to 1 AU that simulates the solar wind acceleration processes by assuming a coronal heating source in contrast to the polytropic approximation adopted in the previous version of the SaP model. We have also presented a new automated method to locate the shock front in the MHD simulation and to extract the plasma parameters at the shock. We have then applied this new SaP model to the simulation of two large SEP events already studied with the previous version of the SaP model.

From the two events modelled, we have shown that the magnetic (BR) and density (r) compression ratios are significantly affected by different solar wind backgrounds as previously demonstrated by Pomoell & Vainio (2012), the new values being in accordance with theoretical MHD predictions. Nevertheless, we show that the evolution of the VR parameter is similar for the two solar wind models. The two methods employed for determining the cobpoint position and the plasma parameters at the shock yield similar values, with maximum differences in the VR parameter of ΔVR = 0.5 occurring when the cobpoint changes position from the central shock front to the shock at the flank that develops in the western wing of the disturbance. The new automated method yields smaller oscillations in the values of the jump parameters than the previous method due to a more robust determination of the shock normal.

We have reproduced the measured proton differential intensity-time profiles for the 2000 April 4 SEP event (19 energy channels, from 0.59 MeV to 80 MeV) and for the 2006 December 13 SEP event (25 energy channels, from 0.3 MeV to 500 MeV). We have obtained the evolution of the injection rate of shock-accelerated particles and we have studied its correlation with the VR parameter at the cobpoint. We have found that linear relations can be established between both parameters, especially for >5 MeV protons, and we have determined the values of the parameters defining these Q(VR) relations. Therefore, in spite of the fact that the values of k are different from those derived using the previous SaP model, we find again linear relations between log Q and VR which enable the computation of synthetic intensity-time profiles for gradual SEP events.

The Q(VR) relation permits the development of software tools for the prediction of peak intensities, fluences and proton intensity-time profiles of SEP events for energies that are relevant to space weather; that is, energies at which protons are capable of traversing typical shielding conditions, E > 10 MeV (Feynman & Gabriel 2000). Tools like SOLPENCO and SOLPENCO2 will be required until physics-based modelling tools that include the processes of particle acceleration by shocks and wave-particle interactions can be routinely used to model SEP events at a similar level of accuracy (when compared with actual data) as provided by the SaP simulations described in this work.

We want to point out that the Q(VR) relation is limited because of its simplicity, since relevant parameters for the acceleration of particles at the shock front, like the θ_Bn angle (e.g., Verkhoglyadova et al. 2012; Battarbee et al. 2013), are not taken into consideration. Other factors, like different seed particle populations (e.g., Li & Zank 2005; Tylka et al. 2005; Battarbee et al. 2011; Vainio et al. 2014) or the inclusion of a variable magnetic field (i.e. non-Parker morphology of the IMF) as studied by different authors (e.g., Giacalone 2005; Kocharov et al. 2009; Guo et al. 2010; Kozarev et al. 2013), have not been included in SaP yet. Our next step, together with the improvement of the modelling of the foreshock region as commented on in the previous section, is to study the correlation between Q and other parameters, particularly with θ_Bn, in order to, eventually, find a multi-parameter dependence of Q.

Acknowledgments

We wish to thank Dr. David Lario for providing the particle anisotropies from ACE/EPAM/LEMS120 and 30 telescopes for the two studied events. The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under Grant Agreement No. 262468 (Spacecast). Angels Aran wishes also to thank the people of the SEPEM Consortium and TEC-EES for their support and collaboration through the years. Angels Aran, Blai Sanahuja and Rosa Rodríguez-Gasén were also partially supported by the Spanish Ministerio de Economía y Competitividad, under the projects AYA2010-17286 and AYA2013-42614. Numerical MHD results were obtained on the HPC cluster VIC of the KU Leuven. Computational support was also provided by the Consorci de Serveis Universitaris de Catalunya (CSUC). SOHO is a project of international cooperation between ESA and NASA. We thank the ACE SWEPAM & MAG instrument teams and the ACE Science Center for providing the ACE data. We also acknowledge the use of STEREO/HET and GOES/SEM data available in the Internet. The editor thanks two anonymous referees for their assistance in evaluating this paper.

Appendix A

Extracting the shock data

In this appendix, the procedures of locating the shock in the MHD simulation data and extracting the plasma parameters of the shock are presented. An overview of the shock extraction process is provided in Section 3. Here we detail the methods used in each subtask of the extraction process and show an example of the result of the algorithm.

A.1. Determination of candidate cobpoint

The task of locating the cobpoint is carried out by using a thresholding procedure in which a set of quantities along the magnetic field line are computed, and if these quantities exceed certain threshold values, a candidate cobpoint is deemed to have been found. The notion candidate is used here, since additional procedures are later on applied in order to assert that the disturbance found by the thresholding is a shock and not a discontinuity or a large-amplitude wave.

The thresholding is performed on perturbations relative to the background solar wind, defined for a quantity q as(A.1)where the SW subscript indicates that the quantity is computed from the steady-state solar wind solution. At the leading edge of the shock front, the relative quantities are ideally all zero in the upstream, and only become non-zero at locations where the CME has perturbed the interplanetary medium. However, due to inherent diffusion in the numerical MHD solver, any structure will be smeared out over a number of grid cells. Moreover, at the flanks of the CME-driven shock, waves excited at the leading edge can be present in the upstream of the shock at the flank, causing the relative quantities to be non-zero. For these reasons, one cannot check simply for the existence of non-zero values. Instead, the quantities should exceed given thresholds. In other words, while following the magnetic field line from the observer towards the Sun, the truth value of the expression(A.2)is computed, and the first occurrence where it is true is flagged as a cobpoint candidate. In the above, S ≡ P/ρ^γ is the entropy and ε_q represents a constant threshold for quantity q that is determined empirically. In practice, we have found it sufficient to test for (Sρ)_rel > ε_Sρ only. In this work we use ε_Sρ = 1.5.

A.2. Computing the shock normal

The direction of the normal of the shock front at the cobpoint is an essential piece of information needed in subsequent calculations. The procedure for determining it is to compute the gradients of the quantities introduced in Section A.1. The rationale is that as the gradient is orthogonal to the level curves of the quantity in question, the gradient will be a good approximation of the normal direction since the quantities exhibit a steep gradient in the direction across the shock.

The gradient of quantity q is computed using standard finite differences,(A.3)where Δr ≡ Δ_grad and Δϕ = 2 arctan(Δ_grad/2r) so that an equal separation in both directions is used. The distance Δ_grad needs to be chosen to be larger than the local grid spacing so as to capture changes in the data.

The gradient is computed at the location of the cobpoint candidate. However, computing the gradient only at one location is error prone, as the vicinity of that one point may not be representative of the gradients occurring at the shock. Therefore, the gradient is computed at three locations, not only at the candidate cobpoint, but also at neighbouring points in the upstream as well as in the downstream direction at a distance d_nc along the field line from the candidate position. The median of the direction slopes evaluated at the three locations is then given as the direction of the shock normal. To increase the robustness of the normal computation even further, this procedure is applied to four quantities separately, i.e. for q = ρ_rel, (Sρ)_rel, v_r;rel, T_rel. As each application gives one estimate for the direction, the median is again computed in order to arrive at the final direction of the normal.

A.3. Determination of the upstream and downstream regions

Once a shock normal has been obtained, the values of the MHD variables along a line parallel to the normal are computed. This one-dimensional line-out of data is used for finding the location of the start of the downstream and upstream regions, and is based on analysing the properties of the derivatives of several quantities along the shock normal.

A.3.1. Shock normal line-out

Before the upstream and downstream regions can be found, the one-dimensional line-out along the shock normal needs to be obtained. Given the slope of the shock normal, computed as described in Section A.2, this is in principle a simple task. However, a subtle issue arises in choosing an appropriate distance between points that constitute the line-out data. If the distance chosen is too small, the bi-linear interpolation from the grid cell data will use the same data as input. This causes the interpolated data in the line-out to behave in a linear fashion, which results in derivatives along the shock normal used in determining the up- and downstream regions becoming inaccurate. A pragmatic solution is to sample the simulation data for the line-out at an interval that is of the order of the local grid spacing.

A.3.2. Locating the upstream region

The start of the upstream region is determined by thresholding the second derivative of the relative temperature. First, a normalized second derivative is obtained:(A.4)

The second derivative q″ ≡ d²q/ds² is computed using a standard centred second-order finite difference formula. The path coordinate s along the one-dimensional line-out runs from −L to L, with s = 0 corresponding to the cobpoint candidate position. The start of the upstream region is given as the location where the inequality(A.5)first holds, where the scan is started from the point corresponding to the most negative values of s and working towards larger values of s, i.e. from the upstream direction towards the downstream. The zero level z is the value of that corresponds to the value of q″ = 0, i.e. z = –min(q″)/(max(q″) – min(q″)). The normalization is done simply for the purposes of plotting since and aids in finding an appropriate value for the constant u_thres.

Note that the procedure is applied to the relative temperature, q = T_rel. This choice for q is somewhat arbitrary, and another quantity such as the relative entropy or relative entropy density could equally well be used. We have used the value u_thres = 0.1 for the threshold in this work.

A.3.3. Locating the downstream region

To find the downstream region of a quantity q (for instance, q = (Sρ)_rel) given a 1D line-out of the data across the shock front, the following procedure is applied:

1. Find all local extrema of q′ = dq/ds that succeed the start of the upstream region, i.e. for s such that s > s _upstream. This is done by computing d² q/ds ² using a standard second-order finite difference approximation and identifying the locations where the second derivative becomes zero by monitoring sign changes.

2. Compute q′ using second-order accurate finite differences, and find the locations beyond the start of the upstream where the derivative becomes zero. This amounts to locating the local extrema of q.

3. Given these locations, the start of the downstream region can be found. This is either given by the first occurrence of q′ = 0 following the first maximum of q′ after the upstream, or by the first minimum of q′ following the first maximum of q′. Discriminating between the two is done by checking which is closer to the first maximum of q′ following the upstream: the closer one is deemed to be the start of the downstream region.

The rationale behind the algorithm is to utilize the general characteristics of (diffused) shock fronts, as revealed by the derivatives, to locate the start of the downstream region: after the initial steep rise in q (a local maximum for q′), a phase where q increases less rapidly must follow (local minimum of q′).

The procedure presented above is carried out for several variables, i.e. q = S_rel, (Sρ)_rel, T_rel and q = ρ_rel. Each application gives generally a different downstream start location. The start of the downstream region is then chosen to be the one that is closest to the start of the upstream region.

A.4. Shock verification

With the upstream and downstream regions determined, a verification computation is performed to ensure that the obtained positions are consistent with a shock solution. This is done by computing the upstream and downstream Mach numbers (fast, slow as well as Alfvénic Mach numbers) in the frame of the shock and discarding solutions where both the upstream and downstream flows are superfast or subslow.

A.5. Example of upstream and downstream localization

Figure A.1 presents an example of the visual output of the shock detection algorithm. In each panel, the horizontal axis is the distance in solar radii along the normal of the shock along which the quantities have been computed, with zero corresponding to the cobpoint candidate position as described in Section A.1. The figure shows the quantities (Sρ)_rel, S_rel, T_rel discussed in the previous sections, as well as the density compression ρ_c = ρ_rel + 1 and normalized radial velocity v_r;rel. All quantities except ρ_c are normalized using Eq. (A.4). The magenta circle (box) indicates the location of the start of the downstream (upstream) region as computed by the algorithm.

Fig. A.1. A selection of plasma parameters along the shock normal for the cobpoint corresponding to the most central observer in Figure 4. See text for details.

The example shown corresponds to the shock crossing at the cobpoint of the most central observer shown in the snapshot in Figure 4.

Appendix B

The Q(VR) relations

In this section the parameters defining the Q(VR) relations (see Eq. (13)) are listed for the 2000 April 4 and 2006 December 13 SEP events in Tables B.1 and B.2, respectively. For each event, the values of the mean energy of the modelled channels are displayed in column 1. For the inner region, columns 2 and 3 show the values of the Q ₀ and k parameters, respectively, whereas column 4 shows the values of the linear correlation coefficient, ξ. Columns 5–7 show the same information for the Q(VR) relation found for the outer region, as indicated in each table.

References

All Tables

All Figures

	Fig. 1. 2000 April 4 event. From left to right, top to bottom: the solar wind density and radial velocity, and the magnetic field strength and the solar wind temperature measured by the ACE spacecraft (red curves). Corresponding values from the simulation are plotted in black. Time is counted from the time of the interplanetary shock passage by the ACE spacecraft.
In the text

	Fig. 2. 2006 December 13 event. From left to right, top to bottom: the solar wind density and radial velocity, and the magnetic field strength and the solar wind temperature measured by the ACE spacecraft (red curves). Corresponding values from the simulation are plotted in black. Time is counted from the time of the interplanetary shock passage by the ACE spacecraft.
In the text

	Fig. 3. The plasma variables as a function of radial distance from the Sun preceding the 2000 April 4 (left) and 2006 December 13 (right) CME-driven shocks.
In the text

Fig. 4. Radial velocity at t = 20.5 h after the CME onset for the 2006 December 13 event. Five observers (dots) with their corresponding magnetic field lines connecting with the shock front are displayed in grey. The observer, corresponding to the near-Earth spacecraft, is the one in the middle. Their corresponding cobpoints are the circles on the shock front. Black lines across them indicate the shock normal whereas magenta lines show the radial direction.

In the text

Fig. 5. The cobpoint position and plasma jumps at the cobpoint for the 2000 April 4 (left) and 2006 December 13 (right) SEP events. From top to bottom: The first and second panels show the radial position and the angular position of the cobpoint with respect to the launch direction of the shock. The third, fourth and fifth panels, respectively, show the normalized downstream-to-upstream radial velocity ratio, VR, the jump in magnetic field strength, BR, and the density compression ratio, r, respectively. The sixth panel displays the Mach number, M; the dashed horizontal line marks M = 1. The bottom panel shows the transit speed of the shock up to each cobpoint position. The black curves correspond to the results applying the methodology described in Section 3. Blue curves display the results obtained with the cobpoint procedure used in SEPEM but with the present MHD simulation model, and the red curves show the same quantities as derived from the model developed during the SEPEM project.

In the text

	Fig. 6. Top: Evolution of the injection rate of shock-accelerated protons at the cobpoint, Q, for a subset of modelled energy channels: left, for the 2000 April 4 SEP event and right for the 2006 December 13 SEP event. The bottom panels display the evolution of the VR parameter (see text for details).
In the text

Fig. 7. Correlation between Q and VR at the cobpoint position for selected modelled energy channels (colour coded crosses): left, for the 2000 April 4 SEP event and right for the 2006 December 13 SEP event. Straight lines show the obtained Q(VR) relations (detailed in Appendix B). Vertical solid lines mark the radial distance boundaries of application of the semi-empirical relations found for each event.

In the text

Fig. 8. Observed (dotted red curves) and calculated (solid black curves) proton differential intensities for 10 out of 19 energy channels simulated for the 2000 April 4 SEP event. The calculated profiles are obtained by using the Q(VR) fits displayed in the left panel of Figure 7 and listed in Appendix B. The two top left panels also show the modelled and observed first-order anisotropy component parallel to the magnetic field, A1/A0. Arrows mark the starting time of the X-ray flare and vertical lines the time of the shock crossing by the ACE spacecraft.

In the text

Fig. 9. Observed (dotted red curves) and calculated (solid black curves) proton differential intensities for 10 out of 25 energy channels simulated for the 2006 December 13 SEP event. The calculated profiles are obtained by using the Q(VR) fits displayed in the right panel of Figure 7 and listed in Appendix B. The top left panel also shows the modelled (black trace) and observed (red dotted trace) first order anisotropy component parallel to the magnetic field, A1/A0. Arrows mark the starting time of the X-ray flare and vertical lines the time of the shock crossing by the ACE spacecraft.

In the text

	Fig. A.1. A selection of plasma parameters along the shock normal for the cobpoint corresponding to the most central observer in Figure 4. See text for details.
In the text