THREE-DIMENSIONAL SOLAR WIND MODELING FROM THE SUN TO EARTH BY A SIP-CESE MHD MODEL WITH A SIX-COMPONENT GRID

The spherical shell is a key computational domain in both computational geophysics and solar–terrestrial physics (Kageyama & Sato 2004; Kageyama 2005; Yoshida & Kageyama 2004; Usmanov & Dryer 1995; Usmanov 1996; Wu et al. 1999; Usmanov & Goldstein 2003; Linker et al. 1999; Odstrcil et al. 2002; Feng et al. 2005; Hayashi 2005; Shen et al. 2007). Simulations in these fields such as geodynamo, mantle convection, global circulation of the atmosphere, and space weather modeling from the Sun to Earth are all performed in spherical shell geometry, such as (r, θ, ϕ) with radius r between some intervals (for example, r₀ ⩽ r ⩽ r₁, 0 ⩽ θ ⩽ π, and 0 ⩽ ϕ ⩽ 2π), if defined in the spherical coordinates.

The most convenient grid system for solar wind modeling is the latitude–longitude grid under spherical shell geometry, which is especially welcomed for numerical schemes for governing equations in the spherical polar coordinates (Usmanov & Dryer 1995; Usmanov 1996; Wu et al. 1999; Usmanov & Goldstein 2003; Linker et al. 1999; Odstrcil et al. 2002; Feng et al. 2005; Hayashi 2005; Shen et al. 2007) since its orthogonality is expected from numerical reasons. The grid mesh of the latitude–longitude grid is uniform when it is seen in the computational space of the spherical surface S = {(θ, ϕ), 0 ⩽ θ ⩽ π, 0 ⩽ ϕ ⩽ 2π} but it is in fact non-uniform when observed in real space. This is a feature common to the grid mesh in the spherical shell. Even if the numerical schemes are implemented in Cartesian or cylindrical coordinates (Tanaka 1994, 1995, 1999, 2000; Groth et al. 2000; Cohen et al. 2008; Feng et al. 2007; Hu et al. 2008; Kleimann et al. 2009; Nakamizo et al. 2009), the computational domain used in solar wind modeling has to approximate the spherical shell geometry.

Intuitively, the Sun's geometry, which is obviously spherical shaped, would suggest the use of spherical coordinates (r, θ, ϕ), especially since the radial convergence of lines of constant θ, ϕ can entail the extra advantage of enhanced spatial resolution near the Sun's surface. But, in such a spherical computational domain, as is well known, there exist two crucial numerical problems: the coordinate singularity and the grid convergence near the poles. In the latitude–longitude grid, the existence of the coordinate singularity arises on the north pole θ = 0 and the south pole θ = π, which by using L'Hospital's rule near the poles leads to three different forms for three local regions of the spherical surface S: $\nabla =(\frac{\partial }{\partial r}, \frac{\partial }{r\partial \theta }, \frac{1}{ r \sin \theta }\frac{\partial }{\partial \phi })$ for the local region (θ, ϕ), 0 < θ < π, 0 ⩽ ϕ ⩽ 2π without the poles, $\nabla =(\frac{\partial }{\partial r}, \frac{\partial }{r \partial \theta }, \frac{1}{ r }\frac{\partial ^2}{\partial \theta \partial \phi })$ for the north pole, and $\nabla =(\frac{\partial }{\partial r}, \frac{\partial }{r \partial \theta }, -\frac{1}{ r }\frac{\partial ^2}{\partial \theta \partial \phi })$ for the south pole. This observation prompts us to devise the numerical scheme for the MHD equations in three parts: a low latitude–longitude region without the poles, the north pole, and the south pole, which results in coding three different subroutines for numerical computation of discretized equations. This may destroy the hyperbolicity of our hyperbolic system such as MHD equations and uniform difference schemes cannot be directly applicable. In addition, this also increases the complexity of the computer programs.

The grid convergence near the poles is another problem in the latitude–longitude grid of spherical coordinates. The grid convergence requires a strict limitation on the time-marching step ▵t in time integration schemes. As θ → 0 or π near the poles, the metric factor 1/(rsin θ) (of the ϕ-component of the spherical gradient operator ∇) becomes extremely large with a constant grid spacing ▵ϕ in the longititudinal direction. In practice, the grid convergence brings us more difficulties than the coordinate singularity due to both the grid redundancy and the severe time-step constraint (especially for large-scale numerical simulations with fine grids) imposed by the Courant–Friedrichs–Levy (CFL) stability condition.

For spherical coordinates, this means that enhancing spatial resolution at small r, however desired for physical reasons, would yield Δt much smaller than what the CFL criterion would require for most parts of the computational domain so that such a small time step seems impractical to be realized in numerical study. While adaptive mesh refinement (AMR) can, in principle, be used with spherical coordinates, its advantage is overcompensated for by the fact that the convergence of grid spacing implies unacceptably low CFL numbers.

The problems of small time steps with coordinate singularity and grid convergence can be avoided by using Cartesian coordinates, which can have equal cell spacing everywhere and thus do save computing time on the larger cells. Many modern numerical schemes are usually established in Cartesian coordinates and can be directly used in solar wind modeling without any further modification. For these reasons, Cartesian coordinates (x, y, z) have also been widely used (Groth et al. 2000; van der Holst & Keppens 2007; Kleimann et al. 2009), for which numerics are faster, simpler (esp ecially with respect to multi-dimensional extension), and more stable. This is especially true for an MHD code built within a framework that allows for Cartesian AMR (e.g., Groth et al. 2000; van der Holst & Keppens 2007; Kleimann et al. 2009).

Evidently, the solar surface, described by a sphere of unit radius located inside the computational domain, cannot be consistent with any of the Cartesian coordinate surfaces. The inner (solar surface) boundary yields the question of how these boundary conditions are well realized on the grid since it contains the surface from which the solar wind emanates, such that numerical artifacts imposed by an improper treatment of this boundary will be quickly propagated through the entire computational domain. As pointed out by Kleimann et al. (2009), simple-minded attempts to implement these inner boundary conditions, such as keeping cell values inside the Sun (i.e., lower boundary) fixed and integrating only those outside, will not work well and will result in block-like artifacts at small radii (essentially tracing the envelope of the set of grid cells considered "inside") where spherical contours would be better required for keeping the problem's symmetry. In order to adequately implement the Sun's spherical surface as an inner boundary on the Cartesian grid, a weighted averaging procedure (Kleimann et al. 2009) was devised in order to handle the huge gradients (most notably of mass density) occurring at this boundary. To deal with the spherical surface inner boundary, an alternative method is to allow for arbitrary embedded boundaries in the Cartesian geometry by employing a cut-cell method (e.g., Zeeuw & Powell 1993; Groth et al. 2000; van der Holst & Keppens 2007). The use of these procedures also contributes to a reduction of spurious departures from the problem's underlying symmetry, which results from the fact that the Sun's spherical boundary surface cannot be mapped to a Cartesian grid of finite cell spacing.

The polyhedron-splitting method seems very efficient for generating the mesh grid on the spherical surface of a unit sphere. In three-dimensional (3D) solar wind simulation, regular octahedron (a Platonic solid composed of eight equilateral triangles, four of which meet at each vertex) and dodecahedron (a Platonic solid composed of 12 regular pentagonal faces, 20 vertices, and 30 edges, with 3 faces meeting at each vertex) splitting grids can be constructed by starting with a regular octahedron (Feng et al. 2007; Hu et al. 2008) and dodecahedron (Nakamizo et al. 2009) inscribed inside the unit spherical surface. These methods divide the sphere homogeneously into triangles of the same size with no directional preference such that grid singularities can be avoided and thus the computation stability can be effectively achieved. In this way, by piling up the triangle mesh along the heliocentric distance direction, a 3D grid system can be generated, which is unstructured in (θ, ϕ) directions and structured in the radial direction. In the computational domain of spherical shell, the "parallel" axis is selected to be the radial direction. This 3D grid system can be efficiently suitable for the finite volume method-like numerical schemes, such as those successfully used formerly (Feng et al. 2007; Hu et al. 2008; Nakamizo et al. 2009). The grid mesh generated by the polyhedron-splitting method is of sphere-surface body fitting that can enable us to easily implement the inner boundary conditions. But, its obvious disadvantage may be its difficulty of paralleling in "(θ, ϕ)" directions.

The general mesh generation technique (Hartmann et al. 2008) seems to be a powerful tool for fitting our spherical shell-shaped computational domain of solar wind modeling from the Sun to Earth. The advantages regarding mesh generation involve fully automated mesh generation based on the triangulated geometry, cell refinement and cell coarsening, automatic solution–AMR/coarsening at run time, and ideally shaped cells in the interior of the computational domain. However, the general mesh generation may seriously deform the cells in body-fitted grids resulting in loss of accuracy and numerical errors due to the quality of the triangulation. Typically, this method is also tedious, and is not able to accurately yield curved 3D boundaries. The grid system from this method may also meet a parallel problem in "(θ, ϕ)" directions. Automatic mesh generation, a relatively new field, seems promising, but at present it has no sound criteria for obtaining the high-quality mesh distribution. Thus, up to now the application of the general mesh generation has not been frequently seen in solar–terrestrial space modeling.

In order to overcome these problems mentioned above, overset grids (Rai 1986; Steger & Benek 1987; Chesshire & Henshaw 1990; Usmanov 1996; Yoshida & Kageyama 2004; Kageyama & Sato 2004; Kageyama 2005) may be one of the best choices. Among these overset grids, the Yin-Yang grid has been applied successfully in simulations of geoscience (Yoshida & Kageyama 2004; Kageyama & Sato 2004; Kageyama 2005). The cubed sphere (Ronchi et al. 1996), which constitutes a spherical surface with six non-overlapping components, is another overset grid recently used for the numerical study of mantle convection, thermal convection, and shallow water (Hernlund & Tackley 2003; Choblet 2005; Chen & Xiao 2008). Based on the concept of an overset grid (Chesshire & Henshaw 1990), Henshaw & Schwendeman (2008) proposed a numerical solution approach for nonlinear evolutionary partial differential equations in complex 3D domains, where the domains are represented by overlapping structured grids, and block-structured AMR is used to locally increase the spatial resolution. Usmanov (1996) also introduced a composite mesh in his 3D solar wind modeling. His composite mesh consists of three overlapping spherical meshes. The first one is the usual spherical mesh with a limited extension in latitude (42° ⩽ θ₁ ⩽ 138°, 0° ⩽ ϕ₁ ⩽ 360°). The polar axis of the mesh is directed along the solar rotation axis. Two other meshes are introduced to cover the polar regions in both hemispheres. These meshes are fragments of spherical coordinates (36° ⩽ θ_2,3 ⩽ 144°, 26° ⩽ ϕ₂ ⩽ 144°, and 216° ⩽ ϕ₃ ⩽ 324°) with the polar axis lying in the equatorial plane of the first coordinate system (θ₁ = 90°, ϕ₁ = 90°).

Here, motivated by the Yin-Yang grid (Kageyama & Sato 2004; Yoshida & Kageyama 2004; Kageyama 2005), the cubed sphere grid (Ronchi et al. 1996), and the successful application of the polyhedron-splitting grid (Feng et al. 2007; Hu et al. 2008; Nakamizo et al. 2009), a six-component grid system for solar wind MHD modeling is proposed for our three-dimensional Solar-InterPlanetary Conservation Element/Solution Element (SIP-CESE) MHD model (Feng et al. 2007, 2009; Hu et al. 2008; Zhou et al. 2008; Zhou & Feng 2008).

The paper is organized as follows. (1) The governing MHD equations with volumetric heating consideration are described. (2) We introduce the six-component grid system and transformation between two components. (3) We give some improvements for a SIP-CESE MHD model with a multiple time-stepping method, Courant-number insensitive method, and multigrid method for the Poisson solver to assure the soloidinal condition (i.e., ∇ · B) is satisfied. (4) Initial and boundary conditions are specified. (5) Numerical results for steady-state solar wind structure of Carrington rotation (CR) 1911 are displayed. (6) We present conclusions and discussions.

2.1. Governing Equations

The solar wind evolution is governed by the modified MHD equations. By splitting the magnetic field into a time-independent potential magnetic field B₀ and a deviation B₁ (Ogino & Walker 1984; Tanaka 1994), we can write the solar–interplanetary governing MHD system as

$$\begin{eqnarray} &&\frac{\partial \rho }{\partial t}+ \nabla \cdot \rho \mathbf {u} = 0\nonumber\\ &&\frac{\partial \rho \mathbf {u}}{\partial t}+ \nabla \cdot \bigg[\rho \mathbf {u} \mathbf {u} + \mathbf {I}\left(p + \frac{1}{2}\mathbf {B}_1^2 + \mathbf {B}_1 \cdot \mathbf {B}_0\right)\nonumber\\ &&\qquad-\mathbf {B}_1 \mathbf {B}_1-\mathbf {B}_1 \mathbf {B}_0 - \mathbf {B}_0 \mathbf {B}_1\bigg] \nonumber\\ &&\quad=\mathbf {j}_0 \times \mathbf {B}_0 + \rho \left[\mathbf {g}- {\bOmega } \times ({\bOmega } \times \mathbf {r})\right]- 2 \rho {\bOmega } \times \mathbf {u} + S_m, \nonumber\\ &&\frac{\partial e_1}{\partial t}+ \nabla \cdot \left[\mathbf {u} \left(e_1 + p + \frac{1}{2}\mathbf {B}_1^2 + \mathbf {B}_1 \cdot \mathbf {B}_0\right)-(\mathbf {u} \cdot \mathbf {B}_1)(\mathbf {B}_1+\mathbf {B}_0)\right] \nonumber\\ &&\quad={-}\mathbf {B}_1 \cdot \frac{\partial \mathbf {B}_0}{\partial t} + \mathbf {E} \cdot \mathbf {j}_0 + \rho \mathbf {u}\cdot [\mathbf {g} - {\bOmega } \times ({\bOmega } \times \mathbf {r})] + Q_e,\nonumber\\ &&\frac{\partial \mathbf {B}_1}{\partial t}+ \nabla \cdot (\mathbf {u}\mathbf {B}-\mathbf {B}\mathbf {u}) = -\frac{\partial \mathbf {B}_0}{\partial t}, \end{eqnarray} \tag{ 1 }$$

where

$$\begin{equation*} \mathbf {E}=\mathbf {u}\times \mathbf {B}, \ \ \mathbf {j}_0 = \nabla \times \mathbf {B}_0,\ \ {\rm and} \ \ e_1=\frac{1}{2} \rho \mathbf {u}^2+\frac{p}{\gamma -1}+\frac{1}{2} \mathbf {B}_1^2 \end{equation*}$$

the latter of which corresponds to the modified total energy density consisting of the kinetic, thermal, and magnetic energy density (written in terms of B₁). A factor of 1/$\sqrt{\mu }$ is included in the definition of B. Since B₀ is constant in time, many terms about B₀ on the right-hand side will vanish.

Here, ρ is the mass density, u = (u, v, w) are the velocities in the x-, y-, and z-directions, 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₁. t and r are time and position vectors originating at the center of the Sun, g = −GM/r³ · r is the solar gravitational force, ${\bOmega }$ is the solar angular speed, and γ is the ratio of specific heats. ρ, u, p, B, r, t, and g are normalized by the characteristic values $\rho _S, a_0, \rho _S a^2_0, \sqrt{\vphantom{A^A}\smash{\hbox{${\rho _S a^2_0}$}}}, R_S, R_S/a_0$, and a²₀/R_S, where ρ_S and a₀ are the density and sound speed at the solar surface. The solar rotation is considered in the present study with angular velocity $| {\bOmega }| = 2\pi /26$ rad day⁻¹ (here normalized by a₀/R_S). For γ, similar to that of Wu et al. (1999), a variable polytropic index is used:

$$\begin{equation*} \gamma =\left\lbrace \begin{array}{@{}ll} 1.05,& r/{R_S}\le 5\\ \ms 1.05+0.03 (r/{R_S}-5),& 5\le r/{R_S}\le 22\\ \ms 1.56,& r/{R_S}>22, \end{array} \right. \end{equation*}$$

S_m and Q_e stand for the momentum and energy source terms, which are responsible for acceleration and heating of the solar wind and will be described later.

The solar wind evolution is calculated in a reference frame of heliographic coordinates corotating with the Sun. For this reference coordinate system we use (r, θ, ϕ) for the position of a point in solar-interplanetary space and (x, y, z) is used to express its corresponding Cartesian coordinates. Sometimes, the analysis of computational results is carried out in the coordinate system at rest to remove these rotational effects.

2.2. Determination of the Momentum Source Term S_m and the Energy Source Term Q_e

Solar wind heating/acceleration has long been debated by scientists for long (Aschwanden 2004; Aschwanden et al. 2008). Since the creative work by Pneuman & Kopp (1971) empirically describing the source terms for the deposition of energy and/or momentum into the solar wind, many works have been devoted to the numerical study of producing solar steady-state conditions (e.g., Usmanov 1993; McKenzie et al. 1997; Mikić et al. 1999; Suess et al. 1996, 1999; Wang et al. 1998; Wu et al. 1999; Groth et al. 2000; Feng et al. 2005, 2007; Zhou et al. 2008; Kleimann et al. 2009; Nakamizo et al. 2009).

To mimic the effects of energy absorption above the transition region—thermal conduction, radiative losses, coronal heating, and solar wind acceleration—these models usually use spatial profiles as the source terms for the deposition of the energy or momentum by exponentials in radial distance, such as by an ad hoc volumetric heating function of exponential decaying form which originated with Hartle & Barnes (1970). The thermal conduction term of Spitzer type $\nabla (\xi T^{\frac{5}{2}} \frac{\nabla T\cdot {\mathbf B}}{B^2})\cdot {\mathbf B}$ for collisionless regimes has also been considered in the heating source term by some authors (e.g., Suess et al. 1996; Wang et al. 1998; Lionello et al 2001; Endeve et al. 2003; Li et al. 2005).

In some models, the energy and momentum exchange between the solar wind plasma and large-scale Alfvén wave (Jacques 1977; Dewar 1970; Barnes 1992; Usmanov et al. 2000; Usmanov & Goldstein 2003) has been considered for heating and accelerating the solar wind. For instance, Usmanov et al. (2000) and Usmanov & Goldstein (2003) used a short-wavelength Wentzel–Kramers–Brillouin approximation of the Alfvén wave turbulence to mimic the energy deposit for solar wind heating and acceleration.

The well-known Wang–Sheeley–Arge (WSA) model (Arge & Pizzo 2000; Arge et al. 2004), a relationship between the flux tube expansion factor f_S and satellite observations of the solar wind speed, has been used to drive the solar wind background simulation. The WSA model has been currently used in the kinematic HAFv.2 model (Fry et al. 2003; Smith et al. 2009) to obtain the solar wind background. Furthermore, Roussev et al. (2003) and Cohen et al. (2007) used the WSA model as an input to present an improved 3D MHD study for the solar wind, where the processes of turbulent heating in the solar wind are characterized by a varied polytropic index. Using this method, Cohen et al. (2007, 2008) are able to reproduce well the fast and slow plasma streams of the solar wind.

In order to obtain the distributions of the coronal density and temperature, by using the observed solar photospheric magnetic field and K-coronal brightness, the self-consistent structure of the distributions of the coronal density and temperature on the source surface at 2.5 R_S is obtained with the help of MHD equations (Wei et al. 2003). Using the self-consistent source surface structures as inputs, Feng et al. (2005) and Shen et al. (2007) develop a 3D MHD code (COIN-TVD MHD model) for the study of the ambient solar wind from the source surface to the Earth or beyond.

Although significant progress has been seen in this area, a sound mathematical description of the 3D topological structure of the interplanetary magnetic field (IMF) used in numerical modeling which affects the energy transfer process in the heliosphere does not exist.

In order to empirically specify the energy and momentum source term to take into account the magnetic field topology effects, we follow the work of Nakamizo et al. (2009) on the expansion factor "f_S" in the volumetric heating function and momentum source term. Thus, the volumetric heating function and momentum source term are

$$\begin{eqnarray} Q_e &=& Q_1 \exp \left({-}\frac{r}{ L_{Q_1}}\right) + Q_2 \left(\frac{r}{R_S}-1\right) \exp \left({-}\frac{r}{ L_{Q_2}}\right)\nonumber\\ &&+\nabla \left(\xi T^{\frac{5}{2}} \frac{\nabla T\cdot {\mathbf B}}{{\mathbf B}^2}\right)\cdot {\mathbf B},\nonumber\\ S_m&=&M\left(\frac{r}{R_S}-1\right)\exp \left({-}\frac{r}{ L_M}\right). \end{eqnarray} \tag{ 2 }$$

Here, r is the heliocentric distance, Q₁, Q₂ and $L_{Q_1}, L_{Q_2}$ are the intensity and decay height of heating, T is the temperature, and M and L_M are the intensity and the decay height of the momentum addition. A Spitzer type thermal conduction has been added in the third term of Q_e and we choose $\xi = 5 \varepsilon _p=1.6 \times 10^{-12} \rm {J} \rm {m}^{-1} \rm {s}^{-1} \hbox{K}^{-7/2}$ according to Endeve et al. (2003).

However, our approach is slightly different from that of Nakamizo et al. (2009). We consider the topology of the magnetic field expansion factor f_S and the minimum angular separation θ_b between an open field foot point and its nearest coronal hole boundary. Specifically, we assume that Q₂ = Q₀C_a, M = M₀C_a, where $C_a =C^{\prime }_a/\rm {max}(C^{\prime }_a)$ with $C^{\prime }_a= \frac{(5.8-1.6e^{[1-(\theta _b/8.5)^3]})^{3.5}}{{(1+f_S)^{2/7}}}$. Some empirical functions with various free parameters, similar to the above formula, have been employed to derive the solar wind speed near the Sun (Arge et al. 2004; Owens et al. 2005). Here, the constant values of Q₁, Q₀, and M₀ are 1.5 × 10⁻⁹ J m⁻³ s⁻¹, 1.18 × 10⁻⁷ J m⁻³ s⁻¹ and 7.9 × 10⁻¹⁴ Nm⁻³, respectively. $L_{Q_1}, L_{Q_2}$, and L_M are set to be 1 R_S. The expansion factor (Wang & Sheeley 1990) reads $f_S=(\frac{R_S}{r})^2\frac{B_{R_S}}{B_r}$ where R_S and r are the solar radius and the distance from the solar center, and $B_{R_S}$ and B_r are magnetic field strength at the solar surface and at r. The idea is motivated by the observation that there exists a close relationship between the solar wind velocity and the inverse of the expansion factor (Levine et al. 1977; Wang & Sheeley 1990) and the angular distance θ_b distinguishes the high-speed solar winds from the low-speed solar wind (high-speed stream originating from the center of an open field region has large θ_b and low-speed stream from the hole boundary has a small θ_b).

Here we use the potential field source surface (PFSS) model (Luhmann et al. 2002; Zhao et al. 2006; Hu et al. 2008) to determine the expansion factor f_S and θ_b and choose the source surface at 2.5 R_S. In the heating source term, we provide f_S and θ_b with additional meanings, that is, we allow them to have a 3D distribution from the solar surface at 1 R_S to the source surface at 2.5 R_S. For a given point O = (r, θ, ϕ) in the open field region from 1 R_S to the source surface at 2.5 R_S, the value of f_S at this point is defined to be $f_S(r, \theta, \phi)=(\frac{R_S}{R_{SS}})^2\frac{B_{R_S}(R_S, \theta _S, \phi _S)}{B_{R_{SS}(R_{SS}, \theta _{SS}, \phi _{SS})}}$ where R_S and R_SS are the solar radius and the source surface radius, and $B_{R_S}$ and $B_{R_{SS}}$ are radial magnetic field strength at the solar surface and at R_SS. Here, (θ_S, ϕ_S) and (θ_SS, ϕ_SS) are the photospheric coordinates and the source surface coordinates of the field line that passes through the point (r, θ, ϕ). During this procedure, the three points (R_S, θ_S, ϕ_S), (r, θ, ϕ), and (R_SS, θ_SS, ϕ_SS) are lying on the same field line. When O = (r, θ, ϕ) is located at the source surface, that is, r = 2.5R_S, f_S has its usual meaning as the expansion factor. By this definition, f_S has the same value along the magnetic field line passing through point O. The value of θ_b at this point is just the minimum angular separation between point O and its nearest open-closed field boundary on the surface r. Keeping these additional meanings in mind, we here have three-dimensional distribution of f_S and θ_b depending on (r, θ, ϕ) in the region 1R_S ⩽ r ⩽ 2.5R_S for our heating source terms. Beyond 2.5 R_S, f_S and θ_b take their distributions on the source surface.

3.1. Six-component Mesh Grid System in the Computational Domain of the Spherical Shell

We introduce a composite mesh that consists of six identical component meshes to envelope a spherical surface with partial overlap on their boundaries (Figure 1). Each component grid is a low-latitude spherical mesh, which is defined in the 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{ 3 }$$

where δ is proportionally dependent on the grid spacing entailed for the minimum overlapping area. Each component is confined by the same region as that in Equation (3) but in different spherical coordinates. We define the red-colored spherical mesh in Figure 1 as our reference heliographic or Cartesian coordinates. The relations between our reference heliographic coordinates C₁ and other coordinates C_j(j = 2, 3, 4, 5, 6) are denoted in the Cartesian coordinates by (x₂, y₂, z₂) = (y₁, − x₁, z₁), (x₃, y₃, z₃) = (−x₁, − y₁, z₁), (x₄, y₄, z₄) = (−y₁, x₁, z₁), (x₅, y₅, z₅) = (z₁, y₁, − x₁), (x₆, y₆, z₆) = (−z₁, y₁, x₁), where (x₁, y₁, z₁) are the Cartesian coordinates of the red-colored part associated with reference Cartesian coordinates; (x₂, y₂, z₂), (x₃, y₃, z₃), (x₄, y₄, z₄), (x₅, y₅, z₅), and (x₆, y₆, z₆) denote those of green-colored, blue-colored, cyan-colored, yellow-colored, and purple-colored parts, respectively, as labeled in Figure 1. From these Cartesian coordinate transforms, we can easily get their relationships in their associated spherical coordinates. For spatial discretization, we define the mesh point on each component C_ℓ(ℓ = 1, ..., 6) as

$$\begin{equation*} \theta _j^{\ell } = \theta _{\rm {min}} + j\Delta \theta,\quad j=0, 1, \ldots, N_\theta +1 \end{equation*}$$

$$\begin{equation*} \phi _k^{\ell } = \phi _{\rm {min}} + k \Delta \phi,\quad k=0, 1, \ldots, N_\phi +1 \end{equation*}$$

and

$$\begin{equation*} \Delta \theta = (\theta _{\rm {max}}-\theta _{\rm {min}})/(N_\theta -1) \end{equation*}$$

$$\begin{equation*} \Delta \phi = (\phi _{\rm {max}}-\phi _{\rm {min}})/(N_\phi -1), \end{equation*}$$

where N_θ and N_ϕ are the mesh numbers of the latitude and longitude, respectively. $\theta _{\rm {\rm {min}}} = \frac{\pi }{4},\theta _{\rm {max}} = \frac{3\pi }{4}, \phi _{\rm {min}} = \frac{3 \pi }{4}, \phi _{\rm {max}} = \frac{5 \pi }{4}$. For the three-dimensional six-component mesh, the construction is straightforward. The mesh grid in the r-direction as shown in Figure 2 is set up by stacking the basic six-component grids shown in Figure 1 in the r-direction by $\emph {r}(1)=1 R_S$, $\emph {r}({i+1})=\emph {r}(\emph {i})+\Delta \emph {r}(\emph {i})$, where i = 1, ..., N_r. The grid spacing on the spherical surface is quasi-uniform, and the ratio of the minimum grid spacing over maximum grid spacing is approximately 0.707. For the spherical surface, the overlapping region between component grids is at least two grid sizes, i.e., δ = Δθ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the present work, the parallel implementation in the whole computational domain of our simulated region is realized by domain decomposition of six-component grids based on the spherical surface and radial direction partition. In order to enhance the resolution in the "sphere" and radial direction we devise six parts among radial parallel decomposition: 1–3.5 R_S, 3.5–10 R_S, 10–25 R_S, 25–75 R_S, 75–166 R_S, and 166–247 R_S with the facility of multiple time stepping integration (to be given in Subsection 4.2 below), and the following grid partitions are employed: for 1–25 R_S, ${N}_\phi ={ N}_\theta = 2\times 2^5-1, \Delta \emph {r}(\emph {i})= 0.01 R_S$ if $\emph {r}(\emph {i}) < 1.1 R_S$; $\Delta \emph {r}(\emph {i}) = {\rm min}(A \times \lg (\emph {r}({i-1})), \Delta \theta \times \emph {r}({i-1}))$ with $A = 0.01/\lg (1.09)$ if $\emph {r}(\emph {i}) < 3.5 R_S$; and $\Delta \emph {r} (\emph {i})= \Delta \theta \times \emph {r}({i-1})$ if $\emph {r}(\emph {i}) > 3.5 R_S$. For 25–75 R_S, N_ϕ = N_θ = 3 × 2⁵ − 1 and $\Delta \emph {r}(\emph {i}) = \Delta \theta \times \emph {r}({i-1})$. For 75–166 R_S, N_ϕ = N_θ = 2 × 2⁶ − 1 and $\Delta \emph {r}(\emph {i}) = \Delta \theta \times \emph {r}({i-1})$. For 166–247 R_S, N_ϕ = N_θ = 3 × 2⁶ − 1 and $\Delta \emph {r} (\emph {i})= \Delta \theta \times \emph {r}({i-1})$. Hereafter, such defined subdomains will be called C^q_ℓ, where ℓ, q = 1, ..., 6. We set the innermost region at the solar surface, and the outermost region corresponds to the sphere with its radius 247 R_S. As a result, for 1–25 R_S, the angular resolution is about 14 in the latitudinal and longitudinal directions; for 25–75 R_S, the angular resolution is about 095 in the latitudinal and longitudinal directions; for 75–166 R_S, the angular resolution is about 071 in the latitudinal and longitudinal directions; for 166–247 R_S, the angular resolution is about 048 in the latitudinal and longitudinal directions, and a spatial resolution in the radial direction gradually varies from 0.01R_S at the inner boundary on the solar surface to 1.8R_S near 1 AU at 215R_S.

Up to this point we have finished the description of our mesh grid system in the computational domain of the spherical shell 0 ⩽ θ ⩽ π, 0 ⩽ ϕ ⩽ 2π and 1R_S ⩽ r ⩽ 215R_S or beyond for solar-interplanetary modeling.

It is evident that the six-component composite mesh has the following advantages: (1) the decomposition of the spherical shell owns six identical regions with the same metric and the component grids can be transformed into each other by coordinate transformation such that we can make efficient and concise programs, (2) each component grid is just a regular low-latitude part of the latitude–longitude grid such that various program resources in the spherical or the Cartesian coordinates grid can be used directly, (3) the boundary or internal border value in the overlapping area can easily be determined by an interpolation from its neighbor component grids according to the related geometrical positions of component grids, (4) the parallelization can be efficiently implemented both in the radial direction and "(θ, ϕ)" directions, and (5) the six-component composite mesh is of sphere-surface body-fitting that enables us to easily implement the inner boundary conditions at the solar surface prescribed for solar wind modeling. These properties have been shared by the Yin-Yang grid (Kageyama & Sato 2004) and the cubed sphere grid (Ronchi et al. 1996).

3.2. Vector Transformation Formulae Between Field Vectors Among Six-component Grids

A novel SIP-CESE MHD model has been developed and used for 3D solar wind modeling and solar transient events from the Sun to Earth (Feng et al. 2007; Hu et al. 2008; Zhou et al. 2008; Zhou & Feng 2008; Feng et al. 2009). Here, the SIP-CESE MHD model will be carried out on our newly introduced six-component grid system using governing Equations (1). Following the SIP-CESE MHD method given by Feng et al. (2007), solving the governing Equations (1) leads us to the discretized integral equation $\Phi ({\mathbf {U}})(\equiv {\mathbf {U}}-\frac{\Delta t}{2} {\mathbf {\eta } }({\mathbf {U}}))={\mathbf {U}}_H$ and its iteration procedure for the corresponding system of linear equations

$$\begin{equation} \frac{\partial \Phi ({\mathbf {U}})}{\partial {\mathbf {U}}^{(i)}}({\mathbf {U}}^{(i+1)}- {\mathbf {U}}^{(i)})=-(\Phi ({\mathbf {U}}^{(i)})-{\mathbf {U}}_H) \end{equation} \tag{ 4 }$$

for the unknown U⁽ⁱ⁺¹⁾ − U⁽ⁱ⁾, i = 0, 1, 2, .... Here, rather than actually computing the inverse of this 8 × 8 Jacobian matrix $\frac{\partial \Phi ({\mathbf {U}})}{\partial {\mathbf {U}}^{(i)}}$ like Equation (17) of Feng et al. (2007), one can directly solve Equation (4) by "intel MKL" with the Fortran library function⁴ according to William et al. (1997).

In code programming, we can directly code the governing equations in each component grid as they are formulated. Meanwhile, we can take the advantage of using various resources available, such as mathematical tools and program libraries in the Cartesian or spherical coordinates. For the parallel computing, it is natural to divide the computational task into six parts, that is, we can realize the parallel in the "(θ, ϕ)" directions. Because the component grids are identical, a balance of the computational loads between the processors can be perfectly achieved. When we have 6N processors, we decompose the radial direction of each component grid into N blocks for N processors each. As done before (Feng et al. 2007), the domain decomposition in the radial direction is simple and straightforward in each component grid since it is a rectangular box in the computational space of (r, θ, ϕ).

In the following subsections, we will give other improvements for the numerical aspects of the SIP-CESE MHD model.

4.1. Courant-number Insensitive Method

As demonstrated by Chang & Wang (2003), the local grid Courant numbers (CFL) ν, generally affect the accuracy of the numerical solution computed from the CESE scheme by using a globally fixed time marching step, and a huge disparity in ν across the non-uniform mesh makes the solutions highly dissipative in regions where ν ≪ 1. In solar wind simulations, the grid points are accumulated heavily near the Sun due to the spherical shell geometry of the computational domain, and the plasma flow and magnetic field vary over many orders of magnitude in the solar–terrestrial space. Such distributions of mesh grids and physical parameters lead to a huge disparity in ν across the mesh. Thus, making the numerical solution less insensitive to variation in ν is necessary to get rid of the excessive numerical dissipation caused by small ν and achieve a higher numerical accuracy. To this end, we employ a multi-dimensional Courant-number insensitive scheme (CNIS) for Euler solvers proposed by Chang & Wang (2003) to our model by re-designing the procedure for the calculation of the spatial derivatives U_mx, U_my, and U_mz as follows. As shown in Figure 3(a), assume that M_A1 is the centroid of octahedron QB₅B₆B₇B₈A₁, and P_A1 is a point on the line segment $\overline{M_{A1}{A_1^*}}$ defined by

$$\begin{equation*} x_{P_{A1}}= x_{M_{A1}} + \nu (x_{A1} - x_{M_{A1}}) \end{equation*}$$

$$\begin{equation*} y_{P_{A1}}= y_{M_{A1}} + \nu (y_{A1} - y_{M_{A1}}) \end{equation*}$$

$$\begin{equation*} z_{P_{A1}}= z_{M_{A1}} + \nu (z_{A1} - z_{M_{A1}}), \end{equation*}$$

where ν represents the grid CFL at Q*. The other five corresponding points P_A2, P_A3, P_A4, P_A5, and P_A6 can be similarly given. The values at P_A1, P_A2, P_A3, P_A4, P_A5, and P_A6 are calculated from the Taylor expansion at A*₁, A*₂, A*₃, A*₄, A*₅, and A*₆, respectively. In CNIS, the spatial derivatives U_mx, U_my, and U_mz at Q* are obtained with a weighted central difference type reconstruction approach by basing the values at P_A1, P_A2, P_A3, P_A4, P_A5, and P_A6, instead of the values at A*₁, A*₂, A*₃, A*₄, A*₅, and A*₆ as done before (Feng et al. 2007).

Zoom In Zoom Out Reset image size

From this procedure, it is easily seen that when ν → 0 (P_A approaches M_A), it will approach the original non-dissipative "a" scheme as shown by Chang & Wang (2003), while ν → 1 (P_A approaches A*), it switches to the so-called α scheme as we used before (Feng et al. 2007). By the introduction of the P_A point, the numerical dissipative generated by small local CFL can be effectively reduced, which is helpful for yielding accurate solutions.

For calculating the local CFL number ν, we use a simplified 3D grid CFL calculation as proposed by Yen et al. (2006). The simplified 3D grid CFL, according to Yen et al. (2006), is defined by enforcing the necessary condition of stability, namely, to require the analytical domain of dependence to be bounded by the numerical domain of dependence, associated with a solution point. Here, both domains of dependence are defined as follows: the numerical domain of dependence (Figure 3(b)) is the octahedron by the neighbor solution points A*₁, A*₂, A*₃, A*₄, A*₅, and A*₆. The analytical domain of dependence, as shown in Figure 3(c), is a sphere with a radius of V_fΔt/2 (V_f is the local speed of fast magnetosonic wave at Q*) and with a center at point P that is displaced by (−uΔt/2, − vΔt/2, − wΔt/2) from point Q*. To satisfy the necessary condition of stability, the sphere, the analytical domain of dependence, has to be bounded by the eight faces of the octahedron, the numerical domain of dependence. Referring to Figure 3(c), the stability requirement is to enforce the geometrical relation between the sphere and the triangle A*₃A*₄A*₆ by $\overline{Q^*T}+\overline{TS}<\overline{Q^*R}$ or $\nu _1 ={(\overline{Q^*T}+\overline{TS})}/{\overline{Q^*R}} <1$, where the line segments Q*R and PG are perpendicular to the triangle A*₃A*₄A*₆, the line segment Q*T is the projection of Q*P on Q*R, and the length of PG equals to the sphere radius V_f▵t/2, TS is its projection on Q*R. The corresponding ν₂, ν₃, ν₄, ν₅, ν₆, ν₇, and ν₈ can be similarly obtained by applying the similar geometrical requirement for the other seven faces of the octahedron. The local CFL ν at the solution point Q* is defined to be the maximum value among ν₁, ν₂, ν₃, ν₄, ν₅, ν₆, ν₇, and ν₈. It should be pointed out that as with done to the tetrahedron grid by Yen et al. (2006), we calculate the simplified grid CFL for only half time step marching instead of the formal grid CFL for a full Δt in order not to avoid incurring extra mesh management effort. The formal grid CFL for a full Δt involves the influencing cells for a solution point larger than the immediate neighboring cells, and the resulting grid CFL calculation is expensive and complex to apply for our large-scale parallel computation application of solar wind modeling with composite six-component grids.

4.2. Multiple Time Stepping

The plasma density, the Alfvén velocity, IMFs, the 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. Typically, the MHD time step is 1–3 s in the corona by the CFL criterion, and it will be 100–300 s in interplanetary space. If applying one uniform time step in the entire solar–terrestrial domain, time accuracy stability limitations on the time step in fine grid parts of the computational domain decrease the performance of numerical methods and it may become much more difficult to achieve efficient load balancing. The multi-time-stepping method (van der Ven et al. 1997; Maurits et al. 1998) is to take different time steps in different parts of the grid with large spatial grid difference, and it avoids the necessity of taking a single time step in the entire computational domain determined by the numerical CFL stability conditions (such as on the smallest grid cells). The application of the multi-time-stepping method can reduce the CFL number disparity that often produces excessive numerical dissipation and saves computation time.

Here, according to van der Ven et al. (1997) and Chang et al. (2005), our multiple time-stepping algorithm is implemented in six subdomains C^q_ℓ with the radial direction partitions: 1–3.5 R_S, 3.5–10 R_S, 10–25 R_S, 25–75 R_S, 75–166 R_S, and 166–247 R_S introduced before. First, we calculate the usual time step Δt_q for each subdomain C^q_ℓ by using the CFL stability condition with the formal Courant number 0.8. Then, the bulk time step Δt_q in the subdomain C^q_ℓ for q>1 is further constrained as follows: Δt_q+1 = M_qΔt_q, M_q = int(Δt_q+1/Δt_q) and q = 1, 2, 3, 4, 5, as shown in Figure 4.

Zoom In Zoom Out Reset image size

When the flow in C^q_ℓ travels M_q iterations of Δt_q, the flow in C^q+1_ℓ advances a single step of Δt_q+1. In this way, the flow in C^q_ℓ advances as many time steps of Δt_q as M_q × M_q+1 × ⋅ ⋅ ⋅ × M₅ for q ⩽ 5 in order to fit in the synchronization time step Δt_s defined as the maximum of these bulk time steps (Δt_s = Δt₆ in our present paper). After each synchronization time level, the bulk time step for each subdomain C^q_ℓ is determined again.

For the overlapping boundary interface between subdomain C^q_ℓ and subdomain C^q+1_ℓ, subdomain C^q_ℓ needs to advance M_q time steps of Δt_q and subdomain C^q+1_ℓ advances one time step of Δt_q+1 to arrive at the same time level. Hence, at the same time level the boundary values are obtained by simple cubic spline spatial interpolation, that is, the boundary values of C^q_ℓ are obtained by interpolating the associated C^q+1_ℓ values and vice versa. Otherwise, at every time step the qth subdomain's values at the boundary points (in subdomain C^q+1_ℓ) are determined by U_(q)(t + mΔt_q) = U_(q+1)(t) + mΔt_q(U_t)_(q+1)(t), where m = 1, ..., M_q − 1 and U_(q+1)(t), (U_t)_(q+1)(t) stand for the values of U and its time derivative in subdomain C^q+1_ℓ at time t. However, in order to keep the space-time conservation of fluxes, when the flow in C^q+1_ℓ is advanced we need to replace the (q + 1)th subdomain's radial flux flux_(q+1)(t)Δt_q+1 at its boundary points (in subdomain C^q_ℓ) by (flux_(q)(t) + flux_(q)(t + Δt_q) + ⋅ ⋅ ⋅ + flux_(q)(t + (M_q − 1)Δt_q))Δt_q, where flux_(q)(t + mΔt_q) is the radial flux in C^q_ℓ at time t + mΔt_q.

It should be noted that in our such grid architecture, in order to enforce flux conservation across an interface separating grid zones of different time step sizes, we follow the concept and methodology proposed by Chang et al. (2005). We have tested the code with dipole or multiple magnetic field as input to find that no spurious reflections are observed from the grid interface and that this procedure is much simpler and more efficient.

4.3. Multigrid Method for Cleaning Magnetic Field Divergence

The numerical error, caused by the nonzero divergence of the magnetic field ∇ · B, can influence the numerical stability and accuracy in multi-dimensional MHD simulations. In order to keep ∇ · B to an allowable numerical error, practioners of MHD simulations have devised some efficient methods, such as the constrained transport (CT) or CT/central difference (CD) approach (Evans & Hawley 1988; Balsara & Spicer 1999; Tóth 2000; Ziegler 2005; Londrillo & Del Zanna 2004; Stone & Gardiner 2009; Li 2008; Cunningham et al. 2009; Lee & Deane 2009), the reconstruction method (Balsara & Kim 2004; Londrillo & Del Zanna 2004; Li & Li 2004; Li 2008; Balsara 2009; Balsara et al. 2009), the eight-wave formulation (Powell et al. 1999), the mixed hyperbolic/parabolic correction procedure with the so-called generalized Lagrange multiplier method (Dedner et al. 2002; Feng et al. 2006; Fragile et al. 2005; Kleimann et al. 2009; Kataoka et al. 2009; Mignone et al. 2010), and the projection method (Brackbill & Baranes 1980; Tóth 2000; Balbás et al. 2004; Balsara & Kim 2004; Crockett et al. 2005; Tóth et al. 2006).

In this paper, we choose the projection method to eliminate ∇ · B, since it does not depend on specific differences and can achieve a better result compared to other methods (Tóth 2000; Balsara & Kim 2004). In this method, the B* field provided by the base scheme in time step n + 1 is projected to a divergence-free Bⁿ⁺¹ field by the following procedure. Since the magnetic field can be expressed by the sum of a curl and a gradient B* = ∇ × A + ∇Φ, we can arrive at a Poisson equation

$$\begin{equation} \nabla ^2\Phi =\nabla \cdot \mathbf {B}^*. \end{equation} \tag{ 5 }$$

With the help of solution Φ to Equation (5), the aim of the projection method is to correct the magnetic field to

$$\begin{equation} \mathbf {B}^{n+1}=\mathbf {B}^*-\nabla \Phi \end{equation} \tag{ 6 }$$

so as to force the numerical divergence of Bⁿ⁺¹ exactly zero. Obviously, this correction does not affect the current density J = ∇ × Bⁿ⁺¹ = ∇ × B*. In order to solve the Poisson Equation (5), Tóth (2000) used a stabilized biconjugate gradient iterative solver, while Tanaka (2003) suggested a conjugate residual method. Tóth (2000) also showed that its order of accuracy cannot be worse than that of his base scheme or a CT/CD type discretization therein. But, the usual method for the Poisson equation is very time consuming, especially for our large-scale simulation of solar–terrestrial space. Thus, a fast solver for the Poisson equation is very necessary. Multigrid methods (Holst & Saied 1993) are such highly efficient numerical techniques for this purpose and can reduce low-frequency error components by changing the grid size such that the global error disappears as a local error on coarse meshes, which is suitable for our projection procedure.

In this paper, we employ the approximate projection operator (Rider 1998) with the full multigrid (F-cycle) method (Holst & Saied 1993) described below to solve the Poisson Equation (5). Usually, the approximate projection operator does not depress high-frequency modes standing for a local decoupling of the magnetic field as occurred for the decoupling of the velocity field in numerical solutions of incompressible flows (Aprovitola & Denaro 2007; Rider 1998). For robust consideration, filtering is necessary, especially for long-term integrations. For this purpose, after carrying out the cleaning procedure through Equation (6), it is also advisable to have a diffusive term of the form (Dedner et al. 2002; Crockett et al. 2005)

$$\begin{equation*} \frac{\partial (\nabla \cdot \mathbf {B})}{\partial t}=\eta \nabla ^2(\nabla \cdot \mathbf {B}), \end{equation*}$$

which acts on the divergence of B(≡Bⁿ⁺¹). This may be described to eliminate a spatial derivative by pulling out a divergence operator:

$$\begin{equation} \frac{\partial \mathbf {B}}{\partial t}=\eta \nabla (\nabla \cdot \mathbf {B}), \end{equation} \tag{ 7 }$$

where $\eta =C \frac{(\Delta x)^2}{\Delta t}$, Δx is the minimal grid spacing, Δt is the time step, and C is a constant. A finite-difference approximation (Rider 1998; Miller & Colella 2001; Crockett et al. 2005) may be used to achieve a simple, single-step filter. In this paper, Δt_q is used in η for C^q_ℓ and the constant C is taken to be 0.05. The application of this filter (7) is twofold: it decreases the cell-centered divergence and effective damping of checkerboard modes (Dedner et al. 2002; Rider 1998; Miller & Colella 2001; Crockett et al. 2005).

Thus, the use of Equations (5), (6), and (7) will conclude our magnetic field divergence cleaning procedure used in the present paper. In what follows, we introduce our F-cycle method for solving the Poisson Equation (5).

In the present solar wind computation for solving Poisson Equation (5) we use the F-cycle iteration by multi-cubic interpolation for both restriction and prolongation between grids of different resolutions.

In Figure 5, numbers 1, 2, 3, 4, 5, and 6 stand for grid levels from coarse to fine grids for subdomain C¹_ℓ. Level 1 is the coarsest grid and Level 6 is the finest grid which is also our grid employed in the present paper. Similar grid level charts follow for C^q_ℓ. In this figure, a star "*" means prolongating the solution from the coarse grid to the fine grid, while a dot "·" denotes restricting the solution from the fine grid to the coarse grid. Figure 5 indicates that we carry out one F-cycle from Level 1 to Level 6 through the following prolongations and restrictions: 1 → 2 → 1 → 2 → 3 → 2 → 1 → 2 → 3 → ⋅ ⋅ ⋅ → 5 → 6. In each grid level, the Gauss–Seidel method is used as a smoother to reduce errors. In the program, two F-cycles have been run when physical time is less than 15 hr, three F-cycles have been performed from 15 to 25 hr, and four F-cycles have been performed after 25 hr. Along with the advancing time, the error of divergence of the magnetic field is accumulated, and more F-cycle runs have to be used in order to achieve a more accurate solution to the Poisson equation and to further reduce the divergence of the magnetic field, which will lead to more computational time.

Zoom In Zoom Out Reset image size

In order to facilitate the F-cycle iteration in our spherical shell computational domain of solar–terrestrial space, Equation (5) is solved on each part C^q_ℓ(ℓ, q = 1, ..., 6) (introduced in Section 3), conforming to the implementation of multiple time stepping. While solving Equation (5) by multigrid, we do not parallelize the multigrid method in each C^q_ℓ. If we use a multi-processor to deal with the F-cycle method, the size of the meshes in each processor at the coarser levels of the multigrid cycle is quite limited (typically 8 or 16 cells). In this case, most of the CPU time will be wasted in data exchanges between processors. Solving Equation (5) in each C^q_ℓ will give us a proper time-cost/performance. In the present work, at each new time step in C^q_ℓ, ∇ · B* is sent to a certain process (assigned to subdomain C^q_ℓ in our parallel implementation of the solar–terrestrial space) to finish the F-cycle iteration, then the obtained Φ through solving Equation (5) is dispatched to every process associated with C^q_ℓ, on which the cleaning procedure (6) is done and then the filter (7) is used.

It should be mentioned that the Neumann boundary condition is used only for the inner boundary at the solar surface, and for all other boundaries the Dirichlet boundary condition is employed. Meanwhile, our code for this F-cycle is modified from free source codes available online⁵. In fact, our cleaning procedure is performed for the time-dependent part of magnetic field B₁.

Since the modified total energy e₁ is an independent unknown in our CESE method, after the implementation of the projection, the only change may be the magnetic energy, which results in the variation or even the negative occurrence of the thermal energy or specifically the temperature. In the extreme case of plasma β ≪ 1, the occurrence of negative pressure may also appear. Thus, after the projection step, in order to avoid numerical difficulty with the possible occurrence of negative thermal pressure, the total energy in our simulation is adjusted as

$$\begin{equation} e_1^{n+1}=e_1^*+\big(\big|\mathbf {B}_1^{n+1}\big|^2-\big|\mathbf {B}_1^*\big|^2\big)/2 \end{equation} \tag{ 8 }$$

for the sake of numerical robustness, while the thermal energy and temperature obtained by the base CESE scheme are kept unchanged. The same procedure has been stated (Balsara & Spicer 1999; Cunningham et al. 2009) for their flux-CT algorithm to preserve thermal energy at the expense of energy conservation.

In order to assess the ∇ · B = 0 constraint numerically, as done in Equation (18) by Feng et al. (2007), we define the following error measurement as:

$$\begin{eqnarray} &&{\rm Error}=\frac{\sum _{k=1}^{M}\int _{V(24)}|(\nabla \cdot \mathbf {B})_k| dV}{\sum \nolimits _{k=1}^{M}\int _{V(24)} dV} \nonumber \end{eqnarray} \tag{ 9 }$$

where M is the number of mesh nodes in the computational domain and V(24), which is similar to V(18) in Feng et al. (2007), is a 24-faced polyhedron involved with the mesh grid in the three-dimensional spatial region.

For numerical simulation of solar-interplanetary space problems in the whole spherical shell between 1 R_S and 215 R_S or beyond, three types of boundary conditions need to be considered: horizontal boundary conditions at the six components' borders in overset parts, lower boundary conditions at r = 1R_S, and top boundary conditions at r = 215R_S or beyond.

The boundary prescription on the inner or lower boundary surface is a challenge: the lower boundary at solar surface 1 R_S is situated at the subsonic/sub-Alfvénic region, thus the numerical inconsistency will appear unless the proper treatment is provided there. The projected normal characteristic method (Nakagawa et al. 1987) imposes the boundary treatment satisfying both the underlying MHD equations and the prescribed boundary constraints. By combining the projected normal characteristic method and the solar wind mass flux limit derived by Neugebauer (1999) from Ulysses, a modified version of this method has been developed (Hayashi 2005; Hayashi et al. 2006). Such a time-dependent inner boundary condition is used here by limiting the mass flux escaping through the solar surface in order to obtain the mathematical and physical consistency for sub-Alfvénic fluid and produce more realistic plasma conditions in the coronal hole and streamer.

Since the top boundary at 215 R_S or beyond is in the supersonic/super-Alfvénic region, here we use a nonreflecting projected normal characteristic boundary condition that gives more stable simulation (Wu et al. 1996, 2006).

For the horizontal (θ, ϕ) boundary or the internal border between two components, according to the general overset technique (e.g., Chesshire & Henshaw 1990), interpolations are applied on the boundary of each component grid to set the boundary values or the internal boundary condition. That is, 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 relative position of the boundary point in the stencils. As Figure 6(a) shows, if a C₅'s boundary point Q(j⁵_hb, k⁵_hb) lies in a neighboring component C₁, we first determine its corresponding coordinates in the C₁'s grid, for example, (j¹₀, k¹₀), then obtain the values at Q(j¹₀, k¹₀) by means of interpolations of the known values at C₁'s grid points (j¹, k¹)'s neighborhood to Q(j¹₀, k¹₀), and finally transform the obtained values at Q(j¹₀, k¹₀) through the rotation matrix, for instance, $\mathbf {u}(Q(j_{hb}^5, k_{hb}^5))=\big(\begin{array}{@{}c c c@{} } \ssty 0& \ssty 0& \ssty 1 \\ [-5pt] \ssty 0& \ssty 1& \ssty 0 \\ [-5pt] \ssty -1& \ssty 0&\ssty 0 \end{array} \big) \mathbf {u} (Q(j_0^1, k_0^1))$ for vector field u, to get the corresponding value at the original C₅'s grid point Q(j⁵_hb, k⁵_hb). To be specific, in this paper, 16 grid points neighboring to Q(j¹₀, k¹₀) are used for our two-dimensional bi-cubic spline interpolation (Figure 6(b)). First, we obtain the values at the cross points A, B, C, and D along the j-direction by using a one-dimensional cubic spline interpolation, then we interpolate the Q(j¹₀, k¹₀) values along the k-direction with the computed values at A, B, C, and D. Finally, for the vector field we employ the rotation transform matrix to get the corresponding value at point Q(j⁵_hb, k⁵_hb). That is, for scalar variables, the interpolation is simply direct. However, for vector fields, we need the rotation transformation between different coordinates because the expressions of a vector in each component grid are different. The boundary points in other component grids could be achieved similarly. After the values at these boundary points are known, we can calculate the whole component region by the numerical procedure. As the numerical tests will later show, using this interpolation to communicate data across boundaries does not degrade the numerical accuracy of the whole model.

Zoom In Zoom Out Reset image size

In order to realistically generate the structured solar wind, we specify the inner initial boundary magnetic field with the observed line-of-sight photospheric magnetic field data. The observed photospheric magnetic field from the Wilcox Solar Observatory at Stanford University is used to deduce a 3D global potential magnetic field as initial magnetic input. Parker's solar wind flow (Parker 1963) provides the initial distributions of the plasma density ρ, gas pressure p, and the plasma velocity v. Here, the initial solar surface temperature and number density are set to be 1.3 × 10⁶ K and 1.5 × 10⁸ cm⁻³, respectively. Then, our code is run in time-relaxation manner until a steady-state equilibrium between flow and magnetic fields is achieved by satisfying some error criteria (Feng et al. 2007).

In this section, we present the 3D numerical results of structured solar wind from the Sun to Earth for CR 1911, which are obtained by executing the techniques introduced in the above sections. Our model is run on a 456 core cluster (19 nodes with each node having 24 cores run at 2.4 GHz). With the use of 156 MPI processes for CR 1911, it takes about 40 hr of wall time to obtain a steady-state solution at the physical time 250 hr.

Figures 7 shows the magnetic field lines, radial velocity, and number density on two different meridional planes at ϕ = 180° − 0° and ϕ = 270° − 90°. Figure 8 is the three-dimensional magnetic field distribution in the solar inner corona drawn by solid black lines while the arrowheads denote the direction of the magnetic field lines, and the color contours represent the number density distribution on the solar surface. The final relaxation solution has revealed the three-dimensional asymmetrical streamer resulting in an arcade structure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From Figures 7 and 8, we can see that a significant dipole component dominates the large-scale coronal magnetic field, which is a typical characteristic of the heliospheric magnetic field during periods of solar minimum. In these field configurations, magnetic field lines at high latitudes extend into interplanetary space to form coronal holes of high speed and low density. As a result, the plasma from both polar coronal holes occupies almost all of the heliospheric latitudes beyond 5 R_S at which the magnetic field is essentially open everywhere. While, at lower latitudes around the equator, a helmet streamer stretched by the solar wind can be observed within about three solar radii. Above the steamer, we can see a thin current sheet between different magnetic polarities. This scenario of the helmet streamer-current sheet system is consistent with those depicted by Pneuman & Kopp (1971) and Gosling et al. (1981).

Due to the proper inner boundary treatment of combining the projected normal characteristic method and the mass flux limit, the model can indeed generate reasonable contrasts of the plasma density (also temperature and velocity omitted here) between the coronal hole and streamer as shown in Figure 8. This point has been made by Hayashi (2005) and Hayashi et al. (2006). The fixed boundary condition will of course smear the density structure on the solar surface.

The coronal polarized brightness measurements are often used to diagnose the coronal density structure. The top panels in Figure 9 present the synthesized polarized brightness images from 2.3 R_S to 6 R_S projected on the meridional planes with ϕ = 180°–0° and ϕ = 270°–90°, respectively, which are computed from the 3D MHD simulation results. As we know, CR 1911 starts on June 28 and ends on 1996 July 24. If we assume that the coronal structure does not change significantly within one CR, which is a reasonable hypothesis during the period of solar minimum, then the coronal polarized brightness images observed on July 5 and June 28 roughly correspond to the top left and top right panels in Figure 9, respectively. Thus, for comparison, two LASCO/SOHO observed polarized brightness images are given in the bottom panels of Figure 9. From these figures, we can see that a qualitative agreement exists between the simulation results and the observations. The most highlighted feature in the right column is the relatively broad latitudinal extent of bright structures in the east limb. This feature results from the projection effect of the high-density regions between ϕ = 230° and ϕ = 280°, which can be seen clearly from the third subfigure in the left column of Figure 10 and from the synoptic polarized brightness images in both the east and the west limbs displayed in Figure 11. These bright structures are also the demonstration of the tilt of the magnetic field neutral line near this region. It should be noted that there is a small latitudinal difference between the bright structures obtained from simulations and observations. The difference results of the fact that the observations made by SOHO are about 7° north of the solar equator and the synthesized images are viewed from 1 AU in the solar equator. The other reason for this difference is due to the imperfect extrapolation of the polar photospheric magnetic field (Arge & Pizzo 2000). In addition, the distributions of high-density streamers in Figure 7 are consistent with the observation of LASCO C2 when considering the little difference mentioned above. In summary, we conclude that our simulation can reproduce the observed coronal density structure well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 10 displays the space weather background synoptic maps for the MHD steady-state solution on the different solar surfaces at 2.5 R_S, 20 R_S, and 215 R_S. This configuration is due to the interaction between the magnetic field and Parker's solar wind flow field. Quantitatively, N and T decrease with heliocentric distance while v_r increases. The flow in the polar regions is faster than that near the equatorial regions, and the heliospheric current sheet (HCS) region is surrounded by higher N. In the top panel of Figure 10, the neutral lines for the MHD solution show significant magnetic shear, and the magnetic field strength falls off along with heliocentric distance, but much less than that of the potential field. The strong warped structure in the streamer belt, located between 230° and 280° longitudes and encompassed by the low plasma pressure region with both low N and T, is spatially coincident with the extensions of the coronal hole boundaries (shown as θ_b distribution in Figure 12(a)). By examining the radial magnetic field, velocity, number density, and temperature distributions at 2.5 R_S and 20 R_S, we recognize that the speed and density patterns around the HCS and in the polar holes are in agreement with the synoptic maps of LASCO C2 in Figure 11. Also, the values of the number density and radial speed at 2.5 R_S and 20 R_S, respectively, are of the same magnitude as those shown in Figures 7 and 8 of Feng et al. (2007), Figure 5 of Hu et al. (2008), and Figures 1 and 2 of Wei et al. (2003). Additionally, the speeds at 2.5 R_S and 20 R_S are supported by the results derived from the observations by LASCO C2 and C3 (Wang et al. 1998; Porfir'eva et al. 2009), by interplanetary scintillation (Breen et al. 2002), and by the Ulysses solar corona experiment (Pätzold et al. 1997).

Zoom In Zoom Out Reset image size

From the distributions of plasmas and the magnetic field at 215 R_S in the right column of Figure 10, we can see that the bimodal structure of the solar wind (a combination of uniform, unipolar, and tenuous fast flows, with slow and relatively dense slow wind for the solar minimum period) and even the width of the slow wind band and the sharpness of velocity transitions are reproduced well in the model. The absence of a significant gradient of B_r with respect to heliolatitude in fast flows is also observed in this simulation.

By comparison with the four panels at 215 R_S in the right column of Figure 10, we would associate the fast solar speed with the presence of an interaction region (high temperature, high density, and intensive magnetic field). The interaction region is approximately aligned with the Parker spiral direction. From the warped structure of the HCS between 230° and 280° at 2.5 R_S and beyond, the features of corotating interaction regions (CIRs) can be seen clearly. At the southern flank of low-speed stream around the HCS, there is a compression region characteristic of steep gradients of density, temperature, and radial magnetic field, which corresponds to the fast solar wind overtaking the slow solar wind. At the northern flank there is a rarefaction region and the physical quantities vary smoothly. The rarefaction region is associated with the low-speed stream following the high-speed stream. The different features of the two regions can be seen more clearly from Figure 16. In addition, we can also see that the maximum latitudinal width of the slow wind flow region decreases approximately from 45° at 2.5 R_S to 25° at 215 R_S. This behavior is a consequence of magnetic pressure gradients.

In Figures 12(a) and (b) we compare the derived coronal hole boundaries mapped on the photosphere and the expansion factors "f_S" at 2.5 R_S with ϕ = 180° versus θ for the MHD model (solid line) and the PFSS model (dash line). Following Linker et al. (1999) and Riley et al. (2006), by identifying the closed-field regions and the open-field regions through tracing magnetic field lines on the lower boundary, a map of the coronal hole boundaries obtained by the MHD model and the PFSS model is presented in Figure 12(a). It is found that the coronal hole boundaries produced by the MHD model and the source-surface model are similar in this calculation. While slight differences exist, the two techniques have achieved the same qualitative characteristics for CR 1911. In particular, the polar coronal holes have the same topological configuration such as overall shape and span and have approximately the same area. In Figure 12(b), the heliolatitude changes of expansion factors "f_S" at 2.5 R_S with ϕ = 180° for the MHD model (solid line) and for the PFSS model (dash line) bear the same patterns, while the expansion factor from the MHD model is higher than that from the PFSS model. Figures 12(c) and (d) also display the distribution of the radial velocity v_r at 5 R_S from the MHD model and the WSA model, the latter of which is calculated by $v_r = 265+\frac{1.5}{(1+f_s)^{2/7}}(5.8-1.6 e^{[1-(\theta _b/7.5)^3]})^{3.5}$ with the help of θ_b and the expansion factor f_s (Owens et al. 2005). By examining this figure, we observe that the expansion factor of the MHD model is higher than that of the PFSS model. This indicates that the MHD model will give a large expansion of the magnetic field and thus present a relatively low speed, which is further confirmed by the radial speed variation in the bottom panel of Figure 12. This can be interpreted physically in terms of the pressure excised by the plasma that further spreads out the magnetic field in the MHD simulation.

Figure 13 displays two cases of the heating coefficient dependence on f_S only and a combination of θ_b and f_S. Figures 13(a) and (b) are the distributions of the corresponding heating coefficients C_a and C_b. They are associated with just our two cases of heating coefficients with θ_b and without θ_b if the exponential decaying factor is neglected from Q_e and S_m.

Zoom In Zoom Out Reset image size

Figures 13(c) and (d) show the MHD steady radial speeds when the heating coefficients C_a and C_b are used in the volumetric heating source terms. From this figure, it is easily seen that θ_b plays an important role in the heating process. C_a acts in a much wider open field region than C_b, and C_a is larger than C_b. Correspondingly, the distributions of the MHD steady radial speed v_r on the surface at 5 R_S with the heating coefficients given by C_a and C_b have some differences such as a narrower low speed region and a higher speed if θ_b is considered.

Figure 14 shows radial variations in the flow speed v_r and number density N from 1 R_S to 20 R_S versus heliocentric distance at ϕ = 0° with θ = 60° by the solid line and θ = 2° by the dashed line, which correspond to the open field region and closed field region. Evidently, we have a high-speed and low-density stream in the polar region (θ = 60°, ϕ = 0°). At 10 R_S and above, the flow is super–Alfvénic and supersonic nearly everywhere and changes very little. However, a low-speed and high-density flow is present in the current sheet region (θ = 2°, ϕ = 0°) and the velocity agrees with the profile of low solar wind speed (Sheeley et al. 1997). The overall distributions of velocity and density are similar to those obtained by Feng et al. (2007), but the fast stream here is much higher due to the addition of the heating and acceleration source terms.

Zoom In Zoom Out Reset image size

It is clear from Figure 14 that heating has much effect on the plasma distribution for an open field, but has little effect for a closed region. That is, by setting the energy and momentum intensity Q, M related to the expansion factor f_s and the minimum angular distance θ_b, the heating and acceleration mostly act on the coronal hole and exert little influence on the closed region of the streamer where the expansion factor goes to infinity.

Figure 15 for temporal evolution of error for ∇ · B defined in Subsection 4.3 is plotted in the region from 1 to 26 R_S, since the magnetic filed in this region is strong. From the evolution of |∇ · B| as shown in Figure 15, the use of the multigrid method can keep the divergence-free condition in an error much smaller than that without any special treatment (Feng et al. 2007). The estimation of error for |∇ · B| is about 10⁻⁶ and continues to be the same even after a long running time, and no obvious large error accumulation from |∇ · B| appears. In our numerical running, the error remains around 10⁻⁶ until the steady state is reached, although large-time marching sees a slight increase as shown by Figure 15(b). More F-cycles can greatly reduce this error but will cost much CPU time.

Zoom In Zoom Out Reset image size

Figure 16 shows the IMF lines of the steady state together with the distribution of the radial solar wind speed in the solar equatorial plane. The arrowheads denote the direction of the magnetic field line. The solar wind stretches the IMF outward into Archimedean spirals due to the solar rotation and the IMF freezing-in effect. The spirals coil more tightly in slow-speed flow than in high-speed flow. For instance, calculated from Figure 16, the angle between the radial direction and the IMF at 1 AU is 44° when the solar wind speed is 320 km s⁻¹ and 38° when the solar wind speed is 560 km s⁻¹. The conclusion that the IMF lines coil more tightly in the slow stream has also been pointed out by Odstrcil & Pizzo (1999).

Zoom In Zoom Out Reset image size

Figure 17 gives the two-state nature of a polar diagram of the computed solar wind speed at 1 AU as a function of heliolatitude. Flow speed profiles are displayed in the meridional plane (ϕ = 180°–0°). From it, we can clearly see a transition from one solar wind region to another at the midlatitudes and the width of transition is to some extent determined by the thickness of the HCS. This simulation result is quite consistent with Ulysses speed observations when normalized to 1 AU (McComas et al. 1998; Groth et al. 2000; Usmanov et al. 2000).

Zoom In Zoom Out Reset image size

In Figure 18, we make an analogous comparison with the 1 hr averaged WIND plasma data. WIND was located in the ecliptic plane near 1 AU (212 ∼ 213 R_S) and observed a more complex pattern of variations than those only partially reproduced by the MHD solution. In Figure 18, there is an overall good agreement between the simulations and the observations. The simulation reproduces well the observed relative high-speed wind part (>500 km s⁻¹), together with structure variations. From these comparisons, we can see that CIR-associated features, such as observed enhancements of density, temperature, and magnetic field magnitude in front of the high-speed wind or at the leading edge of the high-speed wind, are reproduced well. The rate of density enhancement and the slope of the rising solar wind speed in the simulation are closer to those of the 1 hr averaged data. However, some differences exist between the simulation results and the observational data. One difference is the arrival time of the fast wind. In the simulation, the fast wind is shifted late about 1–2 days. The second discrepancy is the structure of the solar wind speed in the last few days. The discrepancies in the last few days partially come from the uncertainties of the photospheric magnetic field measurements, the influence of differential rotation of the photosphere, and the inconsistency of the photospheric magnetic field at the beginning and ending days of a CR due to its evolution.

Zoom In Zoom Out Reset image size

By the way, we can clearly see that the solar wind speed obtained by the MHD simulation is closer to the observation than that derived by the Wang–Sheeley relationship v = 267.5 + 410.0/f^0.4_S (Arge & Pizzo 2000). However, the advantage of the MHD model over the WSA model is that the former can provide all the physical parameters everywhere in the computational domain, while the latter is used to predict the background solar wind speed and the IMF polarity at 1 AU (http://www.swpc.noaa.gov/ws/).

The 3D SIP-CESE MHD model (Feng et al. 2007; Hu et al. 2008; Zhou et al. 2008; Zhou & Feng 2008; Feng et al. 2009) has been greatly improved from the consideration of the following aspects: grid system, CNIS method, time integration, magnetic field divergence cleaning procedure, and volumetric heating. For the grid system, we introduced a composite mesh that consists of six identical component meshes to cover a spherical surface with partial overlap on their boundaries. Like the Yin-Yang grid (Kageyama & Sato 2004) and the cubed sphere grid (Ronchi et al. 1996), the important features of the six-component composite mesh are that these identical component grids with overlapping boundaries can be obtained from each other by coordinate transformation that makes the coding more efficient and concise, and each component with the quasi-uniform grid spacing is just a part of the latitude–longitude grid avoiding mesh convergence and singularities at the pole regions. Meanwhile, the grid system allows easy-paralleling not only in the "(θ, ϕ)" directions but also in the radial direction. The CNIS method can reduce high numerical dissipation in regions with small CFL numbers and thus keep the accuracy of the solution.

In time integration, the multiple time-stepping method is used by dividing solar–terrestrial space into six subdomains. This method enhances the convergence stability and speeds up the calculation by easing CFL number disparity due to spatial grid size variations and the different orders of magnitude of solar wind parameters such as plasma density, the Alfvén velocity, and IMFs from the Sun to Earth. In order to fix the magnetic field divergence error caused by the numerical scheme, the F-cycle iteration can maintain the ∇ · B error to an acceptable level of 10⁻⁶, and it has a proper time-cost/performance of about 40% CPU time in time step. The combination of CNIS and a multigrid Poisson solver is particularly efficient in enhancing solution resolution near the HCS. For the acceleration of solar wind, volumetric heating source terms are considered by using a 3D distribution profile based on the expansion factor f_S and the angular distance θ_b. Although our definition of f_S is slightly different from the original (Wang & Sheeley 1990; Arge & Pizzo 2000; Arge et al. 2004), its role has the same functions as claimed by Nakamizo et al. (2009); that is, this heating is weakened at the location where the magnetic field is in overradial expansion (f_S>1.0), and also it is strengthened at the location where the magnetic field is in underradial expansion (f_S < 1.0).

The use of angular distance θ_b, which can realistically specify solar wind speed at the boundaries and the interiors of coronal holes, can effectively distinguish the high-speed solar wind from the low-speed solar wind (high-speed wind emanating from the center of a coronal hole has large θ_b and low-speed wind from the hole boundary has a small θ_b). Our heating profile configuration can better lead to a non-uniform heating distribution reflecting the topology of the magnetic field in a key region of the solar wind generation. Meanwhile, our numerical results show that the Wang–Sheeley relationship yields higher predicted wind speed values at the source surface than the MHD model does. Also, the MHD model predicts speeds closer to the observation than those predicted by the Wang–Sheeley relationship.

Overall, our model can produce all the physical parameters everywhere within the computation domain. This may imply that the solar surface global heliospheric structure connection can be predicted by the simulation up to a generally acceptable level, although many unsatisfactory points remain as discussed here.

Differences occur due to many factors such as the uncertainties in photospheric magnetic measurements (especially in polar regions), the imperfection of the potential magnetic field approximation, the occurrence of some coronal mass ejections (CMEs) during CR 1911 (1996 June 28–July 25),⁶ the shortage of more sound physical-based and observation-supported heating mechanism, and the neglect of interaction between solar wind and interplanetary discontinuities. Although our heating consideration gives favorable numerical results, we have reasons to believe that volumetric heating alone 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.

Incidentally, our model validations for some other CRs show that the numerical results for the CRs during solar minimum are usually better than those for the CRs during solar maximum, when compared with observations in the solar corona and at 1 AU. This may be due to the fact that the currently used initial input of a magnetic field based on the potential field model cannot take into account the observed heliospheric open flux from active regions and CMEs by way of interchange reconnection frequently occurring at solar maximum (Cohen et al. 2007).

More available spacecraft observations will of course equip us with new findings about the coronal heating mechanism, the causes and mechanisms of CMEs' initiation, CMEs' 3D structure, and the interplanetary evolution process. The recently launched Solar Dynamic Observatory (SDO) will help us understand the Sun's magnetic changes. SDO will determine how the magnetic field is generated and structured, and how the stored magnetic energy is released into the heliosphere and geospace. STEREO observations can provide new insights into the 3D structure of CMEs and their evolution in the heliosphere which can directly be compared with existing models and simulations. Comprehensive data and analysis with multiple spacecraft (such as SDO, STEREO, SOHO, ACE, WIND, or other future missions) will probably help us develop the ability of including physically realistic coronal heating modules into 3D MHD codes, improve the determination of the structure of the ambient solar wind, and further numerically characterize the 3D propagation of CMEs through the heliosphere. Other aspects for space weather event simulations in 3D MHD from the Sun to Earth can follow suggestions made by Dryer (1998, 2007) and Wu et al. (2006). These will be our avenue to future improvements. In particular, this model is being used for the propagation study of transient events from the Sun to Earth.

The work is jointly supported by the National Natural Science Foundation of China (41031066, 40921063, 40874091, 40890162, 40904050, 40874077, and 40536029), the 973 project under grant 2006CB806304, and the Specialized Research Fund for State Key Laboratories. Dr. S. T. Wu is supported by AFOSR (grant FA9550-07-1-0468), AURA Sub-Award C10569A of NSO's Cooperative Agreement AST 0132798, and NSF (grant ATM-0754378). We thank the SOHO/LASCO team for letting us use their data. SOHO is a mission of international collaboration between ESA and NASA. The Wilcox Solar Observatory data used in this study were obtained via the Web site http://wso.stanford.edu at 2009:05:06_17:15:30 PDT courtesy of J.T. Hoeksema. The Wilcox Solar Observatory is currently supported by NASA. We also thank OmniWeb, http://omniweb.gsfc.nasa.gov/, from which we downloaded the hourly average solar wind data by WIND. The numerical calculation was completed on our SIGMA Cluster computing system. Special thanks go to our anonymous reviewer for helpful suggestions that contributed to the improvement of the paper.