© The Authors 2025

1. Introduction

Extreme solar events (Gosling 1993; Schrijver et al. 2015) have been listed as one of the most significant risks coming from space, potentially leading to catastrophic events that forecasting assets must be able to predict as accurately as possible to protect humanity. Particles with high energy, called solar energetic particle (SEPs), can be accelerated during their journey by the solar wind or by magnetic structures in the interplanetary medium (Reames 2013). Flare- and shock-related processes are two of the leading candidates for the efficient acceleration of SEPs observed in situ (see, e.g. reviews by Desai & Giacalone 2016; Klein & Dalla 2017; Vlahos et al. 2019).

Magnetic storms are another type of event caused by coronal mass ejections (CMEs) that are expelled from the Sun, with speeds reaching 1000 to 2000 km/s (Schmieder et al. 2015). Fast and wide CMEs are powerful drivers of shockwaves in the corona and interplanetary space that, in turn, can efficiently accelerate particles to high energies and release them in very distant locations from the flare or eruption site (e.g. Rouillard et al. 2016; Afanasiev et al. 2018; Kouloumvakos et al. 2019, 2022). Magnetic disturbances on Earth can also be initiated by corotating interaction regions (CIRs), which form when high-speed streams (HSS) emanating from a coronal hole interact with the ambient slow solar wind (Hajra et al. 2022). However, the most severe geomagnetic storms arise when the B_z component of the CME’s magnetic field aligns in the opposite direction to the Earth’s magnetopause (Lugaz et al. 2016). This opposite alignment facilitates magnetic reconnection between the magnetic fields of these two structures, thereby channelling energy into the Earth’s inner magnetosphere (Akasofu 1981; Dungey 1961).

These events can cause catastrophic effects on human infrastructure. Geomagnetically induced currents, created due to space weather events, can enter the electricity infrastructure and saturate or damage the equipment (e.g. transformers), causing power blackouts over large areas in the case of severe storms (Eastwood et al. 2017). Additionally, high-energy particles might compromise spacecraft operations, posing a risk to equipment and astronauts’ safety (Bothmer & Daglis 2007; Baker et al. 2004). Understanding the physics of Sun-Earth transients and the mechanisms driving space weather events is crucial for reducing potential risks and associated costs.

In the heliosphere, the properties of interplanetary coronal mass ejections (ICMEs) are measured at various distances from Earth using a suite of spacecraft such as the Advanced Composition Explorer (ACE; Chiu et al. 1998), the Parker Solar Probe (Fox et al. 2016), the Solar Orbiter (Müller et al. 2020), WIND (Harten & Clark 1995), STEREO A and B (Kaiser & Adams 2007). While methods exist to reconstruct the 3D structure of ICMEs from multiple viewpoints (Rodari et al. 2018), the methods are complex due to the limited number of available observations (Démoulin 2010). Consequently, the observations are supplemented by numerical simulations designed to model the properties of the solar wind and/or predict the evolution of the Sun’s magnetic structures as they travel to Earth.

Even though some simulation tools (e.g. Manchester 2004; Mignone et al. 2007; Lugaz et al. 2005; Regnault et al. 2023; Shen et al. 2007, 2013), can be used to track a CME from the Sun to Earth, it is generally preferable to separate the study into different modelling tools due to the varying physical conditions from the solar corona to the heliosphere (Forbes et al. 2006). Additionally, having one single large domain requires numerous grid points, leading to a notable slowdown in the wall-clock times of the simulations.

Therefore, among the various tools available, we find a set of numerical simulations specifically designed to reproduce the solar corona’s complexities accurately. For instance, the Wang-Sheeley-Arge model (WSA; Wang & Sheeley 1990; McGregor et al. 2011; van der Holst et al. 2010) is commonly used to create a solar wind from the observations of photospheric magnetograms (e.g. Verbeke et al. 2019; Scolini et al. 2019; Prete et al. 2024). Other notable examples include the steady global corona model with turbulence transport and heating by Usmanov (1996), Usmanov et al. (2014, 2018), Chhiber et al. (2021); the Wind-Predict-AW code (Réville et al. 2020, 2021; Parenti et al. 2022); and the Alfvén Wave Solar atmosphere Model (AWSoM; van der Holst et al. 2010; van Driel-Gesztelyi & Green 2015; Sachdeva et al. 2019) already integrated within the Space Weather modelling Framework (SWMF; Tóth et al. 2012). A more advanced model describing the ambient plasma density, temperature, and magnetic field of the background corona is provided by the Magneto-hydrodynamics Around a Sphere (MAS) model from the solar surface to 30 R_⊙. This model was developed by Predictive Science Inc. (Riley et al. 2011; Lionello et al. 2009).

In this work, we will focus on the global 3D MHD coronal plasma solver known as COolfluid COroNal UnsTructured (COCONUT; Perri et al. 2022, 2023). COCONUT’s implicit scheme, coupled with its highly scalable parallel architecture, sets it apart from other solvers and bolsters its efficiency in swiftly and precisely addressing MHD challenges. Additionally, the Titov-Démoulin flux rope model (Titov et al. 2014) and the force-free flux rope construction using the Regularised Biot-Savart law (RBSL) method (Titov et al. 2018) in COCONUT accurately describe the initial stages of a CME evolution in the solar corona until O.1 AU (Linan et al. 2023; Guo et al. 2024a).

In addition to solar corona simulations, there is a range of 3D MHD models dedicated to studying the heliosphere from 0.1 AU to the Earth and even beyond 2 AU. Three-dimensional MHD solvers such as ENLIL (Odstrcil 2003), MS-FLUKSS (Singh et al. 2018), ICARUS (Verbeke et al. 2022; Baratashvili et al. 2022), SUSANOO-CME (Shiota & Kataoka 2016), and LFM-helio (Merkin et al. 2011; Pahud et al. 2012), which is time-dependent coupling with the coronal MAS model (Merkin et al. 2016), belong to this class of models. The present study focuses on the European heliospheric forecasting information asset (EUHFORIA; Pomoell & Poedts 2018). This space weather forecasting tool tracks the properties of one or more 3D CME models in the heliosphere after being injected into the domain at 0.1 AU.

The cone model, describing the CME as unmagnetised plasma with a self-similar expansion, was the first CME model implemented in EUHFORIA (Xie et al. 2004). The most often used model is the spheromak, representing a force-free magnetic field with a spherical geometry (Chandrasekhar & Kendall 1957; Verbeke et al. 2019). However, due to a lack of legs, spheromak falls short in accurately representing scenarios where Earth is affected by the flanks of ICMEs. To address this limitation, the ‘Flux Rope in 3D’ model (FRi3D; Isavnin 2016) was implemented by Maharana et al. (2022). This model represents a toroidal flux rope still anchored on the solar surface. Finally, Linan et al. (2024) implemented two new toroidal CMEs: One is based on the modified Miller-Turner (mMT) solution (Romashets & Vandas 2003; Vandas & Romashets 2015), while the other one is derived from the Soloviev equilibrium (Soloviev 1975). The latter two models require considerably fewer computational resources than FRi3D, drastically reducing the overall simulation time.

The accuracy of predictions made by heliospheric simulation tools such as EUHFORIA is significantly influenced by two factors: (1) the coronal model used to initiate the solar wind and (2) the characteristics of the CME model employed to simulate the actual events. For the first point, recent findings by Baratashvili et al. (2024) demonstrate that using COCONUT as the coronal model instead of WSA results in a bi-modal solar wind structure that matches in situ measurements. As for the second point, the closer the CME model aligns with the observational data, the more accurate the predictions are expected to be. However, a major challenge arises when injecting the CME model at 0.1 AU, as it omits the critical physical processes in the lower solar corona. Consequently, CMEs do not experience evolutionary changes within the solar corona before entering the model. Furthermore, tracking CMEs within the currently used WSA model is impossible, as this model relies on empirical relationships and is not a time-dependent 3D simulation.

To account for the crucial interaction between the solar wind and the CME within the solar corona, we introduce the first time-dependent coupling between the coronal model COCONUT and the heliospheric forecasting tool EUHFORIA. Specifically, we perform multiple COCONUT simulations where flux ropes evolve from the solar surface to R_b = 21.5 R_⊙. At each time step, the outer boundary is saved and used as the inner boundary in EUHFORIA. By doing this, the thermodynamic and magnetic properties of CMEs crossing the boundary are directly impacted by their propagation in the lower corona. In the future, we expect predictions made by the space weather chain COCONUT + EUHFORIA will be more accurate than those using CME models independent of the solar corona, such as the spheromak model.

This work presents the first coupling between two simulations for test cases. The application to an actual event is reserved for a future study. The paper is organised as follows: the first section is dedicated to describing COCONUT (cf. Sect. 2), with particular emphasis on the numerical scheme (cf. Sect. 2.1), the implemented flux rope models (cf. Sect. 2.2), and the handling of outputs (cf. Sect. 2.3). Next, we introduce EUHFORIA (cf. Sect. 3), highlighting the standard version of the solver in Section 3.1 as well as the numerical changes made to dynamically adjust the inner boundary based on COCONUT’s outputs (cf. Sect. 3.2). In the subsequent section (cf. Sect. 4), the coupling is tested through a set of six simulations, each varying based on the properties of the flux ropes evolving in COCONUT (cf. Sect. 4.1). After detailing the propagation of CMEs in COCONUT (cf. Sect. 4.2), we illustrate how this information is transferred to EUHFORIA (cf. Sect. 4.3) and how the thermodynamic and magnetic properties continue their evolution in the heliosphere up to Earth (cf. Sects. 4.4 and 4.5). Finally, Section 5 concludes with a presentation of the main findings of our works and discusses potential future developments.

2. COCONUT

2.1. Full MHD coronal model

The first part of our space weather forecasting chain involves tools that generate the solar corona’s magnetic field distribution and thermodynamic parameters. This modelling is achieved through the COCONUT simulation. Unlike some other contemporary solvers, COCONUT employs an implicit second-order finite-volume scheme based on the Computational Object-Oriented Libraries for Fluid Dynamics (COOLFluiD) platform (Lani et al. 2005; Kimpe et al. 2005; Lani et al. 2013). In steady-state simulations, this solver allows for using Courant-Friedrichs-Lewy (CFL) numbers significantly higher than 1 (up to 1000 or more), resulting in up to 30 times faster than state-of-the-art MHD coronal models that use explicit schemes (Perri et al. 2022).

COCONUT’s high performance also results from using an unstructured grid, which facilitates better handling of regions near the poles by avoiding the presence of singularities. Multiple spatial resolutions are possible. In this case, we used a sixth-level subdivision of the geodesic polyhedron, which results in 1.5 million prismatic cells (Brchnelova et al. 2022b). The mesh extends radially from the inner boundary at 1 R_⊙ to the outer boundary at 25 R_⊙. For coupling with EUHFORIA, only the surface at 21.5 R_⊙ is required. However, the upper boundary of COCONUT is set slightly further away to move the possible outer-boundary effects away from the coupling location (Brchnelova et al. 2022a).

A potential-field approximation is performed using the radial component of the magnetic field derived from a magnetic map to establish an initial condition for the magnetic field. Various magnetograms, such as those from the Helioseismic and Magnetic Imager (HMI) or the Global Oscillation Network Group (GONG), can serve as the inner boundary condition (Perri et al. 2023). In every instance, the magnetogram is pre-processed by conducting a spherical harmonics decomposition to smooth the magnetic field. Since it is not yet feasible to use a vector magnetogram as the inner boundary condition in COCONUT, which would provide access to the transverse components of the magnetic field, we apply Neumann (zero gradients) conditions across the inner boundary for the colatitudinal and longitudinal magnetic field components and let them evolve freely with the velocity components in such a way, that the resulting velocity is parallel to the background magnetic field Brchnelova et al. (2022a). The pressure and density are assumed to be constant everywhere (Perri et al. 2022).

After initialisation, the following conservative formulation of the MHD equations is solved:

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\left(\begin{array}{c}\rho \\ \rho \mathit{V}\\ \mathit{B}\\ E\\ \varphi \end{array}\right)+\mathbf{\nabla }·\left(\begin{array}{c}\rho \mathit{V}\\ \rho \mathit{V}\mathit{V}+\mathsf{I}(P+\frac{1}{2}{|\mathit{B}|}^2)-\mathit{B}\mathit{B}\\ \mathit{V}\mathit{B}-\mathit{B}\mathit{V}+\underline{\mathsf{I}\varphi }\\ (E+P+\frac{1}{2}{|\mathit{B}|}^2)\mathit{V}-\mathit{B}(\mathit{V}·\mathit{B})\\ V_{\mathrm{ref}}^2\mathbf{B}\end{array}\right)\\ \hfill =\left(\begin{array}{c}0\\ \rho \mathit{g}\\ 0\\ \rho \mathit{g}·\mathit{V}+\mathbf{S}\\ 0\end{array}\right),\end{array}$$(1)

with ρ the density, E the total energy, P the thermal gas pressure, $\mathit{g}(r)=-(G{M}_{\odot }/r^2){\widehat{\mathit{e}}}_r$ the gravitational acceleration,I$={\widehat{\mathit{e}}}_x\otimes {\widehat{\mathit{e}}}_x+{\widehat{\mathit{e}}}_y\otimes {\widehat{\mathit{e}}}_y+{\widehat{\mathit{e}}}_z\otimes {\widehat{\mathit{e}}}_z$ the identity dyadic and S the energy source term. The solution uses a second-order accurate finite volume discretisation and the artificial compressibility analogy (Chorin 1997; Dedner et al. 2002). The different MHD quantities are normalised by the reference values of ρ_ref = 1.67 ⋅ 10⁻¹³ kg.m⁻³, B_ref = 2.2 ⋅ 10⁻⁴ T, l_ref = 6.9551 ⋅ 10⁸ m, and the corresponding V_ref computed from the values above.

Until now, the vast majority of studies conducted with COCONUT have assumed S = 0, referring to these as polytropic runs. However, recent work by Baratashvili et al. (2024) demonstrates that by incorporating empirical source terms to better approximate the physics of the solar corona, the representation of the solar corona with COCONUT aligns more closely with observations than the purely polytropic version. In light of these findings, we have also decided to use the so-called “full MHD” formulation instead of the polytropic version in this study. The “full MHD” version of COCONUT differs from the polytropic version by adding extra terms on the right-hand side of the energy equation given in Eq. (2). When including the complex physics phenomena in the HMD equations, we switch the adiabatic index to γ = 5/3, which is a realistic polytropic index compared to the artificially reduced one used in the polytropic simulations. In the “full MHD” version, the same inner boundary conditions are used as in the polytropic COCONUT simulations. Furthermore, the MHD equation set is solved uniformly across the entire numerical domain without differentiating between open and closed magnetic field lines.

In detail, the source terms in the energy equation can be described as follows:

$$\begin{array}{c}\hfill S=-\mathrm{\nabla }·\mathbf{q}+Q_{\mathrm{rad}}+Q_H,\end{array}$$(2)

where Q_H is the coronal heating, Q_rad is the radiation loss function and −∇⋅q is the thermal conduction. According to Rosner et al. (1978), the radiative loss can be written as:

$$\begin{array}{c}\hfill Q_{\mathrm{rad}}=-n_e{n}_{p}P(T),\end{array}$$(3)

with n_e and n_p are the electron and proton number densities, and P(T) is a cooling profile depending on the temperature (Rosner et al. 1978).

The thermal conduction follows two regimes depending on whether the plasma is collisional or not (Mikić et al. 1999; Hollweg 1978). Below 10 R_⊙, the plasma is collisional, and the thermal conduction is described by the standard Spitzer-Haärm thermal conduction flux :

$$\begin{array}{c}\hfill {\mathit{q}}_1=-{\kappa }_{\Vert }(\widehat{\mathit{b}}\otimes \widehat{\mathit{b}})·\mathrm{\nabla }T,\end{array}$$(4)

with ${\kappa }_{||}=9\times {10}^{-7}{T}^{\frac{2}{5}}$ in cgs units. Beyond 10 R_⊙, in the regime which is assumed to be collisionless, the thermal conduction flux is then given by:

$$\begin{array}{c}\hfill {\mathit{q}}_2=\alpha n_{e}kT\mathit{v},\end{array}$$(5)

with α a constant defined in Hollweg (1978) and k the Boltzmann constant.

Since the solar coronal heating mechanism is still not fully understood in solar physics, the formulation of the heating term Q_H remains under investigation. Baratashvili et al. (2024), using the same magnetogram date as in this work, tested three widely-used empirical formulations for the heating term in COCONUT. First, they tried a function that depends only on the distance from the Sun’s surface. However, this uniform spherically symmetric heating did not yield a realistic bi-modal solar wind configuration, so we ruled out this option for this article. Baratashvili et al. (2024) also considered a more complex heating function approximation depending on the magnetic field strength and described as:

$$\begin{array}{cc}\hfill Q_H& =Q_H^{\exp }+Q_H^{\mathit{QS}}+Q_H^{\mathit{AR}},\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill Q_H^{\exp }& =H_0{e}^{\frac{-(r-R_{\odot })}{{\lambda }_0}},\hfill \end{array}$$(7)

where H₀ = 4.9128 * 10⁻⁷ erg cm⁻³ s⁻¹, λ₀ = 0.7 R_⊙, Q_H^QS a heating term describing the quiet Sun and Q_H^AR a term describing the effect of active regions (cf. Lionello et al. 2009, for more details). Finally, Baratashvili et al. (2024) tried a heating function directly depending on the strength of the magnetic field, as suggested by Downs et al. (2010) and Pevtsov et al. (2003) :

$$\begin{array}{c}\hfill Q_H=H_0·|\mathbf{B}|·e^{-\frac{r-R_s}{\lambda }},\end{array}$$(8)

where H₀ = 4 ⋅ 10⁻⁵ erg cm⁻³ s⁻¹ G⁻¹ and λ = 0.7 R_⊙. Using either Eq. (6) or Eq. (8) results in a realistic bi-modal solar wind in COCONUT, with a coronal magnetic field distribution that is highly similar between both formulations. To settle on the most suitable expression, Baratashvili et al. (2024) connected the COCONUT output to the heliosphere to examine how these effects propagate to Earth. Specifically, these authors used COCONUT output in the heliospheric simulation ICARUS to initialise the solar wind. This is similar to what is done with EUHFORIA in the current work, but for a single-time instance and without CME. By comparing profiles at Earth with OMNI 1-min data, Baratashvili et al. (2024) found that the modelled heliosphere in ICARUS was more consistent with the OMNI data using the heating introduced in Eq. (8) compared to the heating function in Eq. (6). However, the number density remained significantly overestimated. For this reason, we decided to use the heating term described in Eq. (8) for this study.

2.2. CME implementation

This study track the propagation of multiple CME models in EUHFORIA following their evolution in COCONUT. Tracking a CME’s propagation within COCONUT involves several steps. The solar wind in which the CME will propagate must be modelled initially. To achieve this, COCONUT is run in its relaxation mode, during which the full MHD equations are solved until a quasi-steady state is reached. The convergence is monitored by assessing the global residuals of various physical quantities according to the following formula:

$$\begin{array}{c}\hfill \mathrm{res}(a)=\log \sqrt{{\displaystyle \sum _i}{(a_i^t-a_i^{t+1})}^2},\end{array}$$(9)

where a represents the considered primitive variable, i is the spatial index, and t is the temporal index.

Once the solar wind is established, the magnetic field is modified to superimpose that of a CME model at the solar surface (cf. Linan et al. 2023, for details on the implementation of a flux rope in COCONUT). It is important to note that only the magnetic field is altered during the implementation. The density, velocity, and temperature remain unaffected by adding a CME model.

Currently, two different CME models can be implemented in COCONUT. The first model, whose implementation is described in Linan et al. (2023), is the modified Titov-Démoulin model (hereafter TDm, Titov et al. 2014). It represents a current-carrying and approximately force-free magnetic field along a toroidal geometry. The toroidal segment is partially embedded in the photosphere, with only one circular arc extending into the corona. The magnetic twist is concentrated in a thin layer near the boundary. Moreover, during the implementation, the flux rope’s plane is positioned locally perpendicular to the ambient magnetic fields, which should be potential.

Since no initial velocity is imposed, the flux rope must not be initially in equilibrium to allow for radial evolution. This is the case when the tension of the overlying magnetic field does not balance its magnetic pressure. Practically, this condition is met when the net intensity of the flux rope, I, is set higher than the Shafranov intensity as follows:

$$\begin{array}{c}\hfill I=\zeta I_S,\end{array}$$(10)

with

$$\begin{array}{c}\hfill I_S\approx \frac{4\pi R{B}_{\perp }/\mu }{ln\frac{8R}{a}-\frac{3}{2}+\frac{I_i}{2}}.\end{array}$$(11)

Here, R is the major radius of the torus, a is the minor radius, I_i represents the internal self-inductance per length of the rope, ζ is a positive integer, and B_⊥ is the ambient magnetic field perpendicular to the toroidal axis. When ζ > 1, the ambient magnetic field cannot prevent the expansion of the toroidal segment.

Unlike the TDm model, which follows a toroidal geometry, the second CME model allows for a flux rope with an arbitrary path, the RBSL model. This model was originally proposed by Titov et al. (2018) and implemented in COCONUT by Guo et al. (2024a). The magnetic field can be expressed as:

$$\begin{array}{cc}\hfill \mathit{B}& =\mathrm{\nabla }\times {\mathit{A}}_{\mathit{I}}+\mathrm{\nabla }\times {\mathit{A}}_{\mathit{F}},\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill {\mathit{A}}_{\mathit{I}}\mathbf{(}\mathit{x}\mathbf{)}& =\frac{\mu I}{4\pi }{\displaystyle {\int }_{C\cup C^{\prime }}}{K}_I(r){\mathit{R}}^{\mathbf{\prime }}(l)\frac{\mathit{dl}}{a(l)},\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill {\mathit{A}}_{\mathit{F}}\mathbf{(}\mathit{x}\mathbf{)}& =\frac{F}{4\pi }{\displaystyle {\int }_{C\cup C^{\prime }}}{K}_F(r){\mathit{R}}^{\mathbf{\prime }}(l)\times \mathit{r}\frac{\mathit{dl}}{a{(l)}^2},\hfill \end{array}$$(14)

with a the minor radius, C the axis path, C′ the mirror path of the axis relative to the photosphere, l the arc length, R the radius-vector, R′ the tangential unit vector and r the vector from the source point to the field point. According to Titov et al. (2018), K_I(r) and K_F(r) are the kernels of the RBSL model, which can be written as piecewise functions:

$$\begin{array}{c}\hfill K_I(r)=\{\begin{array}{c}\frac{2}{\pi }(\frac{\arcsin \mathrm{r}}{r}+\frac{5-2r^2}{3}\sqrt{1-r^2})0\le r\le 1\hfill \\ \frac{1}{r}r>1,\hfill \end{array}\end{array}$$(15)

and

$$\begin{array}{c}\hfill K_F(r)=\{\begin{array}{c}\frac{2}{\pi r^2}(\frac{arcsinr}{r}-\sqrt{1-r^2})+\frac{2}{\pi }\sqrt{1-r^2}+\hfill \\ \frac{5-2r^2}{2\sqrt{6}}[1-\frac{2}{\pi }arcsin(\frac{1+2r^2}{5-2r^2})]0\le r\le 1.\hfill \\ \frac{1}{r^3}r>1.\hfill \end{array}\end{array}$$(16)

The magnetic field is force-free with this particular formulation, and the net axial current follows a parabolic distribution.

As mentioned previously, the path of the RBSL model can be controlled. In this work, this path follows the equations defined in Török et al. (2010) and Xu et al. (2020):

$$\begin{array}{c}\hfill f(s)=\{\begin{array}{c}\frac{s(2x_c-s)}{x_c^2}\theta ,0\le s\le x_{\mathrm{c}}\hfill \\ \frac{(s-2x_{{}_{\mathrm{c}}}+1)(1-s)}{{(1-x_c)}^2}\theta ,x_{{}_c}<s\le 1\hfill \end{array}\end{array}$$(17)

with

$$\begin{array}{cc}\hfill x& =(s-x_c)cosf+x_c,\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill y& =(s-x_c)sinf,\hfill \end{array}$$(19)

and

$$\begin{array}{c}\hfill z(x)=\{\begin{array}{c}\frac{x(2x_{\mathrm{h}}-x)}{x_{\mathrm{h}}^2}h,0\le x\le x_{\mathrm{h}}\hfill \\ \frac{(x-2x_{{}_{\mathrm{h}}}+1)(1-x)}{{(1-x_{\mathrm{h}})}^2}h,x_{{}_{\mathrm{h}}}<x\le 1,\hfill \end{array}\end{array}$$(20)

where x_c is the intersection point of the projected curve on the COCONUT surface and the line connecting the two footpoints, θ is the angle defining the orientation of the tangent vector at the position x_c, x_h is the position of the apex, which is at the height h. The implementation of the RBSL model in COCONUT is further detailed in Guo et al. (2024a).

Once the CME is implemented using the TDm or RBSL model, COCONUT is relaunched in its time-accurate mode (Linan et al. 2023). In this mode, the backward Euler scheme used to reach a steady state is replaced by a backward differentiation formula (BDF2) with a time step fixed at 5.10⁻³ (in code units). This particular value was chosen by Linan et al. (2023) and Guo et al. (2024a) as it ensures accurate simulation without excessive computational resource usage. The linearised system was solved using the generalised Minimal RESidual (GMRES) method, complemented with a parallel Additive Schwartz preconditioner from the PETSc library toolkit¹. The convergence of the time-accurate iterative process at each time step was determined either when the residual fell below 1 × 10⁻⁴ or after ten sub-iterations. These settings were chosen after extensive testing because they allow for a fast and accurate simulation. Using more stringent values increased run-time excessively without significantly improving accuracy. Finally, intermediate simulation results (i.e. solution in terms of primitive variables) are saved every 20 iterations, corresponding to approximately 144 s in physical time.

2.3. Output pre-processing

Currently, the process of running the COCONUT and EUHFORIA simulations remains sequential: it is necessary first to run COCONUT entirely, followed by EUHFORIA, as integration into a single execution has not yet been implemented. During the execution of COCONUT, the evolution of various thermodynamic quantities must be regularly saved for later use in EUHFORIA. In this paper, we present for the first time a time-dependent coupling between COCONUT and EUHFORIA.

COCONUT offers the ability to save several types of files. Among these options, we have retrieved the output as VTU files (VTK Unstructured Grid File). These files can be read in the ParaView software (Squillacote et al. 2007), allowing us to visually verify, for example, by tracing magnetic field lines, that the simulation is proceeding correctly, which means that the CME implemented develops radially in the domain over time (e.g. Fig. 5 in Linan et al. 2023).

The coupling between the two simulations requires using the native “CFmesh” files generated by COCONUT. In particular, these ASCII files contain the mesh nodes’ coordinates, the cell-node connectivity, and the values of the thermodynamic quantities at the cell centres. Unlike the VTU files, a single CFmesh file is generated and written in parallel for the entire domain at each iteration, as the solver directly integrates the information processed by the different processors.

As is mentioned in Section 3.2, EUHFORIA requires the values of the magnetic field, velocity, density, and temperature at the surface that is located precisely at 21.5 R_⊙. Since the mesh in COCONUT is unstructured, to extract this specific surface, we proceed to interpolate the unstructured grid into a structured grid corresponding to R_b = 21.5 R_⊙ with 360 points in longitude and 180 in latitude.

The interpolation is performed upstream of EUHFORIA using the radial basis function interpolation method with a linear kernel. To determine the value at a given point, the interpolation is based on the values of the 50 cells nearest to that point on the original grid.

The interpolated surface, illustrated in Fig. 1, is saved in a new ASCII file (hereafter the coupling file or boundary map) that EUHFORIA will read. This step is repeated for all iterations saved by COCONUT, generating many text files representing the surface’s temporal evolution at 21.5 R_⊙. During the creation of the files, a 180-degree rotation in longitude is applied so that the injected CME is directed towards the Earth in EUHFORIA. It should be noted that this interpolation step is computationally costly. Future work is necessary to optimize the data transfer from COCONUT to EUHFORIA.

Fig. 1. Magnetic and thermodynamic quantities derived from the surface at R = 21.5 R⊙ in the COCONUT simulation. The panels represent the density, the temperature, the three components of the magnetic field (BrBclt and Blon), and the three components of the velocity (VrVclt and Vlon). These panels show the surface at the initial moment of the COCONUT simulation, where a CME modelled by the TDm model with ζ = 35 is implemented. These interpolated surfaces, derived from the unstructured mesh of COCONUT and saved in a coupling text file, are used throughout the relaxation phase to build up the solar wind in EUHFORIA.

We also note that the surfaces associated with longitudinal and colatitudinal speed components, cf. the last two panels in Fig. 1, exhibit numerical artefacts. These artefacts are created due to mesh non-orthogonalities causing the formation of spurious fluxes in non-dominant dynamics, as discussed in detail in Brchnelova et al. (2022b). The solar wind is radial at the front boundary before the CME arrives. However, the presence of these artefacts does not affect the results, which is presented later.

3. EUHFORIA

The following subsections describe the functionality and numerical implementation of the heliospheric model used in this work, EUHFORIA. First, we describe the standard version of EUHFORIA (cf. Sect. 3.1). Then, we discuss the updated version of EUHFORIA that enables coupling with COCONUT (cf. Sect. 3.2).

3.1. Standard version of EUHFORIA

3.1.1. Numerical scheme

COCONUT enables tracking the evolution of CMEs in the solar corona up to 21.5 R_⊙. To extend the study beyond this limit and predict the temporal evolution of various plasma quantities at different positions in the inner heliosphere, we couple COCONUT with the physics-based 3D MHD heliosphere model, EUHFORIA (Pomoell & Poedts 2018).

The interface between the two simulations was set at 21.5 R_⊙ because at this distance, we expect a supersonic and super-Alfvénic solar wind, meaning that no information travels towards the Sun. However, for a physical coupling, it is crucial to ensure that the solar wind is not sub-Alfvénic in COCONUT.

In EUHFORIA, the thermodynamic and magnetic quantities are self-consistently computed by solving the set of ideal MHD equations with gravity. The equations solved are similar to those presented for the COCONUT solver, with the critical difference that EUHFORIA is purely polytropic and does not solve the last equation in Eq. (1), which corresponds to the conservation equation for magnetic flux for hyperbolic divergence cleaning. In other words, no source term is explicitly implemented (S = 0). To model the heating that leads to the acceleration of the solar wind, the polytropic index appearing in the definition of the energy, $E=\frac{P}{(\gamma -1)}P+\frac{1}{2}\rho v^2+\frac{B^2}{2{\mu }_0}$ is taken as a non-adiabatic polytropic index (i.e < 5/3; Pneuman & Kopp 1971; Pomoell & Vainio 2012). As in Odstrcil et al. (2004), the reduced index is set to 1.5.

The MHD equations are solved in the Heliocentric Earth Equatorial coordinate system (HEEQ) using a finite-volume method combined with a constrained transport approach to ensure that the magnetic field remains divergence-free. To achieve a robust second-order accurate scheme, standard piece-wise linear reconstruction and an approximate Riemann solver are employed (Kissmann & Pomoell 2012). The Coriolis and centrifugal effects of using a non-inertial frame are neglected, assuming their contributions are negligible (Pomoell & Poedts 2018).

The mesh used in this work is uniform and extends from 0.1 AU (i.e. 21.5 R_⊙) to 2 AU with 256 cells in the radial direction and a 2° angular resolution. While the computational domain spans 360° in longitude, it covers only latitudes from −60° to −60° to avoid polar regions.

3.1.2. Boundary conditions

To establish the background solar wind in the heliosphere, the standard version of EUHFORIA reads the surface at 21.5 R_⊙ which is stored in a single boundary condition map after being generated by a coronal model. The standard EUHFORIA as described by Pomoell & Poedts (2018), only requires the 2D distribution of the radial magnetic field, B_r, radial velocity V_r, density, n, and temperature, T as inner boundary condition. In most of the studies using EUHFORIA (e.g. Maharana et al. 2023; Verbeke et al. 2019; Prete et al. 2024; Linan et al. 2024), the semi-empirical WSA model (McGregor et al. 2011; van der Holst et al. 2010) is used to establish the boundary condition map. This work aims to demonstrate how the solar corona simulation COCONUT can be used in place of the semi-empirical WSA model within EUHFORIA. However, this subsection explains how EUHFORIA currently uses the WSA boundary condition map to provide context for understanding the modifications made to the code to read the series of boundary maps generated by COCONUT (cf. Sect. 3.2).

In the WSA model, a Potential Field Source Surface (PFSS) extrapolation (Altschuler & Newkirk 1969) is performed using the magnetic field derived from a synoptic magnetogram. The magnetic field obtained is then extended up to 0.1 AU following the Schatten Current Sheet model (Schatten et al. 1969). The solar wind speed v_sw = v(f, d) is determined based on the flux tube expansion factor f and the distance d from the footpoint of the flux tube to the nearest coronal hole. Using this solar wind speed, the density and temperature at the outer boundary are calculated through empirical relations. These plasma quantities (i.e. B_r, V_r, n and T) are then saved in an inner boundary file (Pomoell & Poedts 2018). This step is performed upstream of the main EUHFORIA run. The single file thus created is used at each time step in the standard EUHFORIA as the inner boundary condition.

At the inner boundary of EUHFORIA, two layers of ghost cells are used. The values of these ghost cells are determined at each step through linear extrapolation using the boundary condition map and the data from the first in-domain cell. In addition to the information contained in the boundary file when the WSA model is used as a coronal model in EUHFORIA, we assume that the velocity is purely radial at the interface R_b = 21.5 ; R_⊙. Thus, the colatitudinal velocity component, V_clt and longitudinal velocity, V_lon, are set to zero. The same assumption is made for the colatitudinal magnetic field component. The longitudinal magnetic field is determined using the radial magnetic field and the radial velocity according to the Parker Solar Wind model (Parker 1958):

$$\begin{array}{c}\hfill B_{\mathrm{lon}}=R_b\frac{\mathrm{\Omega }{B}_r}{V_r}sin\theta ,\end{array}$$(21)

where θ is the colatitude, and Ω is the solar rotation rate, approximately equal to 2.66622373 × 10⁻⁶ rad s⁻¹.

Open boundary conditions are applied at the outer radial boundary by extrapolating the existing values beyond the domain’s edge. In contrast, a symmetric reflection of values is employed at the latitudinal boundaries to ensure coherent continuity across these boundaries.

In the standard version of EUHFORIA, which uses the semi-empirical WSA coronal model, the same boundary condition file is used throughout the simulation. To model the temporal evolution of the solar corona, a simple rotation of the boundary map is applied, assuming that the coronal model remains static during the run.

3.1.3. Implementation of CME models

In the standard version of EUHFORIA, the initial time, t = 0, corresponds to the observation time of the magnetogram used to create the boundary map via a coronal model. When t is positive, the simulation enters the forecast phase to predict the evolution of one or more CMEs in the heliosphere. However, the solar wind must be already established at t = 0 before tracking CMEs within the domain. Thus, the forecast phase is preceded by a relaxation phase to establish a stable solar wind solution that conforms to the prescribed inner boundary conditions. This phase allows the solar wind to evolve within the simulation until it reaches a steady-state configuration, ensuring the model starts from a consistent and realistic background. During this time, the simulation iterates over several steps, adjusting the plasma properties until the solution no longer changes significantly, creating a balanced initial condition for the heliospheric plasma. Typically, this process takes around 14 days of simulation time, sufficient for a solar wind moving at approximately 250 km/s to travel through 2 AU and reach the outer edge of the simulation domain (Pomoell & Poedts 2018).

During this phase, when EUHFORIA used the WSA model as coronal model the single boundary map is rigidly rotated at each time step by an angle equal to tΩ (t being negative) to maintain the domain in the HEEQ coordinate system. Once the solar wind is established and before the forecast phase, it is possible to insert into the domain all the CMEs that have occurred before the start of the forecast. Since the passage of a CME affects the heliospheric plasma for several days, this ensures that the initial conditions accurately reflect recent perturbations in the spatial environment. By integrating these CMEs, the accuracy of the forecasts is improved by accounting for significant changes in the plasma and magnetic fields caused by these events.

Presently, five different CME models are implemented in EUHFORIA, each with distinct features and limitations. The simplest is the cone model (Xie et al. 2004), which represents an unmagnetised plasma with self-similar expansion, useful for modelling the interplanetary acceleration of low-energy protons during SEP events but unsuitable for studying internal magnetic signatures of ICMEs (Wijsen et al. 2022). The spheromak model depicts a CME with a global spherical shape and a linear force-free magnetic field (Chandrasekhar & Kendall 1957; Verbeke et al. 2019), widely used but inadequate for scenarios where the flanks or legs of ICMEs impact Earth. The flux rope model (FRi3D; Maharana et al. 2022; Isavnin 2016) features an extended flux-rope geometry, including CME legs, closely matching observational data of crescent-shaped CMEs. However, it requires substantial computing time, limiting its operational use. Additionally, two toroidal CME models have been implemented to address the limitations of previous models by incorporating toroidal structures, enhancing accuracy and computational efficiency in representing CME impacts (Linan et al. 2024).

Regardless of the chosen CME model, it is inserted into EUHFORIA by making it pass through the inner boundary at R_b = 21.5 R₀. A function mask identifies the grid points where the CME intersects the inner surface. At these intersection points, the solar wind speed, magnetic field, density, and temperature are replaced by those of the CME. At each time step, the CME centre advances with the initially defined radial speed, updating the 3D mask region as the CME moves through the boundary. Injecting CME models in this manner can potentially lead to an added divergence in the magnetic field.

3.2. Updated EUHFORIA

In this subsection, we describe the code changes to EUHFORIA to enable coupling with the evolving coronal simulation, COCONUT. Unlike a static coronal model such as WSA, the thermodynamic and magnetic properties at the boundary R_b = 21.5 R_⊙ in COCONUT can evolve, reflecting the arrival and passage of a CME that is introduced at the Sun’s surface (cf. Sect. 2.2). As a result, a single file can no longer serve as the inner boundary condition for the entire EUHFORIA run, and the code has been adapted to handle a series of files, each containing the R_b = 21.5 R_⊙ surface of COCONUT at specific time steps. These files are referred to as file_t0, file_t1, and so on, where file_t0 contains the information for the 2D surface at the initial time of the magnetogram, and file_t1 contains the information for the time t₀ + 1 × dt, and so forth.

The new version of EUHFORIA, like the standard version, is still divided into three main phases: the relaxation phase, the CME insertion phase, and the forecast phase (cf. Fig. 2).

Fig. 2. Different phases of an updated EUHFORIA run. The same inner boundary file is used to construct the solar wind throughout the relaxation phase. From the magnetogram date corresponding to the first coupling file, the inner boundary ghost cells of EUHFORIA at a given time t are updated using temporal interpolation of the information in the two files bracketing this particular time t.

As previously described (cf. Sect. 3.1.2), EUHFORIA begins with a relaxation phase to establish the solar wind state across the entire domain at the initial time t = 0, based on information provided at the inner boundary and read from a boundary map. Since the first coupling file, file_t0 contains the magnetic and thermodynamic information at the interface between COCONUT and EUHFORIA at the initial time, this file is used throughout the relaxation phase to initialise the simulation.

In the standard version of EUHFORIA using the semi-empirical WSA model, it is assumed that at the inner boundary of EUHFORIA, the velocity is purely radial, the longitudinal magnetic field is proportional to the radial magnetic field (cf. Eq. (21)), and the colatitudinal magnetic field is zero (cf. Sect. 3.1.2). However, COCONUT provides all velocity and magnetic field components, which are crucial for accurately describing the propagation of a CME. As a result, in the updated version of EUHFORIA that uses COCONUT as the coronal model, no assumptions are made about the magnetic field and velocity components. This information comes directly from COCONUT. The three components of the magnetic field, the three velocity components (cf. Fig. 1), and temperature and density are stored in the coupling files. Temperature is not a primitive variable in COCONUT; it is computed using the ideal gas law for pressure and density. The EUHFORIA code was modified to read these eight variables to update the inner boundary instead of the previous four (cf. Sect. 3.1.2).

Starting from time t = 0, which corresponds to the magnetogram date and the end of the relaxation phase, EUHFORIA enters its forecast phase and updates the inner boundary at each time step using the information from the series of coupling files. However, EUHFORIA’s time step is not fixed; it varies to ensure numerical stability, as it is adjusted according to the CFL condition, which depends on the speed of disturbances within the domain–faster propagation requires smaller time steps. Since the exact time steps of the simulation are not known in advance, it is challenging to have a coupling file, file_ti, corresponding precisely to every simulation time t. Therefore, the code identifies the two files that bracket the current simulation time and performs a linear interpolation of the data in these two boundary maps to determine the required quantities. Having as many coupling files as possible, with minimal temporal gaps between them, is recommended to ensure interpolation accuracy.

The inner boundary of EUHFORIA continues to update at each time step until the simulation time exceeds the time stored in the last coupling file. At that point, and for the remainder of the forecast phase, the last available boundary map is used as the inner boundary condition, similar to how the first coupling file was used during the relaxation phase. A rigid rotation of this final coupling file is performed at each time step to model the temporal evolution of the inner boundary.

Finally, it is important to note that inserting CME models (e.g. spheromak) as done in the standard version of EUHFORIA is still possible. However, additional studies are needed to determine how these models might interact with the CMEs propagated in COCONUT and transmitted to EUHFORIA via the boundary maps.

4. CME propagation from Sun to Earth

4.1. Test cases

To demonstrate the coupling of COCONUT and EUHFORIA, we study the propagation of different CMEs from the Sun to the Earth. We focus on verifying that the information in COCONUT is accurately transmitted to EUHFORIA.

The magnetic field in COCONUT before implementing the CME is obtained via a steady-state run based on the pole-interpolated Helioseismic and Magnetic Imager (HMI) magnetogram product for July 2, 2019. To be used as input to COCONUT, the HMI magnetogram has been processed with spherical harmonics decomposition filtered with l_max = 20 (see Kuźma et al. 2023, for more details about the choice of this particular value and how the projection works.). The pre-processing smooths and homogenizes the magnetic field. However, significant negative and positive polarities exist between longitudes 0 and 60 (cf. Fig. 3).

Fig. 3. Synoptic magnetogram for July 2, 2019. The upper panel corresponds to the original filled HMI magnetogram, while the bottom panel shows the processed magnetogram used as input in COCONUT.

This particular date was chosen because Baratashvili et al. (2024) demonstrate, by comparing with white light images, that the full MHD version of COCONUT effectively reproduces the solar coronal magnetic field distribution. Additionally, 2019 corresponds to a minimum of solar activity, which helps minimize the effects of the background solar wind on CME propagation. Moreover, the simplicity of the magnetic configuration during the solar minimum enhances numerical convergence. Baratashvili et al. (2024) also showed that the output from COCONUT for this particular date could be used to create the solar wind for a heliospheric simulation, ICARUS. In the future, the coupling should also be tested in a maximum solar activity. However, during such intense activity, we expect significant interactions between the solar wind and flux ropes (Linan et al. 2023), adding unnecessary complexity for an initial introduction to the COCONUT-EUHFORIA coupling, as we present here.

Figure 4 shows the radial magnetic field in the meridional plane obtained after the steady-state run for the July 2, 2019, magnetogram. The distribution is bi-modal because the speed is faster at the poles than near the equator. This is due to more intense heating at the poles, where the magnetic field is generally stronger (cf. Baratashvili et al. 2024, for more details about the impact of heating on the solar wind for this particular date).

Fig. 4. Meridional plane of the radial velocity after the steady-state run of the full MHD COCONUT simulation. The black lines represent a sample of the magnetic field lines of the solar wind.

In the background magnetic field, we implement the two different CME models available: the Titov-Démoulin flux rope model and the RBSL model (cf. upper panels in Fig. 5). We aim to demonstrate that the coupling works regardless of the CME model. In detail, we conducted six simulations, three with each of the two CME models (cf. Table 1). In all cases, the CME is positioned at a longitude of −180° and a latitude of 0° to ensure evolution mainly in the equatorial plane. This choice of longitude was made to avoid any possible disturbances caused by a high-speed stream located on the opposite side of the Sun at that time (cf. Sect. 4.2).

Fig. 5. Visualisation of the TDm and RBSL models implemented in COCONUT. The upper panels show the CME models just after their insertion in the solar wind, while the other panels show their expansion 52.8 min after the simulation starts in the top and side views. The coloured lines represent a sample of magnetic field lines from the flux ropes. The origin for tracing these field lines is a sphere with a radius of R = 0.1 R⊙, positioned at the positive polarity. The grey field lines represent the magnetic field from the background solar wind. The radial magnetic field at the solar surface (measured in Gauss) is also displayed.

As in Linan et al. (2023), the flux ropes modelled using the TDm formalism have a major radius of R = 0.3 R_⊙, a minor radius of a = 0.1 R_⊙, and their centre is offset by d = 0.15 R_⊙ from the solar surface. Implementing the model results in the addition of two polarities with an area equal to 4839 Mm² on the solar surface (cf. upper right panel in Fig. 5). The three CMEs implemented using the TDm flux rope model differ only by the parameter ζ used, which takes values of 35, 50, and 70. Increasing the ζ parameter increases the ratio between the magnetic field of the flux rope and the ambient magnetic field. Thus, as described in Section 2.2, setting ζ > 1 ensures that the tension of the overlying magnetic field does not prevent the expansion of the flux rope. The higher the ζ parameter, the higher the magnetic flux. According to Linan et al. (2023) and Regnault et al. (2023), the magnetic flux is directly related to the plane-of-sky speed value measured in the simulation’s initial moments. This relationship is confirmed in the following section (cf. Sect. 4.2).

For the three implemented RBSL flux ropes, we set x_c = x_h = 0.5. The apex of the flux rope is located at a height of h = 0.17 R_⊙. As with the TDm CME models, the minor radius is a = 0.1 R_⊙. The writhe angle θ is 60°. This results in flux ropes with an S-shaped morphology (cf. upper left panel in Fig. 5). In the solar corona, filaments with positive helicity and this particular morphology are frequently observed (Cheng et al. 2017).

The parameter ζ does not appear in the implementation of the RBSL model. Stability is controlled by explicitly modifying the initial magnetic flux. In our work, the three CMEs modelled with the RBSL model have fluxes of 8 × 10²⁰ Mx, 12 × 10²⁰ Mx and 18 × 10²⁰ Mx, respectively. Our tests revealed that such fluxes, in this particular solar wind configuration, allow for the expansion of the flux ropes.

4.2. Initial stage of CME propagation in COCONUT

We initiate the time-dependent COCONUT run once the CMEs are implemented in the steady-state wind. Since the computational resources required for a time-dependent run of COCONUT are significant, we stop the simulation once the magnetic cloud has finished crossing the boundary and there are no significant variations in the different quantities within the domain. Consequently, the simulations stop after 12 000 iterations, corresponding to a physical time span of approximately 24.12 h as the CME has fully exited the domain by this point. The 2D surface at R_b = 21.5 R_⊙ is saved every 20 iterations, resulting in the generation of 600 coupling files.

The two lower panels of Figure 5 show the magnetic field topology after about t = 24.12 minutes in the COCONUT simulations with the CMEs named ‘TDm_3’ and ‘RBSL_3’ (cf. Table 1). The other cases are not represented as they do not show differences in the evolution of the magnetic structure. In Figure 5, we observe that the evolution is very similar to that described in Linan et al. (2023) and Guo et al. (2024a), who studied the propagation of the two flux rope models in detail in the polytropic version of COCONUT. The self-consistent radial expansion of the flux rope creates a pinching of the legs, as described in the 2D “standard model” for CME originally developed by Carmichael (1964), Sturrock (1966), Hirayama (1974), Kopp & Pneuman (1976) and later extended to 3D by Aulanier et al. (2012, 2013), Janvier et al. (2013). The current sheet layer formed by the convergence of upward and downward field lines is a region conducive to magnetic reconnection. However, in COCONUT, magnetic reconnection is purely due to numerical dissipation since there is no prescribed resistivity in the set of MHD equations which are solved. Finally, reconnection leads to the formation of post-flare loops, as seen in observations (Schmieder et al. 1995).

During the propagation, the TDM flux rope remains in the equatorial plane, while the “S-shape” of the RBSL is preserved and found higher in the solar corona (cf. bottom left panel in Fig. 5). This result aligns with the study by Guo et al. (2024b), which focuses on the impact of flux rope morphology on propagation in COCONUT. This indicates that the magnetic distribution of an ICME is directly linked to its progenitor in solar source regions.

Among the features shared with those described in the works using the polytropic version of COCONUT (Linan et al. 2023; Guo et al. 2024a,b), we observe an accumulation of matter ahead of the flux rope that intensifies as the CME expands (cf. Fig. 6). This region of heated and compressed solar matter results from the flux rope travelling without allowing enough time for the solar wind to flow around it. This phenomenon is observed by in situ spacecraft when the CME is sufficiently faster than the surrounding solar wind (Regnault et al. 2020).

Fig. 6. Cross-sections along the equatorial plane of the density for the COCONUT simulation with the TDm flux rope named ‘TDm_3’. The coloured lines represent a sample of the magnetic field lines of the flux rope. The yellow line indicates the Rb = 21.5 R⊙ boundary. A sheath characterised by high density develops ahead of the flux rope during propagation.

Figure 7 presents the equatorial distribution of the heating term, Q_H, the temperature, and the radial velocity at various times for the simulation with the CME named ‘TDm_3’. The other cases have common patterns. Shortly after the simulation begins, at the moment t = 2.4 min (i.e. the first saved time step), the tension of the surrounding magnetic field is insufficient to prevent the expansion of the flux rope. This results in a high radial velocity at the flux rope location (cf. the top right panel in Fig. 7). The initial radial velocities for different simulations can be estimated using this plane. These values are summarised in Table 1. For simulations involving a TDm flux rope, the velocity ranges from 1138 km/s for the CME with the weakest magnetic flux (TDm_1) to 1780 km/s for that with the strongest flux (TDm_3). Similarly, for simulations with an RBSL flux rope, the radial velocity varies from 1263 km/s (TDm_2) to 1804 km/s (TDm_3). As noted by Linan et al. (2023) and Guo et al. (2024a), the greater the magnetic flux, the higher the initial velocity. However, for approximately the same flux, 8 × 10²⁰ Mx, the ‘TDm_2’ case exhibits a lower initial velocity than the ‘RBSL_2’ case. Since the two CMEs do not share the same geometry, the initial resistance created by the surrounding magnetic field is not exactly the same (Guo et al. 2024b), which could explain this difference.

Fig. 7. Equatorial plane of the distribution of the heating term, temperature, and radial velocity at different times for the CME named ‘TDm_3’. The left panels show the heating QH in erg cm−3 s−1 G−1 as defined by Equation (8). The middle panels display the temperature in Kelvin, while the right panels present the radial velocity in km s−1. Each panel includes an overall view of the domain and a zoomed-in view near the Sun’s surface. From top to bottom, the times represented are t = 2.4 min, t = 4 h and t = 24 h.

In the processed magnetogram used as input for COCONUT, we observe strong active regions between O° and 60° longitude (cf. Fig. 7). These regions, characterised by strong magnetic fields, lead to significant heating, as the heating term Q_H scales with magnetic field strength, as described in Section 2.1. Consequently, in the second quadrant (i.e. the region with X < 0 and Y > 0) of Figure 7, we observe high values of Q_H. This heating causes particle acceleration, driving high-speed flows, as reported by Baratashvili et al. (2024), and is visible in the last column of the first row in Figure 7.

Additionally, at the initial stage, the magnetic field of the flux rope we implemented on the solar surface is stronger than the surrounding magnetic field. In the first row of Figure 7, we can logically observe that the heating term, and therefore the temperature, is elevated in the region near the core of the CME.

Later in the simulation, at t = 4 hours, the CME still resides within the domain. Several temperature zones can be identified (cf. middle panel in Fig. 7). To the right, near X = −20 R_⊙, a crescent-shaped area of high temperature is visible, which can be linked to the sheath where temperature rises due to the compression of matter ahead of the flux rope’s progression. Closer to the Sun, around X = −15 R_⊙, temperatures decrease near the flux rope’s core.

During the implementation of the flux rope in COCONUT, we added two opposite polarities on the surface of the Sun. As a result, at t = 4, near the initial insertion region, temperatures remain elevated due to coronal heating at the CME footpoints.

It is also noted that the heating associated with Q_H remains localised near the surface, as it decreases exponentially with increasing radial distance r (cf. Eq. (8)).

Regarding the velocity distribution, two high-speed regions are observed (cf. middle right panel in Fig. 7). The first is located at the bulk of the CME, which continues to expand, while the second is localised in the wake of the CME. This rapid plasma flow is associated with magnetic reconnection occurring where opposite magnetic field lines forming the legs of the CME converge (Ando et al. 2015).

By the last moment of the simulation at t = 24 hours, it is assumed (cf. Sect. 4.4) that the CME has already traversed the domain. Despite this, high-speed streams created during the CME’s wake are observed (cf. bottom right panel in Fig. 7). The simulation was continued for an additional 24 hours with no noticeable changes in the velocity distribution. This is due to a major limitation of the current COCONUT version, where the solar surface’s magnetic field remains fixed throughout the simulation. Consequently, the footpoints of the CME remain on the solar surface even though the flux rope is no longer anchored to these polarities (cf. Linan et al. 2023, for more details). Enhanced magnetic fields in the active region initially implemented lead to significant solar coronal heating. As a result, due to plasma acceleration, the final distribution features two high-speed streams: one in the second quadrant, caused by the active regions initially present on the solar surface as explained before, and the other in the X-direction, generated by the constant heating of the solar corona at the initial CME footpoints, which does not dissipate during the simulation.

4.3. Thermodynamic and magnetic quantities in the coupling files

As detailed in Sections 2.3 and 3.2, the three components of the magnetic field, of the velocity as well as the temperature, and the density at the surface R_b = 21.5 R_⊙ are regularly saved to files. This series of files is used as the inner boundary for EUHFORIA.

Figure 8 shows the magnetic field amplitude computed using the magnetic field components, which are saved in the coupling file at the moment when the flux rope, ‘TDm_3’, begins crossing the R_b = 21.5 R_⊙ surface of COCONUT. This figure shows a weak magnetic field intensity line corresponding to the equatorial streamer belt. We also observe that the increase in the magnetic field associated with the arrival of the flux rope is mainly located at a colatitude of 90° and a longitude of 0°. Since a 180-degree rotation was applied during the creation of the file (cf. Sect. 2.3), we can deduce that the propagation direction remains aligned with the active region of implementation. No major deflection of the CME was observed during its propagation in COCONUT. This could result from our choice to use a solar minimum case for the background solar wind.

Fig. 8. Two-dimensional surface map of the magnetic field amplitude as saved in a coupling file. The magnetic field amplitude is computed using the three magnetic field components Br, Blon, Bclt saved in the 138th ASCII file used for the coupling between COCONUT and EUHFORIA for the case named ‘TDm_3’. The 138th file, marking the time t = 5.5 hours, captures when the flux rope starts crossing the surface at Rb = 21.5 R⊙.

To understand the data transferred to EUHFORIA, we read all 600 coupling files for all six simulations. We extract the values at a colatitude of 90° and a longitude of 0° and plot the temporal evolution of various global quantities in Fig. 9: the magnetic field amplitude B, the density n, the velocity amplitude v, and the temperature T. It is important to note that in EUHFORIA, as described in Section 3.1.2, the inner boundary rotates at each time step by an angle equal to tΩ. This means that the longitude of 0° in the boundary map is perfectly aligned with Earth only at the initial time t = 0 (the magnetogram date). Since the solar rotation rate is slow, Fig. 9 still provides a good estimation of the perturbation transferred towards Earth.

Fig. 9. Velocity, magnetic field, density, and temperature evolution as a function of time, starting from the magnetogram date. Each line corresponds to a different simulation in the background solar wind. These profiles are consistent with the propagation of a flux rope with a sheath ahead. The purple vertical lines approximately mark the start of the sheath, the beginning of the magnetic ejecta (ME), and the end of the ME for the simulation named “TDM_3”.

Figure 9 shows that the different simulations have quite similar profiles. The main differences between all the cases are the amplitude of the peaks reached and the duration of the disturbances generated. In detail, the arrival of the disturbance caused by the propagation of the flux ropes increases all quantities. The first to reach a local maximum is the temperature. At the peak of the temperature, we notice a slight change in the slope of the total magnetic field, |B|. According to Linan et al. (2023), this could indicate the transition between the sheath and the magnetic ejecta. During the crossing of the magnetic ejecta, the temperature decreases and then increases in the wake of the CME. After the flux ropes have passed, the temperature is higher than the initial moment due to constant heating, as explained in the previous section (cf. Sect. 4.2).

The magnetic field profile shows a single increase before returning to its initial value once the magnetic ejecta has finished crossing the boundary R_b = 21.5 R_⊙. For a given CME model (RBSL or TDm), the higher the initial magnetic flux, the higher the maximum value reached since the flux is directly proportional to the strength of the magnetic field within the flux rope.

Similarly, the velocity peaks in Fig. 9 are proportional to the initial magnetic flux. The higher the magnetic flux, the higher the velocity profile values reached. It should be noted, however, that for each simulation, the maximum reached is never more than 50% of the initial speed listed in Table 1. This speed reduction can be attributed to the drag force exerted by the ambient solar wind, which acts to restrain the expanding flux rope. Such a force accounts for the gradual deceleration of CMEs as they travel through space (Sachdeva et al. 2015, 2017). As described in Section 4.2, the two local maxima observed in the velocity profile can be attributed to different phenomena. The first peak corresponds to the passage of the bulk of the CME, while the second is due to the acceleration created by the reconnection outflow in the wake of the CME (Shibata & Tanuma 2001).

Regarding the density in Figure 9, the amplitude of the peaks in the density profiles is directly related to the speed of the CMEs: faster CMEs generate higher peak densities. This is because faster CMEs compress the plasma in front of them more effectively, leading to a denser, more compact region (Winslow et al. 2015). Indeed, the fastest CMEs, ‘TDm_3’ and ‘RBSL_3’, show the highest peak. Conversely, the slower CMEs, ‘TDm_1’ and ‘RBSL_1’, exhibit lower peak densities due to less effective compression. After the passage of the CME, the density is lower than before. This can be attributed to the fact that the CME has swept up and removed a significant amount of solar wind plasma from its path, creating a rarefied region behind it.

Finally, the patterns observed in Figure 9 are consistent with the observations made in the COCONUT domain described in Section 4.2. Therefore, we can assume that the interpolation of the unstructured grid data of COCONUT, which was necessary to write the boundary file, had little impact on the transferred physics.

4.4. Propagation continuation in EUHFORIA

Once the COCONUT simulations are complete and the coupling files generated, these files are used in EUHFORIA as boundary conditions described in Section 3.2. For our EUHFORIA run, eleven days are used for the relaxation phase and seven days for the forecast phase. As a result, there are enough coupling files to cover only 1/7 of the forecast phase. During the remaining six days, the last available file is used.

Figure 10 shows the equatorial distribution of temperature, radial velocity, and density in COCONUT and EUHFORIA at various times. Due to differing time steps in COCONUT and EUHFORIA, there is no exact temporal correspondence between the outputs of COCONUT and those of EUHFORIA. Consequently, the difference in time between the Paraview file produced by EUHFORIA and those made by COCONUT is noted in minutes in the top right corner of each panel. The first black circle marks the boundary of the COCONUT domain.

Fig. 10. Composite visualisation of the equatorial plane from COCONUT and EUHFORIA. Each panel contains two parts. The inner disc is delineated by a black circle derived from COCONUT’s output, while the rest is obtained from EUHFORIA. The region covered by COCONUT is also shown in the bottom left corner of each panel. The left panels represent the temperature distribution in Kelvin, the middle panels represent the radial velocity distribution in km/s, and the right panels represent the density multiplied by the square of the radius. From top to bottom, various time snapshots from EUHFORIA are presented. The number below the date in the top right of each panel indicates the time difference in minutes between the EUHFORIA and COCONUT timestamps. Each panel includes markers at the positions of different planets, virtual satellites, and the Parker spirals connecting them to the Sun.

The top three panels of Figure 10 display the equatorial plane after the relaxation phase in EUHFORIA for the case named ‘TDm_3’. Throughout the 11-day relaxation period, the same coupling file is repeatedly read in EUHFORIA to create the solar wind. In Figure 10, it is apparent that for the three variables (n, v_r, T), the EUHFORIA domain directly extends from the COCONUT domain. The solar wind maintains a coherent structure across both simulations, and the transition between COCONUT and EUHFORIA is seamless, showing no significant changes in the amplitude of these three variables. The information in the coupling file appears to be fully integrated into EUHFORIA. Notably, a high-speed stream opposite Earth is clearly distinguishable in the simulation. As previously discussed (cf. Sect. 4.2), this is due to intense heating in the solar corona.

We can conclude that the steady-state solar wind in EUHFORIA can be constructed from the COCONUT interface at R_b = 21.5 R_⊙.

After the relaxation phase, and for about one day in simulation time, the inner boundary of EUHFORIA is updated through temporal interpolation between the two nearest coupling files to the time t being processed. The second row in Figure 10 captures the approximate moment when the sheath begins to enter EUHFORIA, occurring four and a half hours after the relaxation ends. Since the disturbance caused by the flux rope propagation has just reached the R_b = 21.5 R_⊙ surface, the solar wind in EUHFORIA has only slightly evolved during these four hours. It is particularly noteworthy that despite reading different boundary files and performing temporal interpolations up to this point, no major and apparent numerical artefacts significant enough to impact the large-scale structure of the solar wind have been generated within the domain. For the other quantities saved in the coupling files and not represented in Figure 10 (such as the magnetic field components), a smooth transition at the interface between COCONUT and EUHFORIA can also be observed.

As a further demonstration that information is effectively transferred from COCONUT to EUHFORIA during the reading of boundary files, Figure 11 shows the magnetic field lines of the ‘RBSL_3’ and ‘TDm_3’ flux ropes shortly after they cross the interface at R_b = 21.5 R_⊙. The depicted field lines correspond to the nose of the CMEs, which is corroborated by the distribution of the B_z component of the magnetic field that exhibits a sign change, as expected in the Titov-Démoulin model (Linan et al. 2023). As in COCONUT (cf. Sect. 4.2), the TDm flux rope is mainly in the equatorial plane, while the RBSL presents an “S-shape” as initially implemented. We conclude that the morphology of the CME is well transferred from one simulation to the other.

Fig. 11. Visualisation of the TDm flux rope model and the RBSL model in EUHFORIA. The upper panel shows the CME model for the case named the case named ‘RBSL_3’, while the bottom panel shows the simulation with the case named ‘TDm_3’. The coloured lines represent a sample of magnetic field lines from the flux ropes. The 3D sphere represents the inner boundary of EUHFORIA at Rb = 21.5 R⊙. For the TDm simulation, the equatorial plane shows the Bz magnetic field component in nT.

In EUHFORIA, the coupling files are read and used as boundary conditions as long as the simulation time is less than the date of the last available files, corresponding to about one day in simulation time. Afterwards, the previous available boundary file is utilised until the simulation concludes. The third row in Figure 10 shows the equatorial plane approximately 6 hours after the date corresponding to the last coupling file. By this time, the disturbance has propagated toward Earth with radial expansion. Between the virtual satellites labelled ‘SC01’ and ‘SC02’, a high temperature and density region is observed. Based on COCONUT observations, we identify this region as the sheath, which continues to develop in EUHFORIA. As the flux rope moves forward, material that cannot escape accumulates upstream of the propagation. Like in COCONUT, a higher speed peak is not located in this region but in the wake of the disturbance.

Although the disturbance is primarily visible in the direction of propagation (the X-direction), a circular band of high density covering the entire domain is discernible. We suggest that this is due to a shock created upstream of the flux rope in all directions from the onset of the COCONUT simulation when the CME begins its evolution at speeds exceeding 1700 km/s.

Finally, at the very end of the simulation, as shown in the last row of Figure 10, the solar wind has drastically changed from the initial conditions. Nevertheless, the high-speed stream present after the relaxation phase has shifted according to the co-rotated frame. Another high-speed stream with greater amplitude is also noticeable in the positive X-direction. As explained in Section 4.2, the velocity in the wake of the CME never decreases due to the continuous heating associated with maintaining the active region that originated the flux rope. This rapid plasma flow is constantly transferred into EUHFORIA, resulting in a high-speed region. It is also important to note that the density in this area is lower than initially because the flux rope has pushed material ahead, creating a cavity in its wake (cf. Sect. 4.2).

4.5. Profiles at Earth

EUHFORIA is used mainly to predict the evolution of a CME’s thermodynamic and magnetic properties at Earth. Consequently, Figure 12 shows the evolution of the total velocity, total magnetic field, density, and temperature at Earth for the different simulations performed in this study.

Fig. 12. Evolution of Earth’s magnetic and thermodynamic quantities in EUHFORIA. From top to bottom, the figure shows the profiles of the total velocity in km/s, the total magnetic field in nT, the density in m−3, and the temperature at Earth as predicted by EUHFORIA. Each coloured line represents a different simulation, as in Fig. 9 and the purple vertical lines approximately mark the start of the sheath, the beginning of the ME, and the end of the ME for the simulation named ’TDm_3’.

First, in Figure 12, we observe that the amplitude of the different profiles depends on the simulation, as noted at the interface between the two simulations (Sect. 4.3 and Fig. 9). For example, the initially fastest flux ropes, specifically the cases TDm_3 and RBSL_3, have the highest peak velocities at Earth. Additionally, we notice that the peak values are of the same order of magnitude as those observed at the interface (cf. Fig. 9), indicating that the flux ropes undergo significant deceleration mainly in the solar corona, modelled here by COCONUT. Regarding the velocity profile, the maximum is preceded by a sharp increase in each simulation. In COCONUT, the velocity increase was smoother. However, this may be due to the spatial resolution used. For Figure 12, the values are obtained every 10 minutes, while for Figure 9, the data points are approximately 2.4 minutes apart.

Similarly to the velocity, the total magnetic field profiles’ amplitudes follow the same hierarchy described at the interface between COCONUT and EUHFORIA (cf. Sect. 4.3). The higher the initial magnetic flux of the flux ropes, the higher the amplitude at Earth.

The density and temperature profiles exhibit complex patterns indicating the passage of multiple structures at Earth. Despite a noticeable increase due to the sheath passage, the maximum density ranges between 6 × 10¹⁰ m⁻³ and 4 × 10⁸ m⁻³, which is nearly two orders of magnitude smaller than the maximum values obtained in COCONUT at R_b = 21.5 R_⊙, which vary between 1.6 × 10¹⁰ m⁻³ and 2.2 × 10¹⁰ m⁻³ (cf. Sect. 4.3). The dispersion of matter ahead of a CME as the magnetic structure extends is observed both in the solar corona and the heliosphere. For instance, using a sample of 45 ICMEs observed by Helios 1/2 and the Parker Solar Probe, Temmer & Bothmer (2022) found that the region in front of the magnetic driver becomes denser than the ambient solar wind at about 0.06 AU, with density decreasing as a linear function of distance. A drop in temperature also accompanies this reduction in density.

After the passage of magnetic and thermodynamic perturbations, around July 6 at 10 AM, the temperature is higher than before the event, and the density is lower, consistent with observations made in COCONUT. This consistency again demonstrates a coherence between the changes in the solar wind in the solar corona modelled by COCONUT and those occurring in the heliosphere modelled by EUHFORIA. However, more detailed analyses, which are beyond the scope of this study, are needed to fully explain patterns specific to the evolution in EUHFORIA, such as the sawtooth structure of the temperature and the significant increase in density after the passage of the magnetic structure.

So far, we have described only the evolution of the total magnetic field, |B|. However, the detailed distribution of the magnetic field is crucial since the strongest geomagnetic storms occur when the B_z component of the magnetic field of the CME is opposite in direction to that of the magnetopause (Lugaz et al. 2016). We placed virtual satellites in COCONUT and EUHFORIA at various distances along the propagation direction to track the evolution of the flux rope distributions. Specifically, we examine the time evolution of the magnetic field components, B_x, B_y, and B_z at distances of R = 10 R_⊙ in COCONUT, at the interface between the two simulations (by reading the coupling files), at R = 68 R_⊙ in EUHFORIA, and at Earth. We aim to capture the arrival of the CME towards Earth by positioning the virtual observation points at a colatitude of θ = 90° and a longitude of lon = 0°. It’s worth noting that COCONUT and EUHFORIA separately do not evolve in the same reference frame. To transition to the co-rotated frame of EUHFORIA, the 21.5 R_⊙ boundary of COCONUT is rotated at each time step (cf. Sect. 3.1.2). This means that the virtual satellites in COCONUT and at the interface are slightly offset from those placed in EUHFORIA. Therefore, we do not expect to have perfectly identical profiles between COCONUT and EUHFORIA.

To enhance the figure’s readability, we present only the flux ropes with the highest magnetic flux, the cases named ‘TDm_3’ and ‘RBSL_3’. The other instances exhibit similar profiles but with lower amplitudes.

The first column of Figure 13 shows the evolution of the magnetic field components at ten solar radii. This point is within the domain covered by COCONUT. First, we see that the B_x component is not zero, indicating that the measurement point is slightly offset from the nose of the CME. The profiles of the three magnetic field components depend on the flux rope model used. For example, both models have a B_y profile that exhibits symmetry with respect to the horizontal axis. The B_z component follows a similar evolution in both simulations, with two sign changes. However, the TDm flux rope reaches negative values more than twice as small as the RBSL flux rope, whose negative bump is more discreet. In the simulation with the TDm flux rope model, the first positive bump in B_z, accompanied by a negative peak in the B_x and B_y components, corresponds to the magnetic field of the sheath according to Linan et al. (2023).

At the interface between the two simulations (see the second column in Fig. 13), the profiles for the three magnetic field components are very similar to those obtained closer to the Sun. The main sign changes are present for the B_y and B_z components. According to Linan et al. (2023), Guo et al. (2024a,b), in a minimum activity solar wind in COCONUT, the flux ropes develop radially from the Sun to the outer boundary without deterioration of their overall structure. This explains why we expect to obtain similar profiles at different distances from the Sun in COCONUT.

Fig. 13. Evolution of the magnetic components at different distances from the Sun. The first row depicts the evolution of the Bx component in nT, the second row shows the By component, and the third row presents the Bz component. The columns correspond to different distances from the Sun at a colatitude of θ = 90° and a longitude of lon = 0. Specifically, the first column shows the evolution at R = 10 R⊙ in COCONUT, the second column at Rb = 21.5 R⊙, which is the interface between COCONUT and EUHFORIA, the third column at R = 68 R⊙ in EUHFORIA, and the final column illustrates the evolution at Earth.

In Fig. 13, we can see that while the profile shapes are preserved, the amplitude of the curves has significantly decreased as we move away from the Sun’s surface. For example, the B_z maximum decreases from approximately 2226 nT at 10 R_⊙ to around 729 nT at 21.5 R_⊙, representing a reduction of about 67.25%. The reduction in the magnetic field of a CME due to its radial expansion and interaction with the surrounding solar wind is a well-known phenomenon, thanks to in situ observations and measurements made at different distances in the heliosphere by multiple spacecraft (e.g. Liu et al. 2005; Winslow et al. 2015).

In EUHFORIA, at approximately 68 solar radii, the main patterns of the different magnetic field components are still clearly present, which is consistent with efficient coupling. Nevertheless, for the RBSL₃ simulation, the profile of the B_x component shows some differences compared to the previous profile, such as a noticeable sign change occurring around 25 hours after the start of the simulation. A positive bump is also observed at the rear of the magnetic ejecta in the B_y component, which was not seen at the interface and in COCONUT. Assuming that the differences are not due to poor coupling (which would have been evident in the previous section), profile differences could potentially indicate a slight deflection of the CME during its evolution in EUHFORIA (Kay et al. 2013, 2015). These differences are also a result of the frame discrepancy between COCONUT and EUHFORIA, which prevents us from having perfectly aligned virtual satellites. Moreover, it should be noted that we did not expect to find the same profiles as at the interface because, after travelling more than 40 solar radii, the flux rope has developed and interacted with the solar wind.

At Earth (cf. last column in Fig. 13), the amplitude of the different profiles has drastically decreased compared to those obtained in COCONUT. We can perform a power law fit by taking the maximum magnetic field values at the four distances shown in Fig. 13. We find that the decrease in the magnetic field follows a linear trend with a slope exponent of approximately −1.50 for both the RBSL and TDm cases. This slope is higher than the value reported by Winslow et al. (2015), who derived an exponent of −1.89 ± 0.14 based on observational data, but it is closer to the one obtained by Liu et al. (2005) (−1.40 ± 0.08). However, it is important to note that our fit is highly approximate since we consider four data points. This limitation means the fit’s precision and robustness are constrained, and any conclusions drawn should be interpreted with caution.

Once again, the profiles in EUHFORIA at Earth are consistent with those obtained closer to the Sun. The characteristic pattern in the evolution of the B_y component, notably, is observed at all distances from the Sun. Indeed, the orange curve (RBSL_3) shows a decrease followed by a pronounced positive peak before returning to its initial value. This significant variation is consistently evident across different distances, illustrating the coherent propagation of the magnetic field structure throughout the corona and the heliosphere.

Despite the common patterns observed at previous distances, the profile at Earth shows notable changes compared to those obtained at a distance of 68 solar radii (R = 68 R_⊙). Being far from the inner boundary in EUHFORIA, these differences are not due to coupling issues but rather to the dynamics of traversing the heliosphere. The most significant change is the drastic decrease in the amplitude of the first positive peak of the B_y component in the case named TDm_3, dropping from around 157 nT to 12 nT, with the local minimum shifting from −110 nT to −35 nT. As mentioned earlier, this could be due to a deflection of the CME in the heliosphere, causing Earth to traverse the flank of the CME rather than its nose. Another contributing factor could be the erosion of the CME’s magnetic structure due to interactions with the ambient solar wind.

This erosion results from magnetic reconnection between the CME’s internal magnetic field and the interplanetary magnetic field (IMF) (Hosteaux et al. 2021). However, as the CME travels through the heliosphere, it encounters a solar wind whose polarity and magnetic field components vary significantly. This variability means that magnetic reconnection is not uniform throughout the CME structure. Consequently, reconnection rates and the extent of erosion are influenced by the spatially and temporally fluctuating conditions along the CME’s path, causing certain areas of the CME to undergo more erosion than others.

As a result, the erosion process is complex, varying both spatially across the CME and temporally as it progresses through different solar wind conditions. The observed changes in the magnetic field components at different distances in the heliosphere could also potentially be due to a rotation of the CME structures (e.g. Regnault et al. 2023; Ma et al. 2024). We attempted to highlight this rotation by tracing the magnetic field lines associated with the CME at various stages of the simulation, similar to the approach presented in Fig. 11. However, we could not clearly demonstrate a distinct rotation within the magnetic structure in either the TDm or RBSL simulations. Currently, limitations in spatial resolution make it difficult to conclusively identify rotation in these structures, and future simulations with higher resolution may improve our ability to capture this phenomenon more clearly.

Finally, both rotation and erosion are complex processes that evolve as the CME encounters different solar wind conditions. Further analysis is needed to fully understand the interactions between the solar wind and the CME, including the CME’s deflection, rotation, and spatially non-uniform erosion, all of which impact its evolution through the heliosphere in EUHFORIA. Such an in-depth analysis, however, lies beyond the scope of this paper.

5. Summary and conclusion

In the heliospheric simulations with EUHFORIA, CME models are inserted at the inner boundary at R_b = 21.5 R_⊙. This insertion is done entirely independently of the coronal model used to construct the solar wind. As a result, we miss all the CME dynamics and evolution in the solar corona (cf. Sect. 3.1). To address this, we present the first time-dependent coupling between the coronal simulation with COCONUT and the 3D space weather forecasting simulation with EUHFORIA.

Using a specific magnetogram as input, we performed a set of six runs in COCONUT, where a flux rope is implemented at the solar surface. In three cases, the flux rope was modelled with the TDm CME model, while in the other three, we used the RBSL model (cf. Sect. 4.1). For this study, we used the new full MHD formalism developed by Baratashvili et al. (2024) and detailed in Section 2.1. At regular intervals, the three components of the magnetic field, the three components of the velocity, the temperature, and the density on the 2D surface at R_b = 21.5 R_⊙ are saved in a set of coupling files.

We modified the EUHFORIA code so that the values of the thermodynamic variables and magnetic components at the inner boundary R_b = 21.5 R_⊙ are updated at each time step by reading the created boundary files (cf. Sect. 3.2). After briefly presenting the propagation of the flux ropes in COCONUT and the impact of the full MHD formalism on the properties of the solar wind (cf. Sect. 4.2), we focused on how the simulated coronal data are transmitted into EUHFORIA (cf. Sect. 4.4). We particularly examined the time profiles at Earth generated by the arrival of CMEs initially implemented in COCONUT (cf. Sect. 4.5). Among the various results, our main findings are :

• The full MHD formalism, used for the first time for time-dependent COCONUT runs in a minimum activity solar wind background, does not disrupt the propagation of the initially implemented flux ropes. As in the polytropic version of COCONUT (Linan et al. 2023; Guo et al. 2024a), we observe a self-consistent radial expansion of the flux rope with matter accumulating upstream and the pinching of the CME legs below the core. However, future studies are needed to validate the time-dependent full MHD formalism in a maximum activity scenario.

• Since the inner boundary of COCONUT remains unchanged throughout the simulation, the footpoints of the CME cause continuous heating of the solar corona. This results in a high-speed stream in the wake of the CME. In future studies, we suggest modifying COCONUT to allow for the temporal evolution of the solar surface, for instance, using a series of magnetograms. The goal is to gradually dissipate the active region created by the initial addition of the flux rope.

• In EUHFORIA, the same coupling file from COCONUT is used to initially create the solar wind after a relaxation process. We find that the solar wind thus developed in EUHFORIA is a direct extension of the solar wind created in COCONUT. The transition between the two simulations is smooth, indicating that the creation and reading of the boundary file do not alter the information transmitted by COCONUT.

• The initial perturbations caused by the propagation of a flux rope in COCONUT, visible as significant increases in speed, temperature, and density, continue to develop in EUHFORIA. This validates the time-dependent coupling between the two simulations. Moreover, in the standard version of EUHFORIA, a possible sheath would only develop from 0.1 AU onwards since the coronal model used was not affected by the propagation of a CME. In various studies using EUHFORIA, the size and amplitude of the sheath are often underestimated compared to observational expectations. However, with the coupling to COCONUT, a pre-formed sheath is inserted at the inner boundary. We suggest that this should better represent the sheath in the heliosphere, meaning closer to observations.

• The thermodynamic and magnetic profiles observed at Earth are consistent with those obtained at the interface between the two simulations, reinforcing the idea that simulation data are well transmitted from COCONUT to EUHFORIA. While the peak speed observed at Earth is close to that obtained at the interface of the two simulations, we find a significant decrease in the amplitude of temperature, density, and magnetic field. This result is consistent with various observations using multi-spacecraft in situ measurements, indicating that the density of the region ahead of the flux rope decreases linearly with distance in the heliosphere (e.g. Liu et al. 2005; Winslow et al. 2015).

• By examining the evolution of the three magnetic field components B_x, B_y, B_z at different distances from the Sun, we found that, despite new dynamics due to interaction with the solar wind in the heliosphere, the main structural characteristics of the magnetic field are preserved as the CME propagates outward. This means that the properties of the flux rope, such as its morphology and magnetic flux, need to be carefully controlled at the Sun’s surface in COCONUT since they directly impact the profiles observed at Earth and, consequently, on the geoeffectiveness of an event. The initial conditions set in the coronal model play a crucial role in determining the behaviour of the CME throughout its journey to Earth.

As a primary outcome of the reported implementation and validation efforts, we can conclude that the resulting new simulation chain COCONUT + EUHFORIA is a valuable tool for understanding the dynamics of CMEs from the Sun to the Earth. Now that the coupling has been introduced, the next logical step is to test the effectiveness of this simulation chain in predicting profiles generated by an actual event. For this kind of study, the initial properties of the flux ropes, such as their geometry and magnetic flux, should be derived from observations rather than arbitrarily set, as we have done in our work.

Several improvements can be made to ensure that the COCONUT + EUHFORIA toolchain is fully operational for space weather purposes. In particular, both simulations need to be able to run concurrently so that the outer boundary of COCONUT can be directly processed in EUHFORIA without waiting for the run to finish, similar to the workflows of the interactive tools CORHEL-CME (Linker et al. 2024). Additionally, the parameters should be studied to identify settings (resolution, time step, etc.) that optimize computational time while maintaining a sufficiently good accuracy. In this study, we deliberately chose an excellent temporal resolution in COCONUT (dt = 0.005 in code units) and a high spatial resolution in EUHFORIA. With this setup and the solver configuration used for the simulations, each COCONUT run lasted approximately one day and 19 hours. In contrast, an EUHFORIA run required 5 hours of wall time, plus several additional hours to create the coupling files from the text files generated by COCONUT. Finally, improving the Visualisation tools available for displaying results would also be beneficial, and for this, we can rely on the figures presented in this study as a foundation.

Meanwhile, COCONUT will continue to be developed to enable an increasingly accurate representation of the solar corona that is closer to observations. We are implementing new heating profiles among the current development axes based on the density gradient perpendicular to the background magnetic field. We also intend to consider the heating associated with the wave turbulence-driven heating mechanism.

Acknowledgments

These results were obtained in the framework of the projects C16/24/010 C1 project Internal Funds KU Leuven), G0B5823N and G002523N (WEAVE) (FWO-Vlaanderen), 4000145223 (SIDC Data Exploitation (SIDEX2), ESA Prodex), and Belspo project B2/191/P1/SWiM and the AFOSR basic research initiative project FA9550-18-1-0093. For the computations, we used the VSC–Flemish Supercomputer Center infrastructure, which is funded by the Hercules Foundation and the Flemish Government department EWI. LL acknowledges support from the ERC-2022-Advanced Grant 101095310 (TerraVirtualE, PI: Giovanni Lapenta). SP is funded by the European Union. Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or [name of the granting authority]. Neither the European Union nor the granting authority can be held responsible for them. This project (Open SESAME) has received funding under the Horizon Europe programme (ERC-AdG agreement No 101141362).

References

All Tables

All Figures

Fig. 1. Magnetic and thermodynamic quantities derived from the surface at R = 21.5 R⊙ in the COCONUT simulation. The panels represent the density, the temperature, the three components of the magnetic field (BrBclt and Blon), and the three components of the velocity (VrVclt and Vlon). These panels show the surface at the initial moment of the COCONUT simulation, where a CME modelled by the TDm model with ζ = 35 is implemented. These interpolated surfaces, derived from the unstructured mesh of COCONUT and saved in a coupling text file, are used throughout the relaxation phase to build up the solar wind in EUHFORIA.

In the text

	Fig. 2. Different phases of an updated EUHFORIA run. The same inner boundary file is used to construct the solar wind throughout the relaxation phase. From the magnetogram date corresponding to the first coupling file, the inner boundary ghost cells of EUHFORIA at a given time t are updated using temporal interpolation of the information in the two files bracketing this particular time t.
In the text

	Fig. 3. Synoptic magnetogram for July 2, 2019. The upper panel corresponds to the original filled HMI magnetogram, while the bottom panel shows the processed magnetogram used as input in COCONUT.
In the text

	Fig. 4. Meridional plane of the radial velocity after the steady-state run of the full MHD COCONUT simulation. The black lines represent a sample of the magnetic field lines of the solar wind.
In the text

Fig. 5. Visualisation of the TDm and RBSL models implemented in COCONUT. The upper panels show the CME models just after their insertion in the solar wind, while the other panels show their expansion 52.8 min after the simulation starts in the top and side views. The coloured lines represent a sample of magnetic field lines from the flux ropes. The origin for tracing these field lines is a sphere with a radius of R = 0.1 R⊙, positioned at the positive polarity. The grey field lines represent the magnetic field from the background solar wind. The radial magnetic field at the solar surface (measured in Gauss) is also displayed.

In the text

	Fig. 6. Cross-sections along the equatorial plane of the density for the COCONUT simulation with the TDm flux rope named ‘TDm_3’. The coloured lines represent a sample of the magnetic field lines of the flux rope. The yellow line indicates the Rb = 21.5 R⊙ boundary. A sheath characterised by high density develops ahead of the flux rope during propagation.
In the text

Fig. 7. Equatorial plane of the distribution of the heating term, temperature, and radial velocity at different times for the CME named ‘TDm_3’. The left panels show the heating QH in erg cm−3 s−1 G−1 as defined by Equation (8). The middle panels display the temperature in Kelvin, while the right panels present the radial velocity in km s−1. Each panel includes an overall view of the domain and a zoomed-in view near the Sun’s surface. From top to bottom, the times represented are t = 2.4 min, t = 4 h and t = 24 h.

In the text

Fig. 8. Two-dimensional surface map of the magnetic field amplitude as saved in a coupling file. The magnetic field amplitude is computed using the three magnetic field components Br, Blon, Bclt saved in the 138th ASCII file used for the coupling between COCONUT and EUHFORIA for the case named ‘TDm_3’. The 138th file, marking the time t = 5.5 hours, captures when the flux rope starts crossing the surface at Rb = 21.5 R⊙.

In the text

Fig. 9. Velocity, magnetic field, density, and temperature evolution as a function of time, starting from the magnetogram date. Each line corresponds to a different simulation in the background solar wind. These profiles are consistent with the propagation of a flux rope with a sheath ahead. The purple vertical lines approximately mark the start of the sheath, the beginning of the magnetic ejecta (ME), and the end of the ME for the simulation named “TDM_3”.

In the text

Fig. 10. Composite visualisation of the equatorial plane from COCONUT and EUHFORIA. Each panel contains two parts. The inner disc is delineated by a black circle derived from COCONUT’s output, while the rest is obtained from EUHFORIA. The region covered by COCONUT is also shown in the bottom left corner of each panel. The left panels represent the temperature distribution in Kelvin, the middle panels represent the radial velocity distribution in km/s, and the right panels represent the density multiplied by the square of the radius. From top to bottom, various time snapshots from EUHFORIA are presented. The number below the date in the top right of each panel indicates the time difference in minutes between the EUHFORIA and COCONUT timestamps. Each panel includes markers at the positions of different planets, virtual satellites, and the Parker spirals connecting them to the Sun.

In the text

Fig. 11. Visualisation of the TDm flux rope model and the RBSL model in EUHFORIA. The upper panel shows the CME model for the case named the case named ‘RBSL_3’, while the bottom panel shows the simulation with the case named ‘TDm_3’. The coloured lines represent a sample of magnetic field lines from the flux ropes. The 3D sphere represents the inner boundary of EUHFORIA at Rb = 21.5 R⊙. For the TDm simulation, the equatorial plane shows the Bz magnetic field component in nT.

In the text

Fig. 12. Evolution of Earth’s magnetic and thermodynamic quantities in EUHFORIA. From top to bottom, the figure shows the profiles of the total velocity in km/s, the total magnetic field in nT, the density in m−3, and the temperature at Earth as predicted by EUHFORIA. Each coloured line represents a different simulation, as in Fig. 9 and the purple vertical lines approximately mark the start of the sheath, the beginning of the ME, and the end of the ME for the simulation named ’TDm_3’.

In the text

Fig. 13. Evolution of the magnetic components at different distances from the Sun. The first row depicts the evolution of the Bx component in nT, the second row shows the By component, and the third row presents the Bz component. The columns correspond to different distances from the Sun at a colatitude of θ = 90° and a longitude of lon = 0. Specifically, the first column shows the evolution at R = 10 R⊙ in COCONUT, the second column at Rb = 21.5 R⊙, which is the interface between COCONUT and EUHFORIA, the third column at R = 68 R⊙ in EUHFORIA, and the final column illustrates the evolution at Earth.

In the text