COCONUT: A coronal model with an energy decomposition strategy

TL;DR

This work tackles the numerical instability that arises in low-beta coronal MHD simulations when the thermal pressure is derived from the full energy equation. It introduces an energy-decomposition strategy that evolves a decomposed energy E1 (internal plus kinetic energy) and pairs it with an HLL flux that includes an additional dissipation term to mitigate catastrophic cancellation, all within the time-evolving implicit MHD code COCONUT. The approach achieves near parity with the full-energy results during solar minimum and rising phases while markedly improving stability during solar maximum, demonstrated by quasi-steady tests and a solar-maximum time-evolution that remained stable where the full-energy version failed. The method enables more realistic coronal evolution modeling through solar cycles and opens avenues for data-driven CME modeling, with demonstrated scalability on high-performance computing resources and compatibility with existing boundary-conditions and PF extrapolation practices.

Abstract

In this paper, we propose an energy decomposition method combined with an HLL Riemann solver that includes an additional dissipation term in the energy equation to improve the numerical stability of the fully implicit, time-evolving coronal model COCONUT and extend its applicability to solar-maximum phases. In MHD simulations that evolve conservative variables in time, the thermal pressure is typically computed by subtracting the magnetic and kinetic energies from the total energy. In low-beta (the ratio of thermal to magnetic pressure; $< 10^{-3}$) regions, discretization errors of magnetic energy can be comparable to the thermal pressure, potentially leading to negative thermal pressure and causing the simulation to crash. Therefore, we update the decomposed energy, excluding the magnetic energy, at each time step. It avoids subtracting a large magnetic energy from the total energy to obtain a very small thermal pressure in low-$β$ regions, thereby improving the numerical stability of MHD models. We validate the algorithm using a time-evolving solar-maximum Carrington rotation simulation in 2025, which the previous code failed to run to completion. We also perform quasi-steady-state coronal simulations and 2D benchmark tests to further assess the algorithm's performance. The simulation results show that the algorithm produces results nearly identical to those obtained using the traditional full energy equation during solar minimum, while significantly improving COCONUT's ability to simulate coronal evolution under strong magnetic fields, even including fields exceeding 100 Gauss with $β<10^{-3}$. This method provides a promising approach for performing quasi-realistic coronal simulations during solar maxima.

Figures (11)

• Figure 1: Distribution of the radial plasma speed $V_r \, (\mathrm{km~s^{-1}})$ and proton number density $N \, (10^3~\mathrm{cm^{-3}})$ at 0.1 AU calculated using COCONUT with the decomposed energy equation (first row). The relative differences in radial velocity and plasma density between the results obtained with the decomposed and the traditional full-energy equations at 0.1 AU (second row) and 3 $R_s$ (third row) are also shown. The fourth row presents the white-light pB images observed by LASCO/SOHO and synthesized from electron density which is treated to be the same as the proton density simulated by COCONUT, spanning a range from 2.3 to 6$\;R_s$ with the decomposed energy equation. The black solid and orange dashed lines represent the MNLs calculated using the decomposed and full energy equations, respectively. The orange solid lines denote the magnetic field lines calculated with the decomposed energy equation. These figures correspond to CR 2073.
• Figure 2: Same as Figure \ref{['localdiffCR2073']}, but for CR 2248.
• Figure 3: Synoptic maps of east (top left) and west limb (top right) white-light pB observations from the LASCO C2 instrument onboard SOHO spacecraft at 3 $R_s$ for CRs 2073 and 2248, alongside the simulated electron number density $N \, (10^5~\mathrm{cm^{-3}})$ which is treated to be identical to the proton density
• Figure 4: Full-disk image from EIT/SOHO at 195 Å on 2008 August 3 (top left), the corresponding simulated full-disk image of open- and closed-field regions obtained with the decomposed energy equation (top middle), and the global distributions of open- and closed-field regions derived from COCONUT simulations with the decomposed energy equation for CR 2073 (top right) and CR 2248 (bottom right). The white and black patches denote open-field regions, where magnetic field lines point outward and inward, respectively, relative to the Sun. The grey regions represent closed-field regions. The orange solid lines overlaid on these contours denote the edge of open-field regions derived from simulation results with the full energy equation. The bottom-left panel shows the synoptic EUV map from the 193 Å channel of AIA/SDO for CR 2248.
• Figure 5: Distributions of the inner-boundary radial magnetic field at the 0th, 130th, 260th, 390th, 520th, and 650th hours of the time-evolving coronal simulation during CR 2296. White lines overlaid on the magnetic field contours indicate the boundaries of the open-field regions derived from COCONUT simulation results using the decomposed energy equation.