Modelling large solar proton events with the shock-and-particle model Extraction of the characteristics of the MHD shock front at the cobpoint

We have developed a new version of a model that combines a two-dimensional Sun-to-Earth magnetohydrodynamic (MHD) simulation of the propagation of a CME-driven shock and a simulation of the transport of particles along the interplanetary magnetic field (IMF) line connecting the shock front and the observer. We assume that the shock-accelerated particles are injected at the point along the shock front that intersects this IMF line, i.e. at the cobpoint. Novel features of the model are an improved solar wind model and an enhanced fully automated algorithm to extract the necessary plasma characteristics from the shock simulation. In this work, the new algorithms have been employed to simulate the 2000 April 4 and the 2006 December 13 SEP events. In addition to quantifying the performance of the new model with respect to results obtained using previous versions of the shock-and-particle model, we investigate the semi-empirical relation between the injection rate of shock-accelerated particles, Q, and the jump in speed across the shock, VR, known as the Q(VR) relation. Our results show that while the magnetic field and density compression at the shock front is markedly different than in our previous modeling, the evolution of VR remains largely similar. As a result, we confirm that a simple relation can still be established between Q and VR, which enables the computation of synthetic intensity-time profiles at any location in interplanetary space. Furthermore, the new shock extraction tool is found to yield improved results being in general more robust. These results are important not only with regard to efforts to develop coupled magnetohydrodynamic and particle simulation models, but also to improve space weather related software tools that aim to predict the peak intensities, fluences and proton intensity-time profiles of SEP events (such as the SOLPENCO tool).


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 1983aLee , 1983bLee , 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 CMEdriven 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 tube1 .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) project2 , an upgrade of the SaP model was developed (Aran et al. 2011;Jacobs & Poedts 2011).This version incorporates a twodimensional (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 missions3 .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.

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. 2007Aran et al. , 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: giving the time dependence of the mass density (q), velocity (v), magnetic field (B) and total energy density (E).Furthermore, I is the unit dyad, g ¼ À GM r 2 e r 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 where c is the adiabatic index which is chosen to be constant.
Additionally, an ideal gas law is invoked for computing the temperature, where m is the mean molecular mass taken to be m ¼ m p =2 where m p is the proton mass.

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 allowing us to choose the adiabatic index to be constant and closer to the value of c = 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 q 0 q(r = r 0 ) and temperature T 0 T(r = r 0 ) in the low corona at r 0 = 1.03 R as well as the heating parameters S 0 and L and the adiabatic index c 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 q 0 , T 0 , S 0 , L, c 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 nonmagnetic non-rotating solar wind is a good approximation of the full system.Moreover, in this approximation, the steadystate 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 steadystate 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.

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. (2005Jacobs et al. ( , 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 (q extra ) and radial velocity (v extra ), have the same functional form as Rodríguez-Gasén (2011) and Aran et al. (2011).That is: where The initial disturbance is located at r cme = 1.75 R and centred around u = 180°.The disturbance only affects the region r 2 ½r cme À d cme =2; r cme þ d cme =2 and u 2 ½p À a cme = 2; p þ a cme =2: 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, A ¼ 15 R À1 : The parameter DU 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 DU the shape of the shock front can be varied from parabolic-like (DU = 0.5 a cme ) to circular-like (DU = 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 Dr = 0.003 R at the inner boundary to Dr = 1.09R 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 Du = 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.

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: (a) Compute the magnetic field line passing through the virtual spacecraft.(b) Determine the location of the cobpoint, i.e. the location where the field line computed in the previous step crosses the shock.(c) Compute the shock normal at the cobpoint.(d) Determine the locations of the regions upstream and downstream of the shock.(e) 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.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.

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).

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 À1 (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
Table 1.Overview of the parameters in the shock model.
. Space Weather Space Clim., 5, A12 (2015) 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.

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.As can be seen in these figures, the simulated solar wind (black curves) reproduces the large-scale values of the preevent solar wind conditions for both events.The timeindependent 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.
Note that the Apr00 event developed under slow solar wind conditions (v sw ~380 km s À1 ) whereas the Dec06 event occurred under fast solar wind conditions (v sw ~595 km s À1 ).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 (q cme , v cme and DU) 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 a in lieu of that) flare is used.
First, we keep the value of DU fixed, and we produce different test runs by varying q 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 1 2 q cme v 2 cme (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 q cme and v cme and that of DU, 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 multispacecraft 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).

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: 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 h 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): 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. 2006Aran et al. , 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).
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 DVR = 0.5 in this interval, while for the remainder of the event we have DVR < 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 q cme and a factor of 1.98 for v cme , with the shock being narrower since DU = 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), q 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 (DU 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).

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 focuseddiffusion transport equation used (Ruffolo 1995) accounts for the effects of particle focusing in the diverging IMF, pitchangle 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, l, 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, l, r, p) = A(r) f(t, l, 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, l, r, p).The injection rate of shock-accelerated particles in the phase-space (Heras et al. 1992), Q(cm À6 s À3 /s) (i.e. the distribution function rate), is given by Q(t, l, r, p) = G(t, l, 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(Verkhoglyadova et al. , 2010(Verkhoglyadova et al. , 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, k || .We consider processes of diffusion in pitch-angle, and the mean free path to be given by k || = k ||0 (R/R 0 ) 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.90MeV (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.90MeV 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 0 , that we use as reference.E 0 varies from event to event depending on the available measured differential intensity channels.For example, in the case of the Apr00 event, E 0 = 9.1 MeV and in the case of the Dec06 event, we chose E 0 = 14.3 MeV.This step yields the values of G and k || 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 0 )(E/E 0 ) Àc 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 k || in step 1.The process is iterative, until we find the values of c that may vary in time (two or three values depending on the event) that fit the whole set of differential intensities.The values of c 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 k || 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 k || ) 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.
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, k ||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 k ||c = k ||c0 (R/R 0 ) a , where k ||c = 0.01 AU and a = +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).

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.
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.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   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 n !0.96, except for the 6.4-8.1 MeV channel, for which n = 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).
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 (|n| > 0.86).The synthetic profiles reproduce the observed ones, even for the 5.1-6.4MeV channel for which the correlation coefficient of the fitted Q(VR) relation for the region r > 46 R is close to zero, n = 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 0 (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.

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 DVR = 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 h 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 h Bn , in order to, eventually, find a multi-parameter dependence of Q. A12-p20 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 À1 .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 À1 .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) andVerkhoglyadova et al. (2010).

Fig. 1 .
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 .
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.

Fig. 3 .
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.

Fig. 4 .
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.

Fig. 5 .
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.

Fig. 6 .
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).

-
that are more effective on low-energy protons -have smoothed out the changes of the values of Q. 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 0 , k and n, 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

Fig. 8 .
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 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.

Fig. 9 .
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.

Table 2 .
Solar wind parameters for the 2000 April 4 and 2006 December 13 events.

Table B .
1.The Q(VR) relations for the 2000 April 4 SEP event: values of k and Q 0 per energy channel and the linear correlation coefficient, n.Left: list of the values derived for the 3.5 R r 18 R region.Right: list of the values for the r !70 R region.Table B.2.The Q(VR) relations for the 2006 December 13 SEP event: values of k and Q 0 per energy channel and the linear correlation coefficient, n.Left: list of the values derived for the 3.4 R r 16 R region.Right: list of the values for the r !46 R region.Italic values indicate the energy channels for which a Q(VR) relation is not found (see text for details).