An Entropy-stable Ideal EC-GLM-MHD Model for the Simulation of the Three-dimensional Ambient Solar Wind

Modeling complex structures of solar wind is one of the important aspects in the study of space weather, but laboratory experiments are impractical to carry out as they surpass the capabilities of our technology. Numerical simulation offers us an opportunity to study the evolutionary processes of solar wind that cannot be seen in observations and experiments. The solar wind plasma flow is typically ionized, compressible, and often supersonic and has essentially infinite conductivity (Goldsmith 1970); thus, we can simulate the solar wind plasma flow by solving the ideal MHD equations. The current status of three-dimensional (3D) MHD simulation is comprehensively reviewed in Feng (2020a), including adaptive mesh refinement (AMR) and data-driven MHD modeling for space weather within the framework of cell-centered finite-volume methods. Due to the nonconvexity and the nonlinear hyperbolicity of the ideal compressible MHD equations, discontinuous solutions may develop in the presence of shocks or contact discontinuities even from smooth initial data. Meanwhile, there are many complicated waves included in the solution of the Riemann problem for the MHD system, such as rarefaction waves, overcompressible shocks, compound waves, and the ordinary shocks (Brio & Wu 1988; Miyoshi & Kusano 2005). Therefore, it is a merit to single out the physically relevant solution and obtain stable numerical results for the evolution of the ambient solar wind, particularly for the basic requirements of developing numerical simulation codes with high accuracy and robustness. Fortunately, extra admissibility criteria, such as entropy-stable property, can be supplemented to obtain a unique, stable, and physically relevant weak solution of conservation laws (e.g., Kružkov 1970; Lax 1971; Mock 1980; Harten 1983; Tadmor 1987, 2003; Dafermos 2000). Entropy-stable numerical schemes satisfy a discrete version of the entropy condition and precisely obey the second law of thermodynamics. Thus, constructing and utilizing entropy-stable methods to model large-scale structured solar wind is significant and highly desirable. The focus of the current paper is to show the implementation of the robust entropy-stable numerical scheme for the study of the ambient solar wind from the Sun to Earth in a six-component grid system, by advancing the entropy-consistent nine-wave magnetic field divergence diminishing ideal MHD (i.e., entropy-consistent generalized Lagrangian multiplier MHD; EC-GLM-MHD) equations with the strong stability-preserving (SSP) Runge–Kutta method (Gottlieb 2005), and the combination of the non-negativity-preserving reconstruction and the solenoidality-preserving reconstruction described in Feng et al. (2017, 2019) and Li et al. (2021).

The common numerical methods for an ideal MHD system are not directly coupled with the universally valid second law of thermodynamics and the divergence-free constraint discretization, which may cause spurious numerical oscillations, and/or negative density and/or pressure, leading to unphysical solutions in the computation. Thus, in designing the numerical scheme for the ideal MHD equations, additional admissibility criteria, such as the divergence-free condition of the magnetic field and the universally valid second law of thermodynamics, should be taken into account to single out the physically relevant solutions. The divergence-free condition (also called the involution condition, Barth 2006; Helzel et al. 2013; or solenoidal condition, Jackson & Zia 1977), ∇ · B = 0, means that the magnetic field lines only present in loops or extend to infinity. In extreme environments, numerical discretization errors can cause a nonphysical intermediate state within a numerical discontinuity profile, leading to the violation of the divergence constraint and frequently triggering severe stability problems. Additionally, in contrast to other variables like density, velocity, and energy, a change in magnetic field strength will change the magnetic field morphology and cause serious instability in the numerical simulation. Thus, the development of divergence cleaning numerical approaches never stops, including the nonconservative term correction, the projection approach of Brackbill & Barnes (1980), the constrained transport method introduced by Evans & Hawley (1988), the eight-wave solution method designed by Powell et al. (1999), and the generalized Lagrange multiplier (GLM) hyperbolic divergence error cleaning technique prescribed by Dedner et al. (2002). In the present paper, we particularly focus on the effective built-in hyperbolic divergence cleaning method because it is highly efficient and parallelizable. By introducing a scalar quantity coupled with the momentum and magnetic induction equations, the magnetic field divergence errors can be moved away from the places where they are generated.

Thermodynamics describes the conversion between different energies, explains how and why physical systems behave, and decides how a system cannot behave. For more details, refer to the literature (e.g., Vogel 2002). The second law of thermodynamics must be obeyed for any macroscopic scale in nature, that is to say, the entropy of a closed physical system conserves or tends to increase over time, and never shrinks. The classical ideal MHD model is not directly coupled with the second law of thermodynamics. In a numerical simulation, fluid only equipped with the conservations of mass, momentum, and energy can display strange behaviors and lead to nonphysical phenomenon (Merriam 1989). Fortunately, many recent researches tried to analyze the second law of thermodynamics from a mathematical viewpoint and deduced the entropy-conservative flux function. Consider a general one-dimensional (1D) system of nonlinear conservation laws U _t + f ( U )_x = 0, with U and f being the vector of conserved variables and the flux; the spontaneous formation of shock discontinuities brings in a challenge when designing a numerical scheme. In this situation, the entropy condition plays an important role to single out a physically relevant weak solution in both the theory and the design of the numerical scheme (Lax 1973; Smoller 1983; Dafermos 2000; Tadmor 2003). It is necessary for the variables U to satisfy the additional admissibility equality and inequality constraints: S( U )_t + F( U )_x = 0 for the smooth solution and S( U )_t + F( U )_x ≤ 0 for the discontinuous solution with all admissible entropy pairs (S( U ), F( U )) (Tadmor 2003). So, the second law of thermodynamics can be built into the partial differential equation system in a mathematically consistent way with a special choice of entropy pair (S( U ), F( U )) satisfying the compatibility conditions above. This auxiliary entropy condition provides guidelines for the numerical scheme to single out a physically relevant solution among the many possible weak solutions (Tadmor 2003). A numerical scheme that satisfies the strict equality in discrete form is an entropy-conservative scheme, which is only possible for the smooth solution. However, a scheme that makes the strict inequality a reality is regarded as an entropy-stable scheme, owing to the concentrations along shock discontinuities.

Constructing an entropy-conservative or entropy-stable scheme is a promising method to obtain a more realistic approximate solution and then enhance the robustness of the numerical scheme. Therefore, many researches have been done on constructing entropy-consistent numerical approximations for nonlinear hyperbolic systems, including general numerical frameworks for systems of conservation laws (e.g., Tadmor 1987, 2003, 2016; LeFloch & Rohde 2000; Fjordholm et al. 2012; Fisher & Carpenter 2013; Fisher et al. 2013), scalar conservation laws (e.g., Jiang & Shu 1994; Gassner 2013), shallow water equations (e.g., Fjordholm et al. 2011; Andreas Hiltebrand 2016; Winters & Gassner 2016a; Wintermeyer et al. 2017), Euler equations (e.g., Ismail & Roe 2009; Chandrashekar 2013; Hiltebrand & Mishra 2014; Fjordholm & Ray 2016; Gassner et al. 2016; Ray et al. 2016; Chen & Shu 2017), Navier–Stokes equations (e.g., Carpenter et al. 2014; Parsani et al. 2016; Gassner et al. 2018), ideal MHD models (e.g., Barth 1999; Waagan et al. 2011; Rossmanith 2013; Chandrashekar & Klingenberg 2016; Derigs et al. 2016, 2018; Winters & Gassner 2016b; Winters et al. 2017; Liu et al. 2018), and relativistic MHD equations (e.g., Bohm et al. 2020; Duan & Tang 2020; Wu & Shu 2020). In the following, some recent entropy-stable schemes for MHD equations are briefly outlined.

Nowadays, there are two most popular methods explicitly associated with the property of entropy stability in the area of solving MHD systems, i.e., the discontinuous Galerkin (DG) approach and the finite-volume method. For the entropy-stable DG approach, there are a batch of recent papers emerging. Motivated by the entropy-stable DG approach proposed for systems of conservation laws (Chen & Shu 2017), Liu et al. (2018) extended this scheme to the symmetrizable version of MHD equations with some nonconservative terms proportional to the divergence-free constraint added, which demands special care to ensure the compatibility with the entropy stability and the suitable DG discretization for these nonconservative source terms. At almost the same time, Bohm et al. (2017) constructed and analyzed a nodal DG method that obeys a discrete version of the entropy for the resistive MHD equations with a built-in GLM divergence cleaning mechanism. From several numerical examples, they demonstrated that the high-order-accuracy entropy conservation/stability method with the built-in GLM divergence cleaning mechanism is more robust than the standard high-order DG method without the property of entropy stability for the resistive MHD equations. Afterward, the theory of this novel entropy-stable nodal DG method for the resistive GLM-MHD equations was analyzed and proved in curvilinear unstructured hexahedral meshes (Bohm et al. 2020), and the theoretical analyses are validated through several numerical results. Other similar studies of the entropy-stable DG method for resistive MHD have been carried out (e.g., Duan & Tang 2020; Wu & Shu 2020). In the finite-volume frame, the numerical schemes satisfying the property of entropy stability are also very promising. Winters & Gassner (2016b) derived an affordable entropy-conserving numerical flux function that discretely guarantees the conservation of the entropy for the MHD equations with a particular source term (Janhunen 2000). When a shock is present or the mesh is coarse, some dissipation, such as the scalar artificial dissipation term or entropy-variable-based matrix dissipation term, has to be added to the entropy-conservative numerical flux to construct an entropy-stable scheme. Based on the spatial reconstruction techniques applied to the scaled entropy variables, Chandrashekar & Klingenberg (2016) constructed a high-accuracy entropy-stable finite-volume method for the equations of ideal compressible MHD. Through some standard MHD tests on very fine meshes, the scheme demonstrated that it can treat discontinuous solutions physically and robustly. A nonlinear entropy-stable numerical MHD solver built on the modified ideal MHD system with the Janhunen source term was proposed by Derigs et al. (2016) and was extended into high-order, multiphysics, and multiscale simulation by using state-variable reconstruction techniques on the AMR grid. Their scheme treated the magnetic field divergence-free constraint in a similar manner to the hyperbolic divergence cleaning mechanism. In order to construct an ideal GLM-MHD system consistent with the continuous entropy analysis, recent research considered the Powell nonconservative terms and the nonconservative GLM terms on the right of the MHD equations (Derigs et al. 2018). Meanwhile, they proposed an entropy-stable scheme that obeys the second law of thermodynamics in finite-volume discretization. Thus, motivated by the merits of the entropy-stable scheme presented by former researchers, the current paper tries to explore the numerical study of the ambient solar wind by using the entropy-stable numerical scheme in the finite-volume framework, so as to deal with the high Mach number or low plasma β situations robustly.

Apart from the considerations of the second law of thermodynamics and the divergence constraint, the present paper also takes into account another major robustness issue in the numerical design, i.e., preserving the positivity of pressure. This intention is essential and important for the reason that negative pressure will make the discrete problem ill-posed and cause the breakdown of a modeling (e.g., Zhang & Shu 2010b; Balsara 2012; Xisto et al. 2014; Christlieb et al. 2015b; Fu & Feng 2015; Wu & Shu 2019). Three efforts have been made for the positivity-preserving property in this paper. First, we combine the non-negativity-preserving reconstruction described in Feng et al. (2017) and the solenoidality-preserving least-squares reconstruction used in Li et al. (2021) to ensure that the reconstructed pressures and densities are positive, and the reconstructed magnetic field remains divergence-free to some extent. Second, the current paper adopts the Roe+LF hybrid entropy-stable scheme by automatically adding high amounts of dissipation for the situations at high-pressure gradient regions (Feng et al. 2019), which can balance the robustness and accuracy of the numerical scheme and then play a positive role in preserving the positivity of pressure. Third, we utilize the additional entropy equation to guarantee the positivity of pressure by introducing nonconservative low-pressure limits (Derigs et al. 2016).

The main aim of the current paper is to utilize the robust entropy-stable numerical flux formulation for the physically motivated entropy-consistent GLM-MHD (EC-GLM-MHD) model with an effective built-in divergence cleaning mechanism to simulate the ambient solar wind from the Sun to Earth, accompanied by two simple 1D tests. The remainder of this paper is as follows: Section 2 provides the ideal EC-GLM-MHD solar wind equations and the grid partition. In Section 3, we give the detailed implementation of the entropy-stable scheme, which is the most important part. This section contains two subsections: some basics of entropy, and discretization of the EC-GLM-MHD system. To clarify the organization, Section 3.2 is divided into three subsubsections: entropy-stable numerical flux, discretization of the source terms, and time integration. Some positivity-preserving schemes and the entropy-related positivity-guaranteeing method for pressure are described in Section 4. In Section 5, we give a brief description of the treatments for the projected normal characteristic boundary method at the inner boundary of the corona. Section 6 presents the numerical results of two 1D tests and the large-scale solar wind structures for Carrington Rotation (CR) 2183. In Section 7, the conclusions are made. Finally, some details of the computation for the logarithmic mean are provided in Appendix A. The entries of the discrete symmetric positive-definite matrix ${ \mathcal H }$ for the entropy-stable LF scheme and the average state of the eigenvalues, diagonal scaling matrix, and right eigenvectors for the entropy-stable Roe scheme are presented in Appendices B and C, respectively.

In this paper, the model equations governing solar wind plasma and magnetic field are described as follows:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\boldsymbol{U}}}{\partial t}+{\rm{\nabla }}\cdot {\boldsymbol{F}}=\displaystyle \frac{\partial }{\partial t}\left[\begin{array}{c}\rho \\ \rho {\boldsymbol{v}}\\ E\\ {\boldsymbol{B}}\\ \psi \end{array}\right]\\ +{\rm{\nabla }}\cdot \left[\begin{array}{c}\rho {\boldsymbol{v}}\\ \rho {\boldsymbol{v}}{\boldsymbol{v}}+\left(p+\displaystyle \frac{1}{2}{{\boldsymbol{B}}}^{2}\right){\boldsymbol{I}}-{\boldsymbol{B}}{\boldsymbol{B}}\\ {\boldsymbol{v}}\left(\displaystyle \frac{1}{2}\rho {{\boldsymbol{v}}}^{2}+\displaystyle \frac{\gamma p}{\gamma -1}+{{\boldsymbol{B}}}^{2}\right)-{\boldsymbol{B}}\left({\boldsymbol{v}}\cdot {\boldsymbol{B}}\right)+{c}_{h}\psi {\boldsymbol{B}}\\ {\boldsymbol{v}}{\boldsymbol{B}}-{\boldsymbol{B}}{\boldsymbol{v}}+{c}_{h}\psi {\boldsymbol{I}}\\ {c}_{h}{\boldsymbol{B}}\end{array}\right]={\boldsymbol{S}},\end{array}\end{eqnarray} \tag{ 1 }$$

where the source terms S = S _pow + S _ψ + S _damp + S _sw + S _h , with ${{\boldsymbol{S}}}_{\mathrm{pow}}=-\left({\rm{\nabla }}\cdot {\boldsymbol{B}}\right){\left(0,{\boldsymbol{B}},{\boldsymbol{v}}\cdot {\boldsymbol{B}},{\boldsymbol{v}},0\right)}^{T}$, ${{\boldsymbol{S}}}_{\psi }\,=-{\rm{\nabla }}\psi \cdot {\left({\bf{0}},0,{\boldsymbol{v}}\psi ,0,{\boldsymbol{v}}\right)}^{T}$, ${{\boldsymbol{S}}}_{\mathrm{damp}}=-{\left(0,{\bf{0}},0,{\bf{0}},\tfrac{{c}_{h}^{2}}{{c}_{p}^{2}}\psi \right)}^{T}$, ${{\boldsymbol{S}}}_{\mathrm{sw}}\,=(0,\,-\,\rho \tfrac{{{GM}}_{s}}{{r}^{3}}{\boldsymbol{r}}\,-\,2\rho {\boldsymbol{\Omega }}\,\times \,{\boldsymbol{v}}-\rho {\boldsymbol{\Omega }}\,\times \,\left({\boldsymbol{\Omega }}\,\times \,{\boldsymbol{r}}\right),-\rho {\boldsymbol{v}}\cdot (\tfrac{{{GM}}_{s}}{{r}^{3}}{\boldsymbol{r}}\,+{\boldsymbol{\Omega }}\times ({\boldsymbol{\Omega }}\times {\boldsymbol{r}})),{\bf{0}},0{)}^{T}$ ${{\boldsymbol{S}}}_{h}={\left(0,{S}_{M},\mathrm{0,0},{S}_{E},{\bf{0}},0\right)}^{T}$. The meanings and values of all the related variables can be seen in previous papers (e.g., Feng et al. 2017; Li et al. 2018). We will not repeat them here. The new total energy is $E=\tfrac{p}{\gamma -1}\,+\tfrac{1}{2}\rho {{\boldsymbol{v}}}^{2}+\tfrac{1}{2}{{\boldsymbol{B}}}^{2}+\tfrac{1}{2}{\psi }^{2}$, containing the kinetic energy $\tfrac{1}{2}\rho {{\boldsymbol{v}}}^{2}$, the thermal energy $\tfrac{p}{\gamma -1}$, the magnetic energy $\tfrac{1}{2}{{\boldsymbol{B}}}^{2}$, and the divergence-correcting field energy $\tfrac{1}{2}{\psi }^{2}$, similar to the magnetic field energy. The ideal EC-GLM-MHD equations are not to strictly enforce the divergence-free condition ∇ · B = 0, but rather to evolve toward the divergence-free state. Compared with the definition of the generalized Lagrangian multiplier ψ in Dedner et al. (2002), the scalar field ψ here not only acts as the role of propagating the divergence errors to the domain boundaries but also can reduce the complexity of the thermodynamic analysis, if nonconservative terms proportional to the divergence of magnetic field and the gradient of the Lagrange multiplier ψ are included (e.g., Powell et al. 1999; Derigs et al. 2018). Another important modification of our solar wind model is that the total energy E contains the divergence-correcting field energy term $\tfrac{1}{2}{\psi }^{2}$. That is because the ψ field coupled into the magnetic induction equation alters the energy and the magnetic field stored therein. The modifications in energy conservation equation and additional nonconservative source terms on the right-hand side of Equation (1) can also be attributed to including the magnetic field correction energy in the total energy. The cleaning velocity, c_h , is an important parameter determining the velocity with which the divergence error of the magnetic field is propagated out of the computational region. From this point, the effect of the hyperbolic divergence error cleaning mechanism is best if the cleaning speed c_h is as high as possible. However, the parameter c_h should not be increased to arbitrarily high, as it restricts the time step. If we set the ψ waves slower than the fast magnetoacoustic waves, the same Courant–Friedrichs–Lewy (CFL) stability criterion can be used regardless of whether or not the hyperbolic divergence corrections are applied to the solution. This requirement is satisfied if we choose c_h as ${c}_{h}={\lambda }_{\max }-{v}_{\max ,{\rm{\Omega }}}$, with ${v}_{\max ,{\rm{\Omega }}}={\max }_{{\rm{\Omega }}}(| {v}_{r}| ,| {v}_{\theta }| ,| {v}_{\phi }| )$ being the largest fluid speed in the whole computational region Ω. If ${v}_{\max ,{\rm{\Omega }}}=0$, ${c}_{h}={\lambda }_{\max }$, which means the divergence correction being most effective. If ${v}_{\max ,{\rm{\Omega }}}\ne 0$, then ${c}_{h}\lt {\lambda }_{\max }$. Neither of the eigenvalues corresponding to ψ exceeds the maximum of the remaining ones, which guarantees the time step of the computation not being negatively affected.

As we know, the conservation property will be lost owing to the nonconservative source terms added to the ideal MHD model. These nonconservative source terms also render additional technical challenges on the compatibility of the entropy stability and the suitable discretization. However, they are essential for the continuous entropy analysis and the symmetrization of the ideal MHD equations (Liu et al. 2018). Moreover, the nonconservative source terms ensure Galilean invariance for the ideal MHD system. Galilean invariance is an essential physical law for any well-defined system in nonrelativistic physics, which guarantees that the same physical properties are universally valid and generally true in all inertial reference frames. In spherical coordinates, we write the time-dependent ideal EC-GLM-MHD system (1) in an explicit derivative form:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\boldsymbol{U}}}{\partial t}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial {r}^{2}{{\boldsymbol{F}}}_{r}}{\partial r}+\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial \sin \theta {{\boldsymbol{F}}}_{\theta }}{\partial \theta }\\ \quad +\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial {{\boldsymbol{F}}}_{\phi }}{\partial \phi }={\boldsymbol{S}}+{{\boldsymbol{S}}}_{\mathrm{geo}},\end{array}\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray*}&&{\boldsymbol{U}}={\left(\rho ,\rho {v}_{r},\rho {v}_{\theta },\rho {v}_{\phi },E,{B}_{r},{B}_{\theta },{B}_{\phi },\psi \right)}^{T},\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}{{\boldsymbol{F}}}_{r} & = & \left\{\begin{array}{c}\rho {v}_{r}\\ p+\rho {v}_{r}^{2}+\displaystyle \frac{-{B}_{r}^{2}+{B}_{\theta }^{2}+{B}_{\phi }^{2}}{2}\\ \rho {v}_{r}{v}_{\theta }-{B}_{r}{B}_{\theta }\\ \rho {v}_{r}{v}_{\phi }-{B}_{r}{B}_{\phi }\\ {v}_{r}\hat{E}-{B}_{r}\left({\boldsymbol{v}}\cdot {\boldsymbol{B}}\right)+{c}_{h}\psi {B}_{r}\\ {c}_{h}\psi \\ {v}_{r}{B}_{\theta }-{v}_{\theta }{B}_{r}\\ {v}_{r}{B}_{\phi }-{v}_{\phi }{B}_{r}\\ {c}_{h}{B}_{r}\\ \end{array}\right\},\\ {{\boldsymbol{S}}}_{\mathrm{geo}} & = & \left(\begin{array}{c}0\\ \rho \displaystyle \frac{{v}_{\theta }^{2}+{v}_{\phi }^{2}}{r}+\displaystyle \frac{2p}{r}+\displaystyle \frac{{B}_{r}^{2}}{r}\\ {S}_{\mathrm{geo},3}\\ {S}_{\mathrm{geo},4}\\ 0\\ 2{c}_{h}\displaystyle \frac{\psi }{r}\\ \displaystyle \frac{1}{r}({B}_{\theta }{v}_{r}-{B}_{r}{v}_{\theta })+\displaystyle \frac{{c}_{h}\psi }{r}\cot \theta \\ \displaystyle \frac{1}{r}({v}_{r}{B}_{\phi }-{v}_{\phi }{B}_{r})+({B}_{\phi }{v}_{\theta }-{B}_{\theta }{v}_{\phi })\displaystyle \frac{\cot \theta }{r}\\ 0\end{array}\right),\end{array}\end{eqnarray*}$$

with $\hat{E}=\tfrac{1}{2}\rho {{\boldsymbol{v}}}^{2}+\tfrac{\gamma p}{\gamma -1}+{{\boldsymbol{B}}}^{2}$, ${S}_{\mathrm{geo},3}=\tfrac{1}{r}({B}_{r}{B}_{\theta }-\rho {v}_{r}{v}_{\theta })\,+\left(p+\rho {v}_{\phi }^{2}+\tfrac{{B}_{r}^{2}+{B}_{\theta }^{2}-{B}_{\phi }^{2}}{2}\right)\tfrac{\cot \theta }{r}$, and ${S}_{\mathrm{geo},4}=\tfrac{1}{r}({B}_{r}{B}_{\phi }-\rho {v}_{r}{v}_{\phi })\,+({B}_{\theta }{B}_{\phi }-\rho {v}_{\theta }{v}_{\phi })\tfrac{\cot \theta }{r}$. The flux vectors F _θ and F _ϕ in the θ- and ϕ- directions are omitted here. The variables ρ, v , p, B , Ω, r, and t in the above equations are normalized by the characteristic values before the computation and regain dimensions by the end of the computation (e.g., Feng et al. 2010).

The simulation of the ambient solar wind from the Sun to Earth is coded on the six-component spherical coordinate grid system with the whole computational region divided into two parts: spherical shell I from 1 R_s to 25 R_s , and spherical shell II ranging from 25 R_s to 230 R_s (Feng et al. 2010, 2012a, 2012b, 2012c, 2014a, 2014b; Li et al. 2018). The grid resolutions in the θ- and ϕ- directions are Δθ = Δϕ = 15 and Δθ = Δϕ = 1° for spherical shell I and spherical shell II, respectively. However, in the r-direction, the grid resolutions are highly nonuniform with Δr(i) = 0.01R_s if ${r}_{{im}}\lt 1.1{R}_{s}$ (where r_im is the bottom bound of a finite-volume cell (i, j, k) in the r-direction), ${\rm{\Delta }}r(i)=\min \left(0.01/{\mathrm{log}}_{10}(1.09)\times {\mathrm{log}}_{10}({r}_{(i-1)m}),{r}_{(i-1)m}{\rm{\Delta }}\theta \right)$ with Δθ = 15 if $1.1{R}_{s}\leqslant {r}_{{im}}\lt 3.5{R}_{s}$, Δr(i) = r_(i−1)m Δθ with Δθ = 15 if $3.5{R}_{s}\leqslant {r}_{{im}}\lt 25{R}_{s}$, and Δr(i) = r_(i−1)m Δθ with Δθ = 1° if ${r}_{{im}}\geqslant 25{R}_{s}$. For every component of the six-component grid system, each finite-volume cell (i, j, k) is an approximative hexahedron with $[{r}_{{im}},{r}_{{ip}}]\times [{\theta }_{{jm}},{\theta }_{{jp}}]\times [{\phi }_{{km}},{\phi }_{{kp}}]$, and the six bounding faces of the cell are denoted by indices ${im}=i-1/2$ (jm = j − 1/2, km = k − 1/2) and ip = i + 1/2 (jp = j + 1/2, kp = k + 1/2), where i = 1, ⋯ , N_r , j = 1, ⋯ , N_θ , k = 1, ⋯ , N_ϕ , with N_r = 136, N_θ = N_ϕ = 66 for spherical shell I and N_r = 128, N_θ = N_ϕ = 96 for spherical shell II. Starting from ${r}_{1m}=1{R}_{s},{\theta }_{1m}=\tfrac{\pi }{4}-{\delta }_{1}$, and ${\phi }_{1m}=\tfrac{3\pi }{4}-{\delta }_{2}$ with δ₁ = 3Δθ and δ₂ = 3Δϕ, then the grid mesh for each component in spherical shell I comes into being with ${r}_{{im}}={r}_{(i-1)p},{r}_{{ip}}\,={r}_{{im}}+{\rm{\Delta }}r(i)$, θ_jm = θ_(j−1)p , θ_jp = θ_jm + Δθ, ϕ_km = ϕ_(k−1)p , ϕ_kp = ϕ_km + Δϕ. The same applies for spherical shell II. The volume-averaged coordinates are $({\bar{r}}_{i},{\bar{\theta }}_{j},{\phi }_{k})$ with ${\bar{r}}_{i}=3({r}_{{ip}}^{4}-{r}_{{im}}^{4})/(4({r}_{{ip}}^{3}-{r}_{{im}}^{3}))$, ${\bar{\theta }}_{j}=\left(\sin {\theta }_{{jp}}-\sin {\theta }_{{jm}}-\left({\theta }_{{jp}}\cos {\theta }_{{jp}}-{\theta }_{{jm}}\cos {\theta }_{{jm}}\right)\right)/(\cos {\theta }_{{jm}}-\cos {\theta }_{{jp}})$, and ϕ_k = (ϕ_kp + ϕ_km )/2. The coordinates of the six face centers are $({r}_{{im}},{\theta }_{j}^{r},{\phi }_{k})$, $({r}_{{ip}},{\theta }_{j}^{r},{\phi }_{k})$, $({r}_{i}^{\theta },{\theta }_{{jm}},{\phi }_{k})$, $({r}_{i}^{\theta },{\theta }_{{jp}},{\phi }_{k})$, $({r}_{i}^{\phi },{\theta }_{j}^{\phi },{\phi }_{{km}})$, and $({r}_{i}^{\phi },{\theta }_{j}^{\phi },{\phi }_{{kp}})$, with ${r}_{i}^{\theta }={r}_{i}^{\phi }\,=2({r}_{{ip}}^{3}-{r}_{{im}}^{3})/(3({r}_{{ip}}^{2}-{r}_{{im}}^{2}))$, ${\theta }_{j}^{r}={\bar{\theta }}_{j}$, and ${\theta }_{j}^{\phi }=({\theta }_{{jp}}+{\theta }_{{jm}})/2$.

This section gives some concepts related to the entropy conservation and entropy stability and then introduces the entropy-stable discrete form for the ideal EC-GLM-MHD system.

3.1. Some Basics of Entropy

To construct an entropy-stable scheme, it is essential to introduce the entropy function ${ \mathcal S }$ for the ideal EC-GLM-MHD system (e.g., Winters & Gassner 2016b):

$$\begin{eqnarray}&&{ \mathcal S }=-\frac{\rho s}{\gamma -1},\end{eqnarray} \tag{ 3 }$$

where $s(\rho ,p)=\mathrm{ln}p-\gamma \mathrm{ln}\rho =-(\gamma -1)\mathrm{ln}(\rho )-\mathrm{ln}(\varsigma )-\mathrm{ln}(2)$ is the physical entropy per particle with $\varsigma =\tfrac{\rho }{2p}$. In terms of the second law of thermodynamics, the approximation of a numerical method can be divided into two regimes: entropy conservative and entropy stable. For smooth solutions, if the local changes of the discretized entropy are equal to those predicted by the continuous entropy conservation law, the numerical method is entropy conservative, i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{ \mathcal S }+{\rm{\nabla }}\cdot ({\boldsymbol{v}}{ \mathcal S })=0.\end{eqnarray} \tag{ 4 }$$

However, the approximation of a numerical scheme is entropy-stable if the entropy satisfies the entropy inequality

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{ \mathcal S }+{\rm{\nabla }}\cdot ({\boldsymbol{v}}{ \mathcal S })\leqslant 0\end{eqnarray} \tag{ 5 }$$

for discontinuous solution. We also introduce the entropy variables to switch from conserved space to entropy space:

$$\begin{eqnarray*}\begin{array}{rcl}{ \mathcal V }:= \frac{\partial { \mathcal S }}{\partial {\boldsymbol{U}}} & = & \left(\frac{\gamma -s}{\gamma -1}-\varsigma {{\boldsymbol{v}}}^{2},2\varsigma {v}_{r},2\varsigma {v}_{\theta },2\varsigma {v}_{\phi },\right.\\ & & -{\left.2\varsigma ,2\varsigma {B}_{r},2\varsigma {B}_{\theta },2\varsigma {B}_{\phi },2\varsigma \psi \right)}^{T},\end{array}\end{eqnarray*}$$

which contracts the ideal EC-GLM-MHD system into the entropy equality (4) or entropy inequality (5) according to the smoothness of the solution U . A broader introduction to the concepts associated with entropy has been given by previous literature (e.g., Tadmor 1984; Derigs et al. 2018).

3.2. Discretization of the EC-GLM-MHD System

In spherical coordinates, the finite-volume semidiscretization for system (2) can be written as follows (e.g., Feng et al. 2014b; Feng 2020b):

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}{\overline{{\boldsymbol{U}}}}_{i,j,k}^{n} & = & -\displaystyle \frac{3\left({r}_{{ip}}^{2}{{\boldsymbol{F}}}_{r,{ip},j,k}^{n,\mathrm{ES}}-{r}_{{im}}^{2}{{\boldsymbol{F}}}_{r,{im},j,k}^{n,\mathrm{ES}}\right)}{\left({r}_{{ip}}^{3}-{r}_{{im}}^{3}\right)}\\ \ & & -\displaystyle \frac{\sin {\theta }_{{jp}}{{\boldsymbol{F}}}_{\theta ,i,{jp},k}^{n,\mathrm{ES}}-\sin {\theta }_{{jm}}{{\boldsymbol{F}}}_{\theta ,i,{jm},k}^{n,\mathrm{ES}}}{{r}_{i}^{\theta }\left(\cos {\theta }_{{jm}}-\cos {\theta }_{{jp}}\right)}\\ & & -\displaystyle \frac{{\rm{\Delta }}\theta (j)}{{r}_{i}^{\phi }\left(\cos {\theta }_{{jm}}-\cos {\theta }_{{jp}}\right)}\displaystyle \frac{{{\boldsymbol{F}}}_{\phi ,i,j,{kp}}^{n,\mathrm{ES}}-{{\boldsymbol{F}}}_{\phi ,i,j,{km}}^{n,\mathrm{ES}}}{{\rm{\Delta }}\phi (k)}\\ & & +{\widetilde{{\boldsymbol{S}}}}_{\mathrm{pow},i,j,k}+{\widetilde{{\boldsymbol{S}}}}_{\psi ,i,j,k}\\ & & +{\overline{{\boldsymbol{S}}}}_{\mathrm{damp},i,j,k}+{\overline{{\boldsymbol{S}}}}_{\mathrm{sw},i,j,k}+{\overline{{\boldsymbol{S}}}}_{h,i,j,k}+{\overline{{\boldsymbol{S}}}}_{\mathrm{geo},i,j,k},\end{array}\end{eqnarray} \tag{ 6 }$$

where the cell-averaged variables $\overline{{\boldsymbol{U}}}=\tfrac{1}{{V}_{i,j,k}}{\int }_{{V}_{i,j,k}}{\boldsymbol{U}}{dV}$ are evaluated at the centroid point $({\overline{r}}_{i},{\overline{\theta }}_{j},{\phi }_{k})$ of a cell with geometrical center (i, j, k). The explicit expressions of the numerical fluxes F ^n,ES in Equation (6) will be given in Section 3.2.1. The discretization of the source terms $\overline{{\boldsymbol{S}}}$, ${\widetilde{{\boldsymbol{S}}}}_{\mathrm{pow}}$, and ${\widetilde{{\boldsymbol{S}}}}_{\psi }$ in Equation (6) will be presented in Section 3.2.2. Finally, Section 3.2.3 will provide the time integration of the semidiscretization (6).

3.2.1. Entropy-stable Numerical Flux

It is known that the entropy-conserving approximation is true only if the solution remains smooth. It may suffer spurious oscillations in the presence of discontinuities and result in breakdown eventually. For the solar wind flow, the solution is not purely smooth and an unphysical solution will develop if we only utilize the entropy-conserving approximation without proper treatments during the computation. From the second law of thermodynamics, we find that the process of converting magnetic energy and kinetic energy into heat is irreversible, which can be denoted as dissipation. Therefore, additional numerical dissipation should be added to the entropy-conservative flux function in an entropy-consistent way to ensure that the entropy inequality (5) holds in discrete form for the solar wind simulation. In this paper, we utilize the entropy-conservative (EC) flux (8) to construct an entropy-stable (ES) numerical flux (7) by adding a general form of numerical dissipation (Derigs et al. 2018):

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{r}^{\mathrm{ES}}={{\boldsymbol{F}}}_{r}^{\mathrm{EC}}-\displaystyle \frac{1}{2}{\boldsymbol{D}}[\kern-2pt[ {\boldsymbol{U}}]\kern-2pt] ,\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{r}^{\mathrm{EC}}=\left(\begin{array}{c}{\rho }^{\mathrm{ln}}\{\{{v}_{r}\}\}\\ {\rho }^{\mathrm{ln}}\{\{{v}_{r}\}\}{}^{2}+{\overline{p}}_{\mathrm{tot}}-\{\{{B}_{r}\}\}{}^{2}\\ {\rho }^{\mathrm{ln}}\{\{{v}_{r}\}\}\{\{{v}_{\theta }\}\}-\{\{{B}_{r}\}\}\{\{{B}_{\theta }\}\}\\ {\rho }^{\mathrm{ln}}\{\{{v}_{r}\}\}\{\{{v}_{\phi }\}\}-\{\{{B}_{r}\}\}\{\{{B}_{\phi }\}\}\\ {F}_{r,5}^{* }\\ {c}_{h}\{\{\psi \}\}\\ \{\{{v}_{r}\}\}\{\{{B}_{\theta }\}\}-\{\{{v}_{\theta }\}\}\{\{{B}_{r}\}\}\\ \{\{{v}_{r}\}\}\{\{{B}_{\phi }\}\}-\{\{{v}_{\phi }\}\}\{\{{B}_{r}\}\}\\ {c}_{h}\{\{{B}_{r}\}\}\end{array}\right),\end{eqnarray} \tag{ 8 }$$

with ${\overline{p}}_{\mathrm{tot}}=\tilde{p}+\tfrac{1}{2}\left(\{\{{B}_{r}^{2}\}\}+\{\{{B}_{\theta }^{2}\}\}+\{\{{B}_{\phi }^{2}\}\}\right),$

$$\begin{eqnarray*}\begin{array}{rcl}{F}_{r,5}^{* } & = & {F}_{r,1}^{* }\left[\displaystyle \frac{1}{2(\gamma -1){\varsigma }^{\mathrm{ln}}}-\displaystyle \frac{1}{2}\left(\{\{{v}_{r}^{2}\}\}+\{\{{v}_{\theta }^{2}\}\}\right.\right.\\ & & \left.\left.+\{\{{v}_{\phi }^{2}\}\}\right)\right]\\ & & +{F}_{r,2}^{* }\{\{{v}_{r}\}\}+{F}_{r,3}^{* }\{\{{v}_{\theta }\}\}+{F}_{r,4}^{* }\{\{{v}_{\phi }\}\}\\ & & +{F}_{r,6}^{* }\{\{{B}_{r}\}\}+{F}_{r,7}^{* }\{\{{B}_{\theta }\}\}+{F}_{r,8}^{* }\{\{{B}_{\phi }\}\}\\ & & +{F}_{r,9}^{* }\{\{\psi \}\}\\ & & -\displaystyle \frac{1}{2}\left(\{\{{v}_{r}{B}_{r}^{2}\}\}+\{\{{v}_{r}{B}_{\theta }^{2}\}\}+\{\{{v}_{r}{B}_{\phi }^{2}\}\}\right)\\ & & +\{\{{B}_{r}\}\}\left(\{\{{v}_{r}{B}_{r}\}\}+\{\{{v}_{\theta }{B}_{\theta }\}\}\right.\\ & & \left.+\{\{{v}_{\phi }{B}_{\phi }\}\}\right)-{c}_{h}\{\{{B}_{r}\psi \}\},\end{array}\end{eqnarray*}$$

$\tilde{p}=\tfrac{\{\{\rho \}\}}{2\{\{\varsigma \}\}}$, $\{\{\cdot \}\}=\tfrac{(\cdot )+(\cdot )}{2}$, ${(\cdot )}^{\mathrm{ln}}=\tfrac{[\kern-2pt[ \cdot ]\kern-2pt] }{[\kern-2pt[ \mathrm{ln}(\cdot )]\kern-2pt] }$, 〚 · 〛 = ( · )_R − ( · )_L , and D a dissipation matrix. For details on the logarithmic mean, refer to Appendix A. The entropy-conservative numerical flux (8) conserves the discrete entropy and is consistent with the physical flux satisfying ${{\boldsymbol{F}}}_{\kappa }^{\mathrm{EC}}({\boldsymbol{U}},{\boldsymbol{U}})={{\boldsymbol{F}}}_{\kappa }$ (κ = r, θ, ϕ). Meanwhile, it avoids the problematic discretizations for nonconservative terms (e.g., Winters et al. 2017; Derigs et al. 2018). In order to construct an entropy-stable approximation for the EC-GLM-MHD system, the dissipation term in (7) should be ensured to make negative contributions in the discrete entropy inequality (5). Different forms of the dissipation term correspond to different approximate Riemann solvers. For example, we can construct entropy-stable HLL-family fluxes if the left and right wave speed bounds λ_L and λ_R are properly chosen (Bouchut et al. 2007; Schmidtmann & Winters 2017). We can also obtain Lax–Friedrichs (LF) or Roe entropy-stable numerical fluxes through specific discretization. In our solar wind simulation, we adopt a hybrid entropy-stable approximation that continuously combines LF and Roe entropy-stable numerical fluxes to provide more dissipation in regions containing discontinuities and less dissipation with the flow being well resolved.

If we choose a simple choice of dissipation matrix D as

$$\begin{eqnarray*}&&{{\boldsymbol{D}}}_{\mathrm{LF}}=\left|{\lambda }_{\max }^{\mathrm{global}}\right|{\boldsymbol{I}},\end{eqnarray*}$$

where $| {\lambda }_{\max }^{\mathrm{global}}| =\max \left(| {\lambda }_{i}| \right)$ (i = 1, 2, ⋯ ,9) is the global largest eigenvalue of the ideal EC-GLM-MHD system, the dissipation term can be written as

$$\begin{eqnarray*}&&\begin{array}{l}\displaystyle \frac{1}{2}{{\boldsymbol{D}}}_{\mathrm{LF}}[\kern-2pt[ {\boldsymbol{U}}]\kern-2pt] =\displaystyle \frac{1}{2}\left|{\lambda }_{\max }^{\mathrm{global}}\right|{\boldsymbol{I}}[\kern-2pt[ {\boldsymbol{U}}]\kern-2pt] \\ =\displaystyle \frac{1}{2}\left|{\lambda }_{\max }^{\mathrm{global}}\right|\hat{{ \mathcal H }}[\kern-2pt[ { \mathcal V }]\kern-2pt] \simeq \displaystyle \frac{1}{2}\left|{\lambda }_{\max }^{\mathrm{global}}\right|{ \mathcal H }[\kern-2pt[ { \mathcal V }]\kern-2pt] ,\end{array}\end{eqnarray*}$$

where the entropy Jacobian matrix $\hat{{ \mathcal H }}=\tfrac{\partial {\boldsymbol{U}}}{\partial { \mathcal V }}$ is the transition matrix from the conserved space to the entropy space. Instead of incorporating the jump in conservative variables, the entropy-stable dissipation term incorporates the jump in entropy variables to guarantee negative contributions in the entropy inequality and ensure entropy stability (e.g., Barth 1999; Winters & Gassner 2016b). The entropy Jacobian transition matrix $\hat{{ \mathcal H }}$ can be found easily in continuous space, but it is highly nontrivial to discretize the matrix in a numerical scheme. Moreover, in the presence of strong pressure jump, the commonly used arithmetic mean may lead to an unphysical state (Derigs et al. 2017). So, it is important to seek a specific averaging for the elements of the discrete entropy Jacobian ${ \mathcal H }$ to ensure the relation $[\kern-2pt[ {\boldsymbol{U}}]\kern-2pt] ={ \mathcal H }[\kern-2pt[ { \mathcal V }]\kern-2pt]$ whenever possible. Solving the solutions of the 81 equations of the relation $[\kern-2pt[ {\boldsymbol{U}}]\kern-2pt] ={ \mathcal H }[\kern-2pt[ { \mathcal V }]\kern-2pt]$ step by step, one can obtain the entries of the discrete symmetric positive-definite matrix ${ \mathcal H }$. They are presented in Appendix B. This particular averaging of the matrix ${ \mathcal H }$ holds asymptotically for the jump in total energy ${\left([\kern-2pt[ {\boldsymbol{U}}]\kern-2pt] \right)}_{5}\approx {\left({ \mathcal H }[\kern-2pt[ { \mathcal V }]\kern-2pt] \right)}_{5}$ and holds strictly for other terms ${\left([\kern-2pt[ {\boldsymbol{U}}]\kern-2pt] \right)}_{i}={\left({ \mathcal H }[\kern-2pt[ { \mathcal V }]\kern-2pt] \right)}_{i}\,(i\ne 5)$, which is required to construct a discrete entropy-stable dissipation term (Derigs et al. 2016). Finally, we obtain the LF entropy-stable numerical flux:

$$\begin{eqnarray*}&&{{\boldsymbol{F}}}_{r}^{\mathrm{ES},\mathrm{LF}}={{\boldsymbol{F}}}_{r}^{\mathrm{EC}}-\displaystyle \frac{1}{2}\left|{\lambda }_{\max }^{\mathrm{global}}\right|{ \mathcal H }[\kern-2pt[ { \mathcal V }]\kern-2pt] .\end{eqnarray*}$$

Since the global nature of the LF flux and the wave speed adopted is the maximum eigenvalue in the computational domain, the LF scheme is quite dissipative (especially for slow waves). However, Rusanov (1962) showed a less diffusive and stable scheme by using a local maximum wave speed of the solution, instead of a global wave speed measure. The dissipation term is called the local Lax–Friedrichs (LLF) stabilization term. The form of the dissipation term reads

$$\begin{eqnarray*}&&\displaystyle \frac{1}{2}{{\boldsymbol{D}}}_{\mathrm{LLF}}[\kern-2pt[ {\boldsymbol{U}}]\kern-2pt] =\displaystyle \frac{1}{2}\left|{\lambda }_{\max }^{\mathrm{local}}\right|{ \mathcal H }[\kern-2pt[ { \mathcal V }]\kern-2pt] ,\end{eqnarray*}$$

where $| {\lambda }_{\max }^{\mathrm{local}}| =\max \left(| {\lambda }_{\max ,R}| ,| {\lambda }_{\max ,L}| \right)$. The final entropy-stable LLF numerical flux is

$$\begin{eqnarray*}&&{{\boldsymbol{F}}}_{r}^{\mathrm{ES},\mathrm{LLF}}={{\boldsymbol{F}}}_{r}^{\mathrm{EC}}-\displaystyle \frac{1}{2}\left|{\lambda }_{\max }^{\mathrm{local}}\right|{ \mathcal H }[\kern-2pt[ { \mathcal V }]\kern-2pt] .\end{eqnarray*}$$

The scalar dissipation term will always add dissipation on interfaces separating zones of different densities and cannot resolve contact discontinuities exactly. Besides, unlike a selective dissipative term, this scalar dissipation scheme usually introduces much more dissipation than needed and cannot dissipate selectively based on the waves.

When the dissipation matrix is ${{\boldsymbol{D}}}_{\mathrm{Roe}}={\boldsymbol{R}}\,\left|{\boldsymbol{\Lambda }}\right|{{\boldsymbol{R}}}^{-1}$, we call it the Roe scheme. This scheme is less diffusive compared with LF or LLF methods. Thereinto, R is the matrix consisting of the columns of right eigenvectors, and Λ is the diagonal matrix of the eigenvalues for the Jacobian matrix corresponding to the conservative variables. To cut down the amount of calculation, a particular positive diagonal matrix Z related to the right eigenvector matrix R and the entropy Jacobian ${ \mathcal H }$ is introduced in the entropy-stable scheme, which satisfies

$$\begin{eqnarray}&&{ \mathcal H }={\boldsymbol{R}}{\boldsymbol{Z}}{{\boldsymbol{R}}}^{T}.\end{eqnarray} \tag{ 9 }$$

Fortunately, the ideal EC-GLM-MHD system with nonconservative terms is symmetrizable (Godunov 1972; Barth 1999), which ensures the existence of the entropy scaling (Derigs et al. 2016). Given the relation, the dissipation term of the Roe scheme can be rewritten as

$$\begin{eqnarray*}\begin{array}{rcl}\displaystyle \frac{1}{2}{{\boldsymbol{D}}}_{\mathrm{Roe}}[\kern-2pt[ {\boldsymbol{U}}]\kern-2pt] & = & \displaystyle \frac{1}{2}{\boldsymbol{R}}\left|{\boldsymbol{\Lambda }}\right|{{\boldsymbol{R}}}^{-1}[\kern-2pt[ {\boldsymbol{U}}]\kern-2pt] \simeq \displaystyle \frac{1}{2}{\boldsymbol{R}}\left|{\boldsymbol{\Lambda }}\right|{{\boldsymbol{R}}}^{-1}{ \mathcal H }[\kern-2pt[ { \mathcal V }]\kern-2pt] \\ & = & \displaystyle \frac{1}{2}{\boldsymbol{R}}\left|{\boldsymbol{\Lambda }}\right|{{\boldsymbol{R}}}^{-1}{\boldsymbol{R}}{\boldsymbol{Z}}{{\boldsymbol{R}}}^{T}[\kern-2pt[ { \mathcal V }]\kern-2pt] =\displaystyle \frac{1}{2}{\boldsymbol{R}}\left|{\boldsymbol{\Lambda }}\right|{\boldsymbol{Z}}{{\boldsymbol{R}}}^{T}[\kern-2pt[ { \mathcal V }]\kern-2pt] .\end{array}\end{eqnarray*}$$

The average form of the right eigenvector matrix is ${\boldsymbol{R}}=\left[{{\boldsymbol{r}}}_{+f},{{\boldsymbol{r}}}_{+a},{{\boldsymbol{r}}}_{+s},{{\boldsymbol{r}}}_{+\psi },{{\boldsymbol{r}}}_{E},{{\boldsymbol{r}}}_{-\psi },{{\boldsymbol{r}}}_{-s},{{\boldsymbol{r}}}_{-a},{{\boldsymbol{r}}}_{-f}\right]$. As the MHD system has nine traveling waves, we can also call it the nine-wave formulation denoting the ideal EC-GLM-MHD system. When setting ψ = c_h = 0, one immediately gets an eight-wave scheme. Both GLM waves r _±ψ join to form the commonly known (single) divergence wave with eigenvalue being v_r . The items of each discrete right eigenvector, the average state for the diagonal scaling matrix Z , and the average state of the diagonal matrix Λ can be seen in Appendix C. Finally, the Roe entropy-stable numerical flux is (e.g., Barth 1999; Winters & Gassner 2016b; Derigs et al. 2018)

$$\begin{eqnarray*}&&{{\boldsymbol{F}}}_{r}^{\mathrm{ES},\mathrm{Roe}}={{\boldsymbol{F}}}_{r}^{\mathrm{EC}}-\displaystyle \frac{1}{2}{\boldsymbol{R}}\left|{\boldsymbol{\Lambda }}\right|{\boldsymbol{Z}}{{\boldsymbol{R}}}^{T}[\kern-2pt[ { \mathcal V }]\kern-2pt] .\end{eqnarray*}$$

The complete Roe scheme gives greatly increased robustness, but it is more expensive in terms of computation and needs an entropy fix. According to our numerical tests and previous works (Chandrashekar 2013; Derigs et al. 2016; Winters et al. 2017), we find that the hybrid-type entropy stabilization dissipation term has additional benefits in the robustness of the entropy-stable numerical scheme and does not need an entropy fix. The hybrid dissipation term continuously combines the Lax–Friedrichs numerical flux with the Roe scheme. It applies less dissipation in the smooth region and more near the strong shock to ensure robustness and a more accurate resolution of the features (Quirk 1994). The form of the hybrid dissipation term is

$$\begin{eqnarray*}\begin{array}{rcl}\displaystyle \frac{1}{2}{{\boldsymbol{D}}}_{\mathrm{HD}}[\kern-2pt[ {\boldsymbol{U}}]\kern-2pt] & = & \displaystyle \frac{1}{2}{\boldsymbol{R}}\left|{{\boldsymbol{\Lambda }}}_{\mathrm{HD}}\right|{{\boldsymbol{R}}}^{-1}[\kern-2pt[ {\boldsymbol{U}}]\kern-2pt] \simeq \displaystyle \frac{1}{2}{\boldsymbol{R}}\left|{{\boldsymbol{\Lambda }}}_{\mathrm{HD}}\right|{{\boldsymbol{R}}}^{-1}{ \mathcal H }[\kern-2pt[ { \mathcal V }]\kern-2pt] \\ & = & \displaystyle \frac{1}{2}{\boldsymbol{R}}\left|{{\boldsymbol{\Lambda }}}_{\mathrm{HD}}\right|{{\boldsymbol{R}}}^{-1}{\boldsymbol{R}}{\boldsymbol{Z}}{{\boldsymbol{R}}}^{T}[\kern-2pt[ { \mathcal V }]\kern-2pt] =\displaystyle \frac{1}{2}{\boldsymbol{R}}\left|{{\boldsymbol{\Lambda }}}_{\mathrm{HD}}\right|{\boldsymbol{Z}}{{\boldsymbol{R}}}^{T}[\kern-2pt[ { \mathcal V }]\kern-2pt] ,\end{array}\end{eqnarray*}$$

where ∣Λ_HD∣ = (1 − ω)∣Λ_Roe∣ + ω∣Λ_LF∣, with Λ_Roe being the diagonal matrix corresponding to the Roe dissipation and Λ_LF being the diagonal matrix whose elements are ${\lambda }_{\max }^{\mathrm{global}}$ or ${\lambda }_{\max }^{\mathrm{local}}$. The parameter ω is a simple local pressure indicator $\omega ={\left|\tfrac{{p}_{L}-{p}_{R}}{{p}_{L}+{p}_{R}}\right|}^{\tfrac{1}{2}}$. The case of the parameter ω approaching 1 indicates a strong shock, whereas ω → 0 corresponds to a smooth solution. Finally, the hybrid entropy-stable approximate Riemann solver can be written as

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{r}^{\mathrm{ES}}={{\boldsymbol{F}}}_{r}^{\mathrm{EC}}-\displaystyle \frac{1}{2}{\boldsymbol{R}}\left|{{\boldsymbol{\Lambda }}}_{\mathrm{HD}}\right|{\boldsymbol{Z}}{{\boldsymbol{R}}}^{T}[\kern-2pt[ { \mathcal V }]\kern-2pt] .\end{eqnarray} \tag{ 10 }$$

3.2.2. Discretization of the Source Terms

After constructing the specific numerical flux that is in agreement with the discrete version of thermodynamic entropy for system (2), we give the discretization of the source terms in this section. The source terms $\overline{{\boldsymbol{S}}}$ in Equation (6) are the cell-averaged values of the corresponding variables in the cell (i, j, k). However, the nonconservative terms ${\widetilde{{\boldsymbol{S}}}}_{\mathrm{pow}}$ and ${\widetilde{{\boldsymbol{S}}}}_{\psi }$ in Equation (6) should be given special treatment to ensure that the discretizations do not antagonize the entropy stability of the scheme (e.g., Chandrashekar & Klingenberg 2016; Liu et al. 2018; Bohm et al. 2020; Duan & Tang 2020; Wu & Shu 2020). Assuming that Y _pow denotes ${\left(0,{\boldsymbol{B}},{\boldsymbol{v}}\cdot {\boldsymbol{B}},{\boldsymbol{v}},0\right)}^{T}$ and Y _ψ stands for ${\left({\bf{0}},0,{\boldsymbol{v}}\psi ,0,{\boldsymbol{v}}\right)}^{T}$, the discrete nonconservative terms ${\widetilde{{\boldsymbol{S}}}}_{\mathrm{pow}}$ and ${\widetilde{{\boldsymbol{S}}}}_{\psi }$ in the cell (i, j, k) can be approximated following Bohm et al. (2020),

$$\begin{eqnarray*}\begin{array}{rcl}{\widetilde{{\boldsymbol{S}}}}_{\mathrm{pow},i,j,k} & = & -\displaystyle \frac{1}{{V}_{i,j,k}}{\displaystyle \int }_{{V}_{i,j,k}}\left({\rm{\nabla }}\cdot {\boldsymbol{B}}\right){{\boldsymbol{Y}}}_{\mathrm{pow}}{dV}\\ & \approx & -\displaystyle \frac{1}{{V}_{i,j,k}}\left({\overline{{\boldsymbol{Y}}}}_{\mathrm{pow},i,j,k}{\displaystyle \int }_{\partial V}{B}_{n}{dS}+{\displaystyle \int }_{\partial V}\left({\left({{\boldsymbol{Y}}}_{\mathrm{pow}}{B}_{n}\right)}^{\mathrm{ext}}\right.\right.\\ & & \left.\left.-{\left({{\boldsymbol{Y}}}_{\mathrm{pow}}{B}_{n}\right)}^{\mathrm{int}}\right){dS}\right)\\ & \approx & -\displaystyle \frac{1}{{V}_{i,j,k}}\left({\overline{{\boldsymbol{Y}}}}_{\mathrm{pow},i,j,k}\sum _{{\ell }=1}^{6}{B}_{n,{\ell }}{S}_{{\ell }}+\sum _{{\ell }=1}^{6}\left({\left({{\boldsymbol{Y}}}_{\mathrm{pow}}{B}_{n}\right)}_{{\ell }}^{\mathrm{ext}}\right.\right.\\ & & \left.\left.-{\left({{\boldsymbol{Y}}}_{\mathrm{pow}}{B}_{n}\right)}_{{\ell }}^{\mathrm{int}}\right){S}_{{\ell }}\right)\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}{\widetilde{{\boldsymbol{S}}}}_{\psi ,i,j,k} & = & -\displaystyle \frac{1}{{V}_{i,j,k}}{\displaystyle \int }_{{V}_{i,j,k}}\left({\rm{\nabla }}\psi \right){{\boldsymbol{Y}}}_{\psi }{dV}\\ & \approx & -\displaystyle \frac{1}{{V}_{i,j,k}}\left({\overline{{\boldsymbol{Y}}}}_{\psi ,i,j,k}\cdot {\displaystyle \int }_{\partial V}\psi {\boldsymbol{n}}{dS}\right.\\ & & \left.+{\displaystyle \int }_{\partial V}\left({\left({{\boldsymbol{Y}}}_{\psi }\psi \right)}^{\mathrm{ext}}-{\left({{\boldsymbol{Y}}}_{\psi }\psi \right)}^{\mathrm{int}}\right){\boldsymbol{n}}{dS}\right)\\ & \approx & -\displaystyle \frac{1}{{V}_{i,j,k}}\left({\overline{{\boldsymbol{Y}}}}_{\psi ,i,j,k}\cdot \sum _{{\ell }=1}^{6}{\psi }_{{\ell }}{{\boldsymbol{n}}}_{{\ell }}{S}_{{\ell }}\right.\\ & & \left.+\sum _{{\ell }=1}^{6}\left({\left({{\boldsymbol{Y}}}_{\psi }\psi \right)}_{{\ell }}^{\mathrm{ext}}-{\left({{\boldsymbol{Y}}}_{\psi }\psi \right)}_{{\ell }}^{\mathrm{int}}\right){{\boldsymbol{n}}}_{{\ell }}{S}_{{\ell }}\right),\end{array}\end{eqnarray*}$$

where B_n,ℓ = B · n _ℓ is the magnetic field in the outward normal direction n _ℓ evaluated at the interface ℓ. The vector ${\overline{{\boldsymbol{Y}}}}_{\mathrm{pow},i,j,k}$ is the cell-averaged values of the corresponding variables in the cell (i, j, k). $({{\boldsymbol{Y}}}_{\mathrm{pow}}{B}_{n}{)}_{{\ell }}^{\mathrm{ext}}$ and $({{\boldsymbol{Y}}}_{\mathrm{pow}}{B}_{n}{)}_{{\ell }}^{\mathrm{int}}$ are evaluated at the interface ℓ, with the superscripts "ext" and "int" denoting the values obtained from the reconstructions of the neighbors and the current control volume (i, j, k). The similar symbols in ${\widetilde{{\boldsymbol{S}}}}_{\psi ,i,j,k}$ can be explained by way of analogy.

3.2.3. Time Integration

If the right-hand side of Equation (6) is abbreviated as ${{ \mathcal R }}_{{\boldsymbol{U}}}[\overline{{\boldsymbol{U}}}]$, Equation (6) becomes ${\partial }_{t}\overline{{\boldsymbol{U}}}={ \mathcal R }[\overline{{\boldsymbol{U}}}]$. We can advance the semidiscretization form using the following SSP Runge–Kutta method:

$$\begin{eqnarray*}\begin{array}{rcl}{\overline{{\boldsymbol{U}}}}^{* } & = & {\overline{{\boldsymbol{U}}}}^{n}+{\rm{\Delta }}t\cdot { \mathcal R }[{\overline{{\boldsymbol{U}}}}^{n}]\\ {\overline{{\boldsymbol{U}}}}^{* * } & = & \displaystyle \frac{3}{4}{\overline{{\boldsymbol{U}}}}^{n}+\displaystyle \frac{1}{4}\left({\overline{{\boldsymbol{U}}}}^{* }+{\rm{\Delta }}t\cdot { \mathcal R }[{\overline{{\boldsymbol{U}}}}^{* }]\right)\\ {\overline{{\boldsymbol{U}}}}^{n+1} & = & \displaystyle \frac{1}{3}{\overline{{\boldsymbol{U}}}}^{n}+\displaystyle \frac{2}{3}\left({\overline{{\boldsymbol{U}}}}^{* * }+{\rm{\Delta }}t\cdot { \mathcal R }[{\overline{{\boldsymbol{U}}}}^{* * }]\right),\end{array}\end{eqnarray*}$$

which are the convex combinations of the explicit forward Euler integration. The SSP Runge–Kutta scheme is stable under the same time-step size (Gottlieb 2005).

Near the Sun, it is necessary to enhance the spatial resolution desired for physical reasons, which requires a very small time step limited by the CFL condition and much time cost. Significantly low plasma β and high Alfvén wave speed near the Sun can also lead to a time step much smaller than most parts of the computational domain away from the Sun. Furthermore, the plasma β, the Alfvén velocity, the plasma density, and the spatial resolution span several orders of magnitude from the corona to the interplanetary space, which implies a large variation of the CFL stability condition. If one uniform time step is adopted in the whole computational region from the Sun to Earth, the smallest time step constrained by the stability limitation and high resolution near the Sun will be the choice. It will be very time-consuming and seriously decrease the performance of the entropy-stable numerical scheme. Meanwhile, it is a challenge to implement efficient load balancing.

Fortunately, the adoption of the multiple time-stepping integration in the radial direction can ease this situation (van der Ven et al. 1997; Maurits et al. 1998). The method takes different time steps for different partitions and avoids taking a uniform small time step in the whole computational region severely constrained by the unacceptably low CFL number near the Sun. The application of the multi-time-stepping method can enhance the convergence stability, speed up the calculation by mitigating CFL number disparity, facilitate the power of load balancing for the parallel implementation, and save computation resources. For consistency, the implementation of the multiple time-stepping method is described in the following paragraphs. Readers can also refer to Feng et al. (2010) and Jiang et al. (2010).

In the corona, the multiple time-stepping integration is implemented in three radial spherical shell partitions: 1–3.5 R_s , 3.5–10 R_s , and 10–25 R_s (van der Ven et al. 1997; Chang et al. 2005; Feng et al. 2010). Each partition consists of six components and corresponds to one time step. By using the CFL stability condition with the Courant number being 0.6, the minimum of the usual time step ${\rm{\Delta }}{t}_{q}^{\mathrm{CFL}}(q=1,2,3)$ in each of the three spherical shell partitions can be achieved. For q = 1, the bulk time step ${\rm{\Delta }}{t}_{q}={\rm{\Delta }}{t}_{q}^{\mathrm{CFL}}$. For q = 2, 3, the time step Δt_q is further constrained as follows: Δt₂ = M₁Δt₁ with ${M}_{1}=\mathrm{int}({\rm{\Delta }}{t}_{2}^{\mathrm{CFL}}/{\rm{\Delta }}{t}_{1})$, and then Δt₃ = M₂Δt₂ with ${M}_{2}=\mathrm{int}({\rm{\Delta }}{t}_{3}^{\mathrm{CFL}}/{\rm{\Delta }}{t}_{2})$, as shown in Figure 4 of Feng et al. (2010). For example, the time steps Δt₁, Δt₂, and Δt₃ have the relation 4Δt₁ = 2Δt₂ = Δt₃. Correspondingly, the flow in the first spherical shell partition travels four iterations of Δt₁, and the flow in the second partition advances two steps of Δt₂ to fit in the synchronization time step Δt₃, as shown in Figure 4 of Feng et al. (2010). The bulk time steps for the three spherical shell partitions will be determined again after each synchronization time level. For the overlapping boundary interface between two adjacent spherical shell partitions, the former partition needs to advance M₁ (M₂) steps of Δt₁ (Δt₂) to match one time step of Δt₂ (Δt₃) for the latter partition. Hence, the guard cell filling can depend on a simple time interpolation of solution, that is, the boundary values of the former partition are obtained by interpolating the associated interior values of the latter partition and vice versa. Similarly, the multiple time-stepping integration is also applied to the spherical shell II ranging from 25 R_s to 230 R_s .

In the solar wind simulation, the pressure may become negative owing to numerical errors in the near-Sun, especially in the strong magnetic field and low plasma β region. Recent years have witnessed some progresses in designing positivity-preserving schemes. For example, the decomposition of magnetic field into the potential magnetic field B ₀ and its perturbed component B ₁, such as B = B ₀ + B ₁, may mitigate such an occurrence of negative pressure to some extent. Alternatively, some schemes replace the conservative energy equation with the entropy equation or the internal energy density equation or pressure equation optionally to avoid the occurrence of negative pressure in the regions with low kinematic pressure, low plasma β, small internal energy, and smooth regions without shocks (e.g., Ryu et al. 1993; Bryan et al. 1995; Balsara & Spicer 1999; Cargill et al. 2000; Li et al. 2008; Feng et al. 2010, 2014b; Tóth et al. 2012; Usmanov et al. 2012). These positivity-preserving operations are known as "dual-energy" methods. However, such design is sensitive to the switching conditions. For some comments about this, the reader may refer to Section 2.7 of Usmanov et al. (2018). Besides, some researchers try to avoid the occurrence of negative pressure by establishing hybridized flux for hyperbolic conservation laws. Specifically, by seeking a convex combination of the high-order flux and the first-order Lax–Friedrichs flux, the parameterized maximum-principle-preserving (MPP) flux-corrected limiting procedure can achieve a provably positivity-preserving high-order scheme for both hydrodynamic and MHD flows (e.g., Zhang & Shu 2010b, 2010a, 2011; Xiong et al. 2013, 2016; Christlieb et al. 2015a, 2015b, 2016; Feng 2020b). Also, entropy-stable Riemann solvers have been designed so as to bear the property of positivity-preserving in the finite-volume, finite-difference, and DG method for the hyperbolic conservative system (e.g., Derigs et al. 2016, 2017; Winters & Gassner 2016b; Winters et al. 2017; Liu et al. 2018; Ranocha 2018). To proceed with entropy-stable numerical flux, we propose the EC-GLM-MHD model in the current paper. The designed model aims to handle the numerical challenge of producing negative pressure or negative density in solar wind simulation without the necessity of special treatments described above, such as magnetic field decomposition, dual energy method, or the parameterized MPP flux-corrected limiting procedure. Due to the benefit of the positivity-preserving entropy-stable numerical fluxes, the proposed EC-GLM-MHD solar wind model can further reduce the risk of negative pressure or density by restricting the reconstructed density and pressure (e.g., Waagan 2009; Fuchs et al. 2010, 2011; Balsara 2012; Feng et al. 2017).

Additionally, the entropy pressure equation can be further utilized to nip the negative pressure in the bud. We describe the implementation in the following. Considering the semidiscrete form ${\partial }_{t}\overline{{\boldsymbol{U}}}={ \mathcal R }[\overline{{\boldsymbol{U}}}]$ of the ideal EC-GLM-MHD equations, we can contract it with the entropy vector ${ \mathcal V }$ (Derigs et al. 2016):

$$\begin{eqnarray*}&&{{ \mathcal V }}^{T}\partial {\overline{{\boldsymbol{U}}}}_{t}={{ \mathcal V }}^{T}{ \mathcal R }[\overline{{\boldsymbol{U}}}].\end{eqnarray*}$$

According to the chain rule and the definition of the entropy vector, we obtain ${{ \mathcal V }}^{T}\partial {\overline{{\boldsymbol{U}}}}_{t}={{ \mathcal S }}_{t}$. Then, we deduce the expression of the time evolution of physical entropy density:

$$\begin{eqnarray}&&{{ \mathcal S }}_{t}={{ \mathcal V }}^{T}{ \mathcal R }[\overline{{\boldsymbol{U}}}].\end{eqnarray} \tag{ 11 }$$

The entropy-conservative flux provides us the entropy equality:

$$\begin{eqnarray}&&\displaystyle \frac{\partial { \mathcal S }}{\partial t}+{\rm{\nabla }}\cdot \left({\boldsymbol{v}}{ \mathcal S }\right)=0.\end{eqnarray} \tag{ 12 }$$

From Equations (11) and (12), we obtain a discrete method for the spatial derivative of the entropy flux:

$$\begin{eqnarray*}&&-{\rm{\nabla }}\cdot \left({\boldsymbol{v}}{ \mathcal S }\right)={{ \mathcal V }}^{T}{ \mathcal R }[\overline{{\boldsymbol{U}}}].\end{eqnarray*}$$

Therefore, we can see that Equation (11) provides a consistent, discrete evolution of the entropy. The discrete form of the physical entropy density equation can be update together with the MHD Equations (6) and evolves in time with other fluid variables. From the updated physical entropy density ${ \mathcal S }$, we derive a new expression for the pressure p_s , which guarantees the pressure to be positive:

$$\begin{eqnarray}&&{p}_{s}=\exp \left(\displaystyle \frac{\gamma -1}{\rho }{ \mathcal S }+\gamma \mathrm{ln}\rho \right).\end{eqnarray} \tag{ 13 }$$

The entropy-stable numerical scheme prescribed in Section 3 guarantees that the entropy satisfies the discrete version of the second law of thermodynamics everywhere (Derigs et al. 2016). Thus, the proposed positive pressure preserving method can be applied to any regions of the flow. In the calculation, we update the entropy density ${ \mathcal S }$ through the discrete form (11), together with the conservative values U , and then we obtain ${{ \mathcal S }}^{n+1}$ and U ⁿ⁺¹. If the updated energies $\tfrac{{E}_{\mathrm{int}}^{n+1}}{{E}^{n+1}}\leqslant 0.01$, we use Equation (13) to get the pressure ${p}_{s}^{n+1}$ and recompute the updated total energies by ${E}^{n+1}=\tfrac{{\left(\rho {\boldsymbol{v}}\right)}^{2}}{2{\rho }^{n+1}}+\tfrac{{p}_{s}^{n+1}}{\gamma -1}+\tfrac{{\left({{\boldsymbol{B}}}^{n+1}\right)}^{2}}{2}+\tfrac{{\left({\psi }^{n+1}\right)}^{2}}{2}$ to make the scheme consistent.

As we know, the inner boundary of the corona is situated at the surface of the Sun (1 R_s ), which is the sub-Alfvénic or subsonic region. Thus, proper treatment is a necessity to avoid the numerical inconsistency. The projected normal characteristic method is a good choice (Nakagawa et al. 1987). It can handle both the prescribed physical boundary constraints and the underlying MHD equations. The characteristic compatibility equations correspond to the outgoing waves. These numerical boundaries can give some variables that are not defined by the physical conditions. As usual, we adopt the first-order upwind differencing to discrete the characteristic compatibility equations and solve these values entirely determined by the solution at and within the boundary, because only the grid cells on the interior side of the boundary are available. Through tests and theoretical methods, researchers have verified that the global order of the scheme is not affected if the order at or near the boundary is one order lower than that away from the boundary (Gustafsson 1981; Thompson 1990; Poinsot & Lele 1992; Colonius 2004). They also stated that an unstable solution may appear if the order of the boundary approximation is the same order (second or higher) of the interior scheme. Additionally, we adopt the invariant entropy and the limit of the plasma mass flux escaping from the solar surface to obtain two independent constraints, as previous works did (e.g., Yang et al. 2012; Feng et al. 2015, 2017; Li et al. 2018). The invariant entropy at the inner boundary satisfies the stability of entropy, which is consistent with the EC-GLM-MHD model used in the computational region. If the constraints are not enough to determine all the MHD variables, the nonreflecting boundary condition can be used, meaning no backflow being allowed along the characteristic curve and the plasma permitted to leave the Sun freely. Through the simulation of the solar wind for CR 2183, we did not find any significant evidence of instabilities through the whole run. The current boundary condition treatments do not affect the proposed EC-GLM-MHD model in the interior of the computational domain.

In this section, we present and analyze the numerical results of two simple 1D benchmark tests and the large-scale ambient solar wind for CR 2183 by utilizing the entropy-stable scheme introduced in the previous sections.

6.1. Simple Benchmark Tests

To illustrate the entropy consistency and the increased robustness of the entropy-stable scheme for the ideal EC-GLM-MHD equations, we consider two simple 1D benchmark tests, following Derigs et al. (2017, 2018). These two 1D tests are the modified Brio and Wu shock tube and strong shock traveling through uniform fast-moving and magnetized medium.

6.1.1. Modified Brio and Wu Shock Tube

In this part, the MHD shock tube originally proposed by Brio & Wu (1988) is chosen to assess the entropy consistency of the entropy-stable scheme. Similar to Derigs et al. (2018), a constant velocity is added to the left side to create a rarefaction wave. The detailed initial and boundary conditions are given in Table 1. Figure 1 displays the plot of the fluid density at t = 0.1 achieved from the traditional Roe scheme (red line), the entropy-stable Roe scheme (blue line), and the hybrid entropy-stable scheme (green line). It is clear that the rarefaction wave on the left-hand side is much smaller, which is the result of the compound wave in the MHD shock tube. From the enlarged plot, there is an apparent unphysical oscillation for the traditional Roe scheme at around x = 0.495, while the densities obtained from the two entropy-stable schemes do not have similar oscillation at the same position. This means that the solution obtained from the Roe scheme violates the entropy condition. By comparison, there is no inconsistency visible for the result achieved from the entropy-stable schemes.

Zoom In Zoom Out Reset image size

Table 1. Initial Conditions: Modified Brio and Wu Shock Tube Test

Variables	x < 0.5	x ≥ 0.5
ρ	1.0	0.125
p	1.0	0.1
vx	1.75	0.0
vy , vz	0.0	0.0
Bx	0.75	0.75
By	1.0	−1.0
Bz	0.0	0.0

Domain size	${x}_{\min }=0,\ {x}_{\max }=1$
Adiabatic index	γ = 2.0
Courant number	CFL = 0.5
Simulation end time	${t}_{\max }=0.1$
Boundary conditions	zero gradient

Download table as:  ASCIITypeset image

Watching the details of blue and green lines at around x = 0.495 carefully, we find that the hybrid entropy-stable scheme is a bit more dissipative than the entropy-stable Roe scheme. A similar feature can be detected at around x = 0.670. In brief, the results of the traditional Roe and entropy-stable schemes are basically the same except for the entropy glitch of the traditional Roe scheme in the vicinity of x = 0.495.

6.1.2. A Strong Shock Traveling through a Uniform Fast-moving and Magnetized Medium

To test the robustness of the entropy-stable scheme, we simulate a strong shock traveling through a uniform fast-moving and magnetized medium, which is typical in large-scale astrophysical modeling. We list the initial and boundary conditions in Table 2. From Table 2, we can see that the relative magnitude of the initial pressure pulse is 10⁶ and the initial pressure at ∣x∣ > 0.1 tends to zero. It creates a near-infinite shock strength. Meanwhile, this problem is prone to produce negative pressure during a simulation.

Table 2. Initial Conditions: A Strong Shock Traveling through a Uniform Fast-moving and Magnetized Medium

Variables	∣x∣ ≤ 0.1	∣x∣ > 0.1
ρ	1.0	1.0
p	1.0	1.0 × 10−6
vx	10.0	10.0
vy , vz	0.0	0.0
Bx	1.0 × 10−2	1.0 × 10−2
By , Bz	0.0	0.0

Domain size	[−1, 1]
Adiabatic index	γ = 5/3
Courant number	CFL = 0.8
Simulation end time	${t}_{\max }=5.0\times {10}^{-2}$
Boundary conditions	periodic

Download table as:  ASCIITypeset image

In Figure 2, we present the plots of thermal pressure after the first time step and the results of density at the simulation end time for the fast shock traveling through a magnetized medium. From the figure, we can see that the traditional LLF scheme breaks down during the first time step. This is because the traditional LLF finite-volume method produces negative pressure in the first step as illustrated in the top left panel of Figure 2. However, the entropy-stable LLF and hybrid entropy-stable schemes do not produce negative pressure after the first time step and can run to the simulation end time robustly. Due to just one time step, the pressures achieved from the two entropy-stable schemes overlap exactly, which can be seen from the top right panel of Figure 2. From the bottom right panel of Figure 2, we can see that two shock fronts appear as expected. Note that the amplitude of density obtained from the hybrid entropy-stable scheme is a little larger than that from the entropy-stable LLF scheme. That is because the entropy-stable LLF scheme introduces more dissipation than the hybrid entropy-stable scheme.

Zoom In Zoom Out Reset image size

6.2. Ambient Solar Wind

In this section, we utilize the hybrid entropy-stable EC-GLM-MHD solar wind model to simulate the large-scale structures of the solar wind for CR 2183, with the initial solar surface temperature and number density being 1.3 × 10⁶ K and 1.5 × 10⁸ cm⁻³, respectively, and compare the results with remote sensing observations and in situ data. This solar wind simulation was run on the TianHe-1 (A) supercomputer located at the National Supercomputer Center in Tianjin, China. Using 16 computing nodes with 192 CPU cores, the code for the corona obtained a steady-state equilibrium in about 32.7 hr of wall-clock time, with the minimum time step being 9.22 × 10⁻⁵ and the maximum time step being 5.81 × 10⁻³, while the run for the interplanetary region took about 11.2 hr of wall-clock time to obtain a steady-state solution. Indeed, it is a bit of a waste of time simulating ambient solar wind with a low plasma β region, especially near the Sun. In order to improve the performance of the EC-GLM-MHD model and avoid the Courant limit further, we will adopt an implicit FVM solved by the LU-SGS technique together with the multiple time-stepping integration in our next work.

To test the robustness and positivity-preserving properties of the proposed entropy-stable numerical scheme in dealing with the low plasma β and strong magnetic field problem, the initial potential fields for the input of the entropy-stable EC-GLM-MHD model are obtained from the SDO/Helioseismic and Magnetic Imager (HMI) synoptic map by retaining spherical harmonics of degree 50. In Figure 3, we show the radial magnetic field B_r and the plasma β over the full solar surface, with the left panel denoting the magnetic field obtained from the SDO/HMI full-disk magnetogram, the middle panel standing for the initial radial magnetic field obtained from the potential field model, and the right panel being the initial distributions of plasma β. The color bars of the magnetic field images are saturated between 70 G (red positive) and –50 G (blue negative). The plasma β is as low as a magnitude of 10⁻⁴, which is a challenge for the solar wind simulation. However, the present entropy-stable EC-GLM-MHD model can run robustly and yield satisfactory results without incurring negative pressure or density. Through some simple calculations, we obtain the total unsigned magnetic flux, ${\rm{\Phi }}({R}_{s})={\int }_{0}^{2\pi }{\int }_{0}^{\pi }| {B}_{r}({R}_{s},\theta ,\phi )| {R}_{s}^{2}\sin \theta d\theta d\phi$, for the SDO/HMI magnetogram and the initial magnetic field synoptic map. They are 3.1 × 10²³ Mx and 1.9 × 10²³ Mx, respectively. The difference may be due to the spherical harmonics of the initial potential field model and the lower resolution of the MHD model (Yeates et al. 2018). Accordingly, we evaluate the open magnetic flux at 2.5 R_s when the steady-state solution is reached. It is 4.0 × 10²² Mx. Due to the influence of non-force-free effects in the corona, the drop is not very heavy, which is similar to the CESE-MHD model (Feng et al. 2012c). Through the formula $E={\int }_{{R}_{s}}^{2.5{R}_{s}}{\int }_{0}^{2\pi }{\int }_{0}^{\pi }\tfrac{| {\boldsymbol{B}}(r,\theta ,\phi ){| }^{2}}{8\pi }{r}^{2}\sin \theta d\theta d\phi {dr}$, we also calculate the total magnetic energy E_m of the MHD model from r = R_s to r = 2.5R_s and the potential magnetic energy E_p from the potential field source surface model with the source surface r = 2.5R_s . They are E_m = 4.06 × 10³² erg and E_p = 3.50 × 10³² erg. The excess of nonpotential energy E_m above the potential magnetic energy E_p is the relative free energy. The ratio of the total to the potential magnetic energy is approximately 1.16, which is lower than the results (around 1.4) of the Optimization NLFFF (Tadesse et al. 2014), Linear MHS (Bogdan & Low 1986), Evolving magnetofrictional (Yeates 2014), and MHD-MAS (Mikić et al. 1999) models during the solar eclipse on 2015 March 20 as shown by Yeates et al. (2018). That is because the initial input of the current model is obtained from the potential field model based on the line-of-sight photospheric magnetogram from HMI, which cannot reflect the real information of active regions, although there are several during this declining phase. Also, the present model does not include any mechanism to energize the magnetic field near the Sun. This may be another reason.

Zoom In Zoom Out Reset image size

In the solar declining and minimum phases, the high-latitude open-field regions are usually the source of the fast solar wind, and the low-latitude closed-field regions are closely associated with the slow wind, so it is important to consider the distributions of the open- and closed-field regions. To this end, we draw the distributions of the open- (black) and closed-field (white) regions obtained from the potential field source surface (PFSS) model and the corresponding SDO/AIA 193 Å observations in Figure 4. The red lines overlaid on the right panel stand for the modeled boundaries between the open- and closed-field regions. Despite the assumptions and simplifications of the PFSS model, it is widely used to predict the coronal magnetic field's topology qualitatively, and it has become a rough standard to which our MHD model can be compared. From the left panel in Figure 4, we can see that the southern and northern polar coronal holes occupy the high latitudes of the southern and northern hemispheres. There are a series of equatorward extensions in the southern hemisphere around Carrington longitude 70°–190° and in the northern hemisphere from 230° to 360°, respectively. These robust features can also be seen in the modeled results and the observations as shown in the right panel of Figure 4. However, the observed extending coronal holes extend much farther than the modeled ones in both the southern and northern hemispheres. The observed polar coronal hole in the south pole is much smaller than that obtained from the PFSS model and our MHD model, which can be discerned around Carrington longitude 180°–360°. Given the boundary conditions and the sensitivity of the open-/closed- footpoint regions, the modeled results are in reasonable agreement with the corresponding observations.

Zoom In Zoom Out Reset image size

Since the magnetic field lines are open beyond the source surface (usually 2.5 R_s ) in the PFSS model, we present contours of the radial magnetic field at 2.5 R_s obtained from the entropy-stable EC-GLM-MHD model and the PFSS model in Figure 5. From this figure, we can see that the contour lines of the radial magnetic field achieved from both the MHD model and the PFSS model are basically similar to each other in the sizes, shapes, and locations. However, the contour lines in the low latitudes obtained from the MHD model seem to deflect a little more abruptly than the corresponding ones from the PFSS model. Furthermore, the strengths of the radial magnetic field in the polar regions achieved from the MHD model are slightly weaker than those obtained from the PFSS model. Figure 6 displays the synoptic maps of the east and west limbs at 2.5 R_s recorded from the LASCO C2/SOHO (top panels) and the modeled plasma number density N and radial speed v_r at the same height (bottom panels). The red lines in the top panels and the white ones in the bottom panels are the modeled magnetic neutral lines (MNLs, the bases of the heliospheric current sheets (HCSs) at 2.5 R_s ). From the top images, we can see that the MNLs are wrapped in the white-light polarized brightness (pB) structures, and the regions are getting darker from low latitude to high latitude. For the bottom panels, the MNLs are closely related to the high plasma density and low radial velocity, while both polar regions are characterized by low plasma density and high radial speed. These close associations among the low-speed wind, the high-density flow, and the MNL can also be revealed by combining multistation interplanetary scintillation (IPS) observations with coronagraph and magnetograph observations (Rickett & Coles 1991). At this point, the comparisons of the numerical results and the corresponding observations in Figure 6 also show that the observed white-light bright structures are surrounded by high-density, low-velocity flows, and the dark regions correspond to low-density, high-velocity flows. Moreover, there is a wide flat bulge structure from 40° to 210° in each of the top panels. The modeled results divide the flat hump into two, with the dividing line being ϕ = 160° approximately, as shown by the MNLs in Figure 6. Broadly speaking, the pB synoptic maps from LASCO C2/SOHO provide us a criterion to verify the numerical results to some extent.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For further comparison, we plot the coronal magnetic field configurations and the white-light pB image obtained from the entropy-stable EC-GLM-MHD model with the central meridian plane located at longitude 180°, and the corresponding white-light observations on 2016 October 2 from the LASCO/SOHO in the top panels of Figure 7. The bottom panels of Figure 7 are two zoom-ins of the image in the top left, that is, the enlarged coronal magnetic field configurations of the pseudo-streamers. In the synthesized and observed pB images (middle and right images of the top panels), two broad bright structures stretch outward at the northeast and southwest limbs. Meanwhile, two faint bright structures extend radially at the east and west limbs for the observations and the simulated results. From the left image of the top panels in Figure 7, we find that the inner corona is composed of complex substreamers by a series of closed-loop systems at middle and low latitudes, and the outer corona is mainly composed of the open-field region at high latitudes and two helmet streamers at middle latitudes, around 45°. As shown by the bottom panels, we can see that obvious pseudo-streamers are captured by our model, with the magnetic field lines of the same polarity being separated (Wang et al. 2007). Unlike the helmet streamers giving the origin of the heliospheric current sheet (HCS), these specific configurations of pseudo-streamers are simply plasma sheets. At their bases, pseudo-streamers are the closed loops and usually more confined near the surface of the Sun than the helmet streamers. Owing to their peculiar magnetic configuration, pseudo-streamers do not emit plasma blobs as helmet streamers do (Sheeley et al. 1997). Combining the topological structures of the magnetic field with the pB images, we can see that the helmet streamers correspond to the broad bright structures at the northeast and southwest limbs and end in the magnetic neutral face that extends outward, while the pseudo-streamers result in less bright structures at the east and west limbs with no current sheet (Pinto et al. 2011). Additionally, the structures of the streamers are relatively simple, which is in accordance with the typical features of the solar minimum and declining phases.

Zoom In Zoom Out Reset image size

Figure 8 presents the modeled results for the radial velocities v_r and plasma number density N on the meridional planes of ϕ = 90°–270° and ϕ = 180°–0° from 1 R_s to 20 R_s . Seen from the magnetic field lines, the figure is featured by the large-scale magnetic field, forming a belt around the Sun. Meanwhile, only large helmet streamers can extend to nearly 5 R_s , where the HCS starts. The open magnetic field has negative polarity in the southern polar coronal hole and positive polarity in the northern polar coronal hole. The coronal magnetic topology determines the formation of the solar wind and has a direct influence on the distributions of the fast and slow solar winds, as well as the plasma density, via the sizes, positions, and widths. As a result of the acceleration of the wind flow along its path starting from the Sun's surface, the solar wind flow manifests as latitudinal inhomogeneities as shown by the top panels of Figure 8. The separations between the fast and slow solar winds are easily visible in this figure. These images also demonstrate the fact that the high-density low-speed flows are mostly close to the helmet streamers, while the low-density high-speed flows are closely associated with the open-field areas. These quiescent streamers with the magnetic field organized by low-order multipoles are prevalent in the declining and minimum phase (e.g., Antonucci et al. 2000). Figure 9 shows the modeled plasma temperature T, plasma number density N, and radial velocity v_r at 20 R_s and 215 R_s , respectively. Comparing with the bottom panels in Figure 6, we can see that the number densities decrease as the heliocentric distance increases, while the radial speeds increase at the same time. The solar wind speed at high latitudes is considerably faster than that at the solar equator, which is consistent with the multistation IPS measurements during the declining and minimum phases (Tokumaru et al. 2010). Moreover, the same solar wind structures at 215 R_s shift leftward about 50° compared with those at 20 R_s , and the same structures at 20 R_s are to the left slightly compared with the ones at 2.5 R_s .

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 10 displays the entropy-stable EC-GLM-MHD results for the magnetic field and radial speed on the equatorial plane and the neutral surface of magnetic field with the central meridian located at longitude 90° from 21.5 R_s to 215 R_s . The black arrowheads in the left panel represent the direction of the magnetic field, and the colored contours denote the radial velocity of the solar wind. There are three distinct compression regions in the left panel, where the radial velocity is lower. Meanwhile, the left panel demonstrates three magnetic sectors where the fast solar wind (>500 km s⁻¹) emanates. For the right panel, the neutral surface of magnetic field is almost in the equatorial plane. The left side of the neutral surface stretches northward slightly, and the right side has a little southern excursion. During the declining phases, this southward shift of the HCS is mainly caused by the dominance of the solar activity in the southern hemisphere (Wang 2014), which is pretty obvious from our numerical results in Figure 5. The feature that the open coronal magnetic field is asymmetric in the declining phase is consistent with the former descriptions (e.g., Virtanen & Mursula 2016; Getachew et al. 2017; Koskela et al. 2018; Antonucci et al. 2020).

Zoom In Zoom Out Reset image size

Although the present simulation just produces a quasi-steady solar wind solution and there is no direct time-varying information, we can also roughly create a time series of the solar wind parameters at fixed radial distance of 1 au over a CR (27.27 days) by selecting only the data that are at appropriate radius, latitude, and longitude from the MHD model solution. That is because the solar longitude values correspond to a time span of one CR and can be converted to a time axis with solar longitude 0° corresponding to the end of day 27 of a CR and solar longitude 360° being the start date of the CR.

In Figure 11, we exhibit the comparison of the results obtained from the current entropy-stable EC-GLM-MHD solar wind model and the measurements recorded from the OMNI archive. The blue lines are the modeled results, the gray lines stand for the OMNI in situ observations, and the red ones represent the smoothed observations. By and large, these images reveal that the MHD model has the ability to obtain the analogous temporal variation trends for the plasma parameters at Earth's orbit. For the radial speed, the MHD model produces three peaks centering at days 4, 8, and 24 and three troughs centering at days 2, 6, and 19. By comparing the modeled radial speed and the corresponding in situ observations, we can see that the MHD model captures the peaks and troughs with reasonable accuracy, although the simulated peak centered at day 4 is about 1.5 days earlier than the measurements, the modeled peak at day 8 is approximately 80 km s⁻¹ lower than the observations, and the obtained trough centered at day 19 is lower and narrower than the observed results. Moreover, the entropy-stable EC-GLM-MHD model produces slightly higher plasma number density, a little higher temperature, and a bit weaker radial magnetic field than the observed ones, which are consistent with some previous works (e.g., Lepri et al. 2008; Stevens et al. 2012; Linker et al. 2017; MacNeice et al. 2018). The discrepancies partially come from the simplified input of using the Parker 1D solar wind model, the single fluid model used, the occurrence of coronal mass ejections (CMEs) during CR 2183 (a CME was recorded by LASCO C2/SOHO at 0424 UT on 2016 November 5 according to the "Near-Earth Interplanetary Coronal Mass Ejections Since January 1996" list ³ , problems with the photospheric magnetic field, the missing of more sound physical-based and observation-supported coronal heating and solar wind acceleration mechanism, and other reasons as pointed out by Feng et al. (2010). Some specific explanations are described in the following paragraphs.

Zoom In Zoom Out Reset image size

The current simulation models the solar wind in the quiet declining phase without considering the disturbances triggered by transient solar events. However, the effects of the CME or other disturbances are embedded in the observations and cannot easily be separated from the ambient solar wind. This may influence the comparison of the modeled and observed results. Moreover, this quasi-steady model is designed for the large-scale solar wind simulation, which assumes that the photospheric field evolves slowly and the small-scale fields are trivial for determining the large-scale solar wind structures of corona and heliosphere. It only basically gives the trends of the solar wind plasma parameters over the course of a CR.

Although the present coronal heating and solar wind acceleration mechanism can give favorable background solar wind characteristics, we still have reasons to believe that simple volumetric heating formulae cannot be the only heating and acceleration mechanism, because it is not connected directly to any specific physical process and ignores the energy contained in subgrid-scale turbulence. Therefore, more complete descriptions of the important physical process and some other currently unknown sources should be added to the MHD model. By adding a set of additional equations to the basic MHD system, some global MHD model developers integrated some specific physical processes, such as the energy content, transport, and dissipation of the Alfvén wave fluctuations (e.g., Chandran et al. 2011), into their MHD models (e.g., Usmanov et al. 2014, 2018; van der Holst et al. 2014; Downs et al. 2016; Wiengarten et al. 2016). It is anticipated that coronal heating and solar wind acceleration may be understood in an even better fashion by future theoretical study and comprehensive analysis of more solar observations.

As for the magnetograms, there are a number of well-documented problems. For instance, magnetograph measurements assume that the atmosphere is uniform over a given detector pixel; magnetogram calibration may introduce significant differences because it depends on polarimetric models to account for temporal, spatial, and spectral resolution deficiencies (Leka & Barnes 2012); polar fields are usually inferred by using observations over a solar rotation and interpolation of fields from a fairly limited lower-latitude region because their measurements are poor (e.g., Wang & Sheeley 1995; Riley et al. 2002; Li et al. 2021). Given these issues with the magnetograms, the underestimation of the magnetic field in the heliosphere is hardly surprising (e.g., MacNeice et al. 2018). To remedy the imperfections of the photospheric magnetic field, one possible way is to use magnetograms from different observatories to ensure that the magnetic field can be reliably intercalibrated. Through quantifying the effects of the polar field strengths relative to the equatorial values, the errors in the synoptic magnetic field measurements, the polar field extrapolations, and currently saturation corrections, the most critical points of the photospheric magnetic field can be improved. Additionally, the polar magnetic field is expected to be measured by spacecraft away from the ecliptic plane. The Polarimetric and Helioseismic Imager on board the Solar Orbiter (Müller et al. 2013) and the Full-disk Vector Magnetograph on board the Advanced Space-based Solar Observatory (Gan et al. 2019) will give a more accurate measurement for polar magnetic field. It is expected that a more realistic description of ambient solar wind will be achieved by using better data of photospheric magnetograms as input for MHD models.

From the numerical results and the possible reasons analyzed above, we can see that modeling consistently such an interplay of length scales and forecasting the ambient solar wind near Earth are challenging tasks, even if the solar activity is low, the large-scale magnetic field configuration evolves less rapidly, and the CMEs are less frequent (Owens & Forsyth 2013). Continued efforts are required to improve the capabilities of forecasting more realistic ambient solar wind for the MHD model (Schrijver et al. 2015). In order to further judge whether the differences between the numerical results and the observations are related to using the entropy-stable scheme proposed here, we perform a simulation with the positivity-preserving HLL solver proposed by Wu & Shu (2019) and overlay its results with the green lines in Figure 11. The results that we achieved from the comparison of the two different numerical schemes for solar wind simulation to in situ observations at 1 au are basically consistent with previous researches (Lee et al. 2009; Jian et al. 2011; Gressl et al. 2014). According to the data presented in Figure 11, we quantify the Pearson correlation coefficient of the simulated and observed radial velocity, plasma number density, temperature, and magnetic field magnitude. They are 0.68, 0.48, 0.49, and 0.65 for the hybrid entropy-stable scheme. Meanwhile, for the HLL scheme, they are 0.67, 0.49, 0.52, and 0.62, respectively. The two schemes produce almost the same simulation results for the solar wind parameters with correlation coefficients between 0.48 and 0.68. There is no trend as to which numerical scheme gives systematically more realistic solar wind plasma parameters near Earth than the other. Thus, the differences between the observations and the modeling results are largely the result of the reasons listed above, rather than the entropy-stable numerical flux.

To validate the entropy stability of the hybrid entropy-stable EC-GLM-MHD solar wind model, we plot the evolution of the normalized total entropy

$$\begin{eqnarray*}&&\begin{array}{l}\quad \left({\sum }_{i,j,k}{\int }_{{V}_{i,j,k}}{\rm{\Delta }}{ \mathcal S }{dV}={\sum }_{i,j,k}{\int }_{{V}_{i,j,k}}{ \mathcal S }(i,j,k,t)\right)\\ \qquad \qquad \qquad \qquad \times {dV}-{\sum }_{i,j,k}{\int }_{{V}_{i,j,k}}{ \mathcal S }(i,j,k,0){dV}\end{array}\end{eqnarray*}$$

until a final time, which in principle should numerically decay with respect to time owing to the numerical dissipation if the scheme is entropy stable (e.g., Chandrashekar & Klingenberg 2016; Bohm et al. 2017; Biswas & Dubey 2018; Liu et al. 2018; Friedrich et al. 2019; Gouasmi et al. 2019, 2020; Bhoriya & Kumar 2020; Duan & Tang 2020, 2021; Wu & Shu 2020). The evolution of the differences between the current and the initial total entropy is displayed in the left panel of Figure 12. In Figure 12(a), we clearly observe that the quantity does not increase in time as expected, which confirms that the fully discretized scheme for solar wind simulation is entropy-stable with the Venkatakrishnan limiter applied in the reconstruction of state variables. The discrete total entropy declines rapidly for the first 40 hr and remains almost constant until the end of the simulation. That is because more dissipation should be added to the entropy-stable scheme to make the scheme robust for the earlier simulation. As time goes by, the oscillations are damped by the dissipation of the entropy-stable scheme and the numerical solution becomes smoother. Since the magnetic field in the coronal region is strong, the temporal evolution of the relative magnetic divergence error, $\mathrm{div}B={\sum }_{k=1}^{M}\tfrac{\left|{\int }_{{V}_{k}}{\rm{\nabla }}\cdot {BdV}\right|}{{\int }_{{S}_{k}}\left|{\boldsymbol{B}}\right|{ds}}/M$, with M being the total number of cells, is displayed in this region. From the evolution of the relative magnetic divergence error as shown in Figure 12(b), we can see that the biggest increasing rate occurs around the heading period and then remains broadly constant at around 10⁻⁵ for a long running time. During the whole simulation, there is no obvious large divergence error accumulation. By and large, the built-in GLM divergence cleaning mechanism of the hybrid entropy-stable solar wind model can be comparable to our previous works and some commonly used models (e.g., Duffell & MacFadyen 2011; Gaburov et al. 2012; Pakmor & Springel 2013; Zhang & Feng 2016; Li et al. 2018, 2021).

Zoom In Zoom Out Reset image size

To assess the advantages and disadvantages of the EC-GLM-MHD solar wind model, we perform the same solar wind simulation with the HLL scheme. Without the operations (such as magnetic field decomposition, the dual energy method, and the parameterized MPP flux-corrected limiting procedure) to preserve the positivity of pressure and density, the HLL scheme produces some negative values for pressure and density during the solar wind simulation, with the initial potential magnetic field obtained by retaining spherical harmonics of degree 50. In order to obtain the physical solar wind parameters by using the HLL scheme, we adopt the positivity-preserving HLL solver (Wu & Shu 2019) and add several positivity-preserving methods (magnetic field decomposition, the parameterized MPP flux-corrected limiting procedure, and the non-negativity-preserving reconstruction) to our code. During the simulation, we record the discrete total entropy and then draw the time evolution plot of the normalized total entropy obtained from the HLL scheme to compare with the hybrid entropy-stable scheme. In Figure 13, we can see that the total entropy obtained from the HLL scheme has a temporary gain at around time = 10 hr. Because the entropy glitch does not last long and is quite small, it has little effect on the simulation results for the ambient solar wind. We also plot some general structures of the background solar wind with the data obtained from the HLL scheme. From the comparison, we find that the entropy-stable scheme and the HLL scheme generate almost the same structures morphologically. Through simple computations, the relative errors of all physical variables between the two schemes are around 2% (Feng et al. 2017). For the multiscale nature of the solar wind simulation, assessing the performance of two different solar wind models quantitatively is a challenging task in astrophysical plasma modeling (e.g., Wu & Dryer 2015). The performance of operational frameworks for forecasting the ambient solar wind has been assessed by many studies (e.g., Arge et al. 2003; Owens et al. 2008; Lee et al. 2009; Jian et al. 2011, 2015; Devos et al. 2014; Gressl et al. 2014; Reiss et al. 2016, 2019; Cohen 2017; MacNeice et al. 2018; Kumar et al. 2020). Since each model has its advantages and disadvantages, there is not a winner that performs best. By comparing the solution of the solar wind model at the fixed location of a spacecraft, they reported that almost all models have an acceptable deviation comparable to that obtained from in situ observations for large-scale solar wind structures and for timescales of one to several days. In our work, the time-varying plots of the solar wind speed, temperature, density, and magnetic field obtained from the entropy-stable scheme and the HLL scheme are compared preliminarily in Figure 11. As stated previously, the two schemes have no significant systematic differences and can produce a considerable match with in situ measurements of solar wind plasma parameters obtained from the OMNI database, with temperature and density a little overestimated and magnetic field somewhat underestimated. In order to improve the performance of the EC-GLM-MHD solar wind model, further validation and testing are needed in future.

Zoom In Zoom Out Reset image size

The main mission of this paper is to design the robust entropy-stable numerical approximate Riemann solver to solve the ideal EC-GLM-MHD equations and describe the solar wind plasma distributions and magnetic field topology from one solar radius (1 R_s) to Earth's orbit (1 au) for CR 2183. The developed model employs the nine-wave magnetic field divergence diminishing ideal MHD equations that incorporate the nonconservative Powell source terms to guarantee the Lorentz force being correctly captured when ∇ · B ≠ 0. This EC-GLM-MHD system is also symmetric and Galilean invariant with the nonconservative terms proportional to the gradient of the Lagrange multiplier ψ and the divergence of the magnetic field. On the other hand, the thermodynamical consistency and the entropy conservation (or stability) are accounted for. Due to the hyperbolic nature of the MHD system, discontinuities may develop during the solar wind simulation even if the initial solution is smooth. In the extreme environments, such as the low plasma β region near the Sun, nonphysical intermediate states or severe stability problems may appear within a numerical discontinuity profile. The entropy, introduced by the second law of thermodynamics, separates possible states from impossible ones and should be an auxiliary conservation law for an ideal MHD system. Here, we use the entropy-conservative flux function as a baseline flux and add a hybrid dissipation term to construct an entropy-stable numerical approximation, making it capable to simulate the low plasma β problem. The hybrid dissipation term continuously blends the highly accurate Roe dissipation scheme with the more dissipative LF dissipation scheme to increase the robustness of the entropy-stable numerical scheme in extreme ambient solar wind, while retaining more accurate features in the smooth computational domain. Additionally, this paper makes some efforts to preserve the positivity of pressure in several ways. On the one hand, we combine the non-negativity-preserving reconstruction and the solenoidality-preserving least-squares reconstruction to obtain the reliable interface values (Feng et al. 2017; Li et al. 2021). On the other hand, if the internal energy is lower than the predetermined threshold, the pressure derived from the formula of physical entropy density will be used to replace the one obtained from the energy equation, which is guaranteed to be positive (Derigs et al. 2016). From the numerical results, we can see that the entropy-stable scheme shows its obvious merits in the results of the two simple tests. Moreover, the results of the ambient solar wind from the Sun to Earth for CR 2183 show that the hybrid entropy-stable scheme is able to capture the large-scale structures of solar wind, such as coronal holes, regions characterized by open- and closed-magnetic fields, helmet streamers, pseudo-streamers, high-speed and low-speed solar wind flows, HCS, and so on. Meanwhile, in comparison with the remote measurements and the in situ interplanetary observations, the modeled results for the large-scale solar wind properties show that the entropy-stable EC-GLM-MHD model can efficiently handle complex structured solar wind and solar corona with low plasma β. Overall, based on the numerical results for solar wind simulation presented before, the EC-GLM-MHD solar wind model can be used as a candidate to reproduce the large-scale structures of solar wind to some extent. We believe that models, constrained by more observations of multiple spacecraft and established by taking more physically realistic processes such as coronal heating modules, will be bound to bring an accurate description of the corona−interplanetary space environment.

The work is jointly supported by the National Natural Science Foundation of China (grant Nos. 42030204, 41731067, 41874202, 41861164026, 41531073, and 42104168), the Basic and Applied Basic Research Foundation of Guangdong, China (grant No. 2019A1515111122), and the Specialized Research Fund for State Key Laboratories. The present work is partially supported by the B-type Strategic Priority Program of the Chinese Academy of Sciences, grant No. XDB41000000. This work utilizes data obtained by LASCO C2/SOHO, Solar Dynamics Observatory (SDO)/Atmospheric Imaging Assembly (AIA), and Helioseismic and Magnetic Imager (HMI). We also acknowledge use of NASA/GSFC's Space Physics Data Facility's OMNIWeb (or CDAWeb or ftp) service, and OMNI data. The work was carried out at the National Supercomputer Center in Tianjin, China, and the calculations were performed on TianHe-1 (A).

If the left and right states of positive scalar a are a_L and a_R , then the logarithmic mean of scalar a can be defined as (Ismail & Roe 2009)

$$\begin{eqnarray*}&&{a}^{\mathrm{ln}}(L,R)=\displaystyle \frac{{a}_{L}-{a}_{R}}{\mathrm{ln}{a}_{L}-\mathrm{ln}{a}_{R}}.\end{eqnarray*}$$

Apparently, the expression is not numerically well posed when a_L → a_R . To avoid this case, we make a transformation to write the logarithmic mean in another form, i.e., let $\zeta =\tfrac{{a}_{L}}{{a}_{R}}$, and then the logarithmic mean can be written also,

$$\begin{eqnarray}&&{a}^{\mathrm{ln}}(L,R)=\displaystyle \frac{{a}_{L}+{a}_{R}}{\mathrm{ln}\zeta }\displaystyle \frac{\zeta -1}{\zeta +1}\equiv \displaystyle \frac{{a}_{L}+{a}_{R}}{F}.\end{eqnarray} \tag{ A1 }$$

In order to express the logarithmic mean more clearly, we define $f=\tfrac{\zeta -1}{\zeta +1}$ and u = f².

Case I: If a_L → a_R , then $\zeta =\tfrac{{a}_{L}}{{a}_{R}}\to 1$, $f=\tfrac{\zeta -1}{\zeta +1}\to 0$, u = f² → 0. In this case, expression (A1) cannot be used directly. Utilizing the Taylor series expansion, the natural logarithm function $\mathrm{ln}\zeta$ may be expressed approximately as

$$\begin{eqnarray*}\begin{array}{rcl}\mathrm{ln}\zeta & = & 2\left(\displaystyle \frac{\zeta -1}{\zeta +1}+\displaystyle \frac{1}{3}{\left(\displaystyle \frac{\zeta -1}{\zeta +1}\right)}^{3}\right.\\ & & \left.+\displaystyle \frac{1}{5}{\left(\displaystyle \frac{\zeta -1}{\zeta +1}\right)}^{5}+\displaystyle \frac{1}{7}{\left(\displaystyle \frac{\zeta -1}{\zeta +1}\right)}^{7}+{ \mathcal O }({\zeta }^{9}\right).\end{array}\end{eqnarray*}$$

Then, we have

$$\begin{eqnarray*}&&\displaystyle \frac{1}{F}\equiv \displaystyle \frac{1}{\mathrm{ln}\zeta }\displaystyle \frac{\zeta -1}{\zeta +1}\approx \displaystyle \frac{1}{2\left(1+u/3+{u}^{2}/5+{u}^{3}/7\right)}\ \ \ (u\lt \varepsilon ).\end{eqnarray*}$$

Case II: If a_L ↛ a_R , then u ↛ 0. We obtain

$$\begin{eqnarray*}&&\displaystyle \frac{1}{F}\equiv \displaystyle \frac{1}{\mathrm{ln}\zeta }\displaystyle \frac{\zeta -1}{\zeta +1}=\displaystyle \frac{1}{\mathrm{ln}\zeta }f=\displaystyle \frac{1}{\mathrm{ln}\zeta /f}\ \ \ (\mathrm{else}).\end{eqnarray*}$$

Finally, the logarithmic mean may be written as

$$\begin{eqnarray*}{a}^{\mathrm{ln}}(L,R)=\left\{\begin{array}{cc}\displaystyle \frac{{a}_{L}+{a}_{R}}{2\left(1+u/3+{u}^{2}/5+{u}^{3}/7\right)} & (u\lt \varepsilon )\\ \displaystyle \frac{{a}_{L}+{a}_{R}}{\mathrm{ln}\zeta /f} & (\mathrm{else})\end{array}\right.\end{eqnarray*}$$

with ε = 10⁻³.

In order to ensure the relation $[\kern-2pt[ {\boldsymbol{U}}]\kern-2pt] ={ \mathcal H }[\kern-2pt[ { \mathcal V }]\kern-2pt]$ whenever possible, a choice of the averaging for the discrete entropy Jacobian ${ \mathcal H }$ can be as follows:

$$\begin{eqnarray}{ \mathcal H }=\left[\begin{array}{ccccccccc}{\rho }^{\mathrm{ln}} & {\rho }^{\mathrm{ln}}\{\{{v}_{r}\}\} & {\rho }^{\mathrm{ln}}\{\{{v}_{\theta }\}\} & {\rho }^{\mathrm{ln}}\{\{{v}_{\phi }\}\} & \overline{E} & 0 & 0 & 0 & 0\\ {\rho }^{\mathrm{ln}}\{\{{v}_{r}\}\} & {\rho }^{\mathrm{ln}}\{\{{v}_{r}\}\}{}^{2}+\tilde{p} & {\rho }^{\mathrm{ln}}\{\{{v}_{r}\}\}\{\{{v}_{\theta }\}\} & {\rho }^{\mathrm{ln}}\{\{{v}_{r}\}\}\{\{{v}_{\phi }\}\} & (\overline{E}+\tilde{p})\{\{{v}_{r}\}\} & 0 & 0 & 0 & 0\\ {\rho }^{\mathrm{ln}}\{\{{v}_{\theta }\}\} & {\rho }^{\mathrm{ln}}\{\{{v}_{r}\}\}\{\{{v}_{\theta }\}\} & {\rho }^{\mathrm{ln}}\{\{{v}_{\theta }\}\}{}^{2}+\tilde{p} & {\rho }^{\mathrm{ln}}\{\{{v}_{\theta }\}\}\{\{{v}_{\phi }\}\} & (\overline{E}+\tilde{p})\{\{{v}_{\theta }\}\} & 0 & 0 & 0 & 0\\ {\rho }^{\mathrm{ln}}\{\{{v}_{\phi }\}\} & {\rho }^{\mathrm{ln}}\{\{{v}_{r}\}\}\{\{{v}_{\phi }\}\} & {\rho }^{\mathrm{ln}}\{\{{v}_{\theta }\}\}\{\{{v}_{\phi }\}\} & {\rho }^{\mathrm{ln}}\{\{{v}_{\phi }\}\}{}^{2}+\tilde{p} & (\overline{E}+\tilde{p})\{\{{v}_{\phi }\}\} & 0 & 0 & 0 & 0\\ \overline{E} & (\overline{E}+\tilde{p})\{\{{v}_{r}\}\} & (\overline{E}+\tilde{p})\{\{{v}_{\theta }\}\} & (\overline{E}+\tilde{p})\{\{{v}_{\phi }\}\} & {{ \mathcal H }}_{\mathrm{5,5}} & \tau \{\{{B}_{r}\}\} & \tau \{\{{B}_{\theta }\}\} & \tau \{\{{B}_{\phi }\}\} & \tau \{\{\psi \}\}\\ 0 & 0 & 0 & 0 & \tau \{\{{B}_{r}\}\} & \tau & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \tau \{\{{B}_{\theta }\}\} & 0 & \tau & 0 & 0\\ 0 & 0 & 0 & 0 & \tau \{\{{B}_{\phi }\}\} & 0 & 0 & \tau & 0\\ 0 & 0 & 0 & 0 & \tau \{\{\psi \}\} & 0 & 0 & 0 & \tau \end{array}\right]\end{eqnarray} \tag{ B1 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl}\overline{E} & = & \displaystyle \frac{{p}^{\mathrm{ln}}}{\gamma -1}+\displaystyle \frac{1}{2}{\rho }^{\mathrm{ln}}\overline{{{\boldsymbol{v}}}^{2}},\ {p}^{\mathrm{ln}}=\displaystyle \frac{{\rho }^{\mathrm{ln}}}{2{\varsigma }^{\mathrm{ln}}},\\ \tau & = & \displaystyle \frac{\{\{p\}\}}{\{\{\rho \}\}}=\displaystyle \frac{1}{2\{\{\varsigma \}\}},\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}\overline{{{\boldsymbol{v}}}^{2}} & = & 2\left(\{\{{v}_{r}\}\}{}^{2}+\{\{{v}_{\theta }\}\}{}^{2}+\{\{{v}_{\phi }\}\}{}^{2}\right)\\ & & -\left(\{\{{v}_{r}^{2}\}\}+\{\{{v}_{\theta }^{2}\}\}+\{\{{v}_{\phi }^{2}\}\}\right),\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}{{ \mathcal H }}_{\mathrm{5,5}} & = & \displaystyle \frac{1}{{\rho }^{\mathrm{ln}}}\left(\displaystyle \frac{{\left({p}^{\mathrm{ln}}\right)}^{2}}{\gamma -1}+{\overline{E}}^{2}\right)+\tilde{p}\left(\{\{{v}_{r}\}\}{}^{2}+\{\{{v}_{\theta }\}\}{}^{2}\right.\\ & & \left.+\{\{{v}_{\phi }\}\}{}^{2}\right)+\tau \left(\{\{{B}_{r}\}\}{}^{2}+\{\{{B}_{\theta }\}\}{}^{2}\right.\\ & & \left.+\{\{{B}_{\phi }\}\}{}^{2}+\{\{\psi \}\}{}^{2}\right).\end{array}\end{eqnarray*}$$

The items of each discrete right eigenvectors are as follows:

GLM wave: λ_±ψ

$$\begin{eqnarray*}&&{{\boldsymbol{r}}}_{\pm \psi }={\left[\mathrm{0,0},\mathrm{0,0},\{\{{B}_{r}\}\}\pm \{\{\psi \}\},\mathrm{1,0,0},\pm 1\right]}^{T}.\end{eqnarray*}$$

Entropy wave: λ_E

$$\begin{eqnarray*}&&{{\boldsymbol{r}}}_{E}={\left[1,\{\{{v}_{r}\}\},\{\{{v}_{\theta }\}\},\{\{{v}_{\phi }\}\},\displaystyle \frac{1}{2}\overline{{{\boldsymbol{v}}}^{2}},0,0,0,0\right]}^{T}.\end{eqnarray*}$$

Alfvén waves: λ_±a

$$\begin{eqnarray*}&&{{\boldsymbol{r}}}_{\pm a}=\left[0,0,\pm {\rho }^{\mathrm{ln}}\sqrt{\{\{\rho \}\}}{\bar{\chi }}_{\phi },\mp {\rho }^{\mathrm{ln}}\sqrt{\{\{\rho \}\}}{\bar{\chi }}_{\theta },\right.\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\left.\mp {\rho }^{\mathrm{ln}}\sqrt{\{\{\rho \}\}}\left({\bar{\chi }}_{\theta }\{\{{v}_{\phi }\}\}-{\bar{\chi }}_{\phi }\{\{{v}_{\theta }\}\}\right),0,-{\rho }^{\mathrm{ln}}{\bar{\chi }}_{\phi },{\rho }^{\mathrm{ln}}{\bar{\chi }}_{\theta },0\right]}^{T}.\end{eqnarray*}$$

Fast magnetoacoustic waves: λ_±f

$$\begin{eqnarray*}&&\begin{array}{l}{{\boldsymbol{r}}}_{\pm f}=\left[{\hat{\alpha }}_{f}{\rho }^{\mathrm{ln}},{\hat{\alpha }}_{f}{\rho }^{\mathrm{ln}}\left(\{\{{v}_{r}\}\}\pm {\hat{c}}_{f}\right),\right.\\ {\rho }^{\mathrm{ln}}\left({\hat{\alpha }}_{f}\{\{{v}_{\theta }\}\}\mp {\hat{\alpha }}_{s}{\hat{c}}_{s}{\bar{\chi }}_{\theta }\sigma ({\bar{b}}_{r})\right),\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}&&\begin{array}{l}{\rho }^{\mathrm{ln}}\left({\hat{\alpha }}_{f}\{\{{v}_{\phi }\}\}\mp {\hat{\alpha }}_{s}{\hat{c}}_{s}{\bar{\chi }}_{\phi }\sigma ({\bar{b}}_{r})\right),{\hat{{\rm{\Psi }}}}_{\pm f},0,\\ {\left.{\hat{\alpha }}_{s}{a}^{\varsigma }{\bar{\chi }}_{\theta }\sqrt{{\rho }^{\mathrm{ln}}},{\hat{\alpha }}_{s}{a}^{\varsigma }{\bar{\chi }}_{\phi }\sqrt{{\rho }^{\mathrm{ln}}},0\right]}^{T}.\end{array}\end{eqnarray*}$$

Slow magnetoacoustic waves: λ_±s

$$\begin{eqnarray*}&&\begin{array}{l}{{\boldsymbol{r}}}_{\pm s}=\left[{\hat{\alpha }}_{s}{\rho }^{\mathrm{ln}},{\hat{\alpha }}_{s}{\rho }^{\mathrm{ln}}\left(\{\{{v}_{r}\}\}\pm {\hat{c}}_{s}\right),\right.\\ {\rho }^{\mathrm{ln}}\left({\hat{\alpha }}_{s}\{\{{v}_{\theta }\}\}\pm {\hat{\alpha }}_{f}{\hat{c}}_{f}{\bar{\chi }}_{\theta }\sigma ({\bar{b}}_{r})\right),\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}&&\begin{array}{l}{\rho }^{\mathrm{ln}}\left({\hat{\alpha }}_{s}\{\{{v}_{\phi }\}\}\pm {\hat{\alpha }}_{f}{\hat{c}}_{f}{\bar{\chi }}_{\phi }\sigma ({\bar{b}}_{r})\right),{\hat{{\rm{\Psi }}}}_{\pm s},0,\\ {\left.-{\hat{\alpha }}_{f}{a}^{\varsigma }{\bar{\chi }}_{\theta }\sqrt{{\rho }^{\mathrm{ln}}},-{\hat{\alpha }}_{f}{a}^{\varsigma }{\bar{\chi }}_{\phi }\sqrt{{\rho }^{\mathrm{ln}}},0\right]}^{T},\end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{\hat{{\rm{\Psi }}}}_{\pm s} & = & \displaystyle \frac{{\hat{\alpha }}_{s}{\rho }^{\mathrm{ln}}\overline{{{\boldsymbol{v}}}^{2}}}{2}-{a}^{\varsigma }{\hat{\alpha }}_{f}{\rho }^{\mathrm{ln}}{\bar{b}}_{\perp }+\displaystyle \frac{{\hat{\alpha }}_{s}{\rho }^{\mathrm{ln}}{\left({a}^{\mathrm{ln}}\right)}^{2}}{\gamma -1}\pm {\hat{\alpha }}_{s}{\hat{c}}_{s}{\rho }^{\mathrm{ln}}\{\{{v}_{r}\}\}\\ & & \pm \,{\hat{\alpha }}_{f}{\hat{c}}_{f}{\rho }^{\mathrm{ln}}\sigma ({\bar{b}}_{r})\left(\{\{{v}_{\theta }\}\}{\bar{\chi }}_{\theta }+\{\{{v}_{\phi }\}\}{\bar{\chi }}_{\phi }\right)\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}&&\begin{array}{l}{\hat{{\rm{\Psi }}}}_{\pm f}=\displaystyle \frac{{\hat{\alpha }}_{f}{\rho }^{\mathrm{ln}}\overline{{{\boldsymbol{v}}}^{2}}}{2}+{a}^{\varsigma }{\hat{\alpha }}_{s}{\rho }^{\mathrm{ln}}{\bar{b}}_{\perp }\\ \quad +\displaystyle \frac{{\hat{\alpha }}_{f}{\rho }^{\mathrm{ln}}{\left({a}^{\mathrm{ln}}\right)}^{2}}{\gamma -1}\pm {\hat{\alpha }}_{f}{\hat{c}}_{f}{\rho }^{\mathrm{ln}}\{\{{v}_{r}\}\}\\ \quad \mp {\hat{\alpha }}_{s}{\hat{c}}_{s}{\rho }^{\mathrm{ln}}\sigma ({\bar{b}}_{r})\left(\{\{{v}_{\theta }\}\}{\bar{\chi }}_{\theta }+\{\{{v}_{\phi }\}\}{\bar{\chi }}_{\phi }\right)\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}{\hat{c}}_{a}^{2} & = & {\bar{b}}_{r}^{2},\ {\hat{c}}_{f,s}^{2}=\displaystyle \frac{1}{2}\left(({\bar{a}}^{2}+{\bar{b}}^{2})\pm \sqrt{{\left({\bar{a}}^{2}+{\bar{b}}^{2}\right)}^{2}-4{\bar{a}}^{2}{\bar{b}}_{r}^{2}}\right),\\ \tilde{p} & = & \displaystyle \frac{\{\{\rho \}\}}{2\{\{\varsigma \}\}},\ {\bar{a}}^{2}=\gamma \displaystyle \frac{\tilde{p}}{{\rho }^{\mathrm{ln}}}\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}{\left({a}^{\mathrm{ln}}\right)}^{2} & = & \gamma \displaystyle \frac{{p}^{\mathrm{ln}}}{{\rho }^{\mathrm{ln}}},\ {\left({a}^{\varsigma }\right)}^{2}=\gamma \displaystyle \frac{1}{2\{\{\varsigma \}\}},\\ {\bar{b}}^{2} & = & {\bar{b}}_{r}^{2}+{\bar{b}}_{\theta }^{2}+{\bar{b}}_{\phi }^{2},\\ {\bar{b}}_{\perp }^{2} & = & {\bar{b}}_{\theta }^{2}+{\bar{b}}_{\phi }^{2},\ {\bar{\chi }}_{\theta ,\phi }=\displaystyle \frac{{\bar{b}}_{\theta ,\phi }}{{\bar{b}}_{\perp }}\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}{\bar{b}}_{r,\theta ,\phi }^{2} & = & \displaystyle \frac{\{\{{B}_{r,\theta ,\phi }\}\}{}^{2}}{{\rho }^{\mathrm{ln}}},\\ {\hat{\alpha }}_{f}^{2} & = & \displaystyle \frac{{\bar{a}}^{2}-{\hat{c}}_{s}^{2}}{{\hat{c}}_{f}^{2}-{\hat{c}}_{s}^{2}},\ {\hat{\alpha }}_{s}^{2}=\displaystyle \frac{{\hat{c}}_{f}^{2}-{\bar{a}}^{2}}{{\hat{c}}_{f}^{2}-{\hat{c}}_{s}^{2}},\\ \sigma (\omega ) & = & \left\{\begin{array}{ll}+1 & \mathrm{if}\ \ \omega \geqslant 0,\\ -1 & \mathrm{otherwise}\end{array}\right..\end{array}\end{eqnarray*}$$

The average state for the diagonal scaling matrix Z is

$$\begin{eqnarray*}&&\begin{array}{l}{\boldsymbol{Z}}=\mathrm{diag}\left(\displaystyle \frac{1}{2\gamma {\rho }^{\mathrm{ln}}},\right.\\ \displaystyle \frac{1}{4\{\{\varsigma \}\}{\left({\rho }^{\mathrm{ln}}\right)}^{2}},\displaystyle \frac{1}{2\gamma {\rho }^{\mathrm{ln}}},\displaystyle \frac{1}{4\{\{\varsigma \}\}},\displaystyle \frac{{\rho }^{\mathrm{ln}}(\gamma -1)}{\gamma },\\ \left.\displaystyle \frac{1}{4\{\{\varsigma \}\}},\displaystyle \frac{1}{2\gamma {\rho }^{\mathrm{ln}}},\displaystyle \frac{1}{2\gamma {\rho }^{\mathrm{ln}}},\displaystyle \frac{1}{4\{\{\varsigma \}\}{\left({\rho }^{\mathrm{ln}}\right)}^{2}},\displaystyle \frac{1}{2\gamma {\rho }^{\mathrm{ln}}}\right),\end{array}\end{eqnarray*}$$

and the average state of the diagonal matrix reads

$$\begin{eqnarray}&&{\boldsymbol{\Lambda }}=\mathrm{diag}\left({\hat{\lambda }}_{+f},{\hat{\lambda }}_{+a},{\hat{\lambda }}_{+s},{\hat{\lambda }}_{+\psi },{\hat{\lambda }}_{E},{\hat{\lambda }}_{-\psi },{\hat{\lambda }}_{-s},{\hat{\lambda }}_{-a},{\hat{\lambda }}_{-f}\right)\end{eqnarray} \tag{ C1 }$$

from the nine eigenvalues for the ideal EC-GLM-MHD system evaluated by

$$\begin{eqnarray*}&&\left[\begin{array}{c}{\hat{\lambda }}_{+f}\\ {\hat{\lambda }}_{+a}\\ {\hat{\lambda }}_{+s}\\ {\hat{\lambda }}_{+\psi }\\ {\hat{\lambda }}_{E}\\ {\hat{\lambda }}_{-\psi }\\ {\hat{\lambda }}_{-s}\\ {\hat{\lambda }}_{-a}\\ {\hat{\lambda }}_{-f}\end{array}\right]=\left[\begin{array}{c}\{\{{v}_{r}\}\}+{\hat{c}}_{f}\\ \{\{{v}_{r}\}\}+{\hat{c}}_{a}\\ \{\{{v}_{r}\}\}+{\hat{c}}_{s}\\ \{\{{v}_{r}\}\}+\{\{{c}_{h}\}\}\\ \{\{{v}_{r}\}\}\\ \{\{{v}_{r}\}\}-\{\{{c}_{h}\}\}\\ \{\{{v}_{r}\}\}-{\hat{c}}_{s}\\ \{\{{v}_{r}\}\}-{\hat{c}}_{a}\\ \{\{{v}_{r}\}\}-{\hat{c}}_{f}\end{array}\right].\end{eqnarray*}$$

The form of the discrete wave speeds ${\hat{c}}_{f,a,s}$ is

$$\begin{eqnarray*}&&\begin{array}{l}{\hat{c}}_{f,s}=\frac{\sqrt{2\sqrt{{\hat{b}}_{r}\hat{a}}+\hat{a}+{\hat{b}}_{\phi }+{\hat{b}}_{\theta }+{\hat{b}}_{r}}\pm \sqrt{-2\sqrt{{\hat{b}}_{r}\hat{a}}+\hat{a}+{\hat{b}}_{\phi }+{\hat{b}}_{\theta }+{\hat{b}}_{r}}}{2},\\ {\hat{c}}_{a}=\sqrt{{\hat{b}}_{r}},\end{array}\end{eqnarray*}$$

where $\hat{a}=\{\{p\}\}\{\{{\rho }^{-1}\}\}\gamma$, ${\hat{b}}_{\kappa }=\{\{{B}_{\kappa }\}\}\{\{\tfrac{{B}_{\kappa }}{\rho }\}\}$, with κ = r, θ, ϕ.