ALFVÉN WAVE SOLAR MODEL (AWSoM): CORONAL HEATING

During the last few decades, considerable progress has been made in the understanding of the solar atmosphere due to the increased availability of observational data and the development of analytical and numerical models of the solar wind. One aspect of this development is the construction of complex three-dimensional (3D) models. These models can be validated with observations and further refined to improve the comparison. Important in the further progress of these models is a better understanding of the coronal heating and solar wind acceleration problem. Improvements in the theories of the coronal heating scenarios may result in more realistic models that are more reliable in predicting the solar wind conditions. Eventually, this may further improve the numerical forecasting of space weather events.

Recent Hinode observations suggested that there is more than enough energy in the chromospheric magnetic field fluctuations, which propagate away from the Sun, to heat the solar corona (SC) and maintain the temperatures at 1 MK (De Pontieu et al. 2007). With the Solar Dynamics Observatory (SDO) these waves were shown to be ubiquitously present in the transition region and low corona (McIntosh et al. 2011). These observations suggest that it is therefore appealing to develop a 3D SC and inner heliosphere (IH) model that is based on Alfvén wave turbulence and check if such a theoretical model reproduces emission features seen in the extreme ultraviolet (EUV) images. In our previous corona model (Sokolov et al. 2013), we demonstrated that Alfvén wave turbulence can indeed capture the overall observable EUV emission, with the exception of some of the details such as the emission around active regions.

There is a long history of Alfvén wave turbulence in the solar wind that dates back to the pioneering work of Coleman (1968), who concluded that turbulence is important in the solar wind near 1 AU based on Mariner 2 measurements. The earliest solar wind models that incorporated Alfvén wave turbulence were presented in Belcher & Davis (1971) and Alazraki & Couturier (1971), followed by two-dimensional (2D) global corona models by Usmanov et al. (2000) and Hu et al. (2003). Suzuki & Inutsuka (2006) constructed a self-consistent one-dimensional (1D) Alfvén wave turbulence model from the photosphere to IH that included wave reflection and mode conversion from Alfvén to slow waves. This model was generalized to 2D by Matsumoto & Suzuki (2012) to account for turbulent cascade. Alfvén waves that propagate outward from the Sun are partially reflected by gradients due to stratification, which produces waves propagating in opposite directions, see for instance Heinemann & Olbert 1980, Leroy 1980, Matthaeus et al. 1999, Dmitruk et al. 2002, Verdini & Velli 2007, Cranmer 2010, and Chandran et al. 2011. Counter-propagating waves are essential for the classical incompressible cascade (Velli et al. 1989) and hence the coronal heating. The reflection due to strong magnetic field gradients may be important in close proximity of the active regions, further intensifying the dissipation and the observable EUV emission. In the present paper, we will determine if this enhanced wave reflection does improve our model comparison with the observed EUV images.

Remote observations from the Ultraviolet Coronagraph Spectrometer have shown that the perpendicular ion temperature is much larger than the parallel ion temperature in the corona holes (Kohl et al. 1998; Li et al. 1998). Similarly, Helios observations have shown that the ion temperature in the IH is anisotropic as well (Marsch et al. 1982). Several models have been constructed that take the ion temperature anisotropy into account, see for instance Leer & Axford (1972), Isenberg (1984), Fichtner & Fahr (1991), and Chandran et al. (2011) for 1D analysis and Vásquez et al. (2003) and Li et al. (2004) for 2D analysis. The ion temperature anisotropy has, to the best of our knowledge, not yet been incorporated in a 3D solar wind model.

In recent years, our solar wind modeling efforts in the Space Weather Modeling Framework (SWMF; Tóth et al. 2012) were focused on creating an Alfvén wave turbulence-driven SC and IH model with a two-temperature approach for the electrons and ions. Our goal is to develop a single validated model that can produce realistic synthesized line-of-sight (LOS) images in the multi-wavelength EUV, and accurate 1 AU prediction of the solar wind properties. This model will also produce realistic background solar wind for simulations of coronal mass ejections (CMEs). van der Holst et al. (2010) constructed a solar wind model that has an inner boundary in the 1 MK corona and has separate electron and ion temperatures. It incorporated the collisional electron heat conduction and Alfvén wave transport and dissipation along the open field lines. This model was validated in Jin et al. (2012), while Evans et al. (2012) enhanced the model by including surface Alfvén waves. The two-temperature approach is important in producing a correct CME shock, otherwise a strong and unphysical heat precursor can appear ahead of the CME due to the heat conduction (Manchester et al. 2012; Jin et al. 2013). Meanwhile, Downs et al. (2010) developed a lower corona model that was able to reproduce synthesized EUV images, but for the coronal heating the model did still rely on ad hoc coronal heating functions. Sokolov et al. (2013) combined the aforementioned efforts into a single two-temperature model with Alfvén wave turbulence by incorporating the concept of balanced turbulence at the apex of the closed field lines, and via a concise analysis, the model was able to numerically resolve the upper chromosphere and transition region. This model was able to reproduce the overall morphology of coronal holes and active regions in the EUV images. Landi et al. (2013a) demonstrated that it was also able to capture the charge state evolution, while Oran et al. (2013) compared the model output to in situ Ulysses observations.

Lionello et al. (2009) demonstrated that a 3D magnetohydrodynamic (MHD) lower corona model based on a combination of ad hoc coronal heating functions and with heat conduction was able to well reproduce many features in the observed EUV images. Odstrcil et al. (2005) showed that the ENLIL IH MHD model prescribed by the empirical Wang–Sheeley–Arge (WSA) model (Arge & Pizzo 2000) was able to capture the 1 AU observations, while Cohen et al. (2007), Feng et al. (2010), and van der Holst et al. (2010) used the WSA to prescribe the coupled corona and inner IH simulations to obtain 1 AU results.

In the approach taken in the Alfvén wave solar model (AWSoM) described in the current paper, we no longer rely on ad hoc heating functions. The physics incorporated into the AWSoM model allows us to produce realistic LOS images and 1 AU solar wind predictions in one single model. We present the theoretical approach and results in two companion papers. The present paper is the first one and describes the coronal heating methodology and electron heat conduction and presents the resulting synthesized images of the multi-wavelength EUV emission. We also demonstrate the newly implemented low dissipation MP5 limiter (Suresh & Huynh 1997) in the BATS-R-US model. The reduced numerical dissipation allows us to clearly resolve the fine details in the solution, and allow a better comparison with observations. Our second paper, X. Meng et al. (in preparation), describes the impact of these changes on the solar wind acceleration and the 1 AU comparisons. The outline of the paper is as follows. In Section 2, we present details of the theoretical model and discuss the numerical implementation. We first describe in Section 2.1 the Alfvén wave turbulence approach for the two-temperature model. In Section 2.2 we provide the derivation of the Alfvén wave propagation, reflection, and dissipation. In Section 2.3, we generalize this model to three temperatures by using anisotropic ion temperatures and the isotropic electron temperature. The implementation is described in Section 2.4. In Section 3, we demonstrate the performance of this model for Carrington rotation (CR) 2107 by comparing EUV images with observations. We conclude in Section 4 and provide details of the derivations in the appendices.

We now present the equations of the AWSoM model. We first describe the two-temperature MHD equations for the electrons and ions to demonstrate the incorporation of the Alfvén wave turbulence and collisionless heat conduction. The details of the derivation of the Alfvén wave equations are presented in Section 2.2. In Section 2.3, these equations are generalized to ion temperature anisotropy. In Section 2.4, we give some details of the numerical implementation. We finalize this section with a discussion on the boundary conditions (Section 2.4.1).

2.1. Governing Equations of Two-temperature Model

The starting point of our model is the MHD equations, which well describe the large-scale, low-frequency phenomena of the SC and IH plasma. In the inertial frame, the mass conservation, momentum conservation, and induction equation are

$$\begin{equation} \frac{\partial \rho }{\partial t}+\nabla \cdot (\rho {\bf u})=0, \end{equation} \tag{ 1 }$$

$$\begin{eqnarray} \frac{\partial (\rho {\bf u})}{\partial t}&+&\nabla \cdot \left(\rho {\bf u}{\bf u}-\frac{{\bf B}{\bf B}}{\mu _0}\right)\nonumber\\ &+&\nabla \left(P_i+P_e+\frac{B^2}{2\mu _0}+P_A\right)= -\rho \frac{GM_\odot }{r^3}{\bf r}, \end{eqnarray} \tag{ 2 }$$

$$\begin{equation} \frac{\partial {\bf B}}{\partial t}- \nabla \times ({\bf u}\times {\bf B})=0, \end{equation} \tag{ 3 }$$

respectively. In addition, the initial conditions should satisfy ∇ · B = 0. The notation in these equations is as follows: ρ is the mass density, u is the velocity, assumed to be the same for the ions and electrons, B is the magnetic field, G is the gravitational constant, M_☉ is the solar mass, r is the position vector relative to the center of the Sun, and μ₀ is the permeability of vacuum. The Alfvén wave pressure, P_A, provides additional solar wind acceleration. The isotropic ion pressure P_i and electron pressure P_e are determined by the energy equations:

$$\begin{eqnarray} &\frac{\partial }{\partial t}\left(\frac{P_i}{\gamma -1}+\frac{\rho u^2}{2}+\frac{{\bf B}^2}{2\mu _0}\right)\nonumber\\ &\quad +\nabla \cdot \left[\left(\frac{\rho u^2}{2}+\frac{\gamma P_i}{\gamma -1}+\frac{B^2}{\mu _0}\right){\bf u}-\frac{{\bf B}({\bf u}\cdot {\bf B})}{\mu _0}\right]= \nonumber \\ & \quad -({\bf u}\cdot \nabla)\left(P_e+P_A\right)+ \frac{N_ik_B}{\tau _{ei}}\left(T_e-T_i\right)+Q_i-\rho \frac{GM_\odot }{r^3}{\bf r}\cdot {\bf u},\nonumber\\ \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &\frac{\partial }{\partial t}\left(\frac{P_e}{\gamma -1}\right) +\nabla \cdot \left(\frac{P_e}{\gamma -1}{\bf u}\right)+P_e\nabla \cdot {\bf u}=\nonumber\\ & \quad -\nabla \cdot {\bf q}_e +\frac{N_ik_B}{\tau _{ei}}\left(T_i-T_e\right)-Q_{\rm rad}+ Q_e. \end{eqnarray} \tag{ 5 }$$

in which T_{e, i} are the electron and ion temperatures, N_{e, i} are the electron and ion number densities, and k_B is the Boltzmann constant. We use the simple equation of state P_{e, i} = N_{e, i}k_BT_{e, i} and the polytropic index is γ = 5/3. We have the following additional energy contributions in the plasma energy equations. The optically thin radiative energy loss in the lower corona is given by

$$\begin{equation} Q_{\rm rad} = N_e N_h \Lambda (T_e), \end{equation} \tag{ 6 }$$

where N_h is the hydrogen number density and Λ(T_e) is the radiative cooling curve taken from the CHIANTI version 7.1 database (Landi et al. 2013b). The Coulomb collisional energy exchange between the ions and electrons depends on the relaxation time τ_ei. The electron heat flux q_e consists of two contributions. We use both the collisional formulation of Spitzer:

$$\begin{equation} {\bf q}_{e,S} = -\kappa _e T_e^{5/2}{\bf b}{\bf b}\cdot \nabla T_e, \end{equation} \tag{ 7 }$$

where b = B/B and κ_e ≈ 9.2 × 10⁻¹² W m⁻¹ K^−7/2, as well as the collisionless heat flux as suggested by Hollweg (1978):

$$\begin{equation} {\bf q}_{e,H} = \frac{3}{2} \alpha p_e {\bf u}, \end{equation} \tag{ 8 }$$

in which we assume α = 1.05. We smoothly transition between these formulations

$$\begin{equation} {\bf q}_e = f_S {\bf q}_{e,S} + (1 - f_S){\bf q}_{e,H}. \end{equation} \tag{ 9 }$$

Here, the fraction of Spitzer heat flux is defined as a function of r

$$\begin{equation} f_S = \frac{1}{1 + (r/r_H)^2}, \end{equation} \tag{ 10 }$$

where r_H = 5 R_☉ similar to Chandran et al. (2011). The Spitzer form is in this way the dominant heat flux contributor near the Sun where the density is high enough so that the electron temperature length scale L_T = T_e/|∇T_e| is much larger than the collisional mean free path of the electrons (this is not correct for part of the transition region, but we ignore that fact in the present paper), while far away from the Sun the collisionless heat flux is more significant. In Appendix C, we describe the numerical implementation of the collisionless heat flux. The electron and ion heating functions are denoted by Q_e and Q_i, respectively. Their sum equals the total turbulence dissipation per unit time and per unit volume. The partitioning of the dissipation into the coronal heating of the electrons and ions is obtained from results of linear wave theory and stochastic heating, see Chandran et al. (2011) for details and Appendix B for a brief summary. To determine the Alfvén wave pressure and total wave dissipation, we additionally solve for the propagation, reflection and dissipation of the wave energy densities, w_±, in which the + sign is for waves propagating in the direction parallel to B, while the − sign is for waves propagating antiparallel to B. The turbulence equations are derived in Section 2.2. Here, we summarize the final expressions for the time evolution of w_±:

$$\begin{equation} \frac{\partial w_\pm }{\partial t}+\nabla \cdot \left[({\bf u}\pm {\bf V}_A)w_\pm \right]+\frac{w_\pm }{2}(\nabla \cdot {\bf u})= \mp {\cal R}\sqrt{w_-w_+} -\Gamma _\pm w_\pm, \end{equation} \tag{ 11 }$$

where ${\bf V}_A = {\bf B}/\sqrt{\mu _0\rho }$ is the Alfvén speed. The third term on the left-hand side of Equation (11) represents the energy reduction in an expanding flow due to work done by the Alfvén wave pressure P_A = (w₊ + w₋)/2. The last term in Equation (11) is the wave dissipation per unit time and per unit volume. It is expressed in the form of the phenomenological dissipation rate

$$\begin{equation} \Gamma _\pm =\frac{2}{L_\perp }\sqrt{\frac{w_\mp }{\rho }}, \end{equation} \tag{ 12 }$$

which contains the transverse correlation length of the Alfvén waves in the plane perpendicular to B. Similar to Hollweg (1986), we use a simple scaling law L_⊥∝B^−1/2 with the proportionality constant $L_\perp \sqrt{B}$ as adjustable input parameter of the model. Since Γ_± depends on the returning wave w_∓, it is essential to include the partial reflection of the forward propagating wave w_±. The first term on the right hand side of Equation (11) is the new source term describing the conversion into oppositely propagating waves. The signed reflection rate ${\cal R}$ in this term is derived as

$$\begin{eqnarray} {\cal R} &=& \min \left[{\cal R}_{\rm imb},\max (\Gamma _\pm)\right]\nonumber\\ &&\times \left\lbrace \begin{array}{@{}lll}\left(1-2\sqrt{\frac{w_-}{w_+}}\right)&{\rm if}&4w_-\le w_+\\ 0&{\rm if}&1/4w_-<w_+<4w_-\\ \left(2\sqrt{\frac{w_+}{w_-}}-1\right)&{\rm if}&4w_+\le w_-\end{array}\right.,\qquad \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} && {\cal R}_{\rm imb} = \sqrt{\left[({\bf V}_A\cdot \nabla)\log V_A\right]^2 + \left({\bf b}\cdot [\nabla \times {\bf u}]\right)^2}. \end{eqnarray} \tag{ 14 }$$

The reflection rate consists of three parts: (1) for the strongly imbalanced turbulence, where min (w_±)/max (w_±) ≪ 1, and moderately imbalanced turbulence on open field lines and at the bottom of closed field lines, the wave reflection rate is represented by ${\cal R}_{\rm imb}$. In this case, the reflection is due to Alfvén speed gradients and vorticity along the field lines. (2) We additionally assume that the reflection rate is smaller than the maximum dissipation rate to limit the reflection rate in the transition region, which is accomplished by the min function. (3) The full expression of the reflection rate (13) includes a correction to the right of the curly bracket when the oppositely propagating waves are of equal wave energy density near the apex of the closed field lines. This correction presumes that the waves originating from the two foot points are uncorrelated and as a result, the reflection rate is negligible.

2.2. Propagation, Reflection, and Dissipation of Alfvén Waves

In this section, we derive the Alfvén wave energy density equations in several steps. In Section 2.2.1, we first derive the WKB wave transport equation in the absence of wave reflection and dissipation. The nonlinear cascade rate is discussed in Section 2.2.2. In Section 2.2.3, the wave reflection is presented and we arrive at the final expressions for the oppositely propagating waves. The limit of strongly imbalanced turbulence is further analyzed in Appendix A.

2.2.1. WKB Equation for Alfvén Turbulent Wave Energy Densities

In this section, we derive the equations governing the evolution of the wave energy densities, w_±. Our starting point is reduced MHD, which solves the non-conservative form of Equations (2) and (3):

$$\begin{eqnarray} &&\frac{\partial {\bf u}}{\partial t}+({\bf u}\cdot \nabla){\bf u}+ \frac{\nabla B^2}{2\mu _0\rho }+\frac{\nabla (P_e+P_i)}{\rho }= \frac{({\bf B}\cdot \nabla){\bf B}}{\mu _0\rho }-\frac{GM_\odot }{r^3}{\bf r},\nonumber\\ \end{eqnarray} \tag{ 15 }$$

$$\begin{equation} \frac{\partial {\bf B}}{\partial t}+({\bf u}\cdot \nabla){\bf B}+{\bf B}(\nabla \cdot {\bf u})=({\bf B}\cdot \nabla){\bf u}, \end{equation} \tag{ 16 }$$

as well as the continuity equation (1). We represent the magnetic field and velocity vectors as sums of regular and turbulent parts, ${\bf u}=\tilde{\bf u}+\delta {\bf u}$ and ${\bf B}=\tilde{\bf B}+\delta {\bf B}$ (below tildes are omitted) and simplify the equations for turbulent amplitudes by assuming the incompressibility conditions ∇ · δu = 0 and B · δB = 0:

$$\begin{eqnarray} \frac{\partial \delta {\bf u}}{\partial t}&+& ({\bf u}\cdot \nabla)\delta {\bf u}+(\delta {\bf u}\cdot \nabla)\delta {\bf u}+(\delta {\bf u}\cdot \nabla){\bf u} =\frac{({\bf B}\cdot \nabla)\delta {\bf B}}{\mu _0\rho }\nonumber \\ & +& \frac{(\delta {\bf B}\cdot \nabla)\delta {\bf B}}{\mu _0\rho }+\frac{(\delta {\bf B}\cdot \nabla){\bf B}}{\mu _0\rho }, \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} &\frac{\partial \delta {\bf B}}{\partial t} + ({\bf u}\cdot \nabla)\delta {\bf B}+ (\delta {\bf u}\cdot \nabla)\delta {\bf B}+(\delta {\bf u}\cdot \nabla){\bf B}+\delta {\bf B}(\nabla \cdot {\bf u})\nonumber \\ &\quad = ({\bf B}\cdot \nabla)\delta {\bf u}+ (\delta {\bf B}\cdot \nabla)\delta {\bf u}+(\delta {\bf B}\cdot \nabla){\bf u}, \end{eqnarray} \tag{ 18 }$$

$$\begin{equation} \frac{\partial \rho }{\partial t}+({\bf u}\cdot \nabla)\rho +(\delta {\bf u}\cdot \nabla)\rho +\rho \nabla \cdot {\bf u}=0. \end{equation} \tag{ 19 }$$

The equations for the Elsässer variables, ${\bf z}_\pm =\delta {\bf u}\mp \delta {\bf B}/\sqrt{\mu _0\rho }$ are obtained as the sum Equation (17) $\mp \frac{1}{\sqrt{\mu _0\rho }}\times$ Equation (18) $\pm \frac{\delta {\bf B}}{\rho \sqrt{\mu _0\rho }}\times$ Equation (19):

$$\begin{equation} \frac{d_\pm {\bf z}_\pm }{dt} +{\bf z}_\mp \cdot \nabla {\bf u}\mp \frac{{\bf z}_\mp \cdot \nabla {\bf B}}{\sqrt{\mu _0\rho }}- \frac{{\bf z}_\pm -{\bf z}_\mp }{4\rho } \frac{d_\mp \rho }{dt} =0, \end{equation} \tag{ 20 }$$

where (d_±/dt) = (∂/∂t) + (u ± V_A + z_∓) · ∇ and ${\bf V}_A=({\bf B}/{\sqrt{\mu _0\rho }})$ (see the analogous equation in Velli 1993; Chandran & Hollweg 2009; Chandran et al. 2009). Note, that for the regular plasma velocity and magnetic field the non-conservative form of Equations (2) and (3) is fulfilled, the wave pressure resulting from the turbulent magnetic field being:

$$\begin{equation} P_A^{\rm full}=\frac{(\delta {\bf B})^2}{2\mu _0}=\rho \left(\frac{{\bf z}_+^2}{8}+\frac{{\bf z}_-^2}{8}-\frac{{\bf z}_-\cdot {\bf z}_+}{4}\right). \end{equation} \tag{ 21 }$$

The equations for the wave energy densities, $w_\pm =\rho {\bf z}^2_\pm /4$, is obtained by the sum ρz_±/2 × Equation (20) $+3{\bf z}^2_\pm /8\times$ Equation (19):

$$\begin{eqnarray} &\frac{\partial w_\pm }{\partial t}+\nabla \cdot [({\bf u}\pm {\bf V}_A+{\bf z}_\mp)w_\pm ] +\frac{w_\pm }{2}(\nabla \cdot {\bf u})+\frac{\rho }{2}{\bf z}_\pm \nonumber\\ &\quad\cdot \left[({\bf z}_\mp \cdot \nabla){\bf u}\mp \frac{({\bf z}_\mp \cdot \nabla){\bf B}}{\sqrt{\mu _0\rho }}\right]+\frac{{\bf z}_+\cdot {\bf z}_-}{8}\frac{d_\mp \rho }{dt}=0. \end{eqnarray} \tag{ 22 }$$

As a first approximation we set the oppositely propagating wave z_∓ = 0 in the equations for w_± and obtain:

$$\begin{equation} \frac{\partial w_\pm }{\partial t}+\nabla \cdot \left[({\bf u}\pm {\bf V}_A)w_\pm \right]+\frac{w_\pm }{2}(\nabla \cdot {\bf u})=0. \end{equation} \tag{ 23 }$$

This WKB equation describes Alfvén wave propagation along the magnetic field lines (first two terms) and the wave energy reduction in the expanding plasma (the last term) because of the work done by the wave pressure, the latter in the WKB approximation can be represented as:

$$\begin{equation} P_A=\frac{1}{2}(w_++w_-), \end{equation} \tag{ 24 }$$

This approximation is valid if the waves propagate only in one direction and we neglect reflection in a gradually varying medium so that we have in this case z₊ ≠ 0 and z₋ = 0 or vice versa. Alternatively, one can consider the oppositely propagating waves originating from two footpoints of the closed magnetic field line by assuming that the sources of these waves are not correlated. In this case the quadratic in z_± wave energy densities, w_±, produced by each of the sources do not vanish, while the product of two random and non-correlated amplitudes may be assumed to have a zero average: <z₊ · z₋ > =0

2.2.2. Dissipation Rate

Within the more accurate approximation (but still within the WKB method), the nonlinear term ∇ · (z_∓w_±) results in the turbulent cascade and the wave energy dissipation (see, e.g., Dmitruk et al. 2002; Chandran & Hollweg 2009). The dissipation rate for the wave energy density, w₊, is controlled by the amplitude of the oppositely propagating wave, $|{\bf z}_-| =2\sqrt{w_-/\rho }$, and the correlation length, L_⊥, in the transverse (with respect to the magnetic field) direction, because for the Alfvén wave, propagating along the magnetic field, δu and δB are perpendicular to the magnetic field. Therefore, $\nabla \cdot ({\bf z}_\mp w_\pm)\sim ({2}/{L_\perp })\sqrt{({w_\mp }/{\rho })}w_\pm$ and the WKB equations with an account for nonlinear dissipation read:

$$\begin{equation} \frac{\partial w_\pm }{\partial t}+\nabla \cdot \left[({\bf u}\pm {\bf V}_A)w_\pm \right]+\frac{w_\pm }{2}(\nabla \cdot {\bf u}) =-\Gamma _\pm w_\pm, \end{equation} \tag{ 25 }$$

$$\begin{equation} \Gamma _\pm =\frac{2}{L_\perp }\sqrt{\frac{w_\mp }{\rho }} \end{equation} \tag{ 26 }$$

These equations work for the balanced turbulence (w₊ = w₋ = w). For imbalanced turbulence, they properly reduce the dissipation rate for the dominant wave by expressing this rate in terms of the amplitude of the minor oppositely propagating wave.

The system of Equations (1)–(5) once completed with Equations (25) and (26) is consistent as long as the wave pressure as in Equation (24) is used. The source of the dissipated turbulence energy in Equations (4) and (5) is balanced with the energy sink in Equation (25). The wave pressure is applied in Equation (2) as the momentum source, while the work done by the wave pressure is, accordingly, taken into account in Equation (25). The model consistency is explicitly pronounced in the fact that the sum of energy equations (4), (5), and (25) has a form of a conservation law:

$$\begin{eqnarray} &&\frac{\partial }{\partial t}\left(\frac{P_i+P_e}{\gamma -1}+\frac{\rho u^2}{2}+\frac{B^2}{2\mu _0}+w_++w_-\right)\nonumber\\ &&\quad +\nabla \cdot \left[\left(\frac{\rho u^2}{2}+\frac{\gamma (P_i+P_e)}{\gamma -1}+\frac{B^2}{\mu _0}\right){\bf u}-\frac{{\bf B}({\bf u}\cdot {\bf B})}{\mu _0}\right] \nonumber \\ &&\quad +\nabla \cdot \left[\left(w_++w_-+P_A\right){\bf u}+\left(w_+-w_-\right){\bf V}_A\right]\nonumber\\ &&= -\nabla \cdot {\bf q}_e-Q_{\rm rad} -\rho \frac{GM_\odot }{r^3}{\bf r}\cdot {\bf u}. \end{eqnarray} \tag{ 27 }$$

The non-conservative sources in the right hand side, which do not have a divergence-like form, are due to the energy losses for emission and for the work done against the gravitational force. The conservative property of the governing equations is important, since it ensures the consistency of the model. Below, we carefully keep the conservative property of the model while including the wave reflection.

2.2.3. The Model of Reflection

Equation (20) for z_± demonstrates that even in a linear approximation, the WKB approach dismisses the correlations between inward and outward propagating waves. Indeed, in the Equation (20) for z₊ there are source terms linearly proportional to z₋ and vise versa, while in the WKB approximation these source terms describing the conversion between the oppositely going waves are omitted. The omitted correlations are important for turbulent cascade (Tu & Marsch 1995; Dmitruk et al. 2002). We include them into the solar model as follows. First, in the part of the Elsässer Equation (20) that is responsible for the wave reflection, we use the continuity equation (19) to obtain

$$\begin{equation} \frac{1}{\rho }\frac{d_\mp \rho }{dt} = -\nabla \cdot {\bf u}\mp {\bf V}_A\cdot \nabla \log \rho \mp \frac{1}{\sqrt{\mu _0\rho }}\delta {\bf B}\cdot \nabla \log \rho. \end{equation} \tag{ 28 }$$

The third term on the right hand side is a nonlinear term, which can lead to wave reflection due to very steep density gradients in the direction of the wave amplitude δB. This direction is under the incompressibility condition, B · δB = 0, transverse to the magnetic field. Such surface reflection can for instance arise at coronal loop boundaries and streamer boundaries, ultimately resulting in enhanced heating near these locations. In the present paper, we do not consider these effects, so that (1/ρ)(d_∓ρ/dt) ≈ −∇ · u∓V_A · ∇log ρ. Now, we apply this transformation in Equation (22) and repeat the above procedure to check how Equation (27) is modified. Second, if we admit a non-zero correlator <z₊ · z₋ >, then we should use the full magnetic pressure (21) in the momentum equation. This point is also clear from the observation, that the work done by the wave pressure in Equation (22) now becomes [(w_±/2) − ρ(z₋ · z₊/8)]∇ · u, with their total being equal to $P_A^{\rm full}\nabla \cdot {\bf u}$. One can find that in the equation analogous to Equation (27) the wave pressure is changed for $P_A^{\rm full}$ and the non-divergent term breaking the energy conservation appears in the left hand side as follows:

$$\begin{eqnarray} && \frac{\partial }{\partial t}\left(\frac{P_i+P_e}{\gamma -1}+\frac{\rho u^2}{2}+\frac{B^2}{2\mu _0}+w_++w_-\right)\nonumber\\ &&\quad +\nabla \cdot \left[\left(\frac{\rho u^2}{2}+\frac{\gamma (P_i+P_e)}{\gamma -1}+\frac{B^2}{\mu _0}\right){\bf u}-\frac{{\bf B}({\bf u}\cdot {\bf B})}{\mu _0}\right] \nonumber\\ &&\quad +\nabla \cdot \left[\left(w_++w_-+P_A^{\rm full}\right){\bf u}+\left(w_+-w_-\right){\bf V}_A\right]+Q_{\rm noncons}\nonumber\\ &&= -\nabla \cdot {\bf q}_e-Q_{\rm rad} -\rho \frac{GM_\odot }{r^3}{\bf r}\cdot {\bf u}, \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} &&Q_{\rm noncons}= \frac{\rho }{2}\left[{\bf z}_+\cdot ({\bf z}_-\cdot \nabla){\bf u}+{\bf z}_-\cdot ({\bf z}_+\cdot \nabla){\bf u}\right.\nonumber\\ &&\quad\left. +\,\frac{{\bf z}_-\cdot ({\bf z}_+\cdot \nabla){\bf B}-{\bf z}_+\cdot ({\bf z}_-\cdot \nabla){\bf B}}{\sqrt{\mu _0\rho }}\right]. \end{eqnarray} \tag{ 30 }$$

Since the tensor z₋z₊ + z₊z₋ is symmetric, one can find that the non-conservative energy source involves a contribution proportional to the symmetric part of the deformation velocity tensor, (∂u_i/∂x^j) + (∂u_j/∂xⁱ) and another term proportional to ∇ × B. The energy source can also be expressed as follows:

$$\begin{eqnarray} Q_{\rm noncons}&&= \rho \delta {\bf u}\cdot (\delta {\bf u}\cdot \nabla){\bf u}-\frac{1}{\mu _0}\delta {\bf B}\cdot (\delta {\bf B}\cdot \nabla){\bf u}\nonumber\\ &&\quad +\frac{1}{\mu _0}\delta {\bf B}\cdot (\delta {\bf u}\cdot \nabla){\bf B}-\frac{1}{\mu _0}\delta {\bf u}\cdot (\delta {\bf B}\cdot \nabla){\bf B}. \end{eqnarray} \tag{ 31 }$$

The reason for this energy non-conservation is that the system of Alfvén waves, which are assumed to be transversely polarized (δz_± · B = 0) and interact with each other (via reflection) and with the moving plasma (via the turbulent heating and wave pressure), is not closed. Indeed, one can take the following linear combination: ρu · Equation (17) +(B/μ₀) · Equation (18) + ρδu · u × Equation (19) + δu · Equation (15) +(δB/μ₀) · Equation (16) + (u²/2)  ×  [Equation (19)−Equation (1)] (the latter subtraction is applied because (u²/2) × Equation (1) contributes to Equation (4) and should not be doublecounted), which gives:

$$\begin{eqnarray} &&\frac{\partial }{\partial t}\left(\rho {\bf u}\cdot \delta {\bf u}+\frac{\delta {\bf B}\cdot {\bf B}}{\mu _0}\right)+\nabla \cdot \left[\delta {\bf u}\left(\frac{\rho u^2}{2}+\frac{B^2}{\mu _0}\right)\right.\nonumber\\ &&\quad \left.+\left({\bf u}+\delta {\bf u}\right)\left(\rho \delta {\bf u}\cdot {\bf u}+2\frac{\delta {\bf B}\cdot {\bf B}}{\mu _0}\right)\right] \nonumber \\ &&\quad -\frac{1}{\mu _0}\nabla \cdot \left[\delta {\bf B}\left({\bf u}\cdot {\bf B}\right)+\left({\bf B}+\delta {\bf B}\right)\left(\delta {\bf u}\cdot {\bf B}+\delta {\bf B}\cdot {\bf u}\right) \right]\nonumber\\ &&\quad -Q_{\rm noncons}=0. \end{eqnarray} \tag{ 32 }$$

Here, we did not assume that ∇ · δu = 0 and δB · B = 0, as we did before, and we account for the term, ∇(δB · B), in Equation (17) and the term (B + δB)∇ · δu in Equation (18). We also accounted for the term δB · B∇ · u, which we omitted while deriving Equation (22) under the assumption of ∇ · δu = 0. We omitted all terms related to the ion and electron pressures and their variations. We obtain yet another non-closed conservation law, which governs the change in energy related to the variation in the longitudinal magnetic field. The non-convervative source, Q_noncons appears to come with the opposite sign to the newly derived equation, meaning that this source describes the energy nonlinear conversion from transverse Alfvén waves to the compressible mode with δB · B ≠ 0, i.e., magnetosonic waves. We arrive at an important conclusion: the model of Equation (22) allowing a correlation between the oscillations in the oppositely propagating waves is not closed unless the model for compressible turbulence is included. As long as, here, we do not include such model, we are to eliminate the turbulent energy sink to nowhere and add the term, −Q_noncons/2, to the left hand side of each of Equation (22):

$$\begin{eqnarray} && \frac{\partial w_\pm }{\partial t}+\nabla \cdot \left[({\bf u}\pm {\bf V}_A)w_\pm \right]+\left(\frac{w_\pm }{2}-\frac{\rho {\bf z}_+\cdot {\bf z}_-}{8}\right)(\nabla \cdot {\bf u}) \nonumber \\ &&\quad =\frac{\rho }{4}\left[ [{\bf z}_\pm \times {\bf z}_\mp ]\cdot [\nabla \times {\bf u}]\pm \frac{{\rm Tr}[({\bf z}_\pm {\bf z}_\mp +{\bf z}_\mp {\bf z}_\pm)\cdot \nabla {\bf B}]}{\sqrt{\mu _0\rho }}\right. \nonumber\\ &&\qquad \left. \pm ({\bf z}_+\cdot {\bf z}_-)({\bf V}_A\cdot \nabla)\log \sqrt{\rho }\right] -\Gamma _\pm w_\pm, \end{eqnarray} \tag{ 33 }$$

where we used an easy-to-derive identity, Tr[(ab − ba) · ∇c] = [b × a] · [∇ × c], for any three vectors, a, b, c. We can further simplify this expression by assuming the transverse polarization of the Alfvén waves: z_± · B = 0. For δB part of the Elsässer variables, ${\bf z}_\pm =\delta {\bf u}\mp \delta {\bf B}/\sqrt{\mu _0\rho }$, the condition, δB · B = 0, had been already introduced above, while the requirement for the velocity oscillation δu to be approximately perpendicular to the magnetic field line is a consequence of the assumed incompressibility property of the turbulence, ∇ · δu = 0, together with the expectation that the turbulent wave vectors are directed along the magnetic field, so that to satisfy the incompressibility condition the velocity oscillations should, rather, be transverse. For transverse waves in the 3-by-3 tensor (∇B) as present in the expression for the wave conversion rate, one can leave only 2-by-2 transverse components in the plane perpendicular to the magnetic field, (∇B)_⊥⊥. Only the symmetric part of this tensor matters, which possesses two eigenvalues (for curl-free, hence, potential magnetic field, these eigenvalues can be expressed in terms of the curvature radii of the equipotential surface). If we neglect the difference between two eigenvalues, we can admit that the symmetric part of the tensor is proportional to 2-by-2 unity tensor, the proportionality coefficient can be found by observing that the trace of the 2-by-2 unity tensor equals 2, while Tr[(∇B)_⊥⊥] = ∇_⊥ · B_⊥ = ∇ · B − B · ∇log |B|. Under the specified assumptions, the equation accounting for the Alfvén turbulent wave reflection, that is, the energy exchange between two oppositely propagating waves, which reduces the amplitude of the outgoing wave and amplifies the incoming wave, reads:

$$\begin{eqnarray} & \frac{\partial w_\pm }{\partial t}+\nabla \cdot \left[({\bf u}\pm {\bf V}_A)w_\pm \right] +\frac{w_\pm }{2}(\nabla \cdot {\bf u}) \nonumber \\ &\quad =\mp \frac{\rho }{4}\left\lbrace [{\bf z}_-\times {\bf z}_+]\cdot [\nabla \times {\bf u}]\right.\nonumber\\ &\qquad \left. +\, ({\bf z}_+\cdot {\bf z}_-)({\bf V}_A\cdot \nabla)\log V_A\right\rbrace -\Gamma _\pm w_\pm. \end{eqnarray} \tag{ 34 }$$

Note, that we omitted the contribution from the correlator z₋ · z₊ to the work done by the wave pressure (see the modified multiplier at ∇ · u). Accordingly, we do not use this contribution to the wave pressure in the momentum equation and return to Equation (24). The reason for this is that the correlator z₋ · z₊ is assumed to be non-zero only when the waves of one direction dominate over the opposite ones. In this case, the contribution from this correlator to the wave pressure is negligible compared with that from the dominant wave: $|{\bf z}_-\cdot {\bf z}_+|\le \Vert {\bf z}_-\Vert \Vert {\bf z}_+\Vert \ll \max ({\bf z}_\pm ^2)$. On the other hand, this contribution may be a large loss term in the wave energy equation for the minor wave. The neglect can be justified only if the reflection coefficient is limited in such way, that the maximum admissible value, corresponding to the equality case in the following estimate:

$$\begin{eqnarray} && \frac{\left|[{\bf z}_-\times {\bf z}_+]\cdot [\nabla \times {\bf u}]+ ({\bf z}_+\cdot {\bf z}_-)({\bf V}_A\cdot \nabla)\log V_A\right|}{|{\bf z}_-||{\bf z}_+|}\nonumber\\ &&\quad \le \sqrt{\left({\bf b}\cdot [\nabla \times {\bf u}]\right)^2+[({\bf V}_A\cdot \nabla)\log V_A]^2}, \end{eqnarray} \tag{ 35 }$$

can be only achieved if one of the waves dominates over the opposite, so that min (w_±)/max (w_±) ≪ 1. In this case, the amplitude and polarization of the minor wave are imposed by those for the dominant wave, the sign of the correlator is governed by the requirement that, in the course of reflection, the dominant wave should decrease, the oppositely propagating "reflected" wave should grow. If the oppositely propagating waves are comparable, they are assumed to be non-correlated. The following choice of the final expression for the reflection coefficient satisfies all the listed requirements:

$$\begin{equation} \frac{\partial w_\pm }{\partial t}+\nabla \cdot \left[({\bf u}\pm {\bf V}_A)w_\pm \right]+\frac{w_\pm }{2}(\nabla \cdot {\bf u})= \mp {\cal R}\sqrt{w_-w_+} -\Gamma _\pm w_\pm, \end{equation} \tag{ 36 }$$

where

$$\begin{eqnarray} && {\cal R} = {\cal R}_{\rm imb} \left\lbrace \begin{array}{@{}lcl}\left(1-2\sqrt{\frac{w_-}{w_+}}\right)&{\rm if}&4w_-\le w_+\\ 0&{\rm if}&1/4w_-<w_+<4w_-\\ \left(2\sqrt{\frac{w_+}{w_-}}-1\right)&{\rm if}&4w_+\le w_-\end{array}\right.,\qquad \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} && {\cal R}_{\rm imb} = \sqrt{\left({\bf b}\cdot [\nabla \times {\bf u}]\right)^2+[({\bf V}_A\cdot \nabla)\log V_A]^2}. \end{eqnarray} \tag{ 38 }$$

The above considerations may be limited by the assumption of a small reflection coefficient, which, probably, should be less than the dissipation rate, above the other criteria. Therefore, Equation (37) may be bounded as follows:

$$\begin{eqnarray} {\cal R} &&= \min \left[{\cal R}_{\rm imb},\max (\Gamma _\pm)\right] \left[\max \left(1-2\sqrt{\frac{w_-}{w_+}},0\right)\right.\nonumber\\ &&\quad \left.-\max \left(1-2\sqrt{ \frac{w_+}{w_-}},0\right)\right]. \end{eqnarray} \tag{ 39 }$$

By comparing Equations (34) and (36) we can find the correlators of amplitudes of the counter-propagating waves:

$$\begin{eqnarray} &&\frac{\rho }{4}\left({\bf z}_-\cdot {\bf z}_+\right)=\frac{({\bf V}_A\cdot \nabla)\log V_A}{\left({\bf b}\cdot [\nabla \times {\bf u}]\right)^2+[({\bf V}_A\cdot \nabla)\log V_A]^2}{\cal R}\sqrt{w_-w_+},\nonumber\\ \end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray} &&\frac{\rho }{4}{\bf b}\cdot \left[{\bf z}_-\times {\bf z_+}\right]\nonumber\\ &&\quad =\frac{\left({\bf b}\cdot [\nabla \times {\bf u}]\right)}{\left({\bf b}\cdot [\nabla \times {\bf u}]\right)^2+[({\bf V}_A\cdot \nabla)\log V_A]^2}{\cal R}\sqrt{w_-w_+}.\nonumber\\ \end{eqnarray} \tag{ 41 }$$

In this way, the contribution from the first correlator to the wave pressure, $P^{\rm full}_A$ could be accounted for, but we do not do this in the present paper and use Equation (24). Equations (36) and (37) describing the Alfvén wave propagation, reflection and dissipation close the system of the MHD equation in a physically consistent way, as long as the energy conservation law, Equation (27), is fulfilled.

The reflection model used in Matthaeus et al. (1999) is very similar to ours due to the following common features: (1) the reflection turns to zero in balanced turbulence, (2) the sign of effect is such that the reflection reduces the amplitude of the dominant wave and enhances the counter-propagating minor wave, and (3) the magnitude of the reflection coefficient is controlled by the gradient of the Alfvén speed. Our account of the vorticity is also not quite new, since the effect of sheared flow on the mode conversion in the solar atmosphere is discussed by Hollweg & Kaghashvili (2012) and Hollweg et al. (2013). The important distinction is that, for the sake of the model consistency and energy conservation, we ruled out the non-conservative sources from Equations (36). Now, we can evaluate the neglected terms, Q_noncons/2 and −(ρ/8)(z₋ · z₊)(∇ · u), which might be added to the right hand side of Equations (36) with the opposite sign: (ρ/8)(z₋ · z₊)(∇ · u) − (Q_noncons/2). In the low steady-state SC, the plasma moves mostly along the magnetic field lines: u = uB/B, which allows us to express the transverse components of the velocity derivatives in terms of those of the magnetic field: (∇u)_⊥⊥ = (u/B)(∇B)_⊥⊥. If we neglect the contribution from ∇ × B into the non-conservative source and express approximately ∇ · u ≈ −(u · ∇)log ρ from the continuity equation, the energy source due to the Alfvén wave interaction with the compression mode reads:

$$\begin{eqnarray} &&\frac{\rho }{8} ({\bf z}_-\cdot {\bf z}_+)(\nabla \cdot {\bf u})- \frac{Q_{\rm noncons}}{2}\approx \frac{\rho }{4}({\bf z}_-\cdot {\bf z}_+)({\bf u}\cdot \nabla)\nonumber\\ &&\quad\times ({-}\log \sqrt{\rho }+\log B)=\frac{\rho }{4}({\bf z}_-\cdot {\bf z}_+)({\bf u}\cdot \nabla)\log V_A. \end{eqnarray} \tag{ 42 }$$

Thus, we do not introduce the mode conversion term proportional to (u · ∇)log V_A into our model. The reason for this omission is that including this compressible MHD turbulence term would break the energy conservation.

2.3. Generalization to Ion Temperature Anisotropy

Due to observational evidence of ion temperature anisotropy in the lower corona (Kohl et al. 1998; Li et al. 1998) and in the IH (Marsch et al. 1982), we have generalized our solar wind model to anisotropic ion temperatures. The implementation and global magnetosphere application of the anisotropic ion pressure is presented in Meng et al. (2012a, 2012b). Here we will use the same implementation in the solar context.

The equation for the ion pressure (4) is now decomposed into two equations for both the ion pressure component perpendicular to the magnetic field, P_i⊥, and the ion pressure component parallel to the magnetic field P_i∥. However, for convenience, we solve for the averaged ion pressure P_i = (2P_i⊥ + P_i∥)/3 instead of P_i⊥. The ion pressures are determined by the equations

$$\begin{eqnarray} &&\frac{\partial }{\partial t}\left(\frac{P_i}{\gamma -1} + \frac{\rho u^2}{2} + \frac{B^2}{2\mu _0}\right) + \nabla \cdot \left[ \left(\frac{\rho u^2}{2} + \frac{P_i}{\gamma - 1} + \frac{B^2}{\mu _0}\right){\bf u}\right.\nonumber\\ &&\quad \left.+\, {\bf P}_i\cdot {\bf u} - \frac{{\bf B}({\bf u}\cdot {\bf B})}{\mu _0} \right] = - {\bf u}\cdot \nabla (P_e+P_A)\nonumber \\ &&\quad + \frac{N_ik_B}{\tau _{ei}}(T_e-T_i) + Q_i - \rho \frac{GM_\odot }{r^3} {\bf r}\cdot {\bf u}, \end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray} &&\frac{\partial P_{i\parallel }}{\partial t} + \nabla \cdot (P_{i\parallel }{\bf u}) + 2P_{i\parallel }{\bf b}\cdot (\nabla {\bf u})\cdot {\bf b} = \frac{\delta P_{i\parallel }}{\delta t}\nonumber\\ &&\quad + (\gamma - 1)\frac{N_ik_B}{\tau _{ei}}(T_e -T_{i\parallel }) + (\gamma - 1)Q_{i\parallel }, \end{eqnarray} \tag{ 44 }$$

where T_i∥ is the parallel ion temperature obtained from the equation of state P_i∥ = N_ik_BT_i∥ and P_i = P_i⊥I + (P_i∥ − P_i⊥)bb is the ion pressure tensor. The second term on the right hand sides of Equations (43) and (44) are the collisional energy exchanges with the electrons. The third term on the right hand sides are the heating functions Q_i and Q_i∥ for the averaged ion and parallel ion pressure, respectively. The sum of the electron and averaged ion heating functions, Q_e + Q_i, is equal to the total turbulence dissipation per unit volume per unit time, Γ₊w₊ + Γ₋w₋. The partitioning of the wave dissipation into Q_e, Q_i, and Q_i∥ is described in Appendix B. The first term on the right hand side of Equation (44), δP_i∥/δt, is for the relaxation of the pressure anisotropy by the parallel firehose, mirror, and ion-cyclotron instability constraints. If those instability criteria are met, we reduce the pressure anisotropy so that the plasma is stable again. Details about the stability formulation, implementation and results are given in Meng et al. (2012b). The anisotropic ion pressure also modifies the momentum equation (2):

$$\begin{eqnarray} \frac{\partial \rho {\bf u}}{\partial t} &+& \nabla \cdot \left[ \rho {\bf u}{\bf u} + (P_{i\parallel } - P_{i\perp }){\bf b}{\bf b} - \frac{1}{\mu _0}{\bf B}{\bf B}\right]\nonumber\\ \quad &+& \nabla \left(P_{i\perp } + P_e + \frac{B^2}{2\mu _0} + P_A \right) = -\rho \frac{GM_\odot }{r^3}{\bf r}, \end{eqnarray} \tag{ 45 }$$

in which the second term on the left hand side contains a new contribution due to pressure anisotropy. We further assume that the anisotropic pressure does not significantly change the Alfvén wave turbulence, and hence we use the turbulence as formulated for isotropic temperatures in Section 2.1.

2.4. Model Implementation

In this section we present some details of the implementation of the improved solar wind model. This model uses the numerical schemes of the BATS-R-US MHD solver and the overarching SWMF, see Tóth et al. (2012) for a description of the SWMF and BATS-R-US tools. The SWMF is a software framework for modeling various space physics domains in a single coupled model. It has been used, besides space weather applications for the coupled Sun–Earth system, for many planetary, moon and comet applications as well as the outer heliosphere. It has recently been extended to applications of radiation hydrodynamics in the context of laser-driven high-energy-density physics (van der Holst et al. 2011; van der holst et al. 2013). The new components of the SWMF, presented in this paper, are the SC and IH.

The SC model uses a 3D spherical grid with the radial coordinate ranging from 1 R_☉ to 24 R_☉. The grid is highly stretched toward the Sun with smallest radial cell size Δr = 10⁻³ R_☉ to numerically resolve the steep density gradients in the upper chromosphere. We artificially broaden the transition region similar to that as described in Sokolov et al. (2013) and Lionello et al. (2009) to be able to resolve this region. The grid is block decomposed using the block-adaptive tree library (Tóth et al. 2012). This library is a tool to create, load balance and message pass the adaptive refined mesh and solution data. In the simulations of this paper, the grid blocks consist of 6 × 4 × 4 mesh cells. Inside r = 1.7 R_☉, the angular resolution is 256 cells in longitude and 128 cells in latitude corresponding to an angular cell size of 14, while outside that radius the grid is one level less refined. The system of equations described in Sections 2.1 and 2.3 are solved in the heliographic rotating frame by including centrifugal and Coriolis forces −ρ[Ω × (Ω × r) + 2Ω × u] in the momentum equation and adding the centrifugal contribution −ρu · [Ω × (Ω × r)] to the ion energy equations (4) and (43). Here Ω is the angular velocity of the Sun. We assume a uniform solar rotation with a 25.38 days period so that Ω = 2.865 × 10⁻⁶ rad s⁻¹. For steady-state simulations, we use local time stepping, which speeds up the convergence relative to time accurate simulations. During the steady-state convergence, we apply one additional level of mesh refinement at the heliospheric current sheet (HCS). To resolve the details in the LOS EUV images, we also demonstrate higher resolution in latitude and longitude by using 6 × 6 × 6 grid blocks and hence, an angular cell size of 0.94° near the Sun, in combination with the numerical scheme based on the spatially fifth-order MP5 limiter (Suresh & Huynh 1997) instead of our standard second-order shock-capturing schemes (Tóth et al. 2012).

Details of the IH setup and simulations are provided in X. Meng et al. (in preparation).

2.4.1. Boundary Conditions

Here, we limit our discussion to the pre-specified boundary conditions only, and refer the reader to Sokolov et al. (2013) for a more complete description. The radial magnetic field component B_r is prescribed using synoptic magnetogram data in the following way: this data is first extrapolated to a 3D potential field source surface (PFSS) solution using either spherical harmonics or the finite difference iterative potential-field solver (FDIPS). In the current paper, we use FDIPS since this method avoids the ringing patterns near regions of concentrated magnetic fields to which the spherical harmonics method is susceptible, see Tóth et al. (2011). The PFSS magnetic field is used both as the initial condition and to set the boundary conditions.

The boundary condition for the Alfvén wave energy density is empirically set by prescribing the Poynting flux of the outgoing waves (w is w₊ for positive B_r and w₋ for negative B_r): S_A = V_Aw∝B_☉, where B_☉ is the field strength at the inner boundary and the proportionality constant is estimated in Sokolov et al. (2013) as (S_A/B)_☉ = 1.1 × 10⁶ W m⁻² T⁻¹. Under the assumption of sufficiently small returning flux, this estimate of the Poynting-flux-to-field ratio is equivalent to the following averaged velocity perturbation

$$\begin{equation} \langle \delta {\bf u}_\perp \cdot \delta {\bf u}_\perp \rangle ^{1/2} \approx 15\;{\rm km}\,{\rm s}^{-1}\left(\frac{3\times 10^{-11}\;{\rm kg}\, {\rm m}^{-3}}{\rho }\right)^{1/4}, \end{equation} \tag{ 46 }$$

where the mass density 3 × 10⁻¹¹ kg m⁻³ (ion number density N_i = 2 × 10¹⁶ m⁻³) corresponds to the upper chromosphere. This value is compatible with the Hinode observations of the turbulent velocities of 15 km s⁻¹ (De Pontieu et al. 2007). Hence, the energy density of the outgoing wave is set to $w = (S_A/B)_\odot \sqrt{\mu _0\rho }$. The returning wave energy density is absorbed by setting it to zero.

The temperatures are all set to the same value T_e = T_i = T_i∥ = T_☉ = 50, 000 K uniformly at the inner boundary. In Sokolov et al. (2013) it was demonstrated that the grid spacing of Δr = 10⁻³ R_☉ is, in this case, sufficient to numerically resolve the density scale height. We overestimate the density for this temperature by an order of magnitude with the value of N_e = N_i = N_☉ = 2 × 10¹⁷ m⁻³ at the inner boundary. This overestimate prevents chromospheric evaporation and extends the upper chromosphere to reach the correct lower density, but does not significantly change the global solution as shown in Lionello et al. (2009).

For this paper, we selected CR2107 (2011 February 16–March 16). This rotation was also used in the validation studies of Sokolov et al. (2013). In that paper, we were able to reproduce the overall morphology of the coronal holes and active regions in the LOS EUV images of SDO and STEREO. However, details in those images did not show up. It is the goal of the present section to demonstrate that with the new version of the AWSoM model, we are now able to produce high quality synthesized images that capture details of the EUV observations. The validation of this model with in situ data at 1 AU will be presented in X. Meng et al. (in preparation).

To simulate a background solar coronal solution, we need to specify the radial magnetic field component at the inner boundary, which is located in the upper chromosphere. We obtain this magnetic field in the following way: the synoptic map CR2107 of the SDO/Helioseismic and Magnetic Imager (HMI) is used. The polar field of this map is corrected with a two-dimensional and third-order polynomial fitting of the data above 60° (Sun et al. 2011). We use FDIPS to generate an initial condition for the magnetic field and the boundary condition values for the radial magnetic field component (see Figure 1).

Zoom In Zoom Out Reset image size

In all our results we will use parameter values as summarized in Table 1, unless stated otherwise. The boundary condition values are the same as in Sokolov et al. (2013). The value of $L_\perp \sqrt{B}$ is twice larger compared to the value used in Sokolov et al. (2013) due to a factor two difference in the definition of this parameter. The stochastic heating parameter h_s and the collisionless heat conduction parameter r_H are assigned with the same values as in Chandran et al. (2011).

Table 1. The Model Parameters

Parameter	Value
N☉	2 × 1017 m−3
T☉	50, 000 K
(SA/B)☉	1.1 × 106 W m−2 T−1
$L_\perp \sqrt{B}$	$1.5\times 10^5\;{\rm m}\,\sqrt{\rm T}$
hS	0.17
α	1.05
rH	5 R☉

Download table as:  ASCIITypeset image

3.1. Heat Partitioning

We will first demonstrate the heat partitioning of the turbulence dissipation for our three-temperature model. In this simulation, we used the version of this model with the wave reflection term in Equations (11)–(14). The steady-state solution is obtained with the spatially second-order shock-capturing scheme. In Figure 2, we show the three obtained temperatures T_i⊥, T_i∥, and T_e from top to bottom in the panels on the left. These temperatures are shown in the meridional slice X = 0 along with a few projected field lines to indicate the location of open and closed field lines. Very close to the Sun, the three temperatures are nearly the same. This is due to the high density near the Sun, resulting in a sufficiently high rate of Coulomb collisions that equilibrate the temperatures. The collision rate decreases with the density, so that further away from the Sun the collisions are too infrequent to equilibrate the temperatures. The significant heating of the perpendicular ion temperature in the polar coronal holes is due to the stochastic heating. The heat partitioning fractions of the coronal heating into ion perpendicular, ion parallel, and electron heating are shown in the panels on the right from top to bottom, respectively. The ion perpendicular heating, due to the stochastic heating process, dominates in the lower corona sufficiently far away from both the Sun and the HCS. The parallel ion heating is only significant very close to the HCS due to the high plasma beta β_i, while the electron heating is important very close to the Sun and around the (HCS).

Zoom In Zoom Out Reset image size

3.2. EUV Comparison

From the density and electron temperature distribution, calculated with our solar model, we can produce synthesized EUV that can be compared with the observed ones. Such a comparison serves as a check for the performance of the coronal heating model. The presented model accounts for the partial reflection of the outward propagating waves, which is accompanied by the generation of counter-propagating waves. These oppositely propagating waves are ultimately responsible for the turbulent cascade rate and hence, the coronal heating. The distinct feature in the present model is the enhanced reflection in the presence of strong magnetic fields, such as in close proximity of active regions, that can increase the dissipation and thereby intensify the observable EUV emission.

To better resolve the details in the synthesized LOS images, the latitudinal and longitudinal resolution is increased by using adaptive mesh refinement grid size of 6 × 6 × 6. The first attempt to use the full model as described in Section 2 did not yet provide us the desired LOS image quality. The problem is due to the numerical inaccuracy in the Alfvén speed gradients in the reflection source term. To overcome this issue, we plan to solve the upper chromosphere and transition region semi-analytically in a forthcoming paper. In the present paper we changed for now to the turbulence model that is based on local dissipation in Equations (A4) and (A5) instead of the turbulence model with the wave reflection term in Equations (11)–(14). This model is less susceptible by these numerical errors.

In Figure 3, we computed the STEREO/EUVI emission images in the three coronal bandpasses for Fe emission lines at 171, 195, and 284 Å. These LOS images are produced by assuming that the plasma is optically thin for all the considered wavelengths. In the top row, the images are for the CR2107 steady solution of our previous model (Sokolov et al. 2013) using the spatially fifth-order MP5 limiter. In the middle row, we demonstrate the new model with the MP5 limiter. The new model better captures the active region emissivity as observed by the EUVI imager (Howard et al. 2008) on board STEREO A, as shown in the bottom panels. This enhanced emissivity is due to the increased reflection rate caused by the strong magnetic fields around active regions. We note that the steady-state simulation was performed for a synoptic magnetogram, while the observation is for the time 2011 March 7 20:00 UT, and consequently, the model cannot reproduce time dependent activity during the rotation. Also the polynomial extrapolation toward the pole in the CR2107 magnetogram might distort the high latitudinal region somewhat unfavorably. The observed polar coronal holes are somewhat wider than the coronal holes of the new model. In Figure 4, we similarly plot the results for STEREO B, which shows the other side of the Sun as the two STEREO spacecraft are separated by about 177°. The emissivity of the active regions is, again, improved.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 5, we show the comparison between the model synthesized SDO/AIA images with the images observed by AIA (Lemen et al. 2012) on board SDO. The model results are obtained with the MP5 limiter. The wavelengths indicated at the top of each panel correspond to various characteristic temperatures. Again, the active regions are well captured.

Zoom In Zoom Out Reset image size

We have presented our new AWSoM model. This solar model, which is part of the SWMF, is a 3D Alfvén wave turbulence-driven model ranging from the upper chromosphere to the whole heliosphere. Compared to our previous models, AWSoM includes a generalization of the Alfvén wave turbulence to counter-propagating waves on both open and closed field lines. The outward propagating waves are now partially reflected by the Alfvén speed gradients and field-aligned vorticity. The balanced turbulence at the apex of the closed field lines is accounted for. We have also generalized our separate electron and ion temperature to anisotropic ion temperatures and isotropic electron temperatures. To distribute the turbulence dissipation to the coronal heating of the three temperatures, we use the results of the linear wave theory and nonlinear stochastic heating as presented in Chandran et al. (2011). For the isotropic electron temperature, we have now also incorporated the collisionless heat conduction.

Our new model has no ad hoc coronal heating functions and has only a few adjustable parameters: three to prescribe the boundary conditions (density, temperature, and Poynting flux of the Alfvén waves), a transverse correlation length parameter for the turbulence and heat partitioning, a parameter related to the nonlinear stochastic heating of the ions, and two parameters for the collisionless heat conduction. Some of these parameters could potentially be described more self-consistently. For example, the transverse correlation length could be obtained from a time evolution equation (Breech et al. 2008) instead of the simple scaling with the magnetic field strength.

Since the evolution equations of our model do not assume open or closed field lines, those will develop self-consistently by using the data from photospheric magnetic field observations as boundary conditions for the magnetic field. The correctness of the coronal heating can be tested by comparing the simulated and observed multi-wavelength EUV images. We performed such a validation for CR2107. We demonstrated that our model can reproduce many features seen in the LOS images. The high latitudinal region is somewhat distorted. This might be an artifact due to the polynomial interpolation of the synoptic magnetogram above 60° toward the pole. Future improvements in adapting magnetograms might address this issue. In our companion paper, X. Meng et al. (in preparation), we will showcase the model performance in the IH by comparing the results for two CRs with in situ observations at 1 AU.

This work was supported by the NSF grant AGS 1322543. W.B. Manchester IV was supported by NASA grant NNX13AG25G. The simulations were performed on the NASA Advanced Supercomputing system Pleiades.

We now consider the case that the turbulence is strongly imbalanced, i.e., min (w_±)/max (w_±) ≪ 1. For simplicity, we also assume that w₊ is the dominant wave, and hence, w₋ is the minor wave (w₋ ≪ w₊). We can then simplify the wave equations (36) as

$$\begin{eqnarray} &&\frac{\partial w_\pm }{\partial t} + \nabla \cdot \left[({\bf u}\pm {\bf V}_A)w_\pm \right] + \frac{w_\pm }{2}(\nabla \cdot {\bf u})\nonumber\\ &&\quad= \mp {\cal R}_{\rm imb}\sqrt{w_-w_+} -\Gamma _\pm w_\pm. \end{eqnarray} \tag{ A1 }$$

For the minor wave equation, only the right hand side of Equation (A1) has terms with the dominant wave. Hence, the reflection and dissipation term of the minor wave equation can, in leading order of the small quantity w₋/w₊, be assumed to be balanced. Using the dissipation rate $\Gamma _-= (2/L_\perp) \sqrt{w_+/\rho }$, the minor wave energy density can then analytically be determined as $w_- = ({1}/{4})\rho L_\perp ^2 {\cal R}_{\rm imb}^2$. Using this expression in the dissipation rate $\Gamma _+=(2/L_\perp)\sqrt{w_-/\rho }$ of the dominant wave equation and further noting that, to leading order, the reflection term on the right hand side of Equation (A1) is much smaller than the dissipation term for the dominant wave, we arrive at the dominant wave equation

$$\begin{equation} \frac{\partial w_+}{\partial t} + \nabla \cdot \left[({\bf u}+{\bf V}_A)w_+\right] + \frac{w_+}{2}(\nabla \cdot {\bf u})= -{\cal R}_{\rm imb} w_+, \end{equation} \tag{ A2 }$$

valid for strongly imbalanced turbulence. We note that the resulting equation only depends on the dominant wave energy density. A similar derivation can be performed when the w₋ is the dominant wave, and we combine both cases in a single formulation for strongly imbalanced turbulence:

$$\begin{equation} \frac{\partial w_\pm }{\partial t} + \nabla \cdot \left[({\bf u}\pm {\bf V}_A)w_\pm \right] + \frac{w_\pm }{2}(\nabla \cdot {\bf u})= -{\cal R}_{\rm imb} w_\pm. \end{equation} \tag{ A3 }$$

The exponentially small minor wave energy density is now assumed to be zero due to the absence of a reflection term in these evolution equations (although it could be recovered from the aforementioned analytical expression for the minor wave). Hence, by comparing Equations (A3) and (25), we can conclude that in strongly imbalanced turbulence, the dissipation rate is equal to the reflection rate, i.e., the wave dissipation is local. The dissipation does, in addition, no longer depend on the perpendicular correlation length. However, this derivation still dismisses the case that near the apex of closed field lines the wave energy densities can be of equal amplitude (w₊ ≈ w₋). For the balanced turbulence, the dissipation should still be estimated by the original expression (26). Combining both cases results in the final expression for the wave propagation and local dissipation

$$\begin{equation} \frac{\partial w_\pm }{\partial t}+\nabla \cdot \left[({\bf u}\pm {\bf V}_A)w_\pm \right]+\frac{w_\pm }{2}(\nabla \cdot {\bf u})= -\Gamma _\pm w_\pm, \end{equation} \tag{ A4 }$$

in which the dissipation rate is

$$\begin{equation} \Gamma _\pm = \max \left({\cal R}_{\rm imb}, \frac{2}{L_\perp }\sqrt{\frac{w_\mp }{\rho }} \right). \end{equation} \tag{ A5 }$$

This approximation is valid for strongly imbalanced and balanced turbulence. The applicability to moderately imbalanced turbulence, for example in the transition region and chromosphere, is less certain.

The partial reflection of the Alfvén waves due to Alfvén speed gradients and field-aligned vorticity generates counter-propagating waves. The nonlinear interaction between these oppositely directed waves results in an energy cascade from the large scale L_⊥ through the inertial range to the smaller perpendicular scales, i.e., larger perpendicular wavenumber k_⊥, where it can dissipate. The apportioning of the dissipated energy to the coronal heating functions Q_e, Q_i, and Q_i∥ depends on the microphysics that is involved. In this paper, we follow the partitioning strategy based on the dissipation of kinetic Alfvén waves (KAWs) using the theory described in Chandran et al. (2011). That formalism has the distinct advantage of providing approximated formulas that can readily be implemented in a numerical solar wind model that is based on the turbulent cascade of Alfvén waves. We have implemented those formulas in the present model, and below, we reproduce them for convenience.

In Chandran et al. (2011), the cascading of Alfvén waves transitions into cascading of KAWs. The KAWs can dissipate when k_⊥r_i ∼ 1, where r_i is the ion gyro radius. Among the dissipation mechanisms considered are the linear Landau damping and linear transit time damping of KAWs, which contribute to electron and parallel ion heating. The corresponding damping rates Γ_e and Γ_i∥ are

$$\begin{eqnarray} &&\Gamma _e t_c = 0.01\left(\frac{P_e}{P_i\beta _i}\right)^{1/2}\left[ \frac{1+0.17\beta _i^{1.3}}{1+\left(2800\beta _e\right)^{-1.25}}\right], \end{eqnarray} \tag{ B1 }$$

$$\begin{eqnarray} &&\Gamma _{i\parallel } t_c = 0.08\left(\frac{P_e}{P_i}\right)^{1/4} \beta _i^{0.7} \exp \left(-\frac{1.3}{\beta _i}\right), \end{eqnarray} \tag{ B2 }$$

where β_e = 2μ₀P_e/B² and β_i = 2μ₀P_i/B² are the electron and averaged ion plasma beta. Similar to Chandran et al. (2011), the Alfvén frequency 1/t_c = k_∥V_A for the parallel wavenumber k_∥ can be rewritten as $t_c = \rho \delta v_i^2/(\Gamma _+w_++\Gamma _-w_-)$ under the assumption of the critical-balance condition. The velocity perturbation δv_i of the Alfvén waves and KAWs at k_⊥r_i ∼ 1 is assumed to scale with r_i/L_⊥ via

$$\begin{equation} \rho \delta v_i^2 \approx w_d \sqrt{\frac{r_i}{L_\perp }}, \end{equation} \tag{ B3 }$$

where w_d = max (w₊, w₋) is the dominant wave energy density. The minor wave energy density, w_m = min (w₊, w₋), is assumed to be exponentially small compared to w_d, which is, strictly speaking, not true in the balanced turbulence regime and hence, introduces some uncertainty in the heating partitioning. The above scaling is compatible with the 1 AU observations of Podesta et al. (2007). Furthermore, we assume similar to Chandran et al. (2011), nonlinear damping of KAWs via stochastic heating of ions, resulting in perpendicular ion heating. This energization is effective if δv_i is large enough (Chen et al. 2001; Johnson & Cheng 2001). This form of heating is the result of stochastic ion orbits perpendicular to B in an electrostatic potential. The damping rate for β_i∥ = 2μ₀P_i∥/B² ≲ 1 is

$$\begin{equation} \Gamma _{i\perp } = 0.18 \varepsilon _i\Omega _i \exp \left(-\frac{h_S}{\varepsilon _i}\right), \end{equation} \tag{ B4 }$$

where Ω_i = (e/m_i)B is the ion gyro frequency, $v_{i\perp } = \sqrt{2P_{i\perp }/\rho }$ is the perpendicular ion thermal speed, ε_i = δv_i/v_i⊥, and h_S is an input parameter for the stochastic heating. The heating functions are expressed in terms of the damping rates:

$$\begin{eqnarray} && Q_e = \frac{1+\Gamma _et_c}{1+(\Gamma _e+\Gamma _{i\parallel } +\Gamma _{i\perp })t_c}(\Gamma _+w_++\Gamma _-w_-), \end{eqnarray} \tag{ B5 }$$

$$\begin{eqnarray} && Q_{i\parallel } =\frac{\Gamma _{i\parallel }t_c}{1+(\Gamma _e+\Gamma _{i\parallel }+\Gamma _{i\perp })t_c} (\Gamma _+w_++\Gamma _-w_-), \end{eqnarray} \tag{ B6 }$$

$$\begin{eqnarray} && Q_i = \Gamma _+w_++\Gamma _-w_- - Q_e. \end{eqnarray} \tag{ B7 }$$

The 1 + term in these expressions is for the remaining cascading power that succeeds to cascade to k_⊥r_i ≫ 1, so that it can be assumed to be dissipated via interactions with electrons and hence, contributes to electron heating.

In this Appendix, we will derive the final form of our implemented collisionless electron heat conduction. For convenience, we will limit the derivation to the IH only, where the collisional heat conduction, Coulomb collisional heat exchange with the ions and the radiative cooling, can be neglected in the full electron energy equation (5). We additionally omit the time derivate as we focus in this paper on the steady-state solar wind, and hence, Equation (5) can be simplified as

$$\begin{equation} \nabla \cdot \left(\frac{p_e}{\gamma - 1}{\bf u}\right)+p_e\nabla \cdot {\bf u} = - \nabla \cdot \left[ \frac{3}{2} \alpha p_e{\bf u}\right] + Q_e. \end{equation} \tag{ C1 }$$

The first term on both the left hand side and the right hand side can be combined, resulting in

$$\begin{equation} \nabla \cdot \left(\frac{p_e}{\gamma _H - 1}{\bf u}\right) + p_e\nabla \cdot {\bf u}= Q_e, \end{equation} \tag{ C2 }$$

where we have introduced a new polytropic index γ_H for the electrons in the collisionless regime:

$$\begin{equation} \gamma _H = \frac{\gamma + \frac{3}{2}(\gamma -1)\alpha }{1 + \frac{3}{2}(\gamma -1)\alpha }. \end{equation} \tag{ C3 }$$

For our standard values α = 1.05 (taken from Cranmer et al. 2009) and γ = 5/3, we obtain γ_H ≈ 1.33. By reintroducing the missing terms of Equation (5), we obtain our final evolution equation for the electron pressure:

$$\begin{eqnarray} &&\frac{\partial }{\partial t}\left(\frac{P_e}{\gamma _e-1}\right) +\nabla \cdot \left(\frac{P_e}{\gamma _e-1}{\bf u}\right)+P_e\nabla \cdot {\bf u}\nonumber\\ &&\quad= -\nabla \cdot {\bf q}_e^* +\frac{N_ik_B}{\tau _{ei}}\left(T_i-T_e\right)-Q_{\rm rad}+ Q_e, \end{eqnarray} \tag{ C4 }$$

where

$$\begin{equation} \gamma _e = \gamma f_S + \gamma _H (1-f_S) \end{equation} \tag{ C5 }$$

interpolates the electron polytropic index γ_e between the collisional regime where γ_e = γ and the collisionless regime where γ_e = γ_H and the interpolation function f_S is taken to be the same as Equation (10). The * in ${\bf q}_e^*$ indicates that we set q_{e, H} = 0 in the electron heat flux (9), i.e., ${\bf q}_e^* = f_S {\bf q}_{e,S}$, since it is now parameterized via a spatially varying γ_e.

The main difference between Equation (C4) and (5) is the use of γ_e instead of γ in the time derivative, and hence, the time evolution of both (ad hoc) formulations is different. As a final note, we mention that a spatially varying electron polytropic index does not negatively impact the shock evolution of CMEs. Since the electron speed of sound is much larger than the CME speeds, only the ions will be heated by the CME shock. The ion fluid still uses the standard polytropic index γ = 5/3.