SIP-IFVM: A Time-evolving Coronal Model with an Extended Magnetic Field
Decomposition Strategy
Haopeng Wang
1aa, Liping Yang2aa, Stefaan Poedts1,3aa, Andrea Lani1,4aa, Yuhao Zhou1aa, Yuhang Gao1,5aa, Luis Linan1aa,
Jiakun Lv
6
, Tinatin Baratashvili
1aa, Jinhan Guo1,7aa, Rong Lin1,5aa, Zhan Su2,8, Caixia Li9, Man Zhang2aa, Wenwen Wei10aa,
Yun Yang
11aa, Yucong Li1,2,8, Xinyi Ma2,8, Edin Husidic1,12aa, Hyun-Jin Jeong1,13aa, Mahdi Naja5-Ziyazi1aa, Juan Wang14, and
Brigitte Schmieder
1,15,16aa
1
Centre for Mathematical Plasma-Astrophysics, Department of Mathematics, KU Leuven, Celestijnenlaan 200B, 3001 Leuven, Belgium;
haopeng.wang1@kuleuven.be
2
SIGMA Weather Group, State Key Laboratory of Space Weather, National Space Science Center, Chinese Academy of Sciences, Beijing 100190, People’s
Republic of China; lpyang@swl.ac.cn
3
Institute of Physics, University of Maria Curie-Skłodowska, ul. Radziszewskiego 10, 20-031 Lublin, Poland
4
Von Karman Institute for Fluid Dynamics, Waterloosesteenweg 72, 1640 Sint-Genesius-Rode, Brussels, Belgium
5
School of Earth and Space Sciences, Peking University, Beijing 100871, People’s Republic of China
6
Beijing Institute of Applied Meteorology, Beijing 100029, People’s Republic of China
7
School of Astronomy and Space Science and Key Laboratory of Modern Astronomy and Astrophysics, Nanjing University, Nanjing 210023, People’s Republic of
China
8
College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing 100049, People’s Republic of China
9
Shenzhen Key Laboratory of Ultraintense Laser and Advanced Material Technology, Center for Intense Laser Application Technology, and College of Engineering
Physics, Shenzhen Technology University, Shenzhen 518118, People’s Republic of China
10
Space Sciences Laboratory, University of California, Berkeley, CA 94720, USA
11
School of Mathematical Sciences, Ministry of Education Key Laboratory for NSLSCS, Nanjing Normal University, Nanjing 210023, People’s Republic of China
12
Department of Physics and Astronomy, University of Turku, 20014 Turku, Finland
13
School of Space Research, Kyung Hee University, Yongin, 17104, Republic of Korea
14
School of Systems Science, Beijing Normal University, Beijing 100875, People’s Republic of China
15
Observatoire de Paris, LIRA, UMR8254 (CNRS), F-92195 Meudon Principal Cedex, France
16
SUPA, School of Physics & Astronomy, University of Glasgow, Glasgow G12 8QQ, UK
Received 2025 March 17; revised 2025 April 21; accepted 2025 April 24; published 2025 June 10
Abstract
Time-evolving magnetohydrodynamic (MHD) coronal modeling, driven by a series of time-dependent
photospheric magnetograms, represents a new generation of coronal simulations. This approach offers more
realistic results than traditional steady coronal models constrained by a static magnetogram. However, its practical
application is signi5cantly limited by the low computational ef5ciency and poor numerical stability in solving
low-β issues common in coronal simulations. To address this, we propose an extended magnetic 5eld
decomposition strategy and successfully implement it in an implicit MHD coronal model. The traditional
decomposition strategies split the magnetic 5eld into a time-invariant potential 5eld and a time-dependent
component B1. This works well for quasi-steady-state coronal simulations where |B1| is typically small. However,
when the inner-boundary magnetic 5eld evolves, |B1| can grow signi5cantly, and its discretization errors often
lead to nonphysical negative thermal pressure, ultimately causing the simulation to crash. In the extended
magnetic 5eld decomposition strategy, we split the magnetic 5eld into a temporally piecewise-constant 5eld and a
time-varying component, B1. This effectively keeps |B1| consistently small throughout the simulations and
performs well in solving time-evolving low-β issues, thereby outperforming traditional methods. We incorporate
this improved strategy into our implicit MHD coronal model and apply it to simulate the evolution of coronal
structures within 0.1 au over two solar-maximum Carrington rotations. The results show that this coronal model
effectively captures observational features and performs more than 80 times faster than real-time evolutions using
only 192 CPU cores, making it well suited for practical applications in simulating the time-evolving corona.
Unied Astronomy Thesaurus concepts: The Sun (1693); Magnetohydrodynamics (1964); Magnetohydrodyna-
mical simulations (1966); Solar corona (1483)
Materials only available in the online version of record: animation
1. Introduction
Space weather can impact the performance and reliability of
both space-borne and ground-based technological systems and
pose risks to human health. To enable timely action in
mitigating damage from severe space weather events,
advanced Sun-to-Earth model chains should be developed
(e.g., X. S. Feng et al. 2011, 2013; J. Pomoell &
S. Poedts 2018; S. Poedts et al. 2020; K. Hayashi et al.
2021) to better understand space weather mechanisms and
provide reliable forecasts hours to days in advance (e.g.,
D. N. Baker 1998; X. S. Feng et al. 2013; H. E. J. Koskinen
et al. 2017; X. S. Feng 2020). In the Sun-to-Earth modeling
chain, observed photospheric magnetic 5elds serve as input for
the solar coronal model, providing inner boundary conditions
for the inner heliospheric model, which then supplies boundary
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June https://doi.org/10.3847/1538-4365/add0b1
© 2025. The Author(s). Published by the American Astronomical Society.
aaaaaaa
Original content from this work may be used under the terms
of the Creative Commons Attribution 4.0 licence. Any further
distribution of this work must maintain attribution to the author(s) and the title
of the work, journal citation and DOI.
1
information to the geomagnetic model. The coronal model is
essential for initializing other models and is key in accurately
simulating solar disturbances such as coronal mass ejections
(CMEs; M. Brchnelova et al. 2022; B. Kuźma et al. 2023;
B. Perri et al. 2023). However, physics-based MHD coronal
models are also the most complex and computationally
intensive component (H. P. Wang et al. 2025a, 2025b), and
sometimes encounter problems with low β (the ratio of the
thermal pressure to the magnetic pressure), where β can drop
as low as 10−4 near the solar surface (P.-A. Bourdin 2017),
leading to severe challenges to stability and ef5ciency
(X. S. Feng et al. 2021; H. P. Wang et al. 2022a).
In most MHD coronal simulations, time steps are limited to
a few seconds due to the restriction of the Courant–Friedrichs–
Lewy (CFL) stability condition. Consequently, depending on
the mesh resolution, even state-of-the-art quasi-steady-state
coronal simulations require 10,000–100,000 CPU-hours to
reach a quasi-steady state (X. S. Feng et al. 2019; V. Réville
et al. 2020). In contrast, for most MHD inner heliospheric
models, the CFL-limited time steps are typically of the order of
10 minutes (T. Detman et al. 2006; K. Hayashi 2012). As a
result, coronal models require signi5cantly greater computa-
tional resources to remain synchronized with inner helio-
spheric models. Although empirical solar coronal models
(C. N. Arge et al. 2003) offer higher ef5ciency, they discard
important information and fail to deliver the necessary
accuracy in forecasts (E. Samara et al. 2021). Therefore,
further efforts are needed to develop more ef5cient and
accurate MHD coronal models. Implicit temporal discretiza-
tion strategies can help overcome the limitations imposed by
the CFL condition, thereby enhancing computational ef5-
ciency by allowing longer time steps.
Recently, several successful attempts have been made to
increase the ef5ciency of MHD coronal models by using
implicit solvers (B. Perri et al. 2018, 2022, 2023; Y. Wang
et al. 2019; X. S. Feng et al. 2021; H. P. Wang et al.
2022a, 2022b, 2025a; B. Kuźma et al. 2023). In the implicit
algorithm, the convergence rate is improved by selecting a
considerable time step. Additionally, both the second-order
backward differentiation formula (BDF2) with Newton itera-
tions at each time step and the pseudo–time marching method,
which introduces a pseudotime at each physical time step, are
employed to improve temporal accuracy in implicit coronal
models (J. H. Guo et al. 2023; L. Linan et al. 2023; H. P. Wang
et al. 2025a). While these models have achieved the desired
speedup, they remain constrained by a time-invariant magne-
togram, which contrasts with the reality that the solar coronal
structure evolves (M. J. Owens et al. 2017). This discrepancy
leads to differences between simulation results and coronal
observations (M. D. Cash et al. 2015; V. Réville et al. 2020).
Many researchers have focused on developing time-evol-
ving coronal models to address the limitations of commonly
used quasi-steady-state coronal models, which do not account
for the evolution of coronal structures. These models, typically
driven by hundreds of time-varying observed photospheric
magnetograms, can capture the evolution of coronal structures
with higher 5delity (L. P. Yang et al. 2012; A. Yeates et al.
2018; X. S. Feng et al. 2023). They may also improve the
realism and reliability of solar wind and CME modeling
(R. Lionello et al. 2023). However, time-evolving MHD
coronal simulations remain prohibitively computationally
expensive (A. Yeates et al. 2018). One of the main reasons
is that most state-of-the-art time-evolving coronal models
(L. P. Yang et al. 2012; K. Hayashi et al. 2021; R. Lionello
et al. 2023) still rely on explicit or semi-implicit approaches,
where only speci5c source terms are treated implicitly while
the time step remains constrained by the explicitly treated
terms, leading to extremely low ef5ciency. As a result, explicit
or semi-implicit MHD coronal simulations typically require
thousands of computing cores to achieve faster-than-real-time
execution. A more detailed description of time-evolving
coronal models is given by H. P. Wang et al. (2025b).
Furthermore, H. P. Wang et al. (2025b) extended the quasi-
steady-state CoolLuid Coronal Unstructured (COCONUT)
model, a novel implicit MHD solar corona model based on
the Computational Object-Oriented Libraries for Fluid
Dynamics (COOLFluiD; D. Kimpe et al. 2005; A. Lani
et al. 2005, 2013),17 into a time-evolving coronal model. It is
the 5rst fully implicit time-evolving MHD coronal model,
allowing for long time steps exceeding the CFL condition. It
can simulate the evolution of a full Carrington rotation (CR)
period in just 9 hr of computational time (with a time step of
10 minutes using 1080 CPU cores and approximately 1.5
million cells). However, it currently struggles with resolving
low-β issues. Given that ( )+ +B B B B B22 2 2 2 2,
where B represents the magnetic 5eld discretization error, the
magnetic pressure discretization error can be comparable to
thermal pressure in low-β regions. This can lead to
nonphysical negative thermal pressure when deriving thermal
pressure from energy density. Reducing magnetic 5eld
discretization errors is crucial to mitigating these issues.
By solving decomposed MHD equations, where the
magnetic 5eld B is split into a time-dependent 5eld B1 and a
time-independent potential (e.g., T. Tanaka 1995; K. G. Powell
et al. 1999; F. G. Fuchs et al. 2010; X. C. Guo 2015) or even
an arbitrary background magnetic 5eld B0 (C. Xia et al. 2018),
some previous works (e.g., X. S. Feng et al. 2010, 2021;
C. X. Li et al. 2018; Y. Wang et al. 2019; H. P. Wang et al.
2022a, 2022b, 2025a) have reduced magnetic 5eld discretiza-
tion errors and signi5cantly improved the numerical stability
of quasi-steady-state MHD coronal modeling. However, as B1
increases in the time-evolving coronal simulations and does
not satisfy homogeneous boundary conditions anymore,
discretization errors in B1 are still likely to result in
nonphysical negative thermal pressure and cause the code to
break down. Therefore, L. Griton (2018) incorporated the
rotation of the potential magnetic 5eld B0 in simulations of
planetary magnetospheres, treating B0 as a time-dependent
5eld that corotates with the planet. This approach helps keep
B1 small and satis5es homogeneous boundary conditions
throughout the simulations, thus improving the numerical
stability of the MHD planetary magnetospheric code. Given
that time-evolving coronal simulations involve not only
differential rotation but also magnetic Lux emergence and
cancellation, we further introduce a temporally piecewise-
constant variable in this paper to accommodate part of the
nonpotential 5eld, ensuring that B1 remains consistently small
throughout the simulations of the time-evolving corona.
Based on the above considerations, the paper is organized as
follows. In Section 2, we introduce the numerical formulation
and implementations of the time-evolving MHD coronal
model. The governing equations, the derivation of the
17
https://github.com/andrealani/COOLFluiD/wiki
2
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
extended magnetic 5eld decomposition strategy, the proces-
sing of time-evolving boundary conditions, and the positivity-
preserving (PP) measures applied to enhance the model’s
numerical stability are described in detail. In Section 3, we
present the simulation results, including the evolution of the
corona during CRs 2110 and 2111, and a comparison between
the time-evolving simulation results and observational data.
Finally, in Section 4, we summarize the key features of the
ef5cient and numerically stable fully implicit time-evolving
coronal model and offer concluding remarks.
2. Governing Equations and Numerical Methods
2.1. The Governing Equations
This time-evolving MHD global coronal model is based on
the time-accurate Solar Interplanetary Phenomena Implicit
Finite Volume Method (SIP-IFVM) coronal model
(H. P. Wang et al. 2025a). Additionally, we incorporate the
optically thin radiative loss and adjust the adiabatic index γ
from 1.05 to 5/3 to better represent the adiabatic process.
Furthermore, we develop an extended magnetic 5eld decom-
position strategy, the derivation of which is available in
Section 2.5, and employ it to enhance the numerical stability of
the time-evolving thermodynamic MHD coronal model. The
governing equations are described as follows:
where ( )=B B B B, ,x y z
T and ( )=v u v w, , T denote the
magnetic 5eld and velocity in a Cartesian coordinate system,
respectively, ρ is the plasma density, ( )=B B B B, ,x y z
T
00 00 00 00
is a static potential 5eld, ( )=B B B B, ,x y z
T
01 01 01 01
is a
temporally piecewise-constant 5eld and B0 = B00 + B01,
B1 = B − B0, = + +v BE
p
1
1
1
2
2 1
2
1
2 with the adiabatic
index /= 5 3, and = +p p
B
T1
2
1
2
. During the simulations, as
detailed in Section 2.5, when
B
p
0.5
1
2
falls below a speci5c
threshold, we update B01, B1, and E1 to B01 + B1, 0, and
+ v
p
1
1
2
2, respectively. We then solve Equation (1) with
= 0
B
t
01 to advance the solution in time. This process
continues until
B
p
0.5
1
2
again falls below the threshold, at which
point we repeat the procedure. In the de5nition of the magnetic
5eld, a factor of
µ
1
0
is absorbed, with μ0 = 4 × 10−7πHm−1
denoting the magnetic permeability. Additionally, G is the
universal gravitational constant, Ms denotes the mass of the
Sun, and GMs = 1.327927 × 1020m3 s−2. The thermal
pressure of the plasma is de5ned as p = RρT, where
R = 1.653 × 104m2 s−2K−1 denotes the gas constant, and T
is the temperature of the bulk plasma. In addition, r is the
position vector and | |= rr refers to the heliocentric distance.
SE = Qe + v · Sm + ∇ · q + Qrad is the energy source term,
where Qe and Sm, de5ned the same as by X. S. Feng et al.
(2021) and H. P. Wang et al. (2022a, 2022b), are used to
mimic the effect of coronal heating and solar wind accelera-
tion. Additionally, Qrad and ∇ · q represent radiation loss and
thermal conduction, respectively.
By assuming the radiative losses to be optically thin
(R. Rosner et al. 1978; Y. H. Zhou et al. 2021), the radiative
term Qrad is de5ned as
( ) ( )=Q n n T , 2e prad
where the proton number density np is assumed to be equal to
the electron number density ne for the hydrogen plasma.
Similar to B. van der Holst et al. (2014) and H. P. Wang et al.
(2025a), the radiative cooling function ( )T in this paper is
derived from version 10 of CHIANTI (K. P. Dere et al.
1997, 2023), an atomic database for emission lines. As done
by C. Xia et al. (2011) and H. P. Wang et al. (2025a), we set
( )T to zero when T < 2 × 104K. The heat Lux q is de5ned
in a Spitzer form or collisionless form, mainly according to
the radial distance (J. V. Hollweg 1978):
( )
( )
/
=
>
^ ^
q
b b
v
T T R r R
n k T r R
if 1 10
if 10 ,
3s s
e s
5 2
B
where ˆ =b
B
B
, /= ×9 10 J m s K12 1 1 7 2 (E. Endeve
et al. 2003; X. S. Feng et al. 2010), the parameter α is set to 3/
2 (E. Endeve et al. 2003), ne is the electron number density,
and kB = 1.380649 × 10−23J K−1 is the Boltzmann constant.
The Spitzer form of heat Lux is de5ned along the magnetic
5eld, whereas the collisionless form is de5ned along the
velocity vector. To ensure a smooth transition between these
two forms, we hybridize the Spitzer form qSpitzer and
collisionless form qCollisionless based on the Alfvénic Mach
number over the domain 1 Rs� r� 10 Rs, as described in
Equation (4):
( )
= +
×
q
v
q
v
q
V V
R r R
min 1, 1 min 1,
, if1 10 , 4
a a
s s
Spitzer
Collisionless
( )
( )
( )
[( ) ( )] ( )( )
( ) ( ) ( ) ( )
( ) ( )
( )
( )
+ =
+ + + +
= + +
+ + + = +
+ + +
+ = +
v
v v I B B B B
B B B r S
B B v B v B B B B v B
v B B v B B B B v B v B v r
v B B v B B v
p
E p
S
0
,
1
v B B
B
B B
t
t
G M
r
m
E
t T t
G M
r
E
t t
2 2
00 00
1 01
1
1
1 0 1 1 1 01 1
0 0 0 0
1 01
s
s
2
00
2
3
1 0
3
1 01
3
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
where =V
B
a
0.5
is the Alfvén speed. Meanwhile, we apply a
linear smoothing transition (Z. Mikić et al. 1999; R. Lionello
et al. 2008) between Equation (4) and qCollisionless over the
domain of 7.5 Rs� r� 10 Rs. In our time-evolving coronal
simulation, the plasma velocity occasionally reaches extra-
ordinarily high values, exceeding 1000 km s−1, near low-β
regions. Therefore, we enhance the heat conduction ef5ciency
in low-β regions and ultimately adopt the following plasma β-
dependent heat Lux q
*
:
( )= + ×q q 1 tanh
1 1
100
. 5
2.2. Spatial Discretization
We adopt Godunov’s method to advance cell-averaged
solutions in time by solving a Riemann problem at each cell
interface (S. K. Godunov 1959; B. Einfeldt et al. 1991). The
computational domain is a spherical shell ranging from 1 Rs to
20Rs. The six-component composite mesh (X. S. Feng et al.
2010; X. S. Feng 2020; H. P. Wang et al. 2022a, 2025a) with
each identical component being a low-latitude spherical mesh
con5ned by ( )( )+ × +4
3
4
3
4
5
4
and discretized to 96 × 42 × 42 grid cells, is adopted. Here, δθ
and δf are adjustable parameters proportionally dependent on the
grid spacing entailed for the minimum overlapping area.
As usual, the following discretized integral equations are
obtained by integrating Equation (1) over the hexahedral cell i
and applying Gauss’s theorem to convert the volume integral
of the Lux divergence into a surface integral:
( )+ =
U
RV
t
0, 6i
i
i
where ·=R F n Sd Vi
V
i i
i
denotes the residual operator
over cell i. Here, · ( ) ·=
=
F n F U U nd ,
V
j
ij L R ij ij
1
6
i
; Si =
SPowell,i + Sgra,i + Sheat,i corresponds to the Godunov–Powell
source terms, the gravitational force, and the heating and
radiation loss source terms in cell i, respectively. Ui denotes
the cell-averaged solution variables in cell i. As described in
Section 2.5, when min
B
p
cell 0.5i
i
i1,
2
drops below a threshold, in this
paper we adopt 10, we update B01,i, B1,i, and E1,i to
B01,i + B1,i, 0, and + v
p
i i1
1
2
2i , respectively.
Additionally, in Equation (6), Vi is the volume of cell i, Γij is
the area of interface shared by cell i and its neighboring cell j,
and nij is the unit normal vector of Γij, pointing from cell i to
cell j. Following the approach of X. S. Feng et al. (2021) and
H. P. Wang et al. (2022a), the cell-averaged Powell source
term is calculated as
( )(
( ))
( )
= +
+
=
S T Q U B B B
B B B B B ,
7
i
V
j
ij nL nL nL n i
n i nR nR nL nL ij
Powell,
1
1
6
8
1
1 01 1 ,
01 , 1 01 1 01
i
where =
S
S S
L
R L
(K. L. Wu & C. W. Shu 2019),
=T
T
T
0 0
0 0 0
0 0
0 0 0
1 0
0 1
,ij
ij
ij
8
and ( )=T n t t, ,ij ij ij ij
T
1 2 is a rotation
matrix that transforms the ( )x y z, , coordinate system to the
( )n t t, ,
1 2
coordinate system. Here, t1ij and t2ij are two unit
orthogonal tangential vectors of cell face Γij (e.g.,
X. S. Feng 2020, and references therein), SL and SR are the
velocities of two fast waves as de5ned by H. P. Wang et al.
(2022a), and ( ) ( )=Q U T B v B T v0, , ,nL ij i L i L i L ij i L T, , 1 , , .
The subscripts “i,” “L,” and “R” denote the corresponding
variables at the centroid of cell i, and on Γij extrapolated from cell
i and from cell j, respectively. Moreover, Sgra,i and Qrad,i included
in Sheat,i are de5ned by evaluating their respective formulations
using the corresponding variables at the centroid of cell i. The
inviscid Lux through the interface Γij, described as
( ) ·F U U n,ij L R ij, is computed using the positive-preserving (PP)
Harten–Lax–van Leer (HLL) Riemann solver (X. S. Feng et al.
2021), equipped with a self-adjustable dissipation term
(H. P. Wang et al. 2025b). The cell-averaged heat conduction
term ( · )q i is calculated following Gauss’s theorem, as
described by H. P. Wang et al. (2022a). The volume heating
coef5cients in Qe,i and Sm,i are precomputed from six-hourly
updated magnetograms and are stored as a time series and then
linearly interpolated to the current time step during the time-
evolving coronal simulations.
Following H. P. Wang et al. (2025a), the second-order PP
reconstruction method is used to reconstruct the piecewise
polynomials of primitive variables {ρ, u, v, w, p, T} on the cell
surface Γij,
( ) · ( )
{ } ( )
= +x x xX X X
X u v w p T
,
, , , , , , 8
i i i i i
where X i is the corresponding variable at xi (the centroid of
cell i), ( )=X , ,i Xx
X
y
X
z
i
is the derivative of X at xi, and
Ψi is the limiter used to control spatial oscillation. Further-
more, a globally solenoidality-preserving approach (X. S. Feng
et al. 2019, 2021) is employed to enforce the divergence-free
constraint for the magnetic 5eld. Considering that reducing the
magnetic 5eld discretization error is crucial for enhancing the
PP property of MHD models in low-β regions (H. P. Wang
et al. 2025a), and given that the MHD decomposition method
introduced in this paper has improved the numerical stability
of our code, we opt to abandon the limiter as follows in the
calculation of the piecewise polynomial representations for
B00(x), B0(x), and B1(x) to minimize discretization error in the
magnetic 5eld:
( ) · ( )
{ }
( )
= +x x x
B B B B B B B B B
X X X
X
,
, , , , , , , , .
9
i i i i
x y z x y z x y z00 00 00 0 0 0 1 1 1
2.3. Temporal Integration
In this paper, we 5rst apply the following backward Euler
temporal integration to Equation (6) to perform a quasi-steady-
4
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
state coronal simulation constrained by one time-invariant
magnetogram (X. S. Feng et al. 2021; H. P. Wang et al.
2022a):
( )+ =+
U
RV
t
0. 10i
i
n
i
n 1
The superscripts “n” and “n+1” denote the time level,
=
+
U U Ui
n
i
n
i
n1 is the solution increment between the
nth and (n + 1)th time levels, and Δt = t n+1 − t n is the time
increment. Subsequently, we evolve the magnetograms to
carry out the time-evolving coronal simulation, enabling us to
capture the temporal evolution of coronal structures. To
improve the temporal accuracy for time-evolving simulations,
we further employ the pseudo–time marching method
(H. P. Wang et al. 2025a), which introduced a pseudotime τ to
Equation (10) and updated the solution during each physical
time step Δt by solving a steady-state problem on τ, to solve
Equation (10):
( )+ + =+
U U
RV V
t
0. 11i
i
i
i
n
i
n 1
Here, Δτ is a pseudo–time step, and ΔUi is the solution
increment during Δτ.
Similar to the PP measure employed by H. P. Wang et al.
(2025b), we make the following adjustment on the updated
density and thermal pressure during each pseudotime iteration
in the domain of 1 Rs� r� 2 Rs:
( )
( )
( )
= +
= +
B
B
V
p p
R r R
1
2
1
, if1 2 . 12
i
i
A
o i
i p
i
p o i
s s
2
,max
2 ,
2
min ,
i i
i i
In Equation (12), ( )= + ×0.5 0.5 tanh V VVi A i A, ,maxfac with
| |
=V
B
A i,
i
o i,
0.5
, =V 2
B
A,max
s
max
0.5
, and =V
V
fac
1000
A,max . Bmax repre-
sents the maximum magnetic 5eld strength in the entire
computational domain. Moreover, = +0.5 0.5 tanhpi
×
+B
pi
i
min
0.5 2
fac
with βfac = 10
−7, = 10min
4, and
ε = 10−12 are adopted in this paper. Additionally, we constrain
the plasma velocity in the range 1 Rs� r� 1.1 Rs not to exceed
the speed of sound ( )=Cs i p,
0.5
i
i
, as described in
Equation (13):
( )( )
( )
| |
= × × + ×
13
v v
R r R
min 0.3 tanh 8.68 , 1.0 ,
if 1 1.1 .
v
i o i
C r R
R
s s
,
s i
o i
s
s
,
,
Here, ρ, p, and v with the subscript “o” denote the density,
thermal pressure, and velocity updated during the pseudotime
iterations without adjustment.
During the time-evolving coronal simulation, the time step
is gradually increased from an explicit time step to τLow, a
reference time that is the same as de5ned by X. S. Feng et al.
(2021) and H. P. Wang et al. (2022a). Additionally, we noticed
that during the long-term time-evolving simulations, a
dramatic increase in the magnetograms, where the maximum
magnetic 5eld strength increases by more than 1.5 times, may
occur between two adjacent magnetograms, potentially caus-
ing the code to crash. In this paper, this phenomenon occurred
in the magnetograms at the 882nd and 888th hours of the
approximately 1300 hr time-evolving simulation. Assuming
such cases occur at a moment between tm−1 and tm, where the
subscripts “m−1” and “m” correspond to the ( )m 1 th and mth
magnetograms, we opt to employ a second-order Runge–Kutta
method, with the intermediate states U(1) and U(2) computed
using the backward Euler method (X. S. Feng et al. 2021;
H. P. Wang et al. 2022a, 2022b), to solve Equation (6) over the
time interval from tm−2 to tm+3:
( ) ( )
( ) ( )
( ) ( ) ( )
( )
=
=
= +
+
U U R
U U R
U U U
t
t
,
,
. 14
i i
n
i
i i i
i
n
i
n
i
1 1
2 1 2
1
1
2
2
During this time interval, the inner-boundary magnetic 5elds
between tm−2 and tm+1 are linearly interpolated from the
magnetograms at tm−2 and tm+1.
2.4. Implementation of Inner Boundary Conditions
In this paper, we 5rst perform a quasi-steady-state coronal
simulation constrained by a time-invariant magnetogram
(H. P. Wang et al. 2022a), and then evolve the magnetograms
to drive the subsequent time-evolving coronal simulation
during CRs 2110 and 2111. Around this period, the time
interval between two adjacent GONG-ADAPT magnetograms
is 6 hr. The original GONG-ADAPT magnetograms adopted in
this paper are positioned in the corotating Carrington
heliographic coordinate system. Therefore, we 5rst rotated
these magnetograms to the heliocentric inertial coordinate
system to match the inertial coordinate system. Subsequently,
cubic Hermite interpolation is applied to these input
magnetograms to derive the required inner-boundary magnetic
5elds for each time step during the time-evolving coronal
simulations (H. P. Wang et al. 2025b).
As done by H. P. Wang et al. (2025a), the inner boundary
conditions are speci5ed at the inner boundary face, which
coincides with the solar surface. Four Gaussian points are used
on each inner boundary face. The inner boundary conditions at
these Gaussian points, together with the solutions in the
boundary cell and its 17 neighboring cells that share at least a
vertex with it, are utilized to construct the reconstruction
formulation of primitive variables in the inner boundary cells
(X. S. Feng et al. 2021; H. P. Wang et al. 2025a). For
convenience of description, we take the ith boundary cell,
denoted as cell BD,i, which is a hexahedral cell consisting of
one curved boundary surface and 5ve planar faces (X. S. Feng
et al. 2021; H. P. Wang et al. 2022a), as an example. The
reconstruction formulations of primitive variables within cell
BD,i are presented as follows:
( ) ( ) · ( )
{ } ( )
= +x x x
B B B
X X X
X u v w p T
,
, , , , , , , , . 15
i i i i iBD, BD, BD, BD, BD,
1 00 0
The subscript “BD,i” refers to the variable corresponding to
the ith boundary cell BD,i.
As the inner boundary magnetic 5eld evolves, a tangential
component of the electric 5eld, ( )=E E E,tBD, BD, BD, , arises
at the inner boundary. Consequently, the corresponding
5
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
tangential velocity component, ( )=v v v,tBD, BD, BD, is given
by
( )= =v
E
B
v
E
B
, , 16
r r
BD,
BD,
BD,
BD,
BD,
BD,
ensuring that the Lux evolution match the changes of the
observed radial magnetic 5eld (L. P. Yang et al. 2012;
X. S. Feng et al. 2023). Here, the subscripts “θ” and “f” denote
the respective velocity and magnetic 5eld components along
the θ and f directions. Additionally, following R. Lionello
et al. (2023), we apply the following adjustment to
Equation (16) to regularize the Low near the boundary polarity
inversion line, where BBD,r = 0:
( )
| | | |
( )
/
=
×
+
v
E B
B E BD C
. 17t
t
s
BC,
BC BC
BC
2
BC BC
This means that the magnitude of the tangential velocity will
converge to /C Ds , where D= 3 is adopted in this paper, and Cs
is the local sound speed calculated from the solution in the cell
immediately adjacent to the inner boundary cell in the radial
direction.
During the simulations, the boundary conditions for thermal
pressure pBD, plasma density ρBD, temperature TBD, and radial
velocity vBD,r are classi5ed into two cases based on the radial
velocity vr in the cell immediately adjacent to the inner
boundary cell in the radial direction (C. P. T. Groth et al. 2000;
X. S. Feng et al. 2021; H. P. Wang et al. 2022a, 2022b, 2025a).
(i) vr� 0: The thermal pressure and plasma density at the
inner boundary face are set as /=p 1
BD
and ρBD = 1 in code
units. Following M. Brchnelova et al. (2023) and H. P. Wang
et al. (2025b), the inner boundary density ρBD is further
adjusted based on the local Alfvén speed, as in Equation (12),
to enhance the PP property of the coronal model. Additionally,
the inner boundary temperature is calculated as =T
p
BD
BD
BD
,
and the radial velocity component is constrained as = 0
v
r
rBD, .
(ii) vr� 0: = 0
r
BD , = 0
p
r
BD , = 0
T
r
BD , and vBD,r = 0 are
applied.
Additionally, the variables at the centroid of inner boundary
cells are required to calculate Equation (15). This paper
generates the magnetic 5eld from the potential 5eld solver
with a 15-order spherical harmonic expansion. The tangential
velocity is computed as an average of the values at four
Gaussian points and at the centroid of the neighboring cell
immediately adjacent to the inner boundary cell in the radial
direction. The radial velocity, thermal pressure, plasma
density, and temperature are derived using Parker’s one-
dimensional hydrodynamic isothermal solar wind solution
(E. N. Parker 1963). Additionally, the thermal pressure and
plasma density are adjusted using Equation (12).
2.5. Derivation of the Extended Decomposed MHD Equations
Considering that ( )+ +B B B B B22 2 2 2 2 with
εB denoting the discretization error in the magnetic 5eld B, the
magnetic pressure discretization error can be comparable to the
thermal pressure in low-β regions, and nonphysical negative
thermal pressures are prone to appear when deriving the
thermal pressure from energy density. To avoid such
undesirable situations, traditional decomposed MHD equations
are commonly used, where the magnetic 5eld B is divided into
a time-independent potential 5eld and a time-dependent 5eld
B1, with B1 being used in the derivation of the thermal
pressure. However, discretization errors in B1 are still likely to
result in nonphysical negative thermal pressure as | |B1
increases in the time-evolving simulations, potentially causing
the code to break down. Therefore, we propose the extended
magnetic 5eld decomposition strategy, in which the magnetic
5eld B is divided into a time-independent potential 5eld B00, a
temporally piecewise-constant 5eld B01, and a time-dependent
5eld B1. In the following, we describe how the extended
decomposed MHD equations are derived.
With B00 representing a static potential 5eld, and B01 being
a temporally piecewise-constant 5eld that remains unchanged
during the time interval between t n and t n+ k, where the
subscripts “n” and “n+k” denote the nth and ( )+n k th time
levels in a simulation, the following conditions are satis5ed
within the time interval between t n and t n+ k:
[ ]
( )
( )
= = × =
=
+
+
B B t t t
t t t
0
0
, 0, 0, ,
, , .
18
B
B
t
n n k
t
n n k
00 00
00
01
The original energy equation is described as
· [( ) ( · )] ( · ) · ( )+ + =v B v B B v B
E
t
E p , 19T
where = + +v BE
p
1
1
2
2 1
2
2 and = +p p
B
T 2
2
. Given that
B = B00 + B01 + B1 and = + +v BE
p
1
1
1
2
2 1
2
1
2, the
formulation of E can be described as
· ( )= + +B B BE E
1
2
, 20
1 0
2
0 1
where B0 = B00 + B01.
The original induction equation is described as
· ( ) ( · ) ( )+ =
B
v B B v B v
t
. 21
From Equations (18) and (20) we get
· · · ( )= + + +B
B
B
B
B
BE
t
E
t t t t
, 22
1
0
0
0
1
1
0
( )
( )=
+B B B
t t
. 23
1 01
From Equations (23) and (21), we get
· ( ) · ( )=
B
v B B v v B
B
t t
. 24
1 01
Multiplying both sides of Equation (24) by B0 results in
· · ( ) ·
· · · ( )
=B
B
v B B v B
v B B B
B
t
t
. 25
0
1
0
0 0
01
6
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
Considering that
[ · ( )] · ( · )( · ) ( · ) ·
· [ ( · )] · ( · ) ( · ) ·
= +
= +
v B B v B B v B B
v B B v B B v B B ,
0 0 0
0 0 0
[ · ( )] · ( · )( · ) ( · ) ·
· [ ( · )] · ( · ) ( · ) ·
= +
= +
B v B B v B B v B
B v B B B v B v B ,
0 0 0
0 0 0
Equation (25) is equivalent to
( )
· · [( · ) ( · ) ]
· · ·
=
+ 26
B
B
B B v v B B
v B B B
B
t
t
,
0
1
0 0
0 0
01
where · ( · )= v B B0 –( · )v B · · ( · )B B B v0 0
+ ( · ) ·B v B0.
From Equations (26) and (22) we obtain
· [( · ) ( · ) ]
· · · ·
· ( )
=
+
+ +
B B v v B B
v B B B
B
B
B
B
B
E
t
E
t
t t
t
. 27
1
0 0
0 0
01
0
0
1
0
From Equations (19) and (20) we obtain
· ·
( · ) ] · · ( )
= + + + +
×
B B B B
v v B B v B B
E
t
E p 2
1
2
. 28
1 0
2
0 1 1
2
From Equations (28) and (27) we obtain
· · ( · )
· · ·
· ·
( )
+ + + +
= +
+
B B B v v B B
v B B B
B
B
B
B
B
E
t
E p
t
t t
1
2
.
29
1
1 0 1 1
2
1
1 0
01
0
0
1
0
Since ∇ · B00 = 0 and =
B B
t t
0 01 , we obtain the following
decomposed energy equation:
· [( · ) ( · )]
· · ( ) · ( )
+ + +
= +
B B v B v B
v B B B B
B
E
t
E p
t
, 30
T
1
1
1
1 0 1
1 1 01 1
0
where = +p p
B
T1 2
1
2
. This indicates that the thermal pressure,
p, can be determined from E1. As a result, the accuracy of p is
no longer limited by the discretization error of B but rather by
that of B1. As long as B1 remains small, the discretization error
in B
1
2
1
2 is unlikely to approach the magnitude of the thermal
pressure, reducing the risk of nonphysical negative thermal
pressure.
The original momentum equation is described as
( )
· ( · )
( )
+ + + =
v
v v
B
I B B B B
t
p
2
.
31
2
Given that ∇ × B00 × B00 = ∇ ( )· +B I B B12 002 00 00
( · )B00 B00 and ∇ · B00 = 0, Equation (31) is equivalent
to
( )
( )
·
· ( )
+ + + +
= +
32
v
v v
B B
I B B B B
B B B
t
p
2 2
.
2
00
2
00 00
1 01
Consequently, we obtain the following extended decom-
posed MHD equations:
( )
( )
( )
[( ) ( )] ( )( )
( ) ( ) ( ) ( )
( ) ( )
( )
+ =
+ + + + = +
+ + + = +
+ +
+ = +
v
v
v v
B B
I B B B B B B B
B B v B v B B B v B
B
B
v B B v B B B B v B v B
B
v B B v B B v
B
t
t
p
E
t
E p
t
t t
0
2 2
.
33
T
2
00
2
00 00 1 01
1
1
1
1 0 1 1 01 1
1
0
0 0 0 0
1
1 01
01
7
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
Moreover, we consider a speci5c condition in which B01 is
arti5cially increased by B1 over an in5nitesimally short period
denoted by
+
dtlim
dt 0
. This adjustment is applied when B1
becomes signi5cantly large at time t n, such as when
B
p
0.5
1
2
drops below a certain threshold. Correspondingly, B1, B0, E1,
and pT1 evolve to B1 , B0 , *E1 , and *p
T1
, respectively, ensuring
that the MHD system remains unchanged during the time
interval
+
dtlim
dt 0
. As a result, Equations (24) and (30) can be
described as
· ( )
· ( )
=
+
+
*
*
B B
v B B v
v B
B B
dt
dt
lim
lim , 34
dt
dt
0
1 1
0
01 01
and
( )
· [( · ) ( · )]
· · ( ) ·
+ + +
= +
+
+
+
+
*
* * * * *
* * * *
*
*
*
35
B B v B v B
v B B B B
B B
B B
E E
dt
E p
dt
dt
lim
lim
1
2
lim ,
dt
T
dt
dt
0
1 1
1 1 1 0 1
1 1 01 1
0
0
1
2
1
2
where = +B B B01 01 1, = +B B B0 1 , | |= + Bp p
T1
1
2
1
2,
and [ ( )]= v B B0 ( )v B B B0 [ ( )]B v0
( )+ B v B0 . Consequently, B 01 can be derived from
Equation (34). Keeping in mind that = +B B B0 00 01, B
*
→ B,
and + =B B B0 1 , we can derive | |BE E1 1
1
2
1
2 from
Equation (35). At the end of this in5nitely short period, B01,
B1, and E1 become B01 + B1, 0, and + v
p
1
1
2
2,
respectively.
Therefore, when
B
p
0.5
1
2
drops below a threshold at time t n,
we update B01 to B01 + B1, set B1 = 0, and assign
= + vE
p
1
1
1
2
2 and pT1 = p. Subsequently, we solve
Equation (33) with = 0
B
t
01 to advance the solutions in time.
This process continues until
B
p
0.5
1
2
once again drops below the
threshold, at which moment we repeat the procedure outlined
above.
3. Numerical Results
In this section, the time-evolving SIP-IFVM coronal model,
equipped with the extended magnetic 5eld decomposition
strategy, is employed to simulate the evolution of coronal
structures during CRs 2110 and 2111. This period, spanning
from 2011 May 9 to July 3, corresponds to the solar maximum
of Solar Cycle 24. Approximately 220 GONG-ADAPT
magnetograms,18 with the 5rst realization of the 12-member
ensemble (B. Perri et al. 2023) adopted and updated at a 6 hr
cadence, are used to drive these simulations in an inertial
coordinate system spanning from the solar surface to 20 Rs.
During the time-evolving simulations, the magnetic 5eld near
the solar surface sometimes exceeds 40 G, with the plasma β
reaching a minimum value of around 10−3. Figure 1 displays
the inner boundary magnetic 5eld at various moments. For
convenience of comparison, these magnetic 5eld distributions
are illustrated in a corotating coordinate system. The 5gure
reveals that the inner-boundary magnetic 5eld varies more
signi5cantly at low and middle latitudes, while the regions
beyond 70° in both the north and south poles are predomi-
nantly occupied by magnetic 5elds directed inward and
outward from the Sun. Additionally, a dipole with a strong
magnetic 5eld is observed in the panel at 1042 hr. Video 1
(linked in Figure 1) illustrates the evolution of the ADAPT-
GONG maps during this period, suggesting that this dipole
may have originated on the Sun’s farside and was later updated
with nearside data. This leads to the sudden appearance of a
Figure 1. Distribution of the radial magnetic 5eld used as the inner boundary condition at the solar surface, shown in a corotating coordinate system. Its evolution
throughout the simulation (from t = 3 to 1318) is demonstrated in an animation.
(An animation of this 5gure is available.)
18
https://gong.nso.edu/adapt/maps/gong/2011/
8
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
dipole with a strong magnetic 5eld on the magnetogram at the
888th hour.
All calculations in this paper are performed on the WICE
cluster, part of the Tier-2 supercomputer infrastructure of the
Vlaams Supercomputer Centrum (VSC).19 The simulation
consists of 23,040 time steps, with the average time step being
approximately 3.4 minutes. Using 192 CPU cores, the wall-
clock time for simulating these two solar maximum CR
periods is approximately 16 hr. In addition to the simulation
using the extended magnetic 5eld decomposition strategy, we
also conduct a simulation with the traditional method for
comparison. However, the latter crashes after approximately
20 hr of physical time. When the simulation using the
traditional method fails at a time step between tm−1 and tm,
we restart from tm−2 and apply Equation (14) to solve the
governing equations, as described in Section 2.3. Simulta-
neously, we reduce the time-step limit to half of its previous
value. If the simulation crashes again within the same interval,
we revert to tm−2 and further reduce the time-step by half. If
the failure persists for four consecutive attempts, the reversion
process is discontinued.
Figure 2 illustrates the evolution of physical time versus
computational time for both simulations, using the extended
magnetic 5eld decomposition strategy (black solid lines) and the
traditional method (gray solid lines). The left panel demon-
strates that the evolution of coronal structures over 1322 hr of
physical time is completed within just 16 hr of computational
time, showcasing that our model with the extended magnetic
5eld decomposition strategy operates over 80 times faster than
real-time coronal evolution using only 192 CPU cores. The right
panel, a close-up view of the pro5le between 876 and 912 hr of
physical time, is around the event of a dramatic increase in the
magnetograms described in Section 2.3 and is completed within
0.18 computational hours, demonstrating an ef5ciency approxi-
mately 1.5 times greater than during other periods. Although the
strategy outlined in Equation (14) may overlook some high-
frequency phenomena during this interval and reduce temporal
accuracy, it effectively prevents simulation crashes during
periods of dramatic increases in inner-boundary magnetic 5eld
strength. Furthermore, the middle panel shows that the
simulation using the traditional method crashed after around
20 hr of physical time, even with a small time step. This
highlights the effectiveness of the extended magnetic 5eld
decomposition strategy in handling time-evolving low-β issues.
3.1. Distributions of the Open and Closed Magnetic Field
Regions
Coronal holes (CHs) are dark regions observed in extreme
ultraviolet (EUV) and soft X-ray channels, typically associated
with low plasma density and magnetic 5eld lines that are open
to interplanetary space. CH distributions vary across different
phases of solar activity and are among the most prominent
features of the solar corona (G. J. D. Petrie et al. 2011;
X. S. Feng et al. 2015, 2017, 2019). Generally, three types of
CHs can be identi5ed in EUV and soft X-ray images of the
solar corona. Polar CHs are located at the solar poles and often
extend to lower latitudes, occasionally crossing the solar
equator. Isolated CHs, commonly observed near solar maxima,
are detached from polar CHs and scattered across low and
mid-latitudes. Transient CHs are associated with solar eruptive
events, such as CMEs, solar Lares, and eruptive prominences.
In Figure 3, we derive the simulated open- and closed-5eld
regions by tracing the magnetic 5eld lines from 2.5 Rs to the
solar surface (middle and bottom panels) and compare them
with the synoptic maps of the observations (top panel) from
the 193 Å channel of the Atmospheric Imaging Assembly
(AIA) instrument (J. R. Lemen et al. 2012) on board the Solar
Dynamics Observatory (SDO) satellite (W. D. Pesnell et al.
2012).20 These images are illustrated in a corotating coordinate
system with these synoptic images generated by concatenating
a series of meridional strips extracted from full-disk images
over the duration of a complete CR (A. Hamada et al. 2018).
They reveal that the simulations adequately reproduce the
observed southern polar CH from latitude 70° poleward. The
leading CHs around 220° and 280° in longitude, extending
from the southern pole to near the solar equator, as well as the
isolated CH centered at (θlat, flong) = (40°N, 60°) and
detached from the northern polar CH, are also well captured.
Here, “θlat” represents the heliographic latitude, and “flong”
denotes the Carrington longitude. Additionally, the decreasing
trend of scattered isolated CHs in the low-latitude domain
from CR 2110 to 2111 is also captured in the simulation
results. Notably, the isolated CH centered around
(θlat, flong) = (0°, 150°) rapidly transformed from an upside-
down hook shape at 782 hr to a whirlwind-like structure above
the solar equator at 1042 hr. Together with Figure 1, this rapid
change can be attributed primarily to the emergence of a dipole
centered at (θlat, flong) = (10°N, 160°).
Figure 2. Evolution of physical time vs. computational time in the time-evolving coronal simulation for CRs 2110 and 2111 (left), with close-up views of the pro5le
during the 5rst 1.2 hr (middle) and between 10.24 and 10.65 hr of computational time (right).
19
https://www.vscentrum.be/ 20 https://sdo.gsfc.nasa.gov/data/synoptic/
9
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
In Figure 4, we present the synthesized 171 Å EUV images
in a corotating coordinate system. The synthesized EUV
images are derived by integrating the emission along the radial
direction from 1.02 Rs to 1.5 Rs. The emission at each point x
of our computational domain is calculated via
( ) ( ) ( )=x xI G T n n, ,e e
2
where the wavelength λ = 171Å and Gλ is the temperature-
dependent response function for the 171Å wave band, which is
obtained from the CHIANTI atomic database (K. P. Dere et al.
1997, 2023). The 171 Å EUV images primarily show the solar
corona at temperatures of approximately 0.8–1.0 MK. Since the
lifetime of coronal loops typically ranges from minutes to hours
(F. Reale 2014), and cannot be resolved by the six-hourly updated
synoptic magnetograms, the simulation results cannot accurately
reproduce the bright structures observed in the AIA 171 Å EUV
images. Therefore, we do not compare our results with the AIA
171 Å EUV observations. Alternatively, we compare them with
the distributions of the open magnetic 5eld regions in Figure 3.
It can be seen that dark regions occupy the open-5eld polar
region in Figure 3 in the 171 Å EUV images. At 262 hr, the
dark regions centered at (θlat, flong) = (15°N, 50°),
(θlat, flong) = (15°S, 150°), (20°S, 215°) and (θlat, flong) =
(15°S, 280°) align with open-5eld regions, whereas the dark
region centered at (θlat, flong) = (40°N, 40°) corresponds to a
closed-5eld region. At 522 hr, the dark regions centered at
(θlat, flong) = (20°S, 45°), (θlat, flong) = (20°S, 150°), (25°
N, 250°) and (θlat, flong) = (10°N, 280°) align with open-5eld
regions, whereas the dark region centered at (θlat, flong) =
(45°S, 170°) corresponds to a closed-5eld region. At 782 hr,
the dark regions centered at (θlat, flong) = (20°S, 210°) and
(θlat, flong) = (20°S, 285°) align with open-5eld regions, while
the dark region centered at (θlat, flong) = (25°N, 100°)
corresponds to a closed-5eld region. At 1042 hr, the dark
regions centered at (θlat, flong) = (20°S, 215°),
(θlat, flong) = (10°S, 255°) and (θlat, flong) = (15°N, 320°)
align with open-5eld regions, whereas the dark region centered
at (θlat, flong) = (20°S, 60°) corresponds to a closed-5eld
region. Additionally, the extremely bright region centered at
(θlat, flong) = (15°S, 170°) may be attributed to the transition
of the open-5eld region in this area around 782 hr to a closed-
5eld con5guration at 1042 hr. This demonstrates that most of
the dark regions in these synthesized 171Å EUV images
correspond to open-5eld regions, while the bright regions are
associated with closed-5eld regions. The discrepancies may be
Figure 3. Synoptic maps of the EUV observations from the 193 Å channel of AIA on board SDO (top) for CRs 2110 (left) and 2111 (right), alongside the
distributions of open- and closed-5eld regions modeled by the time-evolving SIP-IFVM coronal model (middle and bottom), all shown in a corotating coordinate
system. In the middle and bottom panels, the white and black patches represent open-5eld regions with magnetic 5eld lines pointing outward and inward relative to
the Sun, respectively, while the gray patches indicate closed-5eld regions.
10
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
attributed to the relatively low resolution of the computational
mesh, the lack of a self-consistent heating mechanism, and the
fact that dark regions in EUV maps do not necessarily
correspond to open-5eld regions (M. A. Reiss et al. 2024).
3.2. Time-evolving Coronal Structures
White-light polarized brightness (pB) images can reveal
various large-scale coronal structures. In these pB images,
bright regions correspond to high-density coronal structures,
such as bipolar and pseudostreamers. In contrast, dark regions
indicate low-density structures, such as coronal holes
(X. S. Feng et al. 2015, 2017, 2019; X. S. Feng 2020). Bipolar
streamers separate CHs with opposite magnetic polarities,
while pseudostreamers separate CHs of the same polarity.
Additionally, bipolar streamers extend outward several solar
radii from the Sun, forming a cusplike structure as they are
drawn into a current sheet above the helmet streamer
(Y. M. N. R. Wang et al. 2007; X. S. Feng et al. 2017, 2019).
In Figure 5, we compare white-light 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) spacecraft21 (R. A. Howard et al. 2008)
with those synthesized from the simulation results ranging
from 2.5 Rs to 15 Rs (middle). The 2D distributions of some
selected simulated magnetic 5eld lines on the meridional plane
in the STEREO-A view, ranging from 1 Rs to 6 Rs, are also
illustrated (bottom). This comparison indicates that the
simulated magnetic 5eld lines and bright structures are
generally consistent with the observed bipolar and pseudos-
treamers. However, some discrepancies exist in the width and
position of these bright structures. The mismatch may come
from the imperfect measurements of photospheric magnetic
5elds and the possible presence of solar disturbances during
the two CRs that are not considered in the MHD modeling.
At 442 hr into the time-evolving coronal simulation,
corresponding to 2011 May 28, both the simulation and
observation exhibit two bright structures centered around 44°S
and 46°N at the east and west limbs, respectively. The magnetic
5eld lines indicate that bipolar streamers form these bright
structures. However, in the simulation, the observed bright
structures centered around 10°N and 20°S at the east and west
limbs are shifted northward by approximately 40° and 20°,
respectively. At 602 hr, corresponding to 2011 June 3, the three
bright structures centered around 22°S at the east limb and 57°N
and 7°N at the west limb are well reproduced in the simulation,
formed by bipolar streamers. However, the pseudostreamers
located around 34°S and 61°N at the west and east limbs appear
10° and 20° farther north than in the observations. At 1102 hr,
corresponding to 2011 June 24, the simulated bipolar streamer
around 12°S at the west limb successfully captured the observed
bright structure, while the simulated bipolar streamer around 22°
S at the east limb is shifted 20° northward compared to
observations. Additionally, the pseudostreamers around 48°N at
both limbs align well with the observed bright structures. At 1262
hr, corresponding to 2011 July 1, the simulated bipolar streamer
around 20°S at the east limb and 6°N at the west limb effectively
reproduced the observed structures. The pseudostreamers around
45°N and 59°N at the west and east limbs also closely match the
observations. However, the observed bright structure around 18°S
is not present in the simulation results despite the existence of a
pseudostreamer in the magnetic 5eld structure.
Additionally, we compare the pB observations at 3 Rs with the
modeled results in Figure 6. The top panels illustrate synoptic
maps of the east-limb observations from the Large Angle and
Spectrometric Coronagraph C2 (LASCO-C2) (G. E. Brueckner
et al. 1995) on board the Solar and Heliospheric Observatory
(SOHO)22 for CRs 2110 (left) and 2111 (right), with the bright
structures representing distributions of the high-density
coronal structures. The middle and bottom panels display the
2D timing diagrams of simulated plasma number density
Figure 4. Synoptic maps of the EUV images derived from the time-evolving simulation results at four different moments.
21
https://stereo-ssc.nascom.nasa.gov/browse/ 22 https://sdo.gsfc.nasa.gov/data/synoptic/
11
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
(105 cm−3) and radial velocity Vr (km s
−1
). These 2D timing
diagrams are synthesized from a series of time-evolving
simulation results with a cadence of one result per 20 hr. The
radial basic function interpolation method (H. P. Wang et al.
2022a) is applied to interpolate the variables to the east-limb
longitude of the Sun in the Earth view. The orange solid lines
overlaid on these counter denote the magnetic neutral lines
(MNLs) modeled by the time-evolving coronal model.
The comparison in Figure 6 demonstrates that the bright
structures in the simulated plasma density are generally
consistent with those in the pB observations. The simulated
high-density, low-speed Lows are primarily distributed around
the MNLs. Additionally, it can be noted that the northernmost
bright structures in the simulation for CR 2110 extend
approximately 20° farther north during the 5rst 320 hr.
Meanwhile, for CR 2111, the bright structures during the 5rst
270 hr of the simulation changed from a trapezoidal cavity to
irregular solid triangles in the observations. These mismatches
may be attributed to the fact that the magnetic 5eld at different
longitudes in each synoptic magnetograph is observed at
different times, as well as the imperfect measurements of the
photospheric magnetic 5elds in the polar regions.
During the time-evolving coronal simulation, two virtual
satellites are placed at 3 Rs and 20 Rs, maintaining the same
latitude as Earth and lagging 60° behind in longitude. In
Figure 7, we present the timing diagrams of radial velocity Vr
(km s−1, left), proton number density (105 cm−3 at 3 Rs and
103 cm−3 at 20 Rs, middle), and decadic logarithms of plasma
β (top right) and radial magnetic 5eld polarities (bottom right),
as monitored by the virtual satellites positioned at 3 Rs (top)
and 20 Rs (bottom).
The top panels of Figure 7 indicate that the plasma β at 3 Rs
around the solar equator usually varies between 0.3 and 100
during the time-evolving simulation, which is consistent with
the values derived by G. A. Gary (2001). It is also observed
that a peak/trough in the plasma density pro5le generally
corresponds to a trough/peak in the velocity pro5le. However,
the peaks in plasma density around 600 and 675 hr are
narrower than the corresponding troughs in the velocity
pro5le, while the trough in plasma density around 850 hr is
Figure 5. White-light pB images observed by COR2/STEREO-A (top), ranging from 2.5 Rs to 15 Rs, alongside corresponding pB images synthesized from
simulation results ranging from 2.5 Rs to 15 Rs (middle) and 2D distributions of some selected simulated magnetic 5eld lines from 1 Rs to 6 Rs (bottom). These
images are shown on the meridional plane in the STEREO-A view.
12
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
broader than the corresponding velocity peak. These mis-
matches reveal the complexity of dynamic mechanisms in the
subsonic/sub-Alfvénic coronal region, highlighting the
requirement for careful consideration of the transformation
between kinetic and magnetic energy in simulations of such
regions.
At 20 Rs, the velocity peaks and troughs align more closely
with plasma density troughs and peaks, indicating that the
dynamics in this supersonic/super-Alfvénic region are sig-
ni5cantly simpler than those in the subsonic/sub-Alfvénic
coronal region. Additionally, we noticed that the two velocity
troughs appearing around 740.5 and 166.5 hr at 3 Rs shifted to
748.2 and 175.0 hr at 20 Rs, respectively. This indicates that
the two troughs propagate from 3 Rs to 20 Rs with an average
velocity of approximately 350 km s−1. By performing a
ballistic propagation, we map the solar wind velocity and
radial magnetic 5eld polarities observed by the Wind
satellite23 (J. H. King & N. E. Papitashvili 2005) at 1 au to 20
Rs (top, gray solid lines). Considering that some velocity peaks
observed at 1 au appear roughly 100 hr earlier when mapped to
20 Rs, we positioned the virtual satellite 60° behind Earth to
compare with in situ observations at 1 au near Earth.
The bottom panels of Figure 7 show that the simulation
captures the observed velocity peaks around 73 hr and 750 hr.
The simulated velocity peak, approximately 580 km s−1,
occurs at 280 hr in the simulation, while it appears at 410 hr
in the observations. Additionally, the velocity exceeding
580 km s−1, which persists between 370 and 410 hr in the
observations, is missing in the simulation. The observed
velocity peak around 1000 hr appears about 70 hr earlier than
in the simulation. Additionally, the simulation approximately
captures the observed velocity troughs around 180 hr and 880
hr. The two observed velocity troughs, centered at 500 and 650
hr, combine into one wide trough in the simulation, while the
observed trough centered at 1100 hr occurs about 150 hr
earlier than in the simulation. As for the radial magnetic
polarities, the simulation captures 83.4% of the observations.
These discrepancies in velocity and radial magnetic polarities
may be attributed to the limitations of the empirical heating
function used to approximate coronal heating and solar wind
acceleration, as well as the constraints of synoptic magneto-
grams, where the magnetic 5eld at different longitudes is
observed at different times. Since this period occurs during
solar maximum, the frequent eruptions, such as coronal mass
ejections, which are not included in this model, may also
contribute to the discrepancies between the simulation results
and observations.
In Figure 8, we further present the 2D timing diagrams of
simulated plasma number density (103 cm−3, top), radial
velocity Vr (km s
−1, middle), and temperature (105K, bottom)
to show whether the code produces the latitudinal structure of
the solar wind, which is basically the latitude distribution of
the fast and slow solar wind. These diagrams are synthesized
from a series of time-evolving simulation results at the same
Figure 6. Synoptic maps of east-limb white-light pB images at 3 Rs observed by the LASCO-C2 instrument on board the SOHO satellite (top) for CRs 2110 (left)
and 2111 (right), alongside the timing diagrams of simulated plasma number density (105 cm−3) and radial velocity Vr (km s
−1 at 3 Rs at the east-limb longitude of
the Sun in the Earth view.
23
https://cdaweb.gsfc.nasa.gov/index.html
13
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
longitude as the virtual satellite, with a cadence of one result
every 20 hr. The solid orange lines overlaid on the contours
represent the MNLs. They show that the MNLs, spanning
between 30°S and 45°S, are Latter than those in Figure 6. The
low-speed solar wind (Vr < 400 km s
−1
), associated with high
plasma density, is primarily concentrated around the MNLs.
Meanwhile, the fast solar wind (Vr > 550 km s
−1
), accom-
panied by low plasma density, dominates the polar regions
south of 50°S and north of 80°N. The simulated plasma
temperature distributions are positively correlated with radial
solar wind speeds (H. A. Elliott et al. 2012; R. F. Pinto &
A. P. Rouillard 2017; C. X. Li et al. 2018), with the fast solar
wind exhibiting temperatures ranging from 1 to 1.5 MK, while
the slow wind corresponds to temperatures between 0.3 and
0.7 MK. This demonstrates that the simulated temperature
range matches observations.
4. Concluding Remarks
In this work, we extended our recently developed SIP-IFVM
model (H. P. Wang et al. 2025a), an implicit MHD coronal
model constrained by a time-invariant magnetogram, into a
time-evolving coronal model. We adjusted the adiabatic index
γ from 1.05 to 5/3 to better represent the adiabatic process in
coronal simulations. Furthermore, we designed a magnetic
5eld decomposition strategy for time-evolving coronal simula-
tions and applied it to the time-evolving SIP-IFVM coronal
model. Our approach introduces a temporally piecewise-
constant variable to accommodate part of the nonpotential
5eld, ensuring that B1 remains consistently small throughout
the simulations. As a result, the occurrence of nonphysical
negative thermal pressure when deriving thermal pressure
from energy density in low-β regions was effectively mitigated
in time-evolving coronal simulations.
The main contribution of this paper is to demonstrate that,
by employing an implicit algorithm and adopting the novel
extended magnetic 5eld decomposition strategy proposed
herein, we can perform 3D global MHD coronal simulations
to evaluate the long-term evolution of coronal structures in
practical applications, achieving a computational speed more
than 80 times faster than the real-time evolution using only
192 CPU cores. By these means, the time-evolving MHD
coronal model, driven by a series of time-evolving magneto-
grams, can robustly and ef5ciently resolve low-β issues.
Retaining more critical time-evolving information than
commonly used quasi-steady-state coronal models do further
enables it to capture the dynamic features of the corona with
higher 5delity, offering the potential for more accurate
simulations of solar disturbances, such as CME propagation.
We employed the thermodynamic time-evolving MHD
coronal model, incorporating the extended magnetic 5eld
decomposition strategy, to simulate the evolution of the global
corona from the solar surface to 20 Rs during two solar
maximum CRs within an inertial coordinate system. Complet-
ing the 1300 hr evolution of the corona within just 16 hr of
computational time validated its ef5ciency. Given that the
maximum magnetic 5eld strength near the solar surface
occasionally exceeds 40 G and the plasma β reaches a
minimum of approximately 10−3, the model demonstrated
Figure 7. Timing diagram of radial velocity Vr (km s
−1, left), proton number density (105 cm−3 at 3 Rs and 10
3 cm−3 at 20 Rs, middle), and decadic logarithms of
plasma β (top right) and radial magnetic 5eld polarities (bottom right) with “−1” denoting the 5eld pointing to the Sun and “1” for the 5eld directed away from the
Sun, as observed by a virtual satellite located at 3 Rs (top) and 20 Rs (bottom), positioned at 60° behind Earth in longitude and at the same latitude. The black solid
lines and black square symbols represent the simulated variables, while the gray solid lines indicate the variables derived from Wind satellite observations, which
have been mapped from 1 au to 20 Rs following the ballistic propagation.
14
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
suf5cient numerical stability for most solar–terrestrial simula-
tion requirements. Furthermore, the results of the time-
evolving simulation basically reproduced remote EUV and
pB observations and captured in situ measurements mapped
from 1 au to 20 Rs, demonstrating its capability to simulate
complex time-evolving coronal structures during solar
maximum.
Given that the computational time for simulating 1300 hr of
physical time is no more than 16 hr using 192 CPU cores, it is
practical to conduct faster-than-real-time CME simulations
from the solar surface to 1 au based on this work. We will
re5ne the mesh and extend the coronal model to 1 au,
establishing a 3D implicit time-evolving MHD Sun-to-Earth
model chain. Additionally, we will use an observation-based
Lux rope to trigger a realistic CME event in the time-evolving
solar–terrestrial MHD model, validating this new generation of
CME simulations. Furthermore, this model will be used to
enhance and streamline the daily Sun-to-Earth forecasting
process and to investigate the dynamic interaction between the
solar wind and the magnetosphere of planets such as Jupiter
and Saturn.
Although this fully implicit time-evolving MHD solar
coronal model, integrated with the extended magnetic 5eld
decomposition strategy, offers signi5cant advantages and is a
promising tool for timely and accurate simulations of the time-
evolving corona in practical space weather forecasting, there is
still considerable room for further improvement. Synchronized
magnetograms (L. Upton & D. H. Hathaway 2014; J. Linker
et al. 2024) are needed to address the limitation of current
synoptic magnetographs, where magnetic 5elds at different
longitudes are observed at different times, and to improve the
accuracy of reproducing inner-boundary magnetic 5eld when
active regions suddenly emerge on the Sun’s farside. More
accurate measurements of the photospheric magnetic 5elds in
the polar regions are required to reproduce more realistic
coronal structures. Additionally, incorporating more physically
consistent heating source terms is necessary to better simulate
coronal heating and solar wind acceleration during time-
evolving simulations. Moreover, self-consistently simulating
the formation and evolution of coronal eruptions is crucial for
enhancing the reliability of space weather forecasting using
numerical models.
In our future work, we plan to incorporate surface Lux
transport models, such as the advective Lux transport model
(L. Upton & D. H. Hathaway 2014), along with Solar Orbiter
Polarimetric and Helioseismic Imager data (P. Loeschl et al.
2024), horizontal velocities inferred from observational data
using the time–distance helioseismology method (J. Zhao
et al. 2012; M. S. Yalim et al. 2017), and arti5cial intelligence-
generated data from STEREO EUV observations (H.-J. Jeong
et al. 2020, 2022), to advect the radial magnetic 5eld with
observed Lows and generate a more realistic real-time
magnetic 5eld evolution at the inner boundary of our coronal
model. Additionally, some models of local active regions (e.g.,
Figure 8. Timing diagrams of simulated plasma number density (103 cm−3, top), radial velocity Vr (km s
−1, middle) and temperature (105 K, bottom) at 20 Rs,
corresponding to the same longitude as in Figure 7.
15
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
C. W. Jiang et al. 2016; T. Amari et al. 2018; Z. Zhong et al.
2021) will be integrated into this global MHD coronal model
to generate energized 5elds corresponding to transient eruption
events. Furthermore, we plan to investigate the wave
turbulence-driven heating mechanism (S. R. Cranmer 2010;
S. Schleich et al. 2023; Y. H. Gao et al. 2024) in our time-
evolving coronal model to gain a better understanding of how
both fast and slow solar wind streams are heated and
accelerated.
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 Inter-
nal Funds KU Leuven), G0B5823N and G002523N (WEAVE)
(FWO-Vlaanderen), and 4000145223 (SIDC Data Exploitation
(SIDEX), ESA Prodex). This work also bene5ts from the
National Natural Science Foundation of China (grant Nos.
42030204, 42204155, 42274213, and 42474216), and the
Specialized Research Fund for State Key Laboratories
managed by the Chinese State Key Laboratory of Space
Weather. The resources and services used in this work were
provided by the Flemish Supercomputer Centre (VSC), funded
by the Research Foundation–Flanders (FWO) and the Flemish
Government. This work utilizes 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). This research is (partially) funded
by the BK21 FOUR program of Graduate School, Kyung Hee
University (GS-1-JO-NON-20242364). E.H. is grateful to the
Space Weather Awareness Training Network (SWATNet)
funded by the European Union’s Horizon 2020 research and
innovation program under the Marie Skłodowska-Curie grant
agreement No. 955620. Y.Z. acknowledges funding from
Research Foundation–Flanders FWO under project No.
1256423N.
ORCID iDs
Haopeng Wangaa https://orcid.org/0000-0002-4217-6990
Liping Yangaa https://orcid.org/0000-0003-4716-2958
Stefaan Poedtsaa https://orcid.org/0000-0002-1743-0651
Andrea Laniaa https://orcid.org/0000-0003-4017-215X
Yuhao Zhouaa https://orcid.org/0000-0002-4391-393X
Yuhang Gaoaa https://orcid.org/0000-0002-6641-8034
Luis Linanaa https://orcid.org/0000-0002-4014-1815
Tinatin Baratashviliaa https://orcid.org/0000-0002-
1986-4496
Jinhan Guoaa https://orcid.org/0000-0002-4205-5566
Rong Linaa https://orcid.org/0000-0001-7655-5000
Man Zhangaa https://orcid.org/0000-0003-3000-2819
Wenwen Weiaa https://orcid.org/0000-0001-8495-9179
Yun Yangaa https://orcid.org/0009-0000-5292-713X
Edin Husidicaa https://orcid.org/0000-0002-1349-3663
Hyun-Jin Jeongaa https://orcid.org/0000-0003-4616-947X
Mahdi Naja5-Ziyaziaa https://orcid.org/0009-0008-
0922-3995
Brigitte Schmiederaa https://orcid.org/0000-0003-3364-9183
References
Amari, T., Canou, A., Aly, J.-J., Delyon, F., & Alauzet, F. 2018, Natur,
554, 211
Arge, C. N., Odstrcil, D., Pizzo, V. J., & Mayer, L. R. 2003, in AIP Conf.
Proc. 679, Solar Wind Ten (Melville, NY: AIP), 190
Baker, D. N. 1998, AdSpR, 22, 7
Bourdin, P.-A. 2017, ApJL, 850, L29
Brchnelova, M., Kuźma, B., Perri, B., Lani, A., & Poedts, S. 2022, ApJS,
263, 18
Brchnelova, M., Kuźma, B., Zhang, F., et al. 2023, SunGe, 15, 59
Brueckner, G. E., Howard, R. A., Koomen, M. J., et al. 1995, SoPh, 162, 357
Cash, M. D., Biesecker, D. A., Pizzo, V., et al. 2015, SpWea, 13, 611
Cranmer, S. R. 2010, ApJ, 710, 676
Dere, K. P., Del Zanna, G., Young, P. R., & Landi, E. 2023, ApJS, 268, 52
Dere, K. P., Landi, E., Mason, B. C., Monsignori, F., & Young, P. R. 1997,
A&AS, 125, 149
Detman, T., Smith, Z., Dryer, M., et al. 2006, JGRA, 111, A07102
Einfeldt, B., Munz, C. D., Roe, P. L., & Sjogreen, B. 1991, JCoPh, 92, 273
Elliott, H. A., Henney, C. J., McComas, D. J., Smith, C. W., & Vasquez, B. J.
2012, JGRA, 117, A09102
Endeve, E., Leer, E., & Holzer, T. E. 2003, ApJ, 589, 1040
Feng, X. S. 2020, Magnetohydrodynamic Modeling of the Solar Corona and
Heliosphere (Singapore: Springer)
Feng, X. S., Li, C. X., Xiang, C. Q., et al. 2017, ApJS, 233, 10
Feng, X. S., Liu, X. J., Xiang, C. Q., Li, H. C., & Wei, F. S. 2019, ApJ,
871, 226
Feng, X. S., Lv, J. K., Xiang, C. Q., & Jiang, C. W. 2023, MNRAS, 519, 6297
Feng, X. S., Ma, X. P., & Xiang, C. Q. 2015, JGRA, 120, 10,159
Feng, X. S., Wang, H. P., Xiang, C. Q., et al. 2013, Sci. Sin-Terrae, 43, 912
Feng, X. S., Wang, H. P., Xiang, C. Q., et al. 2021, ApJS, 257, 34
Feng, X. S., Xiang, C. Q., & Zhong, D. K. 2011, Sci. Sin-Terrae, 41, 1
Feng, X. S., Yang, L. P., Xiang, C. Q., et al. 2010, ApJ, 723, 300
Fuchs, F. G., McMurry, A. D., Mishra, S., Risebro, N. H., & Waagan, K. 2010,
JCoPh, 229, 4033
Gao, Y. H., Van Doorsselaere, T., Tian, H., Guo, M. Z., & Karampelas, K.
2024, A&A, 689, A195
Gary, G. A. 2001, SoPh, 203, 71
Godunov, S. K. 1959, Matematičeskij Sbornik, 47, 271, https://hal.science/
hal-01620642
Griton, L. 2018, PhD Thesis, Université de recherche Paris Sciences et Lettres
https://hal.science/tel-01906490v1/5le/These_Lea_Griton_2018.pdf
Groth, C. P. T., De Zeeuw, D. L., Gombosi, T. I., & Powell, K. G. 2000, JGR,
105, 25053
Guo, J. H., Linan, L., Poedts, S., et al. 2023, A&A, 683, A54
Guo, X. C. 2015, JCoPh, 290, 352
Hamada, A., Asikainen, T., Virtanen, I., & Mursula, K. 2018, SoPh, 293, 71
Hayashi, K. 2012, JGRA, 117, A08105
Hayashi, K., Abbett, W. P., Cheung, M. C. M., & Fisher, G. H. 2021, ApJS,
254, 1
Hollweg, J. V. 1978, RvGeo, 16, 689
Howard, R. A., Moses, J. D., Vourlidas, A., et al. 2008, SSRv, 136, 67
Jeong, H.-J., Moon, Y.-J., Park, E., & Lee, H. 2020, ApJL, 903, L25
Jeong, H.-J., Moon, Y.-J., Park, E., Lee, H., & Baek, J.-H. 2022, ApJS, 262, 50
Jiang, C. W., Wu, S. T., Feng, X. S., & Hu, Q. 2016, NatCo, 7, 11522
Kimpe, D., Lani, A., Quintino, T., Poedts, S., & Vandewalle, S. 2005, in Proc.
12th European Parallel Virtual Machine and Message Passing Interface
Conf., ed. D. K. B. Di Martino & J. J. Dongarra (Berlin: Springer), 520
King, J. H., & Papitashvili, N. E. 2005, JGRA, 110, A02104
Koskinen, H. E. J., Baker, D. N., Balogh, A., et al. 2017, SSRv, 212, 1137
Kuźma, B., Brchnelova, M., Perri, B., et al. 2023, ApJ, 942, 31
Lani, A., Quintino, T., Kimpe, D., et al. 2005, in Computational Science –
ICCS 2005, ed. V. S. Sunderan et al. (Berlin: Springer), 281
Lani, A., Villedieu, N., Bensassi, K., et al. 2013, in 21st AIAA Computational
Fluid Dynamics Conf. (Reston, VA: American Institute of Aeronautics and
Astronautics) AIAA 2013–2589
Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, SoPh, 275, 17
Li, C. X., Feng, X. S., Xiang, C. Q., et al. 2018, ApJ, 867, 42
16
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.
Linan, L., Regnault, F., Perri, B., et al. 2023, A&A, 675, A101
Linker, J., Downs, C., Caplan, R., et al. 2024, in EGU General Assembly 2024
(Munich: European Geosciences Union) EGU24–4200
Lionello, R., Downs, C., Mason, E. I., et al. 2023, ApJ, 959, 77
Lionello, R., Linker, J. A., & Mikić, Z. 2008, ApJ, 690, 902
Loeschl, P., Valori, G., Hirzberger, J., et al. 2024, A&A, 681, A59
Mikić, Z., Linker, J. A., Schnack, D. D., Lionello, R., & Tarditi, A. 1999,
PhPl, 6, 2217
Owens, M. J., Lockwood, M., & Riley, P. 2017, NatSR, 7, 41548
Parker, E. N. 1963, Interplanetary Dynamical Processes (New York:
Interscience)
Perri, B., Brun, A. S., Réville, V., & Strugarek, A. 2018, JPlPh, 84, 765840501
Perri, B., Kuźma, B., Brchnelova, M., et al. 2023, ApJ, 943, 124
Perri, B., Leitner, P., Brchnelova, M., et al. 2022, ApJ, 936, 19
Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, SoPh, 275, 3
Petrie, G. J. D., Canou, A., & Amari, T. 2011, SoPh, 274, 163
Pinto, R. F., & Rouillard, A. P. 2017, ApJ, 838, 89
Poedts, S., Lani, A., Scolini, C., et al. 2020, JSWSC, 10, 57
Pomoell, J., & Poedts, S. 2018, JSWSC, 8, A35
Powell, K. G., Roe, P. L., Linde, T. J., Gombosi, T. I., & de Zeeuw, D. L.
1999, JCoPh, 154, 284
Reale, F. 2014, LRSP, 11, 4
Reiss, M. A., Muglach, K., Mason, E., et al. 2024, ApJS, 271, 6
Réville, V., Velli, M., Panasenco, O., et al. 2020, ApJS, 246, 24
Rosner, R., Tucker, W., & Vaiana, G. 1978, ApJ, 220, 643
Samara, E., Pinto, R. F., Magdaleni, J., et al. 2021, A&A, 648, A35
Schleich, S., Boro Saikia, S., Ziegler, U., Güdel, M., & Bartel, M. 2023, A&A,
672, A64
Tanaka, T. 1995, JGRA, 100, 12057
Upton, L., & Hathaway, D. H. 2014, ApJ, 780, 5
van der Holst, B., Sokolov, I. V., Meng, X., et al. 2014, ApJ, 782, 81
Wang, H. P., Guo, J. H., Yang, L. P., et al. 2025a, A&A, 693, A257
Wang, H. P., Poedts, S., Lani, A., et al. 2025b, A&A, 694, A234
Wang, H. P., Xiang, C. Q., Liu, X. J., Lv, J. K., & Shen, F. 2022a, ApJ,
935, 46
Wang, H. P., Zhao, J. M., Lv, J. K., & Liu, X. J. 2022b, Chin. J. Geophys.,
65, 2779
Wang, Y., Feng, X. S., & Xiang, C. Q. 2019, CF, 179, 67
Wang, Y. M. N. R., Sheeley, J., & Rich, N. B. 2007, ApJ, 658, 1340
Wu, K. L., & Shu, C. W. 2019, NuMat, 142, 995
Xia, C., Chen, P. F., Keppens, R., & van Marle, A. J. 2011, ApJ, 737,
27
Xia, C., Teunissen, J., Mellah, I. E., Chané, E., & Keppens, R. 2018, ApJS,
234, 30
Yalim, M. S., Pogorelov, N., & Liu, Y. 2017, JPhCS, 837, 012015
Yang, L. P., Feng, X. S., Xiang, C. Q., et al. 2012, JGRA, 117, A08110
Yeates, A., Amari, T., Contopoulos, I., et al. 2018, SSRv, 214, 99
Zhao, J., Couvidat, S., Bogart, R. S., et al. 2012, SoPh, 275, 375
Zhong, Z., Guo, Y., & Ding, M. 2021, NatCo, 12, 2734
Zhou, Y. H., Ruan, W. Z., Xia, C., & Keppens, R. 2021, A&A, 648, A29
17
The Astrophysical Journal Supplement Series, 278:59 (17pp), 2025 June Wang et al.