Application of a numerical path-preserving Roe scheme in background solar wind simulation

1 Introduction

Space weather encompasses a series of physical processes driven by solar activity that affect the Sun-Earth space environment, including solar flares, coronal mass ejections (CMEs), solar energetic particle events, and solar wind disturbances. Intense space weather events can impact satellite operations, communication and navigation systems, aviation, power grids, and human spaceflight, making the investigation of their underlying physical mechanisms and the development of predictive numerical methods central challenges in heliophysics and space environment engineering. In contrast to transient solar eruptions such as CMEs, the background solar wind can be regarded as a quasi-steady, magnetized plasma flow that fills the heliosphere and establishes the fundamental plasma and magnetic field environment of interplanetary space. Transient disturbances, including CMEs and interplanetary shocks, propagate through and interact with this background, with their speed, direction, arrival time, and geomagnetic impact being modulated by the solar wind’s velocity distribution, density structure, and magnetic topology. Accurate modeling of the large-scale structure of the background solar wind is therefore essential for understanding Sun-Earth coupling processes and improving space weather forecasting (Feng, 2020).

The solar wind was first explained by Parker as a thermally driven, supersonic expansion of coronal plasma into interplanetary space (Parker, 1958; 1963; Wu et al., 2023). Due to the high conductivity of the solar wind, the solar magnetic field is carried outward with the plasma, forming the interplanetary magnetic field and shaping large-scale heliospheric structures (Burlaga, 1995; Schwenn and Marsch, 1990). The background solar wind is not a uniform radial flow but exhibits structured distributions that vary with source region magnetic topology and solar activity. During solar minimum, high-speed streams originate primarily from polar coronal holes, while low-speed wind is associated with low-latitude regions and streamer flows, and the heliospheric current sheet exhibits a wavy equatorial pattern. During solar maximum, multiple source regions and transient events intertwine, producing more complex spatial structures at all heliolatitudes (Verscharen et al., 2019). These nonuniform background structures collectively establish the heliospheric plasma and magnetic environment and provide the physical context for CME propagation, shock formation, and geomagnetic responses.

Due to the limited constraints provided by remote sensing and in situ observations, three-dimensional numerical simulations have become indispensable for studying and forecasting the background solar wind. MHD models self-consistently describe the interactions between large-scale magnetic fields and plasma flows and offer a practical compromise between physical fidelity and computational feasibility, making them a foundational framework for modeling the background solar wind and Sun-interplanetary space weather (Feng, 2020). Ideal MHD equations, typically neglecting resistive, viscous, and thermal dissipation effects, can be expressed in hyperbolic conservation-law form and must satisfy the divergence-free condition of the magnetic field. More sophisticated coronal-heliospheric MHD models may include additional physical processes such as solar gravity, coronal heating, solar wind acceleration, thermal conduction, and Alfvén-wave turbulence heating (Usmanov et al., 2000; Groth et al., 2000; van der Holst et al., 2010; van der Holst et al., 2014).

Over the years, various three-dimensional, time-dependent MHD modeling frameworks have been developed for the background solar wind and Sun-interplanetary space weather. These representative modeling approaches can be broadly categorized according to their modeling strategies and simulation domains. Coronal and coupled coronal–heliospheric MHD frameworks, such as CISM/CORHEL (Riley et al., 2012a; Goodrich et al., 2004), SWMF/BATS-R-US (Powell et al., 1999; Tóth et al., 2012), and SWIM/SIP–CESE MHD (Feng et al., 2007, Feng et al., 2010; Feng et al., 2014a), solve the MHD equations in the corona and/or couple coronal and heliospheric domains to describe the large-scale solar wind from the near-Sun region into interplanetary space. Hybrid empirical–MHD models, such as SUSANOO (Shiota and Kataoka, 2016) and EUHFORIA (Pomoell and Poedts, 2018), combine empirical or semi-empirical coronal descriptions with heliospheric MHD simulations. Representative heliospheric MHD models include ENLIL, which is widely used to simulate the large-scale solar-wind structure and the propagation of interplanetary disturbances (Odstrcil, 2003), as well as the more recently developed sunRunner3D (González-Avilés et al., 2024) and SWASTi-SW (Mayank et al., 2022), which provide additional numerical frameworks for global heliospheric solar-wind simulations. Together, these representative modeling frameworks provide important tools for reconstructing coronal magnetic structures, simulating the formation and propagation of the background solar wind, modeling heliospheric evolution, and tracking interplanetary disturbances. While existing models reproduce large-scale solar wind structures reasonably well, fine-scale features such as the heliospheric current sheet, fast-slow stream interaction regions, and regions with rapid magnetic field rotations remain sensitive to factors including photospheric magnetic field input, inner boundary conditions, treatment of physical source terms, grid resolution, and numerical discretization (Feng, 2020).

In numerical discretization, the finite volume method is widely employed for ideal MHD and background solar wind simulations due to its local conservation properties and strong capability to capture discontinuities (LeVeque, 2002; Toro, 2009). Within this framework, the calculation of numerical fluxes at cell interfaces is crucial, as it determines wave propagation characteristics and dissipation. Fine-scale structures in the background solar wind, including current sheets and fast-slow stream boundaries, often manifest as thin layers or steep gradients, making the choice of approximate Riemann solvers critical. Common solvers include local Lax-Friedrichs (Rusanov) (Rusanov, 1962), Harten-Lax-van Leer-Einfeldt (HLLE) (Einfeldt et al., 1991; Harten et al., 1983), artificial wind methods (Sokolov et al., 2002), Harten–Lax–van Leer–Discontinuities (HLLD) (Miyoshi and Kusano, 2005), Godunov fluxes (Godunov, 1959), and Roe (Roe, 1981). While Rusanov and HLLE are robust, their limited wave-speed information results in higher numerical dissipation and poor resolution of contact and rotational discontinuities. HLLD improves upon this by restoring key intermediate MHD wave structures, achieving a balance between efficiency and resolution. The Roe solver, with local linearization and complete characteristic decomposition, can resolve fast and slow magnetosonic waves, Alfvén waves, and entropy waves individually, providing low numerical dissipation and high discontinuity resolution.

However, ideal MHD equations are non-strictly hyperbolic and non-convex, rendering accurate discontinuity capture challenging. In MHD Riemann problems with regular solutions but significant tangential magnetic field rotations, conventional solvers may produce nonphysical intermediate shocks absent in the exact solution, which decay only slowly with mesh refinement, resulting in pseudo-convergence (Torrilhon, 2003; Torrilhon and Balsara, 2004). This phenomenon is closely linked to the incorrect capturing of Alfvén waves, whose tangential magnetic field rotates while its magnitude remains constant (Xu et al., 2023, Xu et al., 2025). Numerical methods failing to preserve this property can misrepresent Alfvén waves as intermediate shocks, affecting the propagation of other characteristic waves. Based on this insight, the Numerical Path-Preserving (NPP) approach was proposed and extended to the Roe solver as the NPP-Roe method, which corrects Alfvén wave characteristics to ensure numerical solutions evolve along physically consistent paths, mitigating rotational discontinuity errors and pseudo-convergence (Xu et al., 2025).

Despite its demonstrated effectiveness in canonical MHD problems, the applicability of NPP-Roe to complex background solar wind simulations featuring large-scale magnetic topology, fast and slow stream interactions, current sheets, and rapid magnetic field rotations remains to be assessed. Applying NPP-Roe in this context provides a test of its ability to resolve magnetic structures and discontinuities while preserving physical consistency. Accordingly, this study implements the NPP-Roe solver for background solar wind modeling and evaluates its performance in large-scale solar wind MHD simulations.

The remainder of the paper is organized as follows: Section 2 describes the ideal MHD governing equations, grid system, and boundary conditions; Section 3 presents the finite volume discretization, least-squares reconstruction, traditional Roe flux, and NPP-Roe modifications for Alfvén wave characteristics; Section 4 first compares traditional Roe and NPP-Roe results in one-dimensional MHD Riemann problem and then analyzes three-dimensional simulations of the background solar wind during CR 2203, focusing on large-scale coronal structures, such as solar wind velocity and magnetic topology; finally, conclusions, discussions, and potential directions for further improvement are presented. Supplementary Appendix A provides a simple comparison between the classical Roe and NPP-Roe global background solar-wind solutions obtained under otherwise identical conditions.

2 Numerical methods

2.1 Governing equations

The background solar wind is modeled using the modified MHD equations. To suppress nonphysical phenomena such as negative pressure that may arise during the numerical computation, the magnetic field B = ( B x , B y , B z ) T is decomposed into two components, B 0 and B 1. Here, B 0 = ( B x 0 , B y 0 , B z 0 ) T represents the time-independent background magnetic field, while B 1 = ( B x 1 , B y 1 , B z 1 ) T denotes the perturbation component that evolves with time (Ogino and Walker, 1984; Tanaka, 1994). It should be noted that the magnetic permeability factor is absorbed into B, B 1, and B 0, which are expressed as the corresponding physical magnetic fields divided by μ, where μ = 4 π × 1 0 − 7 H m − 1. This decomposition approach preserves the constraint of the background magnetic field while ensuring improved numerical stability and physical consistency in the evolution of the perturbation magnetic field. Furthermore, the generalized Lagrange multiplier (GLM) divergence-cleaning method proposed by Dedner et al. (2002) is adopted to reduce numerical violations of the magnetic-field divergence-free condition, ∇ ⋅ B = 0. This method introduces an additional constraint variable ψ into the conventional MHD equation system, effectively suppressing nonphysical divergence errors that may arise during numerical computations, thereby improving the stability and consistency of the simulation results. A key advantage of the GLM formulation is that the augmented induction system retains the conservative hyperbolic form of the governing equations, allowing it to be treated directly using existing finite-volume schemes without modifying the reconstruction, numerical-flux, or time-integration procedures. The ideal GLM-MHD governing equations after magnetic field splitting are therefore written as follows (Feng et al., 2014b):where u = ( u x , u y , u z ) T is the plasma velocity, ρ is the plasma mass density, p is the thermal pressure, I is the unit tensor, Ω is the prescribed solar rotation angular velocity, taken as Ω = 2.865 × 1 0 − 6 r a d s − 1, corresponding to a Carrington sidereal rotation period of 25.38 days and ψ is the generalized Lagrange multiplier variable used to control the propagation and damping of magnetic divergence errors. G and M s denote the gravitational constant and the solar mass, respectively. In general,G M s = 1.328 × 1 0 20 m 3 s − 2. The parameter c h controls the propagation speed of the divergence error and is defined as the maximum magnetosonic signal speed within the computational domain in order to avoid additional time-step restrictions. Thus, c h = max i , j , k ( | u x | + c f x , | u y | + c f y , | u z | + c f z ). At the same time, the divergence error decays proportionally to the rate c h 2 c p 2. The damping coefficient c p is determined by the dimensionless parameter α = Δ h c h c p 2 ∈ [ 0,1 ]. This parameter α represents the ratio between the diffusion and advection time scales, i.e.,

∂ ρ ∂ t + ∇ ⋅ ρ u = 0 ∂ ρ u ∂ t + ∇ ⋅ ρ u u + p + B 1 2 2 + B 1 ⋅ B 0 I − B 1 B 1 − B 1 B 0 − B 0 B 1 = − ρ G M s r 3 r − 2 ρ Ω × u − ρ Ω × Ω × r + S M ∂ B 1 ∂ t + ∇ ⋅ u B − B u + ψ I = 0 ∂ E 1 ∂ t + ∇ ⋅ E 1 + p + B 1 2 2 + B 1 B 0 u − B 1 u ⋅ B 1 − B 0 u ⋅ B 1 = − ρ u ⋅ G M s r 2 r + Ω × Ω × r + Q E + u ⋅ S M ∂ ψ ∂ t + ∇ ⋅ c h 2 B 1 = − c h 2 c p 2 ψ (1)

α = Δ t d Δ t a , Δ t d = Δ h 2 c p 2 , Δ t a = Δ h c h .

Here, Δ h denotes the minimum grid spacing. In practical applications, we set α = 0.1. Numerical experiments show that this choice can achieve robust divergence cleaning without introducing excessive diffusion. In addition, the modified total energy density is defined aswhere γ is the effective polytropic index prescribed during the MHD simulation following Wu and Guo (1999) and Feng et al. (2010). It is a time-independent function of heliocentric distance, given by

E 1 = p γ − 1 + 1 2 ρ u 2 + 1 2 B 1 2 ,

γ = 1.05 , r / R s ≤ 5 , 1.05 + 0.03 r / R s − 5 , 5 < r / R s ≤ 22 , 1.56 , r / R s > 22 .

The variable r represents the radial distance from the center of the Sun to a given point, while R s denotes the solar radius. The quantity r / R s is dimensionless and denotes the ratio of the distance from a point to the solar center to the solar radius. Within the momentum source term and the energy source term, S M and Q E represent the components associated with solar-wind acceleration and heating, respectively. They can be written as where Q 2 = Q 0 C a and M = M 0 C a. The topology factor C a is defined aswiththe magnetic expansion factor, R s s = 2.5 R s the source-surface radius, and θ b the minimum angular distance between the footpoint of an open field line and the nearest coronal-hole boundary at the photosphere. We set

S M = M r R s − 1 exp − r L M ,

Q E = Q 1 ⁡ exp − r L Q 1 + Q 2 r R s − 1 exp − r L Q 2 ,

C a = C a ′ max C a ′ , C a ′ = 1 − 0.8 e − θ b / 2.5 2 1 + f s 2 / 9 ,

f s = R s R s s 2 B R s B R s s

and take

M 0 = 3.5 × 1 0 − 13 N m − 3 , Q 0 = 1.18 × 1 0 − 7 J m − 3 s − 1 , Q 1 = 7.5 × 1 0 − 9 J m − 3 s − 1 ,

L Q 1 = L Q 2 = L M = 1 R s .

To improve numerical stability, all variables are nondimensionalized before calculation. In this study, ρ, u, p, B, r, and t are normalized by ρ s, a s, ρ s a s 2, ρ s a s 2, R s, and R s / a s, respectively, where ρ s is the solar-surface density, a s is the solar-surface sound speed. This treatment makes variables dimensionless and comparable in magnitude, and the dimensional results can be recovered through denormalization.

In Cartesian coordinates, Equation 1 can be rewritten in the following conservative form:

where the conserved variable vector U, the flux vectors F x, F y, and F z, and the source term S in Equation 2 are defined as:

∂ U ∂ t + ∂ F x ∂ x + ∂ F y ∂ y + ∂ F z ∂ z = S , (2)

U = ρ , ρ u x , ρ u y , ρ u z , E 1 , B x 1 , B y 1 , B z 1 , ψ T

F x = ρ u x ρ u x 2 + p + 1 2 − B x 1 2 + B y 1 2 + B z 1 2 + − B x 0 B x 1 + B y 0 B y 1 + B z 0 B z 1 ρ u x u y − B x 1 B y 1 − B x 1 B y 0 − B x 0 B y 1 ρ u x u z − B x 1 B z 1 − B x 1 B z 0 − B x 0 B z 1 E 1 + p + B 1 2 2 + B 0 ⋅ B 1 u x − B x u ⋅ B 1 ψ u x B y − B x u y u x B z − B x u z c h 2 B x 1

F y = ρ u y ρ u y u x − B y 1 B x 1 − B y 1 B x 0 − B y 0 B x 1 ρ u y 2 + p + 1 2 B x 1 2 − B y 1 2 + B z 1 2 + B x 0 B x 1 − B y 0 B y 1 + B z 0 B z 1 ρ u y u z − B y 1 B z 1 − B y 1 B z 0 − B y 0 B z 1 E 1 + p + B 1 2 2 + B 0 ⋅ B 1 u y − B y u ⋅ B 1 u y B x − B y u x ψ u y B z − B y u z c h 2 B y 1

F z = ρ u z ρ u z u x − B z 1 B x 1 − B z 1 B x 0 − B z 0 B x 1 ρ u z u y − B z 1 B y 1 − B z 1 B y 0 − B z 0 B y 1 ρ u z 2 + p + 1 2 B x 1 2 + B y 1 2 − B z 1 2 + B x 0 B x 1 + B y 0 B y 1 − B z 0 B z 1 E 1 + p + B 1 2 2 + B 0 ⋅ B 1 u z − B z u ⋅ B 1 u z B x − B z u x u z B y − B z u y ψ c h 2 B z 1

S = 0 − ρ G M s r 2 r x − 2 ρ Ω × u x − ρ Ω × Ω × r x + S M , x − ρ G M s r 2 r y − 2 ρ Ω × u y − ρ Ω × Ω × r y + S M , y − ρ G M s r 2 r z − 2 ρ Ω × u z − ρ Ω × Ω × r z + S M , z − ρ u ⋅ G M s r 2 r − ρ u ⋅ Ω × Ω × r + u ⋅ S M + Q E 0 0 0 − c h 2 c p 2 ψ

To obtain a more compact representation of the governing equations, we introduce the primitive variable vector W = ( ρ , u x , u y , u z , p , B x 1 , B y 1 , B z 1 , ψ ) T. Compared with the conserved variable vector U, these variables provide a more intuitive description of the system. Equation 1 can then be recast in the following quasi-linear form:

∂ W ∂ t + A x ∂ W ∂ x + A y ∂ W ∂ y + A z ∂ W ∂ z = S W

where A x, A y, and A z are the Jacobian matrices in the x-, y-, and z-directions, respectively. For brevity, only the Jacobian matrix in the x-direction is given as

A x = u x ρ 0 0 0 0 0 0 0 0 u x 0 0 1 ρ − B x ρ B y ρ B z ρ 0 0 0 u x 0 0 − B y ρ − B x ρ 0 0 0 0 0 u x 0 − B z ρ 0 − B x ρ 0 0 γ p 0 0 u x γ − 1 u ⋅ B 0 0 1 − γ B x 0 0 0 0 0 0 0 0 1 0 B y − B x 0 0 − u y u x 0 0 0 B z 0 − B x 0 − u z 0 u x 0 0 0 0 0 0 c h 2 0 0 0

It is worth noting that, due to the directional symmetry of the Cartesian MHD equations, A y and A z can be obtained from A x by the corresponding permutation of the velocity and magnetic-field components, together with the associated rows and columns of the Jacobian matrix. Furthermore,the relationship between the conservative and primitive variables can be formally expressed through the following partial derivatives: where k = ( 1 − γ ).

∂ U ∂ W = 1 0 0 0 0 0 0 0 0 u x ρ 0 0 0 0 0 0 0 u y 0 ρ 0 0 0 0 0 0 u z 0 0 ρ 0 0 0 0 0 u ⋅ u 2 ρ u x ρ u y ρ u z − 1 k B x 1 B y 1 B z 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1

∂ W ∂ U = 1 0 0 0 0 0 0 0 0 − u x ρ 1 ρ 0 0 0 0 0 0 0 − u y ρ 0 1 ρ 0 0 0 0 0 0 − u z ρ 0 0 1 ρ 0 0 0 0 0 − k u ⋅ u 2 k u x k u y k u z − k k B x 1 k B y 1 k B z 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1

2.2 Grid system

The computational domain is a spherical shell. To avoid the singularity and grid clustering near the poles that arise in conventional spherical coordinates, the six-patch grid system proposed by Feng et al. (2010), Feng et al., 2014a, Feng et al., 2014b) is adopted, as illustrated in Figure 1. In this grid system, the spherical surface is divided into six identical, partially overlapping patches. The overlap regions facilitate inter-patch communication and boundary treatment. Since all patches have the same size and can be generated by rotating a reference patch, this grid structure is convenient for parallel implementation.

During grid generation, each patch is defined in its own local coordinate system. By transforming the local coordinates to the global coordinate system, the full computational mesh can be constructed. Each patch corresponds to a low-latitude spherical grid, whose coordinate ranges in spherical coordinates are given by

π 4 − δ ≤ θ ≤ 3 π 4 + δ ∩ 3 π 4 − δ ≤ ϕ ≤ 5 π 4 + δ .

The parameter δ is proportional to the grid spacing and is introduced to ensure a minimum overlap between adjacent patches, thereby enabling stable and efficient information exchange. The overlap region corresponds to at least two grid spacings, i.e., δ = 2 Δ θ. For the spatial discretization of ( θ , ϕ ), the grid points on each patch are defined aswhere Δ θ = θ max − θ min N θ , Δ ϕ = ϕ max − ϕ min N ϕ . The quantities N θ and N ϕ denote the number of grid points in the latitudinal and longitudinal directions, respectively, and. θ min = π 4 , θ max = 3 π 4 , ϕ min = 3 π 4 , ϕ max = 5 π 4 .

θ j = θ min + j Δ θ , j = − 1,0,1 , … , N θ + 2 , ϕ k = ϕ min + k Δ ϕ , k = − 1,0,1 , … , N ϕ + 2 ,

For the spatial discretization in the radial direction r, the mesh is generated by stacking layers along the radial direction. The radial grid nodes satisfy r ( 1 ) = 1 R s , r ( i + 1 ) = r ( i ) + Δ r ( i ) , i = 1 , … , N r . The radial grid spacing Δ r ( i ) is defined as

Δ r i = 0.01 R s , r i < 1.1 R s , min A × lg r i − 1 , Δ θ × r i − 1 , r i < 3.5 R s , A = 0.01 lg 1.09 . Δ θ × r i − 1 , r i > 3.5 R s .

On the spherical surface, the grid spacing is quasi-uniform, with the ratio between the minimum and maximum grid spacing being approximately 0.707. In the present simulation, the computational domain covers 1 R s ≤ r ≤ 25 R s, 0 ≤ θ ≤ π, and 0 ≤ ϕ ≤ 2 π. Each angular patch uses N θ = N ϕ = 42, and the radial direction is discretized using N r = 100.

Following the above procedure, the grid system required for the numerical simulation can be constructed. In addition, geometric information, such as cell-centered and face-centered coordinates, is required in the numerical calculations and is evaluated based on the trilinear transformation (Frey et al., 1978; Ivan et al., 2015; Feng et al., 2019).

2.3 Initial and boundary conditions

In the ideal MHD equations, the physical variables include the mass density ρ, temperature T (or gas pressure p), flow velocity u, and magnetic field B. The initial state of the solar wind plasma is obtained from the Parker (Parker, 1963) solar wind model. In this model, the solar-surface temperature T s = 1.3 × 1 0 6 K and density ρ s = 1.5 × 1 0 8 c m − 3 are used as input parameters to solve the corresponding integral equation. The radial velocity v r is then determined as a function of heliocentric distance. Once v r is obtained, the plasma density ρ is derived from the continuity equation, while the temperature T or gas pressure p is calculated using a polytropic equation of state with a constant γ = 1.05 to construct the initial Parker solar-wind solution. The initial tangential velocity components are set to v θ = 0 and v ϕ = 0, indicating that the initial solar wind flows purely in the radial direction.

The initial three-dimensional global magnetic field distribution within the computational domain is constructed using line-of-sight photospheric magnetic field observations provided by the Global Oscillation Network Group (GONG), combined with the potential field source surface (PFSS) model (Altschuler and Newkirk, 1969). The soalr wind simulation in this paper uses the GONG line-of-sight synoptic magnetogram product corresponding to CR 2203, available from the GONG data archive (https://gong.nso.edu/data/magmap/QR/mqs/201805/mrmqs180502/). The magnetogram is preprocessed before being used as the lower boundary of the PFSS model. Specifically, the native GONG synoptic magnetogram used in this study has a grid size of 360 × 180, with a longitudinal sampling interval of 1 ° and a uniform spacing of Δ ( sin ⁡ λ ) = 0.0111111 in the sine of the Carrington latitude, where λ denotes the Carrington latitude. The native GONG synoptic magnetogram is then transformed from the sine-of-latitude representation into a regular Carrington latitude–longitude map through coordinate transformation and linear interpolation. This preprocessing is performed solely to convert the native magnetogram into the coordinate representation required by the PFSS calculation. No additional Gaussian smoothing, spatial filtering, or artificial degradation is applied during this preprocessing. The initial coronal magnetic field is obtained from the PFSS model, in which the magnetic field is assumed to be current-free between the solar surface and the source surface, usually taken to be R ss = 2.5 R s (Riley et al., 2006; Parker, 1958). The corresponding magnetic scalar potential is obtained by solving Laplace’s equation, and the magnetic field is reconstructed using a spherical-harmonic expansion retaining harmonics up to degree 15. Beyond the source surface, the magnetic field is radially extrapolated outward under the assumption of radial magnetic-flux conservation, so that the radial magnetic-field strength decreases with heliocentric distance as r − 2, while B θ and B ϕ are set to zero.

The boundary conditions of the computational domain include both radial and angular boundaries. In the angular directions ( θ , ϕ ), due to the use of a cubed-sphere grid structure, physical quantities at the boundaries are communicated through interpolation between adjacent grid blocks, thereby ensuring the continuity of the numerical solution across the interfaces. In the radial direction, only the inner boundary at r = R s and the outer boundary are specified. At the inner boundary, located at the solar surface, the projected normal characteristic boundary condition is imposed (Nakagawa, 1980; Nakagawa et al., 1987). To ensure that the solar wind flows radially outward, whenever the computed radial velocity satisfies v r < 0, it is forcibly set to zero, while the density and gas pressure are held constant in time. At the outer boundary, an outflow condition is applied, whereby the boundary variables are extrapolated from the nearest interior grid cells, allowing the solar wind to leave the computational domain freely.

3 Spatial discretization of governing equations

3.1 Finite volume discretization

Equation 2 can be spatially discretized using the finite volume method. First, it is rewritten in divergence form as

∂ U ∂ t + ∇ ⋅ F = S .

The computational domain is then partitioned into a set of hexahedral control volumes. Integrating over a control volume V i , j , k yields

∂ ∂ t ∫ V i , j , k U d V + ∫ V i , j , k ∇ ⋅ F d V = ∫ V i , j , k S d V .

Applying the Gauss’s theorem to convert the volume integral into a surface integral leads towhere A denotes the boundary surface of the control volume.

∂ ∂ t ∫ V i , j , k U d V + ∮ A F ⋅ d A = ∫ V i , j , k S d V ,

To obtain a semi-discrete formulation, two approximations are introduced. First, the volume-averaged conservative variables are defined as U ̄ i , j , k = 1 V i , j , k ∫ V i , j , k U d V ≈ U i , j , k, where U i , j , k denotes the cell-centered value of the conservative variables in the control volume V i , j , k. Second, the face-averaged flux is approximated by 1 | A l | ∫ A l F ⋅ d A ≈ F l, where A l represents the l-th cell face. With these approximations, the governing equation becomeswhere S ̄ i , j , k = 1 V i , j , k ∫ V i , j , k S d V representing the volume-averaged source term over the control volume V i , j , k.

∂ U ̄ i , j , k ∂ t = − 1 V i , j , k ∑ l = 1 6 F l ⋅ A l + S ̄ i , j , k , (3)

For the l-th interface, let the outward unit normal vector be n l = ( n l , x , n l , y , n l , z ), and the interface area be | A l |. The flux term can then be expressed as

F l ⋅ A l = F l , x n l , x + F l , y n l , y + F l , z n l , z | A l | .

The evaluation of F l ⋅ A l can exploit the rotational invariance of the Euler equations (Yalim, 2008; Koop, 2008). Specifically, a transformation matrix T l is introduced to rotate the x-axis of the Cartesian coordinate system to align with the local normal direction n l. The velocity and magnetic field components are transformed from the ( x , y , z ) coordinate system to a local ( n , t 1 , t 2 ) coordinate system, where t 1 and t 2 are two orthogonal tangential directions forming a right-handed system with n l. Since the MHD equations retain the same form under orthogonal coordinate transformations, one obtains

T l F l ⋅ A l = T l F l , x n l , x + F l , y n l , y + F l , z n l , z | A l | = F l , x T l U | A l | ,

Using this property, the coordinate system is first rotated to the local frame to compute the normal flux F l , x, and then transformed back to the global coordinate system via T l − 1:

∑ l = 1 6 F l ⋅ A l = ∑ l = 1 6 T l − 1 T l F l ⋅ A l = ∑ l = 1 6 T l − 1 F l , x T l U | A l | . (8)

The transformation matrices are defined as

T = 1 0 0 0 0 0 0 0 0 0 n x n y n z 0 0 0 0 0 0 t 1 x t 1 y t 1 z 0 0 0 0 0 0 t 2 x t 2 y t 2 z 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 n x n y n z 0 0 0 0 0 0 t 1 x t 1 y t 1 z 0 0 0 0 0 0 t 2 x t 2 y t 2 z 0 0 0 0 0 0 0 0 0 1 T − 1 = 1 0 0 0 0 0 0 0 0 0 n x t 1 x t 2 x 0 0 0 0 0 0 n y t 1 y t 2 y 0 0 0 0 0 0 n z t 1 z t 2 z 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 n x t 1 x t 2 x 0 0 0 0 0 0 n y t 1 y t 2 y 0 0 0 0 0 0 n z t 1 z t 2 z 0 0 0 0 0 0 0 0 0 1 .

When n x 2 + n y 2 → 0, the tangential vectors are defined as

t 1 x = n z n x 2 + n z 2 , t 1 y = 0 , t 1 z = − n x n x 2 + n z 2 , t 2 x = − n x n y n x 2 + n z 2 , t 2 y = 1 n x 2 + n z 2 , t 2 z = − n y n z n x 2 + n z 2 .

Otherwise,

t 1 x = n y n x 2 + n y 2 , t 1 y = − n x n x 2 + n y 2 , t 1 z = 0 , t 2 x = n x n z n x 2 + n y 2 , t 2 y = n y n z n x 2 + n y 2 , t 2 z = − n x 2 + n y 2 .

The resulting semi-discrete finite-volume formulation is

∂ U ̄ i , j , k ∂ t = − 1 V i , j , k ∑ l = 1 6 T l − 1 F l , x T l U | A l | + S ̄ i , j , k . (4)

For the GLM-MHD system, the interface flux computation can adopt special treatment of the normal magnetic field component B x and the auxiliary variable ψ. According to the GLM formulation proposed by Dedner et al. (2002), the magnetic field and divergence-cleaning equations are given by

∂ B ∂ t + ∇ ⋅ u B − B u + ψ I = 0 , ∂ ψ ∂ t + ∇ ⋅ c h 2 B = − c h 2 c p 2 ψ .

In the local one-dimensional Riemann problem along the x-direction, the normal magnetic field component B x and the variable ψ form an independent 2 × 2 linear hyperbolic subsystem:

∂ B x ∂ t + ∂ ψ ∂ x = 0 , ∂ ψ ∂ t + c h 2 ∂ B x ∂ x = 0 .

Solving this subsystem yields the interface values

B x , l = B x , L + B x , R 2 − 1 2 c h ψ R − ψ L , ψ l = ψ L + ψ R 2 − c h 2 B x , R − B x , L .

Accordingly, the present work adopts a flux construction strategy that combines a seven-wave MHD solver with a decoupled 2 × 2 GLM subsystem. Specifically, the interface values of B x and ψ are first determined from the above linear subsystem, and their corresponding fluxes are computed as F B x , l = ψ l and F ψ , l = c h 2 B x , l. The numerical fluxes for the remaining conserved variables are then constructed using the seven-wave MHD solver, in which the required interface values of B x and ψ, together with their corresponding fluxes, are directly taken from the GLM subsystem without additional computation. This treatment reduces the contamination of the local physical solution by divergence errors.

After the interface numerical fluxes and source terms have been evaluated, the semi-discrete system given by Equation 9 is advanced in time using a second-order Runge-Kutta (Jameson et al., 1981; Gottlieb and Shu, 1998) scheme:where U ̄ i , j , k n and U ̄ i , j , k n + 1 denote the cell-averaged conservative variables at time levels t n and t n + 1, respectively, U ̄ i , j , k ( 1 ) is the intermediate-stage solution, and R ( U ̄ ) denotes the spatially discretized right-hand-side operator of the semi-discrete system Equation 3. The time-step size Δ t is determined according to the Courant-Friedrichs-Lewy (CFL) (Russell, 1989; Gottlieb and Tadmor, 1991) condition:where c h is taken as the maximum characteristic speed over the computational domain, and Δ h is the minimum spatial scale of the computational grid. In the present study, C F L = 0.5.

U ̄ i , j , k 1 = U ̄ i , j , k n + Δ t R U ̄ i , j , k n , U ̄ i , j , k n + 1 = 1 2 U ̄ i , j , k n + 1 2 U ̄ i , j , k 1 + Δ t R U ̄ i , j , k 1 ,

Δ t = C F L Δ h c h ,

3.2 Least-squares reconstruction

In the finite-volume discretization, the evaluation of interface fluxes requires the left and right states on the two sides of each cell interface, which are obtained through spatial reconstruction of the cell-centered variables. Since only volume-averaged quantities are stored during the computation, it is necessary to recover the spatial variation of variables using information from neighboring cells. In this work, a linear reconstruction based on the least-squares (Yalim, 2008; Mavriplis, 2003; Mandal and Rao, 2011) method is employed. Furthermore, the magnetic field gradients are corrected by incorporating the divergence-free constraint, providing more accurate and physically consistent inputs for the construction of interface states.

The least-squares method is a standard approach for obtaining approximate solutions to overdetermined systems, where the number of equations exceeds the number of unknowns. It establishes an optimal linear relationship between datasets. In the reconstruction procedure, information from the target cell and its 26 neighboring cells is utilized, and the conserved variables are first transformed into primitive variables.

The volume-averaged values of neighboring cells are locally expanded at the target cell using a Taylor expansion:where j = 1,2 , … , 26, and the expansion is applied component-wise to the primitive variable vector W, with W j and W i denoting the corresponding component values at the neighboring and current cell centers, respectively. The component-wise gradients ∂ W ∂ x, ∂ W ∂ y, and ∂ W ∂ z are to be determined. In practice, a weighted least-squares method is applied. A weighting factor ω j is introduced:

W j − W i = ∂ W ∂ x x j − x i + ∂ W ∂ y y j − y i + ∂ W ∂ z z j − z i ,

ω j = 1 x j − x i 2 + y j − y i 2 + z j − z i 2 .

The above system can be written in matrix form as A g = b:

A = ω 1 Δ x 1 ω 1 Δ y 1 ω 1 Δ z 1 ω 2 Δ x 2 ω 2 Δ y 2 ω 2 Δ z 2 ⋮ ⋮ ⋮ ω 26 Δ x 26 ω 26 Δ y 26 ω 26 Δ z 26 , g = ∂ W ∂ x ∂ W ∂ y ∂ W ∂ z i , b = ω 1 W 1 − W i ω 2 W 2 − W i ⋮ ω 26 W 26 − W i .

In the numerical implementation, this reconstruction is applied only to primitive variables except for the magnetic field.

In MHD simulations, the magnetic field must satisfy the divergence-free condition

∇ ⋅ B = ∂ B x ∂ x + ∂ B y ∂ y + ∂ B z ∂ z = 0 .

In the present framework, magnetic-field divergence errors are controlled through complementary treatments at the governing-equation and spatial-discretization levels. The GLM formulation described in Section 2.1 transports and damps divergence errors at the governing-equation level, whereas the divergence-constrained reconstruction introduced below suppresses their generation and accumulation during spatial reconstruction. This constraint is explicitly incorporated into the least-squares reconstruction system. The matrix form is given by A B g B = b B:

A B = ω 1 Δ x 1 ω 1 Δ y 1 ω 1 Δ z 1 0 0 0 0 0 0 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ω 26 Δ x 26 ω 26 Δ y 26 ω 26 Δ z 26 0 0 0 0 0 0 0 0 0 ω 1 Δ x 1 ω 1 Δ y 1 ω 1 Δ z 1 0 0 0 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 0 0 0 ω 26 Δ x 26 ω 26 Δ y 26 ω 26 Δ z 26 0 0 0 0 0 0 0 0 0 ω 1 Δ x 1 ω 1 Δ y 1 ω 1 Δ z 1 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 0 0 0 0 0 0 ω 26 Δ x 26 ω 26 Δ y 26 ω 26 Δ z 26 1 0 0 0 1 0 0 0 1 ,

g B = ∂ B x ∂ x ∂ B x ∂ y ∂ B x ∂ z ∂ B y ∂ x ∂ B y ∂ y ∂ B y ∂ z ∂ B z ∂ x ∂ B z ∂ y ∂ B z ∂ z T , b B = ω 1 B x , 1 − B x , i ⋮ ω 26 B x , 26 − B x , i ω 1 B y , 1 − B y , i ⋮ ω 26 B y , 26 − B y , i ω 1 B z , 1 − B z , i ⋮ ω 26 B z , 26 − B z , i 0 .

Although linear reconstruction achieves second-order accuracy in smooth regions, it may introduce nonphysical oscillations near discontinuities or strong gradients. Therefore, a limiter is applied. In this work, the Venkatakrishnan (Venkatakrishnan, 1993; 1995) limiter is used to control the reconstructed slopes and ensure stability and monotonicity.

For each component of the primitive variable vector W, the maximum and minimum values within the current cell and its 26 neighboring cells are defined as

W max = max W i , max j = 1 , … , 26 W j , W min = min W i , min j = 1 , … , 26 W j ,

and Δ + = W max − W i , Δ − = W min − W i. The reconstructed increment at Gaussian points on each face is

δ W f = W f − W i = ∂ W ∂ x x f − x i + ∂ W ∂ y y f − y i + ∂ W ∂ z z f − z i , f = 1 , … , 6 ,

where ( x f , y f , z f ) denotes the coordinates of the Gaussian point on the f-th face. The corresponding limiter is then defined aswith

φ f = ϕ Δ + δ W f , δ W f > 0 , ϕ Δ − δ W f , δ W f < 0 , 1 , δ W f = 0 ,

ϕ = Δ ± 2 + 2 Δ ± δ W f + ε Δ ± 2 + Δ ± δ W f + 2 δ W f 2 + ε .

Here, ε = ( K h ) 3, where h denotes the local characteristic grid size and K is a user-defined positive parameter. Following the recommendation of Koop (2008), a constant value of K = 0.3 is adopted throughout this work. The final limiter is taken as φ = min f = 1 , … , 6 φ f, and the reconstructed value is W f lim = W i + φ δ W f.

For the magnetic-field variables, the divergence-free condition should be effectively controlled during the numerical calculation. To this end, the gradient-correction method proposed by Hopkins (2015), Hopkins (2016) and Wang et al. (2019) is adopted. A divergence constraint is imposed on the reconstructed magnetic-field gradients to reduce the discrete magnetic-field divergence. Without altering the original least-squares reconstruction framework, this method controls the magnetic-field divergence by constructing local gradient-correction terms. The detailed derivation and associated parameter tests are provided in Wang et al. (2019).

According to the parameter tests of Wang et al. (2019), the divergence error shows little further reduction beyond 60 correction iterations. Therefore, balancing divergence control and computational cost, we set k = 60 in the present simulations. Through the above treatments, the magnetic-field divergence error is maintained at a low level, thereby improving the stability and physical consistency of the numerical computation.

3.3 Roe scheme

As discussed above, in the GLM-MHD system, the normal magnetic-field component B x and the auxiliary variable ψ are solved independently through the 2 × 2 GLM subsystem, yielding the corresponding interface states. The remaining conserved variables therefore constitute a seven-wave MHD subsystem, for which the interface numerical flux can be evaluated using a standard MHD Riemann solver.

In Equation 4, F l , x ( T l U ) denotes the numerical approximation to the normal flux across the l-th interface in the local coordinate system. In practice, this flux is evaluated from the reconstructed left and right states at the interface. Specifically, after spatial reconstruction, the left and right states in the global coordinate system, denoted by U L ( g ) and U R ( g ), are projected onto the local normal coordinate system as

U L = T l U L g , U R = T l U R g .

The local numerical flux is then obtained by solving a one-dimensional approximate Riemann problem in the normal direction, namely,

F l , x T l U ≈ F ̂ l , x U L , U R .

The resulting flux is transformed back to the global coordinate system through T l − 1 and accumulated in the finite-volume discretization. This procedure exploits the rotational invariance of the governing equations and reduces the multidimensional interface-flux evaluation to a set of local one-dimensional problems.

Most approximate Riemann solvers can be written in the following unified form:

F ̂ U L , U R = 1 2 F U L + F U R − 1 2 D U L , U R U R − U L ,

where the first term is the central average of the physical fluxes, and the second term represents the numerical dissipation. The main difference among various approximate Riemann solvers lies in the construction of the dissipation matrix D.

The basic idea of the Roe solver (Roe, 1981) is to approximate the original nonlinear Riemann problem by a locally linearized hyperbolic system. In order for the linearized system to accurately represent the wave-propagation properties of the original nonlinear system near the interface, the Roe-averaged matrix must satisfy several conditions, including consistency. The Roe flux can be written aswhere

F ̂ U L , U R = 1 2 F U L + F U R − | A ̂ | U R − U L ,

The matrix A ̂ is the Roe-averaged Jacobian constructed from the left and right states U L and U R. The diagonal matrix Λ = d i a g ( λ k ) contains the eigenvalues of A ̂, while R and R − 1 denote the matrices of right and left eigenvectors in conservative variables, respectively.

The right and left eigenvectors in conservative variables, denoted by R k and L k, are related to the corresponding eigenvectors in primitive variables, r k and l k, through

L k = l k ∂ W ∂ U , R k = ∂ U ∂ W r k .

The jump in conservative variables across the interface is defined as

Δ U = U R − U L .

It can be expanded in characteristic space as

Δ U = ∑ k = 1 7 η k R k ,

where R k is the k-th right eigenvector and

η k = L k ⋅ U R − U L

is the corresponding wave strength.

According to the wave-propagation properties of hyperbolic systems of conservation laws, each characteristic wave propagates along its corresponding characteristic direction with speed λ k. In order to distinguish wave propagation directions in an upwind manner, the absolute values of the eigenvalues, | λ k |, are introduced into the numerical dissipation term. The Roe numerical flux at the interface is therefore given by

F ̂ U L , U R = 1 2 F U L + F U R − 1 2 R | Λ | R − 1 U R − U L = 1 2 F U L + F U R − 1 2 ∑ k = 1 7 | λ k | η k R k .

3.4 Numerical path-preserving Roe scheme

A correction method based on the NPP approach was proposed by Xu et al. (2023), Xu et al. (2024), Xu et al.(2025). The method is designed to align the numerical path with the characteristic integral curve associated with the Alfvén wave, thereby preserving the magnitude of the tangential magnetic field. After correction, the Alfvén wave characteristics can be better captured, substantially reducing the pseudo-convergence phenomenon while maintaining the conservation properties of the original scheme and the integrity of the characteristic system.

For completeness, the seven eigenvalues of the MHD system, ordered from smallest to largest, are as follows: left fast magnetosonic wave, left Alfvén wave, left slow magnetosonic wave, contact discontinuity, right slow magnetosonic wave, right Alfvén wave, and right fast magnetosonic wave, i.e.,

u x − c f x , u x − c a x , u x − c s x , u x , u x + c s x , u x + c a x , u x + c f x .

Since the correction to the characteristic system only involves the Alfvén wave, only the characteristic vectors, eigenvalues, and wave strengths for the Alfvén wave are provided here. For the characteristic values, eigenvectors, and wave strengths of other waves, please refer to the work of Cargo and Gallice (1997).

In the Roe average matrix, the eigenvalue, right eigenvector, and wave strength corresponding to the Alfvén field are as follows:

where:

λ u ̄ x ± c ̄ a x = u ̄ x ± c ̄ a x ,

R u ̄ x ± c ̄ a x = 0 , 0 , ± ρ ̲ β z , ∓ ρ ̲ β y , ± ρ ̲ u ̄ y β z − u ̄ z β y , − S ρ ̲ β z , S ρ ̲ β y T ,

η u ̄ x ± c ̄ a x = L u ̄ x ± c ̄ a x ⋅ Δ U = 1 2 ∓ β y Δ u z ± β z Δ u y + S ρ ̲ β y Δ B z 1 − β z Δ B y 1

c ̄ a x 2 = B ̲ x 2 ρ ̲ , S = s g n B x , β y , z = B ̲ y , z | B ̲ t | , | B ̲ t | = B ̲ y 2 + B ̲ z 2 .

For the ideal MHD equations, the Roe average is commonly used to construct the intermediate state, where related variables can be given by the following weighted average form:

ξ ̄ = ρ L ξ L + ρ R ξ R ρ L + ρ R , ξ ̲ = ρ L ξ R + ρ R ξ L ρ L + ρ R .

3.4.1 Tangential magnetic field correction

Since the Roe averaging process for the magnetic field involves averaging the components of B t separately, this leads to the traditional Roe average placing the magnetic field on the straight line connecting the left and right states in the tangential plane. This does not satisfy the requirement for maintaining the constancy of the tangential magnetic field magnitude in Alfvén waves. Therefore, the Roe averaging process for the magnetic field is modified as follows:

B y * ̲ = B y ̲ + | B t | * ⁡ cos θ L + θ R 2 , B z * ̲ = B z ̲ + | B t | * ⁡ sin θ L + θ R 2 .

θ L , R = arctan B z B y L , R , θ L , R ∈ − π , π .

δ Δ α x = 1 , | x | ≤ Δ α , 0 , | x | > Δ α .

The additional multiplier δ Δ α π − | θ l − θ r | is introduced to avoid applying Roe averaging to the other characteristic fields, for which the direction of the tangential magnetic field remains unchanged. Δ α can be selected as approximately 0.6 radians based on numerical experiments.

3.4.2 Characteristic system correction

In the traditional Roe method, the evolution of the tangential magnetic field depends on the CFL condition, and its numerical average typically lies on the straight line connecting the left and right states after time stepping. The magnitude of the tangential magnetic field remains unchanged only when CFL = 1. This treatment does not guarantee that the solution strictly satisfies the geometric constraints of the Alfvén wave under general CFL conditions.

We then correct the characteristic speeds, wave strengths, and eigenvector directions of the Alfvén wave characteristic system using the NPP method, ensuring that the grid-averaged state after time-stepping always lies on the Alfvén integration curve, i.e., the circular arc trajectory in the tangential plane, thereby achieving precise preservation of the physical propagation path.

First, the correction to the wave strength is straightforward and simply involves scaling it proportionally to the chord length. Using basic geometric relations, the corrected wave strength is given by:

η u x ± c a x * = η u x ± C a x sin C F L Δ θ 2 sin Δ θ 2

Next, we adjust the characteristic speed so that it propagates the numerical solution onto the target integration path under any CFL condition:

λ u x ± c a x , N P P = sgn λ u x ± c a x Δ x Δ t , sgn x = 1 , x ≥ 0 − 1 , x < 0 (4)

For the nonuniform six-patch spherical grid, Equation (Equation 4) is applied locally at each numerical face l usingwhere x L and x R are the centroids of the adjacent control volumes, x f , l is the face center, and n l is the unit face-normal vector. The corrected Alfvén eigenvalue is then evaluated aswhere Δ t is the synchronized time step. At patch interfaces, the geometric quantities are evaluated in a common physical Cartesian coordinate system. Finally, the correction to the eigenvector is the most critical and intricate part of the NPP-Roe correction. By analyzing the original eigenvector, we find that the direction of the tangential magnetic field is primarily controlled by the parameters β y and β z, which describe the direction of the tangential magnetic field in the plane. In the Roe average, the direction corresponds to the traditional Roe point. However, to meet the NPP requirement, this direction must be rotated to the corresponding NPP-Roe point direction. This rotation is achieved by introducing an angle θ NPP, whose value depends on the wave propagation direction and the angle between the left and right states, and is also related to the CFL number. By doing so, the corrected β y * and β z * will ensure that the tangential magnetic field evolves along the correct circular arc path:

Δ x l = x f , l − x L ⋅ n l + x R − x f , l ⋅ n l ,

λ u x ± c a x , l NPP = sgn λ u x ± c a x , l Δ x l Δ t ,

β y * = cos θ NPP , β z * = sin θ NPP

θ NPP = θ L − C F L θ R − θ L 2 , λ u x ± c a x < 0 θ R − C F L θ R − θ L 2 , λ u x ± c a x > 0 | θ R − θ L | ≤ π

Although the above corrections guarantee that the density, normal momentum, and tangential magnetic field satisfy the jump conditions, the rotational discontinuity also requires that the thermodynamic pressure remains unchanged, a condition that is not automatically satisfied in the energy components. Therefore, we need to apply an additional correction to the energy terms in the eigenvector to maintain the consistency between the corrected eigenvector and the Alfvén field. The enforcement of this constraint yields the corrected energy component as

r E 1 * = ± ρ ̲ β z * u y − β y * u z + S ρ ̲ β y * B z − β z * B y .

The notation { ⋅ } denotes the arithmetic average of the left and right states, i.e., ξ = ξ L + ξ R 2. Consequently, the fully corrected Alfvén eigenvector is expressed as

0 , 0 , ± ρ ̲ β z * , ∓ ρ ̲ β y * , ± ρ ̲ β z * u y − β y * u z + S ρ ̲ β y * B z − β z * B y , − S ρ ̲ β z * , S ρ ̲ β y * .

The corrected Alfvén eigenvector not only satisfies the NPP requirement in terms of direction, but also meets the Rankine-Hugoniot conditions for the Alfvén rotational discontinuity in terms of the energy components, ensuring both geometric and physical consistency of the corrected eigenstructure.

It should be noted that although the notation { ⋅ } still represents the arithmetic average, the left and right states corresponding to this average are not the initial left and right states U L and U R of the entire Riemann problem, but rather the local intermediate states on either side of the Alfvén wave, i.e., the two physical states on either side of this Alfvén wave. Therefore, this average should be calculated across the local jump corresponding to the Alfvén wave, not directly from the initial left and right states at the interface.

At this point, all NPP corrections have been completed. Using the corrected eigenvalues and eigenvectors, numerical fluxes at the grid interfaces can be computed based on the Roe scheme, allowing the time integration of all physical quantities.

4 Numerical results

Background solar wind simulations are primarily used to characterize large-scale flow and magnetic field structures. The resulting large-scale solar wind solution reflects the coupled effects of the magnetic-field topology, boundary conditions, source-term treatment, grid resolution, and interactions among different MHD wave families. Consequently, the background solar wind distribution alone is insufficient to directly evaluate the ability of a numerical scheme to capture Alfvén waves. Based on this consideration, this chapter first adopts a one-dimensional MHD Alfvén rotational discontinuity as a benchmark test case. Numerical simulations are performed using both the classical Roe scheme and the NPP-Roe scheme, which incorporates a numerical path-preserving modification, allowing a direct comparison of their capability in capturing Alfvén wave propagation under controlled conditions.

Following the one-dimensional benchmark test, we further employ GONG observational data to simulate the background solar wind for CR 2203, assessing the applicability of the NPP-Roe scheme in modeling complex coronal-solar wind structures and examining whether the proposed scheme can be incorporated into a realistic global MHD framework without introducing evident distortion into the large-scale background solar wind solution.

4.1 Numerical validation of Alfvén-wave capturing capability

To quantitatively assess the ability of the classical Roe and NPP-Roe schemes to capture Alfvén wave propagation, we construct a one-dimensional MHD finite-angle Alfvén rotational discontinuity as a benchmark problem. In this example, the density, normal velocity, gas pressure, and normal magnetic field are identical on both sides of the discontinuity, while only the tangential magnetic field direction rotates by a finite angle in the transverse plane. The corresponding tangential velocity jump is determined by the Rankine-Hugoniot relation for Alfvén waves. This setup ensures that the initial discontinuity corresponds to a single right-propagating Alfvén wave, minimizing the influence of fast and slow magnetosonic waves and contact discontinuities. Based on this construction, the two schemes can be directly compared in terms of Alfvén wave speed capture, preservation of jump conditions, and maintenance of the tangential magnetic field magnitude. The specific configuration is as follows.

Let the primitive variables be W = ( ρ , u , v , w , p , B y , B z ), where u is the normal velocity, v t = ( v , w ) the tangential velocity, and B t = ( B y , B z ) the tangential magnetic field components. The initial discontinuity is placed at x = 0. The primitive variables are initialized as

with the normal magnetic-field component and the ratio of specific heats set to B x = 0.75 and γ = 1.4, respectively. From this configuration, the tangential magnetic field magnitude on both sides is unity, with a relative rotation angle of 120 °. The tangential velocity jump satisfies the right-propagating Alfvén rotational discontinuity relation, so that the initial discontinuity corresponds to a pure right-propagating Alfvén rotational discontinuity. The theoretical propagation speed is s = u + B x ρ = 0.75 . Hence, the exact solution consists of a single propagating discontinuity located at x / t = 0.75. The computational domain is x ∈ [ − 1.5 , 1.5 ], with a final time t end = 0.8 and a CFL number of 0.5. The theoretical wave front is at x = 0.6, sufficiently far from the domain boundaries to avoid boundary effects on the solution.

W x , 0 = 1,0,0,0,0.5 , 1,0 , x < 0 , 1,0 , 3 2 , − 3 2 , 0.5 , − 1 2 , 3 2 , x > 0 ,

We first adopt a uniform grid with N = 1000 to compare the two schemes under relatively coarse resolution, focusing on the profiles of physical variables, tangential magnetic field phase-plane trajectories, and key error metrics. Subsequently, six grids with N = 500,1000,2000,4000,8000,16000 are employed to examine the error behavior under grid refinement.

Figure 2 presents the numerical results for the one-dimensional Alfvén wave problem with N = 1000 using the classical Roe scheme and the NPP-Roe scheme. For this test case, the density and gas pressure are theoretically expected to remain constant during the propagation of the discontinuity, with ρ = 1 and p = 0.5. Figures 2a,b show that the classical Roe scheme produces noticeable disturbances in density and pressure near the wave front, indicating spurious non-physical perturbations introduced by the scheme and deviations in the characteristic decomposition and propagation of the Alfvén wave. In contrast, the NPP-Roe scheme maintains the density and pressure close to their theoretical values, demonstrating its effectiveness in suppressing spurious perturbations in non-Alfvén wave variables.

Figure 2c depicts the distribution of the tangential magnetic field component B y. The classical Roe scheme exhibits a relatively wide numerical transition region and significant numerical diffusion near the wave front, whereas the NPP-Roe scheme captures the jump structure of B y more sharply, with wave front positions and left/right states closely matching the exact solution. Furthermore, Figure 2d shows the magnitude of the tangential magnetic field | B t |. For an Alfvén wave, the tangential magnetic field primarily undergoes directional rotation, while its magnitude should remain constant, | B t | = 1. The classical Roe scheme causes | B t | to deviate noticeably from the theoretical value within the transition region, whereas the NPP-Roe scheme preserves the magnitude more accurately.

These results indicate that the NPP modification significantly improves the Roe scheme’s capability to capture Alfvén waves, particularly in suppressing non-physical perturbations and maintaining the magnitude of the tangential magnetic field.

To further assess the ability of the two schemes to preserve the tangential magnetic field magnitude, Figure 3 shows the numerical trajectories in the ( B y , B z ) plane. This phase-plane representation illustrates the evolution of the tangential magnetic field within the numerical transition region. For this test case, the exact left and right states lie on the reference circle | B t | = 1, with the tangential field direction rotating from the left to the right state. Therefore, the proximity of the numerical trajectories to the reference circle directly reflects whether each scheme maintains the Alfvén wave’s physical characteristic of changing the field direction without altering its magnitude during propagation.

As shown in Figure 3, the classical Roe scheme’s trajectory shifts noticeably inward from the reference circle, indicating a reduction in | B t | within the transition region. This observation is consistent with the concave dip of | B t | near the wave front in Figure 2d. In contrast, the NPP-Roe scheme’s trajectory closely follows the reference circle and effectively connects the left and right exact states, demonstrating better physical consistency in preserving the rotation path of the tangential magnetic field. The results of Figure 3 corroborate those in Figure 2, further confirming that the NPP-Roe scheme more accurately captures the propagation characteristics of the Alfvén wave.

To quantitatively compare the performance of the two schemes, Table 1 lists key error metrics for N = 1000. For this Alfvén wave test case, the exact solution satisfies ρ = 1, p = 0.5, and | B t | = B y 2 + B z 2 = 1. Accordingly, the maximum errors in density and pressure are defined as

δ ρ ∞ = max i | ρ i − 1 | , δ p ∞ = max i | p i − 0.5 | .

TABLE 1. Quantitative error metrics for the one-dimensional Alfvén rotational discontinuity at N = 1000. Values reported as 0 are below machine precision. Consequently, the corresponding Roe/NPP-Roe ratios are reported as ∞.

The tangential magnetic field magnitude errors are expressed using both the maximum and mean errors:

δ | B t | ∞ = max i B t , i | − 1 , E | B t | 1 = 1 N ∑ i B t , i | − 1 .

To examine whether the numerical solution satisfies the right-propagating Alfvén wave relation between velocity and magnetic field, we define the local residual asand the corresponding metrics

R A , i = v t , i + B t , i ρ i − 1,0 2 ,

R A ∞ = max i R A , i , R A 1 = 1 N ∑ i R A , i .

Additionally, the wave front position error is defined by E s = | ξ front − 0.75 | , where ξ front is determined from the midpoint of B y, and the theoretical wave front is located at ξ = 0.75.

Table 1 shows that the classical Roe scheme introduces noticeable errors in density and pressure, whereas the NPP-Roe scheme maintains these variables close to their theoretical values at the current output precision. This observation is consistent with Figures 2a,b, indicating that the NPP modification effectively suppresses spurious perturbations in non-Alfvén wave variables. For the tangential magnetic field magnitude, the classical Roe scheme yields a maximum and mean error of 2.111 × 1 0 − 1 and 1.882 × 1 0 − 2, respectively, while the NPP-Roe scheme reduces these errors to the order of 1 0 − 7. Similarly, the Alfvén relation residual for the NPP-Roe scheme is roughly 1 0 4 times smaller than that of the classical Roe scheme. These results further demonstrate that the NPP-Roe scheme not only improves the preservation of | B t | but also better satisfies the Alfvén velocity-magnetic field relation. In contrast, the wave front position errors for both schemes are similar and on the order of 1 0 − 2, indicating that at the current resolution, the wave front position does not fully reflect the advantages of the NPP modification. The improvements of NPP-Roe are primarily in maintaining the physical constraints of the Alfvén wave, i.e., significantly reducing the tangential field magnitude error and the Alfvén relation residual.

Based on the single-grid error comparison in Table 1, Figure 4 further shows the error trends of both schemes under different grid resolutions. Here, E | B t | 1 quantifies the ability to preserve the tangential field magnitude, while R A 1 measures the degree to which the numerical solution satisfies the Alfvén velocity-magnetic field relation. Figure 4 indicates that, as the grid resolution increases, both error metrics of the classical Roe scheme decrease, showing that grid refinement can partially reduce its errors. Nevertheless, at all resolutions, the Roe scheme’s errors remain significantly higher than those of the NPP-Roe scheme. In contrast, E | B t | 1 for the NPP-Roe scheme consistently stays on the order of 1 0 − 7, and R A 1 remains around 1 0 − 6, demonstrating that the scheme reliably preserves the main physical constraints of the Alfvén wave across grid resolutions. Slight fluctuations in R A 1 for NPP-Roe at certain grids are observed, but their magnitude remains far below that of the Roe scheme. These minor variations likely arise from the wave front’s position relative to grid points and output precision, and do not affect the overall comparison between the two schemes.

Overall, Figures 2–4 and Table 1 show that, in the present one-dimensional benchmark, the NPP-Roe scheme improves waveform resolution, better preserves the tangential magnetic-field magnitude, and more closely satisfies the Alfvén velocity–magnetic-field relation under both single-grid and multiple-grid conditions.

4.2 NPP-Roe scheme under background solar wind conditions

In this subsection, the NPP-Roe scheme is employed to simulate the background solar wind for CR 2203. The solution is regarded as quasi-steady whenwhere N denotes the total number of computational cells and U i n is the solution vector in cell i at time level n. Otherwise, the time-marching procedure is continued. All simulation results presented below for CR 2203 (14 April–11 May 2018) correspond to the relaxed quasi-steady solution. The calculation is carried out on the supercomputer at the College of Engineering Physics, Shenzhen Technology University. According to the computational configuration, 144 MPI processes are used in the present simulation. Under this configuration, the simulation achieves a quasi-steady solution in approximately 6 h of wall-clock time.

1 N ∑ i = 1 N U i n − U i n − 1 U i n − 1 − ∑ i = 1 N U i n − 1 − U i n − 2 U i n − 2 ≤ ϵ , ϵ = 1 0 − 7 ,

The observational comparisons presented in this section are qualitative and are intended to assess only the consistency of the locations, spatial distributions, and large-scale morphology of the modeled structures. The SDO/AIA 193 Å images are used as observational proxies for the large-scale distributions of coronal holes and bright closed-field regions, whereas the SOHO/LASCO C2 images are used to examine the morphology of streamer-related coronal density structures. No automated boundary extraction, image-correlation analysis, radial-profile comparison, or other quantitative image-validation metric is performed in the present study. Moreover, the synthetic pB image shown in Figure 5 is processed using the Normalizing Radial Graded Filter (NRGF) and is compared with the LASCO C2 image only in terms of the approximate locations, angular extents, extension directions, and overall morphology of the principal brightness-enhanced structures. Therefore, this comparison should be interpreted as a qualitative morphological assessment rather than as a radiometric or pixel-by-pixel validation. A quantitative comparison based on co-registered and consistently processed observational and synthetic images is left for future work.

Figure 5 presents the simulated global coronal magnetic field at the solar surface, r = 1 R s , and compares it with the SDO/AIA 193 Å extreme-ultraviolet (EUV) observations. Figure 5a shows the synoptic map of the radial magnetic field B r obtained using the NPP-Roe method. At low latitudes, the magnetic field exhibits multiple structures with alternating positive and negative polarities, reflecting the concentration and complexity of active-region magnetic flux in the low-latitude regions. In contrast, the magnetic-field distribution at high latitudes is relatively smooth and characterized by large-scale unipolar structures, indicating that the polar magnetic field is more coherent on global scales. The local magnetic-field strength at low latitudes can reach approximately ± 40 G, whereas the field strength at high latitudes is comparatively weaker. This distribution pattern is consistent with observations of the global solar magnetic field and photospheric flux-transport maps, in which active-region magnetic flux is mainly concentrated at low and middle latitudes and enhances the local magnetic complexity, while the high-latitude regions are primarily controlled by the large-scale polar magnetic field (Petrie, 2024).

Figure 5b shows the simulated open-closed magnetic-field topology. The closed magnetic fields, shown in white, are mainly concentrated at low and middle latitudes, with both ends of the field lines anchored at the solar surface, forming closed magnetic loops. The open magnetic fields, shown in black, are primarily distributed at high latitudes and in fragmented regions near the equator. These field lines extend into the heliosphere and correspond to the open magnetic flux in interplanetary space, representing important source regions of the fast solar wind. The local fragments of open magnetic field reflect the complex distribution of open magnetic flux at equatorial and low-to-middle latitudes. This result is consistent with the current theoretical understanding that coronal holes generally correspond to open-field regions, where open magnetic flux extends from the corona into the heliosphere and is closely related to solar-wind source regions (Antonucci et al., 2020; Linker et al., 2021).

Figure 5c presents the SDO/AIA 193 Å EUV Carrington synoptic map. Bright regions generally correspond to closed-field regions with relatively high plasma density, whereas dark regions correspond to open-field regions, namely, coronal holes. The red curves denote the open-closed magnetic-field boundaries derived from the simulation, which show overall agreement with the bright and dark structures observed in EUV emission. Local discrepancies may arise from uncertainties in the input magnetic-field data, uncertainties in the magnetic boundary conditions, and temporal differences between the observations and the simulations. This result is consistent with existing validation approaches for coronal magnetic-field models, in which EUV images of coronal holes, closed-loop systems, and bright-dark emission structures are used to assess the reliability of the simulated open-closed magnetic-field boundaries (Wagner et al., 2022; Parenti et al., 2022).

Overall, the simulation broadly reproduces the principal global topological structures of the solar magnetic field and shows general consistency with the EUV observations, supporting the applicability of the NPP-Roe method to large-scale coronal magnetic-field modeling.

In studies of the large-scale solar coronal structure, the white-light pB serves as a crucial observational parameter for probing the spatial distribution of electron density, as the pB primarily arises from Thomson scattering of sunlight by free electrons, with intensity varying according to electron density and line-of-sight integration (Quémerais and Lamy, 2002). Figure 6 presents the pB observed along the eastern and western limbs of the Sun by the LASCOC2 instrument onboard the SOHO spacecraft. The pB along the eastern and western coronal edges exhibits alternating bright and dark regions. Bright regions typically correspond to high-density structures along the line of sight, manifesting as horizontally elongated, streamer-like bands in low-latitude areas, whereas some dark regions may result from data gaps, detector defects, or image processing masks. Overall, the enhanced pB regions are broadly distributed along the model-calculated 2.5 R s radial magnetic neutral line ( B r = 0 ), indicating a close association between high-density coronal regions and the magnetic neutral line. Moreover, differences in continuity and coherence exist between the bright regions on the eastern and western sides, with some bright structures slightly offset from the neutral line. High-density structures offset from the neutral line are generally associated with non-polar streamers (Riley et al., 2012b).

Figure 7 shows the distributions of radial solar wind speed V r and plasma density N at 2.5 R s, as predicted by the NPP-Roe MHD model. The model indicates that high-density, low-speed flows correspond well to pB regions shown in Figure 6, whereas low-density, high-speed flows appear outside the neutral line (Lamy et al., 2019). Notably, the plasma density N distribution in Figure 7 also reveals high-density structures in mid-to low-latitude regions away from the magnetic neutral line, which are consistently identified in both the MHD simulation and observations (Decraemer et al., 2019). These results suggest that coronal density is not confined solely to the vicinity of the neutral line and exhibits complex asymmetries. Furthermore, they indicate that the model broadly captures the global streamer-belt distribution and the main features of the coronal electron-density structure.

The left panel of Figure 8 shows the synthetic white-light pB image generated from the numerical simulation based on the NPP-Roe method, whereas the right panel presents the corresponding SOHO/LASCO C2 white-light coronal observation. The white-light pB distribution can reflect the electron density structures along the line of sight (Liang et al., 2023). In this study, the synthetic brightness is obtained by integrating the coronal plasma electron density along the line of sight over the radial range of 2.3-6 R s, which enables a morphological comparison with the outer-coronal white-light structures observed by LASCO C2. During the image-processing procedure for the simulated data, the NRGF method was applied. Because the white-light coronal brightness decreases rapidly with increasing heliocentric distance, radial normalization helps suppress the strong radial brightness gradient and thereby enhances angular structures and localized brightness enhancements in the outer corona (Morgan et al., 2006).

As shown in Figure 8, both the synthetic image and the LASCO C2 observation contain bilateral brightness-enhanced structures extending outward from the occulting disk. Their approximate locations, extension directions, and overall large-scale morphologies show qualitative agreement. Because the comparison is based on morphology after image enhancement rather than on a common radiometric and pixel-wise comparison framework, no conclusion regarding quantitative brightness agreement is drawn here.

Nevertheless, some discrepancies remain between the simulation and the observation. Compared with the SOHO/LASCO C2 image, the brightness structures in the synthetic white-light pB image appear generally smoother, as reflected by the insufficient sharpness and clarity of small-scale coronal streamers and their boundaries. Some narrow streamer structures and small-scale brightness variations visible in the observation are not fully reproduced in the simulation. In addition, the angular inclination and streamer-like morphology of the simulated brightness-enhanced regions appear somewhat broader. These differences suggest that, although the model captures the overall large-scale coronal structure, its ability to capture fine local features remains limited by the current spatial resolution and simplified thermodynamic treatment. In practice, the most immediate improvements would be to increase the radial and angular grid resolution, particularly in the low corona and near streamer cusps and open–closed field boundaries. The energy equation will also be extended to include field-aligned thermal conduction, while the present volumetric heating term will be refined to depend on local magnetic properties, such as the field strength or flux-tube expansion. Adaptive mesh refinement (AMR) and time-dependent magnetogram assimilation are not included in the current model and are considered longer-term developments.

Figure 9 presents the latitudinal-longitudinal distributions of the radial velocity V r, number density N, and temperature T of the background solar wind for CR 2203 at heliocentric distances of 2.5, 5, and 20 R s. The results reveal a clear and persistent large-scale structure at all three distances: low latitudes are characterized by slower, denser, and cooler plasma, whereas high latitudes exhibit faster, more rarefied, and relatively hotter solar wind, consistent with the meridional-plane distributions shown in Figure 10. As the heliocentric distance increases, the solar wind undergoes continuous acceleration, accompanied by an overall decrease in density and a progressive smoothing of the temperature distribution, indicating that the background solar wind already exhibits a well-defined radial evolution in the near-Sun region.

In addition, the low-latitude slow and dense structures gradually narrow during outward propagation, while the transition regions between different flow regimes become sharper and exhibit enhanced gradients, with their overall spatial locations remaining largely unchanged. These results suggest that once established in the near-Sun region, the large-scale structure of the background solar wind is largely preserved during radial expansion. Its evolution is primarily manifested as radial variations in the amplitudes of physical quantities rather than a substantial reconfiguration of the spatial structure. This behavior is consistent with recent radial alignment observations from Parker Solar Probe and Solar Orbiter, which show that the same large-scale density structures can be continuously identified during propagation while evolving predominantly in the radial direction (Berriot et al., 2025; Raouafi et al., 2023).

Figure 10 presents the simulated distributions of solar wind parameters during CR 2203 on two representative meridional planes, illustrating the spatial structures of the radial solar wind speed V r and the plasma number density N within the range from 1 to 20 R s. The left column shows the radial solar wind speed V r (km s − 1) on the ϕ = 180 ° − 0 ° and ϕ = 270 ° − 90 ° meridional planes, respectively, while the right column displays the corresponding plasma number density N ( log 10 c m − 3 ). The black curves and arrows denote magnetic field lines and their directions, respectively.

The velocity distributions reveal a pronounced latitudinal stratification of the solar wind. High-latitude regions are predominantly characterized by high-speed outflows, where the magnetic field lines are mainly open and extend radially outward. Open magnetic field regions facilitate the outward release of plasma and the formation of high-speed solar wind streams, consistent with observational and theoretical understanding of coronal holes as the primary sources of fast solar wind (Cranmer, 2009). In contrast, the solar wind speed is significantly lower at low latitudes, particularly near the equatorial region. Such slow-wind structures are typically associated with streamer-dominated flow configurations and the magnetic topology near their cusps, where closed magnetic fields exert stronger confinement on the plasma (Pneuman and Kopp, 1971). From a global perspective, the configuration featuring fast wind at high latitudes, slow wind at low latitudes, and the association of the heliospheric current sheet with the streamer belt represents a typical large-scale structure of the heliospheric magnetic field near solar minimum (Owens and Forsyth, 2013).

The distribution of plasma number density exhibits a clear anti-correlation with the solar wind speed. Higher densities are primarily concentrated in low-latitude streamer regions and in the vicinity of the heliospheric current sheet, forming a band-like enhancement along the equatorial direction. The slow solar wind is generally associated with the streamer belt and regions adjacent to the heliospheric current sheet, and is characterized by relatively higher plasma densities (Abbo et al., 2016; Owens and Forsyth, 2013). In contrast, the high-latitude fast solar wind regions exhibit significantly lower densities, consistent with the low-density, open-field-dominated nature of coronal hole source regions (Cranmer, 2009). Moreover, three-dimensional observations from Ulysses indicate that high-latitude regions are dominated by fast, tenuous, and relatively uniform solar wind, whereas low-latitude regions exhibit slower, more structured, and denser solar wind (McComas et al., 2003).

Figure 11 presents the magnetic-field divergence results after 60 magnetic-field correction iterations, including the spatial distribution of the local divergence and the temporal evolution of the global divergence error. Figure 11a shows the distribution of the divergence of the perturbed magnetic field, B 1, on the ϕ = 270 °-90 ° meridional plane. For each control volume, the local divergence is evaluated using Gauss’s theorem in the finite-volume form, ( ∇ ⋅ B 1 ) I = 1 V I ∮ ∂ V I B 1 ⋅ d A, where V I is the volume of the I-th control volume and ∂ V I denotes its boundary surface. The surface integral is evaluated numerically using Gaussian quadrature over the six faces of each hexahedral control volume. As shown in Figure 11a, the magnitude of ∇ ⋅ B 1 remains on the order of 1 0 − 5 over most of the computational domain. Relatively large values are confined mainly to the vicinity of the inner boundary and regions with strong magnetic-field gradients.

Figure 11b shows the time evolution of the volume-weighted global divergence error of the total magnetic field, B, defined as E r r o r ( ∇ ⋅ B ) = ∑ I = 1 M ( ∇ ⋅ B ) I V I ∑ I = 1 M V I, where M is the total number of control volumes in the computational domain. The normalization by the total computational volume, ∑ I = 1 M V I, makes this quantity the volume-weighted mean of the absolute local divergence over the entire computational domain. The plotted quantity in Figure 11b is log 10 E r r o r ( ∇ ⋅ B ) . Throughout the simulation, this value remains close to − 5.92, corresponding to an average divergence error of approximately 1.2 × 1 0 − 6. Although a slight increase is observed as time advances, the growth remains very limited and does not span even one order of magnitude. These results show that, after 60 correction iterations, the magnetic-field divergence error is maintained at a low level over the entire simulation, supporting the effectiveness and numerical stability of the adopted divergence-control procedure.

5 Discussion and conclusion

For the one-dimensional finite-angle Alfvén rotational-discontinuity benchmark considered in this study, the NPP-Roe scheme modifies the propagation direction, eigenvalues, eigenvectors, and wave strengths associated with the Alfvén characteristic fields so that the numerical solution more closely follows the physically admissible rotational path. Compared with the classical Roe scheme, the NPP-Roe scheme better preserves the tangential magnetic-field magnitude and the Alfvénic velocity-magnetic-field relation, while substantially mitigating the pseudo-convergent behavior observed in this benchmark. Even on relatively coarse grids, its numerical trajectory remains close to the exact rotational path, demonstrating the path-preserving capability of the NPP-Roe scheme for the present one-dimensional test.

Having established this capability under controlled one-dimensional conditions, we next assess whether the NPP-Roe scheme can be incorporated into a realistic three-dimensional background solar-wind model without introducing evident distortion into the large-scale quasi-steady solution. For this purpose, the background solar wind during CR 2203 is simulated by solving the ideal GLM-MHD equations using a finite-volume method. The left and right states at control-volume faces are obtained through least-squares reconstruction, with a divergence constraint imposed on the reconstructed magnetic-field gradients to reduce magnetic-field divergence errors.

The resulting quasi-steady solution broadly captures the principal large-scale coronal and solar-wind features, including open- and closed-field regions, major polar and isolated coronal holes, large-scale streamer-belt structures, and the overall distributions of high- and low-speed solar wind and plasma density. The modeled structures show qualitative consistency with the available SDO/AIA and SOHO/LASCO observations in their approximate locations and large-scale morphology. These comparisons provide qualitative support for the model’s ability to reproduce the principal large-scale coronal structures, rather than a quantitative validation of its observational accuracy.

From the perspective of numerical-method applicability, the stable and physically reasonable quasi-steady solution obtained in the global simulation shows that the NPP modification can be incorporated into a realistic large-scale MHD framework. More directly, the comparison with the classical Roe scheme presented in Supplementary Appendix A shows that, under otherwise identical physical and numerical settings, replacing the classical Roe solver with the NPP-Roe solver does not substantially alter the large-scale quasi-steady background solar-wind solution, although small localized differences remain. This result supports the non-intrusive nature of the NPP modification at the scale examined here. It is also consistent with previous background solar-wind model-comparison and validation studies, which indicate that realistic global solutions are jointly shaped by multiple physical and numerical factors, including the photospheric magnetic-field input, boundary conditions, source-term treatment, grid resolution, and other components of the coupled coronal-heliospheric modeling chain (Gressl et al., 2014; MacNeice et al., 2018; Pomoell and Poedts, 2018; Samara et al., 2021; Reiss et al., 2023). More recent studies further emphasize the influence of magnetic-map inputs and surface flux-transport models on global solar-wind simulations and forecasting performance (Knizhnik et al., 2024; Wu et al., 2024; Heinemann et al., 2025).

Accordingly, the one-dimensional benchmark and the three-dimensional background solar-wind simulation serve complementary purposes. The former directly demonstrates the capability of the NPP-Roe scheme to preserve the Alfvén-wave path under controlled conditions, whereas the latter shows that the scheme can be incorporated into a realistic large-scale MHD framework while maintaining a physically reasonable quasi-steady background solar-wind solution.

The present study further motivates dedicated investigation of localized regions characterized by pronounced magnetic-field rotations. Because the NPP-Roe scheme is designed to preserve the tangential magnetic-field rotational path associated with Alfvén waves, its influence may be more pronounced in such regions. In realistic solar-wind configurations, candidate regions include areas near the heliospheric current sheet and fast-slow stream interaction regions, although these structures are not necessarily pure Alfvén rotational discontinuities. The present global comparison shows only small localized differences between the Roe and NPP-Roe solutions, and the specific contribution of the NPP correction to these structures cannot be quantified without dedicated local diagnostics and controlled analysis. Such an investigation is therefore beyond the scope of the present work.

Although the NPP-Roe scheme is applicable to large-scale background solar wind simulations, certain limitations remain in resolving small-scale and rapidly evolving local features, particularly near streamer cusps, open-closed field boundaries, and localized regions with strong magnetic-field gradients. Brightness-enhanced regions in the simulation tend to appear broader and less sharp compared to observations, indicating that the model underestimates fine-scale coronal structures. These findings suggest that further improvements are needed to enhance the accuracy of the model in representing complex high-latitude magnetic-field structures and small-scale dynamic features.

Several aspects of this study require further investigation and refinement. Future improvements may proceed in the following three directions. First, the capability of the model to capture small-scale features should be enhanced. AMR can be employed to locally refine the grid in regions with strong variations in the flow field and local structures, thereby improving the resolution of fine-scale details (Feng et al., 2014a). Second, the treatment of boundary conditions plays a critical role in determining the quality of numerical results. A self-consistent boundary condition can be adopted and combined with data-driven methods to optimize the boundary conditions in real time, ensuring consistency among physical variables and further improving the simulation accuracy of the model in high-latitude regions and small-scale structures (Yang and Shen, 2021). Third, future work will explore a systematic comparison between the classical Roe and NPP-Roe schemes under identical physical and numerical settings, with particular attention to localized structures involving strong magnetic-field rotations, to quantify the effects of the NPP modification on the global solar wind solution, Alfvénic physical consistency, and numerical stability.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fspas.2026.1909750/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Haopeng Wang, KU Leuven, Belgium

Reviewed by: José Juan González-Avilés, National Autonomous University of Mexico, Mexico

Hao Wu, Harbin Institute of Technology, Shenzhen, China