© The Authors 2023

1. Introduction

Solar flares are sudden bursts of magnetic energy that are released from localized regions of the solar atmosphere called active regions. These emissions of radiation can be accompanied by the expulsion of magnetic plasma from the corona into the interplanetary medium. This is known as a coronal mass ejection (CME) and can be observed with white-light coronagraphs. Solar activity follows an 11-yr cycle with a minimum of activity where only a few CMEs are observed and a maximum of activity characterized by increased magnetic activity levels (Hathaway 2015).

The solar activity can have various impacts on the Earth’s magnetic environment and thus on human technologies (Bothmer & Daglis 2007), and many industrial sectors (e.g., energy, telecommunications, transportation) could be negatively affected (Hapgood & Thomson 2010). In addition to the damage that bursts of energetic particles can cause to space- and ground-based technologies, air crews may also be impacted by solar radiation (Griffiths & Powell 2012). To minimize risks and potential costs, it is necessary to have a thorough understanding of the physics of the Sun-Earth transients and the processes that cause space weather events.

To improve the protection of both on board crews and electronic devices, it is necessary to be able to predict the eruptivity of active regions sufficiently in advance (Schrijver et al. 2015). There are various forecasting methods with varying levels of effectiveness (Barnes et al. 2016), mainly using machine learning or deep learning artificial intelligence algorithms (e.g., Park et al. 2018, 2020). These methods often involve the parameterization of observed solar data in order to characterize the targeted active region (e.g., Bobra & Couvidat 2015; Bobra & Ilonidis 2016; Falconer et al. 2008, 2011; Li et al. 2020). From vector solar magnetograms, various magnetic field parameters are calculated and used as eruptive signatures (Guennou et al. 2017). The quantities at the heart of these algorithms rely upon magnetic field properties (e.g., the free magnetic energy Leka & Barnes 2003), current properties (e.g., the current density Zhang 2001; Török et al. 2014), or the magnetic field polarity inversion line properties (Falconer et al. 2008). However, no algorithm is yet able to predict all solar activity, partly due to the rarity of major X-class events (Aschwanden & Freeland 2012). In recent years, a new quantity defined as the ratio of the current-carrying helicity to the relative helicity (Pariat et al. 2017) has been found to be a good indicator for the eruptivity of numerical magnetic configurations in 3D parametric simulations (Linan et al. 2018; Zuccarello et al. 2018) and for observed active regions (James et al. 2018; Moraitis et al. 2019; Price et al. 2019; Thalmann et al. 2019, 2021; Gupta et al. 2021; Lumme et al. 2022). However, this criterion is not yet used in automatic forecasting models due to the complexity of computing it in the solar atmosphere (Linan et al. 2018, 2020).

Even if we were able to predict an eruptive flare before it occurs, we would also need to be able to track its propagation and evolution throughout the heliosphere because not all CMEs impact Earth and not all CMEs have the same properties nor the same geoeffectiveness (cf. Regnault et al. 2020, for a compilation of CME profiles from 20 yr of in situ observations). Specifically, when the magnetic field lines of a CME and the Earth’s magnetosphere have opposite orientations, magnetic reconnection can occur, resulting in geomagnetic storms with a potentially significant impact on the Earth’s environment (Gonzalez et al. 1994).

In the interplanetary medium, CME signatures are probed at different distances from the Earth by many spacecraft, such as ACE (Chiu et al. 1998), the Parker Solar Probe (Fox et al. 2016), STEREO A and B (Kaiser & Adams 2007), the Solar Orbiter (Müller et al. 2020) and WIND (Harten & Clark 1995). However, there are not enough satellites to perfectly determine the structure of a CME throughout its propagation (Démoulin 2010). One possibility to efficiently reconstruct the complex 3D structure of a CME from a limited set of 1D observations is to use multiple viewpoint reconstruction techniques (e.g., Rodari et al. 2018). However, these techniques rely on approximations of the overall structure of the CME. It is worth noting that the future space mission PUNCH (Polarimeter to UNify the Corona and Heliosphere) will provide four new spacecrafts equipped with narrow and wide field imagers and that the mission is aimed at producing continuous 3D images of the solar corona (DeForest et al. 2022). This should lead to a better understanding and determination of deformation effects that CMEs can undergo as they propagate, including expansion, aging, and erosion (Démoulin et al. 2008).

Observations can be supplemented with numerical simulations by various 3D magnetohydrodynamic (MHD) solvers, such as ENLIL (Odstrcil 2003), EUHFORIA (Pomoell & Poedts 2018), MS-FLUKSS (Singh et al. 2018), and SUSANOO-CME (Shiota & Kataoka 2016), which are all capable of tracking CME propagation in a realistic description of the background solar wind. The latter is usually based on the empirical Wang-Sheeley-Arge (WSA) coronal model, where its properties are derived from the solar wind speed and the expansion of coronal magnetic field lines (Wang & Sheeley 1990). In these simulations, CMEs are inserted at 0.1 AU (21.5 R_⊙) by modifying the solar wind properties (e.g., speed, density, magnetic field components). The parameters of the CMEs are derived from 3D models. Among them, we can mention the cone model describing the CME as an unmagnetized plasma with a self-similar expansion (Xie et al. 2004), the toroidal-like models where the CME is a flux rope connected to the Sun (e.g., “Flux Rope in 3D”, or FRi3D; Isavnin 2016; Maharana et al. 2022), and the linear force-free spheromak model representing a CME with a global spherical shape (Chandrasekhar & Kendall 1957; Verbeke et al. 2019). The initial geometric and magnetic parameters of these models can be partially derived from remote sensing observations. In particular, the magnetic flux and helicity can be deduced from magnetograms, while the geometric and kinematic parameters are obtained using the graduated cylindrical shell model based on white-light images (Scolini et al. 2019; Maharana et al. 2022). Finally, the accuracy of the forecast provided by a simulation depends on the accuracy of the background solar wind, the chosen CME models, and their mutual interaction.

The main limitation of inserting the CME only at 0.1 AU is that it ignores all the physics that occur lower in the solar corona. The CMEs do not evolve in the solar corona before its insertion, and they are impossible to track in the currently used corona model (WSA) since it is not a time-dependent 3D simulation but only a set of empirical relations. Hence, the inserted CME does not interact with the solar wind before the inner boundary. However, this interaction is crucial for accurately determining the geometric and magnetic characteristics of a CME (cf. Green et al. 2018, for a description of the first stages of an eruption). For example, Asvestari et al. (2022), using the spheromak model, found that the CME tilt depends on ambient magnetic field strength and orientation. The interaction with the solar wind can also lead to a deflection of the CME into a streamer (Zuccarello et al. 2011).

To address this limitation, this paper introduces the propagation of magnetic flux ropes in a recently implemented 3D MHD coronal plasma solver called COolfluid COroNal UnsTructured (COCONUT; Perri et al. 2022). The computational MHD model in COCONUT was originally developed to replace the empirical model currently used by the EUropean Heliospheric FORecasting Information Asset (EUHFORIA; Poedts et al. 2020). The current version of COCONUT uses polytropic heating, which means that the thermodynamics of the model are simplified, while the magnetic configuration is realistic. In the polytropic model, the ratio of specific heats, γ, are assumed to be set to 1.05, and the coronal heating, radiation, and conduction terms are neglected (however, their inclusion is the focus of ongoing efforts). The main objective of this work is to demonstrate the potential use of COCONUT to track the propagation of flux ropes in the corona and to provide at 21.5 R_⊙ the necessary information for heliospheric simulations. The boundary conditions imposed on EUHFORIA (or other heliospheric models) should be more accurate than those provided by an independent CME model (e.g., spheromak) and therefore should produce more realistic forecasts. To achieve the aforementioned objective, it was decided not to model a real CME event, which would have required matching the properties of the flux rope with its active region of origin. Instead, we followed the propagation of several “unobserved” flux ropes, which enables us to change the solar wind configuration as well as the properties of the flux rope to determine the advantages and limitations of COCONUT.

Recently, Regnault et al. (2023) studied the propagation of a flux rope from the low corona up to 1 AU using the PLUTO multiphysics code (Mignone et al. 2007). However, in Regnault et al. (2023), the magnetic field configuration of the solar wind is dipolar, while in this study the propagation is examined in two realistic magnetic configurations driven by magnetograms corresponding to a maximum and a minimum of solar activity. Additionally, COCONUT and PLUTO are two different coronal models with different grids, boundary conditions, and numerical schemes. The COCONUT numerical parameters have been optimized to provide the best comparison with observations (Brchnelova et al. 2022a,b). Another difference is that COCONUT uses an implicit solver, while PLUTO uses an explicit solver, which means that the COCONUT code can run within operational times (Perri et al. 2022), making it more suitable for actual forecasting.

The development of coronal MHD models has progressed significantly over the years, with various approaches and techniques being employed to model the solar corona. The first models were based on idealized plasma conditions and did not incorporate observational data into their boundary conditions (e.g., Pneuman & Kopp 1971). However, data-driven models have since been developed, such as the 3D MHD model of Linker et al. (1999), which was validated through white-light images of eclipses and interplanetary observations (see also Mikic et al. 1996).

Among other existing 3D MHD models based on systematic comparisons to both in situ and remote sensing data for validation, some models allow for the study of physical processes that are not yet incorporated in COCONUT. Examples include the steady global corona model with turbulence transport and heating of 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), which allows for a realistic description of the temperature and density distribution thanks to the treatment of wave dissipation, heat conduction, and radiative cooling. This model has already been implemented within the Space Weather Modeling Framework (SWMF; Tóth et al. 2012). The Magnetohydrodynamics Around a Sphere (MAS) solver developed by Predictive Science, Inc. is also one of the most advanced models and is based on time-dependent resistive MHD, including comprehensive energy transport mechanisms from radiation and thermal conduction to Alfvén wave heating (Mikić et al. 1999, 2018).

Although the different models possess apparent advantages, it is worth noting that some of them demand substantial numerical resources to generate a realistic 3D depiction of the solar corona. As a result, the COCONUT solver’s implicit time scheme, combined with its highly scalable parallel architecture, distinguishes it from other solvers and enhances its effectiveness in solving MHD problems quickly and accurately. This justifies our selection of the COCONUT solver for this work.

The structure of our paper is as follows. We begin by introducing the COCONUT code that is used to create two realistic solar wind configurations based on HMI magnetograms from a quiet and an active Sun (cf. Sect. 2). Then, we present the flux rope model used in this study (cf. Sect. 3.1). Special attention is paid to how the flux rope is implemented in two solar wind backgrounds (cf. Sect. 3.2) and to the initial parameters that differentiate the simulations. In Sect. 4, we describe the early stages of propagation and the limitations imposed by high-speed streams, while in Sect. 5 we focus on density, magnetic field, and velocity profiles at 21.5 solar radii. In particular, we discuss the impact of the solar wind and the magnetic flux of the flux ropes on these profiles. Finally, in Sect. 6, we present a comprehensive summary of our work and put forward conclusions for further development.

2. The numerical framework

2.1. COCONUT

In this study, we use the COCONUT code, which is based upon the Computational Object-Oriented Libraries for Fluid Dynamics (COOLFluiD) platform. The COOLFluiD platform was originally developed for multiphysics and, particularly, for the simulations of both flow and plasma (Lani et al. 2005, 2013; Kimpe et al. 2005). The COCONUT solver has been extensively described and validated in Perri et al. (2022) and Kuźma et al. (2023). The second-order finite volume scheme that COCONUT uses is implicit, which means that for steady-state simulations, Courant-Friedrichs-Lewy (CFL) numbers much higher than one can be used, enhancing the performance of the code considerably when compared to other state-of-the-art MHD coronal models.

In addition, the grid is unstructured in COCONUT, which means that more advanced grid refinement techniques (e.g., r-adaptative algorithms Ben Ameur & Lani 2021) and flux reconstruction (high-order Flux Reconstruction algorithms Vandenhoeck & Lani 2019) can be used in the future to further enhance the performance. The mesh used in the code is a sixth-level subdivision of the geodesic polyhedron with 1.5 million cells. More details about the grid and its effects in COCONUT are presented in Brchnelova et al. (2022a).

The code uses a second-order accurate finite volume discretization along with the artificial compressibility analogy (Chorin 1997), which is similar to the hyperbolic divergence cleaning (HDC) of Dedner et al. (2002), to ensure that the divergence of the B-field remains close to zero. The following conservative formulation of the MHD equations was solved by COCONUT (Perri et al. 2022):

(1)

where, ρ is the density, P is the thermal gas pressure, is the gravitational acceleration, and corresponds to the identity dyadic. The last equation in the set above corresponds to the conservation equation for magnetic flux for hyperbolic divergence cleaning. This way, we assured that the divergence of the magnetic field in the domain is negligible. The entire MHD formulation was normalized 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_A computed from the values above.

To create the surrounding field (i.e., without flux rope), the state variables were derived from a one-point backward differentiation scheme. In contrast to the steady-state solver setup introduced in Perri et al. (2022), this study required a transformation into a time-accurate model in order to track the CME propagation. This was achieved by using a three-point backward time discretization with a time step limited to 1 × 10⁻² (in code units; i.e., 14.4 s in physical time). The resulting linearized system was solved using the Generalized Minimal RESidual (GMRES) method with a parallel Additive Schwartz preconditioner from the PETSc library¹. The time-accurate iterative process was considered converged at each time step if the residual is lower than 1 × 10⁻⁴, or after four iterations. We conducted numerous tests to determine the optimal values for the parameters of time step, maximum number of iterations, and targeted residual in order to achieve the best accuracy while minimizing numerical resources. In particular, we conducted tests using time steps that were ten times larger (i.e., 10⁻¹). Unfortunately, this led to convergence issues, resulting in significant alterations in the simulation results during the propagation. We also attempted a simulation using a time step that was ten times smaller (i.e., 10⁻³). There were no notable deviations in the results from what is observed in the simulations that employ a time step of 10⁻² (cf. Sect. 4). The only differences observed were small (< 5%) discrepancies in the magnitude of the profiles presented in Sect. 5 between the 10⁻² and 10³ time steps. However, the simulation that used a time step of 10⁻³ lasted more than three days. Hence, we concluded that the accuracy gains were insufficient to justify the allocation of significantly more numerical resources. It is worth noting that the limitations we encountered, such as the existence of high-speed streams (cf. Sect. 4.2), persisted even when using a time step of 10⁻³.

With the final set of parameters, the simulations on the GENIUS cluster of the Vlaams Supercomputer Centrum² required between 20 and 26 h of computation using 196 cores in parallel to reach 6000 iterations (i.e., the prescribed stop condition). After 6000 iterations, which equates to approximately 24 h in physical time, the magnetic ejecta had successfully traversed the boundary at 0.1 AU in most of the simulations (as discussed in Sect. 5). The only exception is the slowest simulation with ζ = 5 in the minimum activity, as the magnetic ejecta was still in the process of crossing the boundary when the simulation ended. However, this simulation still provided us with valuable information regarding the profiles of various magnetic and thermodynamic quantities, their maximum values, and the CME’s geometry during its propagation (cf. Sects. 4 and 5). It is also noteworthy that shorter simulation durations were sufficient for faster rope flux (e.g., 3750 iterations instead of 6000 in the simulations with ζ > 3 in the maximum activity solar wind). However, we chose to maintain uniformity in our results by employing the same simulation duration for all cases.

The boundary conditions are prescribed as follows: The inner boundary is at 1 R_⊙. The density on the inner boundary is set to a constant value of ρ_ref = 1.67 × 10⁻¹³ kg m⁻³. According to Parker’s solar wind model, a small outflow from the surface is prescribed to follow the magnetic field lines (see Brchnelova et al. 2022b, for details), and the radial component is determined based on the imposed magnetogram. The magnetograms used were obtained by the Helioseismic and Magnetic Imager (HMI) on board the Solar Dynamic Observatory (SDO) satellite (Scherrer 2011), following the study of Perri et al. (2023). To couple the model with EUHFORIA, the outer boundary should be set to at least at 21.5 R_⊙. However, according to the recommendations of Brchnelova et al. (2022b), the grid extends beyond 21.5 R_⊙ to 25 R_⊙ in order to reduce the impacts of the outer boundary on the quantities that will be used in the heliospheric simulations (e.g., EUHFORIA). Figure 17 in Brchnelova et al. (2022b) provides a graphical depiction of the final grid used.

2.2. Solar wind configurations

The first step in our work was to model the solar wind in which the flux rope would propagate. This was achieved by running COCONUT in its relaxation mode, as presented in Perri et al. (2023). As an initial condition, we computed a potential field approximation using a fast finite volume solver implemented within COCONUT for the Laplace equation (cf. Sect. 2.3 in Perri et al. 2023). Then, the global magnetic field distribution was obtained after the MHD relaxation.

Two magnetic configurations were considered for our further discussion. The first test case, based on magnetograms from the solar eclipse of 2 July 2019, corresponds to a minimum of solar activity. The second test case represents a maximum of solar activity and is from the solar eclipse of 20 March 2015. Cases corresponding to total solar eclipses were originally chosen in order to optimize the comparison with observations (cf. Perri et al. 2022, 2023; Kuźma et al. 2023).

The two left panels of Fig. 1 show the original pictures of the solar eclipse as obtained from 128 images showing a solar minimum (top left; Druckmüller and Aniol 2019) and 29 images showing a solar maximum (bottom left; Druckmüller et al. 2015). The right panels illustrate numerical results that were obtained using COCONUT. The sphere color map represents the radial component of the magnetic field at the inner boundary, while white lines represent magnetic field lines of the magnetic coronal structures in the plane of sight.

Fig. 1. Comparison of the solar eclipse images (left) and numerical results (right) for the case of minimum (top) and maximum (bottom) of solar activity. The minimum activity image corresponds to the solar eclipse of 2 July 2019 (© 2019 Miloslav Druckmüller, Peter Aniol: http://www.zam.fme.vutbr.cz/~druck/Eclipse/Ecl2019ch/Tres_Cruses/TC_347mm/0-info.htm), and the maximum activity image corresponds to the solar eclipse of 20 March 2015 (© Miloslav Druckmüller, Shadia Habbal, Peter Aniol, Pavel Starha: http://www.zam.fme.vutbr.cz/~druck/Eclipse/Ecl2015s/800mm/00-info.htm). The red-blue sphere color map represents the radial component of the photospheric magnetic field expressed in Gauss, and the white lines represent magnetic field lines of the magnetic coronal structures in the plane of sight. The reconstructed main streamers align closely with the position and shape of coronal structures observed in white light.

Once the solar wind was modeled and before inserting the flux rope in the global coronal model, we needed to validate the numerical scheme by comparing its output to real features and structures of the solar corona. It was revealed by Yeates et al. (2018), who compared seven different models applied to the same simulation of the coronal magnetic field, that even for the same input dataset, different models can produce considerably different results. In light of these findings, it is essential to compare obtained results with ground-based or in situ observations every time a new model is introduced. This comparison is especially important for simulations related to space weather forecast infrastructure, as even slight changes in the inclination of magnetic structures may greatly affect the propagation of evolving CMEs at 1 AU and thus their geoeffectivity.

The validation process was already performed by Kuźma et al. (2023) for our two solar wind configurations in COCONUT. The validation scheme, originally described in Wagner et al. (2022), is meant to compare fully relaxed, steady-state solutions with the global coronal structures observed in white light. The scheme consists of a five-step process that includes visual classification, feature matching, streamer direction and width determination, brute force matching, and topology classification.

The only difference between the simulations of Kuźma et al. (2023) and those of our study is that we employed a finer grid resolution. For our simulations, we employed a grid that ranged from R = 1 R_⊙ to R = 25 R_⊙ with a variable cell size between 10⁻⁸ and 0.215 R³. In contrast, Kuźma et al. (2023) used a grid that extended radially outward in layers between R = 1 R_⊙ and R = 21.5 R_⊙ with a cell size that varied from 10⁻⁷ to 0.365 R³. Nevertheless, this does not affect the conclusions, as in both solar wind configurations, the equatorial structures are correctly placed, shaped, and inclined (cf. Fig. 1). The apex of streamers matches and the coronal holes are well defined. The only area where a mismatch is observed is the north polar region in the case of the maximum of the activity (cf. bottom panels in Fig. 1). That is, while magnetic structures are present in the eclipse observations, they are absent in the numerical results. However, the aforementioned issue generally arises due to insufficient coverage of the polar regions by photospheric magnetic field observations. This, in turn, leads to inadequate resolution of coronal magnetic structures in that region. We also note that this issue does not usually arise during the minimum of solar activity. At this time, the magnetic configuration tends to be more dipolar in nature and with few small magnetic structures near the polar regions. Such a configuration is therefore easier to reconstruct, even from low-quality observational data. Finally, we concluded that our steady-state model correctly recreates global coronal structures and that this was sufficient for studying CME propagation within its framework.

In addition to the validation conducted by observing the global coronal structures in white light, we propose an additional comparison that is based on computing the electron density of the coronal plasma using a coronal rotational tomography process. Originally developed by Morgan (2015), tomography involves constructing a density map at different heights by using the polarized brightness observations provided by the STEREO A COR2 coronagraph over a half-Carrington rotation period. Tomography is based on advanced data processing and calibration, as detailed in Morgan (2015), followed by a spherical harmonic-based regularized inversion method (Morgan 2019), which was later improved by Morgan & Cook (2020). Different tomography maps covering a broad range of years and heliocentric distance are publicly available online³.

Figure 2 shows a comparison of the electron density map at R = 8 R_⊙ obtained by COCONUT (left) with that obtained by the tomography process (right) for the minimum activity simulation (i.e., the solar eclipse of 2 July 2019). The tomography map shows a more detailed structure than the one of COCONUT, where the distribution of the structure is along nearly constant latitude bands. The highest densities in the tomography map were reached locally at longitudes equal to −50°, 15°, and 96°. In both cases, the structure is close to an equatorial single, narrow streamer sheet. However, the transition between the streamer belt and its surrounding regions is smoother in the COCONUT map than in the reconstructed density map.

Fig. 2. Comparison of the density values derived from COCONUT and those obtained through the tomographic process. Top panels: density maps of the solar corona at height R = 8.0 R⊙ for the total solar eclipse of 2 July 2019 in Carrington longitude-latitude coordinates. The left panel shows the density extracted from the COCONUT model, and the right panel shows the density obtained from the tomography process detailed in the Sect. 2.2. The color bars correspond to the density in normalized units. The dotted lines delineate isocontours for density values of [0.2, 0.4, 0.6, 0.8]. Bottom panel: position of the maximum of density in Carrington longitude-latitude coordinates. The orange curve is the maximum obtained from the COCONUT model, while the blue curve is the maximum from the tomography data. The position of the equatorial streamer belt, which has been reconstructed by COCONUT, is consistent with the one obtained by tomography.

We also observed that the equatorial structure’s width in the COCONUT map is significantly greater than that in the tomography map. Specifically, the width measures approximately 40° in COCONUT, while it only reaches a maximum of approximately 20° in the tomography map. This disparity in widths is directly related to the resolution of the grid used in the simulations. Brchnelova et al. (2022b) found that the current sheet’s shape and behavior are strongly influenced by the numerical dissipation resulting from the finite discretization of the system and the level of magnetic divergence present within the domain. Due to these two factors, premature reconnection occurs at an X-point, and the location of this point depends on the resolution. This early reconnection event causes an increase in the current sheet’s width as it progresses from the X-point. It is noteworthy that even with an increase in the number of cells (from 1.5 million to 2.3 million), this nonphysical behavior persists and is also observed in other solvers, such as Wind-Predict (Perri et al. 2022).

To compare the positions of the streamer belt, the bottom panel of Fig. 2 shows the latitude of the maximum density as a function of longitude for both COCONUT and the tomography process. The curves have been smoothed using a Savitzky-Golay filter (Press & Teukolsky 1990). The discontinuity observed near longitude and latitude 0° is a numerical artifact resulting from the method that was used to extract the maximum density. In Fig. 2, the dynamics of the two curves are almost identical. However, the COCONUT result appears slightly delayed compared to the tomography density before the longitude 0°. Additionally, the variations in the tomography density are wider (cf. Fig. 2, right panel). On average, the difference between the two curves is only 2.1°, with a maximum deviation of 6.7°. Despite the simplifying assumptions of the COCONUT model, which is only polytropic (cf. Sect. 2.1), the density is in good agreement with the observations. This agreement could be improved by incorporating more source terms in the COCONUT MHD equations, but this is beyond the scope of this work.

The comparison between the density of COCONUT and the tomography process was also performed at multiple heights (4, 4.4, 5, 5.4, 5.9, 6.5, 6.9, 7.5, and 8 solar radii) to have more confidence in the validity of the MHD results. However, the structure is almost identical between the different heights, as expected in the inner corona.

For the maximum activity studied in our work, no tomography map is available. The rapid changes related to the large CME occurring during this high period of activity led to a significant error in the “static” tomography method. Nevertheless, the results obtained for the minimum activity, along with the comparison with the white-light images, provided sufficient validation for the coronal magnetic field.

3. Flux rope in COCONUT

3.1. Titov-Démoulin modified model

The CME model implemented in COCONUT is an analytical circular flux rope originally developed by Titov & Démoulin (1999, hereafter TD) and then modified by Titov et al. (2014, hereafter TDm). The model depicts a circular cross section with a current-carrying, approximately force-free magnetic field in the toroidal segment. It is partly immersed in the photosphere such that only one circular arc protrudes into the corona.

In the original TD model, the current density is nearly uniform along the axis of the flux rope. In its modified version, the current density can instead be either (1) mainly localized close to the boundary or (2) parabolically distributed within the cross section, with a maximum at the axis and canceling at the boundary. As a result, the twist can either be concentrated in a thin layer or distributed over the whole torus.

The construction of the flux rope in the model is based on the distribution of the ambient field, which must be potential in the surrounding area. Specifically, the flux rope is placed such that its plane is locally perpendicular to the surrounding fields along one of the approximately circular isocontours. The magnetic field of this isocontour perpendicular to the toroidal axis is denoted as B_⊥. In this configuration, the flux is in equilibrium when the magnetic pressure resulting from the net current I ≠ 0 is balanced with the tension generated by the ambient field. In other words, the intensity must be equal to the Shafranov intensity, such as:

(2)

where R is the radius of the torus, a the minor radius, and I_i the internal self-inductance per unit length of the rope. The curvature of the flux rope is characterized by the ratio between the small radius and the radius of the torus:

(3)

In this configuration, the axial magnetic flux is defined as

(4)

for the density distribution number (1) and as

(5)

for the case (2). In the above equations, the sign is positive if the axial field and the current are counter directed and negative otherwise.

Assuming an axial symmetry of the magnetic field in the torus, it is possible to derive the magnetic vector potential from the current and the magnetic field defined above (cf. Sect. 2 in Titov et al. 2014). Finally, we note that all the characteristics of the flux rope can be determined using only the radius R, the minor radius a and an isocontour of the ambient magnetic field.

3.2. Implementation in the solar wind

As previously discussed (cf. Sect. 3.1), the flux rope is in equilibrium when its intensity is equal to the Shafranov intensity (cf. Eq. (2)). However, in this work, the intensity of the ring current is deliberately set higher than the Shafranov intensity in order to produce an eruptive behavior. Therefore, the intensity of the ring current is defined as:

(6)

with ζ as a positive number. By setting ζ to be greater than one, the flux rope is unstable from the start of the simulation. The overlying background magnetic field cannot prevent the outward expansion of the flux rope. As a result, there is no relaxation phase in which the photospheric patterns (e.g., shearing, twisting, or emergence of polarities) drive the magnetic system until the trigger of eruption due to magnetic reconnection or an ideal instability (Green et al. 2018).

Our simulation process was carried out in two steps. First, the solar wind background was created by a relaxation run (cf. Sect. 2.2). Once this was done, the flux rope was inserted into it, and COCONUT was restarted in its time-dependent version. During the insertion, only the magnetic field in the domain was modified. The density, the velocity, and the pressure were not changed. To insert the flux rope numerically, we used the Python library called pyTDm developed by Regnault et al. (2023), which is available online⁴. Originally developed for the PLUTO solver, the module has been adapted to work with COCONUT data files. The modification is also publicly available online.

The code pyTDm is designed to work only with structured mesh. As a result, an interpolation step was required to convert the unstructured COCONUT grid to the right format. The radial basis function (RBF) interpolation was the method used for this step (Buhmann 2000). This method produces accurate results even for grids with many cells. The choice of the interpolation method and its parameters was optimized to ensure that the solar wind produced by COCONUT after the relaxation run was only slightly modified during this procedure.

The module pyTDm offers the ability to specify the initial position of the flux rope (latitude, lat, and longitude, lon), the radius R, the minor radius a, the depth at which the torus is buried below the photosphere d (i.e., the distance between the surface of the Sun and the center of the torus), the tilt, and the parameter ζ. The current density was distributed along the case (1) described in Sect. 3.1 (i.e., mainly in a thin layer close to the boundary).

In Regnault et al. (2023), they defined specific boundary conditions around the flux rope at R = R_⊙ that were different from those of the solar wind boundary conditions. In contrast, in this study, there is no such distinction. It is also important to note that the insertion of the flux rope structure significantly changes the photospheric magnetic field. This results in the addition of two opposing magnetic polarities, as can be seen in Fig. 3, which illustrates the modification of the photospheric magnetic field after the injection of a flux rope in the minimum activity configuration. In the figure, the TDm is not positioned at the location of an observed active region, meaning that the inserted CME does not correspond to a CME that was possibly observed on the day of the photospheric magnetic field measurement.

Fig. 3. Snapshots of the radial magnetic field Br at R = 1 R⊙ in the COCONUT model (lower radial boundary corresponding to the base of the solar corona). Left panel: magnetic field produced by the initial relaxation run for the maximum activity solar wind. Right panel: same location after the insertion of the flux rope. The insertion of the flux rope results in the addition of two polarities and a modification of the near photospheric field. The blue arrow indicates the negative polarity, while the red one indicates the positive polarity. The magnetic inversion line is demarcated by the dotted line.

In Fig. 3, an inversion line between the two polarities can also be seen. This reflects the local presence of potential magnetic fields that overlay the flux rope. The magnetic field a bit further away from the insertion zone also seems to be modified. This is particularly noticeable in the left sunspot below the positive polarity. The magnetic field of this sunspot is lower after the injection of the flux rope.

Finally, we can mention that the disparities observed in Fig. 3 are the only contrasts between the photospheric magnetic field before and after its insertion. The areas that are not covered by the figure (mainly the opposite side) are identical. This means that the interpolation step that takes place after the execution of the pyTDm module and which puts the magnetic field configuration back into the unstructured grid of COCONUT does not impact the photospheric field and that the modification is therefore mainly related to the insertion of the flux rope.

3.3. Initial parameters

In this study, simulations were conducted to investigate the effect of initial conditions on the properties of flux ropes at 21.5 solar radii. A total of 24 simulations were run, 16 in a low-activity configuration and 8 in a maximum activity solar wind.

All the flux ropes were placed in the equatorial plane at longitude lon = 0° and latitude lat = 90°. Their centers were offset by d = 0.15 R_⊙ from the solar surface. Their major radii were R = 0.3 R_⊙ and their minor radii were a = 0.1 R_⊙. This resulted in a polarity area of A = 4839 Mm² on the inner boundary. The chosen size is deliberately larger than the typical size of an active region to ensure adequate numerical resolution for the description of the CME and to avoid interpreting it as noise. As shown in Fig. 4, after insertion of a flux rope in the maximum activity solution generated by COCONUT with the selected parameter set, only one arc protrudes into the solar corona.

Fig. 4. Visualization of the Titov-Démoulin modified flux rope in COCONUT in the maximum activity solar wind. The colored lines correspond to different magnetic field lines. The seed to trace the field lines is a sphere with a radius of 0.1 R⊙ that is placed at the positive polarity. The gray field lines are a sample of the surrounding field. The radial magnetic field at the solar surface (in Gauss) is also displayed. The flux rope is in the equatorial plane, and initially only an arc extends from the surface of the Sun.

For a given phase of the solar cycle (active or quiet), the difference between all the simulations is only due to variations in the parameter ζ, which affects the intensity of the initial flux rope. Tables 1 and 2 summarize the different cases examined for the minimum and maximum of activity. In the first magnetic configuration, the range of ζ varies from 5 to 20. As the magnetic field of the flux rope at its axis is directly related to the intensity, it increases from 4.2 G to 17.7 G. The magnetic field and the parameter ζ are related by the equation B₀(ζ) = 0.9ζ − 0.3. This linear relationship is a result of using the same surrounding magnetic field as a reference for all runs (cf. Sect. 3.2).

For the maximum activity configuration, the ambient magnetic field is stronger at the same height. Therefore, when the parameter ζ ranges from two to nine, the magnetic field increases from 9.8 G to 44.8 G. The relationship between the two quantities is B₀(ζ) = 5.1ζ − 1.1. As we later explain in Sect. 4.1, the flux ropes with a high magnetic flux lead to a high initial speed. To avoid convergence problems with large velocity gradients, the parameter ζ does not exceed nine in the maximum activity solar wind. In contrast, in the minimum activity cases, the ζ parameter does not go below five in order to maintain a sufficiently high initial velocity.

From the magnetic field, it is possible to compute the magnetic flux according to the Eq. (4). For the minimum activity, the magnetic flux varies from 1.6 × 10²¹ Mx to 7.4 × 10²¹ Mx. While for the maximum configuration, the flux is between 4.1 × 10²¹ Mx to 18.8 × 10²¹ Mx. The high values computed for the maximum activity correspond to major and rare solar events. The magnetic flux in the case of a minimum activity is closer to the values measured for large sunspots (van Driel-Gesztelyi & Green 2015).

4. Propagation in the corona

4.1. Early stages of the evolution

In order to track the propagation of all our flux ropes, we saved the 3D magnetic field and plasma flow outputs of COCONUT every 20 iterations corresponding to a time step of dt = 0.08 h, where t is the physical time in the simulation. All the cases were stopped when the physical time exceeded 24.04 h.

Figures 5 and 6 show the flux ropes at time t = 1.2 h in the solar wind for the minimum and the maximum activity, respectively. The different magnetic field lines of the torus are depicted with different colors and originate from a sphere with a radius of 0.1 R_⊙ situated at the feet of the flux rope with a positive polarity. It is worth mentioning that we tested different choices of seed for the magnetic field lines, but this did not change the main results presented in this section. We also tried adding a sphere at the negative polarity. However, the increased number of magnetic field lines made the figure less clear and did not provide any additional information. Finally, the sphere at the only positive polarity was deemed the most appropriate for portraying the flux rope geometry. Additionally, a sample of the surrounding magnetic field is also shown in white in Figs. 5 and 6. The only variation between the panels in each figure is the ζ parameter chosen during the implementation of the flux rope (cf. Sect. 3.2).

Fig. 5. Visualization of all the flux ropes modeled in the solar wind reconstructed from a minimum of activity at the physical time t = 1.2 h. The magnetic field lines are displayed according to the legend presented in Fig. 4. The small panel in the bottom-right corner of the simulation labeled “ζ = 12” is a zoom in on the polarities, and it highlights the presence of post-flare loops.

Fig. 6. Same as Fig. 5 but for the simulations in the maximum activity solar wind configuration.

The flux ropes in all the simulations exhibit a self-similar evolution, primarily in the equatorial plane. However, they do not reach the same height at the same instant. For a given solar wind, a higher ζ value resulted in a larger travel distance for the front of the flux rope. As expected, flux ropes with the highest ζ (i.e., highest intensity, cf. Eq. (2)) have the highest initial speed. Tables 1 and 2 show the radial velocity of the top part of the CMEs in the minimum and maximum activity, respectively, at the first saved time step, t = 0.08 h. We note that while comparing velocities at t = 0.08 h may not be entirely accurate, as the flux ropes had not reached the same height, it is still a good approximation. It is also important to note that the initial velocities are not numerically prescribed, they are exclusively generated by the imbalance of Lorentz forces that occurs right from the beginning of the simulations (cf. Sect. 2.3 in Scolini et al. 2019, for more details on role of the Lorentz force on CME propagation).

At the minimum of activity (cf. Table 1), the initial radial speed ranges from V₀ = 348 km s⁻¹ for ζ = 5 to V₀ = 1220 km s⁻¹ for ζ = 9. In comparison, the simulations at the maximum of activity cover a wider range of CME speeds, from V₀ = 455 km s⁻¹ for ζ = 2 to V₀ = 2385 km s⁻¹ for ζ = 9 (cf. Table 2). The radial propagation speed of observed CMEs typically varies between ∼200 km s⁻¹ and ∼2000 km s⁻¹ (Chen 2011). The majority of our simulations fall within this range. Only the two most unstable cases in the maximum activity have higher speeds, but they can be considered as rare and fast events that are only occasionally observed.

The higher the magnetic flux of the flux rope is, the higher its initial CME speeds are for a given solar wind. However, the velocity is not solely determined by the magnetic flux, as evidenced by the fact that the simulation with ζ = 11 in the minimum activity and the simulation with ζ = 3 in the maximum activity do not have the same propagation speed despite having a similar magnetic flux (≈4.0 − 4.1 10²¹ Mx). The density distribution and the surrounding field also play a role. There are various formulations to determine the propagation speed based on the characteristics of the solar corona. For instance, Shibata & Tanuma (2001) derived an analytical solution that depends on the Alfvén speed of the corona and the density of both the corona and the flux rope.

Regarding the topology of the flux ropes, there are no significant differences in geometry when comparing a simulation with a given ζ to the simulation with the subsequent ζ value (e.g., ζ = 5 and ζ = 6, ζ = 7 and ζ = 8, etc.) for the same solar wind. The increase in magnetic flux does not have a major impact on the geometry of the flux rope.

Figures 5 and 6 show common points in the topologies of the flux ropes as they grow into the solar corona. As they rise, the potential field lines surrounding them converge and pinch below the core of the flux rope. The convergence of these upward and downward field lines leads to the creation of a current sheet. This region is particularly suitable for magnetic reconnection. However, since there is no resistivity in the set of equations solved by COCONUT, the rearrangement of the field lines was done via numerical dissipation in our cases (cf. panel ζ = 12 in Fig. 5).

Reconnection leads to the formation of post-flare loops below the flux rope, as seen in observations (Schmieder et al. 1995). These field lines are located near the polarities of the Sun and rapidly draw closer to its surface. In addition to the post-flare loops, some field lines wrapping around the flux rope are also created, although they may not be clearly visible in Figs. 5 and 6 because of the choices made for the visualization. Finally, we note that the rearrangement of the line-tied magnetic field lines should remove some constraint on the overlying field and thus facilitate the expansion of the flux rope.

Not all the magnetic field lines become entirely rearranged. One part of the overlying magnetic loop was stretched upward by the advancing flux rope. The dynamics described above (e.g., the presence of post-flare loops) are described in the 2D “standard model” for CME originally developed by Carmichael (1964), Sturrock (1966), Hirayama (1974), Kopp & Pneuman (1976) and later extended in 3D by Aulanier et al. (2012, 2013), Janvier et al. (2013). Their model also accounts for physical processes not included in the COCONUT model, such as the presence of flare ribbons, the conversion of magnetic energy, and particle transport. COCONUT primarily focuses on the global distribution of the magnetic field density in the solar atmosphere, rather than the processes leading to flare appearance.

Comparing Figs. 5 and 6, we saw that the geometry of the flux rope is not the same in the two solar wind configurations. Further emphasizing these differences, Fig. 7 shows two magnetic flux ropes that propagate at approximately the same initial speed in the minimum and maximum activity configurations. For the minimum of activity, the CME has a ζ parameter of 17 and an initial speed of 1067 km s⁻¹, while in the maximum of activity it has a ζ parameter of four and an initial speed of 1060 km s⁻¹. Snapshots were taken at the first four time steps: t = 0 h, t = 0.4 h, t = 0.88 h, and t = 1.28 h.

Fig. 7. Visualization of the propagation of the TDm model in the solar wind reconstructed from a minimum of activity (top row) and from a maximum of activity (bottom row). For the two magnetic configurations, snapshots were taken at the beginning of the simulation and at physical times equal to 0.4 h, 0.88 h, and 1.28 h. The color legend is the same as in Fig. 4. The geometry of the flux ropes evolving in the maximum activity configuration is not identical to the one observed in the minimum activity scenario, meaning the solar wind configuration has an impact on the geometrical properties.

Initially, there is no noticeable difference between the two twisted flux ropes. As noted in Sect. 3.2, they have the same initial size and they are inserted at the same location on the Sun’s surface, with only the background solar wind differing between them. At time t = 0.4 h, the flux ropes have reached almost the same height, which is expected given their nearly identical initial speed. However, their shapes differ.

At the times t = 0.88 h and t = 1.28 h, the geometries of the flux ropes are significantly different. In the case of the flux rope evolving in the minimum of activity solar wind (cf. Fig. 7, top panels), a clear symmetry can be observed along the Z − X plane passing through the center of the torus. The difference between the two flanks can be numerically attributed to the fact that the origin of the magnetic field lines is a sphere placed at only the positive polarity. The geometry of the torus in this case is very close to the one observed by Regnault et al. (2023) for flux ropes evolving in a dipolar, ambient magnetic field. This similarity is expected, as the ambient field during a minimum of activity is typically close to a dipole field (cf. Fig. 1). The results produced by the COCONUT and PLUTO solvers are thus consistent, despite the differences in their numerical scheme and grid resolution.

However, in the maximum of activity (cf. Fig. 7, bottom panels), the symmetry along the meridional plane crossing the front part of the flux rope is less pronounced. The ambient magnetic field surrounding the flux rope is not the same on either side of it. For instance, the positive polarity is particularly close to two important sunspots (cf. Fig. 3). As the background magnetic field varies around the two legs of the flux rope, the interactions via magnetic reconnection between the field lines are also different. Additionally, the magnetic tension exerted by the surrounding field on the flux rope is not uniform along its length. The combination of these two factors leads to asymmetric evolution of the two legs.

We also observed that the magnetic flux rope appears thinner in the maximum activity than in the minimum activity. This difference comes from our choice of seed for tracing the field lines. In both cases and at all times, we only present magnetic field lines crossing a sphere located at the positive polarity. However, for the simulations in the maximum activity, some of the magnetic field of the CME no longer passes through the original sphere. This is something that Regnault et al. (2023) have also observed in their simulations. Indeed, after a few hours, the feet of the flux rope have slightly moved on the photosphere. The rearrangement of the field lines that compose them may lead to a shift in the anchoring area of the flux rope legs.

Further analysis, which is outside the scope of this study, is needed in order to determine more qualitative and quantitative changes of geometry caused by a different configuration of the solar wind, as well as the impact of visualization choices on the observed geometry. However, our first results suggest that the geometry of the flux is impacted by the surrounding magnetic field. In Sect. 5, we investigate the extent to which this difference in geometry affects the dynamics and values of various magnetic and thermodynamic quantities at 21.5 solar radii.

4.2. Known limitations from high initial speed

After the early stages of the evolution described in the previous section (cf. Sect. 4.1), the flux ropes continue to propagate in the solar corona. This propagation impacts all the MHD plasma and magnetic variables in the simulation domain. Figure 8 shows the density and the radial velocity in the equatorial plane for two very distinct simulations when the flux rope reaches around 15 R_⊙. The top panels correspond to the time t = 3.44 h for the flux rope propelling with an initial velocity of approximately 990 km s⁻¹ in the minimum activity, while the bottom panels show the results obtained at the time t = 1.12 h for the flux rope with the initial velocity of around 2385 km s⁻¹ in the maximum activity. In the 2D slice of the density (cf. Fig. 8, left panels), some field lines of the flux rope crossing a sphere of radius R = 1.5 R_⊙ at x = −16 R_⊙ are also displayed. To improve clarity, we have chosen not to display all 24 simulations. The features detailed above concerning density and the magnetic structures are valid for all cases studied for a given state of solar activity. Concerning the radial velocity, the right panels in Fig. 8 show examples for an extreme case (bottom panel) and for a more common case (top panel).

Fig. 8. Cross-sections along the equatorial plane of the density in log scale (left column) and of the radial velocity (right column). The top panels show the flux rope with ζ = 15 evolving in the minimum activity configuration, while the bottom panels correspond to the simulation with ζ = 9 in the maximum activity scenario. In the density panels, some magnetic field lines crossing a sphere of radius 1.5 R⊙ located at x = −16 R⊙ are also displayed. Each color corresponds to one field line. A virtual satellite was placed at 21.5 R⊙ (X = −21.5 R⊙, Y = 0, Z = 0) and is shown as a white sphere. The shape of the sheath is delimited by an ellipse in the density panels. In both solar wind configurations, a sheath is ahead of the flux rope, which is consistent with the observations. When the initial speed of the flux ropes is too high, there are nonphysical high-speed streams in the wake of flux rope legs.

The density distribution is affected in the same way in both solar wind configurations. In the direction opposite to the flux rope’s propagation, the density decreases with greater distance from the Sun. In the direction of the propagation, there is an increase in density by a factor of ten along an almost circular band surrounding the flux rope. As the flux rope travels, it does not allow enough time for the solar wind to flow around it, resulting in an area of heated and compressed solar matter called the sheath. This area is turbulent and characterized by a higher density than the surrounding solar wind. It is worth noting that the accumulation of matter ahead of the flux rope is a process that enhances as the CME expands (Siscoe & Odstrcil 2008). The amount of density reached in the sheath is discussed later, in Sect. 5.3. The accumulation of density can be observed by in situ spacecraft if the CME is faster than the solar wind (Regnault et al. 2020). This confirms the ability of COCONUT to reproduce features observed during the propagation of a magnetic ejecta.

Regarding the magnetic structure (cf. Fig. 8 left panels), no matter the solar wind and the initial velocity, it propagates radially in the equatorial plane. This suggests that there is no rotation of the flux rope. This is also the case for all simulations regardless of the distance from the Sun. Regnault et al. (2023) studied the propagation of the TDm model with the same initial geometry as our cases in a dipolar solar wind, and they found that the flux rope rotates by almost 55° under 15 solar radii. On the contrary, for thinner flux ropes (i.e., flux ropes with a lower minor radius), Regnault et al. (2023) did not observe any particular rotation. Finally, the differences in solar wind configuration between the study of Regnault et al. (2023) and our simulations can explain why our observations on the rotation of the flux rope are different.

Several processes may have caused a rotation to happen near the Sun. For example, Manchester et al. (2017) suggested that the rotation can be due to the kink instability, while for Shiota et al. (2010) the smooth rotation of the structure is caused by the reconnection with the surrounding field. Interactions with the overlying magnetic field can also lead to a deflection of the CME into a streamer. Using a 2.5D MHD simulation of the corona on a specific date, Zuccarello et al. (2011) observed a latitudinal migration of CMEs into a large helmet streamer due to an imbalance in the magnetic pressure and tension force. However, such dynamics are not observed in any of our simulations, despite the presence of important streamers in the maximum of activity configuration. Finally, further analysis is needed to determine whether the absence of rotation was anticipated or if it is a consequence of COCONUT being unable to accurately model the physical phenomena that could have caused it.

The initial magnetic flux, and thus the initial speed, has limited influence on the density perturbation generated by the propagation of the CME. In both cases presented in Fig. 8, a sheath is formed ahead of the flux rope. However, the distribution of the radial velocity is very different (cf. Fig. 8, right panels). In the slower simulation, the maximum of velocity is mainly concentrated near the core of the flux rope. The transition from the solar wind speed to the flux rope velocity is smooth. Behind the flux rope, the region of the current sheet can be distinguished, but there is no particularly high speed at the legs of the CME. This distribution of the speed is similar to the one observed by Maharana et al. (2022) when studying the evolution of a 3D flux rope model in EUHFORIA (see also Verbeke et al. 2019, for the same distribution with the spheromak CME model in the heliosphere).

In the faster simulation (cf. Fig. 8, bottom-right panel), the speed difference between the solar wind and the flux rope is much more pronounced. In this type of propagation, a piston-driven shock wave ahead of it is expected (Chen 2011). When comparing the right panels of Fig. 8, the ellipsoidal area where the solar wind speed has been modified by the propagation of the flux rope is shown to be quite similar in both cases. However, the speed distribution inside this area is not the same. In the bottom panel (the faster CME), the maximum of speed is reached at two high-speed flows starting from the Sun rather than at the bulk of the flux rope. These regions cover the CME legs. The speed inside these streams is typically higher than the initial speed.

According to the “standard model” for CMEs, two other fast-mode shocks could be observed. These shocks occur when the upward reconnection outflow collides with the flux rope or when the downward reconnection outflow collides with the post-flare loops. However, we suggest that the high-speed streams observed in our simulations are mainly numerical errors related to the difficulty encountered by the solver to handle high-speed gradients. Because of these poorly supported speed gradients, the solver produces regions where the pressure is numerically negative. As with Regnault et al. (2023), who faced similar problems in their fastest simulations, the pressure is set to 10⁻¹² in normalized units when it should be negative. We also attempted to mitigate this behavior by increasing the spatial resolution with a 2.3 million-cell grid. However, even with this higher resolution, we still observed the presence of these high-speed streams. It is also important to note that Regnault et al. (2023) encountered the same issue despite using an AMR grid, which offers better resolution at these high-speed stream locations. Finally, to limit the presence of these numerical artifacts, one solution could be to decrease the time step for better resolution, but this would be at the expense of computing time.

In conclusion, while a too high initial velocity does not seem to greatly affect the distribution of the magnetic structure and density, it does generate nonphysical dynamics. Therefore, future users should be aware of this and focus on flux ropes with initial speeds lower than 2000 km s⁻¹. High-speed gradients are only observed in the three fastest simulations (ζ = 7, ζ = 8, and ζ = 9). In all other simulations, regardless of the solar wind configuration, these high-speed streams are not observed. The speed distribution in these simulations is comparable to that of the top-right panel in Fig. 8. Finally, it is worth underscoring that there are only very few CMEs that have a speed higher than 1500 km s⁻¹ (Gopalswamy et al. 2009). COCONUT is therefore suitable for the study of the majority of cases observed.

5. Thermodynamic and magnetic properties at 21.5 solar radii

5.1. Magnetic field components

The simulations were conducted in order to demonstrate the potential use of COCONUT to provide heliospheric models with a realistic description of a CME. Some of these models, such as EUHFORIA, begin at 0.1 AU (R = 21.5 R_⊙). Therefore, in the following sections, the evolution of the plasma and the magnetic properties at a virtual satellite placed at x = −21.5 R_⊙, y = 0 R_⊙, z = 0 R_⊙ are described. This point is directly in the direction of the flux rope’s propagation (cf. the white sphere in Fig. 8).

The first quantity of interest is the magnetic field. Figure 9 shows the evolution of the magnetic field components (B_x, B_y, B_z) as well as the magnetic field amplitude, B, for simulations in the minimum of activity (the top panel) and in the maximum activity (the bottom panel) configurations. The x-direction corresponds to the horizontal axis in Fig. 8, the y-direction is the vertical axis, and the z-direction is perpendicular to both the x- and y-directions. In the solar minimum, the flux rope is initially implemented with ζ = 20 and with ζ = 2 in the solar maximum. These simulations were used to illustrate all the cases, as the dynamics described below can be found for all flux ropes. Only the amplitude of the different components changes between cases, not relative variations (cf. Sect. 5.2).

Fig. 9. Time evolution of the components of the magnetic field at 21.5 solar radii (X = −21.5 R⊙, Y = 0, Z = 0). The top panel corresponds to the flux rope with ζ = 20 in the solar wind background from a minimum of solar activity, while the bottom panel is the case ζ = 2 in the maximum of activity. The dark purple band corresponds to the sheath region, while the purple one is for the magnetic ejecta. The dynamics of the profiles are very similar in the two solar wind configurations: we observed a sheath followed by a magnetic ejecta composed of the flux rope initially inserted.

The two selected flux ropes do not have the same initial velocity. The CME speed is 1220 km s⁻¹ for the flux rope evolving in the minimum activity, while it is 455 km s⁻¹ for the CME in the maximum activity. Therefore, the CME in the minimum activity reaches the 0.1 AU boundary before the one in the maximum activity. This would be the opposite if a faster CME had been taken in the maximum activity.

Regarding the profiles, we emphasize that there is little difference between the simulation evolving in a solar minimum and the one evolving in a solar maximum. This suggests that the differences in shape and orientation observed in Sect. 4.1 do not result in a significant change in the magnetic properties of the flux rope when crossing the distance of 21.5 solar radii.

In particular, during the first hours, the magnetic field remained steady. It was equal to the field of the solar wind. Then, there was a fluctuation of the B_y and B_z components that led to an increase of the magnetic field amplitude (cf. Fig. 9). These fluctuations occurred earlier in the activity minimum since the flux rope evolved faster. This period reflects the crossing of the sheath (see the dark purple band in Fig. 9). At around t = 9 h in the maximum activity configuration, the total magnetic field, B, began to decrease before experiencing a second increase. As shown in Fig. 4 in Regnault et al. (2020), the transition between two bumps indicates the arrival of the magnetic ejecta for relatively fast events with a sheath. However, this pronounced change of slope was not observed in the minimum of activity (cf. top panel in Fig. 9). The observed difference in the sheath between simulations evolving in the minimum and maximum activity configurations can be attributed to the difference in the global magnetic field configuration. The sheath results from the compression of the magnetic field and of plasma by the propagation of the flux rope. As the solar wind configurations are different, the magnetic field profiles in the sheath also differ.

The magnetic field profile within the magnetic ejecta is asymmetric. The maximum is reached in the first half corresponding to the front of the magnetic structure. Because of the speed difference, the maximum is reached more quickly in the simulation with ζ = 20 (cf. Fig. 9, top panel). After the peak, in both cases, the magnetic field decreases before stabilizing. In the solar maximum, we note that the final value for the wake of the CME is higher than before its crossing, as expected. Indeed, studying a set of more than 300 profiles observed at 1 AU, Regnault et al. (2020) found that the most probable value for the total magnetic field is 4.9 ± 0.6 for the solar wind before the CME and 5.5 ± 1.0 when in its wake. This effect is enhanced for speeds that are significantly faster than the solar wind before the CME.

Regarding the magnetic field components (cf. Fig. 9), there is only a minor variation in the B_x component. According to the Lundquist model of flux rope, it should be zero while crossing the front part of the CME (Lundquist 1951). This suggests that the virtual satellite is crossing the flank of the CME or that there may be a very slight rotation not captured by tracing the magnetic field lines. However, this contribution is negligible compared to the other components. The evolution of the magnetic field is primarily dominated by the B_y and B_z components.

In the magnetic ejecta (light purple shaded area in Fig. 9), the B_y starts by increasing to its maximum, which coincides with the maximum of the total magnetic field. Then, the vertical component decreases until it reaches a negative value. In the final hours, B_y increases to reach a slightly higher value than before the event. For B_z (the green lines in Fig. 9), its dynamics are characterized by a change of sign. The evolution of B_z is opposite to that of the B_y component, as it begins by decreasing before rising again once its minimum is reached. In both solar wind configurations, the minimum and maximum peaks reached by B_z are almost identical in absolute value. The profile of B_z is asymmetric. Indeed, in the magnetic ejecta, B_z is more often positive than negative. It is worth noting that the reverse of sign of the B_y and B_z profiles can be measured by an in situ satellite at 1 AU (e.g., Fig. 1 in Regnault et al. 2020).

Finally, the variations of the magnetic components are in agreement with those expected by the TDm model (Titov et al. 2014). This supports the conclusion that the magnetic structure is preserved during its expansion in the solar corona. The rearrangement of the field lines discussed in Sect. 4.1 does not imply a total destruction of the flux rope in either solar wind configuration.

5.2. Density, magnetic field and velocity 1D profiles

In the previous section, we discussed how the magnetic ejecta that crosses the virtual satellite placed at R = 21.5 R_⊙ is indeed the flux rope that traveled from the solar surface. Now, we focus on the impact of the CME passage on the amplitude of the background magnetic field, velocity, and density in all our simulations (cf. Figs. 10 and 11). The components of the magnetic field, and the velocity, and the density are the quantities that are typically injected at the inner boundary of the heliospheric simulation to model a CME (e.g., Verbeke et al. 2019; Maharana et al. 2022).

Fig. 10. Velocity, magnetic field and density evolutions as a function of time at a virtual satellite placed at x ≈ −21.5 R⊙, y ≈ 0 R⊙, and z ≈ 0 R⊙. Each line corresponds to a different simulation in the background solar wind from a solar minimum. These profiles are consistent with the propagation of a flux rope with a sheath ahead of it.

Fig. 11. Same as Fig. 10 but for the simulations in the maximum of activity. The profiles are similar to those obtained in the minimum activity, meaning that COCONUT is effective in describing the propagation of the flux rope even during a solar maximum.

Figure 10 shows the evolution of the three quantities in all the simulations for the solar minimum. Each color corresponds to flux ropes implemented with a different ζ parameter. The figure shows that the changes caused by the passage of the sheath and the flux rope do not start at the same time in all our cases. For the fastest CME, the arrival time is approximately 3.8 h after the beginning of the simulation, while it is around 6 h for the slowest one.

The encounter with the sheath is characterized by an increase in density, magnetic field, and velocity, as expected by the in situ measurement. The increase is more abrupt for simulations with a high ζ parameter, while it is smoother for the slowest cases. After reaching a maximum whose value is discussed in the next section (cf. Sect. 5.3), the three quantities decrease. However, the subsequent evolution differs depending on the parameter.

As described in Sect. 5.1, after the magnetic field reaches its peak, it decreases before stabilizing around a value slightly higher than the pre-event solar wind value. For the density (cf. Fig. 10, last panel), the simulations with a smaller ζ reach a maximum before returning to a value close to the one before the crossing. For the simulations with a higher ζ (and therefore higher initial velocity), during a period of time that we associate with the crossing of the magnetic ejecta, the density is lower than before the arrival of the sheath. At the end of the simulations, the density increases slowly but does not return to its initial value within the next 10 h. Studying the density profiles at 1 AU, Regnault et al. (2020) found the same difference between simulations with relatively fast CMEs (compared to the speed of the solar wind) and slower CMEs. The distinctions observed in COCONUT are thus not surprising and are in agreement with the in situ measurements obtained further in the heliosphere.

We also note that the density bump is thinner than the magnetic field bump. The density takes less time to return to values close to those of the solar wind than the magnetic field. This supports the conclusion that the disturbance encountered at R = 21.5 R_⊙ (0.1 AU) is indeed composed of a sheath where the matter has accumulated and a magnetic ejecta corresponding to the flux rope that was initially implemented.

The evolution of the velocity profile is divided into two bumps (cf. Fig. 10, top panel). The first bump is the shock at the point of contact with the sheath, and the second is the arrival of the magnetic ejecta. The higher the speed of the shock is, the higher is the speed in the wake of the CME. The maximum in the wake is approximately 70 − 75% of the maximum speed of the shock. As for the density, in the relatively fast events, the speed of the solar wind post-CME is expected to be higher than before the CME (Regnault et al. 2020).

Figure 11 show the density, magnetic field, and velocity 1D profiles for the eight simulations performed in a solar maximum configuration. Except for the values of the peaks reached, the profiles of the three quantities for the simulations in a maximum of activity show similar trends to those obtained in the minimum of activity configuration. Most of the observations on profile variations are the same for both wind configurations, indicating that the difference in shape and orientation of the flux ropes (cf. Fig. 7) has little impact on the 1D profiles in the equatorial plane.

However, three differences can be observed between the simulation results in minimum and maximum activity configurations. First, as previously mentioned in Sect. 5.1, there is a change of slope in the magnetic field increase that can be attributed to the transition between the sheath and the magnetic ejecta. This change is not observed in the simulations evolving in a minimum of activity. Secondly, the density profile has a second local maximum that is above the pre-event value. Finally, the main difference concerns the evolution of the velocity. For the two fastest flux ropes (ζ = 8 and ζ = 9 in Fig. 11), the speed in the wake of the CME is higher than during the crossing of the magnetic ejecta. This behavior was already observed in the 2D slices of Fig. 8, where the speed maximum is reached close to the legs of the flux and not ahead of it. Moreover, at the end of the fastest simulations, there are fluctuations in the magnetic field, indicating that the solver has difficulty converging. The issue with convergence might be due to the fact that the modeled phenomena is inherently transient, and thus no immediate steady-state solution that the solver could converge to exists. In conclusion, this reinforces the idea that future studies should be limited to flux rope propagations with a velocity below 2000 km s⁻¹, as suggested in Sect. 4.2, in order to limit the nonphysical behavior.

5.3. Influence of the solar wind configuration

In Figs. 10 and 11, we focus on the dynamics of 1D profiles of density, magnetic field, and velocity at R = 21.5 solar radii. The main conclusion we drew was that the solar wind configuration has little influence on the evolution of the flux ropes considering polytropic solar wind. In this section, we take a closer look at the maximum values reached by these three quantities based on the initial parameters of the flux ropes and the solar wind configuration. Figure 12 shows the maxima of the magnetic field, velocity, and density as a function of the initial magnetic field strength B₀ of the flux ropes in both solar wind background configurations. The maximum values are also summarized in Tables 3 and 4.

Fig. 12. Maximum values of the different temporal profiles presented in in Figs. 10 and 11. From top to bottom: Maximum of the velocity, of the magnetic field amplitude, and of the density reached at x ≈ −21.5 R⊙ as a function of the initial magnetic field of the flux rope. The “+” markers are for the flux rope in the minimum of activity, while the orange dots correspond to the simulation in the maximum of activity. The dotted lines in orange and blue are polynomial regressions between the different values. The slope of the regressions depends on the solar wind configuration, meaning that the properties of the flux rope are related to its initial characteristics and to the surrounding magnetic field.

Our results show that, except for the two slowest flux ropes, the maximum speed is lower than the speed measured at the time t = 0.08 h. For example, in the case where ζ = 9 in the maximum of activity configuration, the speed has decreased by approximately 55%. The fact that two simulations have a higher speed at 21.5 solar radii than near the Sun indicates that the flux ropes accelerate during a given period. This acceleration could be caused by the momentum of the reconnection outflow (Shibata & Tanuma 2001) or by other factors, such as ideal MHD instabilities (Fan & Gibson 2007).

Our findings are consistent with the theory that suggests two forces act on a CME (Sachdeva et al. 2017). The first force is the Lorentz force, which allows for the radial expansion of the CME and whose profile peaks quickly and then declines. The second force is a drag force caused by the surrounding solar wind, which tends to hold down the expanding flux rope. This force explains why CMEs gradually slow down during their propagation. Over the years, the characteristics of these two forces have been studied analytically (e.g., Chen & Kunkel 2010; Subramanian et al. 2012) and by using observational data (e.g., Sachdeva et al. 2015, 2017).

In the top panel of Fig. 12, the maximum of speed, V_max, appears to increase linearly with the initial magnetic strength. To confirm this, we performed linear regressions in both solar wind configurations. In the minimum of activity, the formula we obtained is V_max ≈ 21.9 × B₀ + 277.1 with a coefficient of determination r² = 0.99. In the maximum of activity, the relation is V_max ≈ 23.8 × B₀ + 248.3 with a coefficient of determination of 0.99. For a low B₀ (lower than 15.16 G), the CMEs in the minimum of activity would be faster than CMEs in the maximum of activity. For large B₀ values (higher than 15.16 G), the flux ropes would be faster in the solar maximum. This is more consistent with the observations. Indeed, based on in situ measurements by the ACE satellite at the L1 Lagrangian point, Regnault et al. (2020) found that CMEs during active periods tend to be faster than during quiet periods (cf. also Hundhausen 1999).

The maximum of the magnetic field also appears to increase linearly with the initial field (cf. Fig. 12, middle panel). The correlation is B_max ≈ 28.3 × B₀ + 163.4 with r² = 0.99 in the solar minimum and B_max ≈ 22.3 × B₀ + 322.5 with r² = 0.98 in the solar maximum. Similar to the speed, the higher the magnetic flux of the flux rope is, the higher the magnetic field is at 21.5 solar radii. Comparing the initial magnetic field value to the one measured at the virtual satellite, we observed a decrease in the magnetic field during the propagation (cf. Tables 3 and 4). Different authors, such as Leitner et al. (2007), Liu et al. (2005), Winslow et al. (2015), have derived power laws in order to determine the decrease in magnetic field during the propagation phase. However, these power laws use measurements made at least at 0.3 AU away from the Sun, while our simulation grids stop at 25 R_⊙ (≈0.12 AU). Thus, the comparison between our simulation and the power laws deduced from in situ observation is limited. However, the latter turns out to work well for forecasting the magnetic field at 21.5 solar radii. Leitner et al. (2007), Liu et al. (2005), Winslow et al. (2015) suggested that the magnetic field amplitude of the CMEs is between 10³ and 10⁴ nanotesla at 0.1 AU, which is consistent with the maximum encountered in our set of simulations.

Unlike the other two quantities, the evolution of the maximum density cannot be fitted by a linear regression (cf. Fig. 12, bottom panel). This is simply due to the fact that no material is added during the implementation of the flux rope. The non-addition of material results in an increase in density at the crossing point that is entirely related to the accumulation of material in the sheath created by the propagation of the flux rope. As the mass is conserved in the domain, it cannot have an infinite increase. By fitting a polynomial regression in both solar wind configurations, we observed a horizontal asymptote that seems to be around 1.4 − 1.5 e⁻²⁰ g cm⁻³. Moreover, we observed that for a flux rope with the same magnetic field amplitude, the maximum of density reached is higher in the solar minimum than in the maximum activity. This is due to the fact that after the relaxation runs, the density for the minimum activity solar wind is slightly higher than for the maximum activity solar wind. The density distribution is thus not the same between the two configurations, and the accumulation in the sheath is also different.

The faster the CMEs propagate, the higher is the value of the density. But we also observed that the full width at half maximum of density decreases as the initial magnetic flux increases (cf. Figs. 10 and 11), as expected based on some analytical sheath models such as Takahashi & Shibata (2017). In order to quantify the compression, we derived an approximate size of the sheath by measuring the time between the beginning of the increase of the magnetized field and the first local minimum of the B_y component (cf. Fig. 9). The results are summarized in Tables 3 and 4 as well as in the Fig. 13. We found that the sheath size roughly fits a power law with a slope dependent on the solar wind configurations.

Fig. 13. Approximate size of the sheath (in hours) as a function of the initial magnetic field. The orange dots indicate the simulations at the maximum of activity, while the blue markers correspond to simulations at the minimum of activity. The orange and blue lines are exponential curves fitting the data. As expected, the width of the sheath decreases with increasing initial speed.

To summarize, the properties of the sheath and the flux rope seem to be directly related to the initial parameters in both solar wind configurations. However, the distribution of the surrounding magnetic field and density directly influences the value of the magnetic field, the speed, and the density of the flux ropes, as expected from the observations.

6. Summary and conclusions

The main objective of the present study was to demonstrate how the new 3D MHD coronal model COCONUT can be used to track the evolution of flux ropes in the corona. To this end, we first introduced the time-dependent implicit COCONUT solver (cf. Sect. 2.1). Before implementing a CME, it was necessary to create the surrounding magnetic field in which the flux rope would evolve. We thus ran COCONUT in its relaxation mode to obtain the solar wind from photospheric magnetograms for two opposing solar wind configurations: one corresponding to a solar minimum (2 July 2019, CR2219) and the other representing a solar maximum (20 March 2015, CR2161). In both cases, we used preprocessed HMI magnetograms as initial boundary conditions. We also verified that the configurations created were in agreement with white-light images and tomography measurements (cf. Sect. 2.2).

Once the background MHD corona was created, we could implement flux ropes using the method described in Sect. 3.2. The chosen model for the flux rope was the Titov-Démoulin model (Titov et al. 2014). We performed 24 different simulations in total, with 16 in the minimum solar activity corona and 8 in the maximum activity corona. These simulations varied only in the initial magnetic flux of the implemented flux ropes. In all cases, the propagation of the flux rope was driven by an imbalance caused by the Lorentz force.

We tracked the propagation of magnetic field lines by taking snapshots of them at different instants of our simulations (cf. Sect. 4). A virtual satellite was also placed at R = 21.5 solar radii to extract 1D time profiles of density, magnetic field, and velocity as the structures crossed this specific point in the domain (cf. Sect. 5). Finally, the main results of this study can be summarized as follows:

• During the first stages of the evolution, the shape of the flux rope is strongly affected by the surrounding field. In the minimum activity configuration, the flux rope has a symmetrical propagation, unlike cases evolving in the maximum activity solar wind. However, in both cases, we observed dynamics consistent with the standard eruption model, such as the development of a current sheet layer below the core of the CME. This region is conductive to magnetic reconnection (purely due to numerical dissipation in our cases), which we observed numerically by the presence of post-flare loops in the simulation domain.

• Further out in the corona, the flux rope continues to expand radially without any particular rotation and develops a sheath ahead of it, reflecting the accumulation of matter upstream of a piston-driven shock.

• In 2D slices of the radial speed, we find that artificial high-speed streams develop in the wake of the CME for the fastest flux ropes, leading to convergence problems. We suggest limiting future studies to flux ropes initially propagating at speeds lower than 2000 km s⁻¹, which correspond to most of the events we observed.

• Despite the differences in geometry, the evolution of the magnetic field components at 0.1 AU are similar in both solar wind configurations. They indicate the presence of a sheath followed by a magnetic ejecta composed of the initially implemented flux rope. The values of the solar wind are impacted by the wake of the CME after it passes by.

• The 1D time profiles of the density, magnetic field, and velocity obtained at 0.1 AU match related observations and further in situ measurement in the heliosphere. Only local peaks and the duration of the disturbances differ among the simulations.

• By comparing the maximums of density and of magnetic field strength as a function of the initial magnetic field strength in both solar winds, we find simple linear relationships between quantities. The greater the initial imbalance is of the Lorentz forces, the higher the magnetic field strength will be at 0.1 AU. These relationships also suggest that CME speeds are higher in the maximum activity than in the minimum activity case when the initial magnetic field strength is higher than 15.16 G. Due to the difference in the density distribution, the cases evolving in the minimum activity case have a higher peak of density than those in the maximum activity case.

• As expected, the size of the sheath depends on the initial velocity of the flux rope. The fastest CMEs result in a thinner sheath compared to the slowest cases.

COCONUT is an efficient solver for studying the propagation of a flux rope in a realistic solar wind configuration. Being faster than a state-of-the-art explicit solver, COCONUT appears as a valuable tool for understanding the physical mechanisms that occur during the initial CME propagation and for obtaining a detailed representation of a numerical CME at 0.1 AU after its interaction with the background solar corona reconstructed from observed magnetograms.

Using this work as a foundation, a logical next step could be to model a truly observed event. To do this, it would then be necessary to align the feet of the flux rope with the polarities of the active region. The following steps would then consist of ensuring that both the properties of the flux rope (its shape, intensity, etc.) and the properties of the simulated active region (its position, size, magnetic flux) coincide with what is measured at the photosphere and what is observed in the EUV in the low corona. It would also be necessary to modify the boundary conditions so that the polarities evolve similarly to what is observed at the surface of the Sun. Additionally, our analysis of the impact of the magnetic flux on the properties of the CME at 21.5 solar radii would enable a more precise selection of the flux rope’s initial parameters in order to match the characteristics of the real event. However, the analysis should be expanded to also include other parameters, such as the size and tilt of the flux rope.

We also note that, like the PLUTO solver (Regnault et al. 2023), COCONUT is currently unable to perform the relaxation of the flux rope before provoking its eruption. At this time, the superimposed flux rope is out of equilibrium from the start of the simulations. This is partly due to fixed inner boundary conditions that do not allow for consideration of the emergence and shearing of the magnetic field at the photosphere and the need to extend COCONUT to model the solar plasma below the corona (i.e., to include the transition region and chromosphere).

To achieve better results, particularly in future work, the solver is currently being improved through the addition of more physical terms related to the heat transfer, such as anisotropic thermal conduction, coronal heating, and radiation. Indeed, the current version of COCONUT only includes polytropic MHD, which leads to a biased distribution of density in the domain. In parallel, work is also underway to improve the balance of electromagnetic and gravitational forces in the domain. Combined, these adjustments would allow for a more realistic distribution of the magnetic field and density, as well as a bimodal distribution for the solar wind. In turn, this would lead to a better overall description of CME propagation. Another notable addition would be to consider a non-zero resistivity in the domain to obtain a better description of magnetic reconnection, which is currently entirely due to numerical dissipation. We remark, however, that the resistivity in the hot solar corona is extremely low, and current CPU power does not make it possible to increase the grid resolution such that the numerical resistivity is even lower.

COCONUT was originally developed to replace the simple empirical coronal model used by the space weather simulation model EUHFORIA (Pomoell & Poedts 2018). With improved initial conditions for the solar wind, EUHFORIA should yield more accurate predictions. By using our work as a foundation, we can also use COCONUT to model the evolution of flux rope CMEs closer to the Sun and then insert the CMEs into the inner boundary of EUHFORIA. However, two limitations of COCONUT must be taken into account to facilitate such a coupling of the models. First, the grid resolution of COCONUT decreases as one moves away from the Sun. Thus, at 0.1 AU, the spatial resolution is lower than the standard EUHFORIA resolution at the inner boundary of the heliospheric wind and CME evolution model. A spatial interpolation that could lead to information loss is therefore necessary for the coupling. To remedy this, an increase in resolution of the COCONUT grid can be considered, as well as the use of more advanced grid refinement techniques. This could also mitigate the presence of the high-speed streams mentioned earlier.

The second limitation is the computation time. To obtain our results, it took 2 h to create the background corona and between 20−26 h for the propagation of the flux ropes. With a total time of 22−28 h, COCONUT is a good first step toward a fully operational space weather asset, since CMEs typically take more than two days to reach the Earth. It is worth noting that COCONUT would be much faster than current state-of-the-art coronal models based on explicit solvers. Indeed, Perri et al. (2023) showed that COCONUT (with an implicit solver) is up to 35 times faster than PLUTO (using an explicit solver) to model a realistic solar wind configuration. Currently, the performance of the two solvers with the same initial set up was only made for the reconstruction of the background solar wind (i.e., using the steady-state implicit version of COCONUT). However, according to the results presented here, it is expected that the time-accurate version of COCONUT (with a flux rope inserted) would also be faster than PLUTO for studying the propagation of a CME.

In the future, the simulation model can be made even more efficient, as there are several ways to further reduce the computation time. One way to do so is to increase the resources dedicated to the code execution (e.g., the number of cores), as COOLFluiD scales well on parallel architectures. Additionally, it is possible to optimize the code parameters, such as the time step for the study of a particular event. In our work, we chose the same parameters for all simulations, but many of the flux ropes (the fastest ones) reach the boundary well before 20 h of computation. Moreover, a more detailed study of the impact of the time step combined with the grid resolution needs to be carried out to optimize CPU consumption.

Once the mentioned limitations are addressed and COCONUT is coupled with EUHFORIA, the effectiveness of the latter in predicting CME geoeffectiveness would be improved. Two approaches can be considered for inserting the flux rope that has evolved in COCONUT into the heliospheric wind and CME evolution module of EUHFORIA. The first approach is to extract at 0.1 AU only the properties of the flux rope and then to implement them into EUHFORIA, similar to how it is currently done for the spheromak CME model. However, this requires the ability to identify and isolate the CME structure. The second approach is to include the entire solar wind with the flux rope (i.e., to evolve the entire inlet boundary of EUHFORIA heliosphere) at 0.1 AU in time using COCONUT. In both cases, the inserted flux rope would have interacted with the solar wind before reaching 0.1 AU and should therefore have more realistic magnetic and geometric properties than the current self-similar models that are inserted at 0.1 AU without taking into account the initial evolution and interaction with the background corona. Overall, if coupled with a heliospheric model, COCONUT can be a valuable tool for space weather forecasting.

Acknowledgments

The authors thank the anonymous referee for their valuable comments. The authors acknowledge support from the European Union’s Horizon 2020 research and innovation programme under grant agreement N° 870405 (EUHFORIA 2.0). This work has been granted by the AFOSR basic research initiative project FA9550-18-1-0093. These results were also obtained in the framework of the projects C14/19/089 (C1 project Internal Funds KU Leuven), G.0B58.23N (FWO-Vlaanderen), SIDC Data Exploitation (ESA Prodex-12), and Belspo project B2/191/P1/SWiM. F.R. acknowledges grants 80NSSC20K0431, 80NSSC21K0463 and 80NSSC20K0700. The resources and services used in this work were provided by the VSC (Flemish Supercomputer Centre), funded by the Research Foundation – Flanders (FWO) and the Flemish Government. H.M.I. data are courtesy of the Joint Science Operations Center (JSOC) Science Data Processing team at Stanford University. We thank H. Morgan for the tomography data used.

References

All Tables

All Figures

Fig. 1. Comparison of the solar eclipse images (left) and numerical results (right) for the case of minimum (top) and maximum (bottom) of solar activity. The minimum activity image corresponds to the solar eclipse of 2 July 2019 (© 2019 Miloslav Druckmüller, Peter Aniol: http://www.zam.fme.vutbr.cz/~druck/Eclipse/Ecl2019ch/Tres_Cruses/TC_347mm/0-info.htm), and the maximum activity image corresponds to the solar eclipse of 20 March 2015 (© Miloslav Druckmüller, Shadia Habbal, Peter Aniol, Pavel Starha: http://www.zam.fme.vutbr.cz/~druck/Eclipse/Ecl2015s/800mm/00-info.htm). The red-blue sphere color map represents the radial component of the photospheric magnetic field expressed in Gauss, and the white lines represent magnetic field lines of the magnetic coronal structures in the plane of sight. The reconstructed main streamers align closely with the position and shape of coronal structures observed in white light.

In the text

Fig. 2. Comparison of the density values derived from COCONUT and those obtained through the tomographic process. Top panels: density maps of the solar corona at height R = 8.0 R⊙ for the total solar eclipse of 2 July 2019 in Carrington longitude-latitude coordinates. The left panel shows the density extracted from the COCONUT model, and the right panel shows the density obtained from the tomography process detailed in the Sect. 2.2. The color bars correspond to the density in normalized units. The dotted lines delineate isocontours for density values of [0.2, 0.4, 0.6, 0.8]. Bottom panel: position of the maximum of density in Carrington longitude-latitude coordinates. The orange curve is the maximum obtained from the COCONUT model, while the blue curve is the maximum from the tomography data. The position of the equatorial streamer belt, which has been reconstructed by COCONUT, is consistent with the one obtained by tomography.

In the text

Fig. 3. Snapshots of the radial magnetic field Br at R = 1 R⊙ in the COCONUT model (lower radial boundary corresponding to the base of the solar corona). Left panel: magnetic field produced by the initial relaxation run for the maximum activity solar wind. Right panel: same location after the insertion of the flux rope. The insertion of the flux rope results in the addition of two polarities and a modification of the near photospheric field. The blue arrow indicates the negative polarity, while the red one indicates the positive polarity. The magnetic inversion line is demarcated by the dotted line.

In the text

Fig. 4. Visualization of the Titov-Démoulin modified flux rope in COCONUT in the maximum activity solar wind. The colored lines correspond to different magnetic field lines. The seed to trace the field lines is a sphere with a radius of 0.1 R⊙ that is placed at the positive polarity. The gray field lines are a sample of the surrounding field. The radial magnetic field at the solar surface (in Gauss) is also displayed. The flux rope is in the equatorial plane, and initially only an arc extends from the surface of the Sun.

In the text

	Fig. 5. Visualization of all the flux ropes modeled in the solar wind reconstructed from a minimum of activity at the physical time t = 1.2 h. The magnetic field lines are displayed according to the legend presented in Fig. 4. The small panel in the bottom-right corner of the simulation labeled “ζ = 12” is a zoom in on the polarities, and it highlights the presence of post-flare loops.
In the text

	Fig. 6. Same as Fig. 5 but for the simulations in the maximum activity solar wind configuration.
In the text

Fig. 7. Visualization of the propagation of the TDm model in the solar wind reconstructed from a minimum of activity (top row) and from a maximum of activity (bottom row). For the two magnetic configurations, snapshots were taken at the beginning of the simulation and at physical times equal to 0.4 h, 0.88 h, and 1.28 h. The color legend is the same as in Fig. 4. The geometry of the flux ropes evolving in the maximum activity configuration is not identical to the one observed in the minimum activity scenario, meaning the solar wind configuration has an impact on the geometrical properties.

In the text

Fig. 8. Cross-sections along the equatorial plane of the density in log scale (left column) and of the radial velocity (right column). The top panels show the flux rope with ζ = 15 evolving in the minimum activity configuration, while the bottom panels correspond to the simulation with ζ = 9 in the maximum activity scenario. In the density panels, some magnetic field lines crossing a sphere of radius 1.5 R⊙ located at x = −16 R⊙ are also displayed. Each color corresponds to one field line. A virtual satellite was placed at 21.5 R⊙ (X = −21.5 R⊙, Y = 0, Z = 0) and is shown as a white sphere. The shape of the sheath is delimited by an ellipse in the density panels. In both solar wind configurations, a sheath is ahead of the flux rope, which is consistent with the observations. When the initial speed of the flux ropes is too high, there are nonphysical high-speed streams in the wake of flux rope legs.

In the text

Fig. 9. Time evolution of the components of the magnetic field at 21.5 solar radii (X = −21.5 R⊙, Y = 0, Z = 0). The top panel corresponds to the flux rope with ζ = 20 in the solar wind background from a minimum of solar activity, while the bottom panel is the case ζ = 2 in the maximum of activity. The dark purple band corresponds to the sheath region, while the purple one is for the magnetic ejecta. The dynamics of the profiles are very similar in the two solar wind configurations: we observed a sheath followed by a magnetic ejecta composed of the flux rope initially inserted.

In the text

	Fig. 10. Velocity, magnetic field and density evolutions as a function of time at a virtual satellite placed at x ≈ −21.5 R⊙, y ≈ 0 R⊙, and z ≈ 0 R⊙. Each line corresponds to a different simulation in the background solar wind from a solar minimum. These profiles are consistent with the propagation of a flux rope with a sheath ahead of it.
In the text

	Fig. 11. Same as Fig. 10 but for the simulations in the maximum of activity. The profiles are similar to those obtained in the minimum activity, meaning that COCONUT is effective in describing the propagation of the flux rope even during a solar maximum.
In the text

Fig. 12. Maximum values of the different temporal profiles presented in in Figs. 10 and 11. From top to bottom: Maximum of the velocity, of the magnetic field amplitude, and of the density reached at x ≈ −21.5 R⊙ as a function of the initial magnetic field of the flux rope. The “+” markers are for the flux rope in the minimum of activity, while the orange dots correspond to the simulation in the maximum of activity. The dotted lines in orange and blue are polynomial regressions between the different values. The slope of the regressions depends on the solar wind configuration, meaning that the properties of the flux rope are related to its initial characteristics and to the surrounding magnetic field.

In the text

	Fig. 13. Approximate size of the sheath (in hours) as a function of the initial magnetic field. The orange dots indicate the simulations at the maximum of activity, while the blue markers correspond to simulations at the minimum of activity. The orange and blue lines are exponential curves fitting the data. As expected, the width of the sheath decreases with increasing initial speed.
In the text