A NEW THREE-DIMENSIONAL SOLAR WIND MODEL IN SPHERICAL COORDINATES WITH A SIX-COMPONENT GRID

Magnetohydrodynamics (MHD) equations are the only self-consistent mathematical descriptions currently used to model large-scale space weather phenomena, and numerical MHD simulations are a powerful theoretical approach for retrieving the three-dimensional (3D) structures and dynamics of the solar wind in solar-terrestrial space. Numerical modeling for solar wind in the Sun–Earth system is an important part of space weather. Great progress has been made in recent years, as can be seen from survey articles (e.g., Feng et al. 2011a, 2013). Without exhausting a lengthy list of references, we would like to mention a few of the existing models. The MHD codes for the governing MHD equations of the solar wind plasma and magnetic field are naturally developed in spherical, cylindrical, Cartesian, or even in curvilinear coordinates. For example, some existing 3D MHD codes of solar wind models are written in spherical coordinates (r, θ, ϕ). Usually, in such codes, the model's computational domain is decomposed into two parts along the heliocentric distance: the inner region I (from 1 R_s to certain solar radius at the superAlfvénic region) for the coronal model and the outer region II (from the superAlfvénic region to 1 AU or beyond) for the heliospheric model; different numerical schemes for these two kinds of models can be used or combined, and the solar coronal model provides the inner boundary conditions to drive the heliospheric model. To name a few, the CORonal and HELiospheric (CORHEL) MHD model developed by the Center for Integrated Space weather Modeling (CISM) couples the solar coronal MHD model (MHD-Around-a-Sphere model, MAS model) and the heliospheric MHD model (ENLIL model; Linker et al. 1999; Mikić et al. 1999; Mikić & Linker 1994; Odstrcil et al. 2004, 2005; Riley et al. 2001, 2011, 2012, 2013). The MAS model uses a finite-difference upwind scheme in the (r, θ) directions and a pseudospectral method in the ϕ direction, while the ENLIL model adopts a Total Variation Diminishing Lax-Friedrich (TVDLF) scheme in spherical coordinates with a computational domain covering 30° to 150° in meridional angle and 0° to 360° in azimuthal angle. Another 3D MHD code, the Corona interplanetary MHD (COIN) model (Feng et al. 2002, 2003, 2005; Shen et al. 2007, 2011) is created by combining the TVDLF scheme with the MacCormack II scheme in spherical coordinates for the inner and outer regions. Usmanov (1993, 1996) and Usmanov & Goldstein (2003) established the solar wind model, in which time-dependent MHD equations are solved using a MacCormack or TVDLF scheme for region I and a steady-state MHD equation solved with the numerical scheme of MacCormack in region II (Pizzo 1982). Hayashi (2005, 2013) and Hayashi et al. (2006, 2008) developed a 3D MHD simulation code in spherical coordinates based on the concepts of total variation diminishing (TVD), the monotonic upstream scheme for conservation laws (MUSCL), and the finite volume (FV) method. In cylindrical coordinates, a 3D MHD simulation model has been developed by Nakamizo et al. (2009) for the solar wind system where the MHD equations are integrated in time using the FV method, and a TVD-type MUSCL interpolation with a Van-Leer's differentiating limiter is implemented (Tanaka 1994).

A Block-Adaptive Tree Solarwind Roe-type Upwind Scheme (BATS-R-US) code has been developed for solar coronal and heliospheric studies at the Center for Space Environment Modeling (CSEM) using the FV scheme with upwind Roe approximate Riemann solvers (Tóth et al. 2005; Van der Holst et al. 2010; Tóth et al. 2012), Harten-Lax-van Leer-Einfeldt (HLLE; Harten et al. 1983), and modified HLLE (Linde 2002) approximate Riemann solvers. With the implementation of adaptive mesh refinement, the BATS-R-US model solves the 3D MHD equations on either a Cartesian or a spherical grid. As opposed to the routinely used schemes such as FV and the finite difference schemes, the solar-interplanetary solar wind CESE-MHD (SIP-CESE-MHD) model (Feng et al. 2007, 2010; Feng et al. 2011b; Zhou et al. 2012; Feng et al. 2012a, 2012b) uses a novel numerical conservative element and solution element (CESE) method within a parallel block-based adaptive simulation in a Cartesian or curvilinear coordinates framework by treating space and time as one entity when calculating the flux balance without the Riemann solver. Basically, all of the above-mentioned 3D MHD models are able to generate the large-scale ambient solar wind by using the observational synoptic map of the photospheric and Parker solar wind solution as input.

Additionally, the spherical shell is a natural choice for the computational domain in solar-terrestrial space. For this kind of spherical shell computational domain, we should consider some requirements, such as body-fitting and polar singularity-free properties. In fact, a polyhedron-splitting method with a sphere-fitting property is usually used to construct the grid on the spherical surface. Feng et al. (2007) constructed a non-overlapping polyhedron grid by starting with an ordinary octahedron inscribed inside the spherical surface in order to avoid polar singularities. Nakamizo et al. (2009) used an angularly unstructured grid system with increasing radial spacing based on the dodecahedron-splitting method (see also Tanaka 1994), and covered the spherical surface with triangles of the same size without apparent singularities for spherical coordinates. In order to avoid coordinate singularity and grid convergence near the poles in routine spherical coordinates, various composite meshes have been applied to the solar wind plasma simulation. Usmanov (1996) and Usmanov et al. (2012) adopted a composite grid consisting of three overlapping spherical meshes. In Usmanov's composite grid, the first spherical mesh limits an extension in latitude (42° ⩽ θ₁ ⩽ 138°, 0° ⩽ ϕ₁ ⩽ 360°) with two other meshes (36° ⩽ θ_{2, 3} ⩽ 144°, 26° ⩽ ϕ₂ ⩽ 144° and 216° ⩽ ϕ₃ ⩽ 324°) covering the polar regions where the polar axes of the two polar meshes lie in the equatorial plane of the first coordinate system, that is, θ₁ = 90° and ϕ₁ = 90°. Feng et al. (2012a, 2011b) employed the Yin-Yang overset grid (Kageyama & Sato 2004; Kageyama 2005) with two identical component grids that cover a spherical surface with partial overlap on the component's boundaries. The cubed-sphere grid proposed by Ronchi et al. (1996) constitutes a spherical surface with six non-overlapping identical components where the intersection of two sets of angularly equidistant great circles defines a gridded region, and each of the arcs of the great circles can be associated with either a vertical or horizontal coordinate line. In a parallel block-based adaptive simulation framework of the cubed-sphere mesh, Ivan et al. (2011, 2012, 2013) simulated the hyperbolic conservation laws of both the Euler and ideal MHD equations.

Theoretically, the composite mesh grid with several components is a natural parallel domain decomposition of the whole sphere, since the governing equations have the same form in each component. The three components of the overlapping spherical meshes (Usmanov 1996) are not identical, which makes it difficult to achieve efficient load balancing in parallel components. Although the Yin-Yang overset grid has two identical grid components, the parallel nature of these two components incurs a huge computational cost with each covering 90° in the θ direction and 270° in the ϕ direction. The cubed-sphere grid (Ronchi et al. 1996) is far more uniform than a polar grid in terms of the shape and size of the cells, but there exist eight weak singularities at the corners of the six identical sectors of the grid. In the solar wind simulation, efficiency should also be considered, and composite grids with identical components seem to be a better choice in parallel implementation when considering load balance. For such purposes, Feng et al. (2010) introduced a six-component grid system for solar wind MHD modeling with a composite mesh consisting of six identical components designed to envelope a spherical surface with partial overlap on the component boundaries. The six identical regions have the same metric and the component grids can be transformed into each other by coordinate transformation; parallelization can be efficiently implemented in both the radial direction and components.

In solar wind modeling, a composite grid system has been employed (Usmanov 1996; Usmanov et al. 2012; Feng et al. 2010; Feng et al. 2011b; Ivan et al. 2011) using modern numerical schemes such as CWENO, CESE, and the HLL family schemes. Recently, a semi-discrete central scheme for ideal MHD equations, designed within an FV framework without a Riemann solver or any kind of characteristic decomposition, has been introduced by Ziegler (2011, 2012) on orthogonal-curvilinear grids. In his work, to maintain the ∇ · B = 0 constraint, the constraint transport (CT) technique is used for the magnetic field by utilizing a special discretization on a staggered grid. In the present paper, we use this as our base scheme to establish the 3D solar wind model in spherical coordinates on the six-component grid system.

This paper is organized as follows. In Section 2, model equations for solar wind plasma in spherical coordinates are described. Section 3 is devoted to the mesh grid system and geometrical source term discretization. Section 4 introduces the numerical scheme for the MHD equations. Section 5 gives the initial and boundary conditions in the code. Section 6 presents the numerical results for the steady-state solar wind structure of Carrington rotation (CR) 2064. Finally, we present some conclusions in Section 7.

The magnetic field B = B₁ + B₀ is split into the sum of a time-independent potential magnetic field B₀ and a time-dependent deviation B₁ (Tanaka 1994; Feng et al. 2010). We note that MHD equations can be viewed as a combination of fluid dynamics coupled with magnetic fields. In the present paper, this physical splitting of the MHD equations into fluid and the magnetic parts (Ziegler 2004; Fuchs et al. 2009) is considered in order to design efficient FV schemes with spatial discretization for the fluid equations and the magnetic induction equation adopted from Ziegler (2011). The fluid part of the vector U = (ρ, ρv_r, ρv_θ, ρv_ϕrsin θ, e)^T reads as follows:

$$\begin{equation} \frac{\partial {\mathbf {U}}}{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}r^2\mathbf {F}+ \frac{1}{r\sin \theta }\frac{\partial }{\partial \theta }\sin \theta \mathbf {G}+ \frac{1}{r\sin \theta }\frac{\partial }{\partial \phi }\mathbf {H}=\mathbf {S} \end{equation} \tag{ 1 }$$

$$\begin{equation*} \mathbf {F}=\left(\begin{array}{@{}c@{}}\rho v_r\\ \rho v_r^2+p+\frac{-B_{1r}^2+B_{1\theta }^2+B_{1\phi }^2}{2\mu }+\frac{-B_{1r}B_{0r}+B_{1\theta }B_{0\theta }+B_{1\phi }B_{0\phi }}{\mu }\\ \rho v_r v_\theta -\frac{B_{1r}B_{1\theta }+B_{1r}B_{0\theta }+B_{0r}B_{1\theta }}{\mu }\\ \left(\rho v_r v_{\phi }-\frac{B_{1r}B_{1\phi }+B_{1r}B_{0\phi }+B_{0r}B_{1\phi }}{\mu }\right) r\sin \theta \\ \left(\frac{1}{2}\rho v^2+\frac{\gamma p}{\gamma -1}\right)v_r+\frac{B_{1\theta }}{\mu }(v_rB_\theta -v_\theta B_r)+\frac{B_{1\phi }}{\mu }(v_rB_\phi -v_\phi B_r)\\ \end{array}\right) \end{equation*}$$

$$\begin{equation*} \mathbf {G}=\left(\begin{array}{@{}c@{}}\rho v_\theta \\ \rho v_r v_\theta -\frac{B_{1r}B_{1\theta }+B_{1r}B_{0\theta }+B_{0r}B_{1\theta }}{\mu }\\ \rho v_\theta ^2+p+\frac{B_{1r}^2-B_{1\theta }^2+B_{1\phi }^2}{2\mu }+\frac{B_{1r}B_{0r}-B_{1\theta }B_{0\theta }+B_{1\phi }B_{0\phi }}{\mu }\\ \left(\rho v_\theta v_{\phi }-\frac{B_{1\theta }B_{1\phi }+B_{1\theta }B_{0\phi }+B_{0\theta }B_{1\phi }}{\mu }\right)r\sin \theta \\ \left(\frac{1}{2}\rho v^2+\frac{\gamma p}{\gamma -1}\right)v_\theta +\frac{B_{1r}}{\mu }(v_\theta B_r-v_r B_\theta)+\frac{B_{1\phi }}{\mu }(v_\theta B_\phi -v_\phi B_\theta)\\ \end{array}\right) \end{equation*}$$

$$\begin{equation*} \mathbf {H}=\left(\begin{array}{@{}c@{}}\rho v_\phi \\ \rho v_r v_\phi -\frac{B_{1r}B_{1\phi }+B_{1r}B_{0\phi }+B_{0r}B_{1\phi }}{\mu }\\ \rho v_\theta v_{\phi }-\frac{B_{1\theta }B_{1\phi }+B_{1\theta }B_{0\phi }+B_{0\theta }B_{1\phi }}{\mu }\\ \left(\rho v_\phi ^2+p+\frac{B_{1r}^2+B_{1\theta }^2-B_{1\phi }^2}{2\mu }+\frac{B_{1r}B_{0r}+B_{1\theta }B_{0\theta }-B_{1\phi }B_{0\phi }}{\mu }\right)r\sin \theta \\ \left(\frac{1}{2}\rho v^2+\frac{\gamma p}{\gamma -1}\right)v_\phi +\frac{B_{1r}}{\mu }(v_\phi B_r-v_r B_\phi)+\frac{B_{1\theta }}{\mu }(v_\phi B_\theta -v_\theta B_\phi)\\ \end{array}\right) \end{equation*}$$

$$\begin{equation*} e=\frac{1}{2}\rho \mathbf {v}^2+\frac{p}{\gamma -1}+\frac{1}{2}\mathbf {B}_1^2, \ \ \ \ \mathbf {S}=\mathbf {S}_1+\mathbf {S}_2 \end{equation*}$$

with

$$\begin{equation*} \mathbf {S}_1= \left(\begin{array}{@{}c@{}}0\\ \rho \frac{v_\theta ^2+v_\phi ^2}{r}+\frac{2p}{r}+\frac{B_{1r}^2+2B_{1r}B_{0r}}{r\mu }\\ S_{1,3}\\ 0\\ 0\\ \end{array}\right) \end{equation*}$$

$$\begin{equation*} \mathbf {S}_2=\left(\begin{array}{@{}c@{}}0\\ -\rho \frac{GM_s}{r^2}+\rho \Omega \sin \theta (2v_\phi +\Omega r\sin \theta)+S_M\\ \rho \Omega \cos (2v_\phi +\Omega r\sin \theta)\\ -2\rho \Omega (v_\theta \cos \theta +v_r\sin \theta)r\sin \theta \\ \rho v_r\left(-\frac{GM_s}{r^2}+\Omega ^2r\sin ^2\theta\right)+\rho v_\theta \Omega ^2r\sin \theta \cos \theta +S_E+v_rS_M\\ , \end{array}\right) \end{equation*}$$

where

$$\begin{eqnarray*} S_{1,3}&=&\left(p+\rho v_\phi ^2+\frac{B_{1r}^2+2B_{1r}B_{0r}+B_{1\theta }^2+2B_{1\theta }B_{0\theta }-B_{1\phi }^2-2B_{1\phi }B_{0\phi }}{2\mu }\right)\frac{\cot \theta }{r}\nonumber\\ &&+\frac{1}{r}\left(\frac{B_{1\theta }B_{1r}+B_{1\theta }B_{0r}+B_{0\theta }B_{1r}}{\mu }-\rho v_r v_\theta\right). \end{eqnarray*}$$

Here, e corresponds to the modified total energy density consisting of the kinetic, thermal, and magnetic energy densities (written in terms of B₁). ρ is the mass density, v = (v_r, v_θ, v_ϕ) are the flow velocities in the frame rotating with the Sun, p is the thermal pressure, and B = B₀ + B₁ denotes the total magnetic field consisting of the time-independent potential magnetic field B₀ and its time-dependent derived part B₁. Since B₀ is constant with time, many terms near B₀ on the right-hand side vanish. t and r are time and position vectors originating at the center of the Sun. μ = 4 × 10⁻⁷π, g = −(GM_s/r³)r is the solar gravitational force, G = 6.673 × 10⁻¹¹ m³ s⁻² kg⁻¹ is the gravitational constant, M_s = 1.99 × 10³⁰ kg is the solar mass, and |Ω| = 2π/26.3 rad day⁻¹ is the solar angular speed. In our code, we allow the ratio of specific heats γ to vary from 1.05 to 1.5 along the heliocentric distance r according to Feng et al. (2010), that is, γ = 1.05 for r/R_s ⩽ 5, γ = 1.05 + 0.03(r/R_s − 5) for 5 < r/R_s ⩽ 20, and γ = 1.5 for r/R_s > 20.

In the above formulas, the first source term, S₁ = (S_{1, 1}, ⋅⋅⋅, S_{1, 5})^T, arises from the polar geometrical factors while the second source term, S₂ = (S_{2, 1}, ⋅⋅⋅, S_{2, 5})^T, arises due to the Coriolis, centrifugal, and gravity forces and volumetric heating source terms. S_M and S_E stand for the momentum and energy source terms, which are responsible for the coronal heating and solar wind acceleration. Following Feng et al. (2010), the source terms S_M and S_E are given as follows:

$$\begin{equation*} S_M=M\left(\frac{r}{R_s}-1\right)\exp \left(-\frac{r}{L_M}\right)\qquad\qquad\qquad \end{equation*}$$

$$\begin{equation*} S_E=Q_1 \exp \left(-\frac{r}{L_{Q_1}}\right)+Q_2(r-1)\exp \left(-\frac{r}{L_{Q_2}}\right), \end{equation*}$$

where Q₂ = Q₀C_a, M = M₀C_a, and $C_a=C_a^{\prime }/max(C_a^{\prime })$ with $C_a^{\prime }=({(5.8-1.6e^{[1-(\theta _b/8.5)^3]})^{3.5}}/{(1+f_s)^{2/7}})$, M₀, Q₀, Q₁ in this paper are given as 7.9 × 10⁻¹⁴Nm⁻³, 1.18 × 10⁻⁷Jm⁻³s⁻¹, 1.5 × 10⁻⁹Jm⁻³s⁻¹, respectively. $L_{Q_1}, L_{Q_2}, {\rm and}\ L_M$ are set to be 1 R_s, and f_s is the magnetic expansion factor which reads $f_s=({1}/{R})^2({B_{R_s}}/{B_R})$, where $B_{R_s}$ and B_R are the magnetic field strength at the solar surface and at the heliocentric distance R = 2.5 R_s. θ_b is the minimum angular separation between an open magnetic field foot point and its nearest coronal hole boundary.

It should be noted that we replace the momentum equation ρv_ϕ with the angular momentum equation ρv_ϕrsin θ (Kley 1998). Obviously, the geometrical curvature in spherical coordinates will turn out to produce new terms in the momentum equations, whose discretization must be accounted for and will be discussed in Section 3.2 after the introduction of the grid system. As pointed out by Ziegler (2004) and Fuchs et al. (2009), in our practical realization of numerical code design, the fluid part of the MHD is considered to be an extended Euler system with the magnetic force as a source.

The subsystem for magnetic induction part runs as follows:

$$\begin{equation} \frac{\partial B_{1r}}{\partial t}+\frac{1}{r\sin \theta }\frac{\partial }{\partial \theta }(\sin \theta (v_\theta B_r-v_r B_\theta))- \frac{1}{r\sin \theta }\frac{\partial }{\partial \phi }(v_r B_\phi -v_\phi B_r)=0\quad \end{equation} \tag{ 2 }$$

$$\begin{equation} \qquad\frac{\partial B_{1\theta }}{\partial t}-\frac{1}{r}\frac{\partial }{\partial r}(r(v_\theta B_r-v_r B_\theta)) +\frac{1}{r\sin \theta }\frac{\partial }{\partial \phi }(v_\phi B_\theta -v_\theta B_\phi)=0 \ \end{equation} \tag{ 3 }$$

$$\begin{equation} \qquad\qquad\frac{\partial B_{1\phi }}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r }(r(v_r B_\phi -v_\phi B_r)) -\frac{1}{r}\frac{\partial }{\partial \theta }(v_\phi B_\theta -v_\theta B_\phi)=0. \ \ \end{equation} \tag{ 4 }$$

As usual, ρ, v, p, B, r, and t are normalized by the characteristic values $\rho _s, a_s, \rho _s a_s^2, \sqrt{\mu \rho _s a_s^2}, R_s$, and R_s/a_s, where ρ_s and a_s are the density and sound speed at the solar surface.

This section is divided into two subsections, which are devoted to the introduction of the mesh grid system used in our code and to the consideration of geometrical source term discretization.

3.1. Six-component Mesh Grid

The computational domain is divided into a composite mesh consisting of six identical component meshes designed to envelop a spherical surface with partial overlap on their boundaries (Figure 1). Each component grid is a low-latitude spherical mesh, defined in spherical coordinates by

$$\begin{equation} \left(\frac{\pi }{4}-\delta \le \theta \le \frac{3\pi }{4}+\delta\right)\cap \left(\frac{3\pi }{4}-\delta \le \phi \le \frac{5\pi }{4}+\delta\right), \end{equation} \tag{ 5 }$$

where δ is proportionally dependent on the grid spacing entailed for the minimum overlapping area and is taken to be 2Δθ. Each component is confined to the same region as that in Equation (5) but in different spherical coordinates; vectors can be transformed from each of the six components to another (Feng et al. 2010, 2012a, 2012b, 2012c). The simulation domain is partitioned into the sliding cells i, j, k given by [r_im, r_ip] × [θ_jm, θ_jp] × [ϕ_km, ϕ_kp] with spacings of Δr(i) = r_ip − r_im, Δθ(j) = θ_jp − θ_jm, and Δϕ(k) = ϕ_kp − ϕ_km. Corresponding half-way indices defined by 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) mark the bounding faces of the cell i, j, k, where i = 1, ..., N_r, j = 1, ..., N_θ, k = 1, ..., N_ϕ. The grid mesh is generated as follows: r_im = r_{(i − 1)p}, r_ip = r_im + Δr(i), θ_jm = θ_{(j − 1)p}, θ_jp = θ_jm + Δθ(j), ϕ_km = ϕ_{(k − 1)p}, ϕ_kp = ϕ_km + Δϕ(k) with r_1m = 1R_S, θ_1m = (π/4), and ϕ_1m = (3π/4). The geometrical center of a cell i, j, k is denoted by (r_i, θ_j, ϕ_k) with r_i = ((r_im + r_ip)/2), θ_j = ((θ_jm + θ_jp)/2), and ϕ_k = ((ϕ_km + ϕ_kp)/2). The volume-averaged coordinates $(\overline{r}_i, \overline{\theta }_j,\ {\rm and}\ \phi _k)$ of a cell i, j, k are $\overline{r}_i=(({3r_{ip}^4-3r_{ip}^4})/({4r_{ip}^3-4r_{ip}^3})), \overline{\theta }_j= (({\sin \theta _{jp}-\sin \theta _{jm}-(\theta _{jp}\cos \theta _{jp}-\theta _{jm}\cos \theta _{jm})})/({\cos \theta _{jm}-\cos \theta _{jp}}))$. The coordinates of the six face centers of a cell i, j, k are denoted as $(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}),\ {\rm and}\ (r_i^\phi, \theta _j^\phi, \phi _{kp})$ with $r_i^\theta =r_i^\phi = (({2r_{ip}^3-2r_{im}^3})/({3r_{ip}^2-3r_{im}^2})), \theta _j^r=\overline{\theta }_j, \theta _j^\phi =(({\theta _{jp}+\theta _{jm}})/{2})$. A typical sliding cell is displayed in Figure 2.

Zoom In Zoom Out Reset image size

In the present work, the parallel implementation over the whole computational domain of our simulated region from 1 R_s to 230 R_s is realized by domain decomposition of six-component grids based on the spherical surface and radial direction partition. Due to the radial expansion of this vast spherical shell domain, the MHD governing equations are discretized on highly non-uniform grids in order to resolve the steep gradients in the solution. In order to mitigate this discrete or geometrical stiffness caused by disparate mesh cell widths, the following grid partitions are employed. For 1–25 R_s, N_θ = N_ϕ = 42, Δr(i) = 0.01 R_s if r(i) < 1.1 R_s; Δr(i) = min(A × log₁₀(r(i − 1)), Δθ × r(i − 1)) with A = 0.01/log₁₀(1.09) if r(i) < 3.5 R_s; Δr(i) = Δθ × r(i − 1) if r(i) > 3.5 R_s. For 25-75 R_s, N_θ = N_ϕ = 64 and Δr(i) = Δθ × r(i − 1). For 75-166 R_s, N_θ = N_ϕ = 96 and Δr(i) = Δθ × r(i − 1). For 166-230 R_s, N_θ = N_ϕ = 130 and Δr(i) = Δθ × r(i − 1).

3.2. Geometrical Source Term Discretization

As we know, in polar coordinates, the conservation laws of vector quantities such as the momentum of a fluid parcel yield extra terms in the equations if an additional force is present. Some extra source terms in the momentum equations of the fluid-part Equation (1) will appear if written in the usual conservative form (Usmanov 1993; Feng et al. 2005), compared to the Cartesian coordinate system. Because source terms are always a menace, reducing or avoiding them is a positive thing to do in spatial discretizations. For example, the ϕ momentum equation in the usual conservative form for ρv_ϕ can be written as

$$\begin{eqnarray} &&\frac{\partial \rho v_\phi }{\partial t}+ \frac{1}{r^2}\frac{\partial {r^2 (\rho v_r v_{\phi }-(B_{1r}B_{1\phi }+B_{1r}B_{0\phi }+B_{0r}B_{1\phi })/\mu)}}{\partial r}\nonumber\\ &&\quad+\frac{1}{r\sin \theta }\frac{\partial {\sin \theta (\rho v_\theta v_{\phi }-(B_{1\theta }B_{1\phi }+B_{1\theta }B_{0\phi }+B_{0\theta }B_{1\phi })/\mu )}}{\partial \theta }\nonumber\\ &&\quad+\frac{1}{r\sin \theta }\frac{\partial \left(\rho v_\phi ^2+p+\left(B_{1r}^2+B_{1\theta }^2-B_{1\phi }^2\right)/(2\mu) +(B_{1r}B_{0r}+B_{1\theta }B_{0\theta }-B_{1\phi }B_{0\phi })/\mu \right)}{\partial \phi }\nonumber\\ &&\qquad=S_{1,4}-2\rho \Omega (v_\theta \cos \theta +v_r\sin \theta), \end{eqnarray} \tag{ 6 }$$

where S_{1, 4} = (cos θ/rsin θ)(((B_1θB_1ϕ + B_1θB_0ϕ + B_0θB_1ϕ)/μ) − ρv_θv_ϕ) − (1/r)(ρv_rv_ϕ − ((B_1rB_1ϕ + B_1rB_0ϕ + B_0rB_1ϕ)/μ)) are the extra terms caused by geometrical curvature in the spherical coordinates. In this paper, we replace the ϕ momentum equation with the following angular momentum equation (Kley 1998):

$$\begin{eqnarray} &&\frac{\partial \rho v_\phi r \sin \theta }{\partial t}+ \frac{1}{r^2}\frac{\partial {r^2 (\rho v_r v_{\phi }-(B_{1r}B_{1\phi }+B_{1r}B_{0\phi }+B_{0r}B_{1\phi })/\mu )r \sin \theta }}{\partial r}\nonumber\\ &&\quad+\frac{1}{r\sin \theta }\frac{\partial {\sin \theta (\rho v_\theta v_{\phi }-(B_{1\theta }B_{1\phi }+B_{1\theta }B_{0\phi }+B_{0\theta }B_{1\phi })/\mu )r \sin \theta }}{\partial \theta }\nonumber\\ &&\quad+\frac{1}{r\sin \theta }\frac{\partial \left(\rho v_\phi ^2+p+\left(B_{1r}^2+B_{1\theta }^2-B_{1\phi }^2\right)/(2\mu) +(B_{1r}B_{0r}+B_{1\theta }B_{0\theta }-B_{1\phi }B_{0\phi })/\mu \right)r \sin \theta }{\partial \phi }\nonumber\\ &&\qquad=S_{2,4}\equiv -2\rho \Omega (v_\theta \cos \theta +v_r\sin \theta)r \sin \theta. \end{eqnarray} \tag{ 7 }$$

Then, we can see that the extra terms S_{1, 4} in Equation (6) disappear completely in Equation (7) for the angular momentum ρv_ϕrsin θ, and can easily be incorporated into the momentum flux difference.

In addition, the extra source terms in the r and θ momentum equations should be treated carefully. Now, let us take a look at how we can formulate the conservation Equation (1) in polar coordinates. Equation (1) has the discrete form (∂U_{i, j, k}V_{i, j, k}/∂t) + (Δ_rFS_r)_{i, j, k} + (Δ_θGS_θ)_{i, j, k} + (Δ_ϕHS_ϕ)_{i, j, k} = S_{i, j, k}V_{i, j, k}, and (Δ_rFS_r)_{i, j, k} = F_{ip, j, k}S_{(r), ip, j, k} − F_{im, j, k}S_{(r), im, j, k} denotes the difference of the argument FS_r between the ip and im interfaces, and the same holds for Δ_θ and Δ_ϕ. Here, S_{(r), ip, j, k}, S_{(θ), i, jp, k}, and S_{(ϕ), i, j, kp} denote the surface areas ip, jp, and kp of cell i, j, k, and V_{i, j, k} is the volume of cell i, j, k. Remember that 2p/r in S_{1, 2} for the r momentum equation is generated when we convert ∂p/∂r into (1/r²)(∂r²p/∂r), and the term $({\cot \theta }/{r}) p$ in S_{1, 3} for the θ momentum equation appears when we convert (1/r)(∂p/∂θ) into (1/rsin θ)(∂sin θp/∂θ). In order to preserve discrete consistency, we should compute such source terms following the flux differences in the r and θ directions by replacing (2p/r)V_{i, j, k} by pΔ_rS_r = (S_{(r), ip, j, k} − S_{(r), im, j, k})p_{i, j, k} and $({\cot \theta }/{r}) p V_{i,j,k}$ by pΔ_θS_θ = (S_{(θ), i, jp, k} − S_{(θ), i, jm, k})p_{i, j, k}. Otherwise, one can easily find that the case p(r, θ, ϕ) = constant and (v_r, v_θ, v_ϕ) = 0 in the non-rotating frame cannot be recognized by the solver as a static solution to the pure Euler equations, that is, Equation (1), without considering the magnetic field, solar rotation, solar gravity, and heating source terms.

Following Ziegler (2004) and Fuchs et al. (2009), we design an FV scheme for the governing equations by splitting these equations into a fluid part, Equation (1), and a magnetic induction part, Equations (2)–(4). The fluid part leads to an extended Euler system with magnetic forces as source terms. For spatial discretization of our numerical scheme formulation, we strictly follow those of Ziegler (2011) by using the FV discretization of Equation (1), and by averaging Equations (2)–(4) over facial areas to obtain the semi-integral forms of magnetic induction equations. We used a semi-discrete Godunov-type central scheme for the Euler subsystem. We adopt a second-order accurate linear ansatz reconstruction (Ziegler 2011) for the fluid part, and in the linear reconstruction the derivative terms at volume-averaged coordinates are approximated by using a minmod limiter for oscillation control. The magnetic induction equations are solved using the CT method. For consistency, reconstruction of the magnetic induction part is also of second-order accuracy and cross derivative terms are slope-limited approximations of the exact derivatives at the cell face centers obtained using a minmod limiter. For details, refer to Ziegler (2011).

The full system is integrated over time with a second-order Runge–Kutta scheme. Following the idea of Ziegler (2004) and Fuchs et al. (2009), we simultaneously or sequentially integrate the MHD equations in time for the fluid part and the magnetic induction part. That is, in order to obtain a scheme for the full MHD equations, we move the discretized fluxes to the right-hand sides of the governing Equations (1)–(4), and denote them by their corresponding source terms as $\mathcal {R}_{\mathbf {U}}[\overline{\mathbf {U}},\overline{\mathbf {B}}]$ and $\mathcal {R}_{\mathbf {B}}[\overline{\mathbf {U}}, \overline{\mathbf {B}}]$, respectively, and we can advance using the following procedure:

$$\begin{equation*} \overline{\mathbf {U}}^\ast =\overline{\mathbf {U}}^n+\Delta t\mathcal {R}_{\mathbf {U}}[\overline{\mathbf {U}}^n,\overline{\mathbf {B}}^n] \end{equation*}$$

$$\begin{equation*} \overline{\mathbf {B}}^\ast =\overline{\mathbf {B}}^n+\Delta t\mathcal {R}_{\mathbf {B}}[\overline{\mathbf {U}}^n,\overline{\mathbf {B}}^n] \end{equation*}$$

$$\begin{equation} \overline{\mathbf {U}}^{n+1}=\frac{1}{2}\overline{\mathbf {U}}^n+\frac{1}{2}(\overline{\mathbf {U}}^\ast +\Delta t\mathcal {R}_{\mathbf {U}}[\overline{\mathbf {U}}^\ast, \overline{\mathbf {B}}^\ast ]) \end{equation} \tag{ 8 }$$

$$\begin{equation} \overline{\mathbf {B}}^{n+1}=\frac{1}{2}\overline{\mathbf {B}}^n+\frac{1}{2}(\overline{\mathbf {B}}^\ast +\Delta t\mathcal {R}_{\mathbf {B}}[\overline{\mathbf {U}}^\ast,\overline{\mathbf {B}}^\ast ]). \end{equation} \tag{ 9 }$$

In a cell with a geometrical center of i, j, k, $\overline{\mathbf {U}}=\mathbf {U}(\overline{r}_i, \overline{\theta }_j, \phi _k)$ denotes the volume average for the conservative variables U and $\overline{\mathbf {B}}=(B_r(r_{im},\theta _j^r,\phi _k),B_\theta (r^\theta,\theta _{jm},\phi _k),B_\phi (r^\phi,\theta _j^\phi,\phi _{km})$ denotes the area average for the magnetic field B_r(B_θ, B_ϕ) in face r(θ, ϕ). In the spatial discretization for the magnetic induction equations, all of the reconstruction is based on the staggered magnetic field components without a cell-centered magnetic field involved. Thus, there is only one set of magnetic field updated at the face center. The face-valued magnetic fields $\overline{\mathbf {B}}$ stored on the faces are only averaged to the zone center value B when we solve the pressure equation for positivity consideration in Section 4.1 below.

We can simultaneously piece together the scheme for the fluid part with the scheme for the induction equation. Or sequentially, Equation (9) can be replaced by the update formula

$$\begin{equation*} \overline{\mathbf {B}}^{n+1}=\frac{1}{2}\overline{\mathbf {B}}^n+\frac{1}{2}(\overline{\mathbf {B}}^\ast +\Delta t\mathcal {R}_{\mathbf {B}}[\overline{\mathbf {U}}^{n+1},\overline{\mathbf {B}}^\ast ]). \end{equation*}$$

Furthermore, the order in the above sequential update can also be reversed by first evolving the magnetic part, followed by the evolution of the fluid part according to Fuchs et al. (2009). As usual, the time step length is prescribed by the Courant–Friedrichs–Lewy (CFL) stability condition:

$$\begin{equation*} \Delta t ={\rm CFL}/ {\rm max}\left(\sqrt{\left(\frac{|v_r|+c_{fr}}{\Delta r}\right)^2+\left(\frac{|v_\theta |+c_{f\theta }}{r\Delta \theta }\right)^2+ \left(\frac{|v_\phi |+c_{f\phi }}{r \sin \theta \Delta \phi }\right)^2}\right). \end{equation*}$$

Here, c_fr, c_fθ, andc_fϕ are the fast magnetosonic speeds in the (r, θ, ϕ) directions, defined, respectively, by $c_{fr}=({1}/{2})\sqrt{\vphantom{A^A}\smash{{c_s^2+c_A^2+((c_s^2+c_A^2)^2-4c_s^2({B_r^2}/{\mu \rho }))^{{1}/{2}}}}}$, $c_{f\theta }=({1}/{2})\sqrt{\vphantom{A^A}\smash{{c_s^2+c_A^2+((c_s^2+c_A^2)^2-4c_s^2({B_\theta ^2}/{\mu \rho }))^{{1}/{2}}}}}$, $c_{f\phi }= ({1}/{2})\sqrt{\vphantom{A^A}\smash{{c_s^2+c_A^2+((c_s^2+c_A^2)^2-4c_s^2({B_\phi ^2}/{\mu \rho }))^{{1}/{2}}}}}$, where $c_s=\sqrt{{\gamma p}/{\rho }}\ {\rm and}\ c_A=\sqrt{\vphantom{A^A}\smash{{ ({B_r^2+B_\theta ^2+B_\phi ^2})/{\mu \rho }}} }$ are the sound and Alfvénic speeds. In the following run, we employ a simultaneous time integration with CFL = 0.5.

As pointed out by Feng et al. (2010), the plasma density, Alfvén velocity, interplanetary magnetic fields (IMFs), plasma β, and spatial grid size vary over many orders of magnitude from the Sun to Earth. This implies a large variation of the CFL stability limit from the corona to interplanetary space. If we apply one uniform time step in the entire solar-terrestrial domain, the time accuracy stability limitations on the time step in the fine grid parts of the computational domain decrease the performance of the numerical methods, and it may become much more difficult to achieve efficient load balancing. The present paper uses the multi-time-stepping proposed by Feng et al. (2010) by taking different time steps in different grid parts: 1–3.5 R_s, 3.5–10 R_s, 10–25 R_s, 25–75 R_s, 75–166 R_s, and 166–230 R_s. This multi-time-stepping method somehow avoids the necessity of taking a single time step over the entire computational domain determined by the numerical CFL stability conditions (such as on the smallest grid cells near the Sun). For details, refer to Feng et al. (2010).

4.1. Maintaining Pressure Positivity

When the MHD equations are written in conservative form, the pressure p is derived by the total energy e. One must subtract off the kinetic energy and magnetic energy from the total e. In the low β plasma regions, the plasma kinetic energy and magnetic energy are larger than the pressure, and thus the discretization errors generated when computing the total energy minus kinetic energy and magnetic energy to get the pressure can be larg enough to result in negative pressure. This condition is an unacceptable physical situation in the MHD flows. Therefore, special treatment is necessary in these possible negative regions. Although we use $e=({1}/{2})\rho \mathbf {v}^2+({p}/({\gamma -1}))+({1}/{2})\mathbf {B}_1^2$, this modified energy may mitigate such occurrences to some extent. We still want to enhance this positivity. Here, we refer to Balsara & Spicer (1999) to correct the possible occurrence of negativity by using the pressure equation instead of the energy equation. That is, we use the pressure equation (∂p/∂t) + ∇ · vp = −(γ − 1)p∇ · v + (γ − 1)S_E when the following three conditions are satisfied.

Condition I:

$$\begin{equation*} \overline{p}_{i,j,k}<\alpha _1 \overline{e}_{i,j,k}. \end{equation*}$$

Condition II:

$$\begin{equation*} |\overline{p}_{i+1,j,k}-\overline{p}_{i-1,j,k}|+|\overline{p}_{i,j+1,k}-\overline{p}_{i,j-1,k}| +|\overline{p}_{i,j,k+1}-\overline{p}_{i,j,k-1}| \end{equation*}$$

$$\begin{equation*} <\alpha _2 {\rm min} (\overline{p}_{i,j,k}, \overline{p}_{i+1,j,k}, \overline{p}_{i-1,j,k},\overline{p}_{i,j+1,k},\overline{p}_{i,j-1,k},\overline{p}_{i,j,k+1}, \overline{p}_{i,j,k-1}). \end{equation*}$$

Condition III:

$$\begin{equation*} -\alpha _3 {\rm max} (\overline{C}_{i,j,k}, \overline{C}_{i+1,j,k}, \overline{C}_{i-1,j,k},\overline{C}_{i,j+1,k},\overline{C}_{i,j-1,k},\overline{C}_{i,j,k+1}, \overline{C}_{i, j, k-1}) \end{equation*}$$

$$\begin{equation*} <{\rm min} (\Delta r(i), \overline{r}_i\Delta \theta (j), \overline{r}_i\sin \overline{\theta }_j\Delta \phi (k) )\overline{\nabla \cdot \mathbf {v} }_{i,j,k}, \end{equation*}$$

where

$$\begin{equation*} \overline{C}_{i,j,k}=\left(\frac{\gamma \overline{p}_{i,j,k}}{\overline{\rho }_{i,j,k}} +\frac{{B}_{i,j,k}^2}{\mu \overline{\rho }_{i,j,k}}\right)^{\frac{1}{2}}. \end{equation*}$$

Condition I indicates that the pressure is much smaller than the energy. Condition II is intended to exclude the regions of strong magnetosonic shocks. Condition III acts to exclude the regions where mildly compressive motions (also strong compressive motions) may take place. We note that α₁, α₂, and α₃ can be adjusted in different situations, and α₁ = 0.05, α = 0.3, and α₃ = 0.01 are used in the present paper.

Initially, we specify the magnetic field using the line-of-sight photospheric magnetic data from the Wilcox Solar Observatory to produce a 3D global magnetic field in the computational domain with the potential field source surface (PFSS) model. The initial distributions of plasma density ρ, pressure p, and velocity v are given by Parker's solar wind flow (Parker 1963). The temperature and number density on the solar surface are T_s = 1.3 × 10⁶ K and ρ_s = 1.5 × 10⁸ cm⁻³.

In this paper, the inner boundary at 1 R_s is fixed for simplicity. Since the outer boundary in interplanetary space is the supersonic/super-Alfvénic region, the solar wind parameters at the outer boundary are set equal to the values at their nearest grid points in the computational domain. In the six-component grid system, a horizontal boundary in the (θ, ϕ) directions exists at the six-component's borders in the overlapping parts, and the horizontal boundary or internal border values of each component grid are determined by interpolation from the neighbor stencils lying in its neighboring component grid, with the interpolation coefficients being derived by the corresponding position of the boundary point in the stencils (Feng et al. 2010), and with the second-order Lagrange polynomial interpolation used for the fluid flow and the second-order reconstruction for the magnetic field.

In this section, we present the numerical results from CR 2064 for the steady state solar wind and compare the MHD results at 1 AU with the observational data.

Zoom In Zoom Out Reset image size

Figure 3 shows the model results for the magnetic fields lines, radial speed v_r, and number density N on the meridional plane from 1 to 20 R_s, where the arrowheads on the black lines represent the directions of the magnetic fields. From this figure, it can be seen that the magnetic field lines at the high latitudes extend to interplanetary space where fast solar wind flows and low density appear. However, slow solar wind and high density are shown in the vicinity of the equator or heliospheric current sheet (HCS) region (close to the equator near the solar minimum), and a helmet streamer stretched by the solar wind can be observed within about five solar radii. Above the streamer, a thin current sheet is displayed between different magnetic polarities. These MHD results present the typical structures of the solar wind in the solar minimum.

Zoom In Zoom Out Reset image size

Figure 4 displays the number density, radial speed, and temperature maps for the MHD steady-state solution at 2.5 R_s, 20 R_s, and 215 R_s, respectively. Obviously, the number density decreases with radial distance while radial speed increases. The HCS region has high density, low speed, and low temperature, whereas the inverse condition is seen at the polar or the open field region. From the results at 215 R_s, we find that the fast solar wind from the polar is about 850 km s⁻¹ and the slow solar wind near HCS is ∼350 km s⁻¹ at 1 AU. Figure 5 displays the radial magnetic fields from the MHD and PFSS models, from which the neutral line positions for both models are almost the same.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 is the radial speed v_r at 5 R_s from the MHD and WSA results. The speed from the WSA is specified by $v_r=265+({1.5}/{(1+f_s)^{2/7}})(5.8-1.6e^{[1-(\theta _b/7.5)^3]})^{3.5}$ with the use of θ_b and the expansion factor f_s (Owens et al. 2005). From this figure, it is clear that the speed from the MHD result is smaller than that from WSA, which is because the f_s obtained from the steady-state magnetic fields of the MHD model is bigger than that of the WSA model and will have a large expansion magnetic field and then present a relatively low speed. This can be interpreted physically in terms of the pressure excised by the plasma which further spreads the magnetic field in the MHD simulation (Feng et al. 2010).

Zoom In Zoom Out Reset image size

Figure 7 is the variation of the number density N and radial speed v_r from 1 R_s to 20 R_s at different latitudes. θ = 2° is represented by the solid line and θ = 90° is represented by the dashed line, where θ = 2° corresponds to the open field region while θ = 90° corresponds to the HCS region. As usual, in the open field region, the speed is larger and stands for the fast solar wind, while in the HCS region the speed is smaller and stands for the slow solar wind, and the variation of number density is opposite to that of the speed.

Zoom In Zoom Out Reset image size

Figure 8 is the temporal evolution of error for ∇ · B₁ and ∇ · B in the calculation from 1 R_s to 31 R_s. The error for ∇ · B is defined as Error $=({\Sigma _{k=1}^{M}\int _{V_k}\mid (\nabla \cdot \mathbf {B})\mid dV}/{\Sigma _{k=1}^{M}\int _{V_k}dV})$, where M is the total number of cells in the computational domain and V_k is the kth hexahedron involved with the mesh grids. From Figure 8, we can see that the use of the CT scheme leads to the error for ∇ · B₁ around 10⁻¹⁷ and for ∇ · B of about 10⁻⁶. In practice, this error remains constant after 25 hr and no obvious large error appears after a long run time. The large part of ∇ · B comes from the initial discretizaion of ∇ · B₀. This conforms to the argument provided by Jeltsch & Torrilhon (2005).

Zoom In Zoom Out Reset image size

Figure 9 shows the MHD results for the magnetic fields and radial speed on the equatorial plane from 20 R_s to 215 R_s. The arrowheads denote the direction of the magnetic field and the color contours represent the radial solar wind speed. The solar wind extends the IMF outward into Archimedean spirals due to solar rotation and the IMF freeze-in effect.

Zoom In Zoom Out Reset image size

Figure 10 shows the comparisons between the MHD results and the OMNI observational data. From this figure, we can see that the MHD results are in good agreement with the observations. The speed reaches a first peak at days 11–12 and a second peak at days 20–22, the density peaks at days 8–10 and days 15–18, and the temperature peaks at days 9–11 and days 18–22, which can also be seen from the observations. From this figure, we can also see that the magnetic field has the same trend with the observations. Figure 11 gives the model results for the solar wind speed at 1 AU as a function of heliolatitude in the meridional plane. From this we can clearly see a transition between the fast wind and slow wind at the midlatitudes. The result is in good agreement with Ulysses speed observations when normalized to 1 AU (Usmanov et al. 2000), which states that there is low speed near the heliospheric current sheet region and high speed at both poles comparable to the profiles displayed in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By employing the FV scheme in spherical geometry with a six-component grid for the solar wind plasma MHD system, we establish a new corona-interplanetary solar wind model. In the numerical scheme design, in time integration, we follow the idea of Ziegler (2004) and Fuchs et al. (2009), and in spatial discretization we obey the procedure suggested by Ziegler (2011). The solar-wind-governing MHD equations are divided into two coupled subsystems with a fluid part and a magnetic field part. We adopt the FV discretization (Ziegler 2011) for the extended Euler system with magnetic forces as source terms, while the magnetic field part is treated by CT techniques suggested by Ziegler (2011) in order to maintain the divergence-free nature of the magnetic field. Ziegler's reconstruction (Ziegler 2011), combined with a second-order Runge–Kutta time integration, gives this new numerical scheme second-order accuracy in time and space. The composite mesh with six-component grids proposed by Feng et al. (2010) is used to partition the computational domain of the Sun-to-Earth spherical shell into six identical components with partial overlaps on their boundaries, which is able to solve the problems of coordinate singularity and grid convergence near the poles in spherical polar coordinates. The use of a six-component grid can facilitate parallelization both in the radial direction and in the (θ, ϕ) directions and reduce the computational cost. Following Feng et al. (2010), in order to achieve efficient load balancing, the multiple time-stepping method is adopted in the radial direction, and in the acceleration of solar wind, volumetric heating source terms are considered by using the expansion factor f_s and the angular distance θ_b.

We validate the new model by obtaining the large-scale structure for CR2064. The results show that this model can produce many typical properties of the solar wind, such as an obvious slow speed area near the slightly tilted HCS plane and a fast speed area near the poles, and high density in the slow speed area and vice versa in both poles, all of which makes us very confident in the current MHD code. The numerical results at 1 AU agree well with the Ulysses observation. Both numerical and observational results show that there exists a fast speed flow with a magnitude of 700 ∼ 800 km s⁻¹ at high latitudinal regions, except for a ±20° area near the slightly tilted HCS plane where the velocity is about 350 km s⁻¹. Both also demonstrate that the interface between the high- and low-speed winds is relatively steep due to the addition of volumetric heating source terms. The CT scheme used for the magnetic induction equation can keep the divergence-free condition well and also simplify the computation compared with the projection Poisson divergence cleaning procedure (Feng et al. 2010). Further improvement for the magnetic field divergence error could follow the method suggested by Jeltsch & Torrilhon (2005). Overall, our model can produce all the physical parameters everywhere within the computation domain; however, many unsatisfactory points remain, as we discuss.

Since we only use the volumetric heating to account for the coronal heating/solar wind acceleration, the radial speed (temperature) is a little small compared to OMNI observations and the separation of high flow and low flow is not very satisfactory due to the lack of a physical sound mechanism coupled to our model. However, we have reason to believe that this heating form cannot be the only acceleration process acting on the solar wind and that other presently unknown sources are needed to act within the region between the lower corona and the source surface. Further characterizing and quantifying of the key physical processes/mechanism will clarify an operational route to more physically integrate realistic coronal heating modules into 3D MHD codes. In the future, we will also consider the implementation of an adaptive mesh refinement technique (Ziegler 2005, 2012; Feng et al. 2012a, 2012b, 2014) for this model.

The work was jointly supported by the National Basic Research Program of China (grant No. 2012CB825601), the National Natural Science Foundation of China (grant Nos. 41031066, 41231068, 41274192, 41204127, and 41174150), the Knowledge Innovation Program of the Chinese Academy of Sciences (grant No. KZZD-EW-01-4), and the Specialized Research Fund for State Key Laboratories. The numerical calculation has been completed on our SIGMA Cluster computing system. The Wilcox Solar Observatory is currently supported by NASA. The OMNI data is obtained from the GSFC/SPDF OMNIWeb interface http://omniweb.gsfc.nasa.gov. We appreciate the referee's constructive suggestions to improve this paper.