© The Authors 2024

1 Introduction

Coronal mass ejections (CMEs) are large magnetized plasma clouds erupting from the Sun. They have masses of 10¹³ g up to 10¹⁶ g and erupt with velocities ranging from ~100 km s⁻¹ to ~3000 km s⁻¹ (based on SOHO/LASCO measurements), carrying enormous momentum and often causing strong shocks during their evolution in the inner heliosphere. When CMEs are observed by in situ measurements, they are usually referred to as interplanetary CMEs (ICMES). An ICME has a complex plasma and magnetic structure, consisting of a leading shock front and a turbulent magnetosheath followed by the magnetic plasma cloud, which often resembles the magnetic field configuration it had during the eruption (Kilpua et al. 2017). The topology of the magnetic cloud, as well as the structure and linkage of its magnetic field, is not yet fully understood and probably varies from case to case. Various studies focus on the radial evolution of the CMEs comparing the magnetic field structure at different phases of the evolution (Al-Haddad et al. 2019; Scolini et al. 2021). Overall, ICMEs are the main drivers of the disturbances studied in space weather research. When directed toward Earth, the impact of a CME is notable. Such impacts can cause strong geomagnetic storms, induce electric currents in power grids and even disrupt the transmission. During the “great Quebec blackout,” the Canadian city witnessed the power of a strong CME on March 13, 1989 and they were trying to deal with the consequences for 9 h. The damage caused by regular events, especially during phases of maximum activity of the Sun, has been estimated to accumulate to an economic loss of 10 billion € per year (National Research Council 2008). The potential loss from the future strong Earth-directed events is increasing continuously, since our dependence on telecommunications, navigation, and electronic systems is increasing day-to-day. In 2019, the National Threat and Hazard Identification and Risk Assessment (THIRA) of the US Federal Emergency Management Agency (FEMA) identified space weather as one of the two threats that could potentially disturb our society globally, the other one being a pandemic (THIRA 2019).

During solar minima, CMEs occur only a few times a week, but during solar maxima, the eruptions are far more frequent and several can occur in a single day (Park et al. 2012). To minimize the damage from the arriving ICMEs, realistic, and fast space weather forecasting tools are necessary. Currently, a number of such forecasting tools exist. Examples of physics based forecasting tools are ENLIL (Odstrčil et al. 2004), EUHFORIA (Pomoell & Poedts 2018), SUSANOO-CME (Shiota & Kataoka 2016), and AWSoM (van der Holst et al. 2014). These tools model the solar wind with superposed CMEs on it, propagating them to the Earth and beyond. Currently, ENLIL and EUHFORIA are used in operational settings, performing daily simulation runs to model the solar wind configuration and the CME evolution. The most simple model for a CME is the so-called cone CME model, which represents a homogeneous hydrodynamic plasma cloud that is injected self-similarly and subsequently evolves in the heliosphere. The advantage of this model is that it is simple and quite efficient. Cone CMEs can model the evolution of CME shocks and predict their arrival time at Earth or other locations in the inner heliosphere. However, this model does not take into consideration the internal magnetic field of the CME. As a consequence, the geo-effectiveness of its impact, coming from its interaction with the magnetic field of the Earth, cannot be estimated with this model. In order to predict whether a CME will cause a severe geomagnetic storm upon its arrival at Earth or not, the sign of B_z component of the CME magnetic field, which is the component perpendicular to the equatorial plane, must be modeled accurately. When B_z is negative upon arrival, it is anti-parallel to the magnetic field of the Earth and such a configuration causes substantial magnetic reconnections and severe geomagnetic storms (Kilpua et al. 2017). Therefore, in order to model the B_z component of the magnetic field of the CME, it is important to incorporate magnetized CME models. The simplest magnetized CME model consists of a spheromak magnetic field configuration. Its shape is similar to the cone CME model, but it has an internal magnetic field. In EUH-FORIA, a linear force-free spheromak CME model has been introduced already and it has been used extensively to study different CME events (see e.g., Verbeke et al. 2019, 2022b; Scolini et al. 2019, 2020). In the present paper, we consider the same linear force-free spheromak CME model and implement it in the recently introduced heliospheric model Icarus (Verbeke et al. 2022a).

Icarus is a 3D MHD heliospheric wind and CME evolution model implemented in the framework of MPI-AMRVAC (Keppens et al. 2003; Xia et al. 2018), which is a parallelized architecture, suited for hydro- and magnetohydrodynamic (HD and MHD, respectively) applications. In previous studies, the cone CME was implemented and used for modeling the evolution of CME shocks in the heliosphere. The novelty of Icarus lies in the advanced numerical techniques available in the code that enable better and faster simulations. Baratashvili et al. (2022) examined the effect of radial grid stretching and adaptive mesh refinement (AMR) in the domain on the shock strength of the CMEs and their arrival times, while comparing the wall-clock timings needed for the simulations. Since the interactions of the internal magnetic field of the CME with the magnetic field of the background solar wind play an important role in the CME evolution and its propagation in the heliosphere, the magnetized CME model, namely, the spheromak was now also implemented in Icarus.

The aim of the present paper is to quantify and demonstrate the performance of the new magnetized CME model in Icarus. Different AMR strategies are applied to the CME to refine the grid at the varying location of the CME during its propagation, in order to obtain high resolution simulations and coarsen the grid again after the CME has passed, in order to increase the efficiency of the simulation. For validating the spheromak CME model, the CME that occurred on July 12, 2012 was selected as a typical event to model. This specific event was studied before with the cone CME model in Icarus (Baratashvili et al. 2022) and with both the cone and spheromak CME models in EUHFORIA (Scolini et al. 2019).

The remainder of the present paper is organized as follows. Section 2 describes the numerical model of Icarus. The implementation of the new magnetized CME model is described in detail in Sect. 3. Then we introduce the simulation set-up and the describe the event in Sect. 4. Sections 5 and 6 are devoted to the presentation and analysis of the simulation results and the discussion and conclusions, respectively.

2 Numerical model

Icarus is a 3D MHD heliospheric model recently developed at the Center for mathematical Plasma Astrophysics (CmPA), KU Leuven (Verbeke et al. 2022a). Icarus is implemented in the framework of MPI-AMRVAC (Xia et al. 2018), which has a heavily parallelized architecture, solving different sets of partial differential equations with various numerical methods and lim-iters. Different numerical scheme and flux limiter combinations were examined in the Icarus setting to study the effect on the simulation results and the efficiency (Baratashvili et al. 2023). The experiments were performed for low,- medium-, and high-resolution equidistant grids (see Table 1), in combination with the HLL, HLLC, and TVDLF numerical methods and “van-leer” (van Leer 1974), “koren” (Koren & van der Maarel 1993), “woodward” (van Leer 1977) and “minmod” (Yee & Shinn 1989) flux limiters. Baratashvili et al. (2023) summarized all 36 simulations and compared the timing and the shock sharpness for each simulation. As a result, the combination of the TVDLF numerical scheme and the “woodward” flux limiter was selected as the most efficient. These settings are fixed as default in all the Icarus simulations presented in this paper.

The Icarus heliospheric domain is similar to that in EUHFO-RIA and covers the radial distances from 0.1 AU to 2 AU, 360° in the longitudinal direction and 120° in the latitudinal direction, [−60°, 60°] in co-latitudes, excluding the poles. The output of the EUHFORIA coronal model is used as the input data for the plasma variables at the inner heliospheric boundary for both EUHFORIA and Icarus. In the Icarus simulations, these values are interpolated on the Icarus grid and ideal MHD equations are solved using the obtained values as boundary conditions. In Icarus, unlike EUHFORIA, the reference frame is co-rotating with the Sun, which leads to a true steady solar wind at the end of the relaxation phase (by definition). The time for the relaxation phase is chosen such that the slow wind stream traverses the whole domain from the inner radial boundary to the outward. The parabolic cleaning method is applied to minimize ∇ • B (Dedner et al. 2002).

After the steady background wind is obtained, CMEs are superposed on it in the simulations. So far, in Icarus the coronal mass ejections are introduced with a simple HD cone model, which is a non-magnetized plasma cloud with homogeneous interior. Thus, the internal speed, density, temperature, and pressure are constant. The geometry is extensively described by Scolini et al. (2018). All the mentioned possible shapes are included in Icarus, the same way as in EUHFORIA. The implementation of the cone CME model is described in detail by Verbeke et al. (2022a).

The goal of the new heliospheric model Icarus is to perform simulations accurately, yet effectively. This implies achieving the same results as the original configuration of the heliospheric domain, while spending fewer computational resources and less time. For this purpose, a modified computational grid is used compared to the original EUHFORIA heliospheric simulations. Various advanced techniques are available in the MPI-AMRVAC framework, from which the radial gradual grid stretching and adaptive mesh techniques are addressed in Icarus.

Since the heliospheric domain covers large distances in the radial direction, an equidistant grid is not the most optimal choice. In spherical coordinates, the equidistant grid cells become deformed ever closer to the outer boundaries, and the longitudinal resolution is affected notably. To avoid this issue, gradual radial grid stretching was applied Xia et al. (2018). With this approach, the aspect ratio of the width and the length of the cells is maintained in the whole domain. The choice for the number of cells in the radial direction is the same as given in Baratashvili et al. (2022) and is fixed to N = 60. The number of cells in each resolution for the equidistant and stretched grid simulations are given in Table 1. For the equidistant grid cases, these can directly be translated to the cell sizes, which are constant throughout the domain. For example, in the medium-resolution equidistant grid simulation the size of a cell is ~0.00316(6) AU in the radial direction and 1.875° in the longitudinal and latitudinal directions. For the stretched cases, the longitudinal and latitudinal resolutions are also constant, but in the radial direction the sizes of the cells vary along the radial direction and are not fixed to a constant number.

The second advanced technique exploited in Icarus is called adaptive mesh refinement. The details of block-adaptive AMR in MPI-AMRVAC can be found in Keppens et al. (2003). The AMR can be initiated in the domain by prescribing a criterion for the refinement. Whenever the criterion is met in the simulation, the block is refined to the intended refinement level. Different AMR criteria, suited for heliospheric simulations, are discussed in Baratashvili et al. (2022). The criterion can be aimed at the CME interior, CME and CIR shocks in the solar wind, the combination of the two, or any other regions of interest in the domain. For example, if the CME tracing refinement criterion is prescribed, the CME will be traced along its propagation from the inner heliospheric boundary to the outer one. The whole block containing the portion of the CME is refined, but once the CME passes the block, it is coarsened to the base resolution of the domain, in order to avoid unnecessary high resolution computational grid regions. The base resolution of the simulation is always fixed to the low resolution given in Table 1.

Since the radial resolution of the stretched grid varies, Table 2 summarizes the radial size of the cells at the first Lagrange point L1 with different AMR levels.

The cell size in every AMR level n is 2⁽ⁿ⁻¹⁾ times smaller than in the base resolution of the simulations in every spatial direction, thus, in the next AMR level, every refined grid cell (satisfying the criteria) is cut in eight smaller cells. For example, for a radially stretched grid, the radial cell size at L1 in the AMR level 5 simulation is 0.672 R_⊙, which is 16 times smaller (in every direction) than in the low resolution simulation, corresponding to 10.7 R_⊙ (at L1). Examples for the other AMR levels are given in Table 2.

3 Magnetized CME model

The magnetized CME model implemented in Icarus is a linear force-free (LFF) spheromak model. This is achieved by linking Icarus to the existing spheromak CME model in EUHFORIA (Verbeke et al. 2019). Since the magnetized CME is imported from EUHFORIA and not implemented anew, we discuss the details of the model only briefly. An extensive discussion of the model and its implementation in EUHFORIA is given in Verbeke et al. (2019), while the spheromak CME model is also described and discussed by Kataoka et al. (2009) and Shiota & Kataoka (2016).

The model is in fact a modification of the Gibson-Low magnetized CME model (Gibson & Low 1998), which is not force-free. In the force-free version offered by Shiota & Kataoka (2016), the CME evolves completely through the inner boundary so that its foot points do not stay attached to the inner helio-spheric boundary. Moreover, this LFF Spheromak model implementation pushes it through the inner heliospheric boundary in such a way that positive pressure is guaranteed everywhere, unlike the original Gibson-Low model. The CME is represented by a sphere with a fixed radius when evolving through the inner boundary, thus it is not evolving self-similarly (with growing radius) like the cone model. However, similar to the cone CME model, the plasma variables in the interior are considered to be homogeneous. In order to inject the CME in the heliosphere as a time-dependent boundary condition, each point of the inner heliospheric boundary is checked at every time step in the simulation. When the coordinates of the centre of the CME are given in Cartesian coordinates by (x_CME, z_CME, z_CME), and the coordinates for an arbitrary point at the inner heliospheric boundary are given by (x_IHB , z_IHB , z_IHB), the following condition is checked: $${\left(x_{\text{CME}}-x_{\text{IHB}}\right)}^2+{\left(y_{\text{CME}}-y_{\text{IHB}}\right)}^2+{\left(z_{\text{CME}}-z_{\text{IHB}}\right)}^2\le r_{\text{CME}}^2,$$(1)

where r_CME is the (fixed) radius of the (moving) CME. If this inequality holds for the given point at a certain time, then it belongs to the CME; if not, then it belongs to the background solar wind and should not be adjusted.

The magnetic field is defined in the local spherical coordinate system (r', θ', ϕ'), where the origin is considered the center of the CME, and it exhibits symmetry in the azimuthal direction (ϕ'). The magnetic field vector can be expressed as: $$B=\frac{1}{r^{\prime }\sin {\theta }^{\prime }}\left[\frac{1}{r^{\prime }}\frac{\partial \psi }{\partial {\theta }^{\prime }}{\widehat{r}}^{\prime }-\frac{\partial \psi }{\partial r^{\prime }}{\widehat{\theta }}^{\prime }+Q{\widehat{\varphi }}^{\prime }\right],$$(2)

where ψ and Q are scalar potentials that only depend on r' and θ' (Chandrasekhar 1956), and ${\widehat{r}}^{\prime },{\widehat{\theta }}^{\prime },\text{ and }{\widehat{\varphi }}^{\prime }$denote unit vectors in each of the three spatial coordinate directions. The magnetic field is divergence-free by the given definition, and the solution is determined so that J × B = 0, thus, it is force-free. From the azimuthal force balance, we see that the toroidal field must be the function of a poloidal potential; therefore, we get Q = Q(ψ). In order to obtain the linear force-free model, we assume Q(ψ) = αHψ, where H is the constant describing the helicity of the spheromak, and α is the length scale constant, which is determined from the condition that the magnetic field at the surface is zero. After solving the force-free equation, the poloidal potential is obtained as $$\psi =\frac{B_0}{\alpha }{r}^{\prime }{j}_1\left(\alpha r^{\prime }\right){\sin }^2{\theta }^{\prime },$$(3)

where B₀ is the magnetic field strength and j₁ (x) is the spherical Bessel function of order one. If we expand the magnetic field into the components and use the expression for the poloidal potential, we obtain the magnetic field configuration of the LFF spheromak model as follows: $$B_r^{\prime }=2B_0\frac{j_1\left(\alpha r^{\prime }\right)}{\alpha r^{\prime }}\cos {\theta }^{\prime },$$(4) $$B_{\theta }^{\prime }=-B_0\left[\frac{j_1\left(\alpha r^{\prime }\right)}{\alpha r^{\prime }}+j_1\left(\alpha r^{\prime }\right)\right]\sin {\theta }^{\prime },$$(5) $$B_{\varphi }^{\prime }=H\cdot B_0{j}_1\left(\alpha r^{\prime }\right)\sin {\theta }^{\prime }.$$(6)

The magnetic field is assumed to be zero at the spheromak boundary. If ψ(r = r_CME) = 0, we obtain the condition for the spherical Bessel function of order one $$j_1\left(\alpha r_{\text{CME}}\right)=0$$

using the first zero of spherical Bessel function of order one, we obtain αr_CME ~ 4.4934094579.

The input parameters for the LFF spheromak model also include the parameters characterizing the internal magnetic field configuration, together with the parameters describing the geometry and hydrodynamics characteristics. The overview for the input parameters is given in Table 3, which closely follows the one described by Verbeke et al. (2019).

4 Simulation set-up

4.1 Event description

In order to validate and examine the implemented magnetized LFF spheromak model, we selected a well-studied CME event during which the impact on Earth was observed clearly. The CME originated from the active region, NOAA AR11520. It followed a strong X1.4 class fare, which was observed by LASCO C2 coronagraph at 16:48UT July 12, 2012. The parameters for modeling this CME with the spheromak model, taken from Scolini et al. (2019), are summarized in Table 4.

4.2 Simulation overview

As mentioned before, the magnetized CME model is imported directly from EUHFORIA, it is exactly the same and any differences in the simulations are entirely due to the differences in the numerical methods, including the grid stretching and the mesh adaptivity. Therefore, a link has been set up between Icarus and EUHFORIA to facilitate communication between the models whenever the magnetic field components need to be computed. The simulations are performed with a predefined AMR criterion suited for resolving the complex CME interior (Baratashvili et al. 2022).

4.2.1 Linking Icarus to EUHFORIA

The heliosphere model of EUHFORIA is implemented in Python, while the Spheromak CME model is implemented in C++. The heliosphere model is connected to the spheromak model via a Python-C++ linker to compute magnetic field components at each point belonging to the CME. In order to avoid unnecessary duplication and to ensure exactly the same implementation, Icarus is also connected to the same implemented spheromak model in C++. The visual representation of the connections is plotted in Fig. 1.

Icarus is implemented in a Fortran environment. Hence, a different link has been set up to the spheromak code. First, in EUHFORIA, additional static libraries have been created from the existing spheromak magnetic field computation code, intended for linking to Icarus. In Icarus, on the other hand, a Fortran-C++ linker has been added to the compilation settings and the code has been recompiled. The path to the newly generated static library is indicated for linking to the required library. As a result, the compiled code with these modifications links the spheromak computation code to Icarus. Upon injecting the CME at the inner heliospheric boundary in the simulation, the linked library is called up from Icarus, which first computes and then returns the magnetic field components at the requested points in the computational domain. The parameters characterizing the CME are set locally, as they do not depend on time. Of course, while passing the location of the intended CME point, the difference in the reference frames of the EUHFORIA and Icarus simulations are taken into account.

Fig. 1 Diagram of the procedure to import the spheromak file in Icarus used in EUHFORIA.

4.2.2 AMR criterion

To demonstrate the potential of AMR for resolving the inner magnetic structure of the simulated CMEs, the spheromak CME is injected in Icarus with different resolution grids. The “standard” (for operational use) medium-resolution uniform grid is taken as a reference. To check whether it is possible to improve on the results with advanced techniques, the radial grid stretching can be combined with AMR.

Baratashvili et al. (2022) has examined a range of AMR strategies. The first criterion focuses on the CME interior and is implemented by tracing the density in the heliosphere as the CME propagates toward the outer boundary. The other criterion considered in the paper resolves shocks in the domain by tracing the sign of the divergence of V, hinting at the regions of compression or expansion. The last AMR strategy combines the first two criteria, thus, focuses on the CME interior and its shock-front along its propagation and evolution in the inner heliosphere. The AMR criterion used in all the simulations in this paper is a combination of two criteria, called the “combined” criterion by Baratashvili et al. (2022). This combined criterion is well suited for studies of the evolution of CMEs as it fully resolves the complex CME models. It is therefore also well suited for the considered test case, since the complex internal magnetic field of the CME is modeled in these simulations.

Apart from the CME interior, the CME shock is also maximally resolved to accurately model the shock strength upon arrival at L1. The mesh refinement is limited to a narrow bandwidth around the Earth-Sun line for the space weather forecasting purposes. The modeled time series of solar wind parameters at Earth are compared to each other (for different resolutions) and to observational data. This approach further optimizes the simulations and saves a lot of computational resources. A visual representation of how this combined AMR criterion is applied in the domain is given in Fig. 2. The figure shows cuts of the 3D simulation in the ecliptic plane. A small portion of the ecliptic plane is zoomed in, including the CME, in order to show more details of the CME coverage. The left figure shows the radial velocity values overlaid with the values of the divergence of V (used as the AMR criterion) with corresponding color maps. Thus, in this figure, the shock-front is distinguished with a red color, corresponding to negative div(V) values and it is followed by the CME magnetic cloud, depicted on the basis of V_r values.

The middle figure shows the CME density profile overlaid with the divergence of V values with corresponding color maps. Since the AMR criterion takes into account the CME density and the shock, the combined profile was plotted. The CME density is an independent variable, that only describes the plasma density in the CME interior, and is zero everywhere else in the domain. Hence, the magnetosheath of the CME can be distinguished very easily between the CME density profile and the leading shock.

The right figure shows the same variables as in the middle figure, but with the overlaid mesh. The “combined” refinement criterion has the additional restriction of the bandwidth around the Sun-Earth line, which implies that at every time-step the AMR is applied along the longitude where the Earth lies with a given margin ∆ϕ = ±30°. For example, if Earth is at a longitude ϕ = 60°, the areas meeting the implemented AMR criterion are refined to the indicated refinement level between longitudes 30° < ϕ < 90°. As a result, the whole area including the CME shock and the CME interior in the heliosphere is not refined to the highest level of refinement, but only the regions along the Sun-Earth line.

5 Results and discussion

The event study was performed in order to validate the integration of the magnetized LFF Spheromak CME model in Icarus. As an input for the EUHFORIA coronal model, we selected the standard synoptic magnetogram from the Global Oscillation Network Group (GONG) on July 12, 2012 at 11:54UT. Then, the output of the coronal model, given as the MHD parameters at 0.1 AU, was used as the inner boundary conditions for the onset of the heliospheric simulations in Icarus. They were first extended radially to 2 AU and then Icarus performed an MHD relaxation period of 14 days to obtain a steady background wind. The parameters used for the CME injection from the inner heliospheric boundary are given in Table 4.

Several simulations were performed on different computational grids to compare the performance of the advanced techniques with that of the original uniform grid simulations. The co-latitudinal component of the magnetic field obtained in the simulations is plotted in Fig. 3. The B_clt values are saturated with (−0.3,0.3) nT to emphasize the difference between regions with negative and positive values. For all the simulations with the stretched grid and AMR, the AMR criterion is fixed to the combined refinement criterion described in Sect. 4.2.2. The top-left figure represents the simulation result obtained with AMR level 2, with one refinement level. The top-right figure represents the results obtained with AMR level 3. The bottom-left figure corresponds to the simulation with AMR level 4, and the bottom-right figure shows the results obtained without radial grid stretching, nor AMR. In the latter simulation, the original uniform medium-resolution grid is used. For each simulation, the same snapshot is taken, which is chosen in the time interval in which the spheromak is crossing the sphere with radius of 1 AU.

The interior structure of the spheromak can be well distinguished from the background solar wind at 1 AU. The background wind is not very well resolved in the simulations with the stretched grid and AMR. As a result, more details and smaller-scale structures can be distinguished in the background wind simulated in the uniform grid. The reason for this is that the prescribed AMR criterion does not take into consideration the background solar wind, as it focuses on the CME itself. Therefore, the background wind is resolved to the base resolution of the heliospheric domain, which corresponds to the No AMR resolution from Table 2. Remark that the cell sizes are rather large here, which leads to such smooth profiles. Since the resolution is uniform everywhere in the uniform grid simulation, the solar wind is better resolved at 1 AU in that case. On the contrary, when comparing the spheromak CME interior in the different AMR simulations, more detailed structures are resolved with more AMR levels. The boundaries between opposite polarity regions are larger in the AMR level 2 run than in the AMR level 4 run. The shapes and structures lack details in the AMR level 2 case compared to the AMR level 3 and 4 cases. The overall structure of the magnetized CME is similar in all the simulations, as expected. When increasing the number of AMR levels, a higher resolution is obtained locally as the resolution is doubled for each higher level, which leads to better resolved gradients and local variations between the opposite polarity regions. This is especially important in studying the evolution of the magnetized CMEs in the heliosphere, determined by its interactions also with the magnetized background solar wind, and (in particular) the changes in the complex internal magnetic field structure. The uniform grid simulation result is mostly comparable to the AMR level 3 result and the features are finer and better resolved in AMR level 4 simulation result. This, in fact, is expected, since the longitudinal and latitudinal resolutions are not affected by radial grid stretching. In AMR level 3, the resolution is one level higher than medium uniform simulation in the latitude-longitude plane, but considering the lower radial resolution, AMR level 3 shows slightly smoother profiles at 1 AU.

Figure 4 shows the corresponding result for B_clt component at 1 AU in the EUHFORIA simulation with the medium-resolution uniform grid. A snapshot corresponding to the most similar time as for the shown Icarus results was chosen from the data set. The small differences in the CME interior is caused by the fact that the CME does not arrive exactly at the same time at L1 in Icarus and EUHFORIA; therefore, a slightly different cut through the 3D spheromak structure is captured at the given time in the simulation in EUFHORIA. We note that the background solar wind magnetic field is different from that in the Icarus results. There are more profiles present near the boundaries in the latitudinal direction. The features are well resolved in the background wind, similar to the Icarus uniform medium-resolution grid simulations. Also, the spheromak interior appears similar to that obtained with Icarus results, in terms of negative and positive portions and their alignments. The details and small structures are different in the interior (as expected), since very different numerical methods are used. Details on the EUHFORIA simulation setup with a spheromak CME can be found in Pomoell & Poedts (2018) and Verbeke et al. (2019). The fact that the CME locations are not at the same place in Figs. 3 and 4 is caused by the difference in reference frames (HEEQ in EUHFORIA, while corotating in Icarus). The CME is injected at the same location, but since we are plotting a given snapshot in both simulations, the whole domain is shifted, including the spheromak.

The AMR condition was limited for optimizing forecasting scenarios in these simulations. Thus, the refined area in the longitudinal direction was restricted to ±30° from the Sun-Earth line (as explained in Sect. 4.2.2). Usually, the restriction is also applied to the co-latitudinal direction, ±20° from the equatorial plane to further optimize forecasting simulations. In this case, however, the entire CME was resolved to higher refinement levels with the purpose to show the refined structure better.

Figure 5 shows the evolution of the radial velocity component (upper panel) and the number density (lower panel) at Earth. The black curve represents the corresponding OMNI 1 min data. The simulation was performed for a longer time interval than shown, but the dates were trimmed for better visibility, including only a short period before the arrival of the CME and a longer period until the CME has completely passed Earth. The corresponding EUHFORIA medium-resolution uniform grid simulation result is plotted with a blue dotted line and the Icarus medium- and high-resolution uniform grid simulations with orange and green solid lines. The Icarus AMR level 2, 3, and 4 simulation results are plotted with cyan, purple, and red lines, respectively. The jump in the radial velocity at around noon on July 14 corresponds to the shock arrival at Earth. The shock arrival seems to occur at different times for different simulations with Icarus, at least when the increase of the radial velocity is used as criterion. This is because the shock is more diffuse in the low-resolution simulations than in the higher resolution simulations. As a consequence, it occupies the same amount or more cells in the radial direction (e.g., three to four cells for AMR level 2, while mostly three cells in AMR level 3); moreover, the grid cells become smaller by a factor of 2 for each higher AMR level. As a result, when reading out the temporal evolution of the variables near Earth, the CME shock effectively reaches 1 AU earlier than in higher resolution simulations. This is merely because the shock is more diffuse, while in higher resolution simulations, it is more localized because it is captured with fewer grid cells and, moreover, these cells are smaller. Despite this diffusion effect, the shock strength is more or less independent of the resolution and the peak is reached at the same time in the different AMR level simulations. The shock profile is different in medium- and high-resolution uniform grid simulations, since its grid resolution is different in the entire computational domain, which affects the solar wind too, as can be seen in the figure before the arrival of the shock. The CME shock arrives slightly later in the EUHFORIA simulation than in the Icarus simulations. This can be explained by the difference of the numerical methods used in the two different heliospheric models. The radial velocity values increase from ~300 km s⁻¹ to ~650 km s⁻¹. The number density also jumps from ~5 cm⁻³ to ~30 cm⁻³. We also note that in the different simulations, the steepness of the curve is different, due to the different radial resolutions at 1 AU. The most steep profile is obtained with the AMR level 4 simulation, followed by the Icarus and EUHFORIA medium-resolution uniform grid simulations. As expected, the smoothest profile is modeled by the AMR level 2 simulation, because the radial resolution is the lowest in this case.

Figure 6 shows the time series at Earth for the uniform grid simulations with both EUHFORIA and Icarus together with AMR level 2, 3, and 4 simulation results. All magnetic field components are plotted. Data plotted with black line corresponds to in situ OMNI-1min observational data. As we note earlier in this work, the CME shock arrives slightly later in the EUHFO-RIA simulation than in Icarus simulations. This also leads to different arrival times for the magnetic clouds. The gradients in the time profiles of the plasma variables are the largest for the AMR level 4 simulation, which is showing stronger gradients than the ones obtained with the uniform grid simulation in Icarus. The arrival time is closer to the one in the OMNI data for the EUHFORIA simulation, but the strengths of the magnetic field components are slightly better modeled by Icarus AMR level 4 simulation. Overall, the profiles obtained in the Icarus simulations resemble the OMNI observed data slightly more closely, than with the EUHFORIA simulation. The simulations performed with Icarus are MHD simulations. The aim is to obtain the global picture of the evolution and arrival of the CME cloud and shock in heliosphere. Therefore, by default, our approach gives a macroscopic view. This leaves out microscopic effects such as turbulence and magnetic sheath formation, but it is useful and even necessary to predict the B_z component of the internal magnetic field of the CME and the shock arrival time and strength. Scolini et al. (2022, 2023) showed that the internal structure of the magnetic field of flux ropes changes significantly during their propagation and the magnetic complexity increases during this evolution. Therefore, higher resolution simulations are necessary, thereby enabling us to unravel the intricate internal plasma and magnetic structures and to model the evolution of the CMEs much more accurately. In general, the LFF-Spheromak CME model is successful in introducing the magnetized CME in the heliosphere, as the impact at Earth is mimicked relatively well, that is, the time profiles at Earth (L1) demonstrate a good comparison with the in situ measure profiles. More precisely, the CME shock arrival time, the strength of the magnetic field components and, thus, the direction of the magnetic field, are all in line with the observational data. In the last panel of Fig. 6, the difference between AMR simulation effects are shown clearly. The simulation result with AMR level 2 shows a smooth profile, omitting the variations in the total magnetic field strength before the actual magnetic cloud arrival. The simulation result with AMR level 3 already shows a slight increase in the total magnetic field component, similar to the observational data, along with some variation in the total magnetic field before the shock arrival. For both the AMR level 4 and uniform grid simulations in Icarus and EUHFORIA, this approach is found to model the magnetic cloud in the most similar way to the OMNI data. However, the peak value of the total magnetic field is the largest in the AMR level 4 simulation results, coming closer to the observed value than with the other, lower resolution simulations. We also conducted a simulation with AMR level 5, that is, with four additional levels of refinement, to check how the higher AMR level would perform. The result is given in Appendix A. The peak value in the AMR level 4 result is reached at the time at which the total magnetic field increases to its maximum in the observational data. In the AMR level 5 simulation results, the peak still coincides with the AMR level 4 peak, which suggests that the arrival of the magnetic cloud is modeled rather well in Icarus - and with even higher resolution simulations, the shock arrival would converge to the one seen in the observational data. Since, the peak values in the different magnetic field components are so similar to those in the observed data, the AMR level 5 simulation is not considered necessary, AMR level 4 should be sufficient for forecasting. We note that the considered CME is rather hard to model, since it occurred in an active phase of the Sun and the background wind has a complex structure. During its propagation, the CME interacts with a high speed stream, resulting in high gradients that need to be resolved in the computational domain. Therefore, the simulations with more AMR levels require longer wall-clock times than in less complicated cases. As can be observed in Fig. A.2, the gradients in the time profiles of the magnetic field components are slightly larger in the AMR level 5 results than in the AMR level 4 results, but there are no significant differences in the arrival time of the magnetic cloud nor in the strengths of the magnetic field components, to consider this additional refinement level necessary. The results, including the AMR level 5 simulation time series, are shown in Figs. A.1 and A.2.

Clearly, in the results of simulations using a higher number of AMR levels, we can distinguish more qualitative features, as compared to simulations results from cases with a lower number of AMR levels. In the total magnetic field strength, for instance, we can see that there is some structure being resolved in high AMR level simulation results just before the arrival of the magnetic cloud. The AMR level 4 simulation result shows stronger variations than the lower AMR level simulations. Looking at the AMR level 5 results in Fig. A.2, the profiles of the time evolution at L1 of the magnetic field components, as well as those of the total magnetic field strength, are also different from those obtained with the AMR level 4 simulation. This clearly indicates that higher AMR level simulations resolve smaller scale structures and introduce a more prominent magnetosheath in the simulation. However, elaborating more or deeper on this would require a totally different setting for the simulations, with strongly restricted AMR criteria, focused only on the area of interest in order to enable locally ultra high resolution simulations with very high AMR levels. This is rather important as the magnetosheath region is very turbulent and the oscillations in the observational data are not noise, but they do have a physical meaning. We consider this to be an interesting application for the Icarus heliospheric modeling tool in the future, as standard modeling with uniform grids do not allow for these kinds of studies due to CPU time restrictions. Overall, considering the chosen MHD approach, the synthetic (modeled) data with high resolution Icarus simulations meets the expectations when compared to the OMNI observational data. As a matter of fact, the profiles for the magnetic field components have a similar evolution in time as for the OMNI data, the shock steepness increases with higher resolution simulations and the peak value is reached at a time that is comparable to OMNI data. Also, we can spot in the number density profiles that with the higher resolution simulations, the CME is more localized and less diffuse or spread out, which makes it more comparable to the OMNI data.

Table 5 summarizes the simulation times required for each simulation on 1 node on the Genius cluster with 2 Xeon Gold 6240 CPUs at 2.6 GHz (Cascadelake), 18 cores each, at KU Leu-ven. The medium-resolution uniform grid simulation in EUH-FORIA took 18 h and 34 min, while with Icarus, a very similar setup took only 6 h and 52 min. The simulations with advanced numerical techniques are faster. The simulation on the stretched grid in combination with AMR level 2 takes only 11 min, while with AMR level 3 it takes 27 min and with AMR level 4 it needed 2 h and 34 min. Thus, the AMR level 4 simulation was 2.67 times faster than Icarus uniform grid simulation and 7.23 times faster than the EUHFORIA uniform grid simulation.

Fig. 2 Results in the ecliptic plane for the AMR level 4 simulation with the combined AMR criterion. The left panel shows the radial velocity values overlaid with the divergence of V. The middle panel shows the CME density values overlaid with divergence of V. The right panel shows the mesh superposed on the middle figure.

Fig. 3 Figures given at the 1 AU slice. The Bclt component of the magnetic field is plotted with the same scaling for all the figures for comparison purposes. AMR levels 2, 3, and 4 are used in the first three figures. The last figure shows results for the uniform medium-resolution grid simulation.

Fig. 4 Bclt in the EUHFORIA simulation, also at the 1 AU slice. The magnetic field is scaled similar to the results for the Icarus simulations shown in Fig. 3.

Fig. 5 Time series data at Earth. Time profiles for the radial speed and number density values are given for the uniform grid simulations in Icarus and EUHFORIA, together with Icarus simulations using AMR levels 2, 3, and 4.

Fig. 6 Time series data at Earth. Time profiles of all the magnetic field components are given for the uniform grid simulations in Icarus and EUHFORIA, together with Icarus simulations using AMR levels 2, 3, and 4.

6 Conclusions and outlook

In this study, a magnetized CME model represented by a LFFS was implemented in the new heliospheric wind and CME evolution model Icarus and validated. We considered the effect of advanced numerical techniques such as gradual radial grid stretching and solution-adaptive mesh refinement.

Icarus is a 3D MHD heliospheric wind and CME evolution model, covering the same domain as EUHFORIA, solving ideal MHD equations on uniform spherical or optimized grids with a finite volume numerical scheme. The coordinate frame in Icarus is co-rotating with the Sun, while in EUHFORIA, the equations are solved in the HEEQ coordinate system (in which the longitude of the Earth is fixed, but its latitude and radial distance vary). The model takes the output of the WSA-like coronal model in EUHFORIA and uses it as inner boundary conditions to determine the heliospheric wind. After an MHD relaxation phase, a stationary solar wind is obtained in the simulation, on which CMEs can be superposed by injecting them at the inner boundary, at 0.1 AU.

In this study, the implementation of a magnetized CME model is examined. The LFFS model was imported into Icarus from EUHFORIA via a coupling between the two architectures. The diagram in Fig. 1 shows how the linking was set up between the C++ and Fortran codes. This was achieved by generating a library from the C++ code, calculating the internal magnetic field components of the spheromak CME, which can be accessed from the Fortran code of Icarus. This way, the injected CME in Icarus is guaranteed to be identical to the one in EUHFORIA and, moreover, additional repetitive work is avoided. The LFF-Spheromak CME model in EUHFORIA is represented by the definition for the magnetic field that is divergence-free and force-free (unlike the Gibson-Low sphero-mak model). The equations describing such configuration are only briefly discussed in Sect. 3, since the original model was implemented in EUHFORIA by Verbeke et al. (2019).

Advanced numerical techniques were applied to optimize the computational grid in Icarus. The standard uniform spherical grid is available and it is similar to the one in EUHFORIA. Additionally, gradual radial grid stretching and AMR can be applied. The radial grid stretching preserves the aspect ratio between the width and the length of the cells. Therefore, the cells are not as deformed as in the uniform spherical grid simulations, near the inner and outer boundaries. On the other hand, the block-adaptive AMR guarantees high spatial resolution in the domain only in locations where this is necessary. The refinement con-dition(s) is (are) controlled by the user, leading to substantial flexibility in optimizing the computational grid on which the simulation is performed. Also, the maximum number of refinement levels can be specified. Whenever the implemented AMR criteria are met, the designated block in the computational grid is refined to the prescribed maximum refinement level. Moreover, when the conditions are not met any more, the grid is coarsened again, to avoid unnecessary CPU use.

In this study, the magnetized LFFS CME was injected in Icarus simulations using different grids. The goal of the work is to focus on the macroscopic effects, such as the global structure of the CME, its shape, complex internal magnetic field structure, and associated shock. The results were compared to the similar simulation with EUHFORIA, using the same mag-netogram and inner boundary conditions and with the same settings. Apart from the original uniform grid, several combinations of radial grid stretching and different AMR levels were considered. A real event study was performed in order to validate the model. The CME of July 12, 2012 was modeled with the pre-determined parameters available in the literature. The appropriate GONG magnetogram was chosen as inner boundary condition for the empirical WSA-like coronal model, to provide the plasma conditions at the inner boundary of the heliospheric model (at 0.1 AU). The AMR criteria were the same in all the AMR simulations and consisted of a combination of two refinement criteria, which yield high spatial resolution in both the CME interior and the CME shock. The refinement was also limited to the neighbourhood of the Sun-Earth line for optimization of forecasting simulations. The results were analysed using the time series data of the MHD quantities at Earth and 1 AU slice snapshots. Clearly, the surrounding solar wind is best resolved in uniform grid simulations both with Icarus and EUHFORIA, since this grid is uniform and the resolution is relatively high everywhere in the domain, independently of whether it is a CME region or the solar wind. In simulations using solution AMR, the grid is adjusted depending on the solution and the specified criteria. Since we are focusing on the CME, the resolution is only increased in and around the CME and the surrounding solar wind displays smoother profiles. On the other hand, the CME interior is far better resolved in the AMR simulations. With increasing the maximum number of AMR levels, the observed structures inside the CME are finer and the boundary between inversely polarized regions is narrower. More details can be seen in the AMR level 4 simulation result, compared to the AMR level 2, 3, and uniform grid simulation results. The time series profiles at Earth, as the CME passes through the position of Earth, give a good estimation of how well resolved different regions in the CME are at 1 AU for the different simulations. For assessing the performance of the combined AMR criterion, the radial velocity component and number density profiles are shown in Fig. 5 and the magnetic field components in Fig. 6. The shock arrival is best resolved by the AMR level 4 simulation, since the AMR criterion includes the shock regions in the domain. The profiles are steeper than in the other simulations, converging to the observed data steepness profile with higher numbers of AMR levels. The magnetic field components of the spheromak CME are also best resolved in with AMR level 4, since the combined criterion was also aimed at resolving the CME interior. The strength and profiles of the magnetic field components modeled by Icarus are slightly more similar to the OMNI data than the data modeled by EUHFORIA. The magnetic cloud arrives somewhat later in the EUHFORIA simulation, which is closer to the arrival data in OMNI data, but the strength of the magnetic field in the magnetic cloud is better estimated by the Icarus model.

In conclusion, the simulation results are very similar to those obtained with EUHFORIA, but require much less CPU time. The arrival time was slightly better modeled by EUHFORIA, but the profiles of the different variables considered (radial velocity, number density, and magnetic field components) are modeled slightly better in the Icarus simulation. The AMR level 4 simulation (with 3 refinement levels) produced the best results and when comparing to the AMR level 5 simulation results, with 4 refinement levels, it became clear that AMR level 4 is sufficient for modeling the spheromak CMEs with the purpose of accurately estimating the arrival time of the shock and the magnetic field component configurations. AMR level 5 produces slightly sharper shock results, but not significantly. Therefore, it is not considered worthwhile, in consideration of the longer computational times required.

Table 5 summarizes the wall-clock time each simulation required on 1 node of Genius cluster at KU Leuven. The uniform grid simulations with EUHFORIA and Icarus, having a relatively high spatial resolution everywhere in the domain, require ~18.5 and ~7h, respectively. The simulations combining radial grid stretching and AMR levels 2, 3, and 4 require ~0.2, 0.45, and ~2.5 h, respectively. The AMR level 4 simulation showed the best results, while being ~2.7 times faster than the uniform grid simulation in Icarus and seven times faster than the EUHFORIA simulation. The speed-up obtained by the comparison of these simulations is modest because the most broad refinement criterion is used, which refines both the shock and the CME interior. Simulations only focusing on the CME shock are much faster, since the CME interior area is not taken into consideration for a high resolution in that case. The AMR level 3 model data are the most similar to the uniform grid, standard-resolution simulations in Icarus and EUHFORIA. The AMR level 3 simulation is ~15 and ~41 times faster than the uniform grid with medium-resolution simulations in Icarus and EUHFORIA, respectively. The timings of the simulations are not constant, other real event simulations with different magnetic field configuration of the CME and different solar wind, might require different wall-clock times for computations. However, the case considered here is close to the solar maximum, leading to a quite variable solar wind and the modeled CME has high speed, causing strong shocks along its propagation; in addition, it interacts with a co-rotating interaction region. Therefore, this event study is a good estimation of how the model performs in complex cases.

The newly introduced magnetized CME model in Icarus enables in-depth studies of the effects of the interactions of the CMEs with the magnetic field of the background solar wind, for instance, deformation, deflection and flux, and plasma erosion. With the previous simple cone CME model, only the HD variables could be examined. Therefore, only the arrival time and the shock strength at Earth could be estimated. Since sphero-mak model includes the internal magnetic field configuration, the geo-effectiveness of the CME interaction with the Earth magnetic field can be studied, for instance, by calculating geomagnetic indices. The evolution of the complexity of the internal magnetic field of the CME can also be examined all the way from 0.1 AU to Earth and beyond. A multi-spacecraft study is considered with the magnetized CME model in the future, as it will allow for the model to be validated at multiple points, with the ability to compare the magnetic field can at different heliospheric locations. The combination of advanced numerical techniques with the magnetized CME model provides an opportunity to perform multiple simulations and investigate the effect of different parameters, without using up too many computational resources. This offers deeper insights into the evolution and propagation of the magnetized CMEs all the way to Earth, which will then help in forecasting such events more accurately and mitigate the possible damage more effectively.

Acknowledgements

This research has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 870405 (EUHFORIA 2.0) and the ESA project “Heliospheric modelling techniques” (contract no. 4000133080/20/NL/CRS). These results were also obtained in the framework of the projects C14/19/089 (C1 project Internal Funds KULeuven), G.0B58.23N and G.0025.23N (FWO-Vlaanderen), SIDC Data Exploitation (ESA Prodex-12), and Belspo project B2/191/P1/SWiM. The 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.

Appendix A Time series in Icarus

Appendix A.1 Time series with AMR level 5

Figures A.1 and A.2 show time series profiles at Earth, including the AMR level 5 simulations results. The colors denote the same simulations as in Fig. 5, where here, additionally, the AMR level 5 simulation result is plotted in purple. The AMR level 5 curve is almost identical to the AMR level 4 result, especially regarding the shock arrival time and the peak speed upon arrival. It displays more of a range in the number density values, as the CME is passing through the Earth. This feature is also present in the AMR level 4 result, but they are less pronounced.

Fig. A.1 Time series data extracted at the location of Earth. Radial speed and number density values are given for the uniform grid simulations in Icarus and EUHFORIA, together with Icarus simulations using AMR levels 2, 3, 4, and 5.

The magnetic field components also show more features with highly resolved grids. The AMR level 5 simulation result displays sharper profiles in the magnetic field temporal variations, visible in all the magnetic field components. However, the magnetic field strength time profile is the same as that modeled with the AMR level 4 simulation. The profiles are also very similar most of the part, the difference is mostly present upon the arrival of the magnetic cloud. In the time profile of the B_r-component, there is a more pronounced bump in the AMR level 5 curve before the large increase in the magnitude upon the magnetic cloud arrival. It is also visible in the AMR level 4 curve, but the profile is more smoothed out. In the time profile of the B_clt component, the AMR 4 and 5 curves are rather similar, with slightly higher gradients present upon the magnetic cloud arrival in the AMR level 5 result. In the time profile of the B_lon component, the maximum value of the component is modeled well only by the AMR level 4 and 5 simulations. Similarly, in the time profile of the total magnetic field strength, the strongest magnetic field is obtained in the simulations with the most AMR levels, stronger than in all other simulations, including the uniform grid Icarus and EUHFORIA simulations.

Fig. A.2 Time series data at Earth. All magnetic field components are given for the uniform grid simulations in Icarus and EUHFORIA, together with Icarus simulations using AMR levels 2, 3, 4, and 5.

References

All Tables

All Figures

	Fig. 1 Diagram of the procedure to import the spheromak file in Icarus used in EUHFORIA.
In the text

	Fig. 2 Results in the ecliptic plane for the AMR level 4 simulation with the combined AMR criterion. The left panel shows the radial velocity values overlaid with the divergence of V. The middle panel shows the CME density values overlaid with divergence of V. The right panel shows the mesh superposed on the middle figure.
In the text

	Fig. 3 Figures given at the 1 AU slice. The Bclt component of the magnetic field is plotted with the same scaling for all the figures for comparison purposes. AMR levels 2, 3, and 4 are used in the first three figures. The last figure shows results for the uniform medium-resolution grid simulation.
In the text

	Fig. 4 Bclt in the EUHFORIA simulation, also at the 1 AU slice. The magnetic field is scaled similar to the results for the Icarus simulations shown in Fig. 3.
In the text

	Fig. 5 Time series data at Earth. Time profiles for the radial speed and number density values are given for the uniform grid simulations in Icarus and EUHFORIA, together with Icarus simulations using AMR levels 2, 3, and 4.
In the text

	Fig. 6 Time series data at Earth. Time profiles of all the magnetic field components are given for the uniform grid simulations in Icarus and EUHFORIA, together with Icarus simulations using AMR levels 2, 3, and 4.
In the text

	Fig. A.1 Time series data extracted at the location of Earth. Radial speed and number density values are given for the uniform grid simulations in Icarus and EUHFORIA, together with Icarus simulations using AMR levels 2, 3, 4, and 5.
In the text

	Fig. A.2 Time series data at Earth. All magnetic field components are given for the uniform grid simulations in Icarus and EUHFORIA, together with Icarus simulations using AMR levels 2, 3, 4, and 5.
In the text