Threaded-Field-Lines Model for the Low Solar Corona Powered by the Alfven Wave Turbulence

We present an updated global model of the solar corona, including the transition region. We simulate the realistic tree-dimensional (3D) magnetic field using the data from the photospheric magnetic field measurements and assume the magnetohydrodynamic (MHD) Alfv\'en wave turbulence and its non-linear dissipation to be the only source for heating the coronal plasma and driving the solar wind. In closed field regions the dissipation efficiency in a balanced turbulence is enhanced. In the coronal holes we account for a reflection of the outward propagating waves, which is accompanied by generation of weaker counter-propagating waves. The non-linear cascade rate degrades in strongly imbalanced turbulence, thus resulting in colder coronal holes. The distinctive feature of the presented model is the description of the low corona as almost-steady-state low-beta plasma motion and heat flux transfer along the magnetic field lines. We trace the magnetic field lines through each grid point of the lower boundary of the global corona model, chosen at some heliocentric distance, $R=R_{b}\sim1.1\ R_\odot$ well above the transition region. One can readily solve the plasma parameters along the magnetic field line from 1D equations for the plasma motion and heat transport together with the Alfv\'en wave propagation, which adequately describe physics within the heliocentric distances range, $R_{\odot}<R<R_{b}$, in the low solar corona. By interfacing this threaded-field-lines model with the full MHD global corona model at $r=R_{b}$, we find the global solution and achieve a faster-than-real-time performance of the model on $\sim200$ cores.

1. INTRODUCTION CME initiation by magnetic free energy, these simulations are often performed in a small Cartesian box (e.g. Torok and Kliem 2005), or using global models with no solar wind (e.g. Antiochos et al. (1999) and Fan and Gibson (2004)). So far there have only been a few magnetically driven Sun-to-Earth CME simulations through a realistic interplanetary medium using 3-D MHD (cf. Manchester et al. (2004aManchester et al. ( , 2004bManchester et al. ( , 2005, Lugaz et al. (2007) and Tóth et al. (2007)). The MHD simulation of Tóth et al. (2007) was able to match the CME arrival time to Earth within 1.8 hours and reproduce the magnetic field magnitude of the event.
A simple but convenient way to simulate a magnetically-driven CME is to superimpose a Gibson and Low (1998) (GL) or Titov and Démoulin (1999) (TD) magnetic flux-tube configuration onto the background state of the SC. Specific examples of such CME simulations using the AWSoM model for the SC and IH with a superimposed GL magnetic configuration include Manchester et al. (2012) and Jin et al. (2013Jin et al. ( , 2017a. The GL magnetic configuration describes an erupting magnetic filament filled with excessive plasma density. That filament becomes an expanding flux rope (magnetic cloud) in the ambient solar wind while evolving and propagating outward from the Sun, thus allowing the simulation of the propagation to 1 AU of a magnetically driven CME. In a similar way, by superimposing multiple TD configurations, Linker et al. (2016) have recently modeled the July 2000 CME eruption.
Here we present the development of the AWSOM. A distinctive feature of the presented model is the description of the low SC as almost-steady-state low-beta plasma motion along the magnetic field lines, the heat fluxes also being aligned with the magnetic field. The Low Solar Corona model which ranges from the upper chromosphere to the heliocentric distances about ∼ 1.1 R and includes the transition region at R < R < 1.03 R , is the heart of the global models. In the low SC the Alfvén waves pass from the chromosphere to the solar corona, the plasma temperature increases by two orders of magnitude (from ten thousand to million K), and this is also a place where the solar wind originates. The multi-wavelength observations (in EUV and X-rays) from several satellite locations (SDO, STEREO A,B) may be used to validate the simulation model. Therefore, any global model must account for the processes in this region. On the other hand, for the simulation model to explain the space weather and also have a predictive capability, it should be capable of simulating the dynamic processes faster than they proceed in real time, and the low SC appears to be a bottleneck limiting the computational efficiency and performance.
In numerical simulations of the solar corona, both for the ambient state and especially for dynamical processes, the greatest number of computational resources is spent for maintaining the numerical solution in the low SC and in the transition region, where the temperature gradients are sharp and the magnetic field topology is complicated. The degraded computational efficiency is caused by the need for the highest resolution as well as the use of a fully three-dimensional implicit solver for electron heat conduction. The need to find a numerical method, which would allow us to gain in the computational efficiency, motivates the research presented here.
We benefit from the observation that although the simulations of the low SC are computationally intense, the physical nature of the processes involved is rather simple as long as the heat fluxes and slow plasma motional velocities are mostly aligned with the magnetic field. The Alfvén wave turbulence, is characterized by the wave Poynting flux, which is also aligned with the magnetic field. Therefore, the plasma state at any point within the low SC is controlled by the plasma, particle, and Alfvén wave transport along the magnetic field line, which passes through this point. This physical property is typical for a variety of magnetized plasmas in different astrophysical and laboratory environments and may be used as the base of a new numerical method, which solves the state of plasma in each grid point in the computational domain depth in the following way: (1) by passing the magnetic field line ('thread') through this point and connecting it with the domain boundaries (e.g., with chromosphere and with the global solar corona domain, once the method is applied to the low SC) and (2) by solving a set of one-dimensional transport equations to relate the grid point value to the boundary conditions.
We trace the magnetic field lines through all grid points of the lower boundary of the global coronal model chosen at some heliocentric distance R = R b ∼ 1.1 R well above the transition region. One can readily solve the plasma parameters along the magnetic field line from effectively 1D equations for the plasma motion and heat transfer together with the Alfvén wave propagation, which adequately describe physics within the heliocentric distance range, R < R < R b , i.e. in the low solar corona. By interfacing this Threaded-Field-Line Model (TFLM) for the low corona with full MHD global corona model at R = R b we find the global solution and achieve faster-than-realtime performance of the model with moderate computational resources. Due to the latter feature we call the newly developed model AWSoM-R (AWSoM-Realtime).

Full 3D Governing Equations of the Global Model
The global model within the range of heliocentric distances, R b < R < 1 ÷ 3 AU, R b ∼ 1.1 R employs the standard MHD equations (non-specified denotations are as usually): (herewith, B = |B|), with the full energy equations applied separately to ions ∂ ∂t and to electrons: where, for a hydrogen plasma, N e = N i = ρ/m p , m p being the proton mass. In addition to the standard effects, the above equations account for the radiation energy loss from an optically thin plasma, Q rad = N e N i Λ R (T e ), a possible difference in the electron and ion temperatures, T e,i , the electron heat conduction parallel to the magnetic field lines equals: where m e and e are the electron mass and charge correspondingly, Λ C being the Coulomb logarithm. The energy exchange between electron and ions is parameterized via the energy exchange rate, 3mp(2πk B Te) 3/2 , as this is usually done. The Alfvén wave turbulence pressure, P A = (w − +w + )/2, and dissipation rate, Γ − w − + Γ + w + , are applied in the above equations. Herewith, w ± are the energy densities for the turbulent waves propagating along the magnetic field vector (w + ) or in the opposite direction (w − ). The turbulence energy dissipation (see Eqs.4-5) is split into electron and ion heating. For simplicity, we assume a constant fraction of energy, f p ≈ 0.6, dissipated into ions, the problems caused by this assumption are discussed below. At higher densities as in the low corona, we assume T e = T i and use the single-temperature energy equation for electron and ions, to improve the computational efficiency: where P = P e + P i = 2N i k B T . We use the equation of state, P e,i = N e,i k B T e,i for the coronal plasma with the polytropic index, γ = 5 3 . To complete the model, the equation describing propagation, reflection and dissipation of turbulent waves has been derived in van der Holst et al. (2014) following the approach as adopted in Velli (1993) Tu and Marsch (1995), Dmitruk et al. (2002),  and Chandran and Hollweg (2009)): ρ (note that the definition of L ⊥ and, accordingly the expression for Γ ± used both here and in van der Holst et al. (2014) are by a factor of 2 different from those used in Sokolov et al. (2013)). The reflection coefficient has been found as follows: where I max = 2 is the maximum "imbalance degree". If the "plus" wave strongly dominates, so that w + /w − > I max , the multiplier in the second line tends to +1, in the opposite limiting case of the dominant "minus" wave, it tends to -1. In both these cases the reflection reduces the dominant wave and amplifies the minor one. Otherwise, if the both amplitude ratios do not exceed I max , the turbulence is treated as "balanced" and the reflection coefficient turns to zero. The reflection model used by (Matthaeus et al. 1999) was similar to ours.
An important distinction is that we don't introduce the incompressible-to-compressible mode conversion term proportional to u · ∇ log V A into our model, although it is sometimes accounted for by other authors. The reason for this omission is that this term would break the energy conservation in the model, because it describes the conversion to the compressible MHD turbulence, which is not included (for more detail see van der Holst et al. (2014)).
The boundary condition for the Poynting flux at the top of chromosphere, S A is given by Herewith, we denote with braces the constant parameters of the model, equal to a product or ratio of physical variables. The estimate of the constant Pointing-flux-to-field ratio at the solar surface may be found in Sokolov et al. (2013), Oran et al. (2013) and van der Holst et al.
where the boundary condition for the wave energy density should be applied to the outgoing wave only. The estimate is very close to that which follows from Pevtsov et al. (2003), Suzuki (2006), Abbett (2007), Downs et al. (2010) and Cranmer (2010). To close the model we chose, following Hollweg (1986), the scaling law for the transverse correlation length: Eq.(8) has a form close to the conservation law, which is well suited for solving it numerically within the framework of the global coronal model. However, both for using in the TFLM model and for analytical solution, an alternative form of this equation may be derived based on the substitution, Using the mass conservation law, one obtains: The plasma heating function, Γ + w + + Γ − w − , in these variables equals The dimensionless amplitude, a ± , of the outgoing wave at the lower boundary of the model equals unity. In Eqs.(11) the dimensionless wave amplitudes depend on the plasma dynamical profile only via the plasma velocity as well as the Alfvén speed. In the inner heliosphere, the Alfvén speed is negligible compared to the solar wind speed. Assuming steady state radial solar wind motion with the constant speed (i.e. independent on the heliocentric distance), the dimensionless amplitudes, mass density and the total turbulence energy density decay with the heliocentric distance as follows: , the latter relationship being in agreement with the polytropic index of 3/2, for the Alfvén wave turbulence. In the low SC the Alfvén wave turbulence dynamics is more complicated and discussed below.

Equations of the Threaded-Field-Line Model
Now, the governing AWSoM equations may be applied to simulate the transition region and Low Solar Corona domain at R < R < R b . We present both the simplified 1D model equations for this domain and the way how the model may be interfaced both to the chromosphere at R = R and to the global corona model at R = R b .

Magnetic field
The realistic model for the 3D solar magnetic field includes the boundary condition for the coronal magnetic field taken from the full disc magnetogram incorporating the current and past observation results. The potential magnetic field provides the minimum of magnetic free energy for given boundary condition, therefore, in the "ambient" solution for the solar wind the magnetic field is approximately equal to the potential one in the close proximity of the Sun. Following Ogino and Walker, (1984) and Tanaka (1994) (see also Powell et al. (1999) and Gombosi et al, (2002)) we split the total magnetic field B = B 0 + B 1 , in such way that the potential B 0 field dominates at R = R . If the observable is the radial component of the magnetic field at the photospheric level, then the potential B field may be recovered from the observed magnetogram using the Potential Field Source Surface Method (PFSSM) had been for the first time described in Altschuler et al (1977). The Laplace equation for scalar magnetic potential is solved at R < R < R SS = 2.5R with the given radial gradient of the potential (the observed radial field) at R = R and with vanishing magnetic potential (i.e. purely radial magnetic field) at R = R SS , using the development into a series of spherical harmonics. Accordingly, non-potential B 1 field within the original split field approach (used, particularly, in Sokolov et al. (2013), Oran et al. (2013) and van der Holst et al. (2014)) satisfies zero boundary condition for the radial field component, (B 1 ) R = 0, at R = R , as long as the observed field is fully included into the potential field.
The distinction of the approach presented here is that we neglect non-potential B 1 field in the Low Corona and assume that B 1 ≡ 0 at R < R < R b . Accordingly, the boundary condition, (B 1 ) R = 0, is accepted within the global model at R = R b . In this way we benefit in easily bridging the field observed at R = R to the model starting at R = R b . Second, the lines of the potential, B 0 , field at R < R < R b give us the threads which allow us to bridge the boundary conditions for all other physical quantities from the top of chromosphere to the global model boundary at R = R b .

Magnetic thread and the conservation laws on it.
Now, we introduce a key concept of the Threaded-Field-Line Model (TFLM) -a thread. The boundary conditions for the global model are to be applied at each grid point of the global model boundary at R = R b . The potential magnetic field line, "thread", starting at the grid point can be traced through the Low Corona domain, R < R < R b toward the Sun. To reduce the 3D governing equations to effectively 1D equations, one can integrate Eqs.(1,3,7) over a magnetic flux tube element of a length of ds bounded by two close cross sections of the flux tube, dS 1 and dS 2 , and a bundle of magnetic field lines about the considered thread, which all pass through the contours of these cross sections. The equation, ∇ · B = 0 gives: BdS = const along the flux tube, which allows us to relate the change in the cross-section area along the thread to the magnetic field magnitude. The conservation laws are greatly simplified due to the fact that the velocity of low-beta plasma motion is aligned with the magnetic field. Particularly, the continuity equation (1) where u = (b · u) is the velocity aligned with the magnetic field and ∂ ∂s = (b · ∇). Herewith we denote with braces the combinations of variables, which are constant along the thread (above we did this only for the model parameters, {S A /B} and {L ⊥ √ B} ). As long as the velocity is not solved within the TFLM, the parameter in Eq.(13) should be found from the Global Corona Model (GCM) side: u B T F LM = lim R→R b +0 ((B · u)/B 2 ) GCM . In the momentum equation we neglect u 2 relative to the speed of sound squared, γP/ρ and omit j × B force, vanishing in the potential magnetic field (j ∝ ∇ × B 0 = 0), which gives us the hydrostatic equilibrium equation: The latter can be integrated, if desired, for the given profile of temperature giving the barometric formula, , the values of variables on top of the transition region (TR) are discussed below.
In Eq. (7) we keep the time derivative of temperature as long as the electron heat conduction is a comparatively slow process: Note, that we neglect the Alfvén wave pressure gradients in Eqs. (14,15). Practically in the Low Corona this pressure is small, being proportional to a square root of high density, while the thermal pressure is proportional to the density. Theoretically, keeping this term in Eq.(15) would be inconsistent. Comparing with the Alfvén wave energy deposition (see below) it involves the small ratio of the plasma speed to the Alfvén wave speed, and all such terms are neglected in deriving Eq.(16) below.

Alfvén wave 1D dynamics.
The physical property of the Alfvén waves to have the energy flux aligned with the magnetic field allows us to reduce 3D differential operators in the governing equations to the advective derivatives along the magnetic field lines. In addition, the 1D governing equations, which are obtained in this way, may be further simplified for the low corona environments. Indeed, in the low SC the plasma velocity in Eqs. (11) is negligible compared to the Alfvén wave speed. The steady-state solutions for a ± may be sought for, as long as the non-stationary perturbations propagate with the large Alfvén wave speed across the low corona and quickly converge to an equilibrium, so that Eqs(11) once divided by V A may be written as follows: This equation may be further simplified by substituting As long as the plasma speed is small relative to the Alfvén speed, the velocity curl in the expression for the reflection coefficient is negligible compared with the contribution from the Alfvén speed gradient, hence: ds dξ a max = max(a − , a + ). The formulation of the boundary-value-problem assumes that at the starting point of the magnetic field line, i.e. at minimal ξ = ξ − , the boundary value a +0 should be given for a + wave, propagating in the direction of increasing ξ: a + (ξ = ξ − ) = a +0 . For the oppositely propagating wave the boundary value, a −0 , should be given at the right end point ξ = ξ + , of the magnetic field line section, [ξ − , ξ + ]: a − (ξ = ξ + ) = a −0 . For the case of the closed magnetic field line, starting and ending at the solar surface, a +0 = a −0 = 1 because of our choice of the Boundary Condition (BC) for the Poynting flux. A few examples of the problem formulation for Eq.(16) are delegated to Section 3. In Eq.(15) one can express (Γ − w − + Γ + w + )/B = 2(a − + a + )a − a + dξ/ds{S A /B} = d(a 2 − − a 2 + )/ds{S A /B} and on dividing Eq.(15) by {S A /B}, it can be rewritten as follows: In application to the TFLM it is convenient to denote with "+" and "-" subscripts the waves propagating ourward and inward correspondingly and assume that the variable s along the thread equals zero at the solar surface and is positive at R < R < R b . These assumptions require to re-define the BCss at the interface R = R b between the TFLM and GC models as follows:

BCs for temperature and density
The temperature is governed by Eq. (18) which is of the second order with respect to the spatial coordinate. Hence, at the interface between TFLM and GCM both temperature and its gradient should be continuous, so that the BC for temperature within the TFLM may be taken from the GCM: T T F LM = T GCM at R = R b . Accordingly, once the TFLM equations have been solved with the given temperature at R = R b and the BC at R = R as discussed below, the gradient, ∂T ∂s T F LM , at R = R b is known and can be used to set the radial temperature gradient within the GCM. Assuming the radial component of the temperature gradient to be dominant, one has: ∂T ∂s GCM ≈ ∂T ∂R GCM . Hence, the equation, ∂T ∂R GCM = ∂T ∂s T F LM /|b R |, may be used to close the boundary value problem in the GCM by setting the heating flux through the interface between TFLM and GCM.
The density at the discussed interface is controlled by the direction of u. If u > 0, then All we need now in order to get well-posed problem for the listed 1D governing equations is to close the model with the boundary condition at "low boundary". A good opportunity is to set the boundary conditions on top of the Transition Region (TR), which may be matched with the chromosphere via an analytical model of the transition region. This model has been presented by Lionello et al. (2001) (see also Lionello et al. (2009) and Downs et al. (2010)).
To use this model, we choose at each magnetic thread a short section of the length of, L T R = R T R R ds ∼ 1 Mm to be a width of the TR along the magnetic field line. This is an important distinction from previous works, in which the temperature on top of the TR had been set, rather than the TR width. In the steady-state version of Eq.(18) we keep only the terms which dominate within the TR, i.e. at high density and abrupt temperature gradients: On multiplying Eq.(19) by κ 0 T 5/2 (∂T /∂s) and by integrating from the interface to chromosphere to a given point at a temperature, T , one can obtain: Here the product, {N i T } is assumed to be constant, therefore, it is separated from the integrand, since the temperature gradient scale within the TR is mush shorter than the barometric scale: T /(∂T /∂s) k B T /(GM m p /R 2 ). On the left hand side of the equation the half of the heat flux squared should be taken at the temperature, T , with the positive sign and at temperature T ch with negative sign. We can neglect the contribution from the latter term at the lower boundary of this transition region, if we postulate that the transition region is heated solely by the heat transfer from the corona and the lower boundary of the transition region is where the heat flux from the corona turns to zero. For given T ch ∼ (1 ÷ 2) · 10 4 K and L T R one can solve {N i T } in terms of the radiation loss integral: which also allows us to find the heat flux into the TR from the Low Corona: Here, U heat is a quantity with the dimension of speed, such that 2L T R /((γ − 1)U heat (T T R )) is a typical temperature relaxation time in the TR. We arrive at nonlinear boundary conditions on top of the TR, which for the known width L T R , of the TR and for a given temperature, T T R , on top of the TR allow us to find the heat flux and pressure there using Eqs. (21,22). These BC may be easily implemented if the temperature functions in the right hand side of Eqs.(21,22) as well as the function, Λ R (T ), are all tabulated using the CHIANTI database ) (see Fig. 1). Now, the TFLM is fully described and designed to be solved numerically. There is still a minor uncertainty in the way to distinguish the TR from the top of chromosphere, originating from the fact, that we do not include a consistent chromosphere model. Particularly, the TR solution at some point should be merged to the chromoshere solution with no jump in pressure, at the location, which depends both on {N I k B T } T R and on the pressure barometric distribution in the chromosphere. However, the uncertainty in this location, which also results in some uncertainty in L T R , is negligible, because the barometric scale in the chromosphere is small.  Figure 1. Preprocessed CHIANTI table for radiative cooling, which allows us to formulate the boundary condition at the solar surface, for the TFLM. For a known width, L T R , of the internal transition region and for the input temperature, T , the constant product, {N T k B } may found using the dashed curve. Then, the heat flux into the transition region, U heat {N T k B }, may be found using also the solid curve.

ANALYTICAL SOLUTIONS AND SCALING LAWS FOR THE ALFVÉN WAVE TURBULENCE IN THE LOW SOLAR CORONA
As has been demonstrated above, the main point of the developed approach is to achieve efficient and realistic modeling of the solar atmosphere. On the other hand, the analytical solutions, discussed in the present Section are over-simplified and are relevant only for analyzing the equations describing the Alfvén wave dynamics. Although these solutions are not directly usable for doing simulations, they may give hints on dependencies between the different model parameters. Of a particular interest is the question, which model parameters should be modified to achieve a better agreement with the observations of the solar wind parameters at 1 AU. Therefore, we provide here some solutions describing the wave turbulence in the different regions (closed vs open field lines, lower vs global SC etc).

Solution for Coronal Holes -Weak Reflection
Here, we consider an open magnetic field line, by assuming as we did above, that the wave of amplitude, a + , propagates outward the Sun. First, consider the case when the reflection due to a gradual change in the Alfvén speed magnitude (the latter is assumed to exponentially decay outward the Sun, V A = V A0 exp(−s/L V A )) is small compared with the characteristic dissipation rate: 1 2 ds dξ This assumption is valid at sufficiently high altitudes much above the transition region (where, to the contrary the abrupt density gradients cause large reflection as we discuss below). As the result of a weak reflection, the amplitude of the wave reflected back and propagating toward the Sun is small: a − a + . The governing equations in this limiting case read: For the exponentially decaying profile of the Alfvén speed, by introducing a constant small parameter: one can easily find the solution of Eqs.(23), tending to zero at infinity: We found that, for the exponential profile of the Alfvén speed, the small ratio of the amplitude of the incoming wave to that for the outgoing wave is constant. This observation, in principle, allows us to calculate the dissipation rate for the dominant wave without calculating the amplitude of the reflected wave, since Γ + = 2 The latter comment, although valid only for a particular case of exponentially decaying V A , allows us to link the current model to that described in Sokolov et al. (2013), where we also parameterized the turbulence dissipation within the coronal holes using small dimensionless C refl . The WKB approximation we used in that paper predicted no inward propagating waves originating from the open magnetic field lines (w − = 0). Therefore, we assumed therein a small but finite (due to reflection) amplitude of the inward propagating wave to be parameterized as w − = C 2 refl w + , so that: In Sokolov et al. (2013) we did not discuss the reflection mechanism, so that C refl was an arbitrary and uncertain free parameter. In the model developed here we calculate the reflection coefficient R for realistic distribution of the magnetic field and plasma parameters, to greatly reduce the model uncertainty, however, we see that the choice of C refl = const ≈ 0.01 ÷ 0.1 in Sokolov et al. (2013) was reasonable and might be derived analytically. In a general case of an arbitrary (not necessarily exponential) profile of the Alfvén wave speed, Eqs.(23) still can be solved at small (not necessarily constant) value of C refl . Indeed, the total of two Eqs.(23) is a linear and easy-to-integrate equation, which gives: (a + + a − ) ∝ √ V A , so that a + ≈ V A /V A0 as long as a − a + , with constant V A0 , for a given thread being a characteristic value of the Alfvén wave speed at low altitude. Then, in the second of Eqs.(23) the left hand side is quadratic in small C refl , as the small (∝ C refl ) derivative of a smaller amplitude, a − ∼ C refl a + . Therefore, the two linear in C refl terms on the right hand side should cancel each other, which requirement gives: a − ≈ − 1 2 d log V A dξ . In this way, we arrive at a simple and transparent estimate for the wave energy density within the coronal holes, which linearly scales with the magnetic field: (12), we can now derive the following expressions for the heating function: Note also that for the Alfvén speed profile gradually increasing in the outward direction, the solution for the dominant wave energy density is: Now, we arrive at two important conclusions. First, within the coronal holes both the distribution of the turbulence energy and the heating function do not depend on the dissipation length, L ⊥ , as long as the minor wave amplitude and the dissipation rate are fully controlled by the wave reflection. The latter, in turn, is fully controlled by the field and plasma profile. The model of a self-consistent plasma state within the coronal hole, which is controlled by the heating function dependent only on the plasma state itself, seems to be reasonable and physics-based.
However, the numerical model for the coronal hole at larger heliocentric distances, R ≈ 2 ÷ 10, appears to be vulnerable to quasi-periodic formation of a narrow dip in the plasma density, with sharply peaked ion temperature amounting to 10 7 K. The dip in density corresponds to a local maximum in the Alfvén wave speed, V A = B/ √ µ 0 ρ. By merging the above given solutions for firstincreasing-then-decreasing profile of the Alfvén wave speed and matching the two constant values V A0 to maintain the wave amplitude continuity at the maximum of V A , one can obtain an estimate for the heating function near the density minimum: Γ + w + ∝ √ ρ min /L V A . This estimate becomes large for a sharp density profile, i.e. at L V A → 0. The enhanced heating, which is proportional to ∝ √ ρ may result in the ion temperature growth, because the ion specific heat effect (∝ ρ) and the ion-electron energy exchange rate drop faster as ρ → 0. Finally, the instability onset loop is closed by the ion thermal pressure effect, which sweeps out the locally overheated plasma and further reduces the local plasma density, ρ min , and the plasma specific heat. Eq.(4) allows us to evaluate a threshold for the instability under consideration. The ion heating rate due to the Alfvén wave turbulence dissipation with an account of the above considerations equals ≈ f p Rw + ≈ f p V A w + /L V A . If within the travel time, L V A /C S , for the sound waves escaping with the sound speed, C S , from the vicinity of the density local minimum, the ion temperature increases significantly, i.e. (f p V A w + /L V A )(L V A /C S ) > P i /(γ − 1), the instability is possible. Therefore, the threshold condition becomes as follows: f p > C S P i / [(γ − 1)V A w + ]. Based on this criterion, instability is possible only in a strong turbulence, in which the turbulent energy density (flux) is comparable with the energy flux related to the thermal motion: Since the signatures of such instability (sharp peaks in ion temperature anti-corelated with the plasma density) are not observed in the fast solar wind as originates from the coronal holes, the instability that we noticed in numerous simulations should be considered unphysical and it needs to be suppressed within the framework of numerical method. To achieve this, we limit the fraction of the turbulent energy dissipation absorbed by ions, in the regions of strong turbulence. Specifically, if V A w + ≥ P i C S , then instead of an equation, f p ≈ 0.6, we use the limited expression: f p = min (0.6, C S P i / [(γ − 1)V A w + ]). There is some physical reasoning in favor of such limitation: if the ions are comparatively cold then the number of ions which may be in gyro-resonance with the Alfvén wave is small so the ions cannot be efficiently heated by waves. Including the physics-based model for partitioning energy between ions and electrons (see van der Holst et al. (2014) cited therein), the stable model for the solar wind heating in the coronal loops may be achieved with no artificial limitation.

Scaling Laws for the TFLM
The weakness of any model relying on the solar magnetogram is uncertainty of the solar magnetic field observations. We introduce the parameter of the model, {S A /B} ≈ 1.1 · 10 6 W/(m 2 T) with one digit after a decimal period -is this legitimate? How accurate is the observed magnitude of the solar magnetic field to be multiplied by this model parameter? Numerous observatories provide different values for the measured field. The region in the solar wind which is determined by the polar coronal holes may be large and any realistic model of the solar wind should account these holes, however, the solar magnetic field measured in these holes, by many reasons, may be unrealistically low.
To mitigate the effect of too low and, probably, underestimated magnetogram field, we apply some scaling factor B scale ≥ 1 in our simulations, so that the observed solar magnetic field multiplied by B scale is used as the boundary condition of the model: B TFLM | R=R = B scale · B observed | R=R . We note that the TFLM equations are not affected by this scaling if in accordance with increasing the magnetic field we also decrease the model parameters: because the magnetic field in the TFLM equations is present only in combinations, {S A /B}B and

SIMULATION RESULTS AND DISCUSSION
All simulations were performed with the Space Weather Modeling Framework (SWMF -see Tóth et al. (2004Tóth et al. ( , 2005Tóth et al. ( , 2007Tóth et al. ( , 2012. The SWMF included the models (components) to simulate the Solar Corona and Inner Heliosphere, both models accounting for the contributions from the Alfvén wave turbulence as we described in Sokolov et al. (2013), Oran et al. (2013) and van der Holst et al. (2014) (the AWSoM model). The most important distinction in the current simulations is that we apply the AWSoM model only to the GCM, while the transition region and lower corona are described using the TFLM.
In this way we benefit from saving the computational resources which otherwise should be spent to resolve the true structure of the transition region using a highly refined grid. We start from the observation that the gain in computational efficiency is achieved with no degrade in the accuracy and quality of the numerical results. In Fig. 3 we show a comparison of the numerical results obtained with the ASWoM model (see Sokolov et al. (2013), Oran et al. (2013) and van der Holst et al. (2014)) and presented in the left panel, with the new result obtained with the TFLM (the AWSoM-R model). The grid for the AWSoM run requires much finer grid cells to resolve the transition region. In addition, in time-dependent runs these finest grid cells control and severely reduce the time step making is as short as a few ms. At the same time, in the AWSoM-R result (the right panel) there is no noticeable difference from that obtained with AWSoM, except that it may be obtained much faster and with the time step about one second, which is equivalent to a gain by a factor of several hundreds in the computational efficiency. The capability to do simulations significantly faster than the real time on a couple 100 CPU cores is the main advantage of the AWSoM-R with TFLM.
We did a much longer simulation to test the AWSoM-R model for the inner heliosphere. Specifically, using a simple dipole magnetic field as the boundary condition for the radial magnetic field component, we simulated over 47 days of physical time with the AWSoM-R model. We sampled plasma parameters at the same distances where the Ulysses spacecraft had passed. Figure 4 shows the comparison of radial velocity distribution along different latitudes between the SWOOPS (Bame et al. 1992) plasma measurements (blue) and the simulation results (red). We selected observations of the time frame between the years 1990 and 1997 so that we cover the complete range of latitudes. The observations are during the declining phase of the solar cycle. The solar wind distribution shows the clear difference between slow and fast solar wind regions. Both these wind speeds as well as the transition latitudes are captured by the AWSoM-R model.
The evolution of the solar magnetic field is comparatively slow, therefore, it is usually believed that the steady-state solution with the time independent boundary condition should be a good approximation. However, the advantage of time-dependent solution as we advocate here is in the capability to describe the time-dependent processes that occur even with the steady-state boundary condition for the magnetic field. To demonstrate this capability we simulated 10 days of evolution in the Solar Corona with the steady-state dipole magnetic filed. Such simulation reveals a dynamic helmet streamer structure, which periodically produces plasmoids known as "streamer blobs" (Wang et al. 2000). Our investigation shows that these blobs form as a result of pressure imbalance mainly because of increased ion temperatures at the streamer top. (see Fig.5).
To find a typical time of the blob formation, we choose a radial line in the equatorial plane rotating with the Sun (as an example we used y-axis of the HGR coordinate system) and visualize the distribution of plasma beta along this line as a function of time (x-axis) and radial coordinate(yaxis) as shown in Fig.6. The plasma beta starts increasing at the heliocentric distance of ≈9R s implying the start of the disconnection event and the disconnected blob moves anti-Sunward. These intermittent detachment events occur with a periodicity of about 40 hours (six times within 10 days), which is in a good agreement with the observations. We, thus, show that in the self-consistent Alfvén Wave Turbulence based model the slow solar wind is intermittent even if the solar magnetic field is steady-state and perfectly symmetric.

CONCLUSION
The AWSoM-R model presented here extends the earlier developed AWSoM  and van der Holst et al. (2014)) with the TFLM description for the transition region and Low Solar Corona. It allows us to simulate the Solar-Earth environments on realistic 3D grids faster than real time and with no loss in the results quality.

ACKNOWLEDGMENT
The collaboration between the CCMC and University of Michigan is supported by the NSF SHINE grant 1257519 (PI Aleksandre Taktakishvili). The work performed at the University of Michigan was partially supported by National Science Foundation grants AGS-1322543 and PHY-1513379, NASA grant NNX13AG25G, the European Union's Horizon 2020 research and innovation program under grant agreement No 637302 PROGRESS. We would also like to acknowledge high-performance computing support from: (1) Yellowstone (ark:/85065/d7wd3xhc) provided by NCAR's Computational and Information Systems Laboratory, sponsored by the National Science Foundation, and (2) Pleiades operated by NASA's Advanced Supercomputing Division Figure 5. Variation of plasma beta and magnetic field lines along XY plane in HGR coordinates taken 6 hours 40 minutes apart. As plasma pressure increases, it stretches out the field lines (upper panel) and the reconnected field lines move anti-Sunward (lower panel)