Introduction to SOLPENCO2


The SOLar Particle ENgineering Code, SOLPENCO (Aran et al. [1,2]) was developed to rapidly provide predictions of the proton intensity profiles and fluences for gradual events. These SEP intensity profiles are computed from the onset of the event up to the arrival of the associated interplanetary shock. The core of SOLPENCO contains a database of pre-calculated synthetic flux profiles of gradual proton events for different interplanetary scenarios. These synthetic events are built by combining a MHD model for the description of the propagation of the CME-driven shock and a particle transport model for the simulation of the particle propagation from the shock front to a given location in space, along the connecting IMF line. The interplanetary scenarios depicted are defined by:
  1. the initial speed of the shock perturbation (at 18 solar radii);
  2. the propagation conditions for shock-accelerated particles;
  3. the relative position of the observer in space with respect to the leading direction of the shock (i.e. the heliocentric radial distance of the observer and the heliolongitude of the solar activity that triggers the event).

The SOLPENCO tool is based on the MHD model developed by Wu et al. [3] for the description of shock propagation and in the particle code developed by Lario et al. [4]. SOLPENCO provides upstream intensity profiles, and peak intensity and fluence values for protons with energy between 0.125 MeV and 64 MeV, as seen by observers located at fourteen different heliolongitudes at 0.4 AU or at 1 AU.


The SOLPENCO2 tool is an improved version of SOLPENCO that significantly extends the range of heliocentric radial distances and proton energies. It contains a new database of pre-calculated synthetic SEP profiles, that can predictions of peak fluxes and fluences, and time duration values, for gradual proton (5 MeV < E < 200 MeV) events and for fourteen observers located at seven different radial distances, between 0.2 AU (Solar Orbiter Mission) and 1.6 AU (Mars’ orbit), and for two different solar wind regimes (see Figure 1). Within the frame of the SEPEM project, the outputs of SOLPENCO2 are used as inputs for the statistical model for SEP radiation estimation at distances other than 1 AU, a main aim of the SEPEM project. Therefore, at this first stage of the SEPEM tools release users will not have direct use of this tool.
Figure 1. Location of the observers in the slow (left) and fast (right) wind.

The SEP intensity profiles depend on both how efficiently protons are accelerated, and how the IMF irregularities modulate this population during its journey in space. To build SOLPENCO2 synthetic SEP events database, a new MHD model was developed for the simulation of the CME-driven shocks from near Sun to Mars orbit. The shock propagation model provides values of the MHD variables at the cobpoint [5]; these values are used as inputs to a particle transport model [4]. These models allow computing the injection rate of shock-accelerated particles at a given time, a variable source term in the transport equation describing the propagation of protons along the IMF, towards the observer.

Figure 1 depicts the longitudinal extension covered by SOLPENCO2 for each radial position of the observer. SOLPENCO2 includes simulations of the interplanetary shock propagation on top of two different solar wind regimes (slow and fast solar wind, right and left hand side of Fig. 1).

For instance, the solar source longitudinal span covered by the simulations provided in SOLPENCO2 at 1 AU extends from W85 to E65. In this range 14 different observers (or solar source sites) are considered: W85, W75, W65, W55, W45, W30, W15, W00, E15, E25, E35, E45, E55 and E65.

For each kind of solar wind, SOLPENCO2 includes 8 different shock simulations characterized by different times of shock arrival up to 1 AU: from 17 to 48 hours on top of the fast wind and from 24 to 67 hours on top of the slow wind. Also, to try to cover a larger number of possible scenarios than in SOLPENCO the shock simulations are also performed for two types of shocks according to their longitudinal extension: narrow and wide shocks. Hence, the data base will contain 125,440 SEP flux profiles (7 radial distances x 14 heliolongitudes x 4 particle transport conditions x 16 shocks x 2 solar winds).

In order to obtain intermediate solar-interplanetary scenarios not stored in its data base, the SOLPENCO2 tool interpolates among the four pre-calculated gradual SEPs profiles which have the closest angular positions and transit times of the shock (from the Sun to the observer) to the corresponding values of the given event. In this way, it is possible to produce within minutes a synthetic SEP event for a given observed event at 1 AU or at non-Earth distances.

The MHD simulation of interplanetary shocks

The evolution of collisionless, traveling and expanding interplanetary shocks from Sun to Mars orbit is modelled by means of a new two dimensional MHD model. For the background solar wind plasma a polytropic relation is assumed between density and pressure with a radial varying polytropic index. This model permits the application of two different regimes: fast/low density/high temperature or slow/high density/low temperature solar wind. The MHD model assumes spherical geometry with a radial extent in the equatorial plane from 1.03 solar radii to 1.6 AU. As the shock evolves, the cobpoint slides along the shock front scanning regions of different MHD properties. Figure 2 shows a snapshot of one of these shock simulations with several IMF lines. For the synthesis of each SEP event profile in the database, a simulation for the shock propagation was performed in order to compute the evolution of its strength along the shock front.

Figure 2. Snapshot of an interplanetary shock simulation propagating from Sun to ~0.5 AU. Solar wind plasma speed is colour coded (as indicated) and a set of IMF lines (white lines) is also shown.

The strength of the shock is characterized by means of the different plasma values at the cobpoint, the point of the shock front that magnetically connected to the observer [5]; namely: the normalized downstream-to-upstream plasma velocity jump, VR, the magnetic downstream to upstream ratio, BR, and the local angle between the upstream IMF and the normal to the shock front, θBn [4, 5]. As the shocks expands and evolves in space, the cobpoint slides along the front of the shock and hence, the plasma values change with time, according to the local unperturbed upstream – shocked downstream solar wind conditions. Figure 3 (from [6]) sketches this concept: particles accelerated at the front of the shock are injected in the magnetic field lines connecting the observer (black diamond) with the shock at the cobpoint (red point); in this particular scenario, the position of the cobpoint moves towards the nose of the shock (red arrow), consequently, the efficiency of the injection of shock-accelerated particles increases. Therefore, the values at the cobpoint, for a given time and a given observer considered in the data base of SOLPENCO2 can be derived.

Figure 3. Cartoon showing the evolution of the cobpoint with respect to the front of the CME-driven shock and the position of the observer, for a given solar-interplanetary scenario.

The particle transport equation

The model adopted in [4] for the simulation of the interplanetary propagation of protons uses the transport equation derived by [7]. The transport model requires as input the following variables derived from the shock simulations performed with the MHD model: the position of the cobpoint as the shock progresses and expands from near Sun, the transit speed of the shock and the values of VR, at the cobpoint. The basic parameters of the transport model, as used here, are the mean free path of the protons, λ, the injection rate of shock-accelerated protons in interplanetary space, Q (phase space density) and the dependence with the energy of this injection rate (characterized by a spectral index, γ).

This injection rate changes with time, as it depends on the existing shock jump conditions at the cobpoint. The Q variable is the source term in the transport equation: it gives the efficiency of the shock as injector of accelerated particles in interplanetary medium. The mean free path λ describes the propagation of particles in the diffusive-focused transport model which is the result of the interaction between energetic particles and IMF irregularities in the quasi-linear approximation [8]. For the description of some events, the presence of a turbulent magnetic foreshock region is required to reproduce the flux and the anisotropy values observed in the spiky ESP component at the shock arrival. This foreshock is represented by a region of a given width just ahead of the shock, whose particle mean free path is much smaller than that in the rest of the upstream medium [4]. The transport equation is solved for each shock scenario by using a numerical algorithm that assumes that the inner boundary limit is at the radial distance of the cobpoint, which evolves as the shock propagates.

Deriving the injection rate and its energy dependence

Up to now, the precise details on how the 'acceleration efficiency' of the evolving shock behaves are neither completely clear nor quantified yet. A conceptual point of the shock-and-particle model in which SOLPENCO2 tool relies is the connection of the evolving MHD variables at the cobpoint with the injection rate of shock-accelerated particles; in other words, the Q(VR)-relation derived by Lario et al [4], already applied in different situations (i.e., Aran et al. [9]). This relation implies to assume a functional dependence between Q and VR variables at the cobpoint:
where k is a proportionality factor. This expression allows the modeller to relate the dynamic evolution of the shock strength at the cobpoint to the rate at which shock-accelerated particles are injected into the interplanetary medium. For a given SEP event, the flux and anisotropy profiles are fitted for one energy channel E0. This yields λ0 and Q0, as well as their evolution in the upstream region of the event. The dependence of Q(E) and λ(E) on the energy E is then derived by obtaining the best possible eyeball simultaneous fit to a large number of observed flux profiles (more than twenty) at selected energies, plus simultaneous fit the corresponding first order anisotropy profiles when available [5].

To build the database of SOLPENCO2, averaged values for k and Q0 are assumed. Then, the variation of the injection rate for each scenario is computed by using Eq. (1) because the values of VR are derived from the MHD simulation of the corresponding shock.


SOLPENCO2 consists of multi-dimensional arrays containing either the proton flux or the cumulative fluence value for the whole set of energies and for 12,544 solar-interplanetary scenarios. The parameters defining these scenarios are the following:

Isolated SEP events at 1 AU have been identified in order to calibrate the synthetic flux profiles and to determine the contribution of the downstream region to the fluence of a given event. This has required identifying the solar origin of each event of the list as well as the time of shock arrival and the event duration. It has been also necessary to evaluate how to reproduce compound events as defined in the Statistical Analysis Tool. Finally, we produce the synthetic flux profiles at other locations than 1 AU by assuming that the observers are at the same IMF line as the observer at 1 AU (i.e. the Earth). This allows calculating the peak intensity and fluence values for different energies and at other six radial distances, from 0.2 AU to 1.6 AU. Figure 4 shows an example of the peak flux (dots, left vertical axis) and fluence (squares, right vertical axis) values derived, for two energies: 8.7 MeV (red) and 26.3 MeV (black).

Figure 4. Example of the radial dependences for the peak flux and fluence derived for a specific solar-interplanetary scenario, for two energies.
As can be seen, the derived peak flux and fluence dependences with radial distance are quite different depending on the proton energy considered, as well as with respect the inverse quadratic/cubic law frequently assumed, as observationally verified by and numerically pointed out [10].

This physics based model [11,12] was used to derive the radial variation of the peak intensity and fluence for six SEP event case studies. The SEPEM statistical model is based on 1 AU cleaned data from 1973 to 2009. The solar origin of 204 SEP events (compound and isolated) of the SEPEM Reference Event List (REL) between January 1988 and December 2006 was determined. Their intensity and maximum energy attained were analysed. Then, the REL events were classified into the six studied types in order to assign to them the corresponding radial dependences. This allows to scale peak intensities and fluences to other radial distances and make use of the SEPEM statistical modelling machinery available at 1 AU [13].

Case studies: Comparing SOLPENCO2 synthetic SEP events with 1 AU data from the REL list and the derived heliocentric radial dependences

Six SEP events were selected from the REL list in order to calibrate the synthetic intensity-time profiles generated to reproduce them with SOLPENCO2. These events are:
  1. the 4 April 2000 SEP event,
  2. the 6 and 10 June 2000 SEP events,
  3. the event on 29 March 2001,
  4. and the consecutive events on 13 and 14 December 2006.
Figures 5 and 6 show the comparison of the observed intensity-time profiles for the SEPEM standard differential channels and the synthetic profiles produced by SOLPENCO2 for each SEP event.

Note that in this section, we refer to an "SEP event" following the scientific definition: a particle enhancement generated by a specific solar activity (either flare or/and CME) and often associated with coronal and interplanetary shocks. In the statistical analysis (both at 1 AU and away from 1 AU) SEP events are automatically detected by defining an intensity threshold for a specific energy, and consequently an event defined like this combines different particle enhancements or SEP events according to the first definition. So, in the following, when we refer to compound events, these are particle enhancements automatically detected consisting of a series of SEP events defined by their solar and interplanetary (if identified) sources.

Figure 5. SOLPENCO2 synthetic flux profiles scaled per channel using: peak intensity (solid lines); upstream fluence (dashed lines). Observed and synthetic intensities are displayed for the mean energy, indicated at the right hand side of each panel (in MeV), of the SEPEM Reference Data energy channels. Peak intensities: observed (black triangles), simulated (green circles). Vertical lines mark the passage of an interplanetary shock by the ACE spacecraft at L1.

Figure 6. SOLPENCO2 synthetic flux profiles scaled per channel, as in Figure 5.

Two sets of scaling factors were established per event at 1 AU: for peak fluxes and for upstream fluences, between the observed and synthetic values. Consistency between the total fluence simulated at 1 AU and the observed value was checked by calculating the total-to-upstream observed fluence ratios.

For each simulated energy, the synthetic upstream flux profiles generated at 0.2, 0.4, 0.6, 0.8, 1.3 and 1.6 AU are scaled using these two sets of scaling factors. Note that the virtual observers are placed along the same interplanetary magnetic field line as the observer at 1 AU.

The downstream ratios for duration and fluence calculated from 1 AU REL events are applied to the simulated upstream duration and fluences to get the total event fluence and duration. The ratios applied to the synthetic events away from 1 AU are adopted for the sake of simplicity and due to the scarce amount of data available to verify these assumptions.

Figures 7, 8, and 9 show the intensity-time profiles seen by the virtual observers, for selected proton energies (right panels), and the corresponding radial variations for the peak intensity (circles) and fluences (squares). A power law fit of these quantities with the radial distance (dashed lines) and the derived radial index are indicated. Note that the radial variations found for peak intensity and fluence are independent of the scaling factors but a direct result of the physical model.

Figure 7. Left Panels: Flux profiles at the simulated radial distances for 8.7 MeV and 26.3 MeV protons. Vertical lines mark the time of the shock arrival at each virtual spacecraft. Right panels: Peak intensity and fluence variations with the radial distance of the observers.

Figure 8. Left panels: Flux profiles at the simulated radial distances for 8.7 MeV, 26.3 MeV (and 115 MeV protons for the 13 December 2006 event). Vertical lines mark the time of the shock arrival at each virtual spacecraft. Right panels: Peak intensity and fluence variations with the radial distance of the observers.

Figure 9. Left panels: Flux profiles at the simulated radial distances for 8.7 MeV, 26.3 MeV (and 115 MeV protons for the 14 December 2006 event). Vertical lines mark the time of the shock arrival at each virtual spacecraft. Right panels: Peak intensity and fluence variations with the radial distance of the observers.

Coupling results from SOLPENCO2 with the SEPEM statistical model—Toward a statistical model away from 1 AU

The ultimate goal to reach would be to compare synthetic events provided by SOLPENCO2 with the 204 SEP events contained in the REL event list from January 1988 to December 2006 in order to reproduce the 1 AU REL list at the other six radial distances modelled. Since this is rather a titanic task, in the mean time, it is possible to use the radial dependences found in the studied six events to classify the whole set of events in the REL list during the aforementioned period. This subset of the REL list is referred to as The SEPEM radial dependent reference proton list.

The SEP events were classified according to the following characteristics:

The solar origin and shock passages at 1 AU were identified for each of the SEP event in the REL list from January 1988 to December 2006 (204 events). Figure 10 shows a screen shot of the event classification.

Figure 10. Screen shot of the Excel file containing the 204 SEP events identified and classified in the 1988 to 2007 period. REL uses the SEPEM reference data set based on cleaned GOES and IMP-8 measurements.

From the data the following characteristics were determined:

In compound SEP events the downstream region of one event is overlapped by the onset of the following event. The contributions of any remaining flux resulting from the first shock or of re-acceleration by the following shock are included in the later event. The method has little background physical support. It is aimed to give an operative input for the statistical model.

With this method, the SEP events were classified into six types:

Type 1—Events without observed shock
No shock is detected in the L1 data. Verified that events are of short duration (< 2–3 days) and presence of >66&nsbp;MeV component.
Type 1a. Peak intensity >10 cm-2 sr-1 s-1 MeV-1 for E=8.7 MeV. Use radial dependences from the 14 December 2006 SEP event.
Type 1b. Peak Intensity <10 cm-2 sr-1 s-1 MeV-1 for E=8.7 MeV. Use radial dependences determined for the 10 June 2000 SEP event.

Type 2—Gradual, low energy cases
A shock is detected in the L1 data and there is an absence of particles in the E>66 MeV energy range.
Type 2a. Solar origin between W90 and E15: use radial dependences determined from modelling the 4 April 2000 SEP event.
Type 2b. Solar origin between E15 and E90: use instead the radial dependences from the 6 June 2000 SEP event.

Type 3—Gradual, high energy cases
Type 3a. Peak intensity >50 cm-2 sr-1 s-1 MeV-1 for E=8.7 MeV. Use radial dependences from the 13 December 2006 SEP event.
Type 3b. Peak intensity <50 cm-2 sr-1 s-1 MeV-1 for E=8.7 MeV. Use instead radial dependences determined from the 29 March 2001 SEP event.
The distribution of event types over the events in the REL is given in the table below.

Distribution of event types in the REL
Event TypeFraction in REL
Type 1a3%
Type 1b7%
Type 2a29%
Type 2b21%
Type 3a24%
Type 3b10%
Not classified6%

Once the SEP events in REL are classified into one of the six types, the peak intensity and fluence spectra across radial distances can be computed. Figure 11 shows an example of the resulting spectra for the individual SEP event on 13 May 2005. This event was classified as a Type 2a event. In this case, the spectrum above 66 MeV is an extrapolation from the modelled (synthetic) event. Note that at these energies the event is not included in the statistical model using the recommended default flux and fluence thresholds.

Figure 11. Peak intensity (left) and fluence (right) spectra of the 13 May 2005 SEP event for each modelled radial distance (colour coded as indicated in the legends).

In the case of compound events, the peak intensity of each proton enhancement has been scaled to other radial distances according to its type and the selected value is the highest among the SEP enhancements for each energy and radial distance. To compute the fluences of the compound events, first the fluence of each individual SEP event (or each SEP enhancement) is scaled according to the power law of its type. Next, the scaled fluences are added to account for the fluence of the compound event. Figure 12 shows the example of the Halloween events in October–November 2003. This compound event is formed by 5 SEP events, all of them classified as Type 3a. Note that the flat radial dependence of the fluence in the 12.6 MeV channel determined from the 13 December 2006 SEP event modelling translates into a reversal in the behaviour of the fluence spectra across radial distance from lower to higher energies. For the two lowest energy channels, the farther is the observer from the Sun the larger the measured fluence due to the continuous injection of particles by the interplanetary shock.

Figure 12. Intensity spectra for the Halloween events (26 October–9 November 2003).

Away from 1 AU analysis

The specific method used to include the radial dependent information into the statistical modelling may be outlined as follows:
  1. The total mission period is divided up into successive segments of solar maximum and solar minimum periods, using 7 and 4 year windows around solar maximum (2.5 years before solar maximum to 4.5 years after) and solar minimum, respectively. For missions in the future, the current solar maximum date is repeated with an 11 year cycle.
  2. After a sampled waiting time an SEP event is generated with an associated duration as if the spacecraft were at 1 AU.
  3. The fluence/peak flux of the events is scaled for the spacecraft orbital position at that instance.
  4. This fluence/peak flux is recorded, the next waiting time is generated, the next event size found and scaled, and then added to the total (or compared to the running maximum flux).
  5. This procedure is repeated until the virtual timeline reaches the planned mission end.
  6. The whole process is repeated for 100,000 iterations, for the current solar cycle phase, and the results are stored.
  7. When all solar cycle phase segments have been processed, the results for the individual segments are accumulated.

Future work

The present work emhpasizes the relevance of characterising the downstream region of SEP events, both to quantify its contribution to the total fluence and to determine and separate isolated events.

The present results may be ameliorated by improving also the shock-and-particle model presented here, specifically the way the ESP component of SEP events is modelled. This is currently being done under the EC FP7 project SPACECAST.

More SEP events should be modelled, paying special attention to the effect of intervening structures to enhance SEP intensities [14]. The Solar Orbiter mission may help to improve and verify our models.

The work presented here is documented in a series of internal ESA/SEPEM Technical Reports. Part of these results is still in the process of publication.


[1] A. Aran et al. (2006) Advances in Space Research 37, 1240.

[2] A. Aran et al. (2005) Annales Geophysicae 23, 3047.

[3] S. T. Wu, M. Dryer and S. M. Han (1983) Solar Phys. 84, 395.

[4] D. Lario, B. Sanahuja, and A.M. Heras (1998) Astrophys. Journal 509, 415.

[5] A. M. Heras et al. (1992) Astrophys. Journal 391, 359; A.M. Heras et al. (1995) Astrophys. Journal 445, 497.

[6] A. Aran (2007) Ph. Thesis, Universitat de Barcelona.

[7] E. C. Roelof (1969), in Lectures in High Energy Astrophysics, NASA SP-199, ed. H. Ögelman & J. R. Wayland, 111.

[8] J. R. Jokipii (1966) Astrophys. Journal 146, 480.

[9] A. Aran et al. (2007) Astron. & Astrophys. 469, 1123.

[10] D. Lario et al (2006) Astrophys. Journal 653, 1531; D. Lario et al. (2007) Advances in Space Research 40, 289.

[11] Aran, A., C. Jacobs, B. Sanahuja et al., submitted A&A, 2012.

[12] Jacobs C., S. Poedts, Advances in Space Research, 48, 1958, 2011.

[13] Jiggens. P., S. Gabriel, D. Heynderickx, et al., IEEE Trans. Nuc. Sci., 59(4), 1066, 2012.

[14] Lario, D., A. Aran, R. Gòmez-Herrero, et al., Astrophys. J., 767, 41, 2013.

Last modified: 18 July 2013