© The Authors 2024

1 Introduction

In the solar corona, rapid changes in the magnetic configuration can lead to the release of substantial quantities of magnetic energy, a phenomenon known as a solar flare (Forbes et al. 2006; Schmieder et al. 2007). In some cases, flares can also be accompanied by plasma ejections, such as filaments, jets (Schmieder et al. 1997, 2013). We refer to these as eruptive flares (as opposed to confined flares; Zuccarello et al. 2017). If the ejected structure extends spatially and is subsequently observed in the high corona in white light using a coronagraph, it defines a coronal mass ejection (CME). These spectacular phenomena can impact Earth’s magnetic environment and human technologies (Bothmer & Daglis 2007). Among other effects, satellite communications can be disrupted, GPS signals can be lost due to electronic variations in the ionosphere, and geomagnetically induced currents can damage the power grid and affect the erosion of oil pipelines.

The geoeffectiveness of a solar transient depends on its hydrodynamic and magnetic properties (Dumbović et al. 2015). Signatures of interplanetary coronal mass ejections (ICMEs) can be measured by different spacecraft such as ACE (Smith et al. 1998), the Parker Solar Probe (Fox et al. 2016), STEREO A and B (Kaiser & Adams 2007), the Solar Orbiter (Müller et al. 2020), Helios (Roberts et al. 1987), and Wind (Bougeret et al. 1995). For some events, in situ measurements reveal an interplanetary shock of the ICMEs, followed by a heated region called the sheath as a result of the accumulation of matter upstream of the expansion of the magnetic ejecta (ME; Kaymaz & Siscoe 2006). The ME, characterized by a strong magnetic signature and low temperature, is the direct signature of the twisted structures ejected from the solar corona. This structure is usually described as a flux rope (Démoulin et al. 2008). While interplanetary shocks alone can trigger magnetic disturbances on Earth, the strongest geomagnetic storms occur when the B_z component of the magnetic field of the ME, that is, the component of the magnetic field perpendicular to the equatorial surface, has a direction opposite to that of the magnetopause (Lugaz et al. 2016). Such a configuration allows for reconnection between the magnetic field lines of the two magnetic structures, which leads to an energy transfer towards the inner regions of the magnetosphere (Akasofu 1981; Dungey 1961).

ICMEs present complex signatures that reflect their interaction with the solar wind during their propagation in the heliosphere (Scolini et al. 2022). For example, using a 2.5D magnetohydrodynamics (MHD) simulation, (Zuccarello et al. 2011) found that the CME can be deflected toward the current sheet of a large streamer due to a discrepancy in the magnetic pressure and tension forces. Recently, Asvestari et al. (2022) found that the CME tilt evolves, depending on the strength and orientation of the ambient magnetic field. Interactions with the surrounding field can also lead to rotation of the CME (e.g., Liu et al. 2018; Manchester et al. 2017; Shiota et al. 2010). It should also be noted that the thermodynamic and magnetic properties of successive CMEs are directly affected by their interactions with each other (Scolini et al. 2020; Koehn et al. 2022).

Although multiple viewpoint reconstruction techniques are beginning to reconstruct the complex geometric and magnetic structure of a CME from the current set of observations (e.g., Rodari et al. 2018), the current limited number of satellites does not enable us to fully determine all the properties of ICMEs during their propagation (Démoulin 2010). Therefore, in addition to observations, a set of numerical codes has been developed with the aim of tracking the evolution of CMEs in the heliosphere. In a space weather context, their use is particularly relevant for attempting to predict the geoeffectiveness of a CME before it reaches Earth. Examples of such 3D MHD simulations include ENLIL (Odstrcil 2003), EUHFORIA (Pomoell & Poedts 2018), MS-FLUKSS (Singh et al. 2018), SUSANOO-CME (Shiota & Kataoka 2016), and ICARUS (Verbeke et al. 2022; Baratashvili et al. 2022), where the latter uses the versatile MPI-AMRVAC framework (Keppens et al. 2023) as its MHD solver, enabling radial grid stretching and solution adaptive mesh refinement. In these various simulation models, one or more CMEs propagate through a realistic background solar wind after their injection at the inner boundary of the domain (usually located at 0.1 AU). The predictions are obtained by extracting the temporal evolution of the plasma quantities (velocity, pressure, density, magnetic field, etc.) at different points in the heliosphere. The obtained time profiles depend, especially at L1/Earth, on the inserted CME model and its magnetic and thermodynamic properties (density, speed, magnetic flux, etc.). The characteristics of the injected CME can be deduced from remote sensing observations (e.g., Scolini et al. 2019; Maharana et al. 2022).

In this work, we are particularly interested in the 3D time-dependent magnetohydrodynamic (MHD) model of heliospheric wind and CME evolution, EUHFORIA (EUropean Heliospheric FORecasting Information Asset; Pomoell & Poedts 2018). Before this study, three CME models, each with its own capabilities and limitations, were implemented in the framework of EUHFORIA. The simplest is the cone model, which represents an unmagnetized plasma with a self-similar expansion (Xie et al. 2004). In EUHFORIA in combination with an energetic particle acceleration and transport model called PARADISE, the cone model can be used to model the interplanetary acceleration of low-energy protons during an energetic storm particle (ESP) event (Wijsen et al. 2022). However, this CME model is obviously not suitable for studying the evolution of the internal magnetic signature of an ICME. The second model implemented in EUHFORIA is the spheromak model, representing a CME with a global spherical shape but with an initially linear force-free magnetic field (Chandrasekhar & Kendall 1957; Verbeke et al. 2019). Although the use of this magnetized CME model is widespread (e.g., Verbeke et al. 2019, 2022; Scolini et al. 2019; Asvestari et al. 2022), as a result of its spherical geometry lacking CME legs, the spheromak model does not adequately model scenarios where Earth is impacted by the flanks/legs of ICMEs. To address this limitation, Maharana et al. (2022) implemented in EUHFORIA a flux rope model in which the CME has an extended flux-rope geometry including CME legs (e.g., “Flux Rope in 3D”, or FRi3D; Isavnin 2016). Its geometry is much closer to observations where CMEs are often observed to possess crescent-shaped shapes (Janvier et al. 2013). However, solving the complex equations defining the magnetic distribution of the model requires a substantial amount of computing time, which limits the use of this model in an operational space weather context.

In this study, we introduce the implementation of two new CME models with toroidal geometry in EUHFORIA, with the aim of overcoming the shortcomings of the current CME models. The toroidal geometry of the new CME models is simpler than that of the realistic FRi3D model. However, unlike the latter, we suggest that the torus models are better suited for operational space weather forecasting. Indeed, their magnetic and geometric configurations are based on straightforward analytical formulations. Consequently, they can be applied as time-dependent boundary conditions at 0.1 AU, making them much more time-efficient than the numerically determined FRi3D model.

On the other hand, although the new models do not remain connected to the Sun after injection, just like the spheromak model, they possess a geometry that is more consistent with the observations than the spheromak. Indeed, the toroidal models have a curved front and two ‘legs’, in contrast to the spheromak. Additionally, one of our models is based on a particular solution of the Grad-Shafranov equation. This equation is frequently employed to determine the cross-section of observed magnetic clouds (Möstl et al. 2009; Isavnin et al. 2011). Hence, we anticipate simulations with the new toroidal CME models would closely resemble the cross-sections of the observed ones. Furthermore, in this particular model, it is possible to adjust the ratio of the poloidal and toroidal fluxes. This feature is absent in other models and provides an added degree of freedom to more accurately predict the magnetic profiles observed at Earth.

Lastly, while both the CME models introduced here exhibit a toroidal geometry, they differ in the distribution of the internal magnetic field. This means we can expect varying magnetic field profiles in EUHFORIA, underscoring the rationale for implementing both models. It is important to note that this work focuses on validating the numerical implementation of these two models. In particular, we examine how the models’ distinct free parameters impact the time profiles recorded at Earth (extracted from the 3D time-dependent simulations). While addressing all potential configurations is not feasible, our research seeks to broadly demonstrate the potential of these two new models in EUHFORIA. A subsequent study will be devoted to directly comparing them with the FRi3D and spheromak models to highlight the strengths and weaknesses of all numerical CME models.

The paper is organized as follows: Sect. 2 introduces the two new torus CME models: one with the modified Miller-Turner solution (Vandas & Romashets 2015), and the other with the Soloviev solution (an analytical solution of the Grad-Shafranov equation; Soloviev 1975). In Sect. 3, we provide details on the implementation of the new models in the framework of EUHFORIA. The resulting thermodynamic and magnetic field profiles associated with the CME models in the heliosphere are discussed in Sect. 4. In this section, the influence of various free parameters of the CME models is also investigated. We specifically discuss the impact of initial radial speed (cf. Sect. 4.3), magnetic field strength (cf. Sect. 4.4), density (cf. Sect. 4.5), as well as minor radius (cf. Sect. 4.6). To conclude, Sect. 5 summarizes the main findings and potential implications for future research.

2 CME models

We implement in EUHFORIA two different CME models. In this section, we provide an overview of the two CME models: (1) the modified Miller-Turner solution (Vandas & Romashets 2015), and (2) the Soloviev solution, an analytical solution of the Grad-Shafranov equation (Soloviev 1975).

2.1 The modified Miller-Turner solution

The magnetic configuration of this model was originally derived by Miller & Turner (1981) and represents a constant-a force-free field in a toroidal geometry given by:

$$B_{{\rho }_l}=\frac{B_0}{2{\alpha }_{\text{MT}}{R}_0}{J}_0\left({\rho }_l{\alpha }_{\text{MT}}\right)\sin {\theta }_l,$$(1)

$$B_{{\varphi }_l}=B_0\left(1-\frac{{\rho }_l}{2R_0}\cos {\theta }_l\right)J_0\left({\rho }_l{\alpha }_{\text{MT}}\right),$$(2)

$$\begin{array}{l}{B}_{{\theta }_l}=B_0\frac{\cos {\theta }_l}{2{\alpha }_{\text{MT}}{R}_0}{J}_0\left({\rho }_l{\alpha }_{\text{MT}}\right)\\ -B_0\left(1-\frac{{\rho }_l\cos {\theta }_l}{2R_0}\right)J_1\left({\rho }_l{\alpha }_{\text{MT}}\right),\end{array}$$(3)

where (ρ_l, ϕ_l, θ_l) are the local toroidal curved cylindrical coordinates, α the minor radius of the torus, R₀ the major radius and B₀ the strength of the toroidal component in the center of the torus (cf. Fig. 1). J₀ and J₁ are the Bessel functions of the first kind of order 0 and 1, respectively. The local angles ϕ_l and θ_l range from 0 to 2π, while ρ_l covers the interval [0, a]. Hence, this CME model has a fixed circular poloidal cross-section unlike the other implemented model which contains parameters that determine the shape of the poloidal cross-section (cf. Sect. 2.2). To ensure that the magnetic field is completely poloidal $\text{(i.e., }{B}_{{\varphi }_l}=0\text{ ) }$ at the boundary of the flux rope and that the magnetic field line is confined in the torus, the value α_MT is such that aα_MT is the first root of J₀(ρ_lα_MT):

$${\alpha }_{\text{MT}}\approx C_{\alpha }\frac{2.41}{a},$$(4)

where C_α is defined as the chirality of the magnetic field and takes value +1 (−1) for right (left) handedness. The Miller-Turner solution tends to the Lundquist solution (Lundquist 1951) when R₀ → ∞.

This solution has one main disadvantage: the solenoidal condition ∇ • B = 0 is not exactly fulfilled, but only approximately. To mitigate the effect of the non-solenoidal scenario, in this work, we use a modified Miller-Turner (mMT) solution defined by Romashets & Vandas (2003) such that:

$$B=\nabla \times \frac{B_{\text{MT}}}{{\alpha }_{\text{MT}}},$$(5)

where B_MT is the original Miller-Turner solution (cf. Eqs. (1)− (3)). The new magnetic field can be expressed as:

$$B_{{\rho }_l}=B_0\frac{R_0-2{\rho }_l\cos {\theta }_l}{2{\alpha }_{\text{MT}}{R}_0\left(R_0+{\rho }_l\cos {\theta }_l\right)}{J}_0\left({\rho }_l{\alpha }_{\text{MT}}\right)\sin {\theta }_l,$$(6)

$$B_{{\varphi }_l}=B_0\left(1-\frac{{\rho }_l}{2R_0}\cos {\theta }_l\right)J_0\left({\rho }_l{\alpha }_{\text{MT}}\right),$$(7)

$$\begin{array}{l}{B}_{{\theta }_l}=B_0\frac{R_0-2{\rho }_l\cos {\theta }_l}{2{\alpha }_{\text{MT}}{R}_0\left(R_0+{\rho }_l\cos {\theta }_l\right)}{J}_0\left({\rho }_l{\alpha }_{\text{MT}}\right)\cos {\theta }_l\\ -B_0\left(1-\frac{{\rho }_l}{2R_0}\cos {\theta }_l\right)J_1\left({\rho }_l{\alpha }_{\text{MT}}\right).\end{array}$$(8)

In the mMT configuration, the magnetic field is fully divergence-free and the degree of twisting varies from zero at the center to an infinite value at the outer boundary.

In Fig. 2, the top panels shows the magnetic field distribution in the mTM solution. In this figure, the center of the torus is at the origin (0, 0, 0) of the local coordinate system, assuming that the rotation axis aligns with the z-axis. The major radius, R₀, is equal to 10 R_⊙, while the minor radius is a = 4 R_⊙. The magnetic field in the circular cross-section has a maximum shift toward the torus hole. According to Vandas & Romashets (2015), the shift increases when the aspect ratio, R₀/a, decreases. In theory, the magnetic field is fully force-free only for large aspect ratios (R₀/a ≫ 1). However, Vandas & Romashets (2015) found that the force-free condition in the original Miller-Turner solution (cf. Eqs. (1)−(3)) is held much better than in the mMT model for low aspect ratios (≈2). In this study, we have chosen the mMT solution, which is unique in being completely divergence-free, thus avoiding numerical instability and unphysical solutions. It is worth noting that there is also the Tsuji solution, which describes a magnetic configuration in the same toroidal geometry but is force-free for all aspect ratios (Tsuji 1991). That being said, the expressions for its magnetic field components contain infinite series expansions and the coordinate system transformations required for the implementation in EUHFORIA are too complex for a straightforward and quick numerical implementation.

Fig. 1 Local toroidal curved cylindrical coordinates associated with the Miller-Turner CME model. In the local Cartesian system centered on Ol, the z-axis coincides with the symmetry axis of the torus.

Fig. 2 Normalized magnetic field strength in the two CME models implemented in EUHFORIA. The left panels show a sample of magnetic field lines in a torus with a minor radius a = 4 R⊙ and a major radius R0 = 10 R⊙. The right panels show the dimensionless magnetic field strength B/B0 in a poloidal cross-section. The “+” markers indicate the center of the torus. From top to bottom: magnetic field configuration in the modified Miller-Turner model, in the Soloviev model with αS = 1 and with αS = 10. When the αS parameter is too high, the magnetic field is mainly toroidal.

2.2 The Soloviev solution

2.2.1 Analytical formulation

The second CME model originally corresponds to a global magnetic confinement topology in a tokamak configuration (Goedbloed 1984). This configuration can be described by a particular cylindrical system with R representing the distance from the center to the symmetry axis, Z the vertical coordinate, and ϕ the toroidal angle (cf. Fig. 3). Considering we assume axisymmetry along ϕ, and that the static equilibrium solution satisfies the following MHD equations

$$J\times B=\nabla p,J=\nabla \times B,\text{\hspace{1em} and \hspace{1em}}\nabla \cdot B=0,$$(9)

we find that all solutions are described by the non-linear Grad-Shafranov equation

$$\text{ΔΨ}-\frac{1}{R}\frac{\partial \text{Ψ}}{\partial R}=R\frac{\partial }{\partial R}\left(\frac{1}{R}\frac{\partial \text{Ψ}}{\partial R}\right)+\frac{{\partial }^2\text{Ψ}}{\partial Z^2}=-I{I}^{\prime }-R^2{p}^{\prime }.$$(10)

Here, Ψ is the poloidal flux, which ranges from 0 on the magnetic axis to Ψ₁ on the plasma boundary. The Grad-Shafranov equation, Eq. (10), depends on two arbitrary flux functions: the stream function I(Ψ) of the poloidal current and the pressure p(Ψ) (cf. Sect. 16.2, Goedbloed et al. 2019, for the derivation).

The magnetic configuration used in our work is obtained from an exact analytical solution of the Grad-Shafranov equation, called the Soloviev equilibrium. This particular solution imposes linear profiles for I² and p(Ψ) (i.e., for the case when the right-hand side of Eq. (10) is a constant). The dimensionless flux ψ = Ψ/Ψ₁ is then given by:

$$\psi ={\left[X-\frac{1}{2}ϵ\left(1-X^2\right)\right]}^2+\left(1-\frac{{ϵ}^2}{4}\right)\left[1+ϵ\tau X(2+ϵX)\right]{\left(\frac{Y}{\sigma }\right)}^2,$$(11)

which depends on the dimensionless poloidal coordinates

$$X=\frac{R-R_0}{a},\text{\hspace{1em} and }$$(12)

$$Y=\frac{Z}{a},$$(13)

as well as the inverse aspect ratio ϵ = a/R₀ (≪ 1 in asymptotic expansions), the triangularity τ, and the elongation (or ellipticity) σ. Using the dimensionless quantity ψ, the boundary condition becomes ψ(X, Y) = 1 on the outer boundary of the plasma cross-section.

Fig. 3 Particular cylindrical system associated with the Soloviev solution. The color lines are flux isocontours such as ψ = 1 (in red), ψ = 0.7 (in green) and ψ = 0.2 (in blue).

2.2.2 Geometry

In the Miller-Turner model, the poloidal cross section is circular (cf. the top panels, Fig. 2). However, it is possible to change the shape of the poloidal cross-section in the Soloviev solution by varying the triangularity τ and the elongation σ parameters. Figure 4 shows the isocontour ψ = 1 in the Soloviev model as a function of triangularity and elongation.

The shape of the flux contour can also be modified by changing the triangularity τ (cf. Fig. 4, top panel). The triangularity parameter allows one to pass from an elliptic flux contour when τ = 1 (and σ > 1) to a triangular flux contour whose direction depends on the sign of the triangularity (cf. Fig. 4, top panel). For τ = −0.5 and τ = 1 (cf. Fig. 4, top panel), the isocontour ψ = 1 is closed and delineates the limits of the torus. However, the flux contour is open for τ = −2 and τ = 2.5, which means that the torus is not properly delineated. Analytically, this flux contour shape reflects the presence of hyperbolic points located at:

$$x_s=-\frac{1}{ϵ}\left[1-\frac{1}{\tau }\right]\text{, }$$(14)

$$y_s=\pm \frac{\sigma }{ϵ\tau }\sqrt{\frac{1}{2}\frac{1+{ϵ}^2\tau }{1-\frac{{ϵ}^2}{4}}},$$(15)

that are present in the domain if τ < −1/[ϵ(2 + ϵ)] or τ > 1/[ϵ(2 − ϵ)]. To obtain a closed flux contour, the value of τ must be restricted within these limits. For example, with an inverse aspect ratio of 0.4 as in Fig. 4, the triangularity must lie within the range [−1.04, 1.56] to obtain a closed contour. The values −2 and 2.5, being outside this range, result in an open isocontour for ψ = 1 in both cases, as shown in Fig. 4.

The elongation σ does not change the size in the toroidal plane of the torus, but affects the cross-sectional shape of the torus along the Z-direction, as shown in Fig. 4. Indeed, the outer isocontour moves from a near-circular cross section when σ = 1 to an ellipsoidal cross section when σ > 1 (cf. Fig. 4, bottom panel). The maximum height reached in the Z-direction is σa. For further information on the geometry of the torus, we refer to Goedbloed (1984).

Fig. 4 Isocontours of ψ = 1 in the X − Z plane for the Soloviev solution with different geometries. The upper panel illustrates the impact of the triangularity parameter τ on the isocontour ψ = 1, while the lower panel, characterizes the influence of the elongation parameter σ for a torus with a = 1 R⊙ and R0 = 2.5 R⊙, centered at (0,0,0). The marker “+” indicates the center of the torus, while the dot marker corresponds to x = R0. In the top panel, σ is kept constant at 1.4 and in the bottom plot, τ is constrained to 0.

2.2.3 Magnetic field configuration

Using the analytical solution for the flux (cf. Eq. (11)), the magnetic field components in the particular cylindrical system can be determined such as:

$$B_R=-\frac{B_0}{{\alpha }_{\text{S}}}\frac{ϵ}{1+ϵX}\frac{\partial \psi }{\partial Y},$$(16)

$$B_Z=\frac{B_0}{{\alpha }_{\text{S}}}\frac{ϵ}{1+ϵX}\frac{\partial \psi }{\partial X},$$(17)

$$B_{\varphi }=\frac{B_0}{1+ϵX}{\left[1-2\frac{ϵ\psi }{{\alpha }_{\text{S}}^2}\left(Aϵ-\frac{B}{2}\right)\right]}^{1/2},$$(18)

where,

$$A=2\left[1+\frac{1-\frac{{ϵ}^2}{4}}{{\sigma }^2}\right]\text{, }$$(19)

$$B=4ϵ\left[1+\frac{1-\frac{{ϵ}^2}{4}}{{\sigma }^2}\tau \right],$$(20)

and

$$\frac{\partial \psi }{\partial X}=\left[X-\left(1-X^2\right)\right](1+ϵX)+\left(1-\frac{{ϵ}^2}{4}\right)[2ϵ\tau (1+ϵX)]\left(\frac{Y^2}{{\sigma }^2}\right),$$(21)

while

$$\frac{\partial \psi }{\partial Y}=\left[1-\frac{{ϵ}^2}{4}\right][1+ϵ\tau X(2+ϵX)]\left(\frac{2Y}{{\sigma }^2}\right),$$(22)

and α_S is related to the inverse of the total poloidal flux through the plasma and defined as

$${\alpha }_{\text{S}}=\frac{a^2{B}_0}{{\Psi }_1}.$$(23)

Here, B₀ is a constant magnetic field strength used to scale the different magnetic field components (cf. Eqs. (16)−(18)).

Contrary to the modified Miller-Turner CME model (cf. Sect. 2.1), the Soloviev model allows for varied magnetic field distributions for the same poloidal flux Ψ₁, depending on the value of the parameter α_S. In the Miller-Turner CME model, the magnetic helicity, which signifies the winding of magnetic field lines relative to one another, is a constant parameter. However, the Soloviev solution introduces flexibility in setting helicity by adjusting the parameter α_S (cf. Eq. (23)). Indeed, the local magnetic fields B_R and B_Z (cf. Eqs. (16) and (17), respectively) exhibit variations proportional to 1/α_S while the poloidal field B_ϕ presents a more intricate dependence (cf. Eq. (18)). As a result, the adjustment of α_S modifies the ratio between the toroidal and poloidal fields.

To visualize this, Fig. 2 displays the magnetic field distribution in both the mMT solution and the Soloviev model with α_S = 1 and α_S = 10. A distinct variance can be observed in the distribution of magnetic field lines within the Soloviev model between the α_S = 1 and α_S = 10 cases. In the latter scenario, the high α_S value leads to the dominance of the magnetic field B_ϕ resulting in an almost purely toroidal magnetic field. It is worth noting that, in a tokamak configuration, when the toroidal magnetic field varies as B_ϕ = I(Ψ)/R, this component is always dominant (Goedbloed 1984).

Moreover, the different magnetic helicities within the Soloviev and mMT models indicate their distinct magnetic field distributions. As seen in the right panels in Fig. 2, the isocontours of constant B in the mMT model are circular, while in the Soloviev model they appear as near-vertical bands. It is also worth noting that the maximum of magnetic field is not located at the center, but rather at the boundary in the Soloviev model, which is not usually the case in the magnetic configurations popularly used to model CMEs (Lundquist 1951; Chandrasekhar & Kendall 1957). Given these disparities in magnetic field distributions between the Soloviev and the mMT models, we aim to validate the implementation of both models in EUHFORIA and investigate the differences in their profiles in the heliosphere.

3 Numerical implementation of Miller Turner CME in EUHFORIA

3.1 EUHFORIA

The two magnetic configurations used are implemented as CME models in the physics-based 3D magnetohydrodynamic simulation of the heliosphere wind and CME evolution, EUHFORIA (European Heliospheric Forecasting Information Asset; Pomoell & Poedts 2018). For space weather operational purposes, EUHFORIA can be used to track the propagation of one or more CMEs in the heliosphere and predict the temporal evolution of different plasma quantities (e.g. density, velocity, magnetic field) at different positions in the inner heliosphere.

EUHFORIA is divided into two main parts. The first one is the semi-empirical Wang–Sheeley–Arge model (WSA; McGregor et al. 2011; van der Holst et al. 2010) that inputs a synoptic magnetogram to provide the plasma quantities at 0.1 AU needed to establish the background solar wind in the helio-sphere. The magnetogram can come from different sources, for example from the Global Oscillation Network Group (GONG; Harvey et al. 1996) or from the SDO’s Helioseismic and Magnetic Imager (HMI; Schou et al. 2012). In this coronal part, the 3D coronal magnetic field is first computed from a potential field source surface (PFSS) extrapolation (Altschuler & Newkirk 1969). The magnetic field is then extended to 0.1 AU using the Schatten Current Sheet model (Schatten et al. 1969). The solar wind speed at the 0.1 AU boundary is obtained as a function of the flux tube expansion factor, f, and the distance, d, of the foot point of the flux tube to the nearest coronal hole. From this solar wind speed υ_sw = υ(f, d), the density and the temperature at the outer boundary are computed using empirical relations (cf. Pomoell & Poedts 2018, for more details).

The second part of EUHFORIA is a 3D time-dependent magnetohydrodynamic simulation of the heliosphere. The ther-modynamic and magnetic quantities, computed from the coronal model, are self-consistently determined between 0.1 AU and 2 AU in an MHD relaxation phase in a uniform grid. This evolution is determined by solving the set of ideal MHD equations with a polytropic index equal to 1.5, using a finite volume numerical architecture and a constrained transport scheme ensuring that the solenoidal condition is always satisfied up to machine accuracy.

EUHFORIA can be used to track the propagation of a CME and its possible interaction with the background solar wind by injecting a CME model at the inner boundary as described in the following section (cf. Sect. 3.2).

3.2 CME implementation

The CME injection in EUHFORIA takes place at the inner boundary of the heliospheric model, that is, at 0.1 AU. Starting from a predetermined initiation time, at each time step increment a mask is computed to identify the points on the grid where the injected CME intersects the spherical boundary surface at 0.1 AU. At all these points, the solar wind speed, magnetic field, density, and temperature are substituted with those of the CME. In the subsequent time step, the CME center is advanced with the initially defined purely radial speed, υ_r0, and the 3D mask region is updated as the solution moves through the spherical boundary.

Initially, the center O_l of the torus is located at the position (t_d, θ_i, ϕ_i) in spherical coordinates, where t_d is the distance from the center of the domain O (0,0,0), θ_i is the initial colatitude and ϕ_i is the longitude (cf. Fig. 5). These three coordinates are free parameters of the models and must be set prior to the simulation.

At the intersection point, the components of the solar wind magnetic field that are replaced are in global spherical coordinates (R, θ, ϕ). The equivalent global Cartesian system is (X, Y, Z) centered at O. Hence, it is necessary to transform the magnetic field computed according to Eqs. (6)–(8) or (16)–(18) in a local toroidal coordinate system centered at O_l into magnetic fields (B_r, B_θ, B_ϕ) in the global spherical system centered at O. The specifics of the base change are explained in Appendix A.

Different criteria are used to determine the intersection points for the Soloviev and the mMT models. In the case of the Soloviev solution, the flux Ψ, is computed (cf. Eq. (11)). All points on the boundary where the flux is less than or equal to 1 are considered to be within the CME. As mentioned in Sect. 2.2, depending on the elongation, triangularity, and aspect ratio, the shape of the isocontour Ψ = 1 changes. In some cases, the isocontour is not closed (cf. Fig. 4). Numerically, this will lead to the computation of a magnetic field that is not physical everywhere at the border. Therefore before using the Soloviev model, it must be analytically checked whether condition Ψ < 1 accurately defines a closed contour with the desired geometric parameters (cf. Sect. 2.2.2).

For the Miller-Turner solution, the local toroidal radius ρ_l (cf. Eq. (A.1)) is computed; all points on the boundary where this radius is less than the minor radius of the torus are within the CME. This approach cannot be used for the Soloviev solution, as in that model the poloidal cross section is not (always) circular (cf. Sect. 2.2).

In EUHFORIA, the two toroidal CME models are distinguished by their magnetic field configurations and the free geometrical parameters that define CME kinematics and kinetics. The CME parameters for both models are summarized in Table 1. In the following subsections, we cover various parameters to highlight the differences between the two implemented CME models and how they affect the EUHFORIA simulation obtained at Earth. To limit the number of simulations, we do not discuss the impact of the initial position (i.e., the angles θ_i and ϕ_i) and the internal temperature of the CME (T). We also skip the triangularity and elongation parameters of the Soloviev model as they only influence simulations with high resolution and with large minor radius that enable clearly differentiating between a circular and a “triangular” cross section.

Fig. 5 Relation between the global Cartesian coordinate system of EUHFORIA (X, Y, Z), and the local Cartesian system in which the two CME models are defined (Xl, Yl, Zl). The distance from the center td, the colatitude θi, and the longitude ϕi define the initial position of the center of the torus, Ol. The vector OOl, where O is the center of the domain, corresponds to the Xl axis of the local Cartesian coordinate system.

4 Thermodynamic and magnetic profiles

4.1 Theoretical profiles

Before analyzing the thermodynamic and magnetic profiles obtained in EUHFORIA, it is insightful to examine how the components of the magnetic field evolve in the absence of solar wind. Figure 6 shows the evolution of the components of the magnetic field (in the Geocentric Solar Ecliptic System, GSE) for the two magnetic configurations at the level of three virtual satellites crossed by a torus with a speed of 600 km s⁻¹. It is essential to note that the magnetic field depends on the geometry of the torus, notably through the aspect ratio: R₀/a. However, Fig. 6 provides a general idea of how the tilt of the torus and the in situ measurement point influence the components of the obtained magnetic fields. The minor and major radii of the torus are a = 5 R_⊙ and R = 10 R_⊙. For the Soloviev solution, the elongation, σ, equals 1, the triangularity, σ, is 0.5 and the parameter α_S also equals 1. The center of the torus evolves in the X direction.

In the first row of Fig. 6, the torus crosses the three virtual satellites without any tilt (ω = 0). When the virtual satellite is directly in the propagation path of the torus (i.e., the blue satellite in Fig. 6) there are two sections corresponding to the front and back of the torus, separated by the torus hole. The two parts have the same size, and the distribution of the magnetic field inside depends on the CME model. In both cases, the B_x component is zero, and the B_ɀ component is symmetric in the two sections, with a profile change sign close to the center.

However, the direction of B_ɀ variation is opposite between the two CME models. Another difference is that the B_y component, anti-symmetric between the two magnetic structures, is maximum close to the torus hole for the Soloviev solution, while it reaches its maximum at the magnetic axis for the mMT solution.

If the virtual satellite is in the direction of propagation, it will pass through the torus hole. However, it is also possible to obtain a single magnetic structure by shifting a virtual satellite by an angle, ϕ = 10°, so that it does not pass through the torus hole (cf. Fig. 6, third and last column). As before, the two CME models have different magnetic field distributions, but a double change in the sign of the B_ɀ component is observed in both cases. A change in sign is also observed for the B_y component. Moreover, unlike the case where the virtual satellite is in the direction of propagation, the B_x component is not null. It has a significant influence and is even predominant near the structure’s middle in the Soloviev solution (cf. the top-mid panel in Fig. 6). Also, it is important to note that crossing either the right or left part of the torus (relative to the torus hole) merely changes the sign of B_x in both models.

Asvestari et al. (2022) found that the interaction of the CME with the ambient field in EUHFORIA can modify the inherent CME tilt even in the absence of actual CME rotation in the heliosphere. The middle and bottom rows of Fig. 6 show the magnetic field profiles in the two CME models obtained with a tilt, ω, of 45° and −45°, respectively. The first observation is that at the three virtual satellites, there are two magnetic structures, indicating the passage of a part of the torus hole in all cases. However, the profiles are completely different from those obtained at the same location but without tilt. In particular, symmetry and anti-symmetry is no longer observed for the magnetic field components between the two parts of the torus when the virtual satellite is in the direction of propagation (cf. Fig. 6, second column).

It should also be noted that the total force of the magnetic field (the black dashed and solid line in Fig. 6) at the level of the three satellites is identical between the torus with a tilt of ω = 45° and the torus with an opposite tilt of ω = −45°. However, the distribution of the magnetic field is different, as are the dominant components. Even with tilt, the transition from a virtual satellite offset by an angle ϕ = 10° relative to the direction of propagation to a virtual satellite offset by an angle ϕ = −10° only changes the sign of the component B_x.

Finally, in all configurations, the maximum reached by the total magnetic field is always higher in the Soloviev model than in the mMT model. In the latter, B₀ is the maximum at the magnetic axis, whereas in the Soloviev solution, B can reach approximately 1.7 × B₀ close to the boundary hole. However, this maximum depends on the parameter α_S.

In conclusion, the magnetic field theoretically presents complex profiles that depend on the chosen model, its geometry, its tilt, and also the measurement point. Therefore, we can expect different profiles in EUHFORIA between the two models as well (cf. Sect. 4.2).

Fig. 6 Theoretical profiles of magnetic field components are obtained by advancing a torus at a speed of 600 km s−1 through three virtual satellites. The first column presents an isosurface Ψ = 1 for the Soloviev solution, with three torus inclinations: ω = 0° (top panel), ω = 45° (middle panel), and ω = −45° (bottom panel). The second column represents, for these three configurations, the evolution of magnetic field components (in GSE) at the level of a virtual satellite (blue) placed directly in the direction of torus propagation. The third column corresponds to the passage through the offset virtual satellite (orange) at an angle ϕ = 10° relative to the propagation direction. The fourth column corresponds to the offset green satellite at an angle ϕ = −10°. The dashed lines correspond to the magnetic field components obtained from the Soloviev CME model, whereas the solid lines correspond to the mMT solution.

4.2 Modelled profiles in EUHFORIA

4.2.1 Numerical set-up

For all the simulations we conducted with EUHFORIA for the present paper, the solar wind was reconstructed from the GONG magnetogram obtained at 01:04 UT on November 4, 2015. This date was chosen arbitrarily to illustrate the potential use of EUHFORIA. Hence, as mentioned in Sect. 1, we are not trying to model the CMEs that occurred near this date, but we are interested in how the thermodynamic and magnetic profiles obtained during the propagation of the CME in the solar wind vary depending on the initial parameters.

The mesh used is uniform in all directions with an angular resolution of 4° in the latitudinal and 2° for the longitudinal directions and 256 cells in the radial direction. This resolution is similar to that used by Maharana et al. (2022). Convergence tests were carried out by doubling the number of cells in all three spatial directions, but this did not change our main results while drastically increasing the computation time.

Before studying the impact of different parameters, we carried out a series of simulations that are used as a reference in the following sections. In these simulations, the center of the CME evolves from 12:04 PM on November 11 and from the initial position O_l (t_d = 6.5 R_⊙, θ_i = 81°, ϕ_i = 0°) with a radial speed of 600 km s⁻¹. The torus has a minor radius a = 5 R_⊙, a major radius R_⊙, and a null tilt ω = 0. For the Soloviev solution, the triangularity is τ = 0.5 and the elongation σ = 1. Thus, the geometric properties are identical to the tori presented in the previous section (cf. Sect. 4.1). We chose this elongation and triangularity because they allow a comparable cross section between the mMT and Soloviev models while maintaining a closed Ψ = 1 isocontour (cf. Sect. 2.2).

At the points of intersection between the inner boundary of EUHFORIA and the CME (cf. Sect. 3.2), the injected temperature is constant at 0.8 MK. The density is also constant at n₀ = 1 × 10⁻¹⁷ kg m⁻³. As explained in Sect. 4.5, using both a constant density and temperature results in a constant pressure within the torus. In fact, this contradicts the linear pressure profile we specified to obtain the Soloviev solution (Soloviev 1975). However, we have deliberately chosen to violate this assumption (cf. Sect. 4.5 for more details) and use here only the magnetic distribution of the Soloviev equilibrium.

According to Fig. 6, for the same value of B₀, the magnetic field strength is higher in the Soloviev solution than in the mMT model. Additionally, due to differing magnetic field distributions and slightly different geometries, the two CME models used do not have the same magnetic flux. Bearing this in mind, we still decided to use the same initial magnetic field strength, B₀, in both models, in order to align as closely as possible with the same initial parameters.

Finally, the simulations were run on the wICE cluster of the Vlaams Supercomputer Centrum, a Belgian institution specializing in high-performance computing¹. They have utilized 4 nodes of this supercomputer, with 36 tasks per node, implying a total of 144 parallel processes. Additionally, each CPU allocated to this simulation has 2 GB of memory.

4.2.2 Reference cases

According to the initial magnetic field and the geometrical parameters outlined earlier, Fig. 7 illustrates the evolution of speed, density, magnetic field, and plasma beta (β) at Earth in EUHFORIA for various CME configurations. We simulated three CMEs using the Soloviev solution with differing α_S (α_S = 1, −1, and 10). Additionally, we modeled two CMEs with the mMT solution with C_α = 1 and −1, respectively.

In Fig. 7, we observe that using either the mMT or the Soloviev solution as the CME model results in notably different thermodynamic and magnetic profiles. Specifically, although the CMEs have the same initial speed and the same initial magnetic field strength, B₀, the speed obtained at Earth is higher in simulations with a Soloviev CME than with a mMT CME. For example, the CME computed from the Soloviev model with α_S = 1 reaches a maximum speed that is approximately 9% higher than the mMT CME with C_α = 1. Consequently, the shock arrival time is dependent on the chosen CME model. Additionally, the Soloviev CMEs also demonstrate a higher magnetic field strength, |B|, than the mMT CMEs, as anticipated from the theoretical profiles (cf. Fig. 6). This suggests that the speed of the CME is not solely related to the initially imposed speed but also to the initial magnetic field strength (cf. Sect. 4.3).

At the very end of the simulation, after November 8 in Fig. 7, the speed increases in the two fastest simulations. This increase is a numerical artifact arising from the solver’s challenges in managing steep speed gradients as detailed in Sect. 4.2 of Linan et al. (2023). This can be partially mitigated by enhancing the spatial resolution.

As discussed in Sect. 4.5, even if the CMEs are inserted with the same density, the peak density measured at Earth differs in various simulations. The temporal evolution of the parameter β also depends on the CME model used. However, for all CMEs, an increase in the parameter β is observed at the moment of arrival at Earth, followed by a decrease to values lower than those of the solar wind before the event. The Soloviev CME with α_S = 10 exhibits a distinct evolution from other cases, as there is a second increase in the parameter β a few hours after the first increase, which might indicate the passage of the back part of the CME or some perturbation in the wake of the CME. Except for this case, the profile of β is consistent with real in situ measurements on Earth. Indeed, Regnault et al. (2020), by performing superposed epoch analyzes of 20 yr of ACE data, found that during an active period of the solar cycle, the parameter β increases in the sheath but is lower than before the event during the traversal of the magnetic ejecta.

Considering the magnetic field profiles in Fig. 7, the results align with the theoretical profiles presented in Sect. 6, suggesting that they are influenced by the specific CME model implemented. All CME models feature a non-zero B_x component, implying that Earth is not perfectly in the CME propagation direction (cf. Fig. 6). However, its contribution is minor compared to the B_y and B_ɀ components, whose trends vary from one CME to another.

To explain the trend of the B_y component, Fig. 8 shows the distribution of the longitudinal magnetic field B_lon at two different times for the CME derived from the mMT solution with positive chirality (Figs. 8a,b) and from the Soloviev solution (Figs. 8c,d). The longitudinal magnetic field B_lon in the EUHFORIA coordinate system is directly related to the B_y component in the geocentric solar ecliptic system (B_y = −B_lon). Initially (cf. Fig. 8a), B_lon is positive at the front part of the CME and negative at the back in the mMT model. For the Soloviev CME model with α = 1 (cf. Fig. 8c), the B_lon distribution is almost similar but with opposite signs, as anticipated by the theoretical profiles. We can also note that the longitudinal field in the torus hole is slightly positive in the mMT solution and negative in the Soloviev solution. The magnetic field within the torus closely resembles that of the CME front part. We suggest that it is a consequence of the solver ensuring the magnetic field’s solenoidality and absence of significant gradients.

Upon the arrival of the magnetic ejecta at Earth (cf. Figs. 8b,d), the torus expanded in all directions following a self-similar expansion. The initial toroidal geometry is no longer discernible. Furthermore, ahead of the front part of the CME derived from the mMT model, we can see a thin region where the longitudinal field is negative. Since the longitudinal field of the CME is positive in the front part of the CME (cf. Fig. 6), we associate this negative field with the presence of a sheath that develops as the CME expands (Siscoe & Odstrcil 2008).

With respect to the different time profiles, we suggest a transition between the sheath and the magnetic ejecta at approximately 2015-11-06T19:00, (cf. the vertical green line in Fig. 7) for the mMT simulation with C_α = 1 (the full red curve). In this simulation, the passage through the sheath is characterized by an increase in the B_y-component. On the contrary, the passage through the magnetic ejecta is marked by an extended duration, where B_y remains negative. These two distinct structures result in two main peaks in the total magnetic field strength, |B|. However, for the Soloviev CME cases, we observe mainly an increase in both the B_y-component and |B|. Since the longitudinal field in the front part of the CME is negative, there is no sign change in the B_y-component in this case. Consequently, the presence of a sheath is not as evident as in the mMT CME. We suggest that the rise in plasma β could indicate the passage through a thin sheath (Masías-Meza et al. 2016). However, solely based on the thermodynamic and magnetic profiles shown in Fig. 7, it is not feasible to pinpoint the exact time of the transition between the sheath and the magnetic ejecta in the Soloviev simulations. It is also worth noting that determining this boundary can be complicated due to magnetic reconnection occurring at various points along the ME, which smooth the transition between the different magnetic structures, as discussed by Romashets & Vandas (2003). Moreover, the sheath forms due to turbulence which is not included in the simulation model. Further studies beyond the scope of this work are thus necessary to gain a better understanding of the development and characteristics of sheaths in EUHFORIA. The primary concern is that in EUHFORIA, the sheath only develops after the CME injection at 0.1 AU. By not considering the interactions between the solar wind and the CME closer to the Sun, the size and properties of the sheath can deviate significantly from expectations, potentially leading to inconsistencies when compared to in-situ measurements at L1.

Finally, based on the theoretical profiles (cf. Fig. 6), we expected to see a significant variation of the B_y field at the simulation end, indicating the passage of the torus’s back. In theory, this contribution should be antisymmetric with respect to the torus front (i.e., having the same amplitude but with the opposite sign). The back of the torus, characterized by a positive longitudinal field in the Soloviev case, is visible in Fig. 8d. However, its amplitude is so weak that it is not detectable in the time evolution of the component B_y in Fig. 7.

Figure 9 shows the evolution of the colatitudinal magnetic field, B_clt, for the same two CMEs. In the EUHFORIA coordinate system, B_clt equals −B_z in GSE. Initially (cf. Figs. 9a,c), we distinguish the injected torus through the inner boundary, with a change in the sign of the B_z component as predicted by the theoretical profiles (cf. Fig. 6). Like B_lon, the magnetic field within the torus hole follows the trend of the last point of the CME injected before it. Later in the propagation (cf. Figs. 9b,d), the inversion line in the front part is clearly visible in both cases. As a consequence of this magnetic field distribution, the time evolution of the B_z component presents a significant change in sign. However, as with the B_y component, the temporal profile of B_z does not allow us to distinguish the influence of the back part of the CME (cf. Fig. 7). Further analysis is necessary to understand the exact causes of erosion and/or dispersion of the back part.

Theoretically, the sign of the B_z component is controlled by the sign of chirality, that is, the sign of C_α or α_S. In EUHFO-RIA, changing the sign of chirality indeed reverses the profile of B_z (cf. Fig. 7). However, we also notice that chirality influences the amplitude of the different components, particularly in the simulations with an mMT CME. For example, the maximum of the total magnetic field strength, |B|, is 13.4 nt for the mMT CME with C_α = 1, while |B| reaches only 9.0 nt for the mMT CME with C_α = −1. This leads us to infer that the distribution of the magnetic field affects the nature of the interactions with the ambient solar wind, thus influencing the predicted profiles.

Finally, as expected, increasing the α_S parameter in the simulations with the Soloviev solution allows for a more longitudinal field. Therefore, in the case of the CME with α_S = 10, the predominant magnetic contribution is the B_y component. In the forthcoming sections, we focus exclusively on CMEs with α_S = 1 to facilitate comparisons with CMEs derived from the mMT solution, which show a substantial contribution from the B_z component.

Fig. 7 Magnetic and thermodynamic profiles as function of time obtained at Earth using EUHFORIA with the two CME models. Top to bottom: speed, density, total magnetic field strength, Bx, By, Bɀ components in GSE coordinates, and plasma beta (β in log scale). The solid lines correspond to simulations where the CME is based on the mMT solution, while the dashed lines are obtained from simulations with the Soloviev solution as CME model. The green vertical lines mark the presumed sheath boundaries in the mMT simulation with Cα = 1.

Fig. 8 Visualization of the longitudinal magnetic field component in the EUHFORIA heliosphere domain. For both CME models, two initial mass densities were used: n0 = 1 × 10−16 kg m−3 (panels a and b for the mMT CME model and panels c and d for the Soloviev CME model) and n0 = 1 × 10−16 kg m−3 (panels e and f for the mMT CME model and panels g and h for the Soloviev CME model). The geometrical, magnetic, and thermodynamic parameters of the CMEs are described in Sect. 4.2.1. The left panels (a, c, e, g) show the evolution of the longitudinal field, Blon, on March 4 while the right panels (b, d, f, h) show the simulation on November 6. Each panel is divided into an equatorial cross section, a meridional cross section passing through the Sun and the Earth, and a Mercator projection at 0.1 AU from the Sun. The visualization displays the positions of the different planets in the simulation domain, the location of the current sheet, and the interplanetary magnetic field connecting the Sun and the planets, obtained from Parker spirals.

Fig. 9 Same as Fig. 8 but for the co-latitudinal magnetic field, Bclt.

4.3 Impact of the initial velocity

As outlined in Sect. 3.2, the center of the torus moves in a purely radial direction at each time step, according to the following equation:

$$t_d(t)=t_d(0)+{\upsilon }_{r0}t,$$(24)

where t_d(t) denotes the distance from the domain’s origin at time t, and υ_r0 represents the initial radial speed. Consequently, the torus experiences a translation in the direction of vector OO_l (cf. Fig. 5).

At EUHFORIA’s inner boundary, the solar wind speed is altered to model the CME’s passage. There are two primary strategies to set the injected velocity. The first strategy is to compute the velocity at each point on the torus cross-section, considering its translation and the imposed radial velocity at the center. Therefore, each point in the torus injected at the boundary has a unique velocity, embodied by a uniform motion. This approach is used in the spheromak model (Verbeke et al. 2019). The second strategy models a self-similar evolution at the boundary, as demonstrated in the FRi3D model (Maharana et al. 2022). At 0.1 AU, all intersection points between the torus and the EUHFORIA boundary have only a purely radial speed, equivalent to υ_r0.

The numerical implementation accommodates both methods. However, numerous tests have revealed that simulations exhibit greater stability when using self-similar evolution, as opposed to configurations where each point can have poloidal and toroidal speeds. In the latter scenario, the injected CME is more likely to expand unnaturally upon entering the simulation domain, in turn making the simulation unstable. Therefore, in this work, we consider cases where the injected speed at the boundary is purely radial.

Figure 10 shows the impact of the initial radial speed, υ_r0, on the profile obtained in EUHFORIA. The geometrical and thermodynamic properties of the injected tori are identical to the reference cases presented in the previous section (cf. Sect. 4.2.2) with C_α = 1 for mMT and α_S = 1 for Soloviev models. Only the initial speed varies from 400 km s⁻¹ to 1000 km s⁻¹.

As observed in Fig. 10, a higher initial speed results in a higher shock speed at the time of arrival at Earth. A higher speed, that is, υ_r0 = 1000 km sv, results in a stronger shock, and thus a slightly more abrupt increase in velocity compared to slower cases. However, the maximum speed reached is in all cases lower than the injected speed, indicating that the magnetic structure slows down during its propagation. This deceleration may be caused by a drag force acting on the CME by the surrounding solar wind (cf. Sect. 4.4).

In Fig. 10, the amplitude of the density, magnetic profiles, and plasma β differs according to the initial speed of the CME. The higher the latter, the higher/lower the local minimums/maximums. This result is consistent with the in situ measurement of the ACE spacecraft. Indeed, by comparing a set of fast CMEs and a set of slow CMEs, Masías-Meza et al. (2016) found that the median value of |B| at Earth is higher in the high-speed group than in the slower group (cf. also Regnault et al. 2020). In EUHFORIA, one of the possible reasons for the difference in amplitudes is that, depending on the initial speed, the Earth does not cross the CME at the same location. In Figs. 8 and 9, for example, one can see that the longitudinal and colatitudi-nal fields are not uniformly distributed in the simulation. Thus, a change in the impact position can result in a modification of the amplitude of the different quantities (e.g Lepping et al. 1990).

Furthermore, the time profiles in Fig. 10 seem to indicate substantial compression of the magnetic structure against the strong shock in the cases where υ_r0 = 1000 km s⁻¹. In EUHFORIA, the torus hole does not have the same speed as the torus itself. When the radial speed is high, the back of the torus arrives faster in the domain and catches up more quickly with the front of the torus. The back part of the torus, having had less time to expand, finds itself compressed between the surrounding solar wind and the back. As a result of this compression, the maximum and minimum local values reached by the B_x, B_y, and B_z components at the front of the CME are higher when the initial speed is increased.

The sheath, delineated by the first bump of |B| and B_y in the mMT CMEs, is also more compressed when υ_r0 = 1000 km s⁻¹ due to the fast-moving magnetic ejecta from behind. In contrast, in slow CMEs, the magnetic field and plasma have enough time to approach a near-pressure balance. As a result, an expansion of the sheath can be expected to mirror the trend followed by the CME (Démoulin & Dasso 2009).

The trends of the magnetic field components are also impacted by the change in initial speed (cf. Fig. 10). For the B_x component, there are some differences between the various cases because the magnetic structure is not traversed at exactly the same position at 1 AU. The most significant difference is related to the temporal evolution of B_y. In the fastest cases, that is, with υ_r0 = 1000 km s⁻¹, we see a positive bump in the mMT CME and a negative one in the Soloviev CME at the end of the magnetic ejecta traverse. These bumps are not observed with the slower CMEs (i.e., with υ_r0 = 400 km s⁻¹ and υ_r0 = 600 km s⁻¹). We conclude that these new trends are a trace of the back of the torus, which was characterized by a positive B_y in the mMT model and a negative B_y in the Soloviev solution (cf. Fig. 6).

According to the theoretical profiles (cf. Fig. 6), a change in the sign of the B_z component is expected in the back part of the CME. However, this change in sign occurs only at the very end of the CME. Therefore, we suggest that this part of the CME does not cover a sufficiently large region to be clearly perceptible in the B_z profile and that it is eroded by the interaction with the solar wind in the wake of the front part of the CME.

To circumvent issues of unnatural magnetic field compression and potential inaccuracies in predictions, CMEs might be better modeled using a half-torus, that is, by taking only the front part of the torus into account. In practical terms, this involves injecting only half of the toroidal CME model by halting its passage across the EUHFORIA boundary when the center of the torus reaches 0.1 AU. This approach, which will be explored in a second paper, is employed by Singh et al. (2020) who simulated the 12 July 2012 event by introducing half of a spheromak into a realistic background solar wind. Contrary to a full torus, the injected CME geometry then closely resembles a flux rope featuring a curved front with two legs. However, it is worth noting that in Singh et al. (2020), the CME is initially superimposed onto the solar wind rather than being gradually introduced into the domain via the inner boundary.

Fig. 10 Same as Fig. 7 but the only difference between the different cases computed from the same CME model is the initial radial speed. For both CME solution, the initial speed is υr0 = 300 km s−1, υr0 = 600 km s−1 and υr0 = 900 km s−1.

4.4 Impact of the initial magnetic field strength

To modify the amplitude of the magnetic field profiles, one possibility is to increase the magnetic field B₀, and consequently the magnetic flux. We used the same geometry as our reference CMEs (cf. Sect. 4.3) and performed simulations with a weaker magnetic flux using a magnetic field that is half as strong as in our reference cases, that is, B₀ = 5 × 10⁻⁷ T. We also ran simulations with a magnetic field three times larger, that is, B₀ = 3 × 10⁻⁶ T.

The first significant difference among the simulations is the computation time. The reference simulations require 23 min and 24 s for the mMT case and 37 min and 39 s for the Soloviev one, utilizing the cluster described in Sect. 4.2.1. The Soloviev CME, despite having almost the same geometry, magnetic, and kinetic properties, thus requires more computation time when B₀ equals 1 × 10⁻⁶ T. Since the equations governing the two models are not identical, the time required to calculate them numerically can also vary.

Furthermore, the computation time is intrinsically related to the Courant-Friedrichs-Lewy (CFL) condition, which ensures the stability of the numerical method. Given that the Soloviev CME has slightly more intense temporal profiles than those observed with the mMT CME, it is inferred that more significant gradients exist within the numerical domain. These pronounced gradients can result in swift variations in the domain, increasing the effective propagation speed of disturbances. To maintain stability according to the CFL condition, a rise in the propagation speed leads to a decrease in the time step size. Consequently, a smaller time step requires more steps to cover an equivalent time interval, thus extending the overall computational duration.

The reduction of the value of B₀ by a factor 2 decreases the computation time to 15 min and 19 s for the mMT CME and to 13 min and 34 s for the Soloviev CME. Since the field B₀ is only a constant used to scale the magnetic field obtained in both CME models, it has no influence on the calculation time of the CME properties during its injection at 0.1 AU. Therefore, if changing the value of B₀ affects the simulation time, it is because the resolution of the MHD equations requires more or less time, depending on the amplitude of the magnetic field, once the CME has entered the domain.

Indeed, increasing B₀ requires more computational resources, leading to a considerable increase in computation time. It becomes 27 min and 12 s for Soloviev and 4 h and 14 min for mMT. In the latter cases, the extended simulation time indicates the solver’s difficulty in converging when the gradients between the CME’s magnetic field and the solar wind’s properties are significantly high at the inner boundary.

Examining the thermodynamic and magnetic profiles in Fig. 11, we observe that the larger the magnetic field B₀, the higher the speed. Pal et al. (2018), by studying a set of 30 CMEs, also established a positive correlation between their magnetic field strength and their velocity, as determined from white-light images at 10 R_⊙. They found that, statistically, CMEs with large B₀ expand faster than CMEs with a lower magnetic field strength. This is consistent with our observations. However, their statistical analysis did not definitively determine whether this speed difference is a result of a higher ejection speed influenced by the initial magnetic field strength during the eruption or whether it can also be attributed to physical processes during propagation. In our simulations, the latter scenario is illustrated, as our CMEs are all inserted at 21.5 R_⊙ with the same initial velocity.

Finally, the speed at Earth depends not only on the initial speed (as discussed in Sect. 4.3), but also on the magnetic properties of the model. This phenomenon has also been observed in other numerical simulations. For example, by studying the propagation of different flux ropes in the corona, Linan et al. (2023) discovered that the maximum speed of the CME at 0.1 AU is directly proportional to the initial magnetic field strength (see also Regnault et al. 2023).

The influence of the magnetic field on the propagation speed stems from the balance of forces acting on the CME, which can be described by the following equation:

$$F=m\frac{{\text{d}}^{2}R}{\text{d}{t}^2}$$(25)

$$=F_{\text{Lorentz }}+F_{\text{drag }},$$(26)

where F is the total force, m the CME mass, t the time, and R is the position of the center of mass. The expanding CME experiences two main forces, the first being the Lorentz force F_Lorentz (e.g., Chen & Kunkel 2010; Sachdeva et al. 2015), which is proportional to J × B where J is the current density. This force allows for the radial expansion of the CME. The second, F_drag is the drag force that tends to restrain the CME, explaining why it gradually slows down during propagation (Chen & Kunkel 2010; Subramanian et al. 2012; Sachdeva et al. 2015, 2017). Increasing the magnetic field B₀ alters the balance between these two forces by strengthening the Lorentz force, resulting in the acceleration of the CME. It is also important to note that the expansion of the CME depends not only on its initial velocity and magnetic properties but also on the characteristics of the background solar wind, particularly through the drag force. For instance, by examining a sample of more than a hundred ICMEs, Vršnak & Žic (2007) suggested a linear least squares fit between the Sun-Earth transit time of interplanetary coronal mass ejections and both their speeds and the solar wind velocity.

Increasing the initial magnetic field B₀ (cf. Fig. 11) logically increases the maximum amplitude reached by the total magnetic field |B|. However, in these cases, the maximum magnetic field is weaker than the initial magnetic field B₀. This means that the magnetic field carried by the CME decreases during its propagation through the heliosphere (Leitner et al. 2007; Liu et al. 2005; Winslow et al. 2015).

The increase in B₀ also modifies the trend of |B| (cf. Fig. 11). For the CMEs computed from the mMT solution, we see two distinct local maxima in the cases where B₀ = 5 × 10⁻⁷ T and B₀ = 1 × 10⁻⁶ T, while there is just a slight change in the slope of the magnetic field increase when B₀ = 3 × 10⁻⁶ T. This change is related to the transition between the sheath and the magnetic ejecta (cf. Sect. 4.2.2). As explained in Sect. 4.3, increasing the speed leads to compression of the sheath. This compression is also visible in the profile of plasma β, which is broader in the cases where the initial magnetic field is the lowest.

In both models, the trends of the components B_y, B_z are not significantly impacted by changes in the initial magnetic field (cf. Fig. 11). Only the amplitude and width of the profiles are modified. In the extreme case where B₀ = 3 × 10⁻⁶ T in the mMT model, the B_x component follows the same trend as the B_y component, which was not discernible with a weaker initial magnetic field. As the speed is higher than in other cases, the Earth does not cross at the front of the CME, resulting in a significant B_x component as observed in the theoretical profiles (cf. Fig. 6). However, at each moment, the value of the B_x field remains lower than those of the components B_y and B_z.

It is also worth noting that even though the CMEs with a magnetic field strength of B₀ exhibit speeds close to those obtained in the previous section, where the CMEs were inserted with an initial speed of υ_r0 = 1000 km s⁻¹ (cf. Fig. 10), the various profiles do not display the same characteristics. Especially, all the profiles with B₀ = 3 × 10⁻⁶ T are broad, in contrast to the compression of profiles observed when υ_r0 = 1000 km s⁻¹ in Fig. 10. Furthermore, in Sect. 4.3, the presence of a bump in the B_y component is mentioned, indicating the traversal of the back part of the torus when υ_r0 = 1000 km s⁻¹. We suggest that these features are not observed when B₀ = 3 × 10⁻⁶, since the insertion of the torus still occurs at the speed of υ_r0 = 1000 km s⁻¹. The acceleration of the CME happens gradually during its propagation. Additionally, due to the high magnetic field and the predominance of the Lorentz forces, the front part of the torus is more free to expand spatially and is no longer compressed between the shock and the back part. Finally, in Fig. 11, it is also seen that in both models, increasing the initial field leads to a decrease in the maximum density measured at the Earth’s level. This observation is explained in the following section (cf. Sect. 4.5).

Fig. 11 Same as Fig. 7 but the only difference between the different CMEs is the initial magnetic field strength B0. For CMEs computed from the Turner solution, B0 is equal to 5 × 10−7 T, 1 × 10−6 T and 3 × 10−6 T; while it is 1 × 10−7 T, 2 × 10−7 T and 6 × 10−7 T in the Soloviev cases.

4.5 Impact of the initial mass density

The density and temperature of the solar wind at EUHFORIA’s inner boundary were adjusted to model a CME’s transit. The newly determined density and temperature are uniform at all intersection points, in line with the CME models already implemented in EUHFORIA (e.g., Scolini et al. 2019).

Under the ideal MHD paradigm, upon which the solver is based, imposing a constant density and temperature results in a constant pressure within the torus. However, to achieve the Soloviev equilibrium derived from the Grad-Shafranov equation, the pressure profile is assumed to be linear (cf. Sect. 2.2). Despite this, we chose to adopt constant density and temperature to facilitate numerical computation and comparison between the two CME models employed in our study. Moreover, even if we insert a CME in perfect equilibrium, the equilibrium is immediately broken upon insertion anyway because the external environment changes drastically.

Figure 12 shows the different thermodynamic and magnetic profiles obtained in EUHFORIA using three different initial mass density: n₀ = 1 × 10⁻¹⁸ kg m⁻³, n₀ = 1 × 10⁻¹⁷ kg m⁻³, and n₀ = 1 × 10⁻¹⁶ kg m⁻³. In this figure, it is apparent that all profiles are affected by the change in density.

As mentioned earlier (cf. Sect. 4.4), the total magnetic field of the CME decreases during its propagation. Similarly, due to the expansion of the CME, the density also decreases (Lugaz et al. 2005). Furthermore, when the density is high, specifically n₀ = 1 × 10⁻¹⁶ kg m⁻³, the inertia of the plasma inside the CME counteracts both magnetic expansion and dissipation. At the moment of impact with Earth (cf. Figs. 13e,f), the CME is characterized by a global large, compact magnetic structure with high density and magnetic field. In the meridional cross-section, we can notably distinguish several regions of high intensity corresponding to the sheath, the front part of the torus, the torus hole, and the rear of the torus. These different regions have “pancake” shapes, resulting from the pressure gradients between the CME and the ambient solar wind (Riley & Crooker 2004). Moreover, we can add that since the magnetic field is high within the torus, so is the Lorentz force, leading to a significant velocity (cf. Sect. 4.3).

Conversely, when the density is lower, that is, n₀ = 1 × 10⁻¹⁸ kg m⁻³ (cf. Figs. 13a,b), the CME undergoes over-expansion after its insertion into the domain, leading to a significant decrease in both density and magnetic field. As a result, as illustrated in Fig. 12, the amplitude of the magnetic field, density, and velocity profiles is higher in CMEs with a density of 1 × 10⁻¹⁸ kg m⁻³ compared to CMEs with lower initial mass densities. These trends are also detailed in Maharana et al. (2022), who obtain similar results by varying the initial density of the FRi3D CME model in EUHFORIA.

Altering the initial density of the CME also leads to a change in the distribution of the longitudinal magnetic field (cf. Figs. 8f,h) and the colatitudinal magnetic field (cf. Figs. 9f,h). This results in a modification of the trends of B_y and B_z, as shown in Fig. 12. Specifically, a prolonged phase where B_y is positive in the mMT profile and negative in the Soloviev profile is observed towards the end of the magnetic ejecta when the initial density is high, that is, n₀ = 1 × 10⁻¹⁶ kg m⁻³. As mentioned in Sect. 4.3, and according to the theoretical profile (cf. Fig. 6), this corresponds to the passage through the rear of the torus.

In simulations with the CMEs having the highest initial densities, that is, n₀ = 1 × 10⁻¹⁶ kg m⁻³, the shape of the density profile also reflects the clear traversal of different magnetic structures composing the different parts of the torus as shown in the two dips in plasma beta in Fig. 12 (7 Nov. and 8 Nov.).

Finally, looking at the density distribution in Fig. 13, it is also worth noting that in all the simulations, the maximum density is located at the nose of the CME. Therefore, if Earth were to traverse the flank of the CME, the density encountered would be lower than that if it had passed through the front part. Since the impact position depends on the speed of the CME, the amplitude of the density profile obtained varies according to the initial speed and magnetic field imposed (cf. Figs. 10 and 11).

Fig. 12 Same as Fig. 7 but the only difference between the different cases computed from the same CME model is the initial mass density, n0. For both CME solution, the initial mass density is n0 = 1 × 10−18 kg m−3, n0 = 1 × 10−17 kg m−3 and n0 = 1 × 10−16 kg m−3.

4.6 Impact of major and minor radii

The torus must begin its traversal of the boundary just after the insertion time. In other words, before the insertion time, there are no intersection points between the CME and the inner boundary. To ensure this, the following equation must be satisfied:

$$t_d+a+R_0=21.5R_{\odot }.$$(27)

Therefore, the size of the CME is limited as its maximum value, a + R₀ = 21.5 R_⊙.

Figure 14 shows the time evolution profiles obtained in EUHFORIA for CMEs with varying minor radii a. In these instances, the major radius is fixed at R₀ = 10 R_⊙. Consequently, by modifying the minor radius, we have also adjusted the initial distance from the center of the torus to the center of the domain, that is, t_d, according to Eq. (27).

As the minor radius increases, so does the spatial size of the CME. Therefore, the larger a, the wider the profiles obtained in Fig. 14 will be. On the other hand, increasing the minor radius without altering the initial magnetic field B₀ results in a rise in magnetic flux. Therefore, the magnetic profiles for the case of a = 9 R_⊙ are comparable to those obtained when B₀ = 3 × 10⁻⁶ T in Sect. 4.4. As previously described, the amplitude of the various magnetic profiles and the speed of the CME increase with the minor radius. However, the computation time also increases. It goes from 23 minutes and 24 seconds for the reference mMT CME with a = 5 R_⊙ to 1 hour and 11 minutes for the mMT CME with a = 9 R_⊙. The simulation using a Soloviev model with a = 9 R_⊙ requires approximately three times the computation time as the reference simulation with a = 5 R_⊙.

All different CMEs shown in Fig. 14 have the same initial density n₀. As the minor radius increases, so does the volume of the torus and, hence, the total mass of the CME injected across the boundary. Therefore, the maximum of the obtained density profile is higher for a = 9 R_⊙ than for a = 3 R_⊙, as described in Sect. 4.5. However, it should be recalled that the density profile’s amplitude depends on the impact position between the CME and Earth.

Finally, altering the torus geometry, for example, by adjusting the minor radius, changes the magnetic and thermodynamic properties of the CME and thus the resulting profiles, as described in previous sections. It is important to note that, numerically, the solver (i.e., EUHFORIA) allows the use of a minor radius close to the major radius, that is, R₀/a ≈ 1. Theoretically, the magnetic field is force-free only for large aspect ratios (R₀/a ≫ 1). Therefore, we can conclude that the force-free assumption is not a limiting criterion and that it is not necessary in EUHFORIA to restrict the range of minor and major radii to satisfy R₀/a ≫ 1.

Fig. 13 Distribution of scaled mass density at the moment of arrival at Earth as a function of the initial CME density, n0, and the CME model used. The right panels show simulations with the mMT CME model, while the left panels show simulations with the Soloviev CME model. Three different initial mass densities were used, from top to bottom: n0 = 1 × 10−18 kg m−3 (panels a and b), n0 = 1 × 10−17 kg m−3 (panels c and d) and n0 = 1 × 10−16 kg m−3 (panels e and f). The different legends related to the EUHFORIA vizualisation tool are described in Fig. 8.

5 Conclusion

In order to introduce new CME models in EUHFORIA to fill the gaps in the models currently used, we implemented two magnetic configurations with a toroidal geometry. The first CME model is derived from the modified Miller-Turner (mMT) model (cf. Sect. 2.1). The second corresponds to the Soloviev equilibrium, an analytical solution to the Grad-Shafranov equation (cf. Sect. 2.2). Although both models possess a toroidal geometry, the Soloviev CME model offers more free parameters. Indeed, unlike the mMT model, where the twist is fixed, it is possible to modify the magnetic helicity of the Soloviev solution. On the other hand, the mMT model has a circular cross section, while the poloidal aspect can be modified in the Soloviev CME model by adjusting the triangularity and the elongation (cf. Sect. 2.2.2). The distribution of the magnetic field within the torus is also different in both models, which justifies the implementation of two configurations (cf. Sect. 4.1).

After presenting these differences analytically, we show the numerical implementation in EUHFORIA (cf. Sect. 3.2). This implementation is identical in both models. The torus crosses the inner boundary of EUHFORIA with an initially fixed radial speed. At each intersection point between the torus and the boundary, the speed, magnetic field, density, and temperature of the solar wind are modified to model the passage of the CME.

To illustrate the numerical implementation, we performed a series of simulations in which the CME is injected into a realistic solar wind background derived from the GONG magnetogram obtained at 01:04 UT on November 4, 2015. We examined the impact of different initial parameters on the magnetic and thermodynamic profiles at Earth in EUHFORIA. From these profiles, our main findings were as follows:

• Both models are stable and can be used in EUHFORIA. The magnetic profiles obtained at Earth are consistent with the distribution of the magnetic field injected across the boundary.

• Once within the domain, the CME radially expands with a speed dependent on the radial velocity of the torus center during torus insertion. When the initial velocity is high (approximately 1000 km s⁻¹), we observed a compression of the magnetic structure between the back of the torus and the shock. This leads to narrow profiles. Examining the distribution of the magnetic field, we noticed the development of a sheath ahead of the mMT CME, the thickness of which also depends on the speed of the magnetic ejecta.

• Increasing the magnetic field leads to an increase in the Lorentz force, and therefore increasing the speed of the CME. This also allows us to obtain broad profiles with high amplitudes. However, the computation time is directly related to the intensity of the magnetic field.

• In both models, the CME expansion is directly related to the initial mass density of the torus. When the density is too low (here n₀ = 1 × 10⁻¹⁸ kg m⁻³), the CME experiences over-expansion, leading to dispersion of both density and magnetic field. In contrast, when the density is high (e.g., n₀ = 1 × 10⁻¹⁶ kg m⁻³), the torus maintains a compact structure. It is then possible to discern the influence of the back of the torus on the magnetic profiles.

• The amplitude of the local peaks and the duration of the magnetic disturbances can be modified by adjusting the geometry of the torus. The greater the minor radius, the greater the injected magnetic flux and total mass. Therefore, the amplitude of the magnetic field and density can be altered without changing the initial magnetic field strength and the initial mass density.

Finally, EUHFORIA enables the use of these two new toroidal models. Despite their relatively simple analytical formulation, it is possible to achieve complex profiles, depending not only on the initial parameters of the CME but also on the impact position with Earth.

Now, the next step is to use these two new toroidal CME models to predict the geoeffectiveness of real CME events. Similar to Scolini et al. (2019), the properties of the torus, such as its size, magnetic flux, and speed, must be deduced from various observations to ensure the most realistic simulation possible, thus providing the most accurate predictions (cf. also Maharana et al. 2022). This kind of study will allow us to conclude which model is most suitable for modeling real events. In fact, our study alone cannot determine if one model is superior to the other, since both models provide different magnetic profiles and react almost similarly to changes in the initial parameters. However, we suggest that both the mMT and the Soloviev models can provide profiles that closely align with observations, especially since both models exhibit a change in the B_z profile sign, as can be measured by an in situ satellite at 1 AU (e.g., Regnault et al. 2020).

Moreover, by successively launching two spheromak CMEs into the heliospheric models EUHFORIA and PLUTO, respectively, Scolini et al. (2020) and Koehn et al. (2022) found that the interaction of two moderate CMEs can result in an intense geomagnetic storm. Therefore, it would be insightful to study within EUHFORIA how the thermodynamic and magnetic profiles at Earth are affected by the interaction of multiple CMEs modeled by our mMT and/or Soloviev models.

It may also be interesting to compare the performance of these two models with previously implemented models such as spheromak (Verbeke et al. 2019) and FRi3D (Maharana et al. 2022). Regarding the geometry, given that coronal mass ejections are often approximated as locally cylindrical objects with circular cross sections (Savani et al. 2011), the morphology of these two new toroidal models is likely to be closer to observations than the spheromak model, which is spherical. Moreover, the geometry of our toroidal models can be further refined using the radial stretching approach proposed by Gibson & Low (1998). This enhancement would make them even more closely resembling a realistic flux rope. Finally, we note that the computation time required for the two toroidal CMEs (on the order of a few minutes for our reference runs) is closer to that of the spheromak model than that of FRi3D. Indeed, the latter model, with standard initial parameters and a numerical configuration similar to ours, requires between 7 and 10 h for the calculation, while the spheromak requires between several tens of minutes to a couple of hours (Maharana et al. 2022).

We have found, however, that increasing the magnetic flux can cause a massive increase in computation cost because more time steps are needed to capture the time-dependent evolution on the grid. To speed up the computation, it could be worthwhile to implement these two CME models in the new inner heliospheric simulation ICARUS (Verbeke et al. 2022; Baratashvili et al. 2022). Unlike EUHFORIA, ICARUS is based on MPI-AMRVAC (Keppens et al. 2003, 2012, 2023; Xia et al. 2018), which uses advanced numerical techniques and radial grid stretching together with adaptive mesh refinement. Baratashvili et al. (2022) found that for the same spheromak CME evolving in the same solar wind, ICARUS is up to 17 times faster than EUHFORIA. Furthermore, the adaptive mesh refinement in MPI-AMRVAC could allow the capture of finer details in the simulation that may be missed when using a uniform grid. We plan to validate the implementation of these two new models using observation data in a second paper.

To conclude, the current implementation of the toroidal CMEs in EUHFORIA has laid the groundwork for additional toroidal CME modeling. By addressing all the intricacies specific to a toroidal CME structure, such as determining the insertion points, we have paved the way for a simplified adaptation process. This foundation makes it significantly easier to incorporate other toroidal CME models, such as the Tsuji solution (Tsuji 1991) or the modified Titov-Démoulin CME model (Titov & Démoulin 1999; Titov et al. 2014). These can now be implemented by simply adjusting the magnetic field within the torus.

Fig. 14 Same as Fig. 14 but the only difference between the different CMEs is the minor radius, a = 3 R⊙, a = 5 R⊙ and a = 9 R⊙. The major radius is fixed to R⊙ = 10 R⊙.

Acknowledgements

The authors thank the anonymous referee for their insightful comments. We acknowledge support from the European Union’s Horizon 2020 research and innovation program under the grant agreement no. 870405 (EUHFORIA 2.0) and from the projects C14/19/089 (C1 project Internal Funds KU Leuven), G.0B58.23N and G.0025.23N (FWO-Vlaanderen), 4000134474 (SIDC Data Exploitation, ESA Prodex-12), and Belspo project B2/191/P1/SWiM. R.K. acknowledges support by the C1 project TRACESpace funded by KU Leu-ven and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 833251 PROMINENT ERC-ADG 2018) and a FWO project G0B4521N.

Appendix A Magnetic fields in the global coordinates

Appendix A.1 Miller-Turner transformation system

According to Sect. 2.1, the magnetic field for the modified Miller-Turner CME within the torus is defined in the local (ρ_l, ϕ_l, θ_l) curved cylindrical coordinates. In order to implement in EUHFORIA, we first derive the magnetic fields in the local spherical system (e_ρ, e_θ, e_ϕ) where ρ is the radius, θ is the colatitude, and ϕ is the longitude. The two set of coordinates are related to each other by the relations:

$${\rho }_l=\sqrt{{\rho }^2+R_0^2-2\rho R_0\sin \theta }$$(A.1)

$${\theta }_l=\{\begin{array}{ll}\arcsin \frac{\rho \cos \theta }{\rho },\hfill & \text{ if }\rho \sin \theta -R_0\ge 0\hfill \\ \pi -\arcsin \frac{\rho \cos \theta }{\rho },\hfill & \text{ otherwise. }\hfill \end{array}$$(A.2)

The magnetic field in the local spherical system can be written as:

$$B_{\rho }=\sin \left({\theta }_l+\theta \right)B_{{\rho }_l}+\cos \left({\theta }_l+\theta \right)B_{{\theta }_l}$$(A.3)

$$B_{\varphi }=B_{{\varphi }_l}$$(A.4)

$$B_{\theta }=\cos \left({\theta }_l+\theta \right)B_{{\rho }_l}-\sin \left({\theta }_l+\theta \right)B_{{\theta }_l}.$$(A.5)

In a similar way, the magnetic field can also be obtained in the local Cartesian system (e_Xl, e_Yl, e_Zl) defined by the coordinates (X_l, Y_l, Z_l):

$$\begin{array}{ll}{B}_{Xl}=\hfill & \frac{B_{\rho }{X}_l}{\sqrt{X_l^2+Y_l^2+Z_l^2}}+\frac{B_{\theta }{X}_l{Z}_l}{\sqrt{X_l^2+Y_l^2+Z_l^2}\sqrt{X_l^2+Y_l^2}}\hfill \\ \hfill & -\frac{B_{\varphi }{Y}_l}{\sqrt{X_l^2+Y_l^2}}\hfill \end{array}$$(A.6)

$$\begin{array}{ll}{B}_{Yl}=\hfill & \frac{B_{\rho }{Y}_l}{\sqrt{X_l^2+Y_l^2+Z_l^2}}+\frac{B_{\theta }{Y}_l{Z}_l}{\sqrt{X_l^2+Y_l^2+Z_l^2}\sqrt{X_l^2+Y_l^2}}\hfill \\ \hfill & +\frac{B_{\varphi }{X}_l}{\sqrt{X_l^2+Y_l^2}}\hfill \end{array}$$(A.7)

$$B_{Zl}=\frac{B_{\rho }{Z}_l}{\sqrt{X_l^2+Y_l^2+Z_l^2}}-\frac{B_{\theta }\left(X_l^2+Y_l^2\right)}{\sqrt{X_l^2+Y_l^2+Z_l^2}\sqrt{X_l^2+Y_l^2}}$$(A.8)

with

$$X_l=\rho \sin \theta \cos \varphi$$(A.9)

$$Y_l=\rho \sin \theta \sin \varphi$$(A.10)

$$Z_l=\rho \cos \theta .$$(A.11)

It is worth noting that we use the magnetic field in the spherical coordinate system as an intermediate step rather than directly transforming the magnetic field from the local cylindrical coordinate system to the Cartesian coordinate system.

In order to implement the magnetic field in EUHFORIA, we need to switch from the local Cartesian system (e_Xl, e_Yl, e_Zl) to the global Cartesian system (e_X, e_Y, e_Z) of the simulation. The center of the torus in the EUHFORIA system is the point O_l at ρ = t_d, θ = θ_l, and ϕ = ϕ_l. The original coordinate system is first translated such as:

$$x_{tr}=x-t_d\sin {\theta }_l\cos {\varphi }_l$$(A.12)

$$y_{tr}=y-t_d\sin {\theta }_l\sin {\varphi }_l$$(A.13)

$$z_{tr}=z-t_d\cos {\theta }_l,$$(A.14)

where (x_tr, y_tr, z_tr) define the translated system. Then to obtain the local Cartesian system, the new coordinate system is rotated, so that the new X-axis corresponds to the direction of propagation OO_l, where O is the center of the EUHFORIA domain (cf. Fig. 5).

$$\left(\begin{array}{l}X\hfill \\ Y\hfill \\ Z\hfill \end{array}\right)=\left(\begin{array}{ccc}\sin {\theta }_l& 0& \cos {\theta }_l\\ 0& 1& 0\\ -\cos {\theta }_l& 0& \sin {\theta }_l\end{array}\right)\left(\begin{array}{ccc}\cos {\varphi }_l& \sin {\varphi }_l& 0\\ -\sin {\varphi }_l& \cos {\varphi }_l& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{l}{x}_{tr}\hfill \\ y_{tr}\hfill \\ z_{tr}\hfill \end{array}\right).$$(A.15)

Finally, to apply a tilt defined by an angle ω, the resulting system can be rotated along the new X-axis such as:

$$\left(\begin{array}{l}{X}_{rot}\hfill \\ Y_{rot}\hfill \\ Z_{rot}\hfill \end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \omega & -\sin \omega \\ 0& \sin \omega & \cos \omega \end{array}\right)\left(\begin{array}{l}X\hfill \\ Y\hfill \\ Z\hfill \end{array}\right),$$(A.16)

with (X_rot, Y_rot, Z_rot) define the final coordinate system.

Appendix A.2 Soloviev transformation system

Similarly to the modified Miller-Turner CME model, the Soloviev solution is obtained in a local coordinate system, which differs from the coordinate system used in EUHFORIA. Therefore, a coordinate transformation is necessary. The Soloviev CME model is defined in a specific cylindrical coordinate system, defined by the unit vectors $e_{R_l},e_{{\varphi }_l}$ with R_l the radius, Z_l the vertical axis, and ϕ_l the longitudinal angle. These cylindrical coordinates can be expressed in terms of the local spherical coordinates (ρ, θ, ϕ):

$$R_l=\sqrt{{(\rho \sin \theta \cos \varphi )}^2+{(\rho \sin \theta \sin \varphi )}^2}$$(A.17)

$$Z_l=\rho \cos \theta$$(A.18)

$${\varphi }_l=-\varphi .$$(A.19)

The three components of the magnetic field in the local spherical coordinate system can be obtained from the magnetic field in the cylindrical coordinate system:

$$B_{\rho }=\sin \theta B_{R_l}+\cos \theta B_{Z_l}$$(A.20)

$$B_{\theta }=\cos \theta B_{R_l}-\sin \theta B_{Z_l}$$(A.21)

$$B_{\varphi }=-B_{{\varphi }_l}.$$(A.22)

Once the components of the magnetic field are obtained in the local spherical coordinate system, they allow access to the vectors in the local Cartesian coordinate system before undergoing a translation and rotations to obtain the different vectors in the global system used in EUHFORIA. The transformation from the components B_ρ, B_θ and B_ϕ to the vector in the final EUHFORIA coordinate system is described in Sect. A.1.

References

All Tables

All Figures

	Fig. 1 Local toroidal curved cylindrical coordinates associated with the Miller-Turner CME model. In the local Cartesian system centered on Ol, the z-axis coincides with the symmetry axis of the torus.
In the text

Fig. 2 Normalized magnetic field strength in the two CME models implemented in EUHFORIA. The left panels show a sample of magnetic field lines in a torus with a minor radius a = 4 R⊙ and a major radius R0 = 10 R⊙. The right panels show the dimensionless magnetic field strength B/B0 in a poloidal cross-section. The “+” markers indicate the center of the torus. From top to bottom: magnetic field configuration in the modified Miller-Turner model, in the Soloviev model with αS = 1 and with αS = 10. When the αS parameter is too high, the magnetic field is mainly toroidal.

In the text

	Fig. 3 Particular cylindrical system associated with the Soloviev solution. The color lines are flux isocontours such as ψ = 1 (in red), ψ = 0.7 (in green) and ψ = 0.2 (in blue).
In the text

Fig. 4 Isocontours of ψ = 1 in the X − Z plane for the Soloviev solution with different geometries. The upper panel illustrates the impact of the triangularity parameter τ on the isocontour ψ = 1, while the lower panel, characterizes the influence of the elongation parameter σ for a torus with a = 1 R⊙ and R0 = 2.5 R⊙, centered at (0,0,0). The marker “+” indicates the center of the torus, while the dot marker corresponds to x = R0. In the top panel, σ is kept constant at 1.4 and in the bottom plot, τ is constrained to 0.

In the text

Fig. 5 Relation between the global Cartesian coordinate system of EUHFORIA (X, Y, Z), and the local Cartesian system in which the two CME models are defined (Xl, Yl, Zl). The distance from the center td, the colatitude θi, and the longitude ϕi define the initial position of the center of the torus, Ol. The vector OOl, where O is the center of the domain, corresponds to the Xl axis of the local Cartesian coordinate system.

In the text

Fig. 6 Theoretical profiles of magnetic field components are obtained by advancing a torus at a speed of 600 km s−1 through three virtual satellites. The first column presents an isosurface Ψ = 1 for the Soloviev solution, with three torus inclinations: ω = 0° (top panel), ω = 45° (middle panel), and ω = −45° (bottom panel). The second column represents, for these three configurations, the evolution of magnetic field components (in GSE) at the level of a virtual satellite (blue) placed directly in the direction of torus propagation. The third column corresponds to the passage through the offset virtual satellite (orange) at an angle ϕ = 10° relative to the propagation direction. The fourth column corresponds to the offset green satellite at an angle ϕ = −10°. The dashed lines correspond to the magnetic field components obtained from the Soloviev CME model, whereas the solid lines correspond to the mMT solution.

In the text

Fig. 7 Magnetic and thermodynamic profiles as function of time obtained at Earth using EUHFORIA with the two CME models. Top to bottom: speed, density, total magnetic field strength, Bx, By, Bɀ components in GSE coordinates, and plasma beta (β in log scale). The solid lines correspond to simulations where the CME is based on the mMT solution, while the dashed lines are obtained from simulations with the Soloviev solution as CME model. The green vertical lines mark the presumed sheath boundaries in the mMT simulation with Cα = 1.

In the text

Fig. 8 Visualization of the longitudinal magnetic field component in the EUHFORIA heliosphere domain. For both CME models, two initial mass densities were used: n0 = 1 × 10−16 kg m−3 (panels a and b for the mMT CME model and panels c and d for the Soloviev CME model) and n0 = 1 × 10−16 kg m−3 (panels e and f for the mMT CME model and panels g and h for the Soloviev CME model). The geometrical, magnetic, and thermodynamic parameters of the CMEs are described in Sect. 4.2.1. The left panels (a, c, e, g) show the evolution of the longitudinal field, Blon, on March 4 while the right panels (b, d, f, h) show the simulation on November 6. Each panel is divided into an equatorial cross section, a meridional cross section passing through the Sun and the Earth, and a Mercator projection at 0.1 AU from the Sun. The visualization displays the positions of the different planets in the simulation domain, the location of the current sheet, and the interplanetary magnetic field connecting the Sun and the planets, obtained from Parker spirals.

In the text

	Fig. 9 Same as Fig. 8 but for the co-latitudinal magnetic field, Bclt.
In the text

	Fig. 10 Same as Fig. 7 but the only difference between the different cases computed from the same CME model is the initial radial speed. For both CME solution, the initial speed is υr0 = 300 km s−1, υr0 = 600 km s−1 and υr0 = 900 km s−1.
In the text

	Fig. 11 Same as Fig. 7 but the only difference between the different CMEs is the initial magnetic field strength B0. For CMEs computed from the Turner solution, B0 is equal to 5 × 10−7 T, 1 × 10−6 T and 3 × 10−6 T; while it is 1 × 10−7 T, 2 × 10−7 T and 6 × 10−7 T in the Soloviev cases.
In the text

	Fig. 12 Same as Fig. 7 but the only difference between the different cases computed from the same CME model is the initial mass density, n0. For both CME solution, the initial mass density is n0 = 1 × 10−18 kg m−3, n0 = 1 × 10−17 kg m−3 and n0 = 1 × 10−16 kg m−3.
In the text

Fig. 13 Distribution of scaled mass density at the moment of arrival at Earth as a function of the initial CME density, n0, and the CME model used. The right panels show simulations with the mMT CME model, while the left panels show simulations with the Soloviev CME model. Three different initial mass densities were used, from top to bottom: n0 = 1 × 10−18 kg m−3 (panels a and b), n0 = 1 × 10−17 kg m−3 (panels c and d) and n0 = 1 × 10−16 kg m−3 (panels e and f). The different legends related to the EUHFORIA vizualisation tool are described in Fig. 8.

In the text

	Fig. 14 Same as Fig. 14 but the only difference between the different CMEs is the minor radius, a = 3 R⊙, a = 5 R⊙ and a = 9 R⊙. The major radius is fixed to R⊙ = 10 R⊙.
In the text