MAGNETOHYDRODYNAMIC WAVES AND CORONAL HEATING: UNIFYING EMPIRICAL AND MHD TURBULENCE MODELS

In Usmanov et al. (2000), a three-dimensional (3D) model for the solar wind was suggested in which the Alfvén wave turbulence pressure served to accelerate the solar wind. The solar wind bimodal structure as observed by Ulysses was successfully reproduced in the numerical simulation. However, the long history of quantitative agreement of this model with solar wind observations at 1 AU is insufficient for global space weather simulations. The same criticism seems to be applicable to the more refined and physics-based Alfvén-wave-driven models of the solar wind (Suzuki & Inutsuka 2005; Verdini et al. 2010; Osman et al. 2011).

Therefore, the semi-empirical approach so far has been better suited for global space weather simulations. The most popular parameterization was adopted by Arge & Pizzo (2000) in their Wang–Sheeley–Arge (WSA) model, which well describes the solar wind parameters at 1 AU. The semi-empirical "synoptic" formulae for the solar wind speed employ the solar magnetogram and the properties of the magnetic lines of the potential magnetic field as recovered from the synoptic magnetogram data. The empirical dependence of the solar wind properties from two input parameters, θ, which is the angular distance from the solar wind origin point to the nearest coronal hole boundary, and f_exp, which is the expansion factor, $f_{\rm exp}=|{\bf B}_{R = R_\odot }|R_s^2/[|{\bf B}_{R=2.5\,R_\odot }|(2.5R_s)^2]$, is derived from the following two assumptions and observations. First, the slow solar wind is assumed to originate from the coronal hole boundary (small values of θ), whereas the fast solar wind originates from the central part of the coronal hole (large values of θ). Second, small coronal holes (having large values of f_exp) usually produce slower solar wind compared to large coronal holes (having small values of f_exp).

In Cohen et al. (2007), we coupled the semi-empirical WSA formulae to the global 3D model for the solar corona and inner heliosphere within the Space Weather Modeling Framework (SWMF; Tóth et al. 2005). The WSA formulae were used as the boundary condition for the model that had been coupled to the MHD simulator via the varied polytropic index distribution (see Roussev et al. 2003b). However, the physics of Alfvén wave turbulence has almost no intersection with the semi-empirical model of WSA.

In order to be competitive with semi-empirical models of the solar wind, the turbulence-driven model should quantitatively reproduce the solar wind variation from the coronal hole boundary to the coronal hole central part, similar to that parameterized by the expansion factor.

Proceeding from the solar wind to the coronal heating, we can see again the disconnection between the physics-based models for the turbulent heating in the corona, one of the most advanced models of this kind being recently described in Cranmer (2010) (see also Tu & Marsch 1997; Hu et al. 2000; Li & Habbal 2003; Dmitruk et al. 2002), on one hand and well-established models with semi-empirical heating function, such as that presented in Lionello et al. (2001, 2009), Riley et al. (2006), Titov et al. (2008), and Downs et al. (2010), on the other hand. The heating function is applied to power the plasma in the solar corona with some heating rate, the distribution of which is given by some functions chosen not from some deep physical considerations, but in an ad hoc manner to provide the better agreement with observations.

These heating functions should be applied in conjunction with thermodynamic energy equation(s) in the lower corona (LC) model. A direct accounting is needed of the transition region between the chromosphere and the corona, where non-ideal-MHD thermodynamic terms of energy transport such as electron heat conduction, radiative losses, and coronal heating all become important. In Downs et al. (2010), we used this thermodynamic MHD model to explore empirical parameterizations of coronal heating in the context of realistic 3D magnetic structures observed in the EUV and soft X-rays on 1996 August 27. Through direct comparison of synthetic observables to observations we demonstrated that this model can effectively capture the interplay between coronal heating and electron heat conduction.

The plasma parameter distribution obtained in this way may be successfully applied to generate synthetic EUV and X-ray images, which appear to be in a good agreement with those obtained with the EIT telescope on board Solar and Heliospheric Observatory (SOHO) and SXT on board Yohkoh (up to 2001). The observation synthesis capability has been extended to the major low corona imaging instruments available in space, namely, STEREO/EUVI, SDO/AIA, and Hinode/X-Ray Telescope. The best agreement with observations is achieved with ad hoc heating functions (such as the unsigned-flux-based heating model as discussed in Section 2.1). In order to be competitive with ad hoc coronal heating models, the turbulence-driven model should quantitatively reproduce some successful heating function, that is, the realistic wave dissipation should provide the same heating rate as, for example, the unsigned-flux-based heating model (Abbett 2007).

In addition to the mentioned disconnection between physics-based and observation-driven models, there are two more contradictions to discuss. First, the quantitative Alfvén-wave-driven models for the solar wind and for the coronal heating do not well conform with each other. For example, the coronal heating model of Cranmer (2010), once applied to a realistic 3D global model, does not give realistic plasma parameters in the solar wind region (in the coronal holes). Second, the increased estimates for the Alfvén wave turbulence energetics in the Sun's proximity require us to revisit the models for the evolution of turbulence while the solar wind propagates toward 1 AU. The revisited model should account for both the observed level and frequency spectrum of turbulence at 1 AU and the solar wind ion temperature resulting from the turbulence dissipation.

With the advent of modern computational tools it is now becoming the norm to employ detailed 3D computer models as simulation tools that directly account for the inhomogeneous nature of the Sun–heliosphere environment. The key advantage of this approach is the ability to compare and validate model results through direct comparisons to all kinds of the observational data listed above: for the solar wind and turbulence characteristics at 1 AU (see also Jin et al. 2012) and for EUV and X-ray images of the solar corona.

In the current paper, we get rid of the ad hoc heating functions and parameterize the coronal heating in the LC in terms of the Alfvén wave turbulence dissipation.

The theoretical model for this approach is summarized in Section 2. In Sections 2.1 and 2.3 we demonstrate the possibility to parameterize a popular and successful semi-empirical heating model based on unsigned magnetic flux in terms of the Alfvén wave turbulence dissipation. While choosing the way to parameterize heating, one may try to vary both the boundary condition for the wave energy flux (the Poynting flux) inflowing to the solar corona from the solar surface and the dissipation length for turbulence. In Section 2.2, we show that the boundary condition is strongly restricted, which reduces the model uncertainty. With the pre-specified boundary condition, the drastic difference in the plasma heating mechanism between open and closed field regions should be caused by a difference in the wave dissipation efficiency. In Section 2.4, we suggest and describe the physical mechanism of a nonlinear interaction between oppositely propagating waves. The degraded intensity of the inward propagating waves may be responsible for the reduction in the turbulence dissipation rate in the coronal holes, thus resulting in the bimodal solar wind structure.

In Section 3, we summarize a computational model and the code we use to simulate the state of the solar corona. Section 4 presents the simulation results for CR2107 and their comparison with EUV images. In Section 5 we discuss the plans for a future.

Among the different volumetric heating functions examined in Downs et al. (2010), the method previously adopted by Abbett (2007) was most successful at reproducing both EUV and soft X-ray observations. This heating function is based on the scaling law obtained by Pevtsov et al. (2003). This law establishes the power-law relationship between the total heating power, E = ∫edV, integrated over the plasma volume above 1 R_☉, and the unsigned magnetic flux, $\Phi = \int {|{\bf B}_{R_\odot }\,{\boldsymbol \cdot}\, {\bf dS}|}$, integrated over the photospheric surface, i.e.,

$$\begin{equation} \int _{R \ge R_\odot }{e dV} = {\rm const} \, \Phi , \quad \Phi = \int _{R=R_\odot }{|{\bf B} \,{\boldsymbol \cdot}\, {\bf dS}|}. \end{equation} \tag{ 1 }$$

Here R is the heliocentric distance and R_☉ is the solar radius. In the original work by Pevtsov et al. (2003), the relationship was found not to be linear: E_X ≈ 0.894Φ^1.1488 [CGSE]. Note, however, that the observational data used to derive the scaling law were the X-ray total luminosity, E_X, which is only a small fraction of the total heating power, E. There are also losses due to the electron heat conduction, radiation in spectral ranges other than X-rays, and also solar wind expansion. As a result, a scaling factor that relates E_X to E is required, ζ = E_X/E ≈ (1÷2) × 10⁻². If one assumes that the E_X/E ratio is somewhat elevated during solar maximum, or near active regions where the average magnetic field strength, |B|, is higher, so that ζ∝Φ^0.1488, then we can state that Equation (1) is equivalent to that presented in Abbett (2007) and Pevtsov et al. (2003). On one hand, this relationship quantifies the well-known fact that during solar maximum, or near active regions, the average B is higher. On the other hand, the coronal heating is more intense in this case.

The constraint given by Equation (1) is not sufficient to establish the 3D distribution of the heating function. In Abbett (2007), the heating function was scaled with the magnetic field: e∝|B|, with the constant factor being chosen to satisfy a power-law scaling similar to Equation (1). In Downs et al. (2010), an exponential envelope has been adopted: e∝|B|exp [ − (R − R_☉)/L], where the dissipation length was chosen to be L ≈ 40 Mm (herewith, R is the heliocentric distance). The latter formula, however, appears to be applicable only for closed field regions away from active regions. In coronal holes, in turn, the dissipation length was chosen to be greater by a factor of 10. In all cases, the sophisticated envelope function had to be integrated over the volume in order to obtain the common constant factor satisfying Equation (1).

Below we demonstrate that a model of coronal heating that both satisfies Equation (1) and provides the spatial distribution of the heating envelope similar to that discussed above can be realized assuming Alfvén waves as the main heating source of the solar corona. Without going into details about the wave absorption and nonlinear conversion and transport, we can construct the heating function by (1) imposing a boundary condition for the Poynting flux of Alfvén wave energy and (2) adopting an absorption mechanism for the waves that would result in the desired distribution of energy deposition from waves into the coronal plasma.

Assuming that the Alfvén wave turbulence is the main source for heating the coronal plasma, here we discuss a choice of boundary condition for their Poynting vector, P, at the solar surface. We consider that the local value of the Pointing flux at the photosphere, $P_{R_\odot }$, scales with the magnetic field magnitude, $| {\bf B}|_{R_\odot }$, as

$$\begin{equation} P_{R_\odot } = {\rm const} \,| {\bf B}|_{R_\odot }. \end{equation} \tag{ 2 }$$

With this choice for P, Equation (1) is automatically satisfied if one assumes that the entire energy of the Alfvén waves traveling from the photosphere to the solar corona, $\int _{R = R_\odot }{({P_{R_\odot }}/{|{\bf B}|_{R_\odot }})|{\bf B}\cdot {\bf dS}|}$, is deposited into the coronal plasma, $\int _{R\ge R_\odot }{e\,dV}$. Then

$$\begin{eqnarray*} E&=&\int _{R \ge R_\odot }{e\,dV} = \int _{R=R_\odot } {\frac{P_{R_\odot }}{|{\bf B} |_{R_\odot }}|{\bf B}\,{\boldsymbol \cdot}\, {\bf dS}|} \nonumber\\ &=& \frac{P_{R_\odot }}{|{\bf B}|_{R_\odot }} \int _{R=R_\odot }{|{\bf B}\,{\boldsymbol \cdot}\, {\bf dS}|} = \frac{P_{R_\odot }}{|{\bf B}|_{R_\odot }}\Phi , \end{eqnarray*}$$

meaning that Equation (1) is fulfilled and that the exact same constant is present in Equations (1) and (2):

$$\begin{equation} \frac{E}{\Phi } = \frac{P_{R_\odot }}{|{\bf B} |_{R_\odot }}. \end{equation} \tag{ 3 }$$

However, the relationship given by Equation (1) is not the only way to arrive at Equation (2). It can be also thought of as the percolation boundary condition for the Poynting flux of the waves. Assume that all the magnetic field in the solar corona and all the Alfvén wave turbulence percolate to the solar corona from a depth at the solar interior where both the magnetic field intensity (presumably 50–70 kG) and the Alfvén wave Poynting flux are very high compared to their photospheric values (except in sunspots). In this case, both the magnetic field flux and the wave energy flux "propagate" along the flux tube, and they both vary inversely proportional to the flux tube cross section. (Note that the growth of the flux tube cross section from the solar interior toward the solar surface may be quite large and is hardly predictable.) As a consequence, the local value of the magnetic field at the photosphere is orders of magnitude smaller than that in the solar interior, and its variation along the solar surface may be quite large. The same behavior should be true for the Poynting flux. However, the ratio of the Poynting flux to the magnetic field intensity remains constant along the magnetic tube, and that is exactly what is assumed in Equation (2).

Surprisingly, in addition to these two considerations, we can cite four different groups of authors that arrive at the same condition as in Equation (2) by quite different reasons (and stemming from quite different observations). First, the assumption as in Equation (2) applied to the radial components of the Poynting flux, P_R, and the magnetic field, B_R, is a keystone of the Fisk theory for the solar wind (see Fisk 1996, 2001; Fisk et al. 1999a, 1999b; Fisk & Schwadron 2001). Then, Farrugia et al. (1997) (and the works cited therein) formulated an "abnormal" adiabatic expansion law relating the (turbulent) energy density, w, to the mass density, ρ: w∝ρ^1/2 for the coronal matter, based on numerous observations relating to the coronal plasma and the coronal mass ejection (CME) ejecta. The suggested adiabatic law has an abnormal value of the adiabatic index: w∝ρ^γ, γ = 1/2 < 1. Here, assuming a plasma dominated by turbulence, ρ is the mass density and w could be the energy density of the Alfvén waves (w ∼ (δB)², with δB being the irregular magnetic field). As long as for the Alfvén waves P = V_Aw and the Alfvén wave speed scales with the density and magnetic field as $V_A=({|{\bf B}|}/{\sqrt{\mu _0\rho }}) \propto B\rho ^{-1/2}$, Equation (2) may be also rewritten in the form of the "abnormal" adiabatic law:

$$\begin{equation} w_{R_\odot }={\rm const}\sqrt{\mu _0\rho }=\frac{P}{|{\bf B}|}\sqrt{\mu _0\rho }. \end{equation} \tag{ 4 }$$

As long as media with adiabatic index less than one are not thermodynamically stable, it is hardly instructive to interpret Equation (2) in this manner. However, this unexpected support for Equation (2) is worth mentioning here.

Another example supporting Equation (2) is the paper by Suzuki (2006), which stems from the WSA model and proposes a semi-empirical quantitative model for the solar wind. In their work, it is assumed that the solar wind is mostly powered by the Alfvén wave turbulence, and the input parameter in the model is the constant value of the following average: 〈δB_⊥ · δv_⊥〉 ≈ 0.83 T m s⁻¹ at the solar surface. Here, δB_⊥ and δv_⊥ are the turbulent magnetic field and the turbulent velocity pulsation, respectively. These are both orthogonal to the regular magnetic field. It is straightforward to demonstrate that the above assumption is equivalent to Equation (2):

$$\begin{equation*} \langle \delta {\bf B}_\perp \cdot \delta {\bf v}_\perp \rangle =\frac{|{\bf B}|}{V_A}\frac{w}{\rho }=\frac{|{\bf B}|^2}{\rho V_A^2}\frac{V_A w}{|{\bf B}|}= \mu _0\frac{P_{R_\odot }}{|{\bf B}|_{R_\odot }}. \end{equation*}$$

From here, one can see that the Alfvén wave turbulence may be employed to reproduce the observed solar wind parameters with the use of a boundary condition similar to that given by Equation (2). In this case, the constraint for the constant factor in Equations (1) and (2) is $(P/|{\bf B}|)_{R_\odot } \ge (0.83\,{\rm T\,m\,s^{-1}})/\mu _0\approx 7\times 10^5\,{\rm W\,(m^2\, T)^{-1}}$. This is the lower bound as long as in the model of Suzuki (2006) there is no energy loss mechanism from the solar wind plasma (radiation, heat conduction, etc.). Therefore, the realistic heating rate should be higher in order to balance these losses. Nevertheless, the estimate is very close to that which follows from Pevtsov et al. (2003), Abbett (2007), and Downs et al. (2010). Specifically, for the choice ζ = 2 × 10⁻²〈|B|[G]〉^0.1488, where the average magnetic field intensity over the solar surface is introduced, the estimate of the Pointing-flux-to-field ratio is as follows:

$$\begin{eqnarray*} \frac{E}{\Phi }&=&\frac{E_X}{\Phi \zeta }\!\left(\!{\approx} \frac{0.894 \Phi ^{1.1488}}{2\times 10^{-2}|{\bf B}|^{0.1488}\Phi }\approx 44.7\left(4\pi R_\odot ^2\right)^{0.1488}{\rm [CGSE]}\right)\nonumber\\ &\approx& 1.1\times 10^6\frac{{\rm W}}{{\rm m^2T}}, \end{eqnarray*}$$

or, in all equivalent forms,

$$\begin{eqnarray} P&\approx& 1.1\times 10^6\frac{\rm W}{\rm m^2}\left[\frac{|B|}{1\,T}\right],\nonumber\\ P_R&\approx& 1.1\times 10^6\frac{\rm W}{\rm m^2}\left[\frac{|B_R|}{1\,T}\right],\nonumber\\ w&=&\frac{P}{V_A}\approx 6.6\times 10^{-3}\,{\rm W\,m^{-3}}\sqrt{\frac{\rho }{3\times 10^{-11}\,{\rm kg\,m^{-3}}}},\nonumber\\ \langle \delta {\bf B}_\perp \cdot \delta {\bf v}_\perp \rangle &\approx& 1.4\,{\rm T\,m\,s^{-1}},\nonumber\\ \langle \delta {\bf v}_\perp \cdot \delta {\bf v}_\perp \rangle &\approx& (15\,{\rm km\,s^{-1}})^2\sqrt{\frac{3\times 10^{-11}\,{\rm kg\,m^{-3}}}{\rho }}, \end{eqnarray} \tag{ 5 }$$

where we chose the value of mass density 3 × 10⁻¹¹ kg m⁻³ (N_a ≈ 2 × 10¹⁶ m⁻³) corresponding to the top of the chromosphere (near the transition region), to compare with the Hinode observations of the turbulent velocities ≈15 km s⁻¹. Thus, a variety of models, which are validated against EUV observations, X-ray radiance measurements, and in situ solar wind observations, are all in agreement with Equation (2), with a rather narrow range for the constant factor in the expression.

Next, let us compare this boundary condition with a coronal heating model based on Alfvén wave turbulence, such as the advanced model of Cranmer (2010) (see also the works cited therein). By comparing Equation (2) with Equation (29) from Cranmer (2010), one can see that Equation (2) is adopted in their model not just as a boundary condition at the solar surface, but as a fundamental relationship valid throughout the solar corona (or, at least, for heliocentric distances 1.1 R_☉ ⩽ R ⩽ 1.5 R_☉). The observational constraints for the constant factor in Equation (2) are somewhat weaker than the values discussed here (P/|B| ∼ (0.5÷1.0) × 10⁶ W (m² T)⁻¹). However, the deviation is very small. In addition, there is some uncertainty in these observations (e.g., the Pointing flux is not measured, but rather the oscillating velocity). Also, the solar observations at 1.1 R_☉ cannot be easily converted into a boundary condition at the photosphere. Furthermore, the dissipation length of 40 Mm assumed by Downs et al. (2010) is shorter than the difference between 1.1 R_☉ (where the oscillations are observed) and 1 R_☉ (where the boundary condition is imposed). As a consequence, it is likely that the turbulence observed at 1.1 R_☉ is already weakened by the wave absorption.

In order to proceed to the 3D distribution of the heating function, we represent the volume integral of the heating function as the total of integrals over the flux tubes:

$$\begin{equation*} E=\int {\left(\int {edS(\ell)}\right)d\ell }, \end{equation*}$$

with dℓ being the length differential along the flux tube and dS(ℓ) being the expanding cross section of the flux tube. Along the flux tube, the magnetic flux is conserved, and hence dS(ℓ)|B|(ℓ) = const = |(B · dS)|. Again, assuming that the Alfvén waves propagate only along the flux tube, their energy conservation is given by the following relationship: P(ℓ + dℓ)dS(ℓ + dℓ) − P(ℓ)dS(ℓ) = −edℓdS(ℓ); hence,

$$\begin{equation} \frac{e(\ell)}{|{\bf B}|(\ell)}=-\frac{\partial }{\partial \ell }\left(\frac{P(\ell)}{|{\bf B}|(\ell)}\right). \end{equation} \tag{ 6 }$$

By introducing the absorption length, L, such that e = V_Aw/L = P/L, one can integrate Equation (6) along the magnetic field line:

$$\begin{equation*} P(\ell)=\left(\frac{P}{|{\bf B}|}\right)_{R=R_\odot }|{\bf B}|(\ell)\exp \left({-}\int _0^\ell {\frac{d\ell }{L}}\right), \end{equation*}$$

$$\begin{equation} e(\ell)=\left(\frac{P}{|{\bf B}|}\right)_{R=R_\odot }\frac{|{\bf B}|(\ell)}{L}\exp \left({-}\int _0^\ell {\frac{d\ell }{L}}\right). \end{equation} \tag{ 7 }$$

Depending on the choice of the absorption length, one can obtain the following spatial distribution of the heating function: (1) for l ≪ L, the heating function is the same as that used in Abbett (2007), i.e., e∝|B|; (2) for l ≫ L, the heating function decays exponentially, i.e., e∝exp (− ∫dℓ/L). We now arrive at an important conclusion, that is, the most common heating functions may be parameterized in terms of an effective absorption coefficient for the Alfvén wave turbulence! The desired absorption coefficient, V_A/L, may be directly included in the WKB equation for total wave intensities, w_±, propagating parallel and antiparallel to the mean magnetic field vector, b = B/|B|, correspondingly (Usmanov et al. 2000):

$$\begin{equation} \frac{\partial w_{\pm }}{\partial t}+\nabla \cdot \left({\bf u}w_{\pm }\pm {\bf b}V_Aw_{\pm }\right)+\frac{1}{2}w_{\pm }(\nabla \cdot {\bf u})=-\frac{V_Aw_{\pm }}{L}. \end{equation} \tag{ 8 }$$

In the open field region (in the solar wind) only outward propagation waves were involved in the Usmanov et al. (2000) heliosphere model, that is, the "plus" waves in the regions of the solar wind with the positive radial magnetic field and the "minus" waves elsewhere. In addition to this, the low solar corona model of the present paper also includes the closed field regions, where the waves of both directions are present emitted from the opposite footpoints of a closed magnetic field line. However, Equation (8) can only be applied as a rather crude approximation. In the more refined model as derived in Sokolov et al. (2009), on the left-hand side of Equation (8) the term proportional to the velocity divergence and spectral evolution in the convergent/divergent flows should be properly accounted for. Even more important is the role played by the frequency dependence on the right-hand side of Equation (8):

$$\begin{equation} \frac{\partial I_{\pm }}{\partial t}+\nabla \cdot \left({\bf u}I_{\pm }\pm {\bf b}V_AI_{\pm }\right)-\frac{\omega }{2}\frac{\partial I_{\pm }}{\partial \omega }(\nabla \cdot {\bf u})=-\Gamma (\omega)I_\pm , \end{equation} \tag{ 9 }$$

where the equations are formulated for the wave energy density, I_±, related to the unit volume and to the interval of the wave circular frequency, dω, in the comoving frame of reference, so that w_± = ∫^∞₀I_±dω. Not only is the dependence of the dissipation coefficient, Γ(ω), on the wave frequency significant (practically so that only for the highest wave frequencies is the absorption non-negligible), but also the nonlinear wave–wave interaction may occur in the form of the wave package upshift in frequency. As a result, the wave energy from large-scale perturbations flows through the spectrum toward the short-scale end, and once transferred to the shortest possible scales, the energy dissipates.

2.3.1. Kolmogorov Turbulence and Dissipation

We found above that the unsigned-flux-based heating function with the exponential envelope function can be realized as the dissipation in the Alfvén wave turbulence. In Downs et al. (2010) we saw that the unsigned-flux-based heating function well reproduces EUV and X-ray images of the closed field regions (about the coronal holes; see below) with the choice for the dissipation length of ≈40 Mm. Here we discuss a choice for the dissipation length, L, based on the phenomenological turbulence theory, desiring that in the closed field region L would be about 40 Mm and does not vary too strongly, at least at low altitudes ⩽0.1 R_☉ ≈ 70 Mm.

Consider extra terms in the governing equations for the Alfvén, turbulent energy density, accounting for the wave-wave "cascade" rate. Only the cascade-describing terms are kept. For the closed field region we employ the key assumption about the isotropic (balanced) turbulence:

$$\begin{equation} I_+=I_-=I. \end{equation} \tag{ 10 }$$

For simplicity, assume that Γ_ww(ω) characterizes the rate of conversion of two wave "photons" of the frequency of ∼ω to a single photon with frequency ∼2ω. As long as the wave energy conserves in this processes, this may be thought of as the conservative advection in the frequency space. The positive speed of this advection is ∼ωΓ_ww(ω), as long as with the rate Γ_ww the wave frequency is upshifted by [(∼ 2ω) − (∼ ω)]. So, we have a phenomenological equation as follows:

$$\begin{equation} \frac{\partial I}{\partial t}+ \frac{\partial }{\partial \omega } \left(\omega \Gamma _{ww}(\omega)I\right)=-\Gamma (\omega) I. \end{equation} \tag{ 11 }$$

The dimensionless ratio, Γ_ww/ω, may be parameterized in terms of another dimensionless ratio, ωIμ₀/B², B being the magnetic field magnitude. For two most natural estimates, Γ_ww∝ω · (ωIμ₀/B²)∝ω(δB/B)² or $\Gamma _{ww}\propto \omega \sqrt{\omega I\mu _0/B^2}\propto \omega (\delta B/B)$, we arrive at the Kraichnan spectrum or the Kolmogorov spectrum (see Li et al. 2011) as the steady-state solution of Equation (11) throughout the inertial range of frequencies (we neglect wave absorption in this range):

$$\begin{equation*} \frac{\partial }{\partial \omega }\left[\omega ^2\left(\frac{\omega I\mu _0}{B^2}\right)^{0.5;1}I\right]=0 \end{equation*}$$

for I∝ω^{−5/3; −3/2}. We chose the assumption of the Kolmogorov spectrum, so that

$$\begin{equation*} \Gamma _{ww}\propto \frac{\omega }{V_A}\sqrt{\frac{\omega I}{\rho }}. \end{equation*}$$

On integrating Equation (11) over frequency from some frequency value within the inertial range, ω_in, to infinity, neglecting wave dissipation within the inertial range, we find

$$\begin{equation} \int _{\omega _{{\rm in}}}^\infty {\Gamma (\omega) I d\omega }\propto \left(\frac{\omega ^2I}{V_A}\sqrt{\frac{\omega I}{\rho }}\right)_{\omega =\omega _{{\rm in}}}. \end{equation} \tag{ 12 }$$

The equation is valid for any choice of ω_in within the inertial range, as long as the right-hand side is constant for I∝ω^−5/3 and Γ(ω_in) = 0. The energy of the Kolmogorov spectrum, ∝∫ω^−5/3dω, is concentrated at low frequency; therefore, at the least possible ω_in ∼ ω_min one can estimate $(\omega I)_{\omega =\omega _{{\rm min}}}\sim w$ and

$$\begin{equation*} \int _{\omega _{{\rm in}}}^\infty {\Gamma (\omega) I d\omega }\propto \frac{\omega _{{\rm min}}}{V_A}\sqrt{\frac{w}{\rho }}w. \end{equation*}$$

Following Hollweg (1986), for an isotropic Alfvén wave turbulence we accept

$$\begin{equation} \int _{\omega _{{\rm in}}}^\infty {\Gamma (\omega) I d\omega }\approx \frac{1}{L_\perp }\sqrt{\frac{w}{\rho }}w=\frac{\sqrt{\langle \delta {\bf v}_\perp \cdot \delta {\bf v}_\perp \rangle }}{L_\perp }w, \end{equation} \tag{ 13 }$$

so that the total volumetric heating rate due to the turbulent energy dissipation equals

$$\begin{equation} e=\frac{1}{L_\perp }\sqrt{\frac{w=w_-=w_+}{\rho }}(w_-+w_+). \end{equation} \tag{ 14 }$$

The introduction of L_⊥∝1/k_⊥, which is the transverse correlation length of turbulence, instead of $\omega /V_A=k_\Vert$, which is the wave vector for the Alfvén wave propagating along the magnetic field line, makes the cascade theory better suited to the nature of purely transverse Alfvén waves. For such waves, the ratio of the transverse turbulent velocity to the longitudinal wavelength could hardly characterize the rate of nonlinear wave dissipation. In more refined turbulent theory (which we do not apply and even do not review here), the spectral energy should be introduced separately for parallel and transverse wave vectors. However, within our phenomenological approach we assume that $k_\perp \,{\propto}\, k_\Vert$ and combine, once needed, the spectral energy distribution over $\omega =V_Ak_\Vert$ with the transverse correlation length, to quantify the nonlinear dissipation rate.

To close the model and to compare its predictions with the unsigned-flux-based heating function, we chose, again following Hollweg (1986), the scaling law for the transverse correlation length:

$$\begin{equation} L_\perp \propto |{\bf B}|^{-1/2},\quad L_\perp =\frac{75\,[{\rm km}\,{\rm T}^{1/2}]}{\sqrt{|{\bf B}|}}. \end{equation} \tag{ 15 }$$

In the present work we use the following range of values for the empirical constant, $L_\perp \sqrt{|{\bf B}|}$:

$$\begin{equation*} 50\,[{\rm km}\,{\rm T}^{1/2}]\le L_\perp \sqrt{|{\bf B}|}\le 100\,[{\rm km}\,{\rm T}^{1/2}]. \end{equation*}$$

In the numerical estimates below we use $L_\perp \sqrt{|{\bf B}|}\approx 75\,[{\rm km}\,{\rm T}^{1/2}]$ from Hollweg (1986). The scaling law |B|L²_⊥ = const is well compatible with the "percolation hypothesis" noticed above: while waves are propagating from the solar interior along the flux tube toward the solar surface it is natural to assume that L²_⊥ varies proportionally to the flux tube cross section.

To discuss the value and the variability of the spatial scale in the coronal heating, if the latter is related to the Alfvén wave turbulence dissipation, compare Equations (14) and (15) with the representation of the wave dissipation as in the right-hand side of Equation (8). From this comparison we find $V_A/L\,{=}\,\sqrt{w/\rho }/L_\perp$. At the corona base we can employ Equation (4) to evaluate the boundary values of w, and the estimate for L in this region reads

$$\begin{equation*} L=\frac{V_A(L_\perp \sqrt{|{\bf B}|})}{\sqrt{\frac{|{\bf B}|\left(\frac{P}{|{\bf B}|}\right)\sqrt{\mu _0\rho }}{\rho }}} \approx \sqrt{\frac{V_A}{400\,{\rm km\,s^{-1}}}}40\,{\rm Mm}. \end{equation*}$$

We see that with the choice of the above estimates for the transverse correlation length and for the boundary condition for the Poynting flux the estimate of the dissipation length at the coronal base has only weak scaling (as a square root) with the Alfvén speed. The latter does not vary too much at the coronal base, as long as $V_A\propto |{\bf B}|/\sqrt{\rho }$ and the density correlates with the magnetic field intensity, being larger in the active regions with strong magnetic fields. We also see that for a reasonable estimate for the Alfvén speed at the coronal base, V_A ∼ (200–500) km s⁻¹, the dissipation length is close to the ad hoc value L = 40 Mm as used in Downs et al. (2010). The dissipation length is not used below directly, as long as Equations (14) and (15) do not include it. However, the capability to reproduce the popular unsigned-flux-based heating model with the popular envelop function is important. We can now strengthen the conclusion above in the following manner: the choice of a boundary condition for the Poynting flux together with the estimates of the nonlinear dissipation in an isotropic Kolmogorov turbulence allows us to reproduce in detail the unsigned-flux-based ad hoc heating model for the closed field region.

It has been shown in previous studies that using strongly different heating functions for open and closed field regions is an important component of the ad hoc heating approach. In Downs et al. (2010) and Lionello et al. (2009), in the coronal holes the following heating function was applied:

$$\begin{equation*} e=5\times 10^{-7}\exp [-(R-R_\odot)/(0.7\,R_\odot)]\,{\rm erg\,cm^{-3}\,s^{-1}}. \end{equation*}$$

Comparing with the unsigned-flux-based model for the closed field region, the spatial scale of the heating function in the coronal holes is significantly longer than ∼40 Mm, the value that we used above.

Here, we discuss the possibility of implementing such a heating model via the Alfvén wave absorption. First, we do not assume that in the coronal holes the Poynting flux at the corona base is different from that in the closed field region and attribute the drastic difference in the heating functions to the difference in the dissipation rate only. Second, comparing the integral of the heating function over the coronal hole volume with the accepted value of the Poynting flux for the turbulence, one can find that only a few percent of the Poynting flux is absorbed in the coronal holes within the heliocentric distance range R_☉ ⩽ R ⩽ 2.5 R_☉.

Comparing with the closed field region, where almost all wave energy is absorbed, one can conclude that the wave energy dissipation rate in coronal holes should include a small factor of the order of a few percent compared with Equations (14) and (15).

2.4.1. Imbalanced Kolmogorov Turbulence and Its Dissipation

The drastic decrease in the wave dissipation rate in the coronal holes comes naturally if we take into account a strong turbulence anisotropy in the coronal holes: the outward propagating waves should be much more intense than the waves propagating toward the Sun. Below, having in mind that the anisotropy is a natural property of any MHD turbulence in the directional magnetic field, the turbulence in the coronal holes is referred to as imbalanced, not anisotropic.

Now we consider an imbalanced turbulence, such that I₊ ≠ I₋. As long as the real cascade physics requires the presence of oppositely propagating waves, we require that the wave–wave interaction rate should be the function of minimum of w₋, w₊ in a strongly imbalanced turbulence. To achieve this, instead of Equation (14) we use the following expression for the dissipation rate:

$$\begin{equation} e=\frac{\sqrt{|{\bf B}|}}{(L_\perp \sqrt{|{\bf B}|})} \left(\sqrt{\frac{w_+}{\rho }}w_-+\sqrt{\frac{w_-}{\rho }}w_+\right). \end{equation} \tag{ 16 }$$

For the closed field region this expression is more realistic than that we used before, as long as it correctly captures more intense heating in shorter loops, in which the intensities of the oppositely propagating waves are more uniform, as compared to longer loops, in which the wave energy for counter-propagating waves may be strongly non-uniform and imbalanced. However, for coronal holes this approximation is not sufficient as long as the WKB approximation, under which Equation (11) was derived, is not accurate enough. Within the model we use here, in the coronal holes only outward propagating waves can exist, while in reality the wave outward propagation is accompanied with some reflection, resulting in appearance of waves propagating toward the Sun. The small ratio of the reflected wave amplitude to the bulk wave amplitude is described by the reflection coefficient, C_refl. Therefore, the negligibly small min (w_±) in the expression for the dissipation rate should be floored with C²_reflmax (w_±), giving an ultimate expression for the turbulent dissipation rate we use in the present work:

$$\begin{eqnarray} e&=&\frac{\sqrt{|{\bf B}|}}{(L_\perp \sqrt{|{\bf B}|})}\nonumber\\ &&\times \left(\sqrt{\frac{\max \left(w_+,C^2_{\rm refl}w_-\right)}{\rho }}w_-+\sqrt{\frac{\max \left(w_-,C^2_{\rm refl}w_+\right)}{\rho }}w_+\right).\nonumber\\ \end{eqnarray} \tag{ 17 }$$

2.4.2. How to Derive Equation (16)?

In order to justify Equations (16) and (17), we outline the way to derive them consistently. One can employ the framework of reduced MHD, which solves the equations of motion, induction, and continuity:

$$\begin{equation*} \frac{\partial {\bf u}}{\partial t}+({\bf u}\cdot \nabla){\bf u}+ \frac{\nabla B^2}{\mu _u\rho }+\frac{\nabla (P_e+P_i)}{\rho }=\frac{({\bf B}\cdot \nabla){\bf B}}{\mu _0\rho }, \end{equation*}$$

$$\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*}$$

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

by means of representing the magnetic field and velocity vectors as sums of regular and turbulent parts, ${\bf u}\,{=}\,\tilde{\bf u}\,{+}\,\delta {\bf u},\; {\bf B}= \tilde{\bf B}\,{+}\,\delta {\bf B}$ (below tildes are omitted), and by simplifying the equations for turbulent amplitudes:

$$\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} \nonumber\\ &&\qquad=\frac{({\bf B}\cdot \nabla)\delta {\bf B}}{\mu _0\rho }+\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{ 18 }$$

$$\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\\ &&\qquad=({\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{ 19 }$$

$$\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{ 20 }$$

by assuming the incompressibility conditions: ∇ · δu = 0, B · δB = 0. The equations for the Elsasser variables, ${\bf z}_\pm =\delta {\bf u}\mp \delta {\bf B}/\sqrt{\mu _0\rho }$, are obtained as the sum Equation (18) ${\mp} {1}/{\sqrt{\mu _0\rho }}\times$ Equation (19) ${\pm} {\delta {\bf B}}/{\rho \sqrt{\mu _0\rho }}\times$ Equation (20):

$$\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*}$$

where (d_±/dt) = (∂/∂t)  +  (u  ±  V_A  +  z_∓) · ∇ and ${\bf V}_A=({\bf B}/{\sqrt{\mu _0\rho }})$. The dynamic equations for the wave energy densities, w_± = ρz² _±/4, may be obtained by multiplying the above equation by ρz_±/2 and adding Equation (20) × 3z²_±/8:

$$\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 \cdot \left[({\bf z}_\mp \cdot \nabla){\bf u}\nonumber\\ &&\qquad\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*}$$

As the first (WKB) approximation one can set z_∓ = 0 in the equations for w_± and obtain

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

This 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 P_w = 1/2(w₊ + w₋). Within the more accurate approximation, the nonlinear term ∇ · (z_∓w_±) results in the turbulent cascade and the wave energy dissipation. The dissipation rate for the wave energy density, w₊, is controlled by the amplitude of the oppositely propagating wave, $|{\bf z}_-|\sim \sqrt{w_-/\rho }$, and the correlation length, L_⊥, in the transverse (with respect to the magnetic field) direction, because in the Alfvén wave, propagating along the magnetic field, δu, δB are perpendicular to the magnetic field. Therefore, $\nabla \cdot ({\bf z}_\mp w_\pm)\sim ({1}/{L_\perp })\sqrt{({w_\mp }/{\rho })}w_\pm$ and the WKB equations with an account for nonlinear dissipation read

$$\begin{equation*} \frac{\partial w_+}{\partial t}+\nabla \cdot [({\bf u}+{\bf V}_A)w_+)+\frac{w_+}{2}(\nabla \cdot {\bf u}) =-\frac{1}{L_\perp }\sqrt{\frac{w_-}{\rho }}w_+, \end{equation*}$$

$$\begin{equation*} \frac{\partial w_-}{\partial t}+\nabla \cdot [({\bf u}-{\bf V}_A)w_-)+\frac{w_-}{2}(\nabla \cdot {\bf u})=-\frac{1}{L_\perp }\sqrt{\frac{w_+}{\rho }}w_-. \end{equation*}$$

These equations work both for the balanced and moderately imbalanced turbulence and for the balanced turbulence (w₊ = w₋ = w) and reduce to the equations by Hollweg (1986) as described above. For the 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.

However, the above equations for z_± demonstrate that even in a linear approximation the WKB approach dismisses the correlations between inward and outward propagating waves. Indeed, in the partial differential equation for z₊ the source terms present are linearly proportional to z₋ and vice versa, while in the WKB approximation these source terms describing the conversion between the oppositely going waves are omitted. The omitted correlations are important (see Tu & Marsch 1995; Dmitruk et al. 2002), and we will include them in the model as described in the forthcoming publication.

Here we parameterize this effect in order to incorporate the coronal heating in the coronal holes. The WKB approximation predicts no inward propagating waves originating from the open magnetic field lines, while in reality the inward propagating wave should arise from the partial conversion ("reflection") of the outward propagating wave due to non-WKB effects. The non-WKB generation of the inward propagating waves is parameterized via their amplitude related to that of the outward propagating waves, so that the maximum of the WKB and non-WKB wave amplitudes is used to determine the dissipation rate of the dominant outward propagating waves:

$$\begin{eqnarray} &&\frac{\partial w_+}{\partial t}+\nabla \cdot [({\bf u}+{\bf V}_A)w_+]+\frac{w_+}{2}(\nabla \cdot {\bf u})\nonumber\\ &&\qquad ={-}\frac{1}{L_\perp }\sqrt{\frac{\max \left(w_-,C^2_{\rm refl}w_+\right)}{\rho }}w_+, \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} &&\frac{\partial w_-}{\partial t}+\nabla \cdot [({\bf u}-{\bf V}_A)w_-)+\frac{w_-}{2}(\nabla \cdot {\bf u})\nonumber\\ &&\qquad={-}\frac{1}{L_\perp }\sqrt{\frac{\max \left(w_+,C^2_{\rm refl}w_-\right)}{\rho }}w_-. \end{eqnarray} \tag{ 22 }$$

2.4.3. Estimates for the Reflection Coefficient

In the present research we focus on the study of the LC and do not compare the solar wind generation and its propagation toward 1 AU. Therefore, for our present purpose the estimate

$$\begin{equation*} C_{\rm refl}={\rm const}\sim 0.01\div 0.1 \end{equation*}$$

is sufficient.

Having in mind to develop this approach in the forthcoming publications, we discuss briefly how the idea of the WSA semi-empirical model can be implemented within the Alfvén-wave-turbulence-driven model. As discussed above, for the WSA model the speed of the solar wind originating from the given magnetic field line footpoint may be found from two characteristics of the line—the expansion factor, f_exp, and the distance to the coronal hole boundary, θ:

$$\begin{equation*} u_\infty =u_\infty (\theta ,f_{\rm exp}). \end{equation*}$$

In the original version of the WSA model, for example, the solar wind speed was a function of the expansion factor only: u_∞ = A/f^δ_exp, where the unknown constants, A and δ, were chosen for better fitting the observed solar wind speed at 1 AU.

Once formulated in terms of the Alfvén wave absorption, the model can no longer employ the existing WSA formulae for the solar wind, as long as there is no simple relationship between the wave absorption and the solar wind speed. However, we can keep using the idea of the WSA model and fit the observed solar wind parameters at 1 AU by properly choosing a formula for C_refl that depends on θ and f_exp. The first choice to be tested is

$$\begin{equation} C_{\rm refl}=A\cdot f_{\rm exp}^\delta , \end{equation} \tag{ 23 }$$

where the unknown constants will be chosen to better fit the observations.

Note that there is some reasoning in favor of this approach within the wave-based model for the solar wind, as long as the expansion factor for a given magnetic field line is a good characteristic of the reflection coefficient for the waves, propagating along this line. Indeed, the expansion factor is larger for strongly bent magnetic field lines, but the reflection coefficients at such lines are also larger than on straight lines. At small distances from the coronal boundary the abrupt gradients of the plasma density also can be a reason for the increase in the reflection coefficient (see Evans et al. 2012), resulting in the generally recognized opinion that the slow wind (higher wave absorption) has its origin near the coronal hole boundary (small θ).

The model includes the standard MHD equations (non-specified denotations are as usually):

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

$$\begin{equation} \frac{\partial {\bf B}}{\partial t}+\nabla \cdot \left({\bf u}\otimes {\bf B}-{\bf B}\otimes {\bf u}\right)=0, \end{equation} \tag{ 25 }$$

$$\begin{eqnarray} &&\frac{\partial (\rho {\bf u})}{\partial t}+\nabla \cdot \left(\rho {\bf u}\otimes {\bf u}-\frac{{\bf B}\otimes {\bf B}}{\mu _0}\right)\nonumber\\ &&\qquad+\nabla \left(P_i+P_e+\frac{B^2}{2\mu _0}+\frac{w_-+w_+}{2}\right)=0, \end{eqnarray} \tag{ 26 }$$

with the full energy equations applied separately to ions,

$$\begin{eqnarray} &&\frac{\partial \left(\frac{P_i}{\gamma -1}+\frac{\rho u^2}{2}\right)}{\partial t}+\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}\,{\boldsymbol \cdot}\, {\bf B})}{\mu _0}\right]\nonumber \\ &&\qquad= -({\bf u}\cdot \nabla)P_e+ \frac{N_ik_B}{\tau _{ei}}\left(T_e-T_i\right)+\Gamma _-w_-+\Gamma _+w_+,\quad \end{eqnarray} \tag{ 27 }$$

and to electrons,

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

In addition to the standard effects, the above equations account for a possible difference in the electron and ion temperatures, the electron heat conduction parallel to the magnetic field lines,

$$\begin{equation} \kappa =\frac{{\bf B}\otimes {\bf B}}{B^2}\kappa _\Vert ,\quad \kappa _\Vert \propto T_e^{5/2}, \end{equation} \tag{ 29 }$$

the radiation energy loss from an optically thin plasma,

$$\begin{equation} Q_{\rm rad}=N_eN_i\Lambda (T_e), \end{equation} \tag{ 30 }$$

and the energy exchange between electron and ions parameterized via the relaxation time, τ_ei, as this is usually done. The Alfvén wave turbulence pressure and dissipation rate are applied in the above equations. In particular simulations, the heating due to the turbulence dissipation (see Equation (27)) may be split between electron and ions, and we can also use the total energy equation for electron and ions, to improve the computational efficiency. The turbulence propagation and dissipation are described within the WKB approximation:

$$\begin{equation} \frac{\partial w_{\pm }}{\partial t}+\nabla \cdot \left({\bf u}w_{\pm }\pm {\bf b}V_Aw_{\pm }\right)+\frac{1}{2}w_{\pm }\nabla \cdot {\bf u}=-\Gamma _\pm w_{\pm }. \end{equation} \tag{ 31 }$$

The total wave energy dissipation (i.e., the total of the right-hand side of Equation (31) taken with the opposite signs), e = Γ₋w₋ + Γ₊w₊, is included in the right-hand side of the energy equation (27) as the source term. To summarize the above consideration and derivations regarding the dissipation rate (see Equations (17), (21), and (22)), we apply it in the form

$$\begin{equation} \Gamma _\pm =\frac{\sqrt{|{\bf B}|}}{(L_\perp \cdot \sqrt{|{\bf B}|})}\sqrt{\frac{\max \left(w_\mp ,C^2_{\rm refl}w_\pm\right)}{\rho }}, \end{equation} \tag{ 32 }$$

and with the following ranges for the parameters involved:

$$\begin{equation} (L_\perp \cdot \sqrt{|{\bf B}|})=(20\div 100)\,{\rm km\,T^{1/2}},\quad C_{\rm refl}=0.01\hbox{--}0.1. \end{equation} \tag{ 33 }$$

We use the simplest equation of state for the coronal plasma with the polytropic index, γ = 5/3.

The system of governing equations is solved numerically using the BATS-R-US/SWMF code (see Section 3.3). The boundary condition for the Poynting flux (or for the intensity of the outgoing waves) is given by Equation (5) (see Figure 1). The boundary condition for the coronal magnetic field is taken from the full disk magnetogram. The boundary condition for the density and temperature requires more attention and is discussed in the following subsection.

Zoom In Zoom Out Reset image size

3.2.1. Chromosphere Boundary Condition

Here we discuss the analytical solution of the hydrostatic equilibrium and heat transfer equations with the exponential heating function, Q_h  =  Aexp (− x/L). The solution is as follows:

$$\begin{equation*} T_e=T_i=T_0,\qquad N_e=N_i=N_0 \exp \left(-\frac{m_igx}{k_B(T_e+T_i)}\right), \end{equation*}$$

$$\begin{equation} Q_h=Q_{{\rm rad}}=N_e^2\Lambda (T_0)=N_0^2\Lambda (T_0) \exp \left(-\frac{m_igx}{k_BT_0}\right). \end{equation} \tag{ 34 }$$

Here g = 274 m s⁻² is the gravity acceleration near the solar surface, the direction of this acceleration being antiparallel to the x-axis, and m_i is the proton mass. The two constants in the solution, N₀ and T₀, which are the boundary values for the density and temperature correspondingly, are unambiguously related to the amplitude and decay length of the heating function:

$$\begin{equation} L=\frac{k_BT_0}{m_ig}\approx T_0\cdot (30\,{\rm m\,K^{-1}}),\qquad A=N_0^2\Lambda (T_0). \end{equation} \tag{ 35 }$$

Note a very simple relationship for the exponential decay length for the heating function, which is half of the barometric scale length for density: 2L = L_g = k_B(T_e + T_i)/(m_ig).

The described solution satisfies the equation for the heat conduction as long as the heat transfer in the isothermic solution is absent and heating at each point exactly balances the radiation cooling. The hydrostatic equilibrium is also maintained, as long as

$$\begin{equation*} k_B\frac{\partial (N_eT_e+N_iT_i)}{\partial x}=-gN_im_i. \end{equation*}$$

The suggested solution well describes the chromosphere. The short decay length of the heating function, which is equal to ≈1.5 Mm for T₀ = 5 × 10⁴ K, may be, presumably, related to absorption of (magneto)acoustic turbulent waves, rapidly damping due to the wave-breaking effects, which result in the shock wave formation and energy dissipation. Physically, including this chromosphere heating function would imply that the temperature in the chromosphere is elevated compared to the photospheric temperatures due to some mechanism acting in the chromosphere itself. By no means can this energy be transported from the solar corona as long as the electron heat conduction rate at chromospheric temperatures is very low.

As long as we do not apply such an ad hoc short-scale heating function in the chromosphere, nor include short-scale turbulent dissipation length, we set the boundary condition for the density and the temperature at the top of the chromosphere:

$$\begin{equation} T_{ch}=T_0=5\times 10^4\,{\rm K},\qquad N_{ch}=N_0\approx 2\times 10^{16}\,{\rm m^{-3}}. \end{equation} \tag{ 36 }$$

3.2.2. Transition Region

The analytical solution for the transition region had been published many times. Here we focus on the following issues: (1) merging this solution to that for the chromosphere, (2) the validity of zero-gravity approximation, and (3) the validity of the modified heat conduction model.

The heat transfer equation for a steady-state hydrogen plasma in a uniform magnetic field reads

$$\begin{equation} \frac{\partial }{\partial s}\left(\kappa _0T_e^{5/2}\frac{\partial T_e}{\partial s}\right) +Q_h-N_e^2\Lambda (T_e)=0. \end{equation} \tag{ 37 }$$

Here Q_h = Γ(w₋ + w₊) is the coronal heating function assumed to be constant at the spatial scales typical for the transition region. Note that the coordinate is taken along the magnetic field line, not along the radial direction.

On multiplying Equation (37) by κ₀T^5/2_e(∂T_e/∂s) and by integrating from the interface to chromosphere till a given point at a temperature T_e, one can obtain

$$\begin{eqnarray} &&\left[ \frac{1}{2} \kappa _0^2 T_e^5 \left(\frac{\partial T_e}{\partial s} \right)^2 +\frac{2}{7}\kappa _0 Q_h T_e^{7/2} \right]|^{T_e}_{T_{ch}}= (N_eT_e)^2\nonumber\\ &&\quad\times \int ^{T_e}_{T_{ch}}{\kappa _0T^{1/2}\Lambda (T)dT}. \end{eqnarray} \tag{ 38 }$$

Here the product, N_eT_e, is assumed to be constant; therefore, it is separated from the integrand. For a given T_ch the only parameter in the solution is (N_eT_e). It can be expressed at any point in terms of the local value of the heating flux and the radiation loss integral:

$$\begin{equation} (N_eT_e)=\sqrt{\frac{\frac{1}{2} \kappa _0^2 T_e^5 \left(\frac{\partial T_e}{\partial s} \right)^2 +\frac{2}{7}\left(\kappa _0 Q_h T_e^{7/2}-\kappa _0 Q_hT_{ch}^{7/2}\right)}{\int ^{T_e}_{T_{ch}}{\kappa _0T^{1/2}\Lambda (T)dT}}}. \end{equation} \tag{ 39 }$$

The assumption of constant (N_eT_e) is fulfilled only if the effect of gravity is negligible. Quantitatively this condition is not trivial, as long as both the barometric scale and especially the heat conduction scale are functions of temperature. The barometric scale may be approximated as L_g(T_e) ≈ T_e · (60 m K⁻¹). The heat conduction scale, L_h, can be estimated by noticing that within the large part of the transition region the radiation losses dominate over the heating function; therefore, they are balanced by heat conduction: κ₀T^5/2_e · (Te/L²_h) ∼ Q_rad. Thus, the condition for neglecting gravity is as follows:

$$\begin{equation} L_g(T_e)\approx T_e\cdot (60\ {\rm m\,K^{-1}})\gg L_h\approx \sqrt{ \frac{\kappa _0T_e^{11/2}}{\Lambda (T_e)(N_eT_e)^2}}. \end{equation} \tag{ 40 }$$

In Figure 2 we plot temperature dependencies, L_h(T_e), L_g(T_e), for (N_eT_e) = 10²⁰ K m⁻³. We see that near the chromosphere boundary the approximation (40) works very well as long as the temperature changes with height are very abrupt. The increase in temperature till 10⁵ K occurs at a length shorter than 0.1 Mm. This estimate agrees with the temperature profile recovered from observations of chromospheric lines (see, e.g., Figure 2 and Figure 8 in Avrett & Loeser 2008). However, as the temperature further grows with the height, the gravity effect on the temperature and density profiles becomes more significant. It becomes comparable with the heat conduction effect at T_e ≈ 4.5 × 10⁵ K, which can be accepted as the coronal base temperature, so that the transition region corresponds to the temperature range from T_ch ≈ 5 × 10⁴ K to 4.5 × 10⁵ K, with a typical width being ∼10 Mm ≈ R_S/70. The transition region solution merges to the chromosphere solution with negligible jump in pressure. The merging point in the chromosphere, therefore, is at the density of (N_eT_e)/T_ch ∼ 10¹⁶ m⁻³. The short heat conduction scale at the chromosphere temperature (see Figure 2) ensures that the heat flux from the solar corona across the transition region does not penetrate toward higher densities.

Zoom In Zoom Out Reset image size

It is known that Equation (39) may be used to establish the boundary condition for density at the "coronal base"; however, we use another way to model the transition region by increasing artificially the heat conduction in the lower temperature range (see Abbett 2007). Consider a following transformation of the temperature functions in Equations (37) and (38):

$$\begin{eqnarray} &&\kappa _0\rightarrow f\kappa _0,\qquad ds\rightarrow f\ ds,\nonumber\\ &&\Gamma \rightarrow \Gamma /f,\qquad Q_{{\rm rad}}\rightarrow Q_{{\rm rad}}/f, \end{eqnarray} \tag{ 41 }$$

with a common factor, f ⩾ 1. The equations do not change in this transformation, and the only effect on the solution is that the temperature profile in the transition region becomes a factor of f wider. By applying the factor f = (T_m/T_e)^5/2 at T_ch ⩽ T_e ⩽ T_m, one can achieve that the heat conduction scale in this range is almost constant and is close to ≈2 Mm for a choice of T_m ≈ 2.2 × 10⁵ K (see Figure 2).

It should be emphasized, however, that the choice of the temperature range to apply this transformation is highly restricted by the condition as in Equation (40). In choosing a higher value of T_m the heat conduction scale at the chromospheric temperature exceeds the barometric scale in the chromosphere, resulting in physically meaningless penetration of the coronal heat to the deeper chromosphere. The global model of the solar corona with this unphysical energy sink suffers from the reduced values of the coronal temperature and produces a visible distortion in the EUV and X-ray synthetic images. Thus, in formulating the transition region model we modify the heat conduction, the radiation loss rate, and the wave dissipation rate, and the maximal temperature for this modification does not exceed T_m ≈ 2.2 × 10⁵ K.

The BATS-R-US (Block Adaptive Tree Solar Wind Roe-type Upwind Scheme) code has been developed at the UM. It solves the equations of ideal MHD—a system of eight equations describing the transport of mass, momentum, energy, and magnetic flux (Groth et al. 2000; Powell et al. 1999). This massively parallel code enables Sun-to-Earth simulations to be performed in near real-time when run on hundreds of processors on a supercomputer (Manchester et al. 2004b). The implementation of adaptive mesh refinement (AMR) in BATS-R-US allows orders of magnitude variation in numerical resolution within the computational domain. This is important for a global model of the solar magnetic fields in which one strives to resolve such structures as shocks, volumetric currents of flux ropes, and electric current sheets in a 3D domain, which may extend out to hundreds of R_☉. The use of AMR also enables us to resolve the fine structure of active regions on the Sun, which spawn CMEs. This is vitally necessary for incorporating high-resolution magnetic data into a numerical MHD model. In the context of solar-heliospheric physics, BATS-R-US has been utilized to model the global structure of the solar corona and solar wind (Roussev et al. 2003b; Cohen et al. 2007, 2008), as well as the initiation (Roussev et al. 2003a; Jacobs et al. 2009) and propagation of idealized (Manchester et al. 2004a, 2004b) and not-so-idealized (Roussev et al. 2004, 2007) solar eruptions and associated solar energetic particle (SEP) events (Sokolov et al. 2004).

The SWMF is a high-performance computational tool that has been developed at the University of Michigan to simulate the coupled Sun–Earth system (Tóth et al. 2004, 2005). One of the main modules within the SWMF is the BATS-R-US MHD code. The SWMF is a structured collection of software building blocks to develop Sun–Earth system modeling components, in order to couple and assemble them into applications. The framework was designed to have "plug-and-play" capabilities, and currently it links together nine inter-operating models of physics domains, ranging from the surface of the Sun to the upper atmosphere of the Earth. Tying these models together gives a self-consistent whole in which each region is described by a world-class model, and those models communicate data with each other. The SWMF has achieved faster than real-time performance on massively parallel computers, such as NASA's Columbia supercluster.

To simulate the CR2107, a data-driven boundary condition is used for the inner boundary magnetic field. We use the SDO/HMI synoptic magnetogram in this study (see Figure 3). Given the uncertainty of the polar region field, a correction is made to the polar field (Sun et al. 2011). To extrapolate the initial potential field, Finite Difference Iterative Potential-field Solver (FDIPS) is used (Tóth et al. 2011). The resolution of the magnetogram is 3600×1440.

Zoom In Zoom Out Reset image size

A spherical grid is used in the simulation. The simulation domain inner boundary is at R = R_☉, and outer boundary at R = 20 R_☉. In total, 2.5×10⁵ blocks (6*4*4) are used with 2.1×10⁷ cells. AMR is performed to resolve the heliosphere current sheet (HCS). The smallest cell reaches ∼10⁻³ R_☉, while the largest is ∼0.8 R_☉. The radial velocity at X = 0 plane is shown in Figure 4. The fast and slow winds are evident in that plot. The fast wind speed is ∼750 km s⁻¹. The slow wind speed is ∼300 km s⁻¹.

Zoom In Zoom Out Reset image size

In Figures 5 and 6, zoom-in velocity and plasma beta with magnetic field lines are shown near the Sun. The velocity pattern shows more complex structures compared to the solar minimum condition. The magnetic field shows many structures (e.g., streamers, pseudo-streamers). From the plot of the plasma beta, we can see that the different regions of the corona have very different plasma beta. The polar regions with open field lines have small beta ∼0.01 near the Sun and increase outward. The high beta regions show the location of the current sheets in which the magnetic field is very small.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to evaluate the new steady-state solar wind model near the Sun, we compare the model output to the electron temperature and density derived from EUV images of the Sun by using the DEMT method (Vásquez et al. 2010). In general, the DEMT method uses a time series of EUV images under the assumption of no time variation and uniform solar rotation to derive 3D emissivity distribution in each EUV band. By Local Differential Emission Measure (LDEM) analysis, the 3D distribution of the electron density and temperature can be obtained. The DEMT method assumes that the plasma is optically thin. In this study, we use three bands of EUV observation (171, 193, and 284 A) from SDO/AIA.

In Figures 7 and 8, we show the ratio between the model and DEMT output for the electron density and temperature. The sphere is at r = 1.05 R_☉. The ring is between r = 1.035 R_☉ and r = 1.225 R_☉. For both the density and temperature, we see a good agreement between the model and the DEMT output for most regions. The active regions and polar regions in the model have relatively smaller density (note that the DEMT method has larger uncertainty for the active and polar regions).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 9 we provide a direct comparison of the observed EUV images with those synthesized from simulations. For better visibility we marked the active regions and coronal holes. One can see that our numerical simulations well reproduce the observed morphological structures, thus confirming the physical reasonings of the new global coronal model.

Zoom In Zoom Out Reset image size