© The Authors 2025

1. Introduction

In Sun-to-Earth model chains that couple multiple components to build solar–terrestrial space-weather forecasting frameworks (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), coronal models are critical for initialising the other components (Brchnelova et al. 2022a; Kuźma et al. 2023; Perri et al. 2023), and physics-based magnetohydrodynamic (MHD) coronal models are typically the most complex and computationally intensive component (Wang et al. 2025a,b). Consequently, the overall efficiency and reliability of the entire Sun-to-Earth model chain are significantly affected by the performance of the coronal model.

Generally speaking, coronal and solar wind MHD models can be classified into two types. We classify models that are constrained by a single static magnetogram as quasi-steady-state models, while those that are driven by a sequence of time-evolving magnetograms are referred to as time-evolving models (e.g. Feng et al. 2015; Feng 2020; Lionello et al. 2023; Wang et al. 2025b). Quasi-steady-state models assume that the structures of the coronal and solar wind do not evolve during a specific time interval and are typically used to compute the coronal and solar wind ambience via time-relaxation iterations. In particular, time-evolving coronal and solar wind models are time-accurate and driven by a series of time-evolving magnetograms, which allows us to simulate more realistic and continuously evolving coronal and solar wind structures. Time-evolving modelling is an essential ingredient of coronal mass ejection (CME) modelling.

Recently, implicit temporal discretisation strategies that allow large time steps that exceed the Courant-Friedrichs-Lewy stability condition have significantly increased the computational efficiency of quasi-steady-state (Brchnelova et al. 2023a; Feng et al. 2021; Kuźma et al. 2023; Liu et al. 2023; Linan et al. 2023; Perri et al. 2018, 2022, 2023; Wang et al. 2019, 2022a,b, 2025a,c) and time-evolving (Wang et al. 2025b,d) MHD coronal models. Most of the time-evolving coronal models still rely on explicit or semi-implicit temporal integration methods (Feng et al. 2023; Hayashi et al. 2021; Hoeksema et al. 2020; Linker et al. 2024; Lionello et al. 2023; Mason et al. 2023; Yang et al. 2012), however, where the size of the time-step is constrained by the explicitly treated terms. This makes these simulations prohibitively computationally expensive for practical applications (Yeates et al. 2018). Wang et al. (2025b) successfully extended the quasi-steady-state coronal model called the coolfluid coronal unstructured (COCONUT) model to the first fully implicit time-evolving coronal model for unstructured meshes. COCONUT is a novel implicit MHD solar corona model based on computational object-oriented libraries for fluid dynamics (COOLFluiD) (Kimpe et al. 2005; Lani et al. 2005, 2013)¹. Although it is efficient enough for real-time global coronal simulations during solar minima, its numerical stability is still poor for solar maximum coronal simulations that involve strong complex magnetic fields with low-β values (the ratio of thermal pressure to magnetic pressure), and it risks developing abnormally high-speed streams above active regions (ARs) when the magnetograms are only poorly preprocessed (Kuźma et al. 2023).

The numerical stability of quasi-steady-state MHD coronal models has been significantly improved by adopting decomposed MHD equations (e.g. Feng et al. 2010, 2021; Li et al. 2018; Wang et al. 2019, 2022a,b, 2025a), and the extended magnetic field decomposition method also makes the time-evolving MHD coronal model more stable numerically (Wang et al. 2025d). It remains a non-trivial task to implement the extended MHD decomposition strategy in an implicit global MHD coronal model from scratch, however. Alternatively, adopting a simplified model in the low coronal region, such as solving 1D equations for plasma motion below 1.1 R_s (Sokolov et al. 2021) or using a magneto-frictional model to drive the coronal evolution below 1.15 R_s, can also improve the numerical stability, if at the cost of compromising physical fidelity. In addition, an artificially broadened transition region through an increased plasma density can improve the numerical stability as well (Lionello et al. 2008; Mikić et al. 2013, 2018; Mok et al. 2005).

The methods mentioned above might enhance the numerical stability of the coronal model COCONUT. As a newly developed advanced MHD coronal model that benefits from an unstructured mesh and a fully implicit scheme, COCONUT provides an excellent platform for developing novel positivity-preserving (PP) approaches and for incorporating better established techniques. The time-evolving COCONUT (Wang et al. 2025b) has already employed a PP procedure for the plasma density and adopted a mass-flux limitation strategy (Hayashi 2005; Yang et al. 2012) to constrain the inner boundary plasma velocity that does not significantly exceed the characteristic speeds of supergranulation or sunspots. This paper further introduces a PP procedure for the thermal pressure to improve the capability of COCONUT to address coronal simulations of solar maxima.

During solar maxima, the solar activity evolves significantly more rapidly than during solar minima. This period is characterised by the frequent emergence of ARs, which are visible as sunspots and result from the shearing of magnetic fields by differential rotation and the subsequent rise of new magnetic flux through the convection zone (Brun & Browning 2017; Finley et al. 2024). Additionally, predominantly unipolar magnetic fields near the solar poles undergo a reversal in polarity during solar maxima (Brun et al. 2015). In addition, ARs emerging on the far side of the Sun, which cannot be captured well by single synoptic maps, have significant global effects on the magnetic structure (Perri et al. 2024). Wang et al. (2025b) has demonstrated that time-evolving MHD coronal simulations can be performed efficiently and accurately using an implicit method, thereby offering a more realistic alternative to quasi-steady-state coronal simulations for solar minimum conditions. The next logical step is to compare quasi-steady-state and time-evolving coronal simulations for solar maximum. This comparison illustrates the necessity and applicability of time-evolving MHD coronal simulations for practical space-weather forecasting during solar maximum.

At the onset of solar maximum in cycle 25, the strongest geomagnetic storm since November 2003 was triggered by NOAA AR 13664, the source of numerous CMEs and flares. Visible from Earth between May 2 and 14, 2024, AR 13664 evolved into one of the largest and most flare-productive ARs in recent decades. While these events have attracted considerable interest within the scientific community (Hayakawa et al. 2025; Jarolim et al. 2024; Kwak et al. 2024; Liu et al. 2024; Nedal et al. 2025), quasi-realistic MHD coronal simulations capturing the 3D evolution of global coronal structures during this period remain scarce. To address this gap, we performed coronal simulations that span two Carrington rotations (CRs) around the solar storms in May 2024 to validate the capability of COCONUT to simulate time-evolving coronal structures during solar maximum.

Based on the above considerations, the paper is organised as follows. In Section 2 we introduce the numerical formulation of the time-evolving coronal model and describe the PP measures we used to enhance its numerical stability. In Section 3 we present the simulation results calculated with COCONUT, including the evolution of the corona during two solar maximum CR periods around the solar storms in May 2024. We compare the simulation driven by an artificially rotating static magnetogram with that driven by a sequence of hourly updated time-evolving magnetograms. We also examine the effects of different grid resolutions. Finally, in Section 4, we summarise the key features of solar maximum coronal modelling using the fully implicit time-evolving coronal model, and we conclude.

2. Numerical algorithm

This section mainly describes the governing equations, the grid system, and the inner boundary conditions. It also presents the PP approach applied to the thermal pressure and plasma density, which is used to enhance the numerical stability of the model.

2.1. The governing equations and the grid system

Following Wang et al. (2025b,d), we performed time-evolving coronal simulations driven by a series of evolving magnetograms by solving the thermodynamic MHD equations in an inertial coordinate system. The governing equation is the same as in Wang et al. (2025b) and is described as follows:

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

Here, t is the time, U represents the conservative variable vector, ∇U denotes the spatial derivative of U, F(U) is the inviscid flux vector, and S(U,∇U) = S_gra + S_heat represents the source term vector corresponding to the gravitational force and the heating source terms. The radiative loss term Q_rad, the coronal heating term Q_H, and the thermal conduction term −∇⋅q are included in the heating source term S_heat.

Q_rad and Q_H are defined in the same way as in Wang et al. (2025b). The term −∇⋅q was computed using the Green-Gauss theorem (Wang et al. 2025b), and the formulation of the heat flux q is essentially the same as that used in the solar interplanetary phenomena implicit finite volume method (SIP-IFVM) coronal model (Wang et al. 2025d). The key difference is that in COCONUT, the heat flux is computed during a loop over faces, as a diffusive flux across the interface shared by cells i and j, and its gradient is calculated using the deferred correction approach described by Alvarez Laguna et al. (2016). For SIP-IFVM, the interface gradient in the heat flux is calculated as an averaged gradient combined with a normal direction finite difference correction, as described by Wang et al. (2022a). Additionally, the coronal model COCONUT uses unstructured meshes composed of prismatic elements, while SIP-IFVM uses a six-component composite structured mesh with hexahedral elements.

The Godunov method was adopted to advance the cell-averaged solution in time by solving a Riemann problem at each cell interface (Einfeldt et al. 1991; Godunov 1959). The computational domain was a spherical shell spanning from 1.01 to around 25 R_s. It was discretised into unstructured fifth-, sixth-, and seventh-level subdivided geodesic meshes (Brchnelova et al. 2022b), which includes 378 880, 1 515 520, and 6 062 080 non-overlapping prism cells (each consisting of two triangular and three quadrilateral faces), respectively. There are 74 layers of gradually stretched cells in the radial direction, each containing 5120, 20 480, and 81 920 cells. The finest mesh resolution is approximately 1°, which corresponds to roughly 100 minutes for the Sun to rotate through this angle. We therefore adopted a time step of 5 minutes in all simulations. This is sufficient to capture the evolution of the main coronal structures (Wang et al. 2025b).

2.2. Implementation of the boundary conditions

All time-dependent coronal simulations began with a steady-state solution constrained by a fixed magnetogram. In the steady-state and time-dependent coronal simulations, the observed GONG-zqs photospheric magnetograms provided the inner-boundary magnetic field (Li et al. 2021; Perri et al. 2023). We employed a potential field solver of 25th-order spherical harmonics expansion to extrapolate the photospheric magnetograms, where magnetic field strengths can reach several hundred Gauss near ARs, to the bottom of the low corona (Kuźma et al. 2023; Perri et al. 2022, 2023). The spherical harmonics expansion 𝒰 is originally expressed as

$$\begin{array}{c}\hfill \mathcal{U}\approx {\displaystyle \sum _{l=0}^{l_{\text{max}}}\sum _{m=-l}^l}{I}_l^m{Y}_l^m(\mu ,\varphi ),\end{array}$$(2)

where Y_l^m denotes the spherical harmonic function of degree l and order m, I_l^m represents the corresponding spectral coefficient, ϕ ∈ [0, 2π) is the longitude, and μ ≡ cos θ, with θ ∈ [0, π] denoting the colatitude.

The magnetic field strengths still exceed 50 Gauss and can lead to extremely low plasma β, which results in unphysical negative temperatures or pressures and ultimately causes simulation failures. We therefore applied a filter (McClarren & Hauck 2010) to 𝒰. The filtered spherical harmonics expansion, 𝒰_filtered, we used in the potential field reconstruction is defined as

$$\begin{array}{c}\hfill {\mathcal{U}}_{\mathrm{filtered}}\approx {\displaystyle \sum _{l=0}^{l_{\text{max}}}\sum _{m=-l}^l}\frac{I_l^m{Y}_l^m(\mu ,\varphi )}{1+\xi l^2{(l+1)}^2},\end{array}$$(3)

where a filter strength ξ = 1 × 10⁻⁴ was adopted. This approach effectively limited the maximum magnetic field near the solar surface to approximately 30 Gauss during the simulation.

The temperature and plasma density at the solar surface, denoted by T_s and ρ_s, were set to 1.8 × 10⁶ K and 3.34 × 10⁻¹³ kg m⁻³. Correspondingly, the thermal pressure at the inner boundary was calculated as p_s = 0.01 Pa. The velocity vector, the tangential magnetic field at the inner boundary, and the outer boundary conditions were treated in the same way as in Wang et al. (2025b).

In the time-evolving coronal simulations, we drove the model using a series of hourly updated GONG-zqs magnetograms², which have corrections at the poles to better estimate the global magnetic flux (Perri et al. 2023). The simulations were performed in a quasi-inertial coordinate system, with the Earth permanently positioned at ϕ = 60°. For the quasi-steady-state coronal simulation, we rotated a static magnetogram to generate a series of hourly updated magnetograms, which were then used to drive a rigidly rotating coronal structure. In the quasi-steady-state and time-evolving simulations, the inner boundary magnetic field at each time step was determined by applying a cubic Hermite interpolation to the four nearest hourly updated input magnetograms.

Fig. 1. Distribution of the radial magnetic field used as the inner boundary condition at the solar surface, shown in a co-rotating coordinate system.

2.3. Positivity-preserving measures

Under the presence of low-β regions (β < 10⁻³) and rapid variation in the magnetic field at the inner boundary during the time-evolving coronal simulations, an unphysical negative thermal pressure can easily occur during the calculation when deriving the thermal pressure from energy density. Following our previous work (Wang et al. 2025b,d), we took the PP procedure for the thermal pressure and plasma density in the MHD simulations within the entire computational domain.

We defined the inner boundary plasma density ρ_BC and the thermal pressure p_BC as follows (Brchnelova et al. 2023b; Wang et al. 2025b,d):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\rho }_{\mathit{BC}}=& {\mathrm{{\rm Y}}}_{\rho }\frac{{\mathbf{B}}_{\mathbf{BC}}^{\mathbf{2}}}{V_{A,BC\max }^2}+(1-{\mathrm{{\rm Y}}}_{\rho }){\rho }_s\hfill \\ \hfill p_{\mathit{BC}}=& {\mathrm{{\rm Y}}}_p\frac{{\mathbf{B}}_{\mathbf{BC}}^{\mathbf{2}}}{2}{\beta }_{min}+(1-{\mathrm{{\rm Y}}}_p)p_s\hfill \end{array},\end{array}$$(4)

where ${\mathrm{{\rm Y}}}_{\rho }=0.5+0.5·tanh(\frac{V_A-V_{A,BC\max }}{V_{\mathit{fac}}}·\pi )$ with $V_A=\frac{\left|{\mathbf{B}}_{\mathit{BC}}\right|}{{\rho }_s^{0.5}}$, $V_{A,BC\max }=3000{\mathrm{Km}\mathrm{S}}^{-1}$, and $V_{\mathit{fac}}=2{\mathrm{Km}\mathrm{S}}^{-1}$. For convenience of description, the magnetic field was already divided by $\sqrt{{\mu }_0}$, with μ₀ = 4 × 10⁻⁷π H m⁻¹ denoting the magnetic permeability. ${\mathrm{{\rm Y}}}_p=0.5+0.5·tanh(\frac{{\beta }_{min}-\frac{p_s}{0.5·{\mathbf{B}}_{\mathit{BC}}^2}}{{\beta }_{\mathit{fac}}}·\pi )$ with β_fac = 2 × 10⁻⁶ and β_min = 10⁻³.

During time-dependent simulations, as reported in Wang et al. (2025b), we used the 2nd-order accurate backward Euler (BDF2) method to calculate the temporal integration, and Newton iterations were performed within each time step to update the intermediate solution states. When the updated plasma density or thermal pressure fell beyond a threshold, we adopted 10⁻⁸ρ_s < ρ < 10 ρ_s for the plasma density and 10⁻⁸p_s < p < 100 p_s for the thermal pressure, we reverted to the most recent valid value. This is different to previous works. These thresholds were inspired by the PP reconstruction method by Feng et al. (2021) and Wang et al. (2022a). They limit the reconstructed density and pressure within this permitted interval. Furthermore, we also appropriately adjusted the thermal pressure in addition to the plasma density, as was done in Wang et al. (2025b), during each Newton iteration, as described below,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill p& ={\mathrm{{\rm Y}}}_p\frac{{\mathbf{B}}^{\mathbf{2}}}{2}{\beta }_{min}+(1-{\mathrm{{\rm Y}}}_p)p_o\hfill \\ \hfill \rho & ={\mathrm{{\rm Y}}}_{\rho }\frac{{\mathbf{B}}^{\mathbf{2}}}{V_{A,\max }^2}+(1-{\mathrm{{\rm Y}}}_{\rho }){\rho }_o\hfill \end{array},\end{array}$$(5)

where p_o and ρ_o are the originally updated thermal pressure and plasma density, and $V_{A,\max }=2\frac{{\left|\mathbf{B}\right|}_{\max }}{{\rho }_s^{0.5}}$ with ${\left|\mathbf{B}\right|}_{\max }=\underset{\forall \mathrm{cells}}{max}\left|\mathbf{B}\right|$. Additionally, inspired by Wang et al. (2025d), we constrained the plasma velocity in the range of 1 R_s ≤ r ≤ 1.1 R_s to avoid exceeding the speed of sound $C_s=\frac{\gamma ·p_i}{{\rho }_i}$, as described below,

$$\begin{array}{c}\hfill \mathbf{v}={\mathbf{v}}_o·min(\frac{C_s}{|{\mathbf{v}}_o|}·(0.3+tanh(\frac{r-R_s}{R_s}·8.68)),1.0),\end{array}$$(6)

where v_o is the originally updated velocity in the Newton iteration.

3. Numerical results

In this section, we employ the time-evolving coronal model with enhanced numerical stability to simulate the evolution of coronal structures during CRs 2283 and 2284, which cover the period of the solar storms in May 2024. As mentioned in Section 2.1, the simulations were performed in a quasi-inertial coordinate system with Earth fixed at ϕ = 60°, adopting a 5-minute time step on the fifth-, sixth-, and seventh-level geodesic meshes. Approximately 1300 hourly updated GONG-zqs magnetograms from 5:00 on April 9, 2024, to 16:00 on June 2, 2024, were employed to drive these simulations. Figure 1 presents the inner boundary magnetic field distribution at four different moments in a co-rotating coordinate system.

We performed all calculations on the WICE or Genius cluster, which is part of the Tier-1 and Tier-2 supercomputer infrastructure of Vlaams Supercomputer Centrum (VSC)³. Using 900 CPU cores on the Tier-2 infrastructure and a time step of 5 minutes, the time-evolving coronal simulations ran approximately 48 times and 9 times faster than real-time coronal evolution on the sixth- and seventh-level meshes, respectively. For the fifth-level mesh, the simulation can be approximately 60 times faster than real-time evolution using only 270 CPU cores on Tier-2. In the following subsections, we present the results of the MHD simulations for CRs 2283 and 2284.

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

In this subsection, we compare the results of the time-evolving simulation, driven by a series of hourly updated magnetograms, with those from the quasi-steady-state simulation. As described in Section 2.2, the latter was constrained by a rigidly rotating magnetogram at 11:00 on May 6, 2024, a moment marking the end of CR 2283 and the beginning of CR 2284 and during the solar storms in May 2024, which corresponds to the 654th hour of the time-evolving simulation. All the simulation results presented in this section were calculated on the sixth-level mesh. It shows that the temporal evolution of the magnetic field can lead to pronounced differences in the simulation results.

Fig. 2. Distributions of open- and closed-field regions derived from the time-evolving (contours) and quasi-steady-state (orange lines) simulation results. All shown in a co-rotating coordinate system, and the white and black patches represent open-field regions with magnetic field lines pointing outward and inward relative to the Sun, respectively. The grey patches indicate closed-field regions. The solid orange lines overlaid on these contours denote the edge of close-field regions derived from the quasi-steady-state simulation results.

Coronal holes (CHs) are typically associated with low plasma density and magnetic field lines that are open to interplanetary space and are among the most prominent features of the solar corona (Petrie et al. 2011; Feng et al. 2015, 2017; Wang et al. 2025d). Polar CHs are located at the solar poles and often extend to lower latitudes. They even occasionally cross the solar equator. Isolated CHs are detached from polar CHs, scattered across low- and mid-latitudes, and are commonly observed near solar maxima. Transient CHs are associated with solar eruptive events, such as CMEs, solar flares, and eruptive prominences. In Fig. 2 we trace magnetic field lines to the solar surface to distinguish open- and closed-field regions. This illustrates the pronounced changes in the CHs during a solar maximum CR period.

Compared with the magnetic field distribution between the 258th and 570th hour in Fig. 1, the emergence of positive magnetic field regions, with a polarity opposite to the background, around (75° N, 175° ) and (18° S, 348° ), as well as the disappearance of these regions around (10° N, 7° ), are well captured in the simulation results at the 570th hour. Although the region with an opposite magnetic polarity to the background around (75° N, 175° ) becomes very small by the 876th hour, the time-evolving coronal model still captures it. Additionally, the small-sized dipole that appears around (18° N, 350° ) at the 876th hour and the dipole around (26° N, 37° ) at the 1158th hour are also well captured in the simulation results. These features correspond to a closed-field region within the isolated CH around (26° S, 335° ) at the 876th hour and to the disappearance of the upper left portion of the extended CH around (10° N, 40° ) at the 1158th hour.

Fig. 3. White-light pB images observed from COR2/STEREO-A (top) and synthesised from quasi-steady-state (middle) and time-evolving (bottom) coronal simulation results ranging from 2.5 to 15 Rs on the meridian planes in the STEREO-A view. The orange lines highlight magnetic field lines on these selected meridional planes.

White-light polarised brightness (pB) images can reveal various large-scale coronal structures. High-density features, such as bipolar and pseudo-streamers, typically appear as bright regions in pB images. In contrast, low-density structures, such as CHs, manifest as dark regions (Feng et al. 2015, 2017, 2019; Feng 2020). In Fig. 3 we compare white-light pB images synthesised from the time-evolving (middle) and quasi-steady-state (bottom) simulation results with observed pB images from the innermost coronagraph of the Sun Earth Connection Coronal and Heliospheric Investigation (SECCHI) instrument suite (top) on board the Solar Terrestrial Relations Observatory Ahead (STEREO-A) spacecraft⁴ (Howard et al. 2008).

Figure 3 shows that at the 258th hour, the time-evolving simulation successfully captures the bipolar streamers centred around 37° N, 2° N, and 47° S at the west limb, as well as the bipolar streamer centred around 14° N at the east limb. The quasi-steady-state simulation fails to reproduce the bipolar streamer centred around 2° N at the west limb, however. The observed bright structures near the south pole and centred around 30° S at the east limb are misrepresented in both simulations as a broad bright feature, approximately 40° wide, centred around 60° S at the east limb. At the 570th hour, both simulations produce four bipolar streamers centred around 36° N, 7° S, and 34° S at the east limb and centred around 9° N at the west limb. Additionally, the four simulated bright structures, one at the west limb, two near the two poles, and the bipolar streamer centred around 7° S at the east limb, capture the observed bright structures. At the 876th hour, the time-evolving simulation captures two bipolar streamers centred around 35° N and 49° S at the west limb, as well as one centred around 4° N at the east limb. The quasi-steady-state simulation misses the bipolar streamer centred around 49° S at the west limb, however, and the one at the east limb is displaced significantly northward relative to the observations. Additionally, the quasi-steady-state simulation produces a pseudo-streamer centred around 18° S at the west limb, which is absent in the time-evolving simulation. Nevertheless, both simulations generally reproduce two broad bright structures near the two solar poles, which agrees with the observed features. At the 1158th hour, the bright structure centred around 56° S at the west limb in both simulations is slightly displaced southward in both simulations compared to the observations. The structure centred around 53° S at the east limb is reasonably well reproduced. Neither simulation captures the other three observed bright structures, however, and the broad bright feature around the north pole deviates significantly from what is observed.

This comparison indicates that the simulated plasma density is generally consistent with the observed bright structures, and the time-evolving simulation agrees better than the quasi-steady simulation. The discrepancies in the width and position of these bright structures between observations and simulation may come from the imperfect measurements of photospheric magnetic fields, the lack of direct observation of the polar magnetic field, and the limitation of the synoptic magnetograms in which the magnetic field at different longitudes is observed at different times (Hamada et al. 2018).

Fig. 4. Synoptic maps of east (left) and west limb (right) white-light pB observations from the SOHO instrument LASCO C2 at 3 Rs for CRs 2284 and 2283, respectively. The dashed orange line and solid grey line represent the corresponding MNLs derived from the time-evolving and quasi-steady-state simulations, respectively.

In Fig. 4 we further present synoptic maps of the east- (left) and west-limb (right) observations from the Large Angle and Spectrometric Coronagraph C2 (LASCO-C2) (Brueckner et al. 1995) on board the Solar and Heliospheric Observatory (SOHO)⁵ for CRs 2284 and 2283, respectively. The west- and east-limb white-light pB images for CR 2283 and CR 2284 were synthesised from observations taken between April 16, 2024, and May 13, 2024, and between April 29, 2024, and May 26, 2024, respectively. These correspond to intervals from 162nd to 816th hours and 492nd to 1142nd hours in the time-evolving simulations.

Figure 4 shows that in the west-limb pB image, the magnetic neutral line (MNL) in the time-evolving and quasi-steady-state simulations coincide with the observed bright structures from 10° to 180° and from 210° to 240° in longitude. The MNL from the time-evolving simulation is more consistent with these observed bright structures from 10° to 40° and from 315° to 360°, however. In the east-limb image, the MNL in the time-evolving simulation extends farther north around 10° and 170° in longitude, crosses the region near 5° N from 90° to 115°, and vertically spans from 10° N to 50° S in latitude. This matches the observations better than the quasi-steady-state simulation.

In Fig. 5 we present the 2D timing diagrams of plasma density, velocity, and temperature at 21.5 R_s, derived from the time-evolving simulation results during CRs 2283 and 2284. These images are plotted along the latitudes intersected by the Sun–Earth line. The figure shows that the MNL generally coincides with the latitude distribution of the high-density low-speed solar wind. The distribution of the simulated plasma temperature is positively correlated with the radial speeds of the solar wind (Elliott et al. 2012; Pinto & Rouillard 2017; Li et al. 2018), with the fast solar wind (V_r > 550 km s⁻¹) exhibiting temperatures ranging from 1.1 to 1.8 MK, while the slow wind (V_r < 450 km s⁻¹) corresponds to temperatures between 0.4 and 0.8 MK. Between the 260th and 320th hours, around the 930th hour, and between the 1240th and 1280th hours, however, the MNL is accompanied by low-density high-speed flows. This may be attributed to the absence of a self-consistent heating mechanism in this model.

Fig. 5. Timing diagrams of simulated radial velocity Vr (km s−1, top), plasma number density (103 cm−3, middle), and temperature (105 K, bottom) along the latitudes intersected by the Sun–Earth line at 0.1 AU. The dashed orange and solid grey lines represent the MNLs derived from the time-evolving and quasi-steady-state simulations, respectively.

Additionally, the evolution of MNLs in the time-evolving and quasi-steady-state simulations closely coincides between the 300th and 650th hours, but drifts apart thereafter. This agreement is attributed to the fact that the magnetic fields at latitudes west of the one intersected by the Sun–Earth line at the 654th hour, the reference time for the synoptic magnetogram used in the quasi-steady-state simulation, are gradually updated during this time interval (Hamada et al. 2018) around the 300th and 650th hours. In contrast, the subsequent deviation is due to the absence of the magnetic field evolution beyond this interval in the synoptic magnetogram used in the quasi-steady-state simulation.

3.2. The impact of the grid resolutions on the simulation results

In this subsection, we further evaluate the impact of the grid cell sizes on the time-evolving coronal simulation results. When the cell size is reduced in the tangential direction to one-quarter of that in the fifth-level grid mesh, which yields a resolution that is slightly higher than that of the synchronised magnetograms, the simulated magnetic flux at 0.1 AU increases by more than 40%. The sixth-level grid mesh achieves the desired accuracy in capturing large-scale coronal structures while maintaining high computational efficiency and numerical stability for global coronal simulations near solar maximum. The simulation results were first interpolated onto a structured mesh with a resolution of 150 × 300 cells in the tangential direction to enable a comparison across unstructured meshes with different resolutions. The interpolation was performed using the radial basis function method proposed by Wang et al. (2022a), which is numerically stable, spatially accurate, and computationally efficient. All comparisons were conducted using the interpolated solutions on the structured mesh.

Fig. 6. Contours of open- and closed-field regions at the 876th hour of the time-evolving simulations, performed on the seventh- (left) and fifth-level (right) subdivided geodesic meshes, respectively. The orange lines overlaid on these contours denote the edge of close-field regions derived from the corresponding result on the sixth-level mesh.

In Fig. 6 we display the distribution of open- and closed-field regions derived from the 876th hours of the time-evolving simulation results on the fifth- and seventh-level subdivided geodesic meshes, respectively. The corresponding edges of the simulated open- and closed-field regions on the sixth-level mesh are overlaid. The polar CH regions on the fifth- and seventh-level meshes coincide with those identified on the sixth-level mesh. The seventh-level mesh in addition captures open-field regions at the lower left end of three isolated CHs centred around (21° S, 87° ), (19° S, 245° ) and (29° S, 335° ) and at the northern edge of the extended CH near (10° N, 100° ) and around (50° S, 265° ), which are missing on the sixth-level mesh. In contrast, the open-field regions around these positions and at the lower left end of the isolated CH centred at (7° N, 305° ) appear on the sixth-level mesh but are absent from the fifth-level mesh. The minor bipolar structure around (70° N, 160° ) that is shown in Fig. 1 is well captured on sixth- and seventh-level meshes, but is almost missed on the fifth-level mesh.

Fig. 7. Timing diagram of the radial velocity Vr (km s−1; a, d), proton number density (105 cm−3 at 3 Rs and 103 cm−3 at 21.5 Rs; b, e), and magnetic field strength |B| (G; c, f) observed by two virtual satellites located at 3 Rs (a, b, c) and 21.5 Rs (d, e, f), respectively. These virtual satellites are positioned at the same latitude as Earth, but lag by 60° in longitude. The solid and dashed black lines and diamonds represent the time-evolving simulation results on the fifth-, sixth-, and seventh-level subdivided geodesic meshes, respectively, and the solid grey line indicates the radial velocity observed by the WIND satellite.

In Fig. 7 we present the temporal evolution of radial velocity V_r (km s⁻¹; a, d), proton number density N (10⁵ cm⁻³ and 10³ cm⁻³; b, e), and magnetic field strength |B| (G or nT; c, f) at 3 (a, b, c) and 21.5 R_s (d, e, f). These parameters were monitored at the same solar latitude as Earth, but at a longitude that lags by 60°, assuming that the solar wind takes approximately 100 hours to travel from 0.1 to 1 AU (Wang et al. 2025d). The radial velocity observed by the WIND satellite⁶ (King & Papitashvili 2005) at 1 AU is also presented. It reveals that a peak or trough in the plasma density profile generally corresponds to a trough or peak in the velocity profile in the three meshes. Moreover, the appearance times of the peak and trough in plasma density, radial velocity, and magnetic field strength coincided in different meshes. The magnetic field strength and radial velocity are typically higher in finer meshes, however, whereas the plasma density tends to be lower. The values of the radial velocity and plasma density in the seventh-level mesh closely match those in the sixth-level mesh, except in the vicinity of peaks and troughs.

Additionally, the evolution pattern of the radial velocity at 0.1 AU in the time-evolving simulations between the 400th and 1050th hours generally captures the observed radial velocity at 1 AU. The simulated velocity between the 400th and 550th hours is faster than the observations, however, while the simulated velocity peak between the 760th and 820th hours is noticeably slower. The simulated bi-peak in velocity around the 170th hour may correspond to the observed bi-peak near the 290th hour. The simulation at 0.1 AU fails to reproduce the observed low-speed flow at 1 AU between the 1100th and 1300th hours, however. These discrepancies may be attributed to several factors: the limitations of the empirical heating function used to approximate coronal heating and solar wind acceleration; inaccuracies in the synoptic magnetograms, where magnetic fields at different longitudes are observed at various times and polar magnetic fields are poorly measured, resulting in an incomplete representation of the actual magnetic field distribution; and the absence of transient phenomena such as CMEs, which are not included in the model, but can significantly affect solar wind dynamics. Furthermore, pronounced magnetic flux emergence and cancellation during solar maximum can produce complex solar wind structures that interact as they propagate from 0.1 AU to 1 AU, which substantially alters the solar wind structures.

The radial velocity also increases with higher mesh resolution during most of the simulation period. This may be attributed to the heating source term that we adopted in the model, which is proportional to the magnetic field strength. This simplification does not fully capture the actual heating and acceleration mechanisms. The finer mesh resolution can capture more magnetic flux, which results in stronger magnetic fields at 0.1 AU and consequently higher radial velocities. Additionally, the simulated radial velocity generally exceeds the in situ observations, except during the period between the 600th and 850th hours. This discrepancy is more pronounced on the seventh-level mesh, which also shows improved consistency with observations around the velocity peak between the 760th and 820th hours. These findings suggest that while finer meshes capture more details of the coronal structures, they can also amplify the negative impacts of unphysical or overly simplified assumptions within the model.

In Table 1 we further list the average relative difference in solution variables between the results calculated on the Kth- and fifth-level meshes, denoted as ${\mathrm{RD}}_{\mathrm{ave},\chi }^{{\mathrm{LV}}_K}$,

$$\begin{array}{c}\hfill {\mathrm{RD}}_{\mathrm{ave},\chi }^{{\mathrm{LV}}_K}=\frac{1}{N^{{\mathrm{LV}}_K}}\sum _{i=1}^{N^{{\mathrm{LV}}_K}}\frac{{\chi }^{{\mathrm{LV}}_K}-{\chi }^{{\mathrm{LV}}_5}}{{\chi }^{{\mathrm{LV}}_5}}·\end{array}$$

Here, χ ∈ {|B|,ρ, V_r} means the corresponding variable, and N^LV_K denotes the number of sample points selected from the simulation results on the N_th^LV_K level mesh, as shown in Fig. 7. For the fifth- and sixth-level meshes at 3 R_s and 21.5 R_s, respectively, we selected 218 sample points with a 6-hour interval between two adjacent points. Because the seventh-level mesh output files are extremely large (∼1.7 GB each), we only selected 65 sample points with a 20-hour cadence. The subscript LV_K refers to the corresponding variables in the Kth-level mesh. The simulated magnetic field strength from the sixth- and seventh-level meshes is stronger by more than 20% and 40%, respectively, than the field calculated on the fifth-level mesh at 0.1 AU. This enhancement in the magnetic field strength on finer meshes is also more pronounced at 21.5 R_s than at 3 R_s. Accordingly, the radial velocity increases and the plasma density decreases in the finer mesh, although these fluctuations are less pronounced than the enhancement in the magnetic field strength.

These results indicate that the refined mesh can capture open field regions that are missed by the coarser mesh, which increases the simulated magnetic flux away from the solar surface noticeably. This, in turn, affects the simulated flow field and typically results in an increase in radial velocity and a decrease in plasma density, although the magnitude of the average relative difference is lower by approximately an order of magnitude than that observed in the magnetic field strength.

4. Concluding remarks

Recently, Wang et al. (2025b) established the fully implicit time-evolving MHD coronal model for unstructured grids. It performs more than 70 times faster than real-time global coronal evolutions during a solar minimum (using 1080 CPU cores for 1.5 M grid cells). We further improved its PP property here to enable the efficient computation of time-evolving coronal evolution during solar maximum, with magnetic field strengths that did not significantly exceed 30 Gauss. We drove this model by about 1300 hourly updated GONG-zqs magnetograms to simulate coronal evolutions during two CRs around the solar storms in May 2024.

A comparison between the time-evolving coronal simulation and a quasi-steady-state simulation driven by a rigidly rotating magnetogram demonstrated that the temporal evolution of the magnetic field during a solar maximum CR can lead to significant differences in simulation results. Magnetic flux emergence and cancellation result in noticeable shifts, disappearances, and appearances of high-density streamers in white-light pB images. The MNL in the time-evolving simulation agrees more closely with the observed bright structures in the synoptic maps of pB images, and the evolution of magnetic field distributions on the solar surface are well captured by the simulated time-evolving 3D magnetic field structures. Additionally, the time-evolving simulation produces pronounced differences in the MNL at 0.1 AU compared to the quasi-steady-state simulation. This shows that the evolution of the inner boundary magnetic field within a CR can lead to significant changes in the latitudinal structures of the solar wind at 0.1 AU.

Furthermore, a comparison between the time-evolving simulation results on the fifth, sixth-, and seventh-level meshes demonstrated that mesh refinement enables the model to capture the missed open field regions in coarser meshes. When the resolution is increased by a factor of 16, from the fifth- to the seventh-level mesh, the simulated magnetic flux away from the solar surface increases by more than 40%. Although the changes in the radial velocity and plasma density are less pronounced than those in the magnetic flux, they exhibit corresponding increases and decreases, respectively. This comparison, together with the figures of the evolution of open and closed field regions on the sixth-level mesh, also showed that the sixth-level mesh can capture the evolution of small-sized dipoles well, and also the emergence of small-sized magnetic field regions with a polarity opposite to that of the background.

Additionally, given that both the finest mesh resolution and the distance between adjacent magnetic field pixels in the observed magnetograms are approximately 1°, corresponding to about 100 minutes of solar rotation, a time step of 5 minutes is sufficient to capture the evolution of the coronal structures with an adequate temporal resolution in global coronal simulations. Considering that this model can perform faster by about 48 and 9 times than real-time coronal evolution using 900 CPU cores on the sixth- and seventh-level meshes, respectively, we conclude that a time step of 5 minutes and the sixth-level mesh can generally balance computational efficiency, temporal accuracy, and numerical stability for global coronal simulations near solar maximum. This fully implicit time-evolving MHD coronal model with improved PP property is promising for the timely and accurate simulation of the time-evolving coronal structures in the daily practical space-weather forecasting around solarmaximum.

Several challenges remain to be addressed for further improvement of this novel time-evolving implicit MHD coronal model, however. The model can accurately capture the evolution of even small-scale dipolar structures observed in magnetograms through simulated 3D magnetic field evolution, and it is therefore essential to use synchronised magnetograms, where the magnetic fields across all longitudes are captured simultaneously (Upton & Hathaway 2014; Downs et al. 2025). This is necessary to overcome the limitations of current synoptic magnetograms, in which magnetic fields at different longitudes are observed at different times, which results in significant deviations from the true magnetic field distributions. Extremely finer mesh comparable to the ninth-level mesh, containing 16 times as much as grid cells on the seventh-level mesh each layer, is still required to resolve the much faster small-scale dynamics of ARs, such as sunspot rotation, to self-consistently simulate transient phenomena, including CME events that are triggered by sunspot rotation (Jiang et al. 2021), and to reduce diffusion due to numerical errors for capturing a sufficiently thin current sheet (Jiang & Zhang 2025) in global MHD coronal simulations. Correspondingly, more stable and efficient numerical algorithms are still needed to deal with extremely low-β issues after considering these AR evolutions. Additionally, physically consistent heating source terms are required to better simulate coronal heating and solar wind acceleration during time-evolving coronal simulations around solar maximum. Moreover, accurate measurements of the photospheric magnetic fields in the polar regions are necessary to reproduce more realistic coronal structures during solar maximum.

In our future work, we plan to simplify the novel extended MHD decomposition strategy proposed by Wang et al. (2025d) and implement it in the time-evolving coronal model COCONUT we established here. This improvement aims to further reduce the contamination of magnetic field discretisation errors on thermal pressure derived from the total energy, and consequently, to prevent the occurrence of unphysical negative thermal pressure values, which can cause the code to crash when dealing with low-β issues. Furthermore, we plan to implement a local mesh refinement that increases the resolution to that of the ninth-level mesh only in regions that contain AR evolutions. This approach will enable self-consistent faster-than-real-time simulations of CMEs triggered by sunspot rotation using the implicit time-evolving global MHD coronal model. Additionally, we plan to generate a series of synchronised magnetograms by incorporating surface flux transport models, such as the advective flux transport model (Upton & Hathaway 2014), with data from the Solar Orbiter Polarimetric and Helioseismic Imager (Loeschl et al. 2024), horizontal velocities inferred from observational data using time-distance helioseismology (Zhao et al. 2012; Yalim et al. 2017), and artificial intelligence-generated data from STEREO extreme-UV observations (Jeong et al. 2020, 2022). This method has the potential to generate a more realistic real-time magnetic field evolution at the inner boundary of our coronal model, which will lead to more accurate simulation results that match the observations around solar maximum. We also plan to couple this time-evolving coronal model with EUHFORIA, an inner heliospheric model, to further develop a quasi-realistic faster-than-real-time Sun-to-Earth modelling chain that is suitable for practical daily space-weather forecasting.

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 supported by the National Natural Science Foundation of China (grant No. 42030204) and the BK21 FOUR programme of the Graduate School, Kyung Hee University (GS-1-JO-NON-20242364). 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) programme, 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

All Tables

All Figures

	Fig. 1. Distribution of the radial magnetic field used as the inner boundary condition at the solar surface, shown in a co-rotating coordinate system.
In the text

Fig. 2. Distributions of open- and closed-field regions derived from the time-evolving (contours) and quasi-steady-state (orange lines) simulation results. All shown in a co-rotating coordinate system, and the white and black patches represent open-field regions with magnetic field lines pointing outward and inward relative to the Sun, respectively. The grey patches indicate closed-field regions. The solid orange lines overlaid on these contours denote the edge of close-field regions derived from the quasi-steady-state simulation results.

In the text

	Fig. 3. White-light pB images observed from COR2/STEREO-A (top) and synthesised from quasi-steady-state (middle) and time-evolving (bottom) coronal simulation results ranging from 2.5 to 15 Rs on the meridian planes in the STEREO-A view. The orange lines highlight magnetic field lines on these selected meridional planes.
In the text

	Fig. 4. Synoptic maps of east (left) and west limb (right) white-light pB observations from the SOHO instrument LASCO C2 at 3 Rs for CRs 2284 and 2283, respectively. The dashed orange line and solid grey line represent the corresponding MNLs derived from the time-evolving and quasi-steady-state simulations, respectively.
In the text

	Fig. 5. Timing diagrams of simulated radial velocity Vr (km s−1, top), plasma number density (103 cm−3, middle), and temperature (105 K, bottom) along the latitudes intersected by the Sun–Earth line at 0.1 AU. The dashed orange and solid grey lines represent the MNLs derived from the time-evolving and quasi-steady-state simulations, respectively.
In the text

	Fig. 6. Contours of open- and closed-field regions at the 876th hour of the time-evolving simulations, performed on the seventh- (left) and fifth-level (right) subdivided geodesic meshes, respectively. The orange lines overlaid on these contours denote the edge of close-field regions derived from the corresponding result on the sixth-level mesh.
In the text

Fig. 7. Timing diagram of the radial velocity Vr (km s−1; a, d), proton number density (105 cm−3 at 3 Rs and 103 cm−3 at 21.5 Rs; b, e), and magnetic field strength |B| (G; c, f) observed by two virtual satellites located at 3 Rs (a, b, c) and 21.5 Rs (d, e, f), respectively. These virtual satellites are positioned at the same latitude as Earth, but lag by 60° in longitude. The solid and dashed black lines and diamonds represent the time-evolving simulation results on the fifth-, sixth-, and seventh-level subdivided geodesic meshes, respectively, and the solid grey line indicates the radial velocity observed by the WIND satellite.

In the text