Issue 
A&A
Volume 690, October 2024



Article Number  A184  
Number of page(s)  12  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202449389  
Published online  07 October 2024 
The operationally ready full 3D magnetohydrodynamic model from the Sun to Earth: COCONUT+Icarus
^{1}
Department of Mathematics/Centre for mathematical Plasma Astrophysics, KU Leuven, Celestijnenlaan 200B, 3001 Leuven, Belgium
^{2}
Institute of Physics, University of Maria CurieSkłodowska, ul. Radziszewskiego 10, 20031 Lublin, Poland
Received:
29
January
2024
Accepted:
28
June
2024
Context. Solar wind modelling has become a crucial area of study due to the increased dependence of modern society on technology, navigation, and power systems. Accurate space weather forecasts can predict upcoming threats to Earth’s geospace and allow for harmful socioeconomic impacts to be mitigated. Coronal and heliospheric models must be as realistic as possible to achieve successful predictions. In this study, we examine a novel full magnetohydrodynamic (MHD) chain from the Sun to Earth.
Aims. The goal of this study is to demonstrate the capabilities of the full MHD modelling chain from the Sun to Earth by finalising the implementation of the full MHD coronal model into the COolfluid COroNa UnsTructured (COCONUT) model and coupling it to the MHD heliospheric model Icarus. The resulting coronal model has significant advantages compared to the preexisting polytropic alternative, as it includes more physics and allows for a more realistic modelling of bimodal wind, which is crucial for heliospheric studies. In particular, we examine different empirical formulations for the heating terms in the MHD equations to determine an optimal one that would be able to mimic a realistic solar wind configuration most accurately.
Methods. New heating source terms were implemented into the MHD equations of the preexisting polytropic COCONUT model. A realistic specific heat ratio was applied. In this study, only thermal conduction, radiative losses, and approximated coronal heating function were considered in the energy equation. Multiple approximated heating profiles were examined to see the effect on the solar wind. The output of the coronal model was used to onset the 3D MHD heliospheric model Icarus. A minimum solar activity case was chosen as the first test case for the full MHD model. The numerically simulated data in the corona and the heliosphere were compared to observational products. First, we compared the density data to the available tomography data near the Sun and then the modelled solar wind time series in Icarus was compared to OMNI 1min data at 1 AU.
Results. A range of approximated heating profiles were used in the full MHD coronal model to obtain a realistic solar wind configuration. The bimodal solar wind was obtained for the corona when introducing heating that is dependent upon the magnetic field. The modelled density profiles are in agreement with the tomography data. The modelled wind in the heliosphere is in reasonable agreement with observations. Overall, the density is overestimated, whereas the speed at 1 AU is more similar to OMNI 1min data. The general profile of the magnetic field components is modelled well, but its magnitude is underestimated.
Conclusions. We present a first attempt to obtain the full MHD chain from the Sun to Earth with COCONUT and Icarus. The coronal model has been upgraded to a full MHD model for a realistic bimodal solar wind configuration. The approximated heating functions have modelled the wind reasonably well, but simple approximations are not enough to obtain a realistic densityspeed balance or realistic features in the low corona and farther, near the outer boundary. The full MHD model was computed in 1.06 h on 180 cores of the Genius cluster of the Vlaams Supercomputing Center, which is only 1.8 times longer than the polytropic simulation. The extended model gives the opportunity to experiment with different heating formulations and improves the approximated function to model the real solar wind more accurately.
Key words: magnetohydrodynamics (MHD) / Sun: corona / Sun: heliosphere / solar wind
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Space weather has become a prevailing branch of physics, aimed at studying the environment around the Sun, with a major focus near Earth. Charged particles are continuously emitted from the Sun, following the magnetic field lines of the Sun, forming the socalled solar wind. The latter can interact with the magnetic field of Earth. Moreover, frequently, various eruptions such as coronal mass ejections (CMEs) happen on the solar surface. These CMEs are magnetised clouds that range between 10^{13} g up to 10^{16} g in mass and erupt with velocities ranging from ∼100 km s^{−1} to ∼3000 km s^{−1} (based on SOHO/LASCO measurements). Thus, they end up carrying enormous momentum and often causing strong shocks during their evolution in the inner heliosphere. The CMEs, when directed towards Earth, can cause significant damages, interfere with spaceborne operations and telecommunication and navigation systems, and provoke severe power outages. Similar damages can also be caused on Earth upon the interaction with the corotating interaction region (CIR) shocks (Alves et al. 2006). CMEs usually erupt only a few times during the solar minima, but during solar maxima, a few CMEs can occur during the day (Park et al. 2012). The damages coming from Earthdirected strong CMEs can be partially mitigated only if they are predicted sufficiently in advance. Space weather forecasting has become a standard approach to modelling the solar heliosphere and propagating CMEs in order to avoid disastrous scenarios on Earth. This requires the whole region from the Sun to Earth to be modelled. The dominant physics phenomena along the SunEarth path tend to vary. Modelling the whole domain with a single tool is certainly possible (Regnault et al. 2023), but since the implemented conditions strongly affect the time step in the simulation and vary from the solar corona to the heliosphere, separating them into two modelling tools is preferable. Moreover, since the domain is large and many grid points are necessary to cover the whole domain, the overall slowdown is significant regarding the simulation wallclock times. The common practice is to model the solar corona and heliosphere separately based on the physical conditions, since beyond the Alfvén point, all the information travels radially outwards and nothing goes back to the Sun. The natural separation point between the coronal and heliospheric models is the Alfvén point, which is conventionally assumed to be at 0.1 AU. Currently, a number of such forecasting tools exist. Examples of physicsbased forecasting tools include ENLIL (Odstrčil et al. 2004), EUHFORIA (Pomoell & Poedts 2018), SUSANOOCME (Shiota & Kataoka 2016), and CORHEL (Riley et al. 2012), as well as one of the most advanced magnetohydrodynamic (MHD) coronal model called MAS (Magnetohydrodynamic Algorithm outside a Sphere; Mikić & Linker 1996; Mikić et al. 2018; Downs et al. 2013), AWSoM (van der Holst et al. 2014), and WindPredict (Réville et al. 2016, 2020). These tools model the solar wind with superposed CMEs, propagating them to the Earth and beyond. Currently, ENLIL and EUHFORIA are used in operational settings, performing daily simulation runs to model the solar wind configuration and the CME evolution. The MAS and AWSoM codes have an advanced MHD corona, including a wave turbulence model introduced by Alfvén waves, as well as the transition region in the simulations. Such complex numerical MHD models usually require considerable CPU time to provide a suitable environment for testing different physics mechanisms in investigations of the coronal heating problem, the bimodal nature of the wind, evolution of complex flux rope models, and the acceleration of solar energetic particle (SEP) events (Roussev & Sokolov 2006).
Recently, a novel 3D MHD polytropic coronal model COCONUT (COolfluid COroNa UnsTructured) was developed at the KU Leuven Centre for mathematical Plasma Astrophysics (Perri et al. 2022; Brchnelova et al. 2022a). COCONUT has been implemented within the COOLFluiD (Computational ObjectOriented Libraries for Fluid Dynamics) architecture (Lani et al. 2005; Kimpe et al. 2005; Lani et al. 2013, 2014). COCONUT was developed to become the main MHD coronal model, along with the semiempirical WangSheeleyArge (WSA) model (Arge et al. 2004; Wang & Sheeley 1990), serving for forecasting purposes within the EUropean Heliospheric FORecasting Information Asset (EUHFORIA) 2.0 project (Pomoell & Poedts 2018). Thus, its domain extends from the solar surface to 0.1 AU. Various studies have been performed with COCONUT so far. In particular, the effects of preprocessing of the magnetograms, which are used for the inner boundary conditions, has been investigated in Kuźma et al. (2023), while the influence of the magnetogram types on solar corona simulation results was examined in Perri et al. (2023). In the present study, the results with a full MHD COCONUT model using a realistic adiabatic index (γ = 5/3) are demonstrated. New source terms are implemented into the MHD energy equation, representing the contribution from the radiative losses, thermal conduction and coronal heating. In the first approach, the coronal heating is approximated with different functions and the results are analysed. The output of COCONUT is used as the inner boundary condition for Icarus, a recently developed heliospheric modelling tool (Verbeke et al. 2022; Baratashvili et al. 2022). Icarus is a 3D MHD heliospheric model developed within the framework of MPIAMRVAC (Xia et al. 2018) and its domain covers radial distances from 0.1 AU to 2 AU, where the solar wind is modelled for the given date and CMEs can be injected in the simulation from the inner boundary. COCONUT and Icarus represent the first fully 3D MHD forecasting chain starting from the Sun to Earth within the EUHFORIA project. The produced solar wind profiles by COCONUT+Icarus chain are then validated by comparing to the observational data. The comparison is extensive and focusses on the effects and profiles of the solar wind both in the low corona, close to the Sun, and farther in the domain, near the Earth. Therefore, the profiles obtained by COCONUT are compared to tomography data, white light images, and OMNI data at 1 AU to confirm their results.
The paper is organised as follows. Section 2 describes the standard COCONUT numerical model, including the optimisation techniques to make such computationally expensive simulations feasible for forecasting purposes. Section 3 introduces the full MHD coronal model and discusses the additional source terms in detail. Afterwards, the obtained results with different approximated heating profiles are discussed with a comparison to the observational data in Sect. 3. Section 4 introduces the Icarus model and presents a discussion of the coupling between COCONUT and Icarus. Finally, our conclusions and outlook are presented in Sect. 5.
2. COCONUT Model
COCONUT is a 3D MHD model developed as a polytropic model in its first phase. As it was developed potentially as a part of the forecasting chain, the code is highly optimised to produce results fast enough to forecast space weather. To this end, COCONUT was developed as an implicit finite volume (FV) solver, enabling steadystate solutions up to 30 times faster than explicit solvers, as the Courant–Friedrichs–Lewy (CFL) numbers are not limited to 1, but can go up to 10 000 depending on the complexity of the magnetic field configuration (Perri et al. 2022). The use of unstructured grids also contributes to speeding up convergence since it allows for avoiding singularities near the poles, where local time steps tend to get very small and spurious fluxes might be generated.
In COCONUT, the ideal MHD equations are solved in conservative form in Cartesian coordinates, as described in detail in Perri et al. (2022):
$$\begin{array}{c}\hfill \frac{\partial}{\partial t}\left(\begin{array}{c}\rho \\ \rho \mathit{V}\\ \mathit{B}\\ E\\ \varphi \end{array}\right)+\mathbf{\nabla}\xb7\left(\begin{array}{c}\rho \mathit{V}\\ \rho \mathit{V}\mathit{V}+\mathsf{I}(p+\frac{1}{2}{\mathit{B}}^{2})\mathit{B}\mathit{B}\\ \mathit{V}\mathit{B}\mathit{B}\mathit{V}+\underline{\mathsf{I}\varphi}\\ (E+p+\frac{1}{2}{\mathit{B}}^{2})\mathit{V}\mathit{B}(\mathit{V}\xb7\mathit{B})\\ {V}_{\phantom{\rule{0.333333em}{0ex}}}^{2}\text{ref}\mathbf{B}\end{array}\right)\\ \hfill =\left(\begin{array}{c}0\\ \rho \mathit{g}\\ 0\\ \rho \mathit{g}\xb7\mathit{V}+\mathbf{S}\\ 0\end{array}\right),\end{array}$$(1)
where E is the total energy $\rho \frac{{V}^{2}}{2}+\rho \mathcal{E}+\frac{{B}^{2}}{\mathbf{2}}$, B is the magnetic field, V the velocity, g the gravitational acceleration, ρ the density, and p is the thermal gas pressure. The gravitational acceleration is given by $\mathit{g}(r)=(G{M}_{\odot}/{r}^{2})\phantom{\rule{0.166667em}{0ex}}{\widehat{\mathit{e}}}_{r}$ and the identity dyadic $\mathsf{I}={\widehat{\mathit{e}}}_{x}\otimes {\widehat{\mathit{e}}}_{x}+{\widehat{\mathit{e}}}_{y}\otimes {\widehat{\mathit{e}}}_{y}+{\widehat{\mathit{e}}}_{z}\otimes {\widehat{\mathit{e}}}_{z}$. The closure is given by the ideal equation of state, thus giving for the internal energy density $\rho \mathcal{E}=P/(\gamma 1),$ with a reduced adiabatic index of 1.05. The solar rotation is considered by prescribing a V_{ϕ} component at the inner boundary (Kuźma et al. 2023). Then, S in the energy equation stands for the source terms, which we introduce and explain in detail in Sect. 3. In the polytropic COCONUT simulations, we have S = 0.
We used an unstructured sixthlevel subdivided geodesic polyhedron mesh extended radially outwards in layers between r = R_{⊙} (inner boundary) and r = 30.0 R_{⊙} (outer boundary), as extensively explained in Brchnelova et al. (2022b). The mesh consists of 1.9 million cells with increasing cell size in the radial direction. Unlike in the original description of COCONUT in Perri et al. (2022), in this study, we used more consistent boundary conditions to reduce the generation of the artificial electric field, as described in Brchnelova et al. (2022a).
The typical parameters are prescribed at the solar surface, namely ρ_{⊙} = 1.67 × 10^{−16} g/cm^{3} and similarly for the temperature, T_{⊙} = 1.5 × 10^{6} K, are used for fixedvalue Dirichlet conditions of density and pressure. Then, from the ideal gas law, the surface value of the pressure can be obtained: P_{⊙} = 4.15 × 10^{−2} dyn/cm^{2} with the following formula:
$$\begin{array}{c}\hfill {P}_{\odot}=\frac{2{\rho}_{\odot}{k}_{B}{T}_{\odot}}{\mu {m}_{H}},\end{array}$$(2)
where ρ_{⊙} is the solar surface density, k_{B} is the Boltzmann constant, T_{⊙} is the solar surface temperature, μ is the molecular weight and m_{H} is the mass of hydrogen. In the original approach, μ ∼ 1.27 (Aschwanden 2005). According to Aschwanden (2005), the total mass density in the fully ionised gas consists of electron and ion densities. The contribution from the most abundant elements, hydrogen and helium, is also considered. However, in Aschwanden (2005), the number densities of ions are considered to be equal; whereas, in reality, they have different masses and charges, violating the n_{e} ∼ n_{i} assumption. When taking into consideration that m_{He} = 4 * m_{H}, we see that μ ∼ 1.27 corresponds to Helium abundance of ∼18%, which is a strong overestimation for the solar corona, especially since (in the most recent solar cycles) the Helium abundance is observed to be of the order of 1 − 2% (Yogesh 2021; Moses et al. 2020). Since the helium abundance is observed to be very small, in this approach, we fixed μ = 1, which corresponds to 0 helium abundance in the solar corona.
A Dirichlet boundary condition was prescribed for the radial component of the magnetic field B_{r}, which is directly derived from the magnetic map. The magnetograms represent the magnetic field configuration in the photosphere, where the field is much stronger than at the base of the corona. The magnetogram is preprocessed to smooth the magnetic field, which serves as the inner boundary condition in COCONUT. As suggested in Kuźma et al. (2023), we used spherical harmonics decomposition to filter the high spherical harmonics beyond a certain l_{max} value. This approach decreases the resolution of the magnetic field features and their strength.
Since the currently employed magnetograms provide only the information about the radial magnetic field component, we allowed the magnetic field in the other directions (B_{θ} and B_{ϕ}) to evolve freely throughout the simulation. We have seen that this approach leads to an accurate placement of electromagnetic features in the domain; as, for example, in the work of Kuźma et al. (2023). We also plan to experiment with vector magnetograms that would constrain all three magnetic field components instead. The velocity field at the inner boundary was set so the plasma follows the magnetic field lines, as described in Brchnelova et al. (2022a).
2.1. Choice of magnetogram
The input magnetogram can strongly influence the features of the obtained solar wind in the corona. The effect of using different magnetograms was studied by Perri et al. (2023). As a conclusion of this extensive study, the simulations based on Helioseismic and Magnetic Imager (HMI) produce the best results. Thus, in this study, only HMI magnetograms were used. However, we further looked into HMI products and decided to use the product that uses interpolations near the poles, since these are the regions with the fewest (and worst) observations. This product provides interpolated south and north poles of the Sun instead of data with significant noise.
The left figure in Fig. 1 shows the original HMI product choice. The upper panel shows the HMI magnetogram, while the bottom panel shows the processed input magnetogram for COCONUT at the solar surface. The processing is done as in the original description in Perri et al. (2022). l_{max} = 20 is used for spherical harmonic decomposition. The units are given in Gauss, while in the processed image, the values are code units, meaning the magnetic field is divided by 2.2 G. Positive values can be noticed near the south pole in the processed magnetogram, while that trend is not observed in the original magnetogram. Also, when looking at the original magnetogram, strong noise is present near the poles.
Fig. 1. Synoptic maps for July 2, 2019, Carrington Rotation 2219. The original (left) and pole filled (right) HMI magnetograms are provided in the upper panel. The bottom panel shows the processed corresponding input magnetogram with spherical harmonics decomposition filtered with l_{max} = 20. The units are in Gauss. In the processed magnetogram, the values are given in the code values (divided by 2.2 G). 
The right figure in Fig. 1 represents the pole interpolated HMI magnetogram product. Here, the strength of the magnetic field is higher because it is a higherresolution map than in the original HMI magnetogram case. After processing, we see that the strength of the obtained input boundary condition file is similar to the original one. The main features in the midlatitudes also remain unaffected, where we see notable differences near the poles. The poles here are more homogeneous, without any unexpected positive values in the South Pole.
In comparing the two products, the pole interpolated HMI product in Fig. 1 was chosen for this study to avoid artificially generated artefacts in the modelled corona.
2.2. Polytropic corona for CR 2219
First, we present results for the chosen Carrington Rotation with the polytropic COCONUT model, which was the standard COCONUT model used in all previous works. In this simulation, we used the magnetogram shown in the right panel of Fig. 1 as the initial magnetic field configuration at the solar surface. Figure 2 shows the obtained solution for the radial velocity configuration in the meridional plane. The speed distribution is uniform near the outer boundary, namely, ∼0.1 AU with 350 km s^{−1} values.
Fig. 2. Polytropic solution for the CR2219 corresponding to July 2, 2019. The radial velocity is plotted in the meridional plane. 
Figure 3 shows the same plane as the previous figure but zoomed in on the low corona. The eclipse image is plotted in the background and the magnetic field lines are overlaid by the polytropic COCONUT simulation. The solar surface is coloured with the radial magnetic field values. The large structures seen in the low corona are well mimicked by the magnetic streamlines. This is expected, since plasma is believed to be trapped by the closed loops, enhancing the density in these areas (Lionello et al. 2001).
Fig. 3. Polytropic solution for the CR2219 corresponding to July 2, 2019. The eclipse image is overlaid with the magnetic streamlines from the simulation. The surface of the Sun is coloured with the B_{r} in Gauss. 
Furthermore, we also compared the modelled data to available tomography data. Coronal rotational tomography is used to estimate the electron density in the coronal plasma (Morgan 2015, 2019, and Morgan & Cook 2020). The available data is based on the STEREOA COR2 coronagraph polarised brightness observations. The density profiles are given from 4 to 8 solar radii. When comparing the results from COCONUT simulations to the topography data, the comparison to the similar simulation results from the MAS model (Mikić & Linker 1996; Lionello et al. 2014) was also considered. MAS is a timedependent resistive thermodynamic MHD model, where the MHD equations are solved on a nonuniform, logically rectangular staggered grid using finite differences. MAS simulation results are available freely through the website www.predsci.com. Therefore, we chose the data corresponding to the magnetogram for the 2019 eclipse and visualised the results for the polytropic MAS model from the website.
Figure 4 shows the results from the COCONUT simulation (left) and the MAS model (right). The first two rows show the radial magnetic field and number density at 5 R_{⊙} in the polytropic COCONUT and MAS simulations. The third panel shows the identical tomography reconstruction at 5 R_{⊙} from the observed data. The density is given in normalised units to better emphasise the structures obtained at 5 R_{⊙}. This figure shows that the density profile closely resembles the current sheet in the COCONUT results. This is expected since HCSs are characterised by the slow and highdensity solar wind (Réville et al. 2023; Poirier et al. 2021). The density enhancement position agrees with the tomography data; however, the detailed structure is missing. We spot the same behaviour in the MAS simulation results, where the data is in higher resolution, and the current sheet has a different shape. However, the density enhancements and the current sheet are also similar in the MAS model.
Fig. 4. Horizontal axis shows longitudes in degrees and the vertical axis shows colatitudes in degrees. The left and right panels show the B_{r} and density in normalised units from the COCONUT and MAS models, respectively. The bottom row shows the observed tomography data in normalised units. All quantities are plotted at 5 R_{⊙}. Data from COCONUT are taken from the polytropic simulation. 
MAS and COCONUT are different codes with different numerical setups, namely, using different grids and numerical methods and boundary conditions (Lionello et al. 1999), as mentioned earlier in this paper. This comparison aims to verify whether these two different solvers would produce similar results. When comparing the features present in the radial magnetic field and, consequently, in the number density plots, we see the enhancements and features present at the same places. The current sheet has a different shape in the MAS simulation, which is partially because the processing of the HMI magnetogram is different and their mesh is much finer, which leads to less numerical dissipation, together with the other abovementioned differences.
3. Full MHD coronal model
The real processes happening near the solar surface that propagate towards the Earth are very complex and not fully understood. Since our goal is to include more physics phenomena occurring near the solar surface, we implemented and tested some empirical source terms that are aimed at approximating the physics of the solar corona more closely. To this end, we followed the paper by Mikić et al. (1999), Lionello et al. (2009), and Downs et al. (2010). In the upgraded MHD model, we fixed the adiabatic index to γ = 5/3 and used the same inner boundary conditions as in the polytropic COCONUT simulations. Such source terms are introduced in the energy equation, according to the following formulation:
$$\begin{array}{c}\hfill S=\mathrm{\nabla}\xb7\mathbf{q}+{Q}_{\mathit{rad}}+{Q}_{H},\end{array}$$(3)
where Q_{H} is the coronal heating, Q_{rad} is the radiation loss function Mikić et al. (1999) and −∇⋅q models the thermal conduction. In this study, we neglect resistivity and viscosity and do not introduce these terms in our MHD equations, as their contribution must be small, and we focus on the dominant terms first.
Thermal conduction was implemented the same way as in Mikić et al. (1999) following Hollweg (1978). We defined two regimes where plasma is collisional and collisionless. We define thermal conduction as follows:
$$\begin{array}{c}\hfill {\mathit{q}}_{\mathbf{1}}={\kappa}_{}\widehat{\mathit{b}}\widehat{\mathit{b}}\xb7\mathrm{\nabla}T,\end{array}$$(4)
where q_{1} is the standard SpitzerHärm thermal conduction flux in the collisional regime below 10 R_{⊙} with ${\kappa}_{}=9\times {10}^{7}{T}^{\frac{5}{2}}$ in cgs units, as follows:
$$\begin{array}{c}\hfill {\mathit{q}}_{\mathbf{2}}=\alpha {n}_{e}kT\mathit{v},\end{array}$$(5)
where q_{2} is the thermal conduction flux in the collisionless regime beyond 10 R_{⊙}, α is a constant given in Hollweg (1978) and k is the Boltzmann constant. Radiative loss is defined in the optically thin limit according to Rosner et al. (1978) by:
$$\begin{array}{c}\hfill {Q}_{\mathit{rad}}={n}_{e}{n}_{p}P(T),\end{array}$$(6)
where n_{e} and n_{p} correspond to electron and proton number densities, and it is assumed that n_{e} = n_{p} for the hydrogen plasma, P(T) is defined in Rosner et al. (1978) and is a cooling curve depending on the temperature. The defined radiative loss function profile is given in Fig. 5 in loglog scale.
Fig. 5. Radiative loss cooling curve defined as in Rosner et al. (1978). The horizontal axis denotes temperature and the vertical axis shows the radiative loss cooling curve function in cgs. Loglog scale is applied. 
The mechanism that heats up the corona is still unresolved in solar physics. It is assumed that magnetic energy is transformed into thermal energy, but the exact mechanism still has to be identified. The coronal heating has been approximated with physicsbased models but also with phenomenological approaches (Schrijver et al. 1985; Fisher et al. 1998). In our model, we use the three most common approximations found in the literature. As a first approximation, we used an exponential envelope function similar to that of Lionello et al. (2009) and Downs et al. (2010),
$$\begin{array}{c}\hfill {Q}_{H}={H}_{0}{e}^{\frac{r{R}_{s}}{\lambda}},\end{array}$$(7)
where R_{s} is the solar radius, H_{0} is the local heating rate at r = R_{s} and λ is the scale height. We use H_{0} = 4.9 ⋅ 10^{−5} erg cm^{−3} s^{−1} and λ = 0.7 R_{⊙}.
Pevtsov et al. (2003) established that there is a linear dependence of the magnetic field strength and Xray radiance. Therefore, a function approximating this law is tested in COCONUT with a slight modification. Downs et al. (2010) also suggested using the radial damping of the heating term when using the coronal heating model based on the magnetic field configuration. Thus, the second model considered here is expressed as:
$$\begin{array}{c}\hfill {Q}_{H}={H}_{0}\xb7\mathbf{B}\xb7{e}^{\frac{r{R}_{s}}{\lambda}},\end{array}$$(8)
where H_{0} = 4 ⋅ 10^{−5} erg cm^{−3} s^{−1} G^{−1} and λ = 0.7 R_{⊙}.
Finally, we also considered the most complex heating function approximation that considers exponential heating, which is the contribution describing the quiet Sun and active region heating. The last approximated heating function is taken from Lionello et al. (2009). The function depends on the magnetic field, approximating the heating for the quiet Sun and the active regions:
$$\begin{array}{cc}\hfill {Q}_{H}& ={Q}_{H}^{\mathit{exp}}+{Q}_{H}^{\mathit{QS}}+{Q}_{H}^{\mathit{AR}},\hfill \end{array}$$(9)
$$\begin{array}{cc}\hfill {Q}_{H}^{\mathit{exp}}& ={H}_{0}{e}^{\frac{(r{R}_{\odot})}{{\lambda}_{0}}},\hfill \end{array}$$(10)
where H_{0} = 4.9128 * 10^{−7} erg cm^{−3} s^{−1} and λ_{0} = 0.7 R_{⊙}.
$$\begin{array}{c}\hfill {Q}_{h}^{\mathit{QS}}={H}_{0}^{\mathit{QS}}f(r)\frac{{B}_{t}^{2}}{B({B}_{r}+{B}_{r}^{c})},\end{array}$$(11)
where ${B}_{t}=\sqrt{{B}_{\theta}^{2}+{B}_{\varphi}^{2}}$ is the tangential magnetic field, and
$$\begin{array}{c}\hfill {Q}_{H}^{\mathit{AR}}={H}_{0}^{\mathit{AR}}g(B)(\frac{B}{{B}_{0}}{)}^{1.2}.\end{array}$$(12)
In the performed simulations, the following constants are fixed according to Lionello et al. (2009), H_{0}^{QS} = 1.18 * 10^{−5} erg cm^{−3} s^{−1}, B_{r}^{c} = 0.55 G, H_{0}^{AR} = 1.87 * 10^{−5} erg cm^{−3} s^{−1}, B_{0} = 1 G. The functions f(r) and g(B) are defined as follows:
$$\begin{array}{cc}\hfill f(r)& =\frac{1}{2}(1+tanh\frac{1.7r/{R}_{\odot}}{0.1})exp(\frac{r/{R}_{\odot}1}{0.2}),\hfill \end{array}$$(13)
$$\begin{array}{cc}\hfill g(B)& =\frac{1}{2}(1+tanh\frac{B18.1}{3.97}),\hfill \end{array}$$(14)
The heating profile in Eq. (9) takes into account different heating contributions, but with its initial definition in Lionello et al. (2009), it is specifically tailored for the test case considered in the original paper; therefore, is has required further modification to adjust to the different solar wind configuration. For our case, we changed the cutoff strength for the magnetic field to 4 Gauss, as the magnetic field is weaker. Thus, in COCONUT, we use the following formula for the active region heating:
$$\begin{array}{c}\hfill g(B)=\frac{1}{2}(1+tanh\frac{B4}{4}).\end{array}$$(15)
Figure 6 shows the obtained heating profiles in the meridional plane in the COCONUT simulations for each indicated heating function. The values are given in [J m^{−3} s^{−1}] units. The left figure shows the exponential heating function given in Eq. (7), and, as expected, the obtained heating profile is uniform and decreases radially outward. The total power injected by this function is 7.13 ⋅ 10^{19} J s^{−1}. The heating profile introduced by Eq. (8) is given in the middle panel, where the profile is no longer uniform, as the magnetic field contribution is considered. The total power injected in the corona with this function is 2.88 ⋅ 10^{18} J s^{−1}. The heating profile introduced in Eq. (9) is given in the last panel of Fig. 6. The different heating contributions are separated here, and the heating coming from the strong magnetic field regions is dominant in this case. The total power injected in the corona with this function is 1.44 ⋅ 10^{21} J s^{−1}, which is the highest among the implemented heating profiles. The total power injected by the heating profiles is reported to be of the order of 10^{19} − 10^{21} J s^{−1} in the literature (Lionello et al. 2009; Downs et al. 2010; Réville et al. 2020; Parenti et al. 2022).
Fig. 6. Obtained heating profiles in the meridional plane for the full MHD COCONUT simulations for each heating description given in Eqs. (7)–(9), from left to right, respectively, in [J m^{−3} s^{−1}]. 
Figure 7 shows the obtained heating functions in the meridional plane in the COCONUT simulations for each indicated heating function. The velocity profile near the outer coronal boundary obtained by Eq. (7) is uniform but still very fast. The wind profile obtained with Eq. (8) is bimodal, as the speed is slower near the equator than at the poles. The solar wind speed profile is also almost bimodal in the last panel of Figure 7. However, the slow speed stream is small at 0.1 au in this case.
Fig. 7. Radial velocity values in [km s^{−1}] in the meridional plane for the full MHD COCONUT simulations with the indicated heating functions given in Eqs. (7)–(9), from left to right, respectively. 
Figure 8 shows the overview of the effect of the different heating profiles in the low corona. In all the figures, the observed eclipse image is plotted in the background and magnetic streamlines are overlaid from COCONUT simulations. The solar surface is coloured with the radial magnetic field values.
Fig. 8. Full MHD COCONUT solutions with the indicated heating profiles given in Eqs. (7)–(9), from left to right, respectively. The eclipse image is overlaid with the magnetic streamlines from the COCONUT solution. The radial magnetic field is plotted on the solar surface. 
The left figure shows the profiles for the heating function given in Eq. (7). The heating introduced by this equation with the given scale height values is uniform yet strong. The field lines seem to be blown out and radially stretched – and not collimated towards the equator. The second panel represents the obtained solar wind with Eq. (8). The stream lines are elongated compared to the fieldlines in the polytropic case. Here, the injected heating in the corona was sufficient to speed up the solar wind and obtain the bimodal structure. The right figure shows the obtained coronal configuration for Eq. (9). The obtained heating profile also accelerated solar wind sufficiently. The total power injected by the heating functions in the domain is not overestimated in COCONUT simulations. In the first two cases, it is underestimated. However, we can see that the energy deposition in a specific place plays a crucial role in obtaining a bimodal wind configuration near the outer heliospheric boundary. For example, the total power injected by Eq. (8) is less than in the case of Eq. (9). Nevertheless, the bimodal structure of the solar wind achieved by the former is clearer than by the latter heating profile. We also notice that the deposited heating strongly affects the loops’ shape. A similar effect to the stretched magnetic loops was also obtained in the Wind Predict simulation in Parenti et al. (2022), where the authors observed the stretching of the helmet streamers depending on the location and strength of the deposited thermal energy, together with the number density configuration at the base of the corona.
Figure 9 shows the radial magnetic field, density and tomography data at 5 R_{⊙} from the COCONUT simulation with the full MHD model based on heating functions given in Eq. (7) in the left panel. The middle panel shows the results with the heating from Eq. (8) and the right one includes the heating function from Eq. (9). The structures of the density profiles in all figures show more features compared to the polytropic case (Fig. 4), which again is linked to the current sheet. The density is more enhanced near the equatorial plane in COCONUT simulations than in the tomography data, partially due to the resolution limitation and the prescribed density at the base of the corona. Nevertheless, the distinctive features present in the simulation results (middle row) between [50 ° −180°] degrees in longitude coincide with the observation tomography data (bottom row), as the enhanced density surrounds the lower density region. The distinct contributions from the heating function in Eq. (9) modelled more realistic density profiles in the corona, based on the comparison of the number density profiles between [50 ° −180°] degrees in longitude. The shape of the high values in the density profiles; namely, the highdensity profiles that appear at 10° latitude and 50° longitude and shift towards +20° latitude as the longitude values increase towards 180° degrees are more similar to the tomography observation data plotted in the bottom panel. This signature is absent in the result obtained with the heating profile given in Eq. (8). The number density modelled by the heating profile in Eq. (7) demonstrates a similar profile to the tomography data. The enhanced density is more localised near the equator. However, the prominent feature of the lowerdensity region near the equator surrounded by the higherdensity region is missing. Figure 10 shows the radial magnetic field and number density at 5 R_{⊙} from MAS simulation with the thermodynamic model 1 available on their website. This model includes the transition region and the wave turbulence model, in contrast to the COCONUT simulation. We can see that the features are present in the same place compared to Fig. 9, and the density enhancement is also more prominent than the tomography data. The result of the complex heating function implemented in COCONUT yielded a more similar number density profile to the MAS simulation results, which could be because the additional heating in MAS is introduced with the same formula from Lionello et al. (2009).
Fig. 9. Longitudes in degrees shown on the horizontal axis. The vertical axis gives the colatitudes in degrees. Here, B_{r} is given in [nT], the density from COCONUT in normalised units, and tomography in normalised units, from top to bottom, respectively. All quantities are plotted at 5 R_{⊙}. The left panel shows results for the COCONUT simulation with the heating function from Eq. (7), the middle panel shows the results with the heating from Eq. (8), and the right includes the heating function from Eq. (9). 
Fig. 10. Longitudes in degrees are shown on the horizontal axis. Vertical axis gives the colatitudes in degrees. The panels represent Br in [nT], the density from MAS simulations in normalised units, and tomography in normalised units from top to bottom, respectively. All quantities are plotted at 5 R_{⊙}. The MAS simulation was obtained with thermodynamic model 1. 
Table 1 summarises the computational resources spent on the simulations. All COCONUT simulations were performed on four nodes of the Genius cluster of Vlaams Supercomputing Center. During the iterative process, the run times were taken when the residuals for the density components reach 10^{−3}. The first row represents the polytropic simulation, which takes ∼0.98 h. The secondtolast rows give the results for the full MHD COCONUT simulations with different heating profile approximations. The heating profile introduced by Eq. (7) required 0.7 h to converge to the steadystate solution. The heating profile given by Eq. (8) that resulted in the successful bimodal solar wind took 1.2 h, which is ∼0.2 h longer than the polytropic simulation. The heating profile approximated by Eq. (9) took 1.94 hours to converge, which is only ∼1 h slower than the polytropic simulation. As a result, obtaining the bimodal solar wind by activating the full MHD source terms does not require a significantly increased computational time, making it suitable for space weather forecasting.
Run times (wallclock time) required for the COCONUT simulations to reach convergence –3 in density.
4. Coupling to Icarus
Icarus is a 3D MHD heliospheric modelling tool (Verbeke et al. 2022; Baratashvili et al. 2022) developed in the framework of MPIAMRVAC (Xia et al. 2018). Icarus was developed as an alternative heliospheric model within the EUHFORIA space weather forecasting chain (Pomoell & Poedts 2018). It covers the radial distances starting from 0.1 AU to 2 AU, including the orbit of Mars. It uses the finite volume TVDLF solver and a secondorder slope limiter. The details of the numerical methodology are given in Baratashvili et al. (2023). Icarus is a highly optimised heliospheric tool, as it uses advanced techniques, such as grid stretching and adaptive mesh refinement (AMR), to obtain results quickly at the desired resolution. Different AMR techniques were discussed in Verbeke et al. (2022) and Baratashvili et al. (2022), as the user has the freedom to tailor the AMR condition according to the purpose of the run. Icarus uses the output of the current semiempirical coronal model as inner boundary conditions for driving the solar wind at 0.1 AU. The reference frame is corotating with the Sun; thus, after the relaxation period passes (the time that is required for the slow solar wind to traverse the heliosphere from 0.1 AU to the outer boundary at 2 AU), the solar wind is stationary. After this, the CMEs can be injected from the inner heliospheric boundary. Currently, the simple hydrodynamics cone CME model and the magnetised linear forcefree spheromak models are supported (Baratashvili & Poedts 2024). Thus, Icarus can be used to study the propagation of the CMEs, its interaction with the solar wind and the effect on different satellites. Planets and satellites are included in the heliosphere, where plasma and magnetic field characteristic variables are sampled. Currently, there is Mercury, Venus, Mars, Earth, Stereo A, Stereo B, Solar Orbiter, and Parker Solar Probe sampling included. After performing the simulation, the observed data from the satellites can be compared to the sampled data at the satellite locations.
The coupling of COCONUT and Icarus was not straightforward, as COCONUT uses an unstructured grid and Icarus uses a structured grid. Data were extracted from the converged solution in COCONUT. It was interpolated, and data corresponding to 0.1 AU was extracted. Then, the data were stored similarly as the semiempirical WSA coronal model output, used both by the original EUHFORIA heliosphere and Icarus.
Figure 11 shows the input boundaries generated from the WSA semiempirical coronal model of EUHFORIA, the full MHD COCONUT coronal model with heating from Eq. (8), and the full MHD COCONUT coronal model with heating from Eq. (9). The variables plotted from top to bottom are radial velocity, number density, temperature, and radial magnetic field. The horizontal axis shows longitudes, and the vertical axis shows latitudes in degrees. The values between (30°, 150°) are shown on the vertical axis because the heliospheric models (i.e. the original EUHFORIA heliosphere and Icarus) do not have poles and extend from –60° to 60° in latitudes. In the simulation, the latitudes are transformed to colatitudes, where the equator corresponds to 90° degrees. The results from the WSA model and COCONUT models are quite different; however, the results from two different COCONUT simulations are indeed similar. The density and temperature are strongly overestimated in COCONUT simulations.
Fig. 11. Input boundary files generated from the WSA coronal model (left panel), COCONUT simulation with heating Eq. (8) (middle panel), and COCONUT simulation with heating Eq. (9) (right panel). The horizontal axis shows longitudes, and the vertical axis shows latitudes. The variables are given from top to bottom as follows: radial velocity in m s^{−1}, number density in m^{−3}, temperature in K, and radial magnetic field in T. 
Figure 12 shows the relaxed solutions for the Icarus simulations with the input boundary files demonstrated in Fig. 11. The radial velocity values are plotted in the equatorial plane. The left figure corresponds to the wind modelled by the WSA coronal input file, where the speed ranges between 250–700 km s^{−1}. There are multiple higher and lowerspeed streams. The middle figure corresponds to the wind simulated by the COCONUT model with the heating function in Eq. (8). The speed ranges between 350–700 km s^{−1}. The higher speed streams are more diffused into the lowspeed streams here. The last figure corresponds to the wind modelled by the COCONUT coronal model with the heating function in Eq. (9). The speed is considerably higher here, ranging between 550–1000 km s^{−1}. The time series obtained at Earth by these three different coronal models are demonstrated in Figs. 13 and 14. The former shows the radial velocity and number density values at Earth, whereas the latter shows the magnetic field components. In both figures, the results from the inputs of the WSA coronal model, COCONUT with heating (Eqs. (8) and (9)) are shown in red, orange and green colours, respectively. The black line corresponds to OMNI 1min data. The WSA input boundary slightly underestimates the speed values compared to observations and models the number density well, compared to the overall profile in the observed data; however, the peak in the observed number density profile is missing from the synthetic data. The COCONUT input boundary with the heating profile presented in Eq. (8) models the speed range well compared to the observations; however, it overestimates the overall density profile, although the peak density observed in the density can be noticed in the modelled data. The COCONUT input boundary with the heating profile presented in Eq. (9) strongly overestimates the speed values, but it does model the number density better with respect to the observed data than the COCONUT model with the heating profile from Eq. (8). The magnetic field values modelled by COCONUT simulations are quite similar. They both are different from the one generated with the WSA input boundary file. In the total magnetic field panel, we can see that COCONUT underestimates the total magnetic field strength compared to the WSA model and the observed data. Table 2 shows the times required for the introduced Icarus simulations with different input files. All of them are computed within 6 minutes on four nodes with 36 cores on the Genius cluster of the supercomputing center at KU Leuven.
Fig. 12. Equatorial plane, showing radial velocity values in Icarus simulations using different coronal boundary files from left to right: from the WSA model, from the COCONUT simulation with the heating function defined in Eqs. (8) and (9). Radial velocity values are given in [km s^{−1}]. 
Fig. 13. Time series obtained at Earth. The radial velocity [m s^{−1}] and number density [m^{−3}] values are plotted. The black curve corresponds to the OMNI 1min data. Green, red, and orange curves correspond to modelled data from Icarus with the following input boundary files: COCONUT Eq. (9), WSA model, and COCONUT (Eq. (8)). 
Fig. 14. Time series obtained at Earth. The magnetic field components are plotted in T. The black curve corresponds to the OMNI 1min data. Green, red, and orange curves correspond to modelled data from Icarus with the following input boundary files: COCONUT Eq. (9), WSA model, and COCONUT (Eq. (8)). 
Run times (wallclock time) for the heliospheric simulation for 24 days of simulation time in Icarus.
5. Conclusions and outlook
Our novel full MHD modelling chain was established from the Sun to Earth with COCONUT coronal and Icarus heliospheric models. To use COCONUT for heliospheric modelling, the coronal model was upgraded from the ideal polytropic MHD model to a full MHD model by introducing source and sink terms in the energy equation to account for thermal conduction, radiative cooling and coronal heating. The adiabatic index γ = 5/3 was used instead of γ = 1.05. In the first attempt, the coronal heating was approximated by various functions, depending on the radial distance and the magnetic field strength in the corona. The implemented full MHD coronal model was compared to the coronal model of MAS. This was done by comparing the bimodal structure of the wind in the meridional plane and the density profiles at 5 R_{⊙}. The results of the COCONUT simulations are also compared to tomography data. Different heating profiles generated different coronal configurations. The uniform spherically symmetric heating (introduced via Eq. (7)) did not obtain a realistic bimodal solar wind configuration; however, the obtained solar wind at the outer boundary was fast. The heating profiles, including the magnetic field dependence (i.e. Eqs. (8) and (9)) produced a bimodal solar wind structure. When introducing heating via Eq. (9), the simulated density field was more similar to the one obtained in the MAS model. The next validation mechanism is to connect COCONUT output to the heliosphere and see how these effects propagate to Earth.
The bimodal wind obtained near the outer coronal boundary is used as the onset of the heliospheric model in Icarus. COCONUT and Icarus are coupled by interpolating the plasma variables at 21.5 R_{⊙} from the unstructured grid in COCONUT to the uniform grid used in Icarus. The heliospheric simulation in Icarus is initiated with the coronal boundary file from COCONUT. The initial conditions were relaxed to obtain the steady wind in the Icarus domain that stretched from 21.5 R_{⊙} to 432 R_{⊙}. Furthermore, the time series were extracted at Earth and compared to OMNI 1min data in order to assess the modelling of the various plasma conditions near Earth. The modelled heliosphere is more consistent with the OMNI data, with the heating introduced in Eq. (8) compared to the other heating profile. The radial velocity values were more similar in the case of the heating introduced via Eq. (8); however, the number density was still strongly overestimated. The peak that was present in the number density data before the arrival of the higherspeed stream is better estimated by the same heating profile simulation. The difference in the magnetic field values is small at 1 AU and no strong conclusion can be made on this basis.
A comprehensive examination of the modelled data allowed us to assess the implemented full MHD model and identify the strengths and weaknesses of the model. Different heating approximations have show that it is important to deposit the thermal energy in the correct place to obtain a realistic bimodal wind. Prescribing the uniform heating with exponential damping solely resulted in accelerated wind everywhere near ∼21.5 R_{⊙}, but not in the bimodal structure. Accumulating the thermal energy in the strong magnetic field regions resulted in a more realistic image of the solar wind. Both functions introduced in Eqs. (8) and (9) produced bimodal wind near the outer boundary of the corona. This is the first achievement in the project since, as before, we could only get uniform solar wind distribution. The next question still remains regarding which solar wind is more realistic. The answer is important from the modelling point of view, with the goal to obtain as accurate results as possible, but also from the physics point of view. The main difference between these two approximated functions is that the first one treats the magnetic field only depending on its strength and introduces heating energy depending on the magnitude of the field, whereas the second one takes into account the contribution from the quiet Sun heating and active region heating separately and approximates both regimes with different formulas. In the case of the eclipse of 2019 (considered in this paper), the heating function approximated by Lionello et al. (2009) produced more realistic results when compared to the tomography data.
After comparing the different approximated heating profiles in the solar corona, we can see that we can get our first MHD results comparable to observations and the MAS model; however, crucial physics is still missing in our full MHD model. We intend to experiment with more physicsbased heating profiles, including the gradient of density and magnetic field. We also intend to implement the wave turbulencedriven heating mechanism to model Alfvén waves and the heating associated with them. As a first step, we obtained the full MHD chain from the Sun to Earth that is efficient and capable of operating for space weather purposes. As demonstrated, the advanced techniques, such as a fully implicit solver and an unstructured grid, allow us to include complex physics in our equations without significantly slowing down the code. Furthermore, the propagation of the flux rope CMEs in the full MHD model ought to be considered in order to get a more realistic thermodynamic evolution of the flux rope compared to the polytropic model (Linan et al. 2023; Guo et al. 2024). Finally, we intend to extend our model and include more realistic physics phenomena in subsequent works.
Acknowledgments
TB acknowledges the help and advice received from Jon linker and Cooper Downs while developing the full 3D MHD coronal model in COCONUT. This research has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 870405 (EUHFORIA 2.0) and the ESA project “Heliospheric modelling techniques” (Contract No. 4000133080/20/NL/CRS). These results were also obtained in the framework of the projects C14/19/089 (C1 project Internal Funds KU Leuven), AFOSR FA95501810093, G.0B58.23N and G.0025.23N (FWOVlaanderen), SIDC Data Exploitation (ESA Prodex12), and Belspo project B2/191/P1/SWiM. The Computational resources and services used in this work were provided by the VSCFlemish Supercomputer Center, funded by the Research Foundation Flanders (FWO) and the Flemish GovernmentDepartment EWI.
References
 Alves, M. V., Echer, E., & Gonzalez, W. D. 2006, J. Geophys. Res.: Space Phys., 111, A07S05 [NASA ADS] [CrossRef] [Google Scholar]
 Arge, C. N., Luhmann, J. G., Odstrcil, D., Schrijver, C. J., & Li, Y. 2004, J. Atmos. Sol.Terr. Phys., 66, 1295 [Google Scholar]
 Aschwanden, M. J. 2005, Physics of the Solar Corona. An Introduction with Problems and Solutions, 2nd edn. (New York, Berlin: Springer), 892 [Google Scholar]
 Baratashvili, T., & Poedts, S. 2024, A&A, 683, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baratashvili, T., Verbeke, C., Wijsen, N., & Poedts, S. 2022, A&A, 667, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baratashvili, T., Verbeke, C., Keppens, R., & Poedts, S. 2023, Sun Geosph., 15/2, 49 [NASA ADS] [Google Scholar]
 Brchnelova, M., Kuźma, B., Perri, B., et al. 2022a, ApJS, 263, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Brchnelova, M., Zhang, F., Leitner, P., et al. 2022b, J. Plasma Phys., 88, 905880205 [NASA ADS] [CrossRef] [Google Scholar]
 Downs, C., Roussev, I. I., van der Holst, B., et al. 2010, ApJ, 712, 1219 [CrossRef] [Google Scholar]
 Downs, C., Linker, J. A., Mikić, Z., et al. 2013, Science, 340, 1196 [NASA ADS] [CrossRef] [Google Scholar]
 Fisher, G. H., Longcope, D. W., Metcalf, T. R., & Pevtsov, A. A. 1998, ApJ, 508, 885 [CrossRef] [Google Scholar]
 Guo, J. H., Linan, L., Poedts, S., et al. 2024, A&A, 683, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hollweg, J. V. 1978, Rev. Geophys. Space Phys., 16, 689 [CrossRef] [Google Scholar]
 Kimpe, D., Lani, A., Quintino, T., Poedts, S., & Vandewalle, S. 2005, in Recent Advances in Parallel Virtual Machine and Message Passing Interface, eds. B. Di Martino, D. Kranzlmüller, & J. Dongarra (Berlin, Heidelberg: Springer, Berlin Heidelberg), 520 [Google Scholar]
 Kuźma, B., Brchnelova, M., Perri, B., et al. 2023, ApJ, 942, 31 [CrossRef] [Google Scholar]
 Lani, A., Quintino, T., Kimpe, D., et al. 2005, in Computational Science  ICCS 2005, eds. V. S. Sunderam, G. D. van Albada, P. M. A. Sloot, J. J. Dongarra, (Berlin, Heidelberg: Springer, Berlin Heidelberg), 279 [CrossRef] [Google Scholar]
 Lani, A., Villedieu, N., Bensassi, K., et al. 2013, AIAA 20132589, 21th AIAA CFD Conference, San Diego (CA) [Google Scholar]
 Lani, A., Yalim, M. S., & Poedts, S. 2014, Comput. Phys. Commun., 185, 2538 [NASA ADS] [CrossRef] [Google Scholar]
 Linan, L., Regnault, F., Perri, B., et al. 2023, A&A, 675, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lionello, R., Mikić, Z., & Linker, J. A. 1999, J. Comput. Phys., 152, 346 [NASA ADS] [CrossRef] [Google Scholar]
 Lionello, R., Linker, J. A., & Mikić, Z. 2001, ApJ, 546, 542 [Google Scholar]
 Lionello, R., Linker, J. A., & Mikić, Z. 2009, ApJ, 690, 902 [CrossRef] [Google Scholar]
 Lionello, R., Velli, M., Downs, C., et al. 2014, ApJ, 784, 120 [Google Scholar]
 Mikić, Z., & Linker, J. A. 1996, AIP Conf. Proc., 382, 104 [CrossRef] [Google Scholar]
 Mikić, Z., Linker, J. A., Schnack, D. D., Lionello, R., & Tarditi, A. 1999, Phys. Plasmas, 6, 2217 [Google Scholar]
 Mikić, Z., Downs, C., Linker, J. A., et al. 2018, Nat. Astron., 2, 913 [Google Scholar]
 Morgan, H. 2015, ApJS, 219, 23 [CrossRef] [Google Scholar]
 Morgan, H. 2019, ApJS, 242, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Morgan, H., & Cook, A. C. 2020, ApJ, 893, 57 [CrossRef] [Google Scholar]
 Moses, J. D., Antonucci, E., Newmark, J., et al. 2020, Nat. Astron., 4, 1134 [Google Scholar]
 Odstrčil, D., Riley, P., & Zhao, X. P. 2004, J. Geophys. Res., 109, 2116 [Google Scholar]
 Parenti, S., Réville, V., Brun, A. S., et al. 2022, ApJ, 929, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Park, S.H., Cho, K.S., Bong, S.C., et al. 2012, ApJ, 750, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Perri, B., Leitner, P., Brchnelova, M., et al. 2022, ApJ, 936, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Perri, B., Kuźma, B., Brchnelova, M., et al. 2023, ApJ, 943, 124 [NASA ADS] [CrossRef] [Google Scholar]
 Pevtsov, A. A., Fisher, G. H., Acton, L. W., et al. 2003, ApJ, 598, 1387 [Google Scholar]
 Poirier, N., Rouillard, A. P., Kouloumvakos, A., et al. 2021, Front. Astron. Space Sci., 8, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Pomoell, J., & Poedts, S. 2018, J. Space Weather Space Clim., 8, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Regnault, F., Strugarek, A., Janvier, M., et al. 2023, A&A, 670, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Réville, V., Folsom, C. P., Strugarek, A., & Brun, A. S. 2016, ApJ, 832, 145 [CrossRef] [Google Scholar]
 Réville, V., Velli, M., Panasenco, O., et al. 2020, ApJS, 246, 24 [Google Scholar]
 Réville, V., Poirier, N., Kouloumvakos, A., et al. 2023, J. Space Weather Space Clim., 13, 11 [CrossRef] [EDP Sciences] [Google Scholar]
 Riley, P., Linker, J. A., Lionello, R., & Mikic, Z. 2012, JASTP, 83, 1 [NASA ADS] [Google Scholar]
 Rosner, R., Tucker, W. H., & Vaiana, G. S. 1978, ApJ, 220, 643 [Google Scholar]
 Roussev, I. I., & Sokolov, I. V. 2006, Geophys. Monogr. Ser., 165, 89 [NASA ADS] [Google Scholar]
 Schrijver, C. J., Zwaan, C., Maxson, C. W., & Noyes, R. W. 1985, A&A, 149, 123 [NASA ADS] [Google Scholar]
 Shiota, D., & Kataoka, R. 2016, Space Weather, 14, 56 [CrossRef] [Google Scholar]
 van der Holst, B., Sokolov, I. V., Meng, X., et al. 2014, ApJ, 782, 81 [Google Scholar]
 Verbeke, C., Baratashvili, T., & Poedts, S. 2022, A&A, 662, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wang, Y. M., & Sheeley, N. R., Jr. 1990, ApJ, 355, 726 [NASA ADS] [CrossRef] [Google Scholar]
 Xia, C., Teunissen, J., El Mellah, I., Chané, E., & Keppens, R. 2018, ApJS, 234, 30 [Google Scholar]
 Yogesh, Chakrabarty D., & Srivastava, N., 2021, MNRAS, 503, L17 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Run times (wallclock time) required for the COCONUT simulations to reach convergence –3 in density.
Run times (wallclock time) for the heliospheric simulation for 24 days of simulation time in Icarus.
All Figures
Fig. 1. Synoptic maps for July 2, 2019, Carrington Rotation 2219. The original (left) and pole filled (right) HMI magnetograms are provided in the upper panel. The bottom panel shows the processed corresponding input magnetogram with spherical harmonics decomposition filtered with l_{max} = 20. The units are in Gauss. In the processed magnetogram, the values are given in the code values (divided by 2.2 G). 

In the text 
Fig. 2. Polytropic solution for the CR2219 corresponding to July 2, 2019. The radial velocity is plotted in the meridional plane. 

In the text 
Fig. 3. Polytropic solution for the CR2219 corresponding to July 2, 2019. The eclipse image is overlaid with the magnetic streamlines from the simulation. The surface of the Sun is coloured with the B_{r} in Gauss. 

In the text 
Fig. 4. Horizontal axis shows longitudes in degrees and the vertical axis shows colatitudes in degrees. The left and right panels show the B_{r} and density in normalised units from the COCONUT and MAS models, respectively. The bottom row shows the observed tomography data in normalised units. All quantities are plotted at 5 R_{⊙}. Data from COCONUT are taken from the polytropic simulation. 

In the text 
Fig. 5. Radiative loss cooling curve defined as in Rosner et al. (1978). The horizontal axis denotes temperature and the vertical axis shows the radiative loss cooling curve function in cgs. Loglog scale is applied. 

In the text 
Fig. 6. Obtained heating profiles in the meridional plane for the full MHD COCONUT simulations for each heating description given in Eqs. (7)–(9), from left to right, respectively, in [J m^{−3} s^{−1}]. 

In the text 
Fig. 7. Radial velocity values in [km s^{−1}] in the meridional plane for the full MHD COCONUT simulations with the indicated heating functions given in Eqs. (7)–(9), from left to right, respectively. 

In the text 
Fig. 8. Full MHD COCONUT solutions with the indicated heating profiles given in Eqs. (7)–(9), from left to right, respectively. The eclipse image is overlaid with the magnetic streamlines from the COCONUT solution. The radial magnetic field is plotted on the solar surface. 

In the text 
Fig. 9. Longitudes in degrees shown on the horizontal axis. The vertical axis gives the colatitudes in degrees. Here, B_{r} is given in [nT], the density from COCONUT in normalised units, and tomography in normalised units, from top to bottom, respectively. All quantities are plotted at 5 R_{⊙}. The left panel shows results for the COCONUT simulation with the heating function from Eq. (7), the middle panel shows the results with the heating from Eq. (8), and the right includes the heating function from Eq. (9). 

In the text 
Fig. 10. Longitudes in degrees are shown on the horizontal axis. Vertical axis gives the colatitudes in degrees. The panels represent Br in [nT], the density from MAS simulations in normalised units, and tomography in normalised units from top to bottom, respectively. All quantities are plotted at 5 R_{⊙}. The MAS simulation was obtained with thermodynamic model 1. 

In the text 
Fig. 11. Input boundary files generated from the WSA coronal model (left panel), COCONUT simulation with heating Eq. (8) (middle panel), and COCONUT simulation with heating Eq. (9) (right panel). The horizontal axis shows longitudes, and the vertical axis shows latitudes. The variables are given from top to bottom as follows: radial velocity in m s^{−1}, number density in m^{−3}, temperature in K, and radial magnetic field in T. 

In the text 
Fig. 12. Equatorial plane, showing radial velocity values in Icarus simulations using different coronal boundary files from left to right: from the WSA model, from the COCONUT simulation with the heating function defined in Eqs. (8) and (9). Radial velocity values are given in [km s^{−1}]. 

In the text 
Fig. 13. Time series obtained at Earth. The radial velocity [m s^{−1}] and number density [m^{−3}] values are plotted. The black curve corresponds to the OMNI 1min data. Green, red, and orange curves correspond to modelled data from Icarus with the following input boundary files: COCONUT Eq. (9), WSA model, and COCONUT (Eq. (8)). 

In the text 
Fig. 14. Time series obtained at Earth. The magnetic field components are plotted in T. The black curve corresponds to the OMNI 1min data. Green, red, and orange curves correspond to modelled data from Icarus with the following input boundary files: COCONUT Eq. (9), WSA model, and COCONUT (Eq. (8)). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.