A Data-Driven MHD Model of the Global Solar Corona within Multi-Scale Fluid-Kinetic Simulation Suite (MS-FLUKSS)

TL;DR

This paper introduces a data-driven MHD model of the global solar corona integrated within the MS-FLUKSS simulation suite, capable of utilizing diverse observational data and employing adaptive mesh refinement for improved accuracy.

Contribution

The paper presents a novel data-driven MHD model with characteristic boundary conditions, implemented in MS-FLUKSS, capable of assimilating various solar observational data for the first time.

Findings

Successful benchmark test of dipole field relaxation

Effective relaxation of PFSS field using HMI LOS data

Demonstrated model's capability with different observational inputs

Abstract

We have developed a data-driven magnetohydrodynamic (MHD) model of the global solar corona which uses characteristically-consistent boundary conditions (BCs) at the inner boundary. Our global solar corona model can be driven by different observational data including Solar Dynamics Observatory/Helioseismic and Magnetic Imager (SDO/HMI) synoptic vector magnetograms together with the horizontal velocity data in the photosphere obtained by the time-distance helioseismology method, and the line-of-sight (LOS) magnetogram data obtained by HMI, Solar and Heliospheric Observatory/Michelson Doppler Imager (SOHO/MDI), National Solar Observatory/Global Oscillation Network Group (NSO/GONG) and Wilcox Solar Observatory (WSO). We implemented our model in the Multi-Scale Fluid-Kinetic Simulation Suite (MS-FLUKSS) - a suite of adaptive mesh refinement (AMR) codes built upon the Chombo AMR framework developed at the Lawrence Berkeley National Laboratory. We present an overview of our model, characteristic BCs, and two results we obtained using our model: A benchmark test of relaxation of a dipole field using characteristic BCs, and relaxation of an initial PFSS field driven by HMI LOS magnetogram data, and horizontal velocity data obtained by the time-distance helioseismology method using a set of non-characteristic BCs.

Tables

Table 1. Table 1: Different ways to specify the characteristic BCs. The slow magnetosonic, Alfvénic, and fast magnetosonic velocities are denoted as a s subscript 𝑎 s a_{\mathrm{s}} , a A subscript 𝑎 A a_{\mathrm{A}} , and a f subscript 𝑎 f a_{\mathrm{f}} , respectively.

Case	# of physical BCs	# of numerical BCs	Variables to be imposed as physical BCs	Variables to be imposed as numerical BCs
Subslow inflow (	5	3	$\rho$, $B_r$, $v_{\theta }$, $v_{\varphi }$, and $B_{\theta }$	$v_r$, p, and $B_{\varphi }$
Sub-Alfvénic, but superslow inflow ()	6	2	$\rho$, $v_{\theta }$, $v_{\varphi }$, p, $B_r$, and $B_{\theta }$	$v_r$, and $B_{\varphi }$
Subfast, but super-Alfvénic inflow ()	7	1	$\rho$, $v_{\theta }$, $v_{\varphi }$, p, $B_r$, $B_{\theta }$, and $B_{\varphi }$	$v_r$

Equations

$$\frac{\partial}{\partial t}\left(\begin{array}[]{c}\rho\\ \rho\mathbf{v}\\ \mathbf{B}\\ E\end{array}\right)+\mathbf{\nabla}\cdot\left(\begin{array}[]{c}\rho\mathbf{v}\\ \rho\mathbf{v}\mathbf{v}+\mathbf{I}(p+\frac{B^{2}}{8\pi})-\frac{\mathbf{B}\mathbf{B}}{4\pi}\\ \mathbf{v}\mathbf{B}-\mathbf{B}\mathbf{v}\\ (E+p+\frac{B^{2}}{8\pi})\mathbf{v}-\frac{\mathbf{B}}{4\pi}(\mathbf{v}\cdot\mathbf{B})\end{array}\right)=\left(\begin{array}[]{c}0\\ \rho[\mathbf{g}+(\mathbf{\Omega}\times\mathbf{r})\times\mathbf{\Omega}+2(\mathbf{v}\times\mathbf{\Omega})]+\mathbf{{S}_{M}}\\ 0\\ \rho\mathbf{v}\cdot[\mathbf{g}+(\mathbf{\Omega}\times\mathbf{r})\times\mathbf{\Omega}]+S_{E}\end{array}\right),$$

$$S_{E}=\frac{Q_{0}}{f_{s}}\textrm{exp}\Big{(}-\frac{r}{L_{Q}}\Big{)},\\$$

$$\mathbf{S_{M}}=\frac{M_{0}}{f_{s}}\Big{(}\frac{r}{R_{\odot}}-1\Big{)}\textrm{exp}\Big{(}-\frac{r}{L_{M}}\Big{)},$$

$$f_s=\frac{B(R_{\odot })}{B(R_{SS})}{\left(\frac{R_{\odot }}{R_{SS}}\right)}^2.$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=0$$

$$\frac{\partial \mathbf{v}}{\partial t}+(\mathbf{v}\cdot \nabla )\mathbf{v}=-\frac{\nabla p}{\rho }-\frac{\mathbf{B}\times (\nabla \times \mathbf{B})}{4\pi \rho }+\mathbf{g}+{\mathbf{S}}_{\mathrm{M}}$$

$$\frac{\mathrm{d}p}{\mathrm{d}t}=a^2\frac{\mathrm{d}\rho }{\mathrm{d}t}+S_{\mathrm{E}}$$

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times (\mathbf{v}\times \mathbf{B})$$

$$\nabla \cdot \mathbf{B}=0.$$

$$\frac{\partial \mathbf{u}}{\partial t}+{\bar{A}}_r\frac{\partial \mathbf{u}}{\partial r}+{\bar{A}}_{\theta }\frac{\partial \mathbf{u}}{\partial \theta }+{\bar{A}}_{\varphi }\frac{\partial \mathbf{u}}{\partial \varphi }=\mathbf{S},$$

$$\frac{\partial \mathbf{u}}{\partial t}+{\bar{A}}_r\frac{\partial \mathbf{u}}{\partial r}+\mathbf{H}=0.$$

$$\bar{L}\frac{\partial \mathbf{u}}{\partial t}+\bar{L}{\bar{A}}_r\frac{\partial \mathbf{u}}{\partial r}+\bar{L}\mathbf{H}=0,$$

$${\mathbf{l}}_i\cdot \frac{\partial \mathbf{u}}{\partial t}+{\lambda }_i{\mathbf{l}}_i\cdot \frac{\partial \mathbf{u}}{\partial r}+{\mathbf{l}}_i\cdot \mathbf{H}=0,$$

$$\frac{\partial \mathbf{W}}{\partial t}+\Lambda \frac{\partial \mathbf{W}}{\partial r}+\bar{L}\mathbf{H}=0.$$

$$\frac{\partial}{\partial t}\left|\begin{array}[]{c}\mathbf{W}^{P}\\ \mathbf{W}^{N}\end{array}\right|+\left|\begin{array}[]{c}\bar{L}^{P}\\ \bar{L}^{N}\end{array}\right|\bar{A}_{r}\frac{\partial}{\partial r}\left|\begin{array}[]{c}\mathbf{u}^{P}\\ \mathbf{u}^{N}\end{array}\right|+\bar{L}\mathbf{H}=0,$$

$$\frac{\partial {\mathbf{W}}^P}{\partial t}=({\bar{L}}^P{)}^{-1}{\mathbf{C}}^P(t)\text{ or }\frac{\partial {\mathbf{u}}^P}{\partial t}={\mathbf{C}}^P(t),$$

$$\left|\begin{array}[]{c}\bar{I}^{P}\\ \bar{L}^{N}\end{array}\right|\frac{\partial}{\partial t}\left|\begin{array}[]{c}\mathbf{u}^{P}\\ \mathbf{u}^{N}\end{array}\right|+\left|\begin{array}[]{c}0\\ \bar{L}^{N}\end{array}\right|\bar{A}_{r}\frac{\partial}{\partial r}\left|\begin{array}[]{c}\mathbf{u}^{P}\\ \mathbf{u}^{N}\end{array}\right|+\left|\begin{array}[]{c}-\mathbf{C}\\ \bar{L}^{N}\mathbf{H}\end{array}\right|=0,$$

$$\frac{\partial}{\partial t}\left|\begin{array}[]{c}\mathbf{u}^{P}\\ \mathbf{u}^{N}\end{array}\right|=\mathbf{C}-(\bar{L}_{1})^{-1}\bar{L}_{2}\bar{A}_{r}\frac{\partial}{\partial r}\left|\begin{array}[]{c}\mathbf{u}^{P}\\ \mathbf{u}^{N}\end{array}\right|-(\bar{L}_{1})^{-1}\bar{L}_{2}\mathbf{H}.$$

Full text

1 Introduction

The solar wind (SW) emerging from the Sun is the main driving mechanism of solar events which may lead to geomagnetic storms which are the primary causes of space weather disturbances that affect the magnetic environment of Earth and may have hazardous effects on the space-borne and ground-based technological systems as well as human health. Therefore, accurate modeling of the SW is very important to understand the underlying mechanisms of such storms.

As the number of space-borne and ground-based observatories increased significantly over the past decades, a number of data-driven models of SW that involve formulations of time-dependent boundary conditions (BCs) that incorporate remote and in situ observations of the Sun in a self-consistent way have been developed (see, e.g.  [1, 2, 3, 4, 5, 6]).

The SW emerges from the Sun. In order to model the background SW in a physically-consistent way, it is essential to apply a proper boundary treatment based on observational data on the Sun. For this purpose, we have developed a data-driven magnetohydrodynamic (MHD) model of the global solar corona which applies characteristically-consistent BCs at the inner boundary of our computational domain located in the lower solar corona. Our characteristic BC formulation follows [7], is entirely based on observations, and does not rely on any non-reflection principle, as in [8]. It allows us to specify sufficient number of mathematically admissible BCs. The implementation of characteristic BCs at the lower corona is based not only on the knowledge of how many quantities should be specified, but also on the requirement that time-increments of quantities that are specified as physical BCs should be uniquely expressible in terms of the increments of proper characteristic variables.

With the availability of vector magnetograms and surface velocity components, we can impose our characteristic BC formulation to obtain a realistic time-dependent background SW model. Our model is driven by Solar Dynamics Observatory/Helioseismic and Magnetic Imager (SDO/HMI) synoptic vector magnetograms together with the horizontal velocity data in the photosphere obtained, e.g., with the time-distance helioseismology method [9]. In addition, our model can also be driven by other observational data including the line-of-sight (LOS) magnetogram data obtained by HMI, Solar and Heliospheric Observatory/Michelson Doppler Imager (SOHO/MDI), National Solar Observatory/Global Oscillation Network Group (NSO/GONG), and Wilcox Solar Observatory (WSO).

We obtain the initial distribution of plasma properties in the computational domain from the Potential Field Source Surface (PFSS) model [10, 11] for the magnetic field and Parker’s 1D isothermal SW model [12] for the hydrodynamic variables. We model the SW acceleration by introducing volumetric heating source terms into the momentum and energy equations that incorporate expansion factors to take the coronal magnetic field distribution into account [13, 14].

We have implemented our model in the Multi-Scale Fluid-Kinetic Simulation Suite (MS-FLUKSS) [15] – a suite of adaptive mesh refinement (AMR) codes designed to solve the coupled systems of MHD, gas dynamics Euler, and kinetic Boltzmann equations [16, 17, 18, 19]. MS-FLUKSS is built upon the Chombo AMR framework developed at the Lawrence Berkeley National Laboratory. MS-FLUKSS also has modules that treat pickup ions as a separate fluid and turbulence beyond the Alfvénic surface.

In section 2, we will present an overview of our global solar corona model with particular emphasis given to our characteristic BC formulation and SDO/HMI vector magnetogram and horizontal velocity data. In section 3, we will present results we obtained using our model. Finally, in section 4, we will present our conclusions.

2 Data-driven MHD model of global solar corona

2.1 Governing equations

We solve the set of ideal MHD equations with volumetric heating source terms to model acceleration of SW [13, 14]. They are written in terms of conservative variables, in conservation-law form as follows:

[TABLE]

where $\rho$, $\mathbf{v}$, $\mathbf{B}$, $p$, $E$, and $\mathbf{g}$ are the density, velocity, magnetic field, thermal pressure, specific total energy of the plasma, and gravitational acceleration, respectively. The source terms in the momentum and energy conservation equations include the Coriolis and centrifugal forces which are present when the system is solved in a frame corotating with the Sun. Accordingly, $\Omega$ and $\mathbf{r}$ correspond to the angular velocity of the Sun and position vector, respectively.

In order to model the acceleration of the SW, we introduce a volumetric heating source term, $S_{E}$, into the energy conservation equation, and the corresponding source term, $\mathbf{{S}_{M}}$, into the conservation of momentum equations [13, 14]. They are given as follows:

[TABLE]

[TABLE]

where $L_{\mathrm{M}}$, $L_{\mathrm{Q}}$, $M_{0}$, and $Q_{0}$ are the model constants given as $L_{\mathrm{M}}=L_{Q}=0.9R_{\odot}$, $M_{0}=2.65\times 10^{-14}$ N $\mathrm{m}^{-3}$, and $Q_{0}=1.65\times 10^{-6}$ J \mathrm{m}^{-3}$$\mathrm{s}^{-1}. Additionally, $f_{s}$ is the expansion factor by which a magnetic flux tube expands in solid angle between its footpoint location on the photosphere and the source surface which is typically at $R_{\mathrm{SS}}=2.5R_{\odot}$ [20]:

[TABLE]

These source terms take the coronal magnetic field topology into account by incorporating the expansion factor.

The expansion factor is initially computed in every cell located between the inner boundary and the source surface and kept constant throughout a simulation.

2.2 Observational data

Our model is primarily driven by full disk HMI synoptic vector magnetograms [21] and horizontal velocity data in the photosphere.

HMI [22] provides full disk vector magnetogram data with high cadence ($\sim$ 720 s) and no data gaps [23]. From these data, synoptic vector magnetograms are computed for each Carrington rotation for the entire SDO [24] mission. They have a resolution of 3600$\times$1440 pixels and are provided by the Joint Science Operations Center (JSOC) under the data series name hmi.B_synoptic [21].

The horizontal velocity data are inferred by the time-distance helioseismology method [9]. The acoustic travel times are measured using HMI Dopplergram observations under the data series name hmi.V and this method infers solar interior properties by inverting these measurements. Thus, near-real time full disk maps of subsurface wave-speed perturbations and horizontal flow velocities are produced for depths ranging from 0 to 20 Mm, every 8 hours. Carrington synoptic maps for the subsurface properties are made from these full disk maps. These maps have a resolution of 3000$\times$1000$\times$6 pixels where the last dimension corresponds to the depth beneath the solar surface. We utilize the corresponding data on the photosphere.

In Figures 1 and 2, three components of the magnetic field vector and horizontal velocity components are shown on the solar surface corresponding to Carrington Rotation (CR) 2145, respectively. The horizontal velocity profiles are given within the latitude range of $\pm 60^{\circ}$. Above this latitude, the inferred velocity is not very reliable.

Our model can also be driven by other observational data including NSO/GONG, HMI, SOHO/MDI, and WSO LOS magnetogram data.

2.3 Initial solution

We calculate the initial distribution of magnetic field using the PFSS model [10, 11] based either on the spherical harmonics approach [25, 26, 27] or a finite difference method by incorporating the solution provided by the Finite Difference Iterative Potential-field Solver (FDIPS) code [28] into MS-FLUKSS. FDIPS is a good alternative to the spherical harmonics approach. As the number of spherical harmonics increases, using raw magnetogram data specified on a grid that is uniform in the sine of latitude can result in inaccurate and unreliable results, especially in the polar regions close to the Sun. As the FDIPS code accepts any magnetogram data represented in the Flexible Image Transport System (FITS) file format as input, we can also incorporate potential magnetic field solutions for all such data into MS-FLUKSS. For the rest of the plasma parameters, we compute the initial solution from Parker’s isothermal SW model [12].

2.4 Characteristic boundary conditions

One of the important components of our data-driven MHD model is a genuinely characteristic BC formulation, which became possible only recently because of the availability of vector magnetograms. If the inner simulation boundary is in the vicinity of the photosphere ($R\approx R_{\odot}$), the radial component of the SW velocity is likely to be less than the Alfvén or even slow magnetosonic velocity. Not only this determines the number of physical BCs at $R\approx R_{\odot}$, but the characteristic variables to be specified are also uniquely defined. Since there are no analytic expressions for characteristic variables in MHD, they cannot be found from observational data. On the other hand, any admissible set of physical BCs should be resolvable for the above-mentioned physical characteristic variables. The remaining BCs should be mathematical (numerical). They should be specified, e.g., by extrapolation of the corresponding variables from inside the computational region. This approach is widely applied in aerodynamic simulations [7].

Consider the system of governing equations written in terms of the so-called primitive variables (the abbreviations are conventional).

[TABLE]

Assuming spherical coordinates for definiteness, we rewrite this system in the quasilinear form as follows:

[TABLE]

where $\mathbf{u}$ is the vector of primitive variables, $\bar{A}_{r}$, $\bar{A}_{\theta}$, and $\bar{A}_{\phi}$ are the coefficient matrices in front of the $r$-, $\theta$-, and $\phi$-derivatives, respectively, and $\mathbf{S}$ is the source term vector involving the gravitational force, spherical geometrical factors, and volumetric heating terms.

Let $\mathbf{H}=\bar{A}_{\theta}\frac{\partial\mathbf{u}}{\partial\theta}+\bar{A}_{\phi}\frac{\partial\mathbf{u}}{\partial\phi}-\mathbf{S}$. Then, Eq. (6) becomes

[TABLE]

To write out a set of compatibility relations, we multiply the above equation by the left eigenvector matrix of $\bar{A}_{r}$ [8, 29]:

[TABLE]

or in the component form (note that $\mathbf{l}_{i}\bar{A}_{r}=\lambda_{i}\mathbf{l}_{i}$)

[TABLE]

where $\mathbf{l}_{i}$ and $\lambda_{i}$ are the rows of the left eigenvector matrix $\bar{L}$ and eigenvalues of $\bar{A}_{r}$, respectively.

On introducing a set of characteristic variables $\mathbf{W}$, such that $\delta\mathbf{W}=\bar{L}\delta\mathbf{u}$, we obtain a diagonal system ($\Lambda=\{\lambda_{1},\lambda_{2},\ldots,\lambda_{8}\}$)

[TABLE]

No analytic expressions for $\mathbf{W}$ exist in MHD, but their increments are well defined.

At the inner (SW entrance) boundary, $W_{i}$ ($i=1,\ldots k$) corresponding to $\lambda_{i}>0$ (incoming characteristics) should be specified as physical BCs, whereas the BCs for remaining characteristic variables should be numerical (i.e., they should be derived from the characteristic compatibility relations for ${W}_{j}$ with $\lambda_{j}<0$ for $j=k+1,\ldots n$). For MHD equations, $n=8$.

Following this mathematical principle, we split the compatibility relations (i.e., Eq. (10)) into two sets corresponding to the physical (P) and numerical (N) BCs:

[TABLE]

where $\mathbf{W}^{P}$, $\mathbf{W}^{N}$, $\bar{L}^{P}$, and $\bar{L}^{N}$ are comprised of the characteristic variables and left eigenvectors corresponding to incoming and outgoing characteristics, respectively, and $\mathbf{u}^{P}$ and $\mathbf{u}^{N}$ are the primitive variables to be physically imposed or numerically determined.

The physical BCs are implemented by setting

[TABLE]

where $\mathbf{C}^{P}(t)$ is the vector of time dependent variations of the primitive variables that are physically imposed. At each time step, $\mathbf{u}^{P}$ values are interpolated in time using successive observational data (e.g., vector magnetogram) with the chosen order of accuracy.

Keeping in mind Eq. (12), we can rewrite Eq. (11) in terms of the primitive variables:

[TABLE]

where $\bar{I}^{P}$ is a matrix of size $(P\times(N+P))$ with $I^{P}_{ii}=1$ and $I^{P}_{ij}=0$ for $i\neq j$, and \mathbf{C}=\left|\begin{array}[]{c}\mathbf{C}^{P}(t)\\ 0\end{array}\right|. Let us define \bar{L}_{1}=\left|\begin{array}[]{c}\bar{I}^{P}\\ \bar{L}^{N}\end{array}\right|, and \bar{L}_{2}=\left|\begin{array}[]{c}0\\ \bar{L}^{N}\end{array}\right|.

After making necessary substitutions, the set of primitive variables at the inner boundary can be determined from

[TABLE]

It should be noted that our BCs are evolutionary (i.e., small changes in them lead to small changes in solution, see [30]), iff the left eigenvector matrix of the split system, $\bar{L}_{1}$, is non-singular. This determines the set of admissible physical BCs in terms of the primitive variables. It follows from our analysis that for realistic inflow condition only the plasma density and the radial component of magnetic field should always be specified as physical BCs. The choice of other physical BCs in terms of the primitive variables depends on the availability of appropriate measurements. Table 1 describes reasonable choices, as well as the number of physical and numerical BCs.

Thus, it can be seen that SDO/HMI data allow us to specify the necessary BCs, except for density. It was demonstrated in [31] that underestimating or overestimating density at the inner boundary does not affect the coronal solution. Hence, we can assume uniform density at the inner boundary since we do not have measurements.

We discretize our characteristic compatibility equations (see Eq. (14)) using TVD, Roe-type schemes. For the components of the vector magnetogram and horizontal velocity data which can be imposed as physical BCs (see Table 1), the data values are imposed at the boundary surface and they are used in combination with values in the first and second layers of inner cells to extrapolate the primitive variables into the first and second ghost cell layers, respectively. Since our characteristic BC formulation is in inertial reference frame, the synoptic maps for vector magnetogram and horizontal velocity data are also rotated taking the solar rotation into account. For the rest of the variables, their time variation is determined by solving characteristic compatibility equations (see Eq. (14)) in the ghost cells. These equations involve spatial derivatives of plasma variables which are solved by finite difference approximations taking into account the plasma parameter values in the first and second inner and ghost cell layers.

3 Results

3.1 Relaxation of a dipole magnetic field to a quasi-steady state

We present the results of a benchmark test that we performed to validate our characteristic BC formulations. In this test, we impose a dipole magnetic field with magnetic dipole moment equal to the solar dipole moment, $3.5\times 10^{29}$ N m $\textrm{T}^{-1}$, and Parker’s isothermal SW model as initial conditions for the magnetic field and hydrodynamic plasma parameters, respectively. Then, we relax this configuration to a quasi-steady state using our characteristic BCs. We also performed another simulation with exactly the same grid, numerical schemes and initial conditions, but this time, using a set of non-characteristic BCs (i.e., $n=1.5\times 10^{8}$ $\textrm{cm}^{-3}$, $B_{r}=(B_{r})_{\mathrm{Dipole}}$, $B_{\theta}=(B_{\theta})_{\mathrm{Dipole}}$, $v_{\theta}=v_{\phi}=0$, while $B_{\phi}$, $v_{r}$, and $p$ have zero derivative in the radial direction) instead. Notice that both sets of BCs have the same physical BCs corresponding to the subslow inflow case in Table 1, whereas the rest of the variables are either obtained from the compatibility equations or their radial derivatives are assumed to be zero. The comparison of the obtained solutions is shown in Figure 3.

The computational domain has a span of $1.03R_{\odot}\leq r\leq 20R_{\odot}$, $0\leq\theta\leq\pi$, and $0\leq\phi\leq 2\pi$. The computational grid is 128$\times$56$\times$40 in $r$, $\theta$ and $\phi$ directions, respectively, and is radially stretched. We use a TVD, Roe-type scheme and forward Euler method for spatial and temporal discretizations, respectively. The $\mathbf{\nabla}\cdot\mathbf{B}=0$ constraint is enforced following [32].

In Figure 3, dipole magnetic field lines that relaxed to a quasi-steady state are shown together with the plasma flow speed contours using characteristic and non-characteristic BCs. Moreover, radial velocity variations are also shown along the $y$-axis for both sets of BCs.

As seen from Figure 3, using the characteristic BCs leads to a correct magnetic field configuration whereas the non-characteristic approach shows artificial X-points in the vicinity of the Sun. Furthermore, non-characteristic BCs produce a substantial inflow on the solar equatorial plane which does not seem right.

3.2 Relaxation of the PFSS magnetic field to a quasi-steady state

In this section, we present results related to the relaxation of an initial PFSS magnetic field topology generated by using the HMI LOS magnetogram data observed on 2 January 2014 to a quasi-steady state. We obtained the PFSS solution by applying the spherical harmonics approach for which we took the spherical harmonics coefficients from the pfss_viewer software in IDL SolarSoft, and using 10 modes.

The computational domain has a span of $1.03R_{\odot}\leq r\leq 20R_{\odot}$, $0\leq\theta\leq\pi$, and $0\leq\phi\leq 2\pi$. The computational grid is 128$\times$60$\times$120 in $r$, $\theta$ and $\phi$ directions, respectively, and is radially stretched. The set of non-characteristic BCs applied at the inner boundary is given as follows: $n=1.5\times 10^{8}$ $\textrm{cm}^{-3}$, $T=1.3\times 10^{6}$ K, $B_{r}$ is determined from the HMI LOS magnetogram data, $v_{\theta}$ and $v_{\phi}$ are the horizontal velocity data shown in Figure 2, while $B_{\theta}$, $B_{\phi}$, and $v_{r}$ have zero derivative in the radial direction. The horizontal velocity data inferred by the time-distance helioseismology method have 3000$\times$1000 pixels resolution. We applied bilinear interpolation to interpolate the data to our computational grid at the inner boundary in order to be able to utilize in our BCs. It should be noted that differential rotation and meridional flow profiles are removed from the horizontal velocity data. Therefore, we impose empirical formulae for differential rotation [33] and meridional flow [34] at the entire inner boundary in addition to the data. The simulation is performed in the frame of reference corotating with the Sun.

In Figure 4, we show the magnetic field line distributions in the vicinity of the inner boundary with radial magnetic field line contours for the initial PFSS and quasi-steady solutions, and a full 3D view of the magnetic field line topology. In Figure 5, heliospheric current sheets corresponding to the initial PFSS and quasi-steady solutions are presented.

4 Conclusions

In this paper, we presented an overview of our data-driven MHD model of SW based entirely on observations and not relying on any non-reflection principle in the solar corona.

For this purpose, we developed a set of genuinely characteristic BCs which gives remarkable results already in a benchmark test.

Our model can be driven by a variety of observational data including HMI vector magnetogram and horizontal velocity data, and LOS magnetogram data from HMI, MDI, GONG, and WSO. In order to demonstrate this, we presented results related to the simulation of relaxation of a PFSS magnetic field topology to a quasi-steady state which is driven by HMI LOS magnetogram and horizontal velocity data using a set of non-characteristic BCs.

We will present results related to time-dependent simulations of our global corona model in our next publication. \ackThe authors acknowledge the support from the NASA grant NNX14AF41G and NSF SHINE grant AGS-1358386. We also acknowledge NSF PRAC award OCI-1144120 and related computer resources from the Blue Waters sustained-petascale computing project. Supercomputer allocations were also provided on SGI Pleiades by NASA High-End Computing Program award SMD-15-5860 and on Stampede by NSF XSEDE project MCA07S033.