© A. Afanasiev et al., Published by EDP Sciences 2025

1 Introduction

Solar energetic particle (SEP) events constitute one of the significant aspects of space weather (Baker, 2005; Vainio et al., 2009). Of particular importance are gradual SEP events, which are associated with shock waves driven by coronal mass ejections (CMEs) (e.g., Reames, 1999), as they provide the highest proton intensities and fluences at MeV energies (Desai & Giacalone, 2016). Substantial effort has been addressed worldwide to develop forecasting methods of such SEP events, using empirical and physics-based modeling approaches (a recent review by Whitman et al., 2023, and references therein). The latter approach includes modeling of particle acceleration at shock waves. These simulations are computationally demanding if done self-consistently, i.e., if the effect of the turbulence amplification by accelerated particles is included. This high computational cost currently prevents such models from being directly used in an operational SEP forecasting framework. Consequently, alternative models have been developed, which approximate the particle acceleration processes through semi-analytical treatments.

An example of such model is the improved Particle Acceleration and Transport in the Heliosphere (iPATH; e.g., Hu et al., 2017). This model includes a particle acceleration module that, at each time step in the simulation, relies on steady-state analytical solutions for diffusive shock acceleration (DSA; e.g., Drury, 1983) at the shock wave to inject SEPs into the simulation. The utilized DSA solution relies on the local shock properties, which are obtained from a magnetohydrodynamic (MHD) simulation of a CME. The injected particles diffuse within an onion-shell model situated behind the CME-driven shock. Those particles that succeed in escaping upstream of the shock are tracked by iPATH’s transport module, which solves the focused transport equation (FTE; e.g., van den Berg et al., 2020) under the assumption of a uniform Parker solar wind configuration. Over the years, different versions of the iPATH model have been successfully employed to model several SEP events (see e.g., Li et al., 2021; Ding et al., 2022, for some recent examples).

A semi-analytical model of the foreshock region of an interplanetary shock driven by a CME was developed by Vainio et al. (2014). Its primary objective was to provide a rapid self-consistent characterization of the foreshock region, including particle intensities and mean free paths across a range of particle energies. This model is built upon a steady-state DSA theory given by Bell (1978), which accounts for self-generated Alfvén waves upstream of the shock, and the Coronal Shock Acceleration (CSA) simulation model (Vainio & Laitinen, 2007; Battarbee et al., 2013) capable of simulating coupled ion acceleration and Alfvén wave generation in a time-dependent manner in a spatially varying solar wind plasma. CSA simulations agree with Bell’s theory near the shock upon reaching the steady state (Vainio & Laitinen, 2007; Afanasiev et al., 2015). Therefore, the semi-analytical foreshock model employs, as an analytical framework, Bell’s theory with some modifications to account for the finite time span of the acceleration process and extent of the foreshock region. The self-consistent CSA simulations were fit to the obtained analytical framework in order to calibrate its parameters. The model provides a description of the foreshock region in terms of the energetic particle intensity distribution and the pitch-angle scattering mean free path as functions of energy and distance from the shock (along a magnetic field line), provided that the energy spectrum cutoff energy, the ambient plasma, and the MHD shock parameters are given.

In this paper, we introduce a new particle acceleration and transport model called PARASOL, which is based on the SOLar Particle Acceleration in Coronal Shocks (SOLPACS) model simulating proton acceleration at shocks, coupled with Alfvén wave generation (Afanasiev et al., 2015), and the PArticle Radiation Asset Directed at Interplanetary Space Exploration (PARADISE) model simulating three-dimensional (3-D) particle transport through the solar wind (Wijsen, 2020). SOLPACS treats the wave-particle interactions more accurately than CSA by using the full resonance condition of pitch-angle scattering, and has been recently successfully validated against energetic storm particle (ESP) event observations at 1 au (Afanasiev et al., 2023). PARADISE tracks the evolution of energetic particle distributions through solar wind configurations generated by the time-dependent 3-D MHD model known as the EUropean Heliospheric FORecasting Information Asset (EUHFORIA; Pomoell & Poedts, 2018). Moreover, EUHFORIA can be used to simulate the propagation of CMEs through the solar wind. When sufficiently fast, such a CME will drive a shock wave whose properties serve as input to the PARASOL model.

PARASOL includes a new semi-analytical foreshock model because directly integrating SOLPACS with PARADISE would result in a computationally heavy simulation model that is hardly usable in the operational space-weather context. Our approach to construct such a model is different from the approach undertaken in Vainio et al. (2014). This is because Bell’s theory can not be used in full since SOLPACS simulations give an inverse (compared to Bell’s theory and CSA) dependence of the particle mean free path on energy (Afanasiev et al., 2015), which results from the full resonance condition of pitch-angle scattering utilized in the model. Furthermore, SOLPACS is a local model as it performs simulations under assumption of spatially constant ambient plasma and shock parameters, so no focusing effects are included. As a result, the semi-analytical model developed here describes only the inner part of the foreshock (Fig. 1). While the inner foreshock is characterized analytically, the outer foreshock is handled numerically by means of simulations with the PARADISE model, which propagates the accelerated SEP distributions through the EUHFORIA-generated solar wind. Figure 2 presents a flow chart of PARASOL indicating the input required by the constituting models. It should be noted in this connection that PARASOL can be used not only with EUHFORIA, but with any other model that can provide the necessary bulk plasma and shock parameters.

Figure 1 Diagram of the PARASOL model. The thin black line upstream of the shock schematically depicts the interface between the inner foreshock modeled analytically based on self-consistent SOLPACS simulations and the outer foreshock simulated in PARADISE. It should be noted that the actual shock-normal distance from the shock front to the model matching interface depends on the shock magnetic geometry (see Sect. 2 for details).

Figure 2 Flow chart of PARASOL indicating the input required by the constituting models.

This paper is structured as follows. In Section 2, we outline the derivation of the semi-analytical model and its integration with PARADISE. In Section 3, we present the results of EUHFORIA + PARASOL simulations of an SEP event that was observed in-situ in July 2012. In Section 4, we discuss in detail some deviations of the simulation results from observations, as well as further potential developments of the model. Finally, in Section 5, we provide our conclusions.

2 Model development

2.1 General consideration

In SEP transport equations that do not address the particle acceleration process at the shock explicitly, the energetic particle emission from the shock is usually described by the source term, which represents the resulting accelerated particle population at the shock (e.g., Lario et al., 1998). Assuming the shock to be an MHD discontinuity, the source term can be given as

$$S(\mathbf{x},p,t)=Q(p,t)S_{\mathrm{sh}}(\mathbf{x},t)$$(1)

where Q(p, t) is the particle emission rate as a function of the particle momentum p and time t, and S_sh(x, t) describes the 3-D shock surface at time t. For instance, introducing heliocentric spherical coordinates, the latter function can be given as S_sh(r, θ, ϕ, t) = δ(r − r_sh(θ, ϕ, t))/(4πr² sinθ), where r_sh is the radial coordinate of the shock surface at the co-latitude θ and the longitude ϕ, and δ stands for the Dirac’s delta function (this particular form of S_sh assumes that the shock front can intersect a given line drawn from the origin (ϕ = const., θ = const.) only at one point at a given time).

The diffusive acceleration of particles at the shock can be described by the 1-D steady-state Parker equation (Drury, 1983):

$$u\frac{\partial f}{\partial z}-\frac{p}{3}\frac{\partial u}{\partial z}\frac{\partial f}{\partial p}=\frac{\partial }{\partial z}\left(\kappa \frac{\partial f}{\partial z}\right)+\frac{q}{4\pi p_{\mathrm{inj}}^2}\delta \left(z\right)\delta (p-p_{\mathrm{inj}}),$$(2)

where f is the isotropic part of the phase-space distribution function of particles, u is the scattering-center speed with respect to the shock and κ = κ(z, p) is the spatial diffusion coefficient, both measured along the z direction normal to the shock (z < 0 in the upstream). It is further assumed that seed particles are injected into the acceleration process at a rate q (per unit area of the shock and unit time) at a momentum p_inj. The resulting distribution function at the shock f_sh = f(z = 0) is then equal to

$$f_{\mathrm{sh}}\left(p\right)=\frac{3}{4\pi }\frac{q}{\Delta u{p}_{\mathrm{inj}}^3}{\left(\frac{p}{p_{\mathrm{inj}}}\right)}^{-\sigma },$$(3)

where Δu = u₁ – u₂, σ = 3u₁/Δu, and u₁ and u₂ refer to the shock-normal scattering-center speeds upstream and downstream of the shock, respectively.

Equation (2) can be reformulated to include instead of the seed-particle source term the particle emission term describing the resulting spectrum of accelerated particles:

$$u\frac{\partial f}{\partial z}-\frac{\partial }{\partial z}\left(\kappa \frac{\partial f}{\partial z}\right)=Q\left(p\right)\delta \left(z\right),$$(4)

where

$$Q\left(p\right)=-\frac{\Delta u}{3}p\frac{\partial f_{\mathrm{sh}}\left(p\right)}{\partial p}$$(5)

and f_sh(p) is given by equation (3).

We employ equation (5) to model the particle emission rate at the shock in PARADISE, taking into account that the transport equation solved in PARADISE is formulated in terms of the particle differential intensity j = p²f (Wijsen, 2020, also see Sect. 2.3):

$$\begin{array}{ll}{Q}^{\left(P\right)}\left(p\right)& =-\frac{\Delta u}{3}{p}^3\frac{\partial f_{\mathrm{sh}}^{\left(P\right)}\left(p\right)}{\partial p}\\ & =-U_1\frac{r_{\mathrm{sc}}-1}{r}\frac{p^3}{3}\frac{\partial }{\partial p}\left(p^{-2}{j}_{\mathrm{sh}}^{\left(P\right)}\right),\end{array}$$(6)

where U₁ is the normal component of the upstream bulk plasma velocity in the de Hoffmann-Teller (HT) frame, r is the conventional gas compression ratio of the shock, $r_{\mathrm{sc}}=u_1/u_2=r\left(1-M_{\mathrm{A}}^{-1}\right)$ is the scattering-center compression ratio of the shock, M_A = U₁/(V_A cosθ_Bn) is the Alfvénic Mach number of the shock as defined in the HT frame, θ_Bn is the shock-normal angle, and V_A is the Alfvén speed in the unshocked plasma. The above expression for the scattering-center compression ratio, r_sc, is obtained under the assumption that upstream of the shock particles are scattered by parallel forward (anti-sunward) propagating Alfvén waves (thus, u₁ = U₁ − V_A cosθ_Bn) and downstream of the shock particles are scattered by turbulence frozen-in to the bulk plasma (e.g., Vainio et al., 2014). In equation (6) and below, we use the subscript (P) to explicitly acknowledge quantities related to PARADISE (and the subscript (S) to SOLPACS).

According to equation (6), the particle emission rate Q^(P) is determined by the particle intensity spectrum at the shock, $j_{\mathrm{sh}}^{\left(P\right)}\left(p\right)$. To obtain the particle intensity in the foreshock, we need to specify the diffusion coefficient (or the scattering mean free path) in the shock vicinity as a function of momentum (or energy) and distance from the shock. We do this by making use of self-consistent simulations of particle acceleration at shocks performed with the SOLPACS code (Afanasiev et al., 2015, 2023, also see Sect. 2.2.1). SOLPACS provides accelerated proton distributions at various energies upstream of the shock, along with their corresponding scattering mean free paths. Due to the self-generated turbulence, the SOLPACS mean free paths in a close vicinity of the shock can be several orders of magnitude smaller than those (∼0.1 au) typically used in PARADISE simulations of particle transport in the interplanetary medium. The direct usage of such small mean free paths would slow down a PARADISE simulation substantially. We deal with this problem in the following way. We require that the particle omnidirectional flux in PARADISE matches the particle omnidirectional flux in SOLPACS only at a certain distance x_M (along the magnetic field line) upstream of the shock: $j_{\mathrm{M}}^{\left(P\right)}\left(p\right)=j_{\mathrm{M}}^{\left(S\right)}\left(p\right)$. The PARADISE flux at the shock, $j_{\mathrm{sh}}^{\left(P\right)}\left(p\right)$, is then obtained by “demodulating” the SOLPACS flux at the matching distance x_M, $j_{\mathrm{M}}^{\left(S\right)}\left(p\right)$, by using the solution of the steady-state Parker equation:

$$j_{\mathrm{sh}}^{\left(P\right)}\left(p\right)=j_{\mathrm{M}}^{\left(S\right)}\left(p\right)\exp \left\{{\int }_0^{x_{\mathrm{M}}}\mathrm{}\frac{u_{\parallel }}{{\kappa }_{\parallel }^{\left(P\right)}(x,p)}\mathrm{dx}\right\},$$(7)

where x = −z/cosθ_Bn is the distance along the magnetic field line with respect to the shock (x > 0 in the upstream region), u_∥ = u₁/cosθ_Bn and

$${\kappa }_{\parallel }^{\left(P\right)}(x,p)=\frac{{\kappa }^{\left(P\right)}}{{\cos }^2{\theta }_{\mathrm{Bn}}}=\frac{1}{3}v{\lambda }_{1 \mathrm{au}}^{\left(S\right)}(x,E)$$(8)

is the spatial diffusion coefficient along the magnetic field, determined by the parallel mean free path ${\lambda }_{1 \mathrm{au}}^{\left(S\right)}$ and the particle speed v (to simplify the notation we drop here and below the subscript “∥” in ${\lambda }_{\parallel }^{\left(S\right)}$ as in SOLPACS simulation the cross-field diffusion is not considered). For the mean free path in equation (8), we use an analytical function of distance and energy E that we infer by fitting the mean free paths resulting from SOLPACS simulations for the shock assumed to be located at 1 au. In this way, we can expect PARASOL to be able to predict ESP events measured at 1 au and, at the same time, to be more computationally efficient while simulating particle propagation in the vicinity of the shock (i.e., in the foreshock) when the shock is at distances <1 au.

The SOLPACS flux at the matching distance x_M is obtained via

$$j_{\mathrm{M}}^{\left(S\right)}\left(p\right)=j_{\mathrm{sh}}^{\left(S\right)}\left(p\right)\exp \left\{-{\int }_0^{x_{\mathrm{M}}}\mathrm{}\frac{u_{\parallel }}{{\kappa }_{\parallel }^{\left(S\right)}(x,p)}\mathrm{d}x\right\},$$(9)

with

$${\kappa }_{\parallel }^{\left(S\right)}(x,p)=\frac{1}{3}v{\lambda }^{\left(S\right)}(x,E),$$(10)

and $j_{\mathrm{sh}}^{\left(S\right)}\left(p\right)$ is the particle flux at the shock as modeled by SOLPACS.

One way to define the matching distance x_M is to require the residual modulation of the upstream particle flux as modeled in SOLPACS to be small, that is:

$$M_{\mathrm{res}}={\int }_{x_{\mathrm{M}}}^{x^0}\mathrm{}\frac{u_{\parallel }}{{\kappa }_{\parallel }^{\left(S\right)}}\mathrm{d}x<<1,$$(11)

where x⁰ is the distance along the magnetic field line to the point where the SOLPACS-modeled mean free path becomes equal to the mean free path λ⁰ used in PARADISE beyond the foreshock (see Eq. (43)). This allows one to use the larger mean free paths of PARADISE already at x > x_M in PARASOL, and at the same time to ensure that the escaping particle flux is modeled correctly. As an approximation, we calculate M_res assuming that x⁰ = x_feb, where x_feb is the distance along the magnetic field line from the shock to the free escape boundary in the upstream, which is equal to the simulation box size in SOLPACS (50 R_⊙ in all simulations in this work). Checking the approximated M_res in our simulations demonstrated that x_M > 10 R_⊙ at energies at which particle intensities and mean free paths in a close vicinity of the shock (x ≲ 0.1 R_⊙) are in a quasi-steady state. Moreover, x_M is energy-dependent. In this work, we shall treat x_M as a free parameter of the PARASOL model (and fix it to 10 R_⊙).

To summarize this section, in order to construct an analytical model of the foreshock, we need to obtain analytical functions for the particle intensity distribution j^(S)(x, E) and mean free path λ^(S)(x, E) in the upstream region (here and below we use particle energy E instead of momentum p). As detailed in the next section, we obtain those by fitting the mean free paths resulting from SOLPACS simulations with analytical functions that represent a modification of the mean free path function of Bell’s theory. These fitting functions are then used to evaluate equations (7) and (9). This, along with an analytical representation of the simulated intensity spectrum at the shock, $j_{\mathrm{sh}}^{\left(S\right)}\left(E\right)$, allows us to construct the foreshock model based solely on fitting the mean free.

2.2 Fitting of SOLPACS simulation results

2.2.1 SOLPACS simulation runs

In the SOLPACS simulation model (Afanasiev et al., 2015, 2023), protons and Alfvén waves are traced in a spatially 1-D simulation box (along the magnetic field line) between the shock front and the free-escape boundary placed in the upstream region. The ambient plasma and the magnetic field are assumed homogeneous in the simulation box, and the shock velocity is taken to be constant in a given simulation. The equations solved by SOLPACS can be given in the following form:

$$\frac{\partial f}{\partial t}+\left[\mathrm{v\mu }+\left(1-M_{\mathrm{A}}\right)V_{\mathrm{A}}\right]\frac{\partial f}{\partial x}=\frac{\partial }{\partial \mu }\left(D_{\mu \mu }\frac{\partial f}{\partial \mu }\right)$$(12)

$$\frac{\partial I}{\partial t}+\left(1-M_{\mathrm{A}}\right)V_{\mathrm{A}}\frac{\partial I}{\partial x}=\mathrm{\Gamma } I,$$(13)

where equations (12) and (13) describe the evolution of the gyro-averaged distribution function f(x, v, μ, t) of particles and the evolution of the wave intensity I(x, k, t), correspondingly. Here t is time, x is the spatial coordinate measured along the magnetic field line from the shock toward upstream, v and μ are the particle speed and the pitch-angle cosine as measured in the rest frame of Alfvén waves, and k is the wavenumber. In equation (12), the second term on the left-hand side describes particle streaming and convection, and the right-hand side term provides the pitch-angle scattering of particles. In equation (13), the second term on the left-hand side describes the wave advection toward the shock, and the term on the right-hand side describes the wave growth. The coefficients D_μμ and Γ are the quasi-linear pitch-angle diffusion coefficient and the wave growth rate, correspondingly:

$$D_{\mu \mu }\left(\mu \right)=\frac{\pi }{2}{\mathrm{\Omega }}_0\frac{\left|k_{\mathrm{r}}\right|I\left(k_{\mathrm{r}}\right)}{\gamma B^2}\left(1-{\mu }^2\right)$$(14)

$$\mathrm{\Gamma }\left(k\right)=\frac{\pi }{2}\frac{{\mathrm{\Omega }}_0}{n{V}_{\mathrm{A}}}\int \mathrm{}{\mathrm{d}}^{3}p v\left(1-{\mu }^2\right)\left|k\right|\delta \left(k-k_{\mathrm{r}}\right)\frac{\partial f}{\partial \mu },$$(15)

where B is the mean magnetic field magnitude, n is the plasma density, Ω₀ = eB/m_p is the (non-relativistic) proton cyclotron frequency (e is the electron charge, m_p is the proton mass), γ is the relativistic Lorentz-factor, k_r is the resonant wavenumber given by the quasi-linear resonance condition

$$k_{\mathrm{r}}=\frac{{\mathrm{\Omega }}_0}{\gamma \mathrm{v\mu }},$$(16)

and I(k_r) is the wave intensity at the resonant wavenumber. In equation (15), δ(k − k_r) is the Dirac delta-function of its argument, and the integration is performed over the particle momentum space.

The input parameter list of the SOLPACS model consists of five main parameters: the Alfvén speed V_A and density n of the ambient plasma, the Alfvénic Mach number of the shock M_A, the shock-normal angle θ_Bn and the particle injection efficiency of the shock ϵ_inj (Afanasiev et al., 2023).

In this study, we used fourteen SOLPACS simulation runs. The parameter values in each run are given in Table 1. The scaling property of the SOLPACS equations, α^1/2ϵ_inj = const., where α is the density factor (see Appendix A for details), allows us to narrow down the parameter space. To emphasize this feature, instead of the plasma density n in Table 1, we provide the density factor α, which is related to the density via n = αn_ref, where n_ref is the reference value of 5 cm⁻³.

2.2.2 Fitting particle mean free paths and deriving particle intensity distributions

As outlined above, our approach consists of fitting the simulated mean free paths and applying the steady-state solution of the Parker equation (e.g., Eq. (9)) to calculate the particle intensity distribution in the upstream region of the shock.

We found that the following function fits well the simulated mean free paths:

$${\lambda }^{\left(S\right)}(x,E)=\mathrm{\Lambda }\left(E\right)X_{\lambda }(x,E),$$(17)

where

$$\mathrm{\Lambda }\left(E\right)=K{\left(\frac{E}{E_{\mathrm{inj}}}\right)}^{{\beta }_{\lambda }}\exp \left[{\left(\frac{E}{E_{\mathrm{b}}}\right)}^{{\delta }_{\lambda }}\right]$$(18)

gives the mean free path at the shock as a function of energy E, and K, β_λ, E_b and δ_λ are the fitting parameters, and

$$X_{\lambda }(x,E)={\left[1+{\left(\frac{x}{\Delta x}\right)}^q\right]}^{\frac{1}{q}}$$(19)

characterizes the spatial extent of the foreshock, and Δx(E) and q(E) are fitting parameters, which are functions of energy. Figure 3 shows example of the fits.

Figure 3 Examples of the fits of the spatial (left panel) and energy (right panel) dependencies of the mean free path in the simulation run 6. The former is obtained at E = 5.27 MeV. The crosses indicate points that were not considered during the fitting process.

Equation (18) represents a modification of the energy dependence of the Bell’s mean free path at the shock (see Appendix B). We see that the applied fitting function Λ(E) describes well the lower-energy part of the mean free path at the shock, which corresponds to the quasi-steady-state power-law part of the particle energy spectrum at the shock. In turn, equation (19) presents a modification of the spatial dependence of the Bell’s mean free path. The introduced parameter q < 1 reflects the pitch-angle dependence in the quasi-linear resonance condition of particle scattering. Interestingly, at x ≫ Δx the mean free path is a linear function of x, identical to Bell’s theory.

Using equations (17) and (19), the particle intensity spectrum at the matching point x_M in SOLPACS can be calculated as

$$\begin{array}{ll}{j}_{\mathrm{M}}^{\left(S\right)}& =j_{\mathrm{sh}}^{\left(S\right)}\exp \left\{-\frac{3V_{\mathrm{A}}(M_{\mathrm{A}}-1)}{v}{\int }_0^{x_{\mathrm{M}}}\mathrm{}\frac{\mathrm{d}x}{{\lambda }^{\left(S\right)}}\right\}\\ & =j_{\mathrm{sh}}^{\left(S\right)}\exp \left\{-\frac{3V_{\mathrm{A}}(M_{\mathrm{A}}-1)}{v\mathrm{\Lambda }}{ }_2{F}_1\left[\frac{1}{q};\frac{1}{q};1+\frac{1}{q};-{\left(\frac{x_{\mathrm{M}}}{\Delta x}\right)}^q\right] x_{\mathrm{M}}\right\},\end{array}$$

where ₂F₁ [·] is the hypergeometric function (e.g., Abramowitz & Stegun, 1964). The particle spectrum at the shock to be used in PARADISE is

$$j_{\mathrm{sh}}^{\left(P\right)}=j_{\mathrm{M}}^{\left(S\right)}\exp \left\{\frac{3V_{\mathrm{A}}(M_{\mathrm{A}}-1)}{v}{\int }_0^{x_{\mathrm{M}}}\mathrm{}\frac{\mathrm{d}x}{{\lambda }^{\left(P\right)}}\right\},$$(20)

where

$${\lambda }^{\left(P\right)}={\lambda }_{1 \mathrm{au}}^{\left(S\right)}={\mathrm{\Lambda }}_0{\left[1+{\left(\frac{x}{\Delta x_0}\right)}^{q_0}\right]}^{\frac{1}{q_0}},$$(21)

the subscript “0” indicates that a quantity corresponds to the shock at 1 au. Substituting, we get

$$j_{\mathrm{sh}}^{\left(P\right)}=j_{\mathrm{M}}^{\left(S\right)}\exp \left\{\frac{3V_{\mathrm{A}}(M_{\mathrm{A}}-1)}{v{\mathrm{\Lambda }}_0}{ }_2{F}_1\left[\frac{1}{q_0};\frac{1}{q_0};1+\frac{1}{q_0};-{\left(\frac{x_M}{\Delta x_0}\right)}^{q_0}\right] x_{\mathrm{M}}\right\}.$$(22)

2.2.3 Particle energy spectrum at the shock

The particle energy spectrum at the shock simulated in SOLPACS can be reproduced by using the following functional form (Afanasiev et al., 2015):

$$j_{\mathrm{sh}}^{\left(S\right)}\left(E\right)=J{\left(\frac{E}{E_{\mathrm{inj}}}\right)}^{-\beta }\exp \left[-{\left(\frac{E}{E_{\mathrm{c}}}\right)}^{\delta }\right],$$(23)

where the particle injection energy E_inj, the amplitude factor J, and the cutoff (i.e., roll-over) energy E_c are parameterized as follows:

$$E_{\mathrm{inj}}=\frac{1}{2}{m}_{\mathrm{p}}{M}_{\mathrm{A}}^2{V}_{\mathrm{A}}^2{\left[\left(3-\frac{2r_{\mathrm{mag}}}{r}\right)\left(1-{\mu }_{\mathrm{c}}\right)+{\alpha }_{\mathrm{inj}}{\mu }_{\mathrm{c}}\right]}^2,$$(24)

where ${\mu }_{\mathrm{c}}={\left(1-r_{\mathrm{mag}}^{-1}\right)}^{1/2}$ is the critical value of the pitch-angle cosine, r_mag is the magnetic compression ratio of the shock, α_inj = 4 is the correction factor,

$$J=\frac{\beta +1}{\pi \sqrt{8m_{\mathrm{p}}{E}_{\mathrm{inj}}}}{ϵ}_{\mathrm{inj}}{n}_{\mathrm{p}},$$(25)

$$E_{\mathrm{c}}=8.25\cdot 10^{-2} {\left(\frac{E_{\mathrm{c},\mathrm{th}}}{E_{\mathrm{ref}}}\right)}^{0.86} \left[\mathrm{MeV}\right],$$(26)

where

$$E_{\mathrm{c},\mathrm{th}}=E_{\mathrm{inj}}{\left[1+2\left({\beta }_{\lambda }+\frac{1}{2}\right)M_{\mathrm{A}}^2\frac{r_{\mathrm{sc}}(r_{\mathrm{sc}}-1)}{r^2}\frac{V_{\mathrm{A}}}{v_{\mathrm{inj}}}\frac{V_{\mathrm{A}}\Delta t}{K}\right]}^{1/({\beta }_{\lambda }+\frac{1}{2})},$$(27)

and E_ref = 1 MeV.

Equation (24) is obtained by an ad-hoc modification of the expression for the injection momentum, p_inj = 2m_pM_AV_A (1 − r_magr⁻¹), used by Vainio et al. (2014), whereas Equation (25) is derived in the same manner as equation (29) in Vainio et al. (2014).

Equation (26) is obtained by using an approach similar to Vainio et al. (2014). Namely, we first integrate the analytical expression for the particle acceleration rate (Drury, 1983) over time, assuming a purely power-law dependence for the mean free path at the shock, i.e., $\mathrm{\Lambda }\left(E\right)=K(E/E_{\mathrm{inj}}{)}^{{\beta }_{\lambda }}$ (cf. Eq. (15) in Vainio et al., 2014). This results in equation (27). Then we fit the simulated energy spectra, using equation (23) as a fitting function with $\tilde{J}$ ($=J{E}_{\mathrm{inj}}^{\beta }$), β, E_c and δ treated as independent fitting parameters. Finally, we compare the theoretical values of the cutoff energy E_c,th with the fitted ones and fit a power-law function to the resulting dependence (see the left panel of Fig. 4), which results in equation (26).

Figure 4 Left panel: Cutoff energy obtained in Runs 1–7 versus the cutoff energy calculated based on equation (27). The red line is a power-law fit given by equation (26). Right panel: Example of the simulated particle spectrum at the shock (blue line) and the reproduced spectrum based on the analytical model (orange line). The case presented is run 14.

In equation (27), Δt is time equivalent to the simulation time in SOLPACS. We note, however, that when these equations are implemented in PARADISE, it becomes a free parameter that has to be constrained based on test simulations. The right panel of Figure 4 shows, as an example, a simulated particle energy spectrum at the shock and the reproduced spectrum obtained based on equations (23)–(27).

2.2.4 Obtaining semi-analytical functions for the parameters of the mean free path fitting function

To obtain analytical representations for the fitting parameters K, β_λ, E_b and δ_λ of the function Λ(E) given by equation (18), we studied their evolution with time in several simulation runs. The parameter values were obtained by fitting the energy dependence of the mean free path at the shock at the output times. Figure 5 shows the typical temporal evolution of the parameters in a SOLPACS simulation. One can see that all parameters but one, δ_λ, exhibit clear systematic changes during a simulation.

Figure 5 Example of temporal evolution of the parameters of the fitting function Λ(E) in a SOLPACS simulation. The plots show each parameter versus a simulation snapshot number. The case presented is run 13.

Clearly, the evolution of the amplitude factor K and the power-law index β_λ is due to the contribution of particles being accelerated to progressively higher energies, which then generate waves that resonate with lower-energy particles. Figure 5 shows that, if a simulation runs long enough, the temporal evolution of these parameters slows down. This is expected because there are progressively smaller number of higher-energy particles, modifying the mean free paths of lower-energy particles. Therefore, as the first step, we neglect this temporal evolution at later times and treat these parameters as constants. Specifically, for β_λ we take a value of 0.2. This value is motivated by our earlier finding that the spectrum of self-generated Alfvén waves in SOLPACS simulations approaches asymptotically the k⁻² form as a function of wavenumber k (Afanasiev et al., 2015). Taking into account that in the quasi-linear theory the energy dependence of the mean free path is given as λ(E) ∝ E^1 − q/2, where q (>0) is the power-law index of the wave spectrum, we conclude that the power-law-like part of the energy dependence of the mean free path at the shock should be weak, at least weaker than for the initial k^−3/2 spectrum, i.e., ∝E^1/4 (cf. Fig. 7 in Afanasiev et al., 2015).

Based on Bell’s theory (see Eqs. (B2) and (B4) in Appendix B), one would expect the amplitude factor K to scale with V_A(ϵ_injΩ₀)⁻¹, where Ω₀ is the proton cyclotron frequency. We verified this using the full set of simulation runs and found that this relationship holds true (Fig. 6). Fitting the dependence with a linear function gives:

$$K=0.22 \frac{V_{\mathrm{A}}}{{ϵ}_{\mathrm{inj}}{\mathrm{\Omega }}_0}.$$(28)

Figure 6 Fit to the amplitude factor K of the mean free path Λ at the shock versus VA/(ϵinjΩ0).

Regarding the parameter δ_λ, we fix it at 1.4. It should be noted that δ_λ in fact varies at least with ϵ_inj (see, e.g., Fig. 7 in Afanasiev et al., 2015). We will leave this for a detailed investigation in a future study.

The spatial dependence of the mean free path in the upstream region is fitted using equation (19), in which Δx and q are fitting parameters. These parameters are energy-dependent. We found that the following functions provide good fits for the obtained Δx and q as functions of energy:

$$\Delta x\left(E\right)=C_1\exp {\left(\frac{E}{\Delta E_1}\right)}^{{\alpha }_1},$$(29)

where C₁, ΔE₁, and α₁ are fitting parameters, and

$$q\left(E\right)=C_2\exp {\left(\frac{E}{\Delta E_2}\right)}^{{\alpha }_2},$$(30)

where C₂, ΔE₂, and α₂ are fitting parameters as well. For the two sets of fitting parameters, we have obtained the following representations:

$$C_1=0.35{\left(\frac{M_{\mathrm{A}}}{{ϵ}_{\mathrm{inj}}{x}_{\mathrm{ref}}}\frac{V_{\mathrm{A}}}{{\mathrm{\Omega }}_0}\right)}^{0.71} \left[R_{\odot }\right],$$(31)

$${\alpha }_1=0.83-0.02\frac{E_{\mathrm{b}}}{E_{\mathrm{ref}}},$$(32)

$$\Delta E_1=1.64{\left(\frac{E_{\mathrm{b}}}{E_{\mathrm{ref}}}\right)}^{1.07} \left[\mathrm{MeV}\right],$$(33)

$$C_2=0.66,$$(34)

$${\alpha }_2=0.47,$$(35)

$$\Delta E_2=9.45{\left(\frac{E_{\mathrm{b}}}{E_{\mathrm{ref}}}\right)}^{0.82} \left[\mathrm{MeV}\right],$$(36)

where x_ref = 1 R_⊙ and E_ref = 1 MeV are the reference values.

The idea underlying the derivation of equations (31)–(36) for the fitting parameters of Δx(E) and q(E) is to make them related to the plasma and shock parameters (and to the model input parameters). To obtain equation (31) for the fitting parameter C₁, which characterizes the foreshock spatial extent Δx, we used Bell’s theory as a starting point since the mean free path at E ≪ E_c is in a quasi-steady state by the end of all the SOLPACS simulations used in this work. In Bell’s theory, C₁ ∝ V_A/(σϵ_injΩ₀) (see Appendix B). This scaling, however, does not work well with the values of C₁ obtained from fitting Δx(E) in different SOLPACS simulations. We found that including the Alfveńic Mach number into the Bell’s scaling organizes the fitted values of C₁ into a monotonic dependence that is quite close to a linear dependence (Fig. 7, left panel).

Figure 7 Fits to the parameters of the energy dependence of the foreshock spatial scale Δx.

It is more difficult to find a suitable combination of plasma and shock parameters for C₂ in Eq. (30) that would organize its values into a monotonic dependence, because q does not exist in Bell’s theory. Figure 8 (left panel) shows, as an example, the values of C₂ versus cosθ_Bn/(ϵ_injαM_A). The values, however, vary only within ∼20% of their average. Therefore, for C₂, we take its average over all the simulations performed in this work.

Figure 8 Fits to the parameters of the energy dependence of the foreshock spatial parameter q.

The parameters ΔE₁, α₁ and ΔE₂ were found to scale with the roll-over energy E_b of Λ(E) (Figs. 7 and 8, central and right panels). We used power-law functions to fit ΔE₁ and ΔE₂ versus E_b and a linear function to fit α₁ versus E_b. The parameter α₂ was found to be approximately the same in all the simulation runs (Fig. 8, central panel).

2.3 Implementation in PARADISE

The PARADISE model solves the focused transport equation (e.g., van den Berg et al., 2020), that is,

$$\frac{\partial j}{\partial t}+\nabla \cdot \left(\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}j\right)+\frac{\partial }{\partial \mu }\left(\frac{\mathrm{d}\mu }{\mathrm{d}t}j\right)+\frac{\partial }{\partial p}\left(\frac{\mathrm{d}p}{\mathrm{d}t}j\right)=\frac{\partial }{\partial \mu }\left(D_{\mu \mu }\frac{\partial j}{\partial \mu }\right)+S,$$(37)

with

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}={\mathbf{V}}_{\mathrm{s}w}+\mu v\widehat{\mathbf{b}},$$(38)

$$\frac{\mathrm{d}\mu }{\mathrm{d}t}=\frac{1-{\mu }^2}{2}\left(v\nabla \cdot \widehat{\mathbf{b}}+\mu \nabla \cdot {\mathbf{V}}_{\mathrm{s}w}-3\mu \widehat{\mathbf{b}}\widehat{\mathbf{b}}:\nabla {\mathbf{V}}_{\mathrm{s}w}-\frac{2}{v}\widehat{\mathbf{b}}\cdot \frac{\mathrm{d}{\mathbf{V}}_{\mathrm{s}w}}{\mathrm{d}t}\right),$$(39)

$$\frac{\mathrm{d}p}{\mathrm{d}t}=\left(\frac{1-3{\mu }^2}{2}(\widehat{\mathbf{b}}\widehat{\mathbf{b}}:\nabla {\mathbf{V}}_{\mathrm{s}w})-\frac{1-{\mu }^2}{2}\nabla \cdot {\mathbf{V}}_{\mathrm{s}w}-\frac{\mu }{v}\widehat{\mathbf{b}}\cdot \frac{\mathrm{d}{\mathbf{V}}_{\mathrm{s}w}}{\mathrm{d}t}\right)p,$$(40)

where j denotes the differential particle intensity, $\widehat{\mathbf{b}}$ the unit vector in the magnetic field direction, and V_sw is the solar wind velocity vector. These two vector fields are taken from simulations preformed with the time-dependent MHD model EUHFORIA (Pomoell & Poedts, 2018). The source function S is nonzero exclusively at the shock surface, as defined in equation (1). The value of S in each point of the shock is obtained by substituting equations (22) and (6). Moreover, for the PARADISE simulations presented in this work, we only incorporate a pitch-angle diffusion process, as described by D_μμ in equation (37), while excluding perpendicular diffusion (e.g., Wijsen et al., 2019) and guiding center drifts (e.g., Wijsen et al., 2020). The impact of these processes on the PARASOL simulations will be explored in future investigations.

PARADISE uses the results of quasi-linear theory (QLT; Jokipii, 1966) to prescribe the following pitch-angle diffusion coefficient (see Agueda & Vainio, 2013; Wijsen et al., 2019):

$$D_{\mu \mu }=D_0\left(\frac{\left|\mu \right|}{1+\left|\mu \right|}+ϵ\right)\left(1-{\mu }^2\right),$$(41)

where μ denotes the cosine of the pitch-angle, ϵ = 0.048 is a parameter bridging the resonance gap at μ = 0 (e.g., Klimas & Sandri, 1971). The parameter D₀ is determined by prescribing the particles’ parallel mean free path λ_||, which relates to D_μμ through (Hasselmann & Wibberenz, 1970)

$${\lambda }_{\parallel }=\frac{3v}{8}{\int }_{-1}^1\mathrm{}\frac{{\left(1-{\mu }^2\right)}^2}{D_{\mu \mu }}\mathrm{d}\mu .$$(42)

Upstream of the shock wave, we choose the parallel mean free path

$${\lambda }_{\parallel }=\min \left[{\lambda }^0,{\lambda }^{\left(P\right)}\right],$$(43)

where λ^(P) is defined in equation (21), and ${\lambda }^0=(0.1 \mathrm{au}){\left(R/R_{\mathrm{ref}}\right)}^{2-q_0}$. Here, q₀ = 5/3 represents the spectral index of a Kolmogorov turbulence spectrum, and R is the particle rigidity, with R_ref set at 43 MV, which corresponds to the rigidity of a 1 MeV proton.

Downstream of the shock, we choose the parallel mean free path as

$${\lambda }_{\parallel }=\min \left[{\lambda }^0,{\mathrm{\Lambda }}_0{e}^{\left(-d/d_0\right)}\right],$$(44)

where d represents the distance from the shock wave, and we set d₀ = 1 R_⊙. To summarize, equations (43) and (44) provide the particles with a constant parallel mean free path λ⁰ when they are located far from the shock. As a particle approaches the shock from the upstream side, it enters a foreshock region where the parallel mean free path, now given by λ^(P), gradually decreases based on the fittings to the SOLPACS simulations. On the downstream side of the shock, λ_|| exhibits an exponential increase towards λ⁰.

3 The SEP event of July 12, 2012

3.1 Observations

To evaluate the performance of our new model, PARASOL, we conducted simulations of an SEP event observed in situ. Since our solar wind model, EUHFORIA, only starts at 21.5 R_s, this paper will primarily focus on the energetic storm particle (ESP) component of the SEP event. Specifically, we simulated the SEP event associated with a halo CME that occurred on July 12, 2012, and was observed in situ by near-Earth spacecraft. This CME originated from NOAA Active Region 11520, located near the center of the solar disk at the time of the eruption.

In Figure 9, panel (1a) displays proton intensities in various ion energy channels, as recorded by the Low-Energy Magnetic Spectrometer (LEMS120) of the Electron, Proton, and Alpha Monitor (EPAM; Gold et al., 1998) aboard the Advanced Composition Explorer (ACE) spacecraft. The panel also includes energy channels from the Solar Energetic Particle Environment Modelling Project (SEPEM) reference data set (RDS) version 2.0 (Jiggens et al., 2018). The figure illustrates that this SEP event was characterized by an energetic storm particle (ESP) event, which commenced on July 14, 2012, and lasted for approximately a day. The particle intensity peak in the low energy channels coincides with the arrival of the CME-driven shock at Earth, as indicated by the solid line in the figure.

Figure 9 Comparison of observed and simulated proton intensities and solar wind plasma parameters for the 12 July 2012 SEP event. The panels correspond to (a) omnidirectional particle intensities, (b) IMF magnitude, (c) solar wind speed and (d) number density. The left figure displays in-situ observational data obtained from the ACE and WIND spacecraft (refer to the text for further details). In contrast, the central and right figures showcase simulation results for Earth’s position and the W25 observer. The vertical lines indicate the arrival of the observed and simulated shock waves at the corresponding observers.

Figure 9a, panel 1b, displays the interplanetary magnetic field (IMF) magnitude measured by the Magnetic Field Investigation (MFI) on Wind (Lepping et al., 1995), while panels (1c) and (1d) show, respectively, the solar wind speed and proton number density, as measured by the Solar Wind Experiment (SWE) on Wind (Ogilvie et al., 1995). The shock arrival is clearly evident as a sudden jump in all the solar wind variables at 17:26 UT on July 14.

3.2 EUHFORIA simulation

The CME responsible for this SEP event was previously simulated using the linear force-free spheromak model within EUHFORIA, with parameters derived from remote-sensing observations (Scolini et al., 2019). Additionally, the ESP event was recently modeled by Wijsen et al. (2022) using the PARADISE and EUHFORIA model. In this study, we will utilize the same EUHFORIA simulation as presented in Wijsen et al. (2022) to describe the solar wind conditions and the propagation of the CME. For an in-depth analysis of the CME and SEP event, as well as the EUHFORIA simulation, we refer to these two papers. To recap the essential parameters, the spheromak CME was introduced into the simulation on July 12, 2012, at 19:24 UT, originating from the inner boundary of EUHFORIA’s computational domain, positioned at a heliocentric distance of 0.1 au (21.5 R_s). The CME was directed towards Earth and the insertion speed and radius of the spheromak were assumed to be 763 km s⁻¹ and 16.8 R_⊙, respectively. A snapshot of the EUHFORIA simulation can be seen in Figure 3g of Wijsen et al. (2022).

3.3 PARASOL simulation

The PARASOL model comprises three free parameters: the particle injection efficiency ϵ_inj, the timescale Δt, and the matching distance x_m. As an illustrative example, we present a PARASOL simulation with ϵ_inj = 5 × 10⁻⁴, Δt = 10 h, and x_m = 10 R_s. A comprehensive parametric investigation into the impact of these three parameters will be detailed in a forthcoming study.

The omnidirectional time-intensity profiles, as modeled by PARASOL, are presented in panels 2a and 3a of Figure 9 for Earth’s location and 25° west of Earth, referred to as the W25 observer, respectively. Panels 2b–2d and 3b–3d display the time profiles of the EUHFORIA solar wind at the locations of these two spacecraft. Observing panel 2a, we note that PARASOL simulates an ESP event. Similar to the observations, the intensity profiles of the lower energy channels peak at the arrival of the modeled CME-driven shock. Moreover, when comparing panels 1a and 2a, it becomes evident that PARASOL slightly underestimates the peak intensity levels in the lowest energy channels and overestimates the intensities in the highest energy channels. In other words, the modeled energy spectrum appears to be harder than the observed one. This is illustrated in Figure 10a, where the observed and simulated energy spectra at Earth are shown in blue and orange, respectively. The figure highlights that, although PARASOL underestimates intensity at lower energies, the observed and simulated spectra exhibit similar slopes.

Figure 10 Energy spectra of observed and simulated proton intensities for the 12 July 2012 ESP event. Panel (a) presents the energy spectrum at the shock arrival, while panel (b) displays it eight hours before the shock arrival. Solid lines represent simulated intensities, with the orange line indicating Earth and the green line representing the W25 observer. Dashed lines correspond to in-situ data from ACE and SEPEM, as also depicted in panel 1a of Figure 9.

Figure 9 also shows that the onset of the SEP event is not accurately reproduced. This discrepancy is a consequence of the spheromak CME not being wide enough to adequately capture the shock wings and the IMF being significantly more radial than the modeled IMF, as elaborated in Wijsen et al. (2022). Consequently, in the simulation Earth establishes a magnetic connection with the shock only on July 13, around 22:00 UT. In contrast, the observed onset of the SEP event suggests that Earth likely had a direct magnetic field connection to the shock wave shortly after the CME eruption on July 12, around 17:00 UT.

Due to its more westward position relative to Earth, the W25 observer establishes a magnetic connection to the Earth-directed CME soon after it is inserted into the EUHFORIA simulation at 0.1 au. As a result, the modelled onset of the SEP event for the W25 observer is closer to the observed one. It is still slightly delayed, however, since the computational domain of EUHFORIA, and hence PARASOL, does not include the solar corona.

Figure 10a displays the energy spectrum for the W25 observer at the time of the (simulated) shock arrival. At lower energies, the spectrum closely resembles that of the simulated Earth observer, while at higher energies, the spectrum softens, similar to the observed energy spectrum.

As previously mentioned, the intensities in different energy channels are sensitive to the choice of the injection efficiency parameter, ϵ_inj. This sensitivity is demonstrated in Figure 11, which showcases the particle intensities for the W25 observer in a PARASOL simulation with ϵ_inj = 10⁻⁴. A direct comparison between Figure 11 and panel 3a of Figure 9 reveals that a reduced ϵ_inj results in lower particle intensities across all energy channels and a less pronounced ESP event.

Figure 11 Modeled omnidirectional particle intensities for the W25 observer in a PARASOL simulation with ϵinj = 10−4.

To address this sensitivity and determine the optimal value for ϵ_inj, a future study will involve modeling several different SEP events. This approach will help quantify the most suitable choice for the injection efficiency parameter. In the remainder of this section, we will focus our analysis on the PARASOL simulation with ϵ_inj =5 × 10⁻⁴.

Figure 12 shows the omnidirectional proton intensities for two different energy channels within the plane of constant latitude encompassing Earth, 27 h after the insertion of the CME in the simulation domain. Figure 12 offers valuable insights into the spatial distribution of proton intensities resulting from the shock’s interaction with the surrounding solar wind. In comparing the two panels, we observe that low-energy particles are more tightly confined near the CME than high-energy particles, due to their lower speed (relative to the CME) and their smaller mean free path in the foreshock region, as further discussed in the next section. Furthermore, the right panel reveals that the highest intensities in the high-energy range are concentrated at the CME’s center and west flank, whereas in the low-energy range (left panel), the highest intensities are concentrated along the CME’s east flank.

Figure 12 The modeled omnidirectional proton intensities within the plane of constant latitude encompassing Earth, 27 h after the insertion of the CME in the simulation domain. The left panel displays proton intensities within the 310–580 keV range, while the right panel presents data for the 15–31 MeV range. The white line denotes the position of the CME-driven shock wave. Earth and the W25 observer are marked by green and blue dots, respectively, and their magnetic connectivity to the shock wave is also illustrated.

These energy-dependent variations of the SEP intensities along the shock surface can be better understood by examining the properties of the shock, as illustrated in Figure 13, which corresponds to the same date and time as the intensity snapshots in Figure 12. Figures 13a and 13b depict the shock angle θ_Bn and the Alfvenic Mach number M_A, respectively. One key observation is that the eastern flank of the CME exhibits quasi-parallel shock geometry (θ_Bn < 45°), while the center and the western flank of the CME display quasi-perpendicular shock geometry (θ_Bn > 45°). This spatial variation in the shock geometry is a direct consequence of the spiral geometry of the upstream IMF. Additionally, the Alfvenic Mach number M_A reaches its maximum value at the center and the western flank of the CME shock. This is attributed partly to the shock being quasi-perpendicular in these regions (recall that M_A is measured in the HT-frame) and partly due to variations in the non-uniform upstream density and solar wind speed (not shown in the figure).

Figure 13 The CME-driven shock, observed 27 h after its introduction into the simulation domain, is depicted. The shock surface is color-coded to represent various parameters: the shock obliquity θBn (a); the Mach Alfvén number MA (b); the intensity J and defined in equation (23) (c); and the cutoff energy Ec, as defined in equation (26) (d). Green and blue dots indicate Earth and the W25 observer, respectively, and their magnetic connections to the shock wave are also illustrated. The white lines on the shock indicate the location of the solar equatorial plane and the meridional slice containing Earth.

These spatial variations translate into different levels of particle emission and energy spectra in the PARASOL model at the shock wave. To illustrate this, we present the intensity parameter J and the cutoff energy E_c in Figures 13c and 13d. As indicated in equation (23), J represents the energy-independent part of the particle intensity at the shock wave. Figure 13c clearly illustrates that the highest values for J are concentrated in the east flank of the shock wave. However, the same eastern flank exhibits relatively lower values for the cutoff energy E_c, compared to the nose and the west flank of the CME (see Fig. 13d).

In summary, the varying shock conditions along the shock surface explain why the eastern shock flank generates the highest intensities in the low-energy channels, while the west flank produces the highest intensities in the high-energy channels. The intricate interplay of shock geometry, Alfvenic Mach number, and upstream conditions plays thus a pivotal role in shaping the energy-dependent patterns in the distribution of proton intensities.

4 Discussion

A test PARASOL simulation of the July 12, 2012 SEP event using the EUHFORIA MHD simulation of the solar wind and the CME have reproduced the observed ESP event (E ≲ 5 MeV) in the close vicinity of the shock reasonably well (within one order of magnitude in intensity). The overestimated intensities at higher energies are due to an overestimated cutoff energy at 1 au in the simulation, which is governed by equation (27) and largely results from the assumed Δt = 10 h. It should be noted that this parameter is not necessarily constant and can be a function of heliocentric distance.

In the particle acceleration model where the shock upstream turbulence is self-generated, one of the possibilities is that the cutoff energy is ultimately governed by the flux of protons injected into the acceleration process. In this case, a decrease in the plasma density in the interplanetary medium should lead to a decrease in the injected proton flux (assuming that ϵ_inj does not change much). This should lower the cutoff energy while the shock propagates further away from the Sun, if the maximum cutoff energy was achieved when the shock was in the corona. In this case, to constrain the parameter Δt as a function of distance from the Sun one needs to use a SOLPACS simulation in a global setup, i.e., taking into consideration heliospheric variations of the solar wind plasma density and the magnetic field. Another option (perhaps more applicable for strong shocks) is that the cutoff energy is controlled by the particle escape from the foreshock due to magnetic focusing (Vainio et al., 2014). In this case, it can be shown that Δt does not enter the equations determining the cutoff energy. This will be addressed in a future study.

As we mentioned above, the ESP event under consideration was modeled earlier by Wijsen et al. (2022), in which PARADISE was used not only to propagate but also to accelerate particles (50 keV seed protons) at the EUHFORIA-modeled CME shock. The study employed in-situ observations of the ESP event to estimate the spatial variation of the mean free path in the foreshock region. In contrast, the new PARASOL model does not rely on SEP observations for the event, making it more suitable for forecasting purposes.

The energy spectrum at the shock in Wijsen et al. (2022) showed a strong agreement with observations up to 1 MeV. However, above 1 MeV, the spectrum was too soft. The authors attributed this softening to low SEP acceleration efficiency at the modeled MHD shock as a result of the finite resolution of the MHD simulation leading to an overly thick shock. With PARASOL, instead of modeling particle acceleration with PARADISE, we inject an accelerated particle population based on MHD shock properties and the SOLPACS-based model. This approach circumvents the issue of shock thickness due to finite resolution, but makes the cutoff energy to depend on a free parameter (Δt) of the model.

A striking discrepancy between the observed and simulated time-intensity profiles in Figure 9, not related to the EUHFORIA simulation of the solar wind and CME, can be seen in the low energy channels (up to ∼0.5 MeV) before the shock arrival. The observed time-intensity profiles at these energies nearly overlap, which corresponds to a nearly flat energy spectrum of particles, starting from some distance from the shock in the upstream region. This feature is further illustrated in Figure 10b, which shows the energy spectra eight hours before the shock arrival, at the onset of the ESP event. Here, the observed energy spectrum is particularly hard (nearly constant) below 1 MeV. In contrast, in the simulations, intensities at low energies decrease rapidly with an increasing distance from the shock in the upstream region, resulting in an “inverted” energy spectrum (in which intensity increases with energy) at a given location in the upstream region. The formation of such a spectrum in the PARASOL simulation stems from the SOLPACS self-consistent simulations. Namely, in SOLPACS, higher-energy protons are able to generate Alfvén waves in the upstream region at frequencies resonant with lower-energy protons, thus leading to a decrease of the scattering mean free path at the lower energies. This leads to an increasingly efficient trapping of lower-energy protons close to the shock, and may result (depending on how efficiently the waves are produced) in an “inverted” energy spectrum of protons in the upstream region.

Overlapping particle intensities upstream of ESP events is not uncommon, as recently pointed out by Lario et al. (2018). However, the mechanism is responsible for such flat, rather than inverted, energy spectrum at low-energy particle intensities is currently not understood, but some ideas have been proposed. Those include a balance between the particle acceleration at the shock and adiabatic cooling (Prinsloo et al., 2019) and a velocity filter mechanism (Perri et al., 2022), which suggests that the lower the particle energy, the smaller the fraction of such particles that can escape from the shock, assuming an isotropic distribution of particles. Clearly, the inverted spectrum situation obtained in our simulations is a result of too efficient trapping of particles near the shock due to too small mean free paths of low-energy particles. If the mean free paths at low energies were larger, the intensities of particles at those energies would be larger (the “inverted” part of the spectrum would be flatter). Importantly, the mean free path in SOLPACS (and PARASOL) is an increasing function of particle energy.

Earlier, Wijsen et al. (2022) calculated the mean free paths in this ESP event from the observed particle intensities and found that they increase with energy only in a close vicinity of the shock (≲1.7 R_⊙ along the normal to the shock), where the intensities are not overlapped, but further away, where the intensities become overlapped, the energy dependence of the mean free path turns over. Moreover, the mean free paths at the lowest energies further away from the shock are substantially (about one order of magnitude) larger than close to the shock.

The rigidity dependence of the mean free path in the quasi-linear theory is determined by the power-law index q of the turbulence spectrum I(k), λ ∝ P^2−q (e.g., Bieber et al., 1994). In SOLPACS, the spectrum of self-generated Alfvén waves is characterized by power-law indices q ≲ 2, which explains the resulting energy dependence of the mean free path. However, a steeper turbulence spectrum (q > 2) would result in an opposite energy dependence of the mean free path (the mean free path decreases with particle energy). One candidate process, not accounted for in the SOLPACS model, that can produce such a steeper spectrum at the high end of the MHD frequency range is the Alfvén wave damping on thermal protons (e.g., Berezhko & Taneev, 2016). The spectrum steepening, however, produces a resonance gap (Smith et al., 1990) and makes it more difficult for particles to scatter across the pitch angle of 90°, which should lead toward reduction of the particle acceleration efficiency at the shock. The effect of resonance broadening is usually invoked to deal with this problem, which allows a particle to interact not only with the wave spectral component strictly determined by the quasi-linear resonance condition (e.g., Vainio, 2003), but with spectral components from a wider range of wavenumbers (see, e.g., Bieber et al., 1994; Ng & Reames, 1995, and references therein).

Ng & Reames (1995) proposed a theory where the resonance broadening was attributed to large-amplitude medium-scale fluctuations of the magnetic field, non-resonant with the particles at the considered energies. Bearing that theory in mind, we hypothesize that the resonance broadening can be very efficient close to the shock (stronger large-amplitude magnetic field variations), so that it allows large pitch-angle particles to interact with lower-frequency components of the self-generated Alfvén wave spectrum not affected significantly by the damping process. It can be shown that this can lead to a mean free path increasing with energy. Further away from the shock, the resonance broadening becomes less efficient, so the effect of a steeper spectrum due to the wave damping becomes more important, thus leading to the opposite dependence of the mean free path on energy. Also, one can expect in this case that the mean free path values, especially at lower energies, should be larger in magnitude further away from the shock than close to the shock. This scenario qualitatively agrees with the theory of Ng & Reames (1995, see their Fig. 9).

It is also interesting to note that in some events where the presence of the flat spectrum was reported, the spectra are in fact slightly inverted than flat (see, e.g., Fig. 4 in Ng et al., 2012). Besides, the significantly inverted particle spectrum in a gradual SEP event has been recently reported by Lario et al. (2021) and Ding et al. (2024). In the latter work, such a spectrum was attributed to the presence of a pre-existing enhanced turbulence in the ambient solar wind.

It would be interesting to test the described scenario of the formation of overlapped intensities (nearly flat spectra) in the shock upstream region. This can be done by including the effect of Alfvén wave damping on thermal protons to the SOLPACS model and considering the resonance broadening as decreasing with the distance from the shock (towards upstream). A successful reproduction of this phenomenon in SOLPACS simulation would allow us to fine-tune our semi-analytical model of the inner foreshock in order to model (and eventually predict) this behavior of particle intensities at lower energies.

5 Conclusions and outlook

In this paper, we introduced a new physics-based simulation model, PARASOL, for simulating proton acceleration in a CME-driven shock and their transport in the interplanetary medium. PARASOL combines a semi-analytical model of the inner foreshock region, derived from self-consistent simulations of particle acceleration using the SOLPACS model, with the test-particle simulation model, PARADISE, which handles particle transport in the outer foreshock and the ambient solar wind, where wave growth is not significant. PARASOL requires a MHD model of the solar wind and the shock in order to obtain the necessary plasma and shock parameters. In particular, the semi-analytical model of the foreshock requires the ambient plasma density, n, the Alfvén speed, V_A, the shock-normal angle, θ_Bn, all at the shock location, and the Alfvénic Mach number (in the HT frame) M_A of the shock. For a test PARASOL simulation, we used a EUHFORIA simulation of the CME associated in July 12, 2012 SEP event.

To constrain better the free parameters of PARASOL, the focus of a future modeling study should be on ESP events having the so-called classic component, i.e., a slow increase of the particle intensity starting a few hours before the shock arrival (Lario et al., 2003; Giacalone, 2012). Such a rise of intensity is consistent with the presence of the foreshock.

The ongoing development of PARASOL aims to create an efficient simulation model for operational SEP forecasting that also provides an accurate description of the physical processes involved in particle acceleration, particularly Alfvén wave growth upstream of the shock. The first version of the model presented here focuses mainly on the ESP component of an SEP event. Currently, the entire simulation chain takes around 1500 CPU hours. On a supercomputer, this allows us to produce a forecast within roughly 2 h, which is adequate for predicting an ESP event in a timely manner. However, there are several potential optimizations that could significantly reduce both the computational and prediction times. For instance, EUHFORIA and PARASOL are currently run sequentially, but they could be executed in parallel to save time. Additionally, in the simulations presented here, particles were injected along the entire shock wave, which adds considerable computational cost. In a near-Earth forecasting scenario, it would be sufficient to inject particles only in the vicinity of a shock location that is magnetically connected to Earth, thus reducing the necessary CPU hours drastically.

To address a full SEP event, such as capturing the event onset, an MHD model of the propagating shock in the corona is required. We are working on such a simulation framework by combining EUHFORIA with the recently developed coronal model COolfluid COroNal UnsTructured (COCONUT) (Perri et al., 2022), which can perform self-consistent flux-rope simulations (see e.g., Linan et al., 2023). This integrated approach promises to enhance the accuracy and reliability of SEP event predictions significantly.

Acknowledgments

This research has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreements No 870405 (EUHFORIA 2.0). A.A. and R.V. also acknowledge funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101004159 (SERPENTINE) and from the Finnish Centre of Excellence in Research of Sustainable Space (FORESAIL) funded by the Research Council of Finland (grant no. 352847). This study has been performed in the framework of the European Union’s Horizon 2020 research and innovation programme (EU/H2020) under the Marie Skłodowska-Curie grant agreement No. 955620 (SWATNet). N.W. acknowledges funding from the Research Foundation – Flanders (FWO – Vlaanderen, fellowship no. 1184319N) and the KU Leuven project 3E241013. Computational resources and services used in this work were provided by the Finnish IT Center for Science (CSC) and the FGCI project (Finland) and by the VSC (Flemish Supercomputer Centre), funded by the FWO and the Flemish Government-Department EWI. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

Appendix A

Scaling of SOLPACS equations

In equation (12), the particle distribution function f scales as ϵ_injn, where ϵ_inj is the model parameter describing the efficiency of the shock to inject low-energy particles to the acceleration process (for details, see, e.g., Afanasiev et al., 2023). Taking this into account and using a scaled wavenumber k* = k/B in equation (15), we infer that $\mathrm{\Gamma }\left(k\right)=\mathrm{\Gamma }\left(k^{\mathrm{*}}\right)\propto {ϵ}_{\mathrm{inj}}\sqrt{n}$. Also, we note that Γ does not depend explicitly on B.

Furthermore, one can write for a power-law wave spectrum (or for a portion of the wave spectrum approximated by a power law):

$$I\left(k\right)=I_0{\left(\frac{k}{k_0}\right)}^{-q_0}=I_0{\left(\frac{k^{\mathrm{*}}}{k_0^{\mathrm{*}}}\right)}^{-q_0}\equiv I^{\mathrm{*}}\left(k^{\mathrm{*}}\right).$$(A1)

Therefore, using the scaled wavenumber k*, the pitch-angle diffusion coefficient can be given as

$$D_{\mu \mu }=\frac{\pi }{2}{\mathrm{\Omega }}_0^{\mathrm{*}}\frac{\left|k_{\mathrm{r}}^{\mathrm{*}}\right|I^{\mathrm{*}}\left(k_{\mathrm{r}}^{\mathrm{*}}\right)}{\gamma }\left(1-{\mu }^2\right),$$(A2)

where ${\mathrm{\Omega }}_0^{\mathrm{*}}={\mathrm{\Omega }}_0/B=e/m_{\mathrm{p}}$, so it does not depend on the magnetic field magnitude. Thus, under any transformation of the plasma and shock parameters such that the parameters M_A, V_A and ${ϵ}_{\mathrm{inj}}\sqrt{n}$ are constant, equations (12) and (13) are invariant. We also have to require the constancy of θ_Bn in order to have invariant boundary conditions for particles at the shock. This transformation can also be summarised as follows:

$$\begin{array}{c}{M}_{\mathrm{A}}=\mathrm{const}.\\ V_{\mathrm{A}}=\mathrm{const}.\\ {\theta }_{\mathrm{Bn}}=\mathrm{const}.\\ n\to \alpha n\\ {ϵ}_{\mathrm{inj}}\to {\alpha }^{-1/2}{ϵ}_{\mathrm{inj}},\end{array}$$

where α is a constant. Note that while the particle equation is invariant, the distribution function of accelerated particle being proportional to ϵ_injn scales under this transformation as α^1/2.

We have performed test simulations to confirm the discovered scaling property. Figure A1 shows an example of testing this transformation property for B = 0.335 G, n = 3.55 · 10⁶ cm⁻³, V_sh = 1500 km/s, θ_Bn = 0, and α = 0.01. It can be seen that, as expected, the wave spectrum as a function k/B is invariant, and the particle intensity scales as α^1/2 = 0.1 while the spectral shape is preserved.

Figure A1 Comparison of the Alfvén wave spectra (the blue and red lines depict the opposite circular polarizations) in the vicinity of the shock (top row) and the proton spectra at the shock (bottom row) resulting from two simulation runs to test the scaling property of the SOLPACS equations.

Appendix B

Bell’s theory

The 1D steady-state theory of particle acceleration in a parallel shock accounting for self-generated Alfvén waves (Bell, 1978) can be formulated as follows (see also Vainio et al., 2014; Afanasiev et al., 2015):

$$f(x,p)=f_{\mathrm{s}}\frac{x_0}{x+x_0},$$(B1)

$$\lambda (x,v)=\frac{C_{\lambda }}{v}(x+x_0),$$(B2)

$$f_{\mathrm{s}}\left(p\right)=C_f{\left(\frac{p}{p_{\mathrm{inj}}}\right)}^{-\sigma },$$(B3)

$$x_0\left(p\right)=C_x{\left(\frac{p}{p_{\mathrm{inj}}}\right)}^{\sigma -3},$$(B4)

where f is the particle distribution function, λ is the scattering mean free path, v and p are the particle speed and momentum, σ = 3r_c/(r_c − 1) is the spectral index of the particle distribution function at shock, $r_{\mathrm{c}}=r_{\mathrm{g}}\left(1-M_{\mathrm{A}}^{-1}\right)$ is the scattering-centre compression ratio, r_g is the gas compression ratio of the shock, p_inj is the effective injection momentum (see Afanasiev et al., 2015, for details), and

$$C_{\lambda }=3\left(u_1-V_{\mathrm{A}}\right), C_f=\frac{\sigma {ϵ}_{\mathrm{inj}}n}{4\pi p_{\mathrm{inj}}^3}, C_x=\frac{2}{\pi {ϵ}_{\mathrm{inj}}\sigma }\frac{V_{\mathrm{A}}}{{\mathrm{\Omega }}_0},$$

are the functions of the plasma and shock parameters, which are constant in the context of local SOLPACS simulations.

References

All Tables

All Figures

Figure 1 Diagram of the PARASOL model. The thin black line upstream of the shock schematically depicts the interface between the inner foreshock modeled analytically based on self-consistent SOLPACS simulations and the outer foreshock simulated in PARADISE. It should be noted that the actual shock-normal distance from the shock front to the model matching interface depends on the shock magnetic geometry (see Sect. 2 for details).

In the text

	Figure 2 Flow chart of PARASOL indicating the input required by the constituting models.
In the text

	Figure 3 Examples of the fits of the spatial (left panel) and energy (right panel) dependencies of the mean free path in the simulation run 6. The former is obtained at E = 5.27 MeV. The crosses indicate points that were not considered during the fitting process.
In the text

	Figure 4 Left panel: Cutoff energy obtained in Runs 1–7 versus the cutoff energy calculated based on equation (27). The red line is a power-law fit given by equation (26). Right panel: Example of the simulated particle spectrum at the shock (blue line) and the reproduced spectrum based on the analytical model (orange line). The case presented is run 14.
In the text

	Figure 5 Example of temporal evolution of the parameters of the fitting function Λ(E) in a SOLPACS simulation. The plots show each parameter versus a simulation snapshot number. The case presented is run 13.
In the text

	Figure 6 Fit to the amplitude factor K of the mean free path Λ at the shock versus VA/(ϵinjΩ0).
In the text

	Figure 7 Fits to the parameters of the energy dependence of the foreshock spatial scale Δx.
In the text

	Figure 8 Fits to the parameters of the energy dependence of the foreshock spatial parameter q.
In the text

Figure 9 Comparison of observed and simulated proton intensities and solar wind plasma parameters for the 12 July 2012 SEP event. The panels correspond to (a) omnidirectional particle intensities, (b) IMF magnitude, (c) solar wind speed and (d) number density. The left figure displays in-situ observational data obtained from the ACE and WIND spacecraft (refer to the text for further details). In contrast, the central and right figures showcase simulation results for Earth’s position and the W25 observer. The vertical lines indicate the arrival of the observed and simulated shock waves at the corresponding observers.

In the text

Figure 10 Energy spectra of observed and simulated proton intensities for the 12 July 2012 ESP event. Panel (a) presents the energy spectrum at the shock arrival, while panel (b) displays it eight hours before the shock arrival. Solid lines represent simulated intensities, with the orange line indicating Earth and the green line representing the W25 observer. Dashed lines correspond to in-situ data from ACE and SEPEM, as also depicted in panel 1a of Figure 9.

In the text

	Figure 11 Modeled omnidirectional particle intensities for the W25 observer in a PARASOL simulation with ϵinj = 10−4.
In the text

Figure 12 The modeled omnidirectional proton intensities within the plane of constant latitude encompassing Earth, 27 h after the insertion of the CME in the simulation domain. The left panel displays proton intensities within the 310–580 keV range, while the right panel presents data for the 15–31 MeV range. The white line denotes the position of the CME-driven shock wave. Earth and the W25 observer are marked by green and blue dots, respectively, and their magnetic connectivity to the shock wave is also illustrated.

In the text

Figure 13 The CME-driven shock, observed 27 h after its introduction into the simulation domain, is depicted. The shock surface is color-coded to represent various parameters: the shock obliquity θBn (a); the Mach Alfvén number MA (b); the intensity J and defined in equation (23) (c); and the cutoff energy Ec, as defined in equation (26) (d). Green and blue dots indicate Earth and the W25 observer, respectively, and their magnetic connections to the shock wave are also illustrated. The white lines on the shock indicate the location of the solar equatorial plane and the meridional slice containing Earth.

In the text

	Figure A1 Comparison of the Alfvén wave spectra (the blue and red lines depict the opposite circular polarizations) in the vicinity of the shock (top row) and the proton spectra at the shock (bottom row) resulting from two simulation runs to test the scaling property of the SOLPACS equations.
In the text