© The Authors 2025

1. Introduction

Space weather refers to the variable physical conditions on the Sun and in the solar wind, the magnetosphere, and the ionosphere. It can influence the performance and reliability of space- and ground-based technological systems and affect human life and health. There is thus a need to develop advanced Sun-to-Earth model chains to understand the mechanisms of space weather and ultimately provide reliable space weather forecasts hours to days in advance (e.g. Baker 1998; Feng et al. 2011a, 2013; Feng 2020; Koskinen et al. 2017). Physics-based magnetohydrodynamic (MHD) modelling is the first principal method capable of self-consistently bridging large heliocentric distances from near the Sun to well beyond Earth’s orbit (e.g. Detman et al. 2006; Feng et al. 2007, 2014b, 2010, 2011b, 2012a,b, 2014a, 2017; Feng 2020; Gombosi et al. 2018; Hayashi et al. 2006; Li & Feng 2018; Li et al. 2020; Lugaz & Roussev 2011; Mikić et al. 1999; Nakamizo et al. 2009; Riley et al. 2012; Shen et al. 2021; Tóth et al. 2012; Usmanov 1993; Usmanov & Goldstein 2003; Wu & Dryer 2015; Yang et al. 2021; Zhou et al. 2012; Zhou & Feng 2017). However, realistic MHD simulations of the solar-terrestrial system are complex, involving various physical phenomena across diverse spatiotemporal scales, and are very computationally intensive. We need to develop more efficient and reliable MHD models to improve space weather forecasting (e.g. Feng 2020; Owens et al. 2017, and references therein).

Since solar-terrestrial space involves diverse spatiotemporal scales and phenomena, using a single model to simulate the entire solar-terrestrial space is inefficient. Coupling different models, respectively dedicated to specific regions and physical problems, has become a preferred approach for establishing a space weather forecasting framework (e.g. Feng et al. 2013; Goodrich et al. 2004; Hayashi et al. 2021; Odstrcil et al. 2004; Pomoell & Poedts 2018; Poedts et al. 2020; Tóth et al. 2012). In a coupled Sun-to-Earth modelling chain, observed photospheric magnetic fields serve as input data for the solar corona model. The solar corona model provides the inner heliospheric model’s inner boundary conditions. The inner heliospheric model gives boundary information to the geomagnetic model. The solar corona model is crucial for determining the initialisation of the modelling chain. Thus, it is also a key factor affecting the simulation results of solar disturbance propagation and evolution (Brchnelova et al. 2022a; Perri et al. 2023). Specifically, the solar wind speed increases from subsonic and/or sub-Alfvénic to supersonic and/or super-Alfvénic in the solar corona, and solar disturbances such as coronal mass ejections (CMEs) and solar proton events also propagate through this layer (Feng 2020; Kuźma et al. 2023).

Though the solar corona is a crucial link in the Sun-to-Earth modelling chain and significantly impacts the ultimate effects of space weather, physics-based MHD corona models are also the most complex and computationally intensive component. Specifically, depending on the mesh resolution, even the currently state-of-the-art quasi-steady-state corona simulations take 10 ∼ 100 k CPU-hours to reach a quasi-steady state (Feng et al. 2019; Réville et al. 2020). In these simulations, time steps are limited to a few seconds due to the restriction of the Courant-Friedrichs-Lewy (CFL) stability condition. In contrast, the time steps limited by the CFL condition are typically of the order of 10 minutes for most MHD inner heliospheric models (Detman et al. 2006; Hayashi 2012). As a result, MHD corona models require significantly more computational resources to remain synchronised with inner heliospheric models. Significant simplifications in solar corona modelling are frequently needed to improve efficiency. However, empirical solar corona models (Arge et al. 2003; Yang et al. 2018) discard important information, and it has been demonstrated that even a relatively simple MHD model provides better forecasts (Samara et al. 2021). Therefore, more efforts are required to establish more efficient and accurate MHD solar corona models.

Sokolov et al. (2021) solved 1D equations for the plasma motion and heat transport along fixed potential magnetic field lines together with the Alfvénic wave propagation in the low coronal region and interfaced this threaded-field-line model with the full MHD global corona model at 1.1 R_s. This enabled the updated Alfvén wave solar atmosphere model (AWSoM) model, AWSoM-R (AWSoM-realtime), to achieve a faster-than-real-time performance on approximately 200 CPU cores. In contrast, the original AWSoM model required 1000–2000 CPU cores to reach a similar computing time (Jin et al. 2017). More generally, implicit temporal discretisation strategies can be used to overcome the limitation imposed by the CFL condition, thereby improving the overall computational efficiency by utilising a larger time step.

Recently, several successful attempts have been made to increase the efficiency of MHD corona models by using implicit solvers. For instance, Wang et al. (2019b) achieved speedup ratios of 31.27 and 28.05 in MHD solar corona simulations with an effective matrix-free implicit scheme. Feng et al. (2021) and Wang et al. (2022a,b) further improved the implicit method and developed an efficient parallel lower-upper symmetric Gauss-Seidel solver, improving computational efficiency even with a plasma β smaller than 10⁻⁴. COolfluid COroNal UnsTructured (COCONUT), a novel MHD solar corona model based on the Computational Object-Oriented Libraries for Fluid Dynamics (COOLFluiD)¹, gained a speed up of 35 compared to the state-of-the-art time-explicit wind-predict model (Perri et al. 2018) while achieving the same level of accuracy in steady-state simulations (Perri et al. 2022, 2023). In addition, CME simulations (Guo et al. 2023; Linan et al. 2023) demonstrate that the COCONUT model has the potential to be faster than the explicit MHD SC models in time-dependent simulations while maintaining time-accurate using a relatively large time step. Wang et al. (2025) proposed an efficient time-accurate implicit MHD model for the solar corona and CMEs, which is capable of timely and accurately simulating time-varying events in the solar corona, even with low plasma beta. The fully implicit scheme of the block-adaptive tree solarwind Roe-type upwind scheme (BATS-R-US) code (Tóth et al. 2012), the explicit scheme of which was used by AWSoM corona model, can also produce speedup ratios of the order of 10–20 compared to the explicit version in simulations of some geophysical applications (Tóth et al. 2006; Tóth & Ma 2008).

However, it is still necessary to further develop coronal modelling to capture time-evolving solar coronal structures more accurately, as current models frequently assume a (quasi-)steady corona during one Carrington rotation (CR) and use a time-invariant photospheric magnetic field as the inner boundary. The (quasi-)steady simplification differs from the reality that the solar coronal structure evolves (Owens et al. 2017) and thus leads to discrepancies between the simulation results and observations of the solar coronal structures (Cash et al. 2015; Réville et al. 2020). Time-evolving corona models typically driven by hundreds of time-evolving observed photospheric magnetograms can capture time-evolving coronal structures with higher fidelity (Feng et al. 2023; Yang et al. 2012) and thus may also improve solar wind and CME modelling (Lionello et al. 2023).

For instance, Detman et al. (2006), Linker et al. (2016), and Merkin et al. (2016) used a time series of photospheric magnetograms to drive an empirical source surface current sheet corona model to provide time-evolving lower boundary conditions for the inner heliosphere MHD solar wind model. These hybrid Sun-to-Earth modelling systems are efficient for real-time operations but fail to produce coronal transient events or model closed loops that rise beyond the source surface (usually between 2.0 and 2.5 R_s). Hoeksema et al. (2020) and Hayashi et al. (2021) employed the electric field derived at the Sun’s photosphere from a sequence of vector magnetogram and Doppler velocity measurements (Fisher et al. 2020) to drive a magneto-frictional (MF) model, and then used the MF model to produce time-evolving boundary magnetic fields at 1.15 R_s for their global coronal-heliospheric MHD (GHM) model. Although the MF model can be numerically stable in generating time-dependent 3D coronal magnetic structures, it failed to provide the required initial dynamic states of the plasma in the low atmosphere. To ensure the boundary plasma quantities at 1.15 R_s evolve consistently with both the variations of the magnetic field specified from the MF model and the governing MHD equations of the GHM model, they employed the projected normal characteristics method (e.g. Sauerwein 1966) at the interface boundary. Feng et al. (2023) and Yang et al. (2012) specified the tangential boundary electric field to make the flux evolution match the changes of the observed radial magnetic field and employed the projected normal characteristic method to make boundary conditions self-consistent. Feng et al. (2012a); Lionello et al. (2023), and Mason et al. (2023) employed the surface flux transport model (DeVore et al. 1984; Schrijver & DeRosa 2003) to obtain input maps that incorporate magnetic flux emergence and surface flows for their MHD corona models.

Another time-evolving corona model, magnetohydrodynamic algorithm outside a sphere (MAS; Lionello et al. 2023; Mason et al. 2023), adopts a semi-implicit approach in which only some source terms are treated implicitly, requiring the time step to be selected according to the limitations imposed by the explicitly treated terms. Feng et al. (2023) used the fully implicit method (Feng et al. 2021; Wang et al. 2022a) only at the six-layer grid close to the solar surface in the radial direction and reported that it helped avoid the occurrence of negative density and pressure caused by the strong magnetic field near the Sun. The remaining models still use explicit methods.

We have extended the quasi-steady-state COCONUT to a time-evolving corona model and uses cubic Hermite interpolation to derive the time-evolving magnetograms at each physical time step. This extension demonstrates that the time-evolving corona model can achieve significantly better efficiency and numerical stability by adopting fully implicit algorithms appropriately.

For robust solar corona modelling, ensuring the positivity-preserving (PP) property of thermal pressure and density in MHD simulations is vital. For example, a self-adjusting PP reconstruction method, initially proposed by Balsara (2012) for solving hydrodynamic and MHD equations, was implemented in solar corona simulations by Feng et al. (2017). It was applied to conservative variables via a flattener function defined according to the flow’s rarefactive and compressive motions. Furthermore, Feng et al. (2021) and Wang et al. (2022a) extended this method to primitive variables and implemented it in implicit MHD corona models. Additionally, some PP Harten-Lax-van Leer (HLL) Riemann solvers (Wu & Shu 2019) have been used to design PP MHD corona models (Feng et al. 2021; Wang et al. 2022a). It has been shown that the approximate reconstruction method and flux solver are beneficial to maintaining the PP property of MHD models. In this study, under the cell-centred finite-volume framework of COCONUT, we implemented the Venkatakrishnan limiter (Venkatakrishnan 1993) in the reconstruction formula of primitive variables as before to control spatial oscillation. Furthermore, we designed an HLL Riemann solver with a self-adjustable dissipation term to accommodate low- and high-speed flows.

Brchnelova et al. (2023) improved the PP property of COCONUT by manipulating the inner-boundary density according to the local Alfvénic velocity. In this method, the active region density was carefully increased to avoid an abnormally high Alfvénic speed when the local inner-boundary Alfvénic velocity exceeds a prescribed maximum Alfvénic speed. This method can prevent some non-physical negative thermal pressures and very high-speed streams developed in the domain above the active regions (Kuźma et al. 2023). Consequently, the performance, in terms of both convergence and physical accuracy, can be improved. To enhance the PP property of this time-evolving MHD corona model, we further extended this method to all the computational domains and applied a smooth hyperbolic tangent function to gradually increase the plasma density when the local Alfvénic speed exceeds 2000 km/s. In the future, we can take some extra measures, such as applying the self-adjusting PP reconstruction method (Feng et al. 2021; Wang et al. 2022a) and incorporating a mass flux limitation (Hayashi 2005; Yang et al. 2012), to enhance the PP property of COCONUT further.

Moreover, reducing the magnetic field discretisation error is also helpful in maintaining the PP property of MHD models in low β regions. Considering that (B+ϵB)² − B² ≡ 2 ϵB² + ϵ²B², with ϵB denoting the magnetic field discretisation error, the magnetic pressure discretisation error can be comparable to thermal pressure in low β (the ratio of the thermal pressure to the magnetic pressure) regions, and non-physical negative thermal pressure is prone to appear when deriving thermal pressure from energy density. Some researchers try to decrease the magnetic field discretisation error by adopting fine meshes near the solar surface to avoid such an undesirable situation. In AWSoM, a spherical grid ranging from 1 R_s to 24 R_s is used for the coronal component. The grid is highly refined in the radial direction near the Sun, with the smallest radial grid spacing being approximately 700 km and the total number of cells being 29.7 M in van der Holst et al. (2022). For the MAS model, a non-uniform spherical grid consisting of 27.3 M cells and covering a radial range from 1 R_s to 30 R_s is used in Caplan et al. (2017). It is highly non-uniform in the radial direction, and the smallest cell size in the radial direction is 340 km. This paper adopts the unstructured geodesic mesh (Brchnelova et al. 2022b; Perri et al. 2022), which consists of approximately 1.5 M cells. The grid covers a radial range from 1 R_s to 25 R_s and gradually stretches in the radial direction with an initial grid spacing of about 170 km at the solar surface.

The paper is organised as follows. In Sect. 2 we introduce the governing equations, grid system, and initial conditions adopted in the solar corona simulations. In Sect. 3 and Appendix A, the numerical formulation of the time-evolving corona simulations is described in detail. These sections mainly describe the discretisation of the MHD equations, the implementation of time-evolving boundary conditions, the HLL Riemann solver with a self-adjustable dissipation term, and the PP measures used to enhance the corona model’s numerical stability. We demonstrate the simulation results in Sect. 4: the evolution of the corona during two CR periods around the 2019 eclipse simulated by the time-evolving COCONUT is demonstrated, a comparison of the 2019 eclipse simulations performed by both the time-evolving COCONUT and its quasi-steady-state version is illustrated, and the simulation results calculated with different time-step sizes are compared. In Sect. 5 we summarise the main features of the efficient, fully implicit, time-evolving corona model and give some concluding remarks.

2. Governing equations and grid system

This section mainly describes the governing equations, grid system, and initial conditions for simulating quasi-steady-state corona and time-evolving corona simulations.

2.1. The governing equations

In this study we developed the time-evolving version of COCONUT for time-evolving corona simulations. First, we initiated a quasi-steady-state corona simulation constrained by a fixed magnetogram to get the background corona. Once the steady-state simulation converged, we drove the subsequent evolution of the dynamic corona by a series of time-evolving photospheric magnetograms.

The governing MHD equations are calculated in the heliocentric inertial (HCI) coordinate system (Burlaga 1984; Fränz & Harper 2002):

$$\begin{array}{c}\hfill \frac{\partial \mathit{U}}{\partial t}+\mathrm{\nabla }·\mathit{F}\left(\mathit{U}\right)=\mathit{S}(\mathit{U},\mathrm{\nabla }\mathit{U}).\end{array}$$(1)

Here U = (ρ,ρv,B,E,ψ)^T denotes the conservative variable vector, ∇U corresponds to the spatial derivative of U, the inviscid flux vector F(U) is

$$\begin{array}{c}\hfill \mathit{F}\left(\mathit{U}\right)=\left(\begin{array}{c}\rho \mathit{v}\\ \rho \mathit{v}\mathit{v}+p_{\mathrm{T}}\mathit{I}-\mathit{B}\mathit{B}\\ \mathit{vB}-\mathit{Bv}+\psi \mathit{I}\\ (E+p_{\mathrm{T}})\mathit{v}-\mathit{B}(\mathit{v}·\mathit{B})\\ V_{\mathrm{ref}}^2\mathit{B}\end{array}\right),\end{array}$$

and S(U,∇U) = S_gra + S_heat represents the source term vector corresponding to the gravitational force and the heating source terms defined as follows:

$$\begin{array}{c}\hfill {\mathit{S}}_{\mathrm{gra}}=-\frac{\rho G{M}_{\mathrm{s}}}{{\left|\mathit{r}\right|}^3}\left(\begin{array}{c}0\\ \mathit{r}\\ \mathbf{0}\\ \mathit{r}·\mathit{v}\\ 0\end{array}\right),{\mathit{S}}_{\mathrm{heat}}=\left(\begin{array}{c}0\\ \mathbf{0}\\ \mathbf{0}\\ -\mathrm{\nabla }·\mathit{q}+Q_{\mathrm{rad}}+Q_{\mathrm{H}}\\ 0\end{array}\right).\end{array}$$(2)

In these formulations mentioned above, B = (B_x,B_y,B_z) and v = (u,v,w) denote the magnetic field and velocity in Cartesian coordinate system, $E=\frac{p}{\gamma -1}+\frac{1}{2}\rho {\mathit{v}}^2+\frac{1}{2}{\mathit{B}}^2$ means the total energy density with the adiabatic index $\gamma =\frac{5}{3}$, $p_{\mathrm{T}}=p+\frac{{\mathit{B}}^2}{2}$ is the total pressure, ρ and p represent the density and thermal pressure of the plasma, r is the position vector, r = |r| denotes the heliocentric distance, and t represents the time. For convenience of description, in the definition of the magnetic field, a factor of $\frac{1}{\sqrt{{\mu }_0}}$ is absorbed with μ₀ = 4 × 10⁻⁷π H m⁻¹ denoting the magnetic permeability. As usual, G means the universal gravitational constant, M_s means the mass of the Sun, and GM_s = 1.327474512 × 10²⁰ m³ s⁻². The thermal pressure of the plasma is defined as p = ℜρT, where T is the temperature of the bulk plasma, ℜ = 1.299 × 10⁴ m² s⁻² K⁻¹ denotes the gas constant and is calculate by $\mathfrak{R}=\frac{2\ast k_{\mathrm{B}}}{m_{\mathrm{cor}}\ast m_{\mathrm{H}}}$ with k_B = 1.3806503 × 10⁻²³ J K⁻¹ denoting the Boltzmann constant, the molecular weight setting to m_cor = 1.27 (Aschwanden 2005) and m_H = 1.67262158 × 10⁻²⁷ kg representing the mass of hydrogen. We adopted the hyperbolic generalised Lagrange multiplier method (Dedner et al. 2002; Yalim et al. 2011) to constrain the divergence error; ψ and V_ref denote the Lagrange multiplier and the propagation speed of the numerical divergence error. In the energy source term, Q_H, Q_rad, and −∇⋅q are the coronal heating, radiation loss, and thermal conduction, respectively.

As done in Baratashvili et al. (2024), the heat flux, q, included in S_heat is defined in a Spitzer or collisionless form according to the radial distance (Hollweg 1978; Mikić et al. 1999):

$$\begin{array}{c}\hfill \mathit{q}=\{\begin{array}{c}-\xi T^{5/2}(\widehat{\mathit{b}}·\mathrm{\nabla }T)\widehat{\mathit{b}},\text{if}1\le r\le 10R_{\mathrm{s}}\\ \alpha n_{\mathrm{e}}{k}_{\mathrm{B}}T\mathit{v},\text{if}r>10R_{\mathrm{s}}\end{array}.\end{array}$$(3)

Here $\widehat{\mathit{b}}=\frac{\mathit{B}}{\left|\mathit{B}\right|}$, $\xi =9.0\times {10}^{-12}\mathrm{J}{\mathrm{m}}^{-1}{\mathrm{s}}^{-1}{\mathrm{K}}^{-\frac{7}{2}}$, α is set to $\frac{3}{2}$ (Lionello et al. 2008) and n_e is the electron number density. By assuming the radiative loss to be optically thin (Rosner et al. 1978; Zhou et al. 2021), the radiative term is

$$\begin{array}{c}\hfill Q_{\mathrm{rad}}=-n_{\mathrm{e}}{n}_{\mathrm{p}}\mathrm{\Lambda }\left(T\right),\end{array}$$(4)

where the proton number density, n_p, is assumed to be equal to the electric number density, n_e, for the hydrogen plasma. Λ(T) is a temperature-dependent radiative cooling curve function. Similar to van der Holst et al. (2014), Λ(T) in this paper is derived from version 9 of CHIANTI (Dere et al. 2019), an atomic database for emission lines. The profile of log(Λ(T)) along log(T), with Λ(T) and T in unite of erg s⁻¹ cm³ and K, is shown in Fig. 1. As did in Xia et al. (2011), we set Λ(T) to zero when T < 2 × 10⁴ K, which means the plasma has become optically thick and is no longer fully ionized. We adopted the following empirical heating source term discussed and recommended in Baratashvili et al. (2024). It is proportional to the local magnetic field strength (Mok et al. 2005) while also exhibiting an exponential decay in the radial direction (Downs et al. 2010):

$$\begin{array}{c}\hfill Q_{\mathrm{H}}=H_0·\left|\mathit{B}\right|·e^{-\frac{r-R_{\mathrm{s}}}{\lambda }},\end{array}$$(5)

Fig. 1. Radiative cooling curve profile derived from version 9 of the CHIANTI atomic database (Dere et al. 2019). The horizontal axis denotes the decadic logarithm of temperature, and the vertical axis the logarithm of the radiative cooling curve function value. The temperature and radiative cooling curve function units are K and erg s−1 cm3, respectively.

where H₀ = 4 × 10⁻² J m⁻³ s⁻¹ T⁻¹ and λ = 0.7R_s.

The variables r, ρ, v, p, B, and t are normalized by R_s, ρ_s, V_a, s, ρ_sV_a, s², B_s, and $\frac{R_{\mathrm{s}}}{V_{\mathrm{a},\mathrm{s}}}$, respectively. Here R_s = 6.95 × 10⁸ m is the solar radius and ρ_s = 1.67 × 10⁻¹³ kg m⁻³, B_s = 2.2 × 10⁻⁴ T, and $V_{\mathrm{a},\mathrm{s}}=\frac{B_{\mathrm{s}}}{{\rho }_{\mathrm{s}}^{0.5}}$ denote the characteristic plasma density, magnetic field strength, and Alfvénic velocity at the solar surface, respectively.

2.2. Computational grid system and initialisations

The governing equations are solved using the cell-centred finite-volume method. The computational domain is a spherical shell ranging from 1.01 to 25 R_s and discreted into an unstructured sixth-level subdivided geodesic mesh (Brchnelova et al. 2022b), which includes 1495040 non-overlapped pyramid cells, each consisting of 2 triangular faces and 3 quadrilateral faces. There are 73 layers of gradually stretched cells in the radial direction, each containing 20480 cells.

The observed GONG-zqs (zeropoint-corrected) hourly-updated synoptic maps of the photospheric magnetic field provide inner-boundary conditions of the magnetic field (Li et al. 2021; Perri et al. 2023). Considering that the original photospheric magnetic field strength can be hundreds of Gauss near active regions and may cause non-physical negative temperatures or pressures, and the magnetic maps reconstructed by a potential-field (PF) solver of 20-order spherical harmonic expansion are already fully sufficient for space weather forecasting (Kuźma et al. 2023), we used a PF solver that filters out spherical harmonic expansions above the 20th order to eliminate very small-scale structures from the original magnetograms. These preprocessed photospheric magnetograms constrain the corona simulations in COCONUT (Kuźma et al. 2023; Perri et al. 2022, 2023). Meanwhile, the temperature and plasma density at the solar surface are set to 1.5 × 10⁶ K and 1.67 × 10⁻¹³ kg m⁻³, respectively.

3. Numerical methods

This section introduces the discretisation of the MHD equations and the implementation of inner-boundary conditions.

3.1. Spatial discretisation and temporal integration

Godunov’s method is used in COCONUT to advance cell-averaged solutions in time by solving a Riemann problem at each cell interface (Godunov 1959). By integrating Eq. (1) over the prismatic cell cell_i and applying Gauss’s theorem to calculate the volume integral of the flux divergence, we obtained the following discretized equations:

$$\begin{array}{c}\hfill V_i\frac{\partial {\mathit{U}}_i}{\partial t}=-{\displaystyle {\oint }_{\partial V_i}\mathit{F}·\mathit{n}d\mathrm{\Gamma }+V_i{\mathit{S}}_i,}\end{array}$$(6)

where ${\oint }_{\partial V_i}\mathit{F}·\mathit{n}d\mathrm{\Gamma }={\displaystyle \sum _{j=1}^5}{\mathit{F}}_{\mathit{ij}}·{\mathit{n}}_{\mathit{ij}}{\mathrm{\Gamma }}_{\mathit{ij}}$ and S_i = S_gra, i + S_heat, i. U_i and S_i denote the cell-averaged solution variables and source terms in cell_i, V_i is the volume of cell_i, Γ_ij means the interface shared by cell_i and its neighbouring cell_j, and also denotes the area of this interface, n_ij = (n_x, ij,n_y, ij,n_z, ij) is the unit normal vector of Γ_ij, oriented from cell_i to cell_j, and F_ij ⋅ n_ij is described as F_n, ij, denoting the inviscid flux along normal direction of Γ_ij.

We employed an HLL Riemann solver (Einfeldt et al. 1991) with a self-adjustable dissipation term, as described in Appendix A, to compute the inviscid flux F_n, ij. The cell-averaged source terms S_gra, i and Q_rad, i and Q_H, i were calculated by substituting the corresponding variables at the centroid of cell_i into formulations of S_gra and Q_rad and Q_H. (∇⋅q)_i was calculated using the Green-Gauss theorem as follows:

$$\begin{array}{c}\hfill {(\mathrm{\nabla }·\mathit{q})}_i=\frac{1}{V_i}{\displaystyle \sum _{j=1}^5{\mathit{q}}_{\mathit{ij}}·{\mathit{n}}_{\mathit{ij}}{\mathrm{\Gamma }}_{\mathit{ij}},}\end{array}$$

where q_ij = q(T_ij,(∇T)_ij,U_ij) represent the heat flux through Γ_ij. In this formulation, T_ij is derived from the equation of state $T_{\mathit{ij}}=\frac{p_{\mathit{ij}}}{\mathfrak{R}{\rho }_{\mathit{ij}}}$. The plasma states at the cell surface Γ_ij are required for the calculation of the inviscid flux F_n, ij and heat flux q_ij ⋅ n_ij. We used a second-order accurate reconstruction method to calculate the piecewise linear polynomial of the primitive variables:

$$\begin{array}{c}\hfill P_i(\mathit{x})={P|}_i+{\varphi }_i{\left(\mathrm{\nabla }P\right)|}_i·(\mathit{x}-{\mathit{x}}_i),\end{array}$$(7)

where P ∈ {ρ, u, v, w, B_x, B_y, B_z, p, ψ}, P|_i is the corresponding primitive variables at the centroid of cell_ix_i, ${\left(\mathrm{\nabla }P\right)|}_i={(\frac{\partial P}{\partial x},\frac{\partial P}{\partial \mathit{y}},\frac{\partial P}{\partial z})|}_i$ is the derivative of P at x_i calculated by least square method (Barth 1993; Lani 2008), and ϕ_i is the Venkatakrishnan limiter (Venkatakrishnan 1993) used to control spatial oscillation.

In the quasi-steady corona simulations, we applied a backward Euler temporal integration to Eq. (6) and obtained the following equation:

$$\begin{array}{c}\hfill V_i\frac{\mathrm{\Delta }{\mathit{U}}_i^n}{\mathrm{\Delta }t}+{\mathit{R}}_i\left({\mathit{U}}^{n+1}\right)=\mathbf{0}.\end{array}$$(8)

The superscripts ‘ⁿ’ and ‘^n + 1’ denote the time level, ${\mathit{R}}_i\left({\mathit{U}}^{n+1}\right)=\sum _{j=1}^5{\mathit{F}}_{n,ij}({\mathit{U}}_L^{n+1},{\mathit{U}}_R^{n+1}){\mathrm{\Gamma }}_{\mathit{ij}}-V_i{\mathit{S}}_i^{n+1}$ means the residual operator over cell_i at the (n + 1)-th time level, ΔU_iⁿ = U_i^n + 1 − U_iⁿ represents the solution increment between the n-th and (n + 1)-th time levels, and Δt = t^n + 1 − tⁿ is a user-defined time increment. In the quasi-steady corona simulations, the time variable t does not refer to a physical time but a relaxation time used to get a quasi-steady coronal structure.

After reaching a steady-state coronal structure, we utilised a series of time-varying magnetograms to drive the following evolution of the dynamic corona. Along with spatial accuracy, temporal accuracy is also crucial in time-evolving coronal simulations. To enhance the temporal accuracy of the implicit solver with time step lengths exceeding the CFL condition, we adopted the second-order accurate backward differentiation formula (BDF2) for temporal integration. Newton iterations are used within each time step of the implicit algorithm to update the intermediate state until the L₂ norm of the differences in state variables between two consecutive iterations decreases to 10⁻³, or until ten iterations are reached.

To enhance the PP property of COCONUT, the density updated during the Newton iterations is appropriately adjusted according to the local Alfvénic velocity. This manipulation ensures that the local Alfvénic velocities calculated from the updated intermediate coronal state remain within a reasonable range, thereby enhancing the PP property of COCONUT. In this adjustment, we applied a smooth hyperbolic tangent function to gradually increase the plasma density when the local Alfvénic speed reaches a prescribed maximum Alfvénic speed, V_A, max, as follows:

$$\begin{array}{c}\hfill \rho ={\mathrm{{\rm Y}}}_{\rho }\frac{{\mathit{B}}^2}{V_{\mathrm{A},\max }^2}+(1-{\mathrm{{\rm Y}}}_{\rho }){\rho }_{\mathrm{o}},\end{array}$$(9)

where ${\mathrm{{\rm Y}}}_{\rho }=0.5+0.5·tanh(\frac{V_{\mathrm{A}}-V_{\mathrm{A},\max }}{V_{\mathrm{fac}}}·\pi )$ with $V_{\mathrm{A}}=\frac{\left|\mathit{B}\right|}{{\rho }_{\mathrm{o}}^{0.5}}$, V_A, max = 2000 km s⁻¹, and V_fac = 2 km s⁻¹. Here, the subscript ‘_o’ on ρ refers to the density updated in the Newton iteration without adjustment.

3.2. Implementation of boundary conditions

We first performed a quasi-steady-state coronal simulation constrained by a time-invariant magnetogram (Baratashvili et al. 2024), and then evolved the magnetograms to drive the following time-evolving coronal simulation.

We adopted one layer of ghost cells near the inner and outer boundaries (Lani 2008). The primitive variables in ghost cells near the inner boundary are defined as

$$\begin{array}{c}\hfill P_{\mathrm{G}}=2P_{\mathrm{BC}}-P_{\mathrm{inner}},\end{array}$$(10)

where the subscripts ‘_G’ and ‘_inner’ denote the corresponding variables at the ghost cell’s and the nearest inner cell’s centroid, respectively. In this study we set ψ_BC = 0 (Perri et al. 2022), and ρ_BC and p_BC were defined following Brchnelova et al. (2023) to improve the PP property of COCONUT:

$$\begin{array}{c}\hfill p_{\mathrm{BC}}={\mathrm{{\rm Y}}}_p\frac{{\mathit{B}}_{\mathrm{BC}}^2}{2}{\beta }_{min}+(1-{\mathrm{{\rm Y}}}_p)p_{\mathrm{s}},\end{array}$$(11)

where ${\mathrm{{\rm Y}}}_{\mathrm{p}}=0.5+0.5·tanh(\frac{{\beta }_{min}-\frac{p_{\mathrm{s}}}{0.5·{\mathit{B}}_{\mathrm{BC}}^2}}{{\beta }_{\mathit{fac}}}·\pi )$ with β_fac = 2 × 10⁻⁶ and β_min = 0.02. ρ_BC is defined in the same manner as in Eq. (9), but with ρ_o and B replaced by ρ_s and B_BC.

Considering the coronal simulations in this paper are conducted in an inertial coordinate system, we accordingly rotate the GONG-zqs magnetograms to the HCI coordinate system. In quasi-steady-state coronal simulations, the radial component of the inner-boundary magnetic field, B_BC, r, is linearly interpolated from a time-invariant magnetogram. We use cubic Hermite interpolation for time-evolving coronal simulations, which uses four magnetograms as the stencil. We generate an intermediate magnetogram at each time step and then linearly interpolate B_BC, r from this intermediate magnetogram. Meanwhile, the tangential components of the inner-boundary magnetic field, B_BC, θ and B_BC, ϕ, are defined as the corresponding values at the centroid of the nearest inner cell. Although this simplification may reduce accuracy by one order at the inner boundary, the overall accuracy of the numerical solution can still be maintained at the desired order (Gustafsson 1981; Wang et al. 2016).

Around the 2019 eclipse, the interval between two adjacent GONG-zqs magnetograms is approximately 1 hour, with fluctuations of several minutes. Considering that it is more practical to handle a batch of input magnetogram files with a uniform interval rather than non-uniform ones, we first apply cubic Hermite interpolation to the observed magnetograms to generate a series of magnetograms with a consistent 1-hour interval, which are then used to drive the time-evolving coronal simulations. In the time-evolving coronal simulations, step sizes are set to be a few minutes – shorter than the interval between adjacent magnetograms. Therefore, as detailed below, we use cubic Hermite interpolation on the input magnetograms to obtain the required inner-boundary magnetic field for each time step:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill B_{\mathrm{BC},\mathrm{r}}(t,\theta ,\varphi )& =h_{00}(t\prime )B_{\mathrm{r}}{(\theta ,\varphi )}_m\hfill \\ \hfill & +h_{10}(t\prime )(t_{m+1}-t_m){\left(\frac{\partial B_{\mathrm{r}}(\theta ,\varphi )}{\partial t}\right)}_m\hfill \\ \hfill & +h_{01}(t\prime )B_{\mathrm{r}}{(\theta ,\varphi )}_{m+1}\hfill \\ \hfill & +h_{11}(t\prime )(t_{m+1}-t_m){\left(\frac{\partial B_{\mathrm{r}}(\theta ,\varphi )}{\partial t}\right)}_{m+1}\hfill \\ \hfill \end{array}\end{array}$$(12)

with $t\prime =\frac{t-t_m}{t_{m+1}-t_m}$.

Here B_BC, r(t,θ,ϕ) is the interpolated magnetic field at the inner-boundary facial centroid located at (R_s,θ,ϕ), t denotes the physical time, and the subscripts ‘_m’ and ‘_m + 1’ refer to the m-th and (m + 1)-th magnetograms, two adjacent magnetograms nearby t, with t_m < t ≤ t_m + 1, ${\left(\frac{\partial B_{\mathrm{r}}(\theta ,\varphi )}{\partial t}\right)}_m=\frac{1}{2}(\frac{B_{\mathrm{r}}{(\theta ,\varphi )}_{m+1}-B_{\mathrm{r}}{(\theta ,\varphi )}_m}{t_{m+1}-t_m}+\frac{B_{\mathrm{r}}{(\theta ,\varphi )}_m-B_{\mathrm{r}}{(\theta ,\varphi )}_{m-1}}{t_m-t_{m-1}})$ and ${\left(\frac{\partial B_{\mathrm{r}}(\theta ,\varphi )}{\partial t}\right)}_{m+1}=\frac{1}{2}(\frac{B_{\mathrm{r}}{(\theta ,\varphi )}_{m+2}-B_{\mathrm{r}}{(\theta ,\varphi )}_{m+1}}{t_{m+2}-t_{m+1}}+\frac{B_{\mathrm{r}}{(\theta ,\varphi )}_{m+1}-B_{\mathrm{r}}{(\theta ,\varphi )}_m}{t_{m+1}-t_m})$. Besides, h₀₀(t′) = 2t′³ − 3t′² + 1, h₁₀(t′) = t′³ − 2t′² + t′, h₀₁(t′) = −2t′³ + 3t′² and h₁₁(t′) = t′³ − t′² are used to normalize the cubic Hermite interpolation formula.

Additionally, the inner-boundary velocity v = (u_BC,v_BC,w_BC) is defined as the velocity at the centroid of the nearest cell in the inner domain and constrained by a predefined speed V_BD, max. Given that an inner-boundary face, a patch of the solar surface, can encompass features like a supergranulation or a sunspot, where the typical average velocity is generally below 1 km s⁻¹, we further constrained the facial average plasma velocity at inner-boundary to not significantly exceed this velocity during our simulations:

$$\begin{array}{c}\hfill {\mathit{v}}_{\mathrm{BC}}=\frac{{\mathit{v}}_{\mathrm{inner}}}{\left|{\mathit{v}}_{\mathrm{inner}}\right|}({\mathrm{{\rm Y}}}_{\mathit{v}}{V}_{\mathrm{BD},max}+(1-{\mathrm{{\rm Y}}}_{\mathit{v}})\left|{\mathit{v}}_{\mathrm{inner}}\right|),\end{array}$$(13)

where ${\mathrm{{\rm Y}}}_{\mathit{v}}=0.5+0.5·tanh(\frac{\left|{\mathit{v}}_{\mathrm{inner}}\right|-V_{\mathrm{BD},max}}{V_{\mathrm{BD},\mathrm{fac}}}·\pi )$ with V_BD, max = 1 km s⁻¹ and V_BD, fac = 2 m s⁻¹. Here v_inner refers to the velocity of plasma flow at the centroid of the nearest cell in the inner domain.

The outer boundary of the present model is set to 25 R_s. Since this distance is where the plasma flow has become supersonic and super-Alfvénic, we prescribed the Neumann boundary conditions at the outer boundary (Brchnelova et al. 2022b). Therefore, we extrapolated r²B_r, B_θ, B_ϕ, ρ, u, v, wp, and ψ from the outermost cell centres in the computational domain to the ghost cells with a zero gradient (Brchnelova et al. 2022a; Perri et al. 2022).

4. Numerical results

In this section, the time-evolving coronal model developed in this paper is employed to simulate the evolution of the coronal structures during CRs 2219 and 2220, encompassing the solar eclipse on July 2, 2019. About 1300 GONG-zqs magnetograms², spanning the period from 11:00 on June 29, 2019, to 21:00 on August 22, 2019, are used to drive these simulations. As mentioned in Sect. 2.2, the initial magnetic fields for the quasi-steady-state coronal simulations are achieved from a PF model, with the bottom boundary condition specified by the synoptic map of the radial photospheric magnetic field centred on 11:00 on June 29, 2019. Following the quasi-steady-state coronal simulation, the time-evolving photospheric magnetograms are then used to drive the time-evolving coronal simulations, extending from the solar surface to 25 R_s throughout these two CRs within an inertial coordinate system.

In Sect. 4.1, we perform quasi-steady-state and time-evolving MHD coronal simulations and compare their simulation results to validate the time-evolving approach. It demonstrates that, compared with the quasi-steady-state simulation constrained by only one time-invariant magnetogram during a CR period, there are pronounced differences in the time-evolving simulation results. In Sect. 4.2, we evaluate the impact of time-step sizes on the time-evolving coronal simulation results. It indicates that using a moderately larger time-step size can significantly enhance computational efficiency with minimal impact on accuracy.

In this paper, all the calculations are performed on the KU Leuven/UHasselt Tier-2 cluster wICE of the Vlaams Supercomputer Centrum (VSC)³. Utilising 1080 CPU cores, the wall-clock time for the time-evolving coronal simulations covering two CRs was approximately 17.54 hours with a 10-minute time-step size and 39.06 hours with a 2-minute time-step size. In the following subsections, we present the results of the MHD coronal simulations for CRs 2219 and 2220.

4.1. Time-evolving versus quasi-steady-state coronal simulation results

In this subsection we compare the time-evolving simulation results at different moments with the quasi-steady-state simulation, which is constrained by a magnetogram at 19:00 on July 2, 2019 – the day of the 2019 eclipse – corresponding to the magnetogram at 82 hours of the time-evolving simulation. The time-evolving coronal simulation presented here adopts a constant time-step size of 10 minutes.

Movie 2 (linked in Fig. 4), spanning from the 82nd to the 735th hour of the time-evolving simulation, depicts the evolution of some selected magnetic field lines within the HCI coordinate system. Movie 1 (linked in Fig. 2) illustrates the evolution of polarisation brightness (pB) images synthesised from the simulation results, as viewed from the perspective of STEREO-A. Furthermore, we transformed the simulation results from the HCI coordinate system into the heliographic co-rotating coordinate system to evaluate the relative differences between the time-evolving and quasi-steady-state simulations. The evolution of these differences is visualised in movies 3 and 4 (linked in Fig. 5). Additionally, in Fig. 6, we present the variation patterns of some simulated plasma parameters in the high-density, low-speed (HDLS) and low-density, high-speed (LDHS) regions. The figure also compares the simulated velocity and observed data, mapped from 1 AU to 20 R_s using ballistic propagation.

Fig. 2. White-light polarisation brightness (pB) images observed from COR2 (STEREO-A) on July 2, 2019 (a), corresponding pB images synthesised from quasi-steady-state coronal simulation results (b), and synthesised from simulation results at 82 (c) and 735 (d) hours of the time-evolving coronal simulations. These synthesised images span from 2.5 to 15 Rs on the meridian planes perpendicular to the line of sight from STEREO-A. Orange lines highlight magnetic field lines on these selected meridional planes. The evolution of simulated pB images over this time interval is showcased in online movie 1.

In Fig. 2 we compare white-light pB images from 2.5 to 15 R_s that are observed from the outer coronagraph of the Sun-Earth Connection Coronal and Heliospheric Investigation (SECCHI) instrument suite (panel a) on board the STEREO-A spacecraft (Frazin et al. 2012; Howard et al. 2008; Kaiser et al. 2008; Thompson et al. 2003; Thompson & Reginald 2008) and the white-light pB images synthesised from the quasi-steady-state b and time-evolving (c and d) coronal simulation results. The top panel demonstrates the pB image observed on July 2, 2019 a, alongside the image synthesised from the simulation result of quasi-steady-state simulation results b. The bottom panel presents pB images synthesised from the time-evolving simulation results at the 82nd hour c and the 735th hour d, with an interval of one CR period between the two images. The observed and simulated images reveal four bright structures near the solar equator. However, compared to the observation, the locations of these simulated bright structures in the northern hemisphere are farther north, and the simulated bright structures at the south-east limb extend farther outwards. This discrepancy may be attributed to inaccurate observations for both polar photospheric magnetic fields.

Additionally, it can be seen that when confined by the same magnetogram at the inner boundary, the pB images synthesised from the time-evolving and quasi-steady-state coronal simulation results show only minor differences. After a CR period, the bright structure at the south-west limb shifts about 10° towards the equator, aligning more closely with the observed results. It demonstrates that bright structures are typically associated with closed-field regions. Additionally, interested readers can refer to online movie 1 to see the evolution of simulated pB images observed from the COR2 (Coronagraph 2 instrument) on STEREO-A field of view between the 82nd and 735th hours of the time-evolving coronal simulation.

In Fig. 3 we present the synthesised white-light pB images at 3 R_s. These images were synthesised from observations during CRs 2219 (panel a) and 2210 (panel b) and displayed in the co-rotating coordinate system overlaid with magnetic field neutral lines (MNLs) from the quasi-steady-state simulation result and from the 82nd and 735th hours of the time-evolving simulation results, respectively. It can be seen that when constrained by the same magnetogram at the inner boundary, the MNLs calculated from the quasi-steady-state and time-evolving simulations are nearly identical. After one CR period, the observed bright structure between 240° and 310° in longitude becomes more concave; the MNLs also capture this variation from the time-evolving coronal simulation. Also, it reveals that between 20° and 240° in longitude, the segments of these MNLs calculated by MHD models are approximately 10° farther north than the central axis of the observed bright structures. These discrepancies may be attributed to the empirical, rather than self-consistent, formulation of certain mechanisms in this MHD coronal model – such as coronal heating and solar wind acceleration – which affect magnetic pressure and, consequently, the distribution of MNLs.

Fig. 3. Synoptic maps of white-light pB observations from the SOHO instrument LASCO C2 at 3 Rs for CRs 2219 (left) and 2220 (right). The dashed yellow, dashed black, and solid grey lines represent the MNLs derived from the quasi-steady-state simulations and time-evolving simulations at 82nd and 735th hours, respectively.

In Fig. 4 we shows the distributions of 2D magnetic field lines. These radial velocity contours range from 1 to 20 R_s on the same meridians as in Fig. 2. The top-left panel a denotes the quasi-steady-state simulation result, the top-right and bottom-left panels present the time-evolving simulation results at the 82nd and 735th hours, respectively, and the bottom-right panel d demonstrates the 3D distribution of some selected magnetic field lines. Together with Fig. 2, this figure reveals that close-field structures in low-latitude regions dominate the HDLS flow. In contrast, LDHS flows are prevalent in the high-latitude areas, characterised by open-field structures. This is a crucial characteristic of solar minima.

Fig. 4. 2D distribution of magnetic field lines overlaid on contours of the radial plasma speed, Vr (km s−1), illustrated on meridional planes spanning from 1 to 20 Rs and perpendicular to the STEREO-A line of sight (a, b, and c), and 3D distribution of magnetic field lines (d) at the 82nd hour of the time-evolving simulation. Panel a shows the results of the quasi-steady-state simulation, panels b and c present the time-evolving simulation results at the 82nd and 735th hours, respectively, and panel d the 3D distribution of some selected magnetic field lines. The evolution of these selected magnetic field lines is shown in online movie 2.

Besides, it can be seen that after a CR period, the open magnetic field lines in the eastern hemisphere, spanning from 20 ° S to 15 ° N, change their direction from pointing outwards from the Sun to pointing towards the Sun. Significant variations are also evident in the closed field lines near the Sun. Furthermore, interested readers can refer to online movie 2 to see the evolution of some simulated magnetic field lines in a 3D representation within the HCI coordinate system over the 82nd and 735th hours of the time-evolving coronal simulation. The continuous formation and disappearance of closed field lines, a phenomenon absent in quasi-steady-state coronal simulations, is evident during this time-evolving process.

In Fig. 5 we present the distribution of radial velocity (panel a) and proton number density b at 20 R_s calculated at the 82nd hour of the time-evolving coronal simulation. We also illustrate relative differences in radial velocity ${\mathrm{RD}}_{V_{\mathrm{r}}}=\frac{V_{\mathrm{r}}^{\mathrm{TE}}-V_{\mathrm{r}}^{\mathrm{QSS}}}{V_{\mathrm{r}}^{\mathrm{QSS}}}$ c and plasma density ${\mathrm{RD}}_{\rho }=\frac{{\rho }^{\mathrm{TE}}-{\rho }^{\mathrm{QSS}}}{{\rho }^{\mathrm{QSS}}}$ d between the time-evolving and the quasi-steady-state coronal simulations. Furthermore, movies 3 and 4 illustrate the evolution of the relative differences of these parameters. In this comparison, the simulation results at 20 R_s during 82 and 735 hours of the time-evolving simulation are compared with the quasi-steady-state coronal simulation constrained by a time-invariant magnetogram. It can be seen that the relative differences in radial velocity RD_Vr and plasma density RD_ρ can reach 15% and 25%. Additionally, in Table 1 we list the average relative differences in plasma density RD_ave, ρ, velocity RD_ave, |v|, and magnetic field strength RD_ave, |B| between the steady-state and the time-evolving coronal simulation results at the 82nd and 735th hours, respectively. Here

$$\begin{array}{cc}\hfill {\mathrm{RD}}_{\mathrm{ave},\rho }& =\frac{\sum _{i=1}^N|{\rho }_{i,\mathrm{QSS}}-{\rho }_{i,\mathrm{TE}}|}{\sum _{i=1}^N{\rho }_{i,\mathrm{QSS}}},\hfill \\ \hfill {\mathrm{RD}}_{\mathrm{ave},\left|\mathit{v}\right|}& =\frac{\sum _{i=1}^N|{\mathit{v}}_{i,\mathrm{QSS}}-{\mathit{v}}_{i,\mathrm{TE}}|}{\sum _{i=1}^N\left|{\mathit{v}}_{i,\mathrm{QSS}}\right|},\hfill \\ \hfill {\mathrm{RD}}_{\mathrm{ave},\left|\mathit{B}\right|}& =\frac{\sum _{i=1}^N|{\mathit{B}}_{i,\mathrm{QSS}}-{\mathit{B}}_{i,\mathrm{TE}}|}{\sum _{i=1}^N\left|{\mathit{B}}_{i,\mathrm{QSS}}\right|}.\hfill \end{array}$$

Fig. 5. Distribution of the radial plasma speeds, Vr, in units of km s−1 (a) and proton number density in units of 103 cm−3 (b) at 20 Rs derived from the 82nd hour of the time-evolving coronal simulation. Panels c and d show the relative differences in radial velocity and plasma density between the quasi-steady-state and time-evolving coronal simulations. Online movies 3 and 4 illustrate the evolution of these relative differences. The solid black lines and the dashed orange lines represent the MNLs calculated by the time-evolving and quasi-steady-state coronal models.

The superscript ‘^t’ means the t-th hour of the time-evolving simulation, the subscript ‘_i, QSS’ or ‘_i, TE’ denotes the corresponding variable in cell i calculated by the quasi-steady-state or time-evolving coronal models, and N is the number of cells in the computational domain. It can be seen that the overall differences in the simulation plasma densities, velocities, and magnetic field strength between the quasi-steady-state and time-evolving coronal models are relatively small when constrained by the same magnetogram. However, these differences can increase significantly during a CR period, particularly in magnetic field strength, reaching up to 35%.

In Fig. 6 we display the timing diagrams of radial velocity (panels a and d), proton number density (b and e) and magnetic field strength (c and f) during the 82nd and 735th hours of the time-evolving coronal simulation by solid black lines. These parameters are observed by two virtual satellites placed at (r,θ,ϕ) = (20 R_s,0,201°) and (r,θ,ϕ) = (20 R_s,−70°,201°), which located in the HDLS (q, b, c) and LDHS (d, e, f) regions, respectively. Additionally, we mapped the heliolongitude in quasi-steady-state coronal simulations to a CR period of physical time as

$$\begin{array}{c}\hfill t=\{\begin{array}{c}82-653\frac{\varphi -{\varphi }_0}{{360}^{°}},\text{if}\varphi \le {\varphi }_0\\ 735-653\frac{\varphi -{\varphi }_0}{{360}^{°}},\text{if}\varphi >{\varphi }_0\end{array},\end{array}$$

Fig. 6. Timing diagram of simulated radial velocity, Vr, in km s−1 (a, d), proton number density in units of 103 cm−3 (b, e), and magnetic field strength in unit of Gauss (c, f) observed by two virtual satellites located at the HDLS (a, b, c) and LDHS (d, e, f) regions. The solid black lines denote variables observed at (r,θ,ϕ) = (20 Rs,0,201°) and (r,θ,ϕ) = (20 Rs,−70°,201°) during the 82nd and 735th hours of the time-evolving coronal simulation. For the quasi-steady-state coronal simulations constrained by magnetograms at 82 (denoted by dashed black lines) and 735 (denoted by dash-dot black lines) of this time interval, the heliolongitude is mapped to a CR period corresponding to the time series on the horizontal axis of the timing diagram. The solid grey line (a) denotes the velocity observed by the WIND satellite and mapped from 1 AU to 20 Rs following the ballistic propagation.

where ϕ₀ = 201°. The selected variables were calculated using the quasi-steady coronal models, constrained by magnetograms at 19:00 on July 2, 2019, and 00:00 on July 30, 2019, along the mapped time at (r,θ) = (20 R_s,0). It reveals that the peaks and troughs in the timing diagrams from the time-evolving simulation differ from those in the quasi-steady-state coronal simulations, with these differences being more pronounced in the LDHS regions than in the HDLS regions. In the LDHS regions, the peaks and troughs observed around the 160th, 220th, 320th, and 660th hours in the time-evolving coronal model are absent in the quasi-steady-state simulations. Additionally, in the HDLS regions, the troughs in the magnetic field strength timing diagrams occur earlier in the time-evolving simulation than those from the quasi-steady-state models.

Furthermore, we made a ballistic propagation to map the observed velocity at 1 AU around the equatorial plane to 20 R_s and compared them with the simulation results, as shown in the left-top panel of Fig. 6. Both quasi-steady and time-evolving simulations capture the velocity peak between 250 and 490 hours of the coronal simulation. However, between 300 and 330 hours and 370 and 470 hours, the time-evolving simulation result aligns more closely with the observation. It shows that the time-evolving coronal model has the potential to make simulation results coincide more closely with observations.

4.2. The impact of time-step sizes on simulation results

In this subsection we evaluate the impact of time-step sizes on the time-evolving coronal simulation results. The model performs well with a time-step size of 10 minutes or less. However, increasing the time steps to larger intervals, such as 20 minutes, can cause the simulation to crash under certain conditions. To further evaluate the impact of time-step sizes on accuracy and efficiency, we conducted a time-evolving coronal simulation with a smaller time-step size of 2 minutes and compared the simulated plasma density, velocity, and magnetic field strength with those from the previously discussed simulation, which used a 10-minute time step.

In the time-evolving coronal simulation with time steps of 10 minutes, approximately 7–9 Newtonian iterations, as described in Sect. 3.1, are typically performed in each time step. However, instances exist where the L₂ norm of the differences between the state variables updated in two consecutive Newton iterations fluctuates around a value larger than the threshold used to stop the Newtonian iterations. In such cases, we terminated the Newton iterations at the tenth iteration. For the time-evolving simulation with 2-minute time steps, approximately five to six Newtonian iterations are typically performed per time step; the computational time was about 2.23 times longer than in the simulation with a 10-minute time step.

In Fig. 7 we show the relative differences in plasma density RD_ρ (panels a and b) and radial velocity RD_Vr (c and d) between the results calculated with time steps of 2 and 10 minutes at 20 R_s. Here, RD_ρ and RD_Vr are defined similar to those in Fig. 5, but with the variables from the quasi-steady-state simulation replaced by their counterparts from the time-evolving simulation with a time step of 2 minutes. Overlaid on these contours are the MNLs calculated with time steps of 2 and 10 minutes. It can be seen that the MNLs calculated with different time-step sizes coincide with each other. In most regions, the relative differences in plasma density and radial velocity are less than 0.2% and 0.1%. Additionally, it reveals that the most significant relative differences mainly appear around the MNLs.

Fig. 7. Distribution of the relative differences in plasma density (a, b) and radial velocity (c, d) at 20 Rs between the results calculated with time steps of 2 and 10 minutes. The solid black and dashed orange lines denote the MNLs calculated with time steps of 2 and 10 minutes, respectively. The left and right panels denote the results at the 82nd and 735th hour of the time-evolving coronal simulations.

Furthermore, Table 2 presents the average relative differences in plasma density ${\mathrm{RD}}_{\mathrm{ave},\rho }^t$, velocity ${\mathrm{RD}}_{\mathrm{ave},\left|\mathit{v}\right|}^t$ and magnetic field ${\mathrm{RD}}_{\mathrm{ave},\left|\mathit{B}\right|}^t$ between the simulation results calculated with time steps of 2 and 10 minutes. Here, ${\mathrm{RD}}_{\mathrm{ave},\rho }^t$, ${\mathrm{RD}}_{\mathrm{ave},\left|\mathit{v}\right|}^t$ and ${\mathrm{RD}}_{\mathrm{ave},\left|\mathit{B}\right|}^t$ are defined similarly to those in Table 1, except that the variables from the quasi-steady-state simulation are replaced by their counterparts from the time-evolving simulation with a time-step size of 2 minutes. Additionally, the table includes the wall-clock time required for both time-evolving coronal simulations. It indicates that the overall differences between the simulation results with dt = 10 min and dt = 2 min are minimal, remaining below 1%. In contrast, the computational time with dt = 2 min is approximately 2.23 times as long as that with dt = 10 min. This demonstrates that the computational efficiency of the implicit time-evolving coronal model can be significantly enhanced by increasing the time-step size while still maintaining the required level of temporal accuracy.

5. Concluding remarks

The commonly used steady-state or time-dependent coronal models are typically constrained by a time-invariant magnetogram, which fails to reflect that the corona evolves. This simplification usually discards crucial information needed to simulate solar disturbances accurately, such as CME propagation. It also leads to discrepancies compared with 1 AU observations, resulting in an overly simplistic representation of magnetic connectivity. In contrast, time-evolving coronal models, driven by a series of more realistic time-evolving magnetograms, can capture the dynamic features of the corona with higher fidelity, offering the potential to address these shortcomings.

In this study we extended the recently developed COCONUT coronal model, constrained initially by a time-invariant magnetogram, to a more realistic time-evolving coronal model that is driven by a sequence of time-varying magnetograms, enabling it to capture the dynamic features of the solar corona. It is the first fully implicit time-evolving MHD coronal model, offering the flexibility to use large time-step sizes exceeding the CFL condition. This paper demonstrates that accurate real-time simulations of the time-evolving global corona can be achieved using no more than 20 CPU cores. This approach is considerably more efficient than explicit or semi-implicit coronal models.

A series of hourly updated GONG-zqs magnetograms around the 2019 eclipse drive the time-evolving coronal simulations in an inertial coordinate system. The cubic Hermite interpolation method reconstructs the inner-boundary magnetic field at each physical time step. Furthermore, we designed an HLL Riemann solver with a self-adjustable dissipation term, which reduces the dissipation term for low-speed flows and reverts to the original one for high-speed flows. This method minimises the degradation of the solution accuracy in low-speed regions and improves the numerical stability of our coronal model. Additionally, during the time-evolving coronal simulations, the model adaptively adjusts the plasma density to avoid the abnormally high Alfvénic speed. This strategy enhances the PP property of the present coronal model.

The simulation results demonstrate that this implicit time-evolving coronal model successfully replicates observed coronal structures, including those captured in pB images and in situ observations. Additionally, the comparison with a quasi-steady-state coronal simulation, constrained by a single time-invariant magnetogram during a CR period, demonstrates that the time-evolving coronal model with a 10-minute time step can produce coronal structures comparable to those generated by the quasi-steady-state model at the moment when both models use the same magnetogram to define their inner boundaries. The average relative differences in plasma density, velocity, and magnetic field strength between the time-evolving and quasi-steady-state coronal simulations are less than 0.6%, 1.3%, and 4.3%, respectively. However, these differences can increase significantly during a CR period, with the magnetic field strength differences reaching as high as 35%. Consequently, we can conclude that the time-evolving coronal model with a significant time-step size not only reproduces coronal structures comparable to those generated by the quasi-steady-state model at the moment confined by the same magnetogram at the inner boundary; it also captures the continuous temporal evolution of these structures, a capability absent in the quasi-steady-state model.

Additionally, we evaluated the impact of time-step size on the simulation results. The average relative differences in plasma density, velocity, and magnetic field strength between simulations with time steps of 2 minutes and 10 minutes are no more than 0.09%, 0.10%, and 0.60%, respectively. Meanwhile, the computational time required for a 10-minute time-step size is less than half of that needed for a 2-minute time-step size. We can conclude that increasing the time-step size can significantly enhance the computational efficiency of the fully implicit time-evolving coronal model while still preserving the required temporal accuracy. We also notice that excessively large time-step sizes, such as 20 minutes, can occasionally cause the code to crash. This time-evolving coronal model is designed to be coupled with an inner heliosphere model in the future. The implicit model’s flexibility in selecting significant time steps will allow us to determine the time steps of the MHD coronal model via the inner heliosphere model. It is well known that most MHD inner heliosphere models use explicit numerical schemes, with time steps limited by the CFL stability criterion to around 10 minutes. Thus, a 10-minute time step is efficient and numerically stable enough for our research.

Given that the computational time of 1281 hours of physical time is no more than 18 hours when adopting a time-step size of 10 minutes, it is practical to adopt a smaller time step to achieve more accurate simulation results for higher-frequency events such as CMEs, albeit at the cost of a manageable reduction in computational efficiency. The physical time-step size can be adjusted based on the accuracy required for specific research or applications. This allows for optimised time-step sizes that balance computational efficiency with the necessary numerical stability and accuracy. This approach enabled us to model the entire solar-terrestrial domain with a single model, thereby streamlining the Sun-to-Earth forecasting process.

Although this established fully implicit time-evolving MHD solar coronal model has many advantages and is a promising tool for timely and accurately simulating the time-evolving corona in practical space weather forecasting, there remains considerable room for improvement. In future work, we plan to couple this highly efficient time-evolving coronal model with an inner heliosphere model, such as EUHFORIA, to further validate its potential for enhancing the performance of the current Sun-to-Earth model chain in providing timely and accurate space weather forecasting. Additionally, we plan to incorporate surface flux transport models, such as the advective flux transport model (Upton & Hathaway 2014), to advect the radial magnetic field with observed flows. This approach may help generate a more realistic magnetic field evolution at the inner boundary of our coronal model. Furthermore, we may also try to incorporate additional observations into the inner boundary conditions, such as inferring horizontal velocities from observational data using the time-distance helioseismology method (Yalim et al. 2017; Zhao et al. 2012). Additionally, reducing magnetic field discretisation errors by solving decomposed MHD equations (e.g. Feng et al. 2010; Li et al. 2018; Wang et al. 2019a) could enhance the numerical stability of COCONUT, even when addressing extremely low-β problems.

Data availability

Movies associated to Figs. 2, 4, and 5 are available at https://www.aanda.org

Acknowledgments

This project has received funding from the European Research Council Executive Agency (ERCEA) under the ERC-AdG agreement No. 101141362 (Open SESAME). These results were also obtained in the framework of the projects FA9550-18-1-0093 (AFOSR), C16/24/010 (C1 project Internal Funds KU Leuven), G0B5823N and G002523N (WEAVE) (FWO-Vlaanderen), and 4000145223 (SIDC Data Exploitation (SIDEX), ESA Prodex). This work is also part of a project supported by the Specialized Research Fund for State Key Laboratories, managed by the Chinese State Key Laboratory of Space Weather. It is also supported by the National Natural Science Foundation of China (grant Nos. 42030204, 42074208 and 42104168). 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. The Research Council of Norway supports F.Z. through its Centres of Excellence scheme, project No. 262622. This work utilises data obtained by the Global Oscillation Network Group (GONG) program, managed by the National Solar Observatory and operated by AURA, Inc., under a cooperative agreement with the National Science Foundation. The data were acquired by instruments operated by the Big Bear Solar Observatory, High Altitude Observatory, Learmonth Solar Observatory, Udaipur Solar Observatory, Instituto de Astrofísica de Canarias, and Cerro Tololo Inter-American Observatory. The authors also acknowledge the use of the STEREO/SECCHI data produced by a consortium of the NRL (US), LMSAL (US), NASA/GSFC (US), RAL (UK), UBHAM (UK), MPS (Germany), CSL (Belgium), IOTA (France), and IAS (France).

References

Appendix A: HLL Riemann solver with a self-adjustable dissipation term

Approximate Riemann solvers, such as HLL- (Feng et al. 2021) and AUSM-type (Kitamura 2020; Wang et al. 2022a) solvers, can be used to compute the inviscid flux F_n, ij between two cells i and j. The subscript _i, j is omitted hereafter for simplicity. In this study we adopted the HLL Riemann solver, which can be described as follows:

$$\begin{array}{c}\hfill {\mathit{F}}_n({\mathit{U}}_L,{\mathit{U}}_R)=\frac{{\mathit{F}}_n\left({\mathit{U}}_L\right)+{\mathit{F}}_n\left({\mathit{U}}_R\right)}{2}-\frac{1}{2}{\mathit{D}}_{\mathrm{HLL}}({\mathit{U}}_L,{\mathit{U}}_R),\end{array}$$(A.1)

with

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathit{D}}_{\mathrm{HLL}}({\mathit{U}}_L,{\mathit{U}}_R)=& \frac{(S_L+S_R)({\mathit{F}}_n\left({\mathit{U}}_R\right)-{\mathit{F}}_n\left({\mathit{U}}_L\right))}{S_R-S_L}\hfill \\ \hfill & -\frac{2S_R{S}_L}{S_R-S_L}({\mathit{U}}_R-{\mathit{U}}_L),\hfill \end{array}\end{array}$$(A.2)

and

$$\begin{array}{c}\hfill {\mathit{F}}_n\left(\mathit{U}\right)=\left(\begin{array}{c}{V}_n\rho \\ V_n\rho u+p_T{n}_x-B_n{B}_x\\ V_n\rho v+p_T{n}_y-B_n{B}_y\\ V_n\rho w+p_T{n}_z-B_n{B}_z\\ (v{B}_x-B_{y}u)n_y+(w{B}_x-B_{z}u)n_z+\psi n_x\\ (u{B}_y-B_{x}v)n_x+(w{B}_y-B_{z}v)n_z+\psi n_y\\ (u{B}_z-B_{x}w)n_x+(v{B}_z-B_{y}w)n_y+\psi n_z\\ V_n(E+p_T)-B_n\mathit{B}·\mathit{v}\\ B_n{V}_{\mathrm{ref}}^2\end{array}\right).\end{array}$$

Here, V_n(U) = v ⋅ n_ij and B_n(B) = B ⋅ n_ij denote the velocity and magnetic field along the normal direction of Γ_ij, S_L and S_R usually stand for the conventional two fast waves in the HLL Riemann solver (Einfeldt et al. 1991). Unless otherwise stated, the subscripts ‘_L’ and ‘_R’ denote the corresponding left and right variables evaluated on the centroid of Γ_ij. In this paper, we set S_L = min(0,λ_min(U_L),λ_min(U_R)) and S_R = max(0,λ_max(U_L),λ_max(U_R)) with

$${\lambda }_{min}\left(\mathit{U}\right)=min(V_n\left(\mathit{U}\right),V_n\left(\mathit{U}\right)\pm c_{f\mathrm{or}s}\left(\mathit{U}\right),V_n\left(\mathit{U}\right)\pm c_A\left(\mathit{U}\right)),$$

and

$${\lambda }_{max}\left(\mathit{U}\right)=max(V_n\left(\mathit{U}\right),V_n\left(\mathit{U}\right)\pm c_{f\mathrm{or}s}\left(\mathit{U}\right),V_n\left(\mathit{U}\right)\pm c_A\left(\mathit{U}\right)),$$

where $c_{f\mathrm{or}s}\left(\mathit{U}\right)=\sqrt{0.5(\frac{\gamma p}{\rho }+\frac{{\mathit{B}}^2}{\rho }\pm \sqrt{{(\frac{\gamma p}{\rho }+\frac{{\mathit{B}}^2}{\rho })}^2-4\frac{\gamma p}{\rho }\frac{B_n^2}{\rho }})}$, and $c_A\left(\mathit{U}\right)=\frac{\left|B_n\right|}{{\rho }^{0.5}}$.

Although the HLL Riemann solver performs well in preserving the PP property for compressible high-speed flows with large Mach numbers, its diffusion is excessive for incompressible low-speed flows with small Mach numbers. Due to excessive diffusion in the original solver, the accumulation of solution accuracy degradation in low-speed flow regions may eventually lead to code crashes. A suitable low-dissipation numerical flux solver is required for low-speed flows to avoid a degradation of the solution accuracy (e.g. Esquivel et al. 2010; Kitamura et al. 2011; Kitamura & Balsara 2018; Minoshima et al. 2020; Minoshima & Miyoshi 2021; Wang et al. 2022a). Therefore, we introduced the factor

$$\begin{array}{c}\hfill \phi =\frac{max(\left|S_L\right|,\left|S_R\right|)}{S_R-S_L}\end{array}$$(A.3)

to the second term on the right-hand side of Eq. (A.2), which plays a major role in low Mach number regions and decreases to zero in high Mach number regions to reduce the diffusion of HLL Riemann solver for low-speed flow. Consequently, we get the following HLL Riemann solver,

$$\begin{array}{c}\hfill {\mathit{F}}_{n,ij}({\mathit{U}}_L,{\mathit{U}}_R)=\frac{{\mathit{F}}_n\left({\mathit{U}}_L\right)+{\mathit{F}}_n\left({\mathit{U}}_R\right)}{2}-\frac{1}{2}{\mathit{D}}_{\mathrm{HLL}}\prime ({\mathit{U}}_L,{\mathit{U}}_R),\end{array}$$(A.4)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathit{D}}_{\mathrm{HLL}}\prime ({\mathit{U}}_L,{\mathit{U}}_R)=& \frac{(S_L+S_R)({\mathit{F}}_n\left({\mathit{U}}_R\right)-{\mathit{F}}_n\left({\mathit{U}}_L\right))}{S_R-S_L}\hfill \\ \hfill & -\phi \frac{2S_R{S}_L}{S_R-S_L}({\mathit{U}}_R-{\mathit{U}}_L).\hfill \end{array}\end{array}$$(A.5)

The factor φ will reduce the dissipation term of the HLL Riemann solver by half for low-speed flows and recover to the original one for high-speed flows. We find that incorporating this factor enhances the numerical stability of our coronal model, as the model would be prone to crashing without this adjustment.

All Tables

All Figures

	Fig. 1. Radiative cooling curve profile derived from version 9 of the CHIANTI atomic database (Dere et al. 2019). The horizontal axis denotes the decadic logarithm of temperature, and the vertical axis the logarithm of the radiative cooling curve function value. The temperature and radiative cooling curve function units are K and erg s−1 cm3, respectively.
In the text

Fig. 2. White-light polarisation brightness (pB) images observed from COR2 (STEREO-A) on July 2, 2019 (a), corresponding pB images synthesised from quasi-steady-state coronal simulation results (b), and synthesised from simulation results at 82 (c) and 735 (d) hours of the time-evolving coronal simulations. These synthesised images span from 2.5 to 15 Rs on the meridian planes perpendicular to the line of sight from STEREO-A. Orange lines highlight magnetic field lines on these selected meridional planes. The evolution of simulated pB images over this time interval is showcased in online movie 1.

In the text

	Fig. 3. Synoptic maps of white-light pB observations from the SOHO instrument LASCO C2 at 3 Rs for CRs 2219 (left) and 2220 (right). The dashed yellow, dashed black, and solid grey lines represent the MNLs derived from the quasi-steady-state simulations and time-evolving simulations at 82nd and 735th hours, respectively.
In the text

Fig. 4. 2D distribution of magnetic field lines overlaid on contours of the radial plasma speed, Vr (km s−1), illustrated on meridional planes spanning from 1 to 20 Rs and perpendicular to the STEREO-A line of sight (a, b, and c), and 3D distribution of magnetic field lines (d) at the 82nd hour of the time-evolving simulation. Panel a shows the results of the quasi-steady-state simulation, panels b and c present the time-evolving simulation results at the 82nd and 735th hours, respectively, and panel d the 3D distribution of some selected magnetic field lines. The evolution of these selected magnetic field lines is shown in online movie 2.

In the text

Fig. 5. Distribution of the radial plasma speeds, Vr, in units of km s−1 (a) and proton number density in units of 103 cm−3 (b) at 20 Rs derived from the 82nd hour of the time-evolving coronal simulation. Panels c and d show the relative differences in radial velocity and plasma density between the quasi-steady-state and time-evolving coronal simulations. Online movies 3 and 4 illustrate the evolution of these relative differences. The solid black lines and the dashed orange lines represent the MNLs calculated by the time-evolving and quasi-steady-state coronal models.

In the text

Fig. 6. Timing diagram of simulated radial velocity, Vr, in km s−1 (a, d), proton number density in units of 103 cm−3 (b, e), and magnetic field strength in unit of Gauss (c, f) observed by two virtual satellites located at the HDLS (a, b, c) and LDHS (d, e, f) regions. The solid black lines denote variables observed at (r,θ,ϕ) = (20 Rs,0,201°) and (r,θ,ϕ) = (20 Rs,−70°,201°) during the 82nd and 735th hours of the time-evolving coronal simulation. For the quasi-steady-state coronal simulations constrained by magnetograms at 82 (denoted by dashed black lines) and 735 (denoted by dash-dot black lines) of this time interval, the heliolongitude is mapped to a CR period corresponding to the time series on the horizontal axis of the timing diagram. The solid grey line (a) denotes the velocity observed by the WIND satellite and mapped from 1 AU to 20 Rs following the ballistic propagation.

In the text

Fig. 7. Distribution of the relative differences in plasma density (a, b) and radial velocity (c, d) at 20 Rs between the results calculated with time steps of 2 and 10 minutes. The solid black and dashed orange lines denote the MNLs calculated with time steps of 2 and 10 minutes, respectively. The left and right panels denote the results at the 82nd and 735th hour of the time-evolving coronal simulations.

In the text