© E. Husidic et al., Published by EDP Sciences 2024

1 Introduction

Space weather events are having an ever-growing socio-economic impact on today’s societies that are heavily reliant on nanotechnology, Global Navigation Satellite Systems (GNSS), and satellites for communication (Hapgood, 2011; Schrijver et al., 2015). Solar eruptive events, especially coronal mass ejections (CMEs), stand out as the primary drivers of hazardous conditions in the inner heliosphere that spacecraft and astronauts are exposed to (Gopalswamy, 2018). Moreover, fast CMEs directed towards Earth can interact with the terrestrial magnetic field, causing geomagnetic storms and potentially affecting electrical power grids, pipelines, and telecommunication (Hapgood, 2011; Schrijver et al., 2014, 2015).

A frequent consequence of solar eruptions is the acceleration of charged particles that are collectively referred to as solar energetic particles (SEPs), comprising electrons, protons, and ions (Reames, 2017). Historically, two types of SEP events were distinguished, namely impulsive and gradual SEP events, where the first type was linked to solar flares and the second type to CMEs (Reames, 1999). However, Kallenrode (2003), and later also (Cane et al., 2010), showed that the largest SEP events include contributions from both flare and CME shock acceleration. Besides solar eruptive phenomena, the shocks in corotating interaction regions (CIRs) function occasionally as particle accelerators to generate energetic particles in the inner heliosphere, and also accelerate more efficiently particles in the outer heliosphere (Richardson, 2004; Wijsen et al., 2019). A CIR is a region where a fast solar wind stream catches up with a slow solar wind stream and compresses it. Some authors (Jian et al., 2006; Richardson, 2018) reserve the term CIR for interaction regions with a lifetime greater than one solar rotation and refer to events of shorter duration as stream interaction regions (SIRs). In the present paper we refer to all such interaction regions as CIRs. At radial distances well beyond 1 au, CIRs are often bounded by two shock waves, a forward shock wave accelerating the slow wind, and a reverse shock decelerating the fast wind. However, in the majority of cases observed at 1 au, CIRs are only bounded by two compression waves, and the shock waves form further out (Richardson, 2018). In observational data, the suprathermal particle intensity profiles often show a double-peaked profile as a result of particles being accelerated at both the forward and reverse shock/compression waves (Smith & Wolfe, 1976; Richardson, 2018; Allen et al., 2021; Wijsen et al., 2021).

In large SEP events, accelerated particles can achieve energies of several hundreds of MeV/nuc or even up to a few GeV/nuc, posing a dangerous threat to satellites and astronauts. Therefore, there is significant interest in understanding the acceleration and transport processes behind SEP events. Despite huge and decades-long efforts, there is still no consensus in the scientific community about the exact acceleration mechanisms for SEPs (Petrosian, 2016; Perri et al., 2022). A frequently proposed candidate for the acceleration mechanism in impulsive events at solar flare sites during magnetic reconnection is stochastic acceleration (SA) due to turbulence (Petrosian & Liu, 2004). Possible explanations for the acceleration of SEPs at shock fronts include diffusive shock acceleration (DSA) and shock drift acceleration (SDA), accompanied by SA (Desai & Giacalone, 2016; Vainio & Afanasiev, 2018). In the SDA mechanism that occurs when the charged particle interacts with a quasi-perpendicular shock front, the particle travels multiple times upstream and downstream of the shock. The difference in magnetic field strengths in the two regions leads to a change in the cyclotron radius of the particle which shifts the particle’s guiding centre along the electric field, resulting in a continuous acceleration. The DSA mechanism dominates at quasi-parallel shocks. Here, the charged particle also crosses multiple times the shock front to gain energy but is scattered at frozen-in magnetic turbulence on both sides of the shock, altering the particle’s pitch angle. The SA occurs due to resonant wave-particle interactions between the SEP and magnetic field fluctuations (for more details on the acceleration mechanisms, see e.g. Kallenrode, 2004; Vainio & Afanasiev, 2018).

In recent years, several simulation models have been developed both for theoretical studies and the forecasting of SEP and CIR events, characterised by a common feature, the coupling of a heliospheric solar wind model with a particle transport code. A few of these models shall be mentioned in the following. The SEPMOD (Solar Energetic Particle MODel) developed by Luhmann et al. (2017) uses background solar wind configurations generated by the heliospheric three-dimensional (3D) magnetohydrodynamic (MHD) model Enlil (Odstrcil et al., 2004, 2005) to simulate the scatter-free transport of SEPs. Another SEP propagation model is the Energetic Particle Radiation Environment Module (EPREM), which includes the diffusion and drift effects of the particles, introduced by Kozarev et al. (2010). EPREM was initially coupled to the Enlil model (Kozarev et al., 2010), while later Kozarev et al. (2013) advanced EPREM to take input from the Block Adaptive Tree Solar-Wind Roe Upwind Scheme (BATS-R-US) model (Manchester IV et al., 2012). Recently, the PArticle Radiation Asset Directed at Interplanetary Space Exploration (PARADISE) code (Wijsen, 2020) has been introduced that uses solar wind background configurations by the heliospheric 3D MHD model European Heliospheric Forecasting Information Asset (EUHFORIA; Pomoell & Poedts, 2018) to propagate SEPs as test particles and calculate the intensity distributions of the energetic particles. Similar to Enlil, EUHFORIA incorporates a uniform equidistant grid. EUHFORIA + PARADISE has been applied successfully both in theoretical studies with synthetic solar wind maps (Wijsen et al., 2019) and by reproducing real SEP events (Wijsen et al., 2021, 2023).

In the present work, we introduce our new model, Icarus + PARADISE. Similar to EUHFORIA, Icarus is a 3D MHD heliospheric wind model that simulates solar wind configurations from 0.1 au onward and models the propagation and evolution of transient structures such as CMEs and CIRs (Verbeke et al., 2022). When Icarus uses an equidistant grid, it functions similarly to EUHFORIA or Enlil. However, unlike the latter models, Icarus can incorporate advanced numerical techniques such as radial grid stretching and (solution) adaptive mesh refinement (AMR), to increase or lower the spatial resolution in selected regions of the computational domain. For demonstration purposes, we choose in this work a synthetic idealised solar wind map to provide the inner boundary conditions at 0.1 au. In the resulting heliospheric solar wind model that contains CIRs, we use our particle transport code PARADISE to inject 1 MeV protons in the forward and reverse shock of the CIR. By solving the focused transport equation (FTE; see e.g. van den Berg et al., 2020 and references therein), intensities of the energetic particles are obtained. To validate our new model, we compare simulation results to the EUHFORIA + PARADISE model that has been successfully tested with real data. Finally, we apply different levels of AMR to the background solar wind configurations of Icarus and show the significance of an increased resolution at the shock waves for the obtained particle intensities.

The paper is structured as follows. In Section 2 we describe the energetic particle code PARADISE (Sect. 2.1), the heliospheric models Icarus and EUHFORIA (Sect. 2.2), the setup for the subsequent simulations (Sect. 2.3) as well as considerations made for the coupling of PARADISE to Icarus (Sect. 2.4). Results of the simulations are displayed in Section 3, where the Icarus + PARADISE model is validated by reproducing results of the established EUHFORIA + PARADISE model (Sect. 3.1), and subsequently applied with different levels of AMR (Sect. 3.2). The paper concludes with a summary and discussion in Section 4.

2 Numerical models and setup

The Icarus + PARADISE model represents a coupling of two models. The PARADISE code is used to calculate the spatio-temporal evolution of energetic particles, and its functionality and setup are described in Section 2.1. PARADISE requires as input background solar wind configurations which are provided by the Icarus code. Together with the previously coupled EUHFORIA model that is used here to validate the Icarus + PARADISE, the two heliospheric models are illustrated in Section 2.2. The setup for the simulations in Section 3 is described in Section 2.3, while Section 2.4 outlines considerations and challenges regarding the coupling of PARADISE to Icarus.

2.1 The PARADISE code

PARADISE is a particle transport code utilised for solving the focused transport equation (FTE) to derive spatio-temporal distributions of energetic particles (Wijsen, 2020). The FTE (Roelof, 1969; Skilling, 1971, 1975; Ruffolo, 1995; Isenberg, 1997; le Roux & Webb, 2009) can be derived from the Vlasov equation under the assumption that the energetic particles propagate as test particles through a prescribed solar wind background. It describes the evolution of the five-dimensional gyrotropic particle distribution function f(x, p, μ, t), which depends on various parameters: spatial coordinates x, momentum magnitude p, pitch angle cosine μ ≡ cos(α) with pitch angle α = arctan(v_⊥/v_∥) where, respectively, v_⊥ and v_∥ denote the velocity components perpendicular and parallel with respect to the ambient mean magnetic field, and time t. PARADISE implements the FTE in the form

$$\frac{\partial f}{\partial t}+\frac{\mathrm{d}\mathit{x}}{\mathrm{d}t}\cdot \nabla f+\frac{\mathrm{d}\mu }{\mathrm{d}t}\frac{\partial f}{\partial \mu }+\frac{\mathrm{d}p}{\mathrm{d}t}\frac{\partial f}{\partial p}=\frac{\partial }{\partial \mu }\left(D_{\mu \mu }\frac{\partial f}{\partial \mu }\right)+\nabla \cdot \left({\mathit{\kappa }}_{\perp }\cdot \nabla f\right)$$(1)

with D_μμ being the pitch angle diffusion coefficient, and κ_⊥ the spatial cross-field diffusion tensor (see below for details). The total derivatives of x, μ, and p in equation (1) are given by

$$\frac{\mathrm{d}\mathit{x}}{\mathrm{d}t}={\mathit{V}}_{\mathrm{sw}}+{\mathit{V}}_{\mathrm{d}}+\mu v \mathit{b},$$(2)

$$\frac{\mathrm{d}\mu }{\mathrm{d}t}=\frac{1-{\mu }^2}{2}\left(v\nabla \cdot \mathit{b}+\mu \nabla \cdot {\mathit{V}}_{\mathrm{sw}}-3\mu \mathit{bb}:\nabla {\mathit{V}}_{\mathrm{sw}}-\frac{2}{v}\mathit{b}\cdot \frac{\mathrm{d}{\mathit{V}}_{\mathrm{sw}}}{\mathrm{d}t}\right),$$(3)

$$\frac{\mathrm{d}p}{\mathrm{d}t}=\left[\frac{1-3{\mu }^2}{2}\left(\mathit{bb}:\nabla {\mathit{V}}_{\mathrm{sw}}\right)-\frac{1-{\mu }^2}{2}\nabla \cdot {\mathit{V}}_{\mathrm{sw}}-\frac{\mu }{v}\mathit{b}\cdot \frac{\mathrm{d}{\mathit{V}}_{\mathrm{sw}}}{\mathrm{d}t}\right]p .$$(4)

Here, V_sw denotes the solar wind bulk velocity, V_d the drift velocity of a particle’s guiding centre (GC) accounting for magnetic drifts (in the present work we are not concerned with GC drifts and thus set V_d = 0), v ≡ p/(γm) is the particle speed with γ being the Lorentz factor and m the particle rest mass, and b the unit vector in direction of the mean magnetic field. Furthermore, the colon in equations (3) and (4) denote the Frobenius inner product, which can be expressed by using Einstein’s summation convention as bb : ∇V_sw = b_ij ∂v_i/∂x_j. The FTE in equations (1)–(4) is written in mixed coordinates, that is, the spatial coordinates are measured in an observer’s inertial frame, while the momentum magnitude and the pitch angle are given in a frame co-moving with the solar wind. We note that equations (1)–(4) ignore terms of the order of V_sw/c, where c denotes the speed of light. This is a reasonable approximation, given that V_sw/c ∼ 10⁻³. For a derivation of the FTE retaining terms of higher order in V_sw/c, see, for instance, Chapter 2 of Wijsen (2020), and references therein. Inside PARADISE, equations (1)–(4) are rewritten to an equivalent set of stochastic differential equations (SDEs) that are integrated forward in time using Itô calculus (Wijsen, 2020; Strauss & Effenberger, 2017).

Global heliospheric MHD models, such as EUHFORIA and Icarus (see Sect. 2.2), derive large-scale solar wind configurations, where the interplanetary magnetic field (IMF) follows a nominal spiral pattern unless it is perturbed by transient structures like CMEs and CIRs. However, even in the absence of such structures, on small spatial and temporal scales, the IMF exhibits turbulent fluctuations that are not captured by the MHD model and which can, for instance, scatter the energetic particles. In the FTE, the effect of such turbulent fluctuations on the energetic particle transport is modelled through a set of diffusion processes in phase space. An expression of these diffusion coefficients can be derived from, for instance, quasi-linear theory (Jokipii, 1966).

In PARADISE, the diffusion processes are treated in a modular fashion, that is, the diffusion terms in the FTE describing turbulence can be added, removed, or exchanged. PARADISE contains different expressions for the pitch angle diffusion coefficient which are derived under the assumption of quasi-linear theory (QLT) and magnetic slab turbulence (see e.g. Jokipii, 1966 or Jackel & Schlickeiser, 1992). In the present work, we use an expression similar to those of Agueda et al. (2008) and Agueda & Vainio (2013), given by

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

Here, ϵ describes a parameter to resolve the resonance gap at μ = 0 (see e.g. Klimas & Sandri, 1971; in the present work, ϵ = 0.048), and the scaling factor D₀ is determined by the particle’s rigidity and assuming a parallel mean free path (MFP) that is constant in space. In particular, we take the parallel MFP as λ_∥ = 0.3 au for a 1 MeV proton. For details, we refer to Wijsen (2020).

Besides pitch angle diffusion, we consider also cross-field spatial diffusion (CFD) described by the diffusion tensor κ_⊥. Spatial CFD is a proposed mechanism to explain, for instance, SEP events that are observed by spacecraft which are seemingly not magnetically connected to the particle source (e.g. Klassen et al., 2015). In this work, we account for CFD in the PARADISE simulations by prescribing a constant perpendicular MFP length of λ_⊥ = 3 × 10⁻⁴ au. Thus, with λ_⊥/λ_∥ = 10⁻³, the particles will predominantly propagate along the IMF lines in our simulations.

2.2 Heliospheric solar wind models

Icarus, as detailed in Verbeke et al. (2022) and Baratashvili et al. (2022), operates within a reference frame corotating with the Sun. Icarus solves the ideal 3D MHD equations to obtain heliospheric solar wind configurations from 0.1 au onward. A simple Icarus solar wind and a cone CME model is maintained as a test case in the Message Passing Interface-Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC), which constitutes a collection of various numerical schemes tailored for handling (near-)conservation laws and shock problems in hydrodynamics and magnetohydrodynamics (Keppens et al., 2012; Xia et al., 2018). One notable feature of Icarus is its flexibility regarding the choice of a coronal model that provides the inner boundary conditions, as long as the solar wind is considered super-Alfvénic at 0.1 au (in fact, the speed has to be above the fast magnetosonic speed). While this is a common approach in MHD solar wind modelling, recent observations by the Parker Solar Probe indicate that the solar wind can still be sub-Alfvénic at distances beyond 0.1 au (Bandyopadhyay et al., 2022; Jiao et al., 2024). Moreover, Icarus is able to exploit the output of EUHFORIA’s coronal model (see next paragraph), providing it with the inner boundary conditions. The time advantage of radial grid stretching and AMR applications in Icarus compared to simulations of EUHFORIA employing an equidistant grid has been demonstrated in studies conducted by Verbeke et al. (2022) and Baratashvili et al. (2022).

EUHFORIA is a physics-based space weather forecasting tool that generates solar wind configurations similar to Icarus from 0.1 au onward, including transient structures such as CMEs (Pomoell & Poedts, 2018). The model comprises two primary components, a coronal and a heliospheric part. The semi-empirical coronal model accepts as input synoptic magnetograms by the Global Oscillation Network Group (GONG) and employs the Wang-Sheeley-Arge model in order to obtain the solar wind parameters at the inner boundary. Using these solar wind parameters as inner boundary conditions, the heliospheric part of EUHFORIA solves the ideal 3D MHD equations to generate solar wind configurations from 0.1 au onwards. Unlike Icarus, EUHFORIA solves the ideal MHD equations using Heliocentric Earth EQuatiorial (HEEQ) coordinates with Earth’s longitude fixed at 0° on an equidistant grid without the ability of AMR.

2.3 Simulation setup

For the simulations outlined in Section 3 we use the following setup. The background solar wind of Icarus is displayed in the left panel of Figure 1 and was obtained by using a synthetic solar wind map at the inlet boundary, 0.1 au, for the heliospheric wind model. The map is displayed in Figure 2 and shows the radial solar wind velocity component of a tilted slow solar wind stream embedded in a fast wind stream which is based on the work of Pizzo (1991). The band of slow solar wind represents a wind coming from a helmet streamer around the magnetic equator of the Sun. Such a configuration produces two CIRs as shown in Figure 1. Following Pizzo (1991), in all simulations it is assumed that the magnetic axis of the Sun is tilted by 30° with respect to its rotation axis.

Figure 1 Contour plots showing the radial velocity component in the ecliptic plane from 0.1 au to 2 au by Icarus (left panel), and the corresponding divergence of the solar wind velocity in the Icarus simulation.

Figure 2 Contour plot of the radial speed at the inner boundary both for Icarus and EUHFORIA.

To determine the injection region for the test particles, we compute the divergence of the solar wind velocity of the heliospheric wind. As an example, in the right panel of Figure 1, we present a cross-section at θ = 90° of the solar wind velocity divergence in the Icarus wind. The forward and reverse shocks of the CIR are clearly identifiable where the velocity divergence is negative. Given the perfect symmetry of this structure, we arbitrarily select the CIR that extends to the lower half-plane of the right panel of Figure 1 for injecting 1 MeV protons within a region spanning from 0.12 to 2 au, and limited in colatitude in the interval 90° ± 2° at the shock wave. We note that injection energy on the order of 1 MeV is significantly higher than the anticipated energies of seed particle populations at CIR shocks. This divergence arises from the primary objective of our present study, which is to validate our model and illustrate the impact of AMR on acceleration efficiency within the model, rather than replicating a real CIR event.

The particles are injected according to the differential intensity formula, represented as

$$j(r,E)=C\delta (E-1 \mathrm{MeV})r^{-2}$$(6)

with C being a proportionality constant and δ the Dirac delta distribution, and where div V_SW < 0 and j = p²f. In practical terms, this is accomplished within PARADISE by uniformly distributing 1 MeV protons throughout the region where div V_SW < 0 (tracing the CIR shock) and assigning them a statistical weight of 1/r². To trace the shock, we used a modified version of the shock tracer code by Wijsen et al. (2022). Furthermore, we assume that j is independent of time. To achieve this, we first calculate the Green’s function solution for the FTE, which emerges when all particles are injected at time t = 0. Subsequently, we derive the steady-state solution for particle injection with constant time dependence by convolving Green’s function solution with a time-independent function. This is done for all simulations in Section 3.

For the application of AMR in Section 3.2, we use a region criterion to determine the refinement location, constraining the injection region to one CIR. Following the approach outlined in Wijsen et al. (2021) and Verbeke et al. (2022), we estimate the longitudinal position $\widehat{\varphi }$ of the CIR that we wish to refine by

$$\widehat{\varphi }=\varphi +\frac{r-r_{\mathrm{B}}}{U}\mathrm{\Omega } .$$(7)

In equation (7), [r, ϕ] denote coordinates in the domain, r_B is the radial distance at the inner boundary (i.e., 0.1 au), U is the characteristic speed of the fast stream (in the present work, U = 500 km s⁻¹), and Ω is the solar synodic rotation rate (Ω ≈ 2.67 ∙ 10⁻⁶ s⁻¹). To approximate the CIR spiral, equation (7) must satisfy the condition

$${\varphi }_{\mathrm{lower}}<\widehat{\varphi }<{\varphi }_{\mathrm{upper}},$$(8)

with ϕ_lower and ϕ_upper denoting the upper and lower longitudinal limits, respectively. We note that the standard setting of Icarus uses the solar sidereal rotation rate, but since EUHFORIA uses the solar synodic rotation rate, we adjusted it for the present study.

The MPI-AMRVAC framework contains a number of numerical schemes and shock-capturing methods that can be selected. In this particular study, we choose the second-order (in time and space) total variation diminishing Lax-Friedrichs (TVDLF) scheme as our shock-capturing method (Tóth & Odstrcil, 1996) combined with the Woodward limiter (Woodward & Colella, 1984). The TVDLF is a robust scheme, but more diffusive in comparison to other solvers, such as the Harten-Lax-van Leer (HLL) Riemann solver. A key distinction between EUHFORIA and Icarus is that in EUHFORIA ∇∙B is prescribed to be 0. While MPI-AMRVAC has the option to also use a constraint transport approach, in Icarus, small non-zero values of ∇∙B are locally allowed. The parabolic cleaning method by Keppens et al. (2003) is included, which diffuses any local non-zero divergence of the magnetic field.

To validate the Icarus + PARADISE model in Section 3.1, we conducted simulations replicating outcomes of the EUHFORIA + PARADISE model. To achieve this, we generated a similar heliospheric wind in EUHFORIA, employing the same coronal map utilised for Icarus, as presented in Figure 2. In Figure 3, we provide a comparison of the longitudinal profiles of the magnetic field magnitude (top panel), the solar wind speed (middle panel), and the particle number density (bottom panel) at 1.0 au and θ = 90°. The results are in reasonable agreement with slight differences in all three solar wind parameters which are due to the entirely different numerical schemes used in EUHFORIA and Icarus. The spatial domain of both winds covers a radial distance of 0.1–2 au since particle acceleration in CIRs often occurs beyond 1.0 au, and a colatitudinal and longitudinal range of 30°–150° and 0°–360°, respectively. The resolution of both winds is the same with 600 cells in radial direction or Δr ≈ 1.2 R_S, 60 cells in colatitude and 180 cells in longitude, which corresponds to an angular resolution of Δθ = Δϕ = 2°.

Figure 3 Plots of the magnetic field magnitude (top panel), the solar wind speed (middle panel), and the solar wind particle number density (bottom panel) versus longitude ϕ at radial distance r = 1 au and colatitude θ = 90° obtained by Icarus (blue solid line) and EUHFORIA (red dashed line).

The simulations for the comparison of different levels of AMR in Section 3.2 utilises the same Icarus solar wind as in Section 3.1. However, on AMR level 1, which corresponds to a uniform initial grid, we choose a low resolution of only 300 cells in radial direction, 40 cells in colatitude and 96 cells in longitude. With each subsequent AMR level, the number of cells in each direction doubles, at least locally, i.e., where the given AMR criteria are satisfied. In the highest AMR level achievable in the chosen simulation setup, AMR level 5, there are “effectively” 4800 cells in the radial direction, 640 “effective” cells in colatitude, and 1536 “effective” cells in longitude, reflecting the highest refinement level.

Finally, we note that we made two modifications in the inner boundary conditions in comparison to the standard Icarus version available in MPI-AMRVAC. Since larger discrepancies in v_ϕ and B_ϕ in Icarus and EUHFORIA at and near the inner boundary were noticed resulting from the use of different coordinate systems, we changed the inner boundary in Icarus to match those in EUHFORIA. The longitudinal velocity component at the inner boundary was altered to v_ϕ = −r_BΩ sin(θ) in order to obtain v_ϕ ≈ 0 at the inner boundary in the inertial frame. According to Pomoell & Poedts (2018), we also adjusted the longitudinal magnetic field component to B_ϕ = B_r v_ϕ/v_r, where B_r and v_r are the radial magnetic and velocity components, respectively. This change in B_ϕ was made to ensure that the electric field is zero in the corotating frame, a necessary condition for obtaining a steady-state solution of the solar wind in the corotating frame, as implemented in EUHFORIA.

2.4 Coupling of PARADISE to Icarus

The PARADISE particle transport code employs a grid-free numerical approach for the propagation of energetic test particles. However, the input background solar wind configurations required by PARADISE are provided on a numerical grid, containing parameters of the solar wind plasma. To integrate the set of equations constituting the FTE for a specific particle, it becomes necessary to locate the test particle on the numerical grid of the particular heliospheric solar wind model, read out the solar wind plasma variables of the grid points surrounding the particle, and interpolate them at each computational step. For this purpose, the method of trilinear interpolation that has been applied for the previously coupled EUHFORIA model has been adopted also for the implementation of the Icarus model due to the rectangular arrangement of solar wind variables within Icarus.

While EUHFORIA utilises a uniform and equidistant grid contained in a single large rectangular block, the Icarus grid adopts a more complex structure comprising multiple blocks. In the case of Icarus, each level of refinement divides the target block into halves along each coordinate direction. As higher levels of AMR are applied, these blocks with different AMR levels follow a strict nesting protocol, that is, two blocks touching each other either orthogonally or diagonally, will never exhibit an AMR level difference exceeding one (Nool & Keppens, 2002).

The primary challenge encountered when coupling PARADISE with Icarus arose from the block-based structure of the Icarus grid when multiple levels of AMR are applied. When a particle is close to the boundary of a block, that is, when it is in some form outside the internal grid structure of the current block, communication between multiple neighbouring blocks becomes necessary to provide requisite data for numerical computation processes such as trilinear interpolation or finite volume differencing. This issue is illustrated in Figure 4 which has been adapted from the MPI-AMRVAC documentation where further information about the AMR implementation can be found¹. In cases involving more than two levels of AMR, it becomes evident that a block of some level L must eventually communicate with either a so-called parent block of level L − 1 or a so-called children block of level L + 1. Communication between blocks of the same level occurs both with and without applied AMR.

Figure 4 Illustration of the AMR nesting and communication between blocks. Using AMR likely requires blocks of some level L to communicate both to coarse blocks (level L − 1) and further refined blocks (level L + 1).

The solution of having ghost cells for each block, as employed by the built-in particle tracer in AMRVAC, was not feasible for PARADISE due to resulting significantly larger data files storing these ghost cell values. As an alternative, we chose to individually address the various scenarios based on a particle’s position relative to the numerical grid and acquire all necessary data from neighbouring blocks as needed.

3 Simulation results

In this section, we showcase simulation results obtained from the Icarus + PARADISE model. Implementing the setup described in Section 2.3, we first validate our new model by reproducing the results of the previously established EUHFORIA + PARADISE model using in Icarus a uniform grid (Sect. 3.1). Subsequently, we explore the impact of higher levels of AMR to illustrate its effects (Sect. 3.2).

3.1 Validation

Figures 5 and 6 show longitudinal profiles of the particle intensities obtained by PARADISE at radial distances r = 1 au (Fig. 5) and r = 1.5 au (Fig. 6). In each figure, the top panel shows the intensity profiles of three selected energy channels that are defined in the legend below the figures. Additionally, we provide the longitudinal profiles of solar wind speed (middle panel) and the particle number density (bottom panel). The Icarus-based results are shown in dotted lines, while the results obtained in the EUHFORIA wind are shown in the same figures in the same colours, but in solid lines. Similarly, the solar wind parameter profiles are shown in dotted format for Icarus and in solid for EUHFORIA.

Figure 5 Longitudinal profiles of intensities (top panel) obtained for 1 MeV protons at θ = 90° and r = 1.0 au. The intensities obtained in the Icarus wind are represented by dotted lines, while the EUHFORIA-based trends are depicted by solid lines, each in their respective colours. Similarly, the solar wind speed and solar wind particle number density (middle and bottom panel, respectively) are shown in dotted lines for Icarus and solid lines for EUHFORIA.

Figure 6 Longitudinal profiles of intensities (top panel) obtained for 1 MeV protons at θ = 90° and r = 1.5 au. The intensities obtained in the Icarus wind are represented by dotted lines, while the EUHFORIA-based trends are depicted by solid lines, each in their respective colours. Similarly, the solar wind speed and solar wind particle number density (middle and bottom panel, respectively) are shown in dotted lines for Icarus and solid lines for EUHFORIA.

For both models, we observe the result of adiabatic deceleration (channel below 1 MeV) and shock acceleration (channels above 1 MeV). Furthermore, the intensity profiles reveal a characteristic two-peaked structure typically associated with a CIR. The right intensity peak corresponds to particles accelerated at the forward shock of the CIR, whereas the valley between the two peaks coincides with the stream interface (SI) of the CIR, that is, the transition between the compressed slow and fast wind. The left intensity peak is attributed to particles originating from the reverse shock.

In the two lowest energy channels, the two models exhibit good agreement, resembling in shape and following a similar trend. It can be noted that, depending on the longitude, the Icarus-based intensities exceed the EUHFORIA-based results, and vice versa. These small differences can be attributed to the small differences between the underlying solar wind configurations. This variability is expected, given that the underlying models are not identical, and energetic particles traverse a significant distance through these two solar winds before reaching observers. Cumulatively, small differences in both solar winds can impact the particles’ properties. Moreover, in the highest energy channel, it can be seen that the intensities in EUFHORIA peak are higher than in Icarus. This difference can be attributed to the fact that, at the same resolution, Icarus exhibits a slightly higher level of diffusion compared to EUHFORIA. Consequently, the shock in Icarus is less steep than in the EUHFORIA solar wind, resulting in a reduction of particle acceleration efficiency. However, as demonstrated in the following section, the implementation of AMR enables Icarus to achieve a resolution at the shock significantly surpassing the one depicted in Figures 5 and 6.

3.2 Effect of AMR application

In contrast to similar heliospheric MHD models such as Enlil or EUHFORIA, which employ uniform equidistant grids, the MPI-AMRVAC framework empowers Icarus to utilise AMR to refine regions of interest while maintaining a low resolution in other regions. We illustrate the application of five levels of AMR in the panels of Figure 7 that are similarly structured as Figures 5 and 6 with the longitudinal intensity profiles in the top panels, the solar wind speed profiles in the middle panels, and the solar wind particle number density in the bottom panels in each subfigure. Here, AMR level 1 corresponds to a uniform grid with the resolution as described in Section 2.3. The subfigures in Figure 7 depict the intensities at r = 1.8 au.

Figure 7 Longitudinal profiles of intensities obtained for 1 MeV protons in the Icarus wind at θ = 90° for five levels of AMR (panels a through e). Each plot displays intensities at r = 1.8 au together with the background solar wind speed (middle panels) and solar wind particle number density (bottom panels).

Similarly to Figures 5 and 6, we observe two peaks in the intensities, the right peak corresponding to the forward shock, and the left peak corresponding to the reverse shock. For the higher energy channels, the intensity peak associated with the forward shock becomes only visible by increasing the resolution (i.e., higher AMR level). With increasing levels of AMR, the channels of higher energies, that is the channels measuring the accelerated particles, are increasingly populated. It is crucial to note that shock waves simulated by 3D MHD models are typically much wider than actual shock waves in interplanetary space. Since the ratio of the MFP length to the width of the shock plays a pivotal role, simulations with a more refined shock (here, due to higher levels of AMR) result in more efficient particle acceleration for the same MFP length, see e.g. Wijsen et al. (2022).

To elucidate the impact of higher levels of AMR on particle acceleration, Figure 8 displays intensity-energy profiles at different radial distances at the reverse shock. The intensity curves consist of the peaks of the corresponding energy channel at the particular radial distance, where panel (a) shows the results at 0.5 au, panel (b) at 1 au, panel (c) at 1.5 au, and panel (d) at 1.8 au. All panels consistently indicate an increase in the intensity at higher energies with increasing levels of AMR.

Figure 8 Plots of intensities versus energy at r = 0.5 au (panel a), r = 1.0 au (panel b), r = 1.5 au (panel c), and r = 1.8 au (panel d), using different levels of AMR. The panels show the intensity peaks of the different energy channels at the reverse shock from a simulation up to 2 au.

With increasing resolution, the particle intensities and their energy spectra should converge if the underlying solar wind converges. However, it should be kept in mind that by increasing the resolution of the MHD simulation, we are resolving smaller-scaled solar wind structures, which can subsequently affect the energetic particles. While the results in Figure 8 do not show convergence yet, the differences between consecutive AMR levels can be seen to progressively decrease.

4 Summary and discussion

The significance of energetic particles generated in the heliosphere is growing due to humanity’s dependence on space technology and its renewed interest in space exploration. Consequently, there is a pressing need for the development of numerical tools to reliably and timely predict harsh events. In order to make progress in our space weather forecasting capabilities and our theoretical understanding of the transport and acceleration of solar energetic particles, we have introduced the Icarus + PARADISE model. This two-part model comprises the 3D MHD code Icarus, integrated into the MPI-AMRVAC framework, which generates solar wind background configurations. Accompanying Icarus is the energetic particle transport code PARADISE, which takes the Icarus wind as input to propagate energetic particles as test particles through it and tracks their acceleration and deceleration by solving the FTE in a stochastic manner.

To introduce the new model, we presented a theoretical study using a synthetic solar wind map for the coronal model providing Icarus with inner boundary conditions to generate a solar wind configuration that contains two CIRs. As part of the validation process, we reproduced outcomes from the earlier EUHFORIA + PARADISE model and showcased longitudinal intensity profiles at different radial distances. The simulation results demonstrated overall good agreement, exhibiting a strong similarity in the lower energy channels. However, some disparities surfaced in the higher channels, where the intensities derived from EUHFORIA surpassed those obtained in the Icarus simulation for the same underlying MHD resolution. We attribute these variations in intensity profiles to differences between the two solar wind models, stemming from the distinct numerical schemes employed by EUHFORIA and Icarus. For instance, the TVDLF scheme used in Icarus for shock capturing was found to be more diffusive compared to the one utilised in EUHFORIA. Alternative shock-capturing schemes, such as the Harten-Lax-van Leer (HLL) Riemann solver or Roe-type approximate Riemann solvers, are available in the MPI-AMRVAC framework (Keppens et al., 2023) and are known to be less diffusive than the TVDLF scheme. They will be explored for evaluation in future studies due to their potential advantages.

The MPI-AMRVAC framework employed by Icarus offers the flexibility to include AMR. In our study, we strategically utilised a region criterion to refine only the shock region where we injected the test particles while maintaining a lower resolution for the remaining solar wind. Using five levels of AMR, we displayed longitudinal profiles of the intensities at 1.8 au. Furthermore, we extracted the intensity peaks for the different energy channels and illustrated the corresponding intensity-energy plots at the reverse shock at different radial distances. Our findings consistently demonstrated that as the level of AMR increased, a more pronounced particle acceleration was observed. This result holds significant implications for models that simulate particle acceleration at shock waves generated by MHD models, for instance, Wijsen et al. (2022) using EUHFORIA + PARADISE to study proton acceleration in the inner heliosphere, or Young et al. (2021) using EPREM and the CORona-HELiosphere (CORHEL) model to investigate proton acceleration in the corona.

In our future work, we anticipate further applications of the Icarus + PARADISE model. These potential applications encompass, among others, the examination of solar wind configurations of real events involving CMEs and observed SEP events, mirroring previous studies with EUHFORIA + PARADISE. Furthermore, Icarus is also capable of implementing grid stretching, that is the ability to keep grid cells in a more cubic shape that in an otherwise equidistant grid in spherical coordinates would become more elongated. The performance of Icarus with applied radial grid stretching with a combination of AMR is demonstrated in Baratashvili et al. (2022). Such an approach can be tested for SEP studies with sufficient AMR levels in order to gain significant speed-up in computation time while obtaining high-resolution data where needed. Furthermore, a major update for Icarus + PARADISE is foreseen, for instance, a modification of PARADISE such to load only the necessary blocks of the Icarus snapshots, rather than loading entire data files. This and other enhancements are designed to accelerate simulations inside PARADISE, enable the use of higher levels of AMR, and augment the already significant acceleration of simulation time achieved by Icarus.

Acknowledgments

The authors thank Rony Keppens and Fabio Bacchini from KU Leuven, Belgium, for fruitful discussions regarding the MPI-AMRVAC framework. In addition, the authors thank the anonymous referees for their helpful comments that improved the quality of the manuscript. E.H. is grateful to the Space Weather Awareness Training Network (SWATNet) funded by the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 955620. Furthermore, E.H. acknowledges the travel grant V477923N from the Fonds voor Wetenschappelijk Onderzoek – Vlaanderen (FWO). N.W. acknowledges support from the Research Foundation – Flanders (FWO-Vlaanderen, fellowship no. 1184319N). R.V. acknowledges the support of the Research Council of Finland (grant no. 352847). We also acknowledge support from the European Space Agency project “Heliospheric modelling techniques” (ESA Contract No. 4000133080/20/NL/CRS) as part of the GSTP (General Support Technology Programme), the projects C14/19/089 (C1 project Internal Funds KU Leuven), G.0B58.23N and G.0025.23N (WEAVE) (FWO-Vlaanderen), 4000134474 (SIDC Data Exploitation, ESA Prodex-12), and Belspo project B2/191/P1/SWiM. Computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation – Flanders (FWO) and the Flemish Government – department EWI. The editor thanks Pablo Dmitruk, Ian Richardson and David Ruffolo for their assistance in evaluating this paper.

References

All Figures

	Figure 1 Contour plots showing the radial velocity component in the ecliptic plane from 0.1 au to 2 au by Icarus (left panel), and the corresponding divergence of the solar wind velocity in the Icarus simulation.
In the text

	Figure 2 Contour plot of the radial speed at the inner boundary both for Icarus and EUHFORIA.
In the text

	Figure 3 Plots of the magnetic field magnitude (top panel), the solar wind speed (middle panel), and the solar wind particle number density (bottom panel) versus longitude ϕ at radial distance r = 1 au and colatitude θ = 90° obtained by Icarus (blue solid line) and EUHFORIA (red dashed line).
In the text

	Figure 4 Illustration of the AMR nesting and communication between blocks. Using AMR likely requires blocks of some level L to communicate both to coarse blocks (level L − 1) and further refined blocks (level L + 1).
In the text

Figure 5 Longitudinal profiles of intensities (top panel) obtained for 1 MeV protons at θ = 90° and r = 1.0 au. The intensities obtained in the Icarus wind are represented by dotted lines, while the EUHFORIA-based trends are depicted by solid lines, each in their respective colours. Similarly, the solar wind speed and solar wind particle number density (middle and bottom panel, respectively) are shown in dotted lines for Icarus and solid lines for EUHFORIA.

In the text

Figure 6 Longitudinal profiles of intensities (top panel) obtained for 1 MeV protons at θ = 90° and r = 1.5 au. The intensities obtained in the Icarus wind are represented by dotted lines, while the EUHFORIA-based trends are depicted by solid lines, each in their respective colours. Similarly, the solar wind speed and solar wind particle number density (middle and bottom panel, respectively) are shown in dotted lines for Icarus and solid lines for EUHFORIA.

In the text

	Figure 7 Longitudinal profiles of intensities obtained for 1 MeV protons in the Icarus wind at θ = 90° for five levels of AMR (panels a through e). Each plot displays intensities at r = 1.8 au together with the background solar wind speed (middle panels) and solar wind particle number density (bottom panels).
In the text

	Figure 8 Plots of intensities versus energy at r = 0.5 au (panel a), r = 1.0 au (panel b), r = 1.5 au (panel c), and r = 1.8 au (panel d), using different levels of AMR. The panels show the intensity peaks of the different energy channels at the reverse shock from a simulation up to 2 au.
In the text