CORONAL HEATING BY SURFACE ALFVÉN WAVE DAMPING: IMPLEMENTATION IN A GLOBAL MAGNETOHYDRODYNAMICS MODEL OF THE SOLAR WIND

The ability to accurately forecast space weather events at Earth is limited, in part, by the lack of a realistic background solar wind model in which to simulate the propagation of solar disturbances such as coronal mass ejections (CMEs). Global magnetohydrodynamics (MHD) models require sources of energy and/or momentum in order to match both coronal and interplanetary solar wind properties (see discussion in Usmanov & Goldstein 2003), but the physical origins of these sources are not well understood. Approaches to modeling the solar corona without specifying the physical processes responsible for the extra energy include adding an empirical heating function in the MHD energy equation (Suess et al. 1999; Mikić et al. 1999; Sittler & Guhathakurta 1999; Groth et al. 2000; Lionello et al. 2009; Downs et al. 2010; Airapetian et al. 2011), and defining the polytropic index to be less than the adiabatic value and sometimes spatially varying its value (Riley et al. 2001; Roussev et al. 2003; Cohen et al. 2007). These models successfully benchmark with near-Earth quiet-time data (for example, by the Advanced Composition Explorer, Wind, and Ulysses). However, the successes achieved for quiet-time solar wind do not extend well to transient cases. First, Evans et al. (2008) showed that the resulting Alfvén speed profile in the lower corona (R < 8 R_☉) is too shallow, allowing very slow disturbances to drive strong shocks, in contrast with observations of radio bursts and solar energetic particle events. Second, Loesch et al. (2011) found that the thermodynamics of a modeled CME is sensitive to the coronal heating treatment in the solar wind. Third, the variable polytropic index approach distorts the basic thermodynamics of the solar wind and influences shock compression ratios in a modeled CME.

Therefore, a more physics-based model is needed in order to simulate CMEs in the lower corona. Our knowledge of the properties of the solar corona and origins of the solar wind are largely due to remote sensing of the corona by instruments like the Ultraviolet Coronagraph Spectrometer (Kohl et al. 1995) and the Large Angle Spectroscopic Coronagraph (Brueckner et al. 1995) onboard the Solar and Heliospheric Observatory (SOHO/UVCS and SOHO/LASCO, respectively). The observations inspire and provide tests for theories of coronal heating and solar wind acceleration. Although the applicability of these theories (to be discussed more in the next paragraph) to the solar corona is strongly debated, it is likely that many different physical processes are contributing to the heating and acceleration problems. The structure of the corona, with open and closed field lines and fast and slow solar wind, suggests that different processes might be dominant in the different regions.

Magnetic reconnection in nanoflares (Parker 1988; Klimchuk 2006) could play a dominant role in the heating of the slow solar wind, thought to originate in streamer cusps and the boundaries of streamers and active regions (Edmondson 2011; Brooks & Warren 2011). Alfvén waves are also a leading candidate: the wave pressure gradient has been suggested as a force to explain the acceleration of the fast solar wind above 700 km s⁻¹, and wave dissipation could provide a source of heat (Parker 1965; Belcher 1971). This theory is supported by observations of Alfvénic signatures in the heliosphere (Belcher & Davis 1971). Properties of these waves inferred from observations suggest that there is sufficient energy in the waves to heat the corona (Tomczyk et al. 2007; De Pontieu et al. 2007; McIntosh et al. 2010). The frequency range of interest for the damping of Alfvén waves in the lower corona is 10⁻⁴–1 Hz: waves above 1 Hz are strongly damped in the chromosphere (Leake et al. 2005), and the dominating Alfvén waves measured in the heliosphere have frequencies near 10⁻⁴ Hz (Belcher & Davis 1971). There are several damping mechanisms relevant to these low-frequency waves in the corona—see review by Narain & Ulmschneider (1996). Generally, studies of the dissipation of Alfvén waves are specific to either open or closed magnetic field regions. The dissipation of counterstreaming Alfvén waves via MHD turbulence has been proposed as a dominant mechanism along closed loops (Rappazzo et al. 2007; van Ballegooijen et al. 2011). For open field regions, commonly discussed mechanisms include ion-cyclotron resonance (Cranmer et al. 1999; Isenberg et al. 2001; Chandran et al. 2010) and turbulent dissipation of counterstreaming waves (Verdini & Velli 2007; Cranmer et al. 2007).

MHD models which include the interaction of the background with Alfvén wave energy often formulate wave dissipation with an empirical length scale (Allen et al. 1998; Ofman & Davila 1998; Usmanov et al. 2000); for discussion of the challenges facing wave-driven models, see Ofman (2010) and Cranmer (2010). Global models which specify physically motivated wave dissipation are currently being developed (Lionello et al. 2009, 2010; Oran et al. 2010; van der Holst et al. 2010). However, none of these models focus on heating near the boundaries of open and closed fields, even though the heating signatures and region characteristics used to constrain heating mechanisms are different from those in fast solar wind.

For example, the preferential heating of heavy ions observed by SOHO/UVCS in the center of coronal holes suggests that ion-cyclotron resonance is important, because its heating is proportional to the ratio A/Z (A is the atomic mass and Z is electric charge; Cranmer et al. 1999). However, SOHO/UVCS data showed that the temperature of Si xii was larger than the temperature of O vi at the streamer center and at the border with open field lines (Zangrilli et al. 1999). There is also some conflict with turbulent dissipation, another mechanism proposed to be important in the fast solar wind. This mechanism requires incoming (Sunward) waves, thought to be generated by the reflection of outgoing waves due to gradients in the Alfvén speed along the direction of the magnetic field (Velli 1994; Cranmer & van Ballegooijen 2005). However, at boundary regions, the gradient in the Alfvén speed transverse to the magnetic field could be comparable or stronger to that along the magnetic field (Vásquez et al. 2010).

Additionally, there are observations which suggest that strong heating occurs near the boundaries of open and closed fields. A technique called differential emission measure tomography (DEMT), which produces the local difference emission measure from extreme-ultraviolet images, was recently developed Frazin et al. (2009). For the solar minimum periods corresponding to Carrington rotations (CRs) 2068 and 2077, Vásquez et al. (2010, 2011) applied the DEMT method to Solar Terrestrial Relations Observatory Extreme UltraViolet Imager (STEREO/EUVI) and produced three-dimensional (3D) reconstructions of the electron density and temperature in the low corona (R < 1.225 R_☉). They found temperature enhancements within the streamer belt closed field region. By comparing the DEMT maps with global potential magnetic field models of those periods, they found that the temperature-enhanced regions were located along and around polarity inversion lines that are inside the streamer belt closed regions, near the boundary with the surrounding open regions. The largest temperatures are found inside the closed region, but enhanced temperatures are also present in the open surrounding regions. This result suggests that an efficient heating mechanism is operating near the interface of open and closed field lines and near polarity inversion lines. They also found a density enhancement in those regions with respect to the streamer core equatorial latitudes. In those regions, the density falls off more slowly with height than in the streamer equator. As a result, for the height range 1.1–1.2 R_☉, the largest streamer densities were located at intermediate latitudes between the equator and the open and closed magnetic boundaries. Because local electron densities are determined by the integral of the flux of dissipated Alfvén waves and the radiative heat flux along a field line, the observed electron enhancement could indicate a downward heat flux due to wave damping above. Lastly, Vásquez et al. (2011) found that small closed loops in which the temperature increased from the footprint to apex were located near the boundaries of open and closed fields, also corresponding to the temperature-enhanced regions mentioned above, which indicates that the heating in this region does not distinguish between open and closed field line geometries.

Combining all of these points suggests that some other process (aside from ion cyclotron and turbulent dissipation) could be dominant near the boundary regions of open and closed magnetic fields and along polarity inversion lines. We suggest that wave dissipation in this region is driven by transverse gradients in the local Alfvén speed. One such mechanism is phase mixing (Heyvaerts & Priest 1983; Ruderman et al. 1999; Ofman & Aschwanden 2002). Magnetic field lines, initially at rest in a plasma that is nonuniform in the transverse direction, are perturbed. These oscillations travel outward along the magnetic field lines with different phase speeds, which leads to frictional damping. It has been argued that phase mixing is not likely to be strong along closed field loops or in the region of heating in coronal holes (Parker 1991; Ruderman et al. 1998; De Moortel et al. 1999; Aschwanden 2004; Erdélyi 2005). For these reasons, we do not consider it as the most likely mechanism to explain the previously discussed observations.

Another wave dissipation that could be dominant at the boundary of the fast/slow solar wind is resonant absorption (Ionson 1978; Hollweg 1987, 1997; Parker 1991; Goossens et al. 2002, 2011). In a nonuniform plasma, such as a discontinuity separating two uniform but different plasmas, surface waves can be supported. If instead a smoothly varying layer separates the two plasmas, then the propagating wave will resonate with the layer at the location where the wave phase speed is the same as the local Alfvén speed (Hasegawa & Uberoi 1982; Lee & Roberts 1986). The gradients are steepened until dissipation of the wave energy occurs. The physical dissipation mechanism must involve viscosity and/or resistivity, but the method of energy conversion (from wave energy to heating of the plasma) does not need to be specified in order to develop analytic expressions of the dissipation (Goossens & Ruderman 1995).

Resonant absorption has been applied to the problem of heating small-scale loops (Ionson 1978; Poedts et al. 1990; Sakai et al. 2001). Aschwanden et al. (2003) compared resonant absorption rates that were measured with Transition Region and Coronal Explorer observations to theoretical approximations, and found general agreement (a ratio of the damping time to oscillation period of 2.0 and 1.8, respectively). In Verth et al. (2010), Coronal Multi-channel Polarimeter (CoMP) data (Tomczyk & McIntosh 2009) in the 1–4 mHz range along closed loops was used to estimate resonant absorption due to gradients in density across the loops. They found that resonant absorption is an efficient dissipation process in this region and frequency range.

Resonant absorption of surface Alfvén waves (SAWs) has also been proposed as a heating mechanism along open flux tubes in the corona. During solar minimum, the equatorial region contains large-scale streamers (closed magnetic loops). The edges of coronal hole boundaries must expand much more rapidly than radially (so-called superradially) in order to fill the space above the streamers. This superradial expansion of field lines at heights of 1–2 solar radii from the surface has been verified observationally with SOHO/UVCS observations (Dobrzycka et al. 1999; Strachan et al. 2012). Strachan et al. (2012) confirm that, for solar minimum conditions, the largest flux tube expansion occurs near the interface of coronal holes and streamers (i.e., the open/closed magnetic field boundaries). Previous modeling efforts of one- and two-dimensional (1D/2D) open flux tubes with prescribed magnetic field geometry and solar wind parameters demonstrated that SAW damping could be a significant contribution to heating in these regions (Jatenco-Pereira & Opher 1989; Narain & Sharma 1998).

The first study of SAW damping along open field lines using a three dimensional (3D) solar wind model was performed by Evans et al. (2009). They estimated how much heat could be deposited along open field lines due to SAW dissipation in a given solar wind steady state solution using the 3D MHD model of Cohen et al. (2007). They found that waves with frequencies larger than 0.28 mHz would be appreciably damped in the low corona, and would provide a total energy contribution to the solar wind that was comparable to the ad hoc heating term in the Cohen et al. MHD model. All of these theoretical, numerical, and observational studies support the notion that wave dissipation could replace ad hoc heating in global MHD models of the solar corona.

Evans et al. (2009) estimated what the contribution of SAW would be to coronal heating given a self-consistent MHD steady state solution for the solar corona and wind that did not include waves, i.e., the results were not derived self-consistently because the model did not treat waves. Thus, the interaction between the heating mechanism and the plasma could not be studied. In order to better understand the role of SAW dissipation in the solar corona and wind evolution, it is necessary to determine if it remains an important effect when it interacts self-consistently with the plasma; doing so is the objective of this study.

The numerical tool we use is a 3D data-driven solar wind model within the Space Weather Modeling Framework (SWMF; van der Holst et al. 2010; Toth et al. 2011). The wave pressure gradient accelerates the plasma, and coronal heating is achieved through wave dissipation. In van der Holst et al. (2010), the wave dissipation was set in a manner set to mimic turbulence (Hollweg 1986). In this work, we implement SAW damping into the van der Holst et al. model, and show that the addition of SAW dissipation improves agreement with observations from the low corona to the inner heliosphere. We also discuss implications for shock formation.

This paper is organized as follows: Section 2 presents the theory and simulation details, Section 3 provides results and connections to observations which support the presence of SAW damping, and Section 4 contains discussion and conclusions.

2.1. Wave-driven Solar Wind

In the low-frequency, Wenztel–Kramers–Brillouin (WKB) limit, MHD waves are treated as wave train packets (Dewar 1970; Jacques 1977; Sokolov et al. 2009), and therefore their evolution can be described by a wave transport equation written as

$$\begin{equation} \frac{dE_W}{dt} + \nabla \cdot {\bf F} = {\bf u} \cdot \left(\nabla \cdot \bf {P_W} \right) - \Gamma E_W, \end{equation} \tag{ 1 }$$

where F is the wave energy flux,

$$\begin{equation} {\bf F}= {{\bf V{}}_g} E_W + {\bf u} \cdot {{\bf P{}}_W}, \end{equation} \tag{ 2 }$$

and E_W is the wave energy density, u is the bulk fluid velocity, P_W is the wave pressure tensor, and Γ is the wave energy dissipation rate. The group velocity V_g is defined as

$$\begin{equation} {{\bf V{}}_g} = \frac{ \partial \omega }{\partial {\bf k}} = \frac{ \partial \omega _0}{\partial {\bf k}} + {\bf u} = {{\bf V}_{g0}} + {\bf u}, \end{equation} \tag{ 3 }$$

where ω is the frequency in the moving frame (subscript "0" indicates inertial frame) and k is the wave vector. In the case of a transverse Alfvén wave, V_g0 is equal to the Alfvén velocity (v_A = ${\textbf {B}}/{\sqrt{4\pi \rho }}$, where B is the local magnetic field and ρ is the plasma mass density), and E_W = (δB)²/8π.

The wave pressure tensor is isotropic for Alfvén waves in the low frequency limit and can be written as p_W = (1/2)E_W. A spectrum of waves (E_W = ∫dωE'_W(ω)) will evolve according to (Sokolov et al. 2009)

$$\begin{eqnarray} \frac{\partial E^{\prime }_W (\omega)}{\partial t} &+& \nabla \cdot ({{\bf V{}}_g} E^{\prime }_W (\omega)) - \frac{1}{2} (\nabla \cdot {\bf u})\nonumber\\ &\times &\frac{\partial E^{\prime }_W (\omega)}{\partial {\rm log} (\omega)} = - \Gamma (\omega) E^{\prime }_W (\omega). \end{eqnarray} \tag{ 4 }$$

In this model, we consider a flat (or gray) spectrum of exclusively outgoing (anti-sunward) waves launched from the inner boundary. (Boundary conditions will be described in detail in Section 2.3.) This approach is equivalent to performing a frequency-integrated average to Equation (4). The resulting wave transport equation is (van der Holst et al. 2010)

$$\begin{equation} \frac{\partial E_W^{\pm }}{\partial t} + \nabla \cdot \left[ {{\bf V{}}_g} E_W^{\pm } \right] - \frac{1}{2}E_W^{\pm } \left(\nabla \cdot {\bf u} \right) = -\Gamma _{\pm } E_W^{\pm }. \end{equation} \tag{ 5 }$$

In the numerical implementation, Equation (5) is solved for two Alfvén wave energy densities, E⁺_W and E⁻_W. E⁺_W represents waves that propagate in the direction of B (and therefore having a velocity +v_A). E⁻_W represents waves that propagate in the direction antiparallel to B (and therefore having a velocity −v_A). In both cases, all of the waves propagate away from the Sun.

There is no wave reflection in the WKB limit; therefore there are no counterpropagating waves in the model. Additionally, the wave intensities E⁺_W and E⁻_W are averaged over one wave period and not a frequency-resolved spectrum, so we cannot treat phase evolution. Therefore, wave–wave interactions (like mode conversion or a turbulent cascade) cannot be modeled self-consistently; doing so would require a treatment with Elsasser variables, as in Cranmer et al. (2007), Matthaeus et al. (2003), and Verdini et al. (2010). We emphasize that the wave dissipation mechanisms (to be discussed in Section 2.2) are self-consistently implemented but phenomenological by construction.

The waves are coupled to the plasma via the MHD energy equation,

$$\begin{equation} \frac{\partial E_b}{\partial t} + \nabla \cdot {E_b \bf u} + (\gamma -1)E_b\nabla \cdot {\bf u} = \Gamma _+ E_W^+ + \Gamma _- E_W^- \end{equation} \tag{ 6 }$$

and the momentum equation,

$$\begin{equation} \rho \frac{d{\bf u}}{dt} = -\nabla p -\nabla \frac{B^2}{8 \pi } + \frac{\left({\bf B} \cdot \nabla \right){\bf B}}{4 \pi } + \rho {\bf g} -\nabla p_W^+ -\nabla p_W^-, \end{equation} \tag{ 7 }$$

where E_b is the internal energy density of the background plasma and γ is the polytropic index (set to the adiabatic value of 5/3). The equations of mass continuity, magnetic induction (for ideal MHD), and the divergenceless of B close the system. The reader is referred to van der Holst et al. (2010) for more details.

For this work, the crucial term is the one which describes the wave energy density dissipation rate, Γ. It can be expressed generally as

$$\begin{equation} \Gamma = \frac{V_g}{L}, \end{equation} \tag{ 8 }$$

where L is the length scale associated with the dissipation mechanism. In the following section, we describe how we calculate Γ for each dissipation mechanism.

2.2. Wave Dissipation Mechanisms

Here, we describe the physics of the wave dissipation utilized in this study. We emphasize that the dissipation is implemented self-consistently, but the forms of the dissipation rates are ad hoc by construction.

2.2.1. Surface Alfvén Wave (SAW) Damping

The environment of the solar corona is highly structured, consisting of regions in which local plasma conditions can strongly vary. The presence of variations in the local plasma environment, such as the Alfvén speed, across structures like open flux tubes can affect the propagation of waves traveling along these field lines. For example, the interface between flux tubes can host surface waves (Ionson 1978; Wentzel 1979). A schematic of this interface is shown in Figure 1: two uniform media, separated by a region in which the local Alfvén frequency varies in the x-direction from ω₁ to ω₂. The evolution of a surface wave with frequency ω₀ (between ω₁ and ω₂), which propagates in this inhomogeneous layer, has been derived in the case where the plasma varies linearly in the interface as an initial boundary value problem (Lee & Roberts 1986). Starting from the ideal MHD equations and prescribing an exponential function to describe the waves, it can be shown that the Fourier-transformed equation for the wave evolution contains a critical point when the wave frequency equals the local Alfvén frequency. The damping can be categorized as strong if the layer width is small compared with the wavelength, i.e., k₀a ≪ 1 (k₀ is the wavenumber), and in this limit the damping rate for SAWs can be expressed as in Ionson (1978), Hasegawa & Uberoi (1982), and Lee & Roberts (1986):

$$\begin{equation} \Gamma _{{\rm SW}} = \pi (\bar{k}a) \left(\frac{\omega _{2}^2-\omega _{1}^2}{8 \omega _0}\right), \end{equation} \tag{ 9 }$$

where $\bar{k}$ is the average wavenumber. The subscript "SW" stands for surface wave.

This decay, in the context of ideal MHD, is not strictly a heating rate, as there is no dissipative process available to transfer the wave energy. Instead, this coupling between the SAW and the bulk Alfvén waves can be considered a conversion process of the SAW to other wave modes. These wave modes can then undergo dissipative processes (viscosity, resistivity, and other kinetic effects), and this ultimately leads to plasma heating. In the ideal MHD model presented in this paper, we do not model any dissipative processes, but instead assume that all wave dissipation leads directly to plasma heating. We propose the following form for the damping rate and length to be implemented in a global MHD model. We consider only gradients in density and write Equation (9) as

$$\begin{equation} \Gamma _{{\rm SW},i} = \frac{V_{g,i}}{L_{{\rm SW},i} } \end{equation} \tag{ 10 }$$

$$\begin{equation} L_{{\rm SW},i} = \frac{C_{{\rm SW}}^{\prime }}{\omega } \frac{\rho }{\sqrt{ \left(\frac{\partial \rho }{\partial x_j}\right)^2 + \left(\frac{\partial \rho }{\partial x_k}\right)^2 } }, \end{equation} \tag{ 11 }$$

where C'_SW controls the damping strength, and i, j, and k represent the X-, Y-, and Z-directions. It is assumed that the wave polarization is in the direction perpendicular to the resonant layer. This formalism captures the behavior of the mechanism: linear dependence on frequency and characterized by transverse gradients. It should be emphasized that this phenomenological damping rate is in agreement with those formally derived in the thin tube (strong gradient) limit (Terradas et al. 2010; Soler et al. 2011).

The density gradients calculated in the simulation were compared to values calculated from the DEMT reconstruction of the mean electron density (Vásquez et al. 2010) from R = 1.035 to 1.225 R_☉ for the same time period. Although the simulation does not resolve small-scale features present in the DEMT reconstructions, the order of magnitude and distribution of the gradients in the simulation are in agreement with those calculated using the tomographic results.

In this paper, Alfvén waves are modeled as a flat (or gray) spectrum, so the frequency dependence in Equation (11) is absorbed into the free parameter (C_SW = C'_SW/ω). In reality, C'_SW is not a free parameter: it depends on the background magnetic field intensity and geometry (through the expansion factor of the flux tube), and the density at the base of the corona, and it is also a function of frequency. All of these considerations were taken into account in Evans et al. (2009), where the analysis was performed only along selected magnetic field lines. However, in the implementation in the 3D model with a flat/gray spectrum, the frequency dependence is absorbed, and so it is truly a scaling factor. It is beyond the scope of this work to constrain the free parameter by the amplitude of the waves, background field intensity, or other physical parameters. In the future, when we can model a frequency-resolved spectrum, C'_SW can again be constrained with the local plasma parameters.

The damping rate Γ_SW is calculated at each spatial location for all waves, and at each iteration the waves are dissipated. The dissipated wave energy is passed to the plasma (as in Equation (6)), so that the total energy (waves and plasma) is conserved.

Although we do not treat a spectrum of waves, we can estimate the frequency range for which this treatment is valid. The condition of strong damping (ka < 0.1) is estimated by

$$\begin{equation} ka \approx \frac{\omega }{v_a} \frac{\rho }{\frac{d\rho }{dx}}, \end{equation} \tag{ 12 }$$

and is satisfied in these simulations for ω ⩽ 13 mHz. This estimated frequency limit overlaps both the range that Evans et al. (2009) found to be important for SAW damping (ω ⩾ 0.28 mHz), and the range that Verth et al. (2010) used to make observational constraints of resonant absorption along closed loops (ω = 1–4 mHz).

2.2.2. Hollweg Dissipation

Solar wind spectra display power laws, a known property of turbulent media (Coleman 1968). Turbulent cascades have been suggested as a method for transferring energy from large to small scales, where dissipation and plasma heating can occur. This process is rooted in the interaction of two counterstreaming waves. As discussed in the Introduction, it is assumed that the sunward waves are generated by reflection due to a sharp gradient in the Alfvén speed profile along the magnetic field (Dmitruk et al. 2002; Cranmer & van Ballegooijen 2005; Chandran & Hollweg 2009).

In this model, we are working in the WKB approximation where wave reflection is not possible, and we consider only the frequency-averaged wave energy density. Therefore, turbulent cascades cannot be treated explicitly. However, it is still possible to apply a turbulence-like dissipation as developed by Hollweg (1986) for the purpose of heating the corona. The motivation for Hollweg's equation was empirically based. The Kolmogorov cascade rate derived for isotropic turbulence was used to write a phenomenological description of the dissipation of Alfvén waves in the solar wind. The volumetric heating rate is expressed as

$$\begin{equation} \Gamma _{{\rm HOL}\pm } = \frac{1}{C_{{\rm HOL}}} \sqrt{\frac{E_W^{\pm }B}{\rho }}, \end{equation} \tag{ 13 }$$

where $L=C_{{\rm HOL}}\sqrt{(B)}$ is a correlation length of the Alfvén wave fluctuations, representative of the local plasma properties, and C_HOL is a free parameter to adjust the strength of the dissipation. As this formulation does not account for the loss of outgoing waves due to reflection, in general applying this term to a specific wave spectrum will result in an overestimation of the energy dissipation compared to a more strict calculation (as in Cranmer & van Ballegooijen (2005). The formalism in Equation (14) was used as the source for coronal heating in van der Holst et al. (2010). It resulted in strong heating in the center of the coronal holes (see Figure 6 of that paper).

The purpose of including Hollweg dissipation in this study is twofold: first, we compared the damping rate to determine the regions where SAW damping is important, and second, to demonstrate that multiple wave dissipation mechanisms are required in a global MHD model to achieve the best match with known properties of the solar wind.

2.3. Simulation Details and Data-driven Boundary Conditions

We obtain a steady state solar wind solution using the wave-driven solar corona (SC) component of the SWMF (van der Holst et al. 2010). The inner boundary condition for the radial component of the magnetic field is specified from the National Solar Observatory's Global Oscillation Network Group (GONG) magnetogram data for CR2077 (UT 2008 November 20 through December 17), as shown in Figure 2(a). The magnetogram data were multiplied by a factor 1.8 due to the low spatial resolution of GONG. The initial conditions for the magnetic field are calculated with the Potential Field Source Surface (PFSS) model (Altschuler & Newkirk 1969). The inner boundary conditions for the solar wind density and temperature (shown in Figures 2(b) and (c)) are determined from the DEMT reconstructions of mean electron temperature and density (Vásquez et al. 2010) using STEREO/EUVI data for CR2077. The initial conditions for the single-fluid plasma density and solar wind speed are given by the isothermal Parker solar wind solution (Parker 1958).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The wave energy density at the inner boundary is set to zero for all closed magnetic field lines, as determined by the PFSS model and the magnetogram data. The formalism we apply for the evolution of Alfvén waves considers the Alfvén wave energy averaged over a wave period in the WKB approximation. Therefore, there is no wave reflection in the model. Additionally, this formalism does not treat counterpropagating waves. In closed loops, counterstreaming waves occur naturally (from waves launched along both footpoints, wave bouncing, and reflection), and the interactions of these waves are important to describing the heating of small loops and streamers (Rappazzo et al. 2007). However, in this model, wave energy may leak onto closed field lines due to numerical diffusion, or the evolution of the field during the MHD steady state relaxation. In regions where this occurs, the dissipation is calculated according to Equations (13) and (14) (e.g., we have not defined limitations on the region where damping occurs).

For open field lines, a purely anti-Sunward (outgoing) wave energy density is specified at each boundary location by assuming conservation of energy along a flux tube and applying the Bernoulli integral from the surface to 1 AU (Suzuki & Inutsuka 2006). At 1 AU, the solar wind energy is assumed to be purely kinetic, with a velocity u_WSA calculated using the Wang–Sheeley–Arge (WSA) model (Arge et al. 2003). This approach allows the wave energy input at the boundary to be calculated according to

$$\begin{equation} E_W^{\pm }= \frac{\left(\rho u r^2\right)_{1 \,{\rm AU}}}{r^2} \frac{f_{{\rm exp}}}{|v_{A,r}|} \left(\frac{u_{{\rm WSA}}^2}{2} + \frac{GM_{\odot }}{R_{_\odot }} -\frac{\gamma }{\gamma -1}\frac{p}{\rho } \right), \end{equation} \tag{ 14 }$$

where f_exp is the expansion factor for the magnetic field line and ρu is the mass flux at 1 AU. The predicted wave pressure at the inner boundary is shown in Figure 2(d), and ranges from 0 to 0.0015 dyn cm⁻².

The goal of this paper is to determine in which region SAW damping is important in the solar corona. To achieve this goal, we explore two damping mechanisms, both individually, and in combination. Damping rates for each mechanism can be written as

$$\begin{equation} \Gamma _{{\rm SW}\pm } = \frac{V_{g,\pm }}{ C_{{\rm SW}} \frac{\rho }{\nabla _{\perp }\rho } } \end{equation} \tag{ 15 }$$

$$\begin{equation} \Gamma _{{\rm HOL}\pm } = \frac{1}{C_{{\rm HOL}}} \sqrt{\frac{E_W^{\pm }B}{\rho }}. \end{equation} \tag{ 16 }$$

The free parameters used to set the strength of the dissipation are C_HOL and C_SW. In all three simulations, the values were chosen to obtain a maximum fast solar wind velocity of 750 km s⁻¹ at 0.1 AU. In the two simulations containing only one dissipation, this selection was straightforward: in the simulation with only Hollweg dissipation, $C_{{\rm HOL}}=2.1\,R_{\odot } \sqrt{G}$, and in the simulation with only SAW, C_SW = 10. In the combined case, for each value of C_SW there is a unique C_HOL that results in a maximum wind speed of 750 km s⁻¹. In this work, the goal is to determine the role of SAW damping in corona heating, and so we chose to keep that damping strong, and introduce a weak Hollweg dissipation. We diminished SAW dissipation by 20%, and then determined the value of C_HOL with the solar wind velocity ($C_{{\rm HOL}}=10.3\,R_{\odot } \sqrt{G}$ and C_SW = 12). The damping rate for Hollweg dissipation depends on the square root of the wave energy density, whereas SAW damping is independent of the wave energy density. Therefore, the simulation is more sensitive to the free parameter of one dissipation than the other. This explains why the free parameters changed by different factors (1.2 as compared to 5). To see how the behavior of the solution changes with the free parameter, see Appendix C.

We could have presented several combination simulations in which both C_SW and C_HOL were varied. However, we do not because we are not trying to constrain realistic values for these parameters in this study; this is outside the scope of this work for several reasons. We do not have a frequency-resolved spectrum (so we cannot treat ion-cyclotron damping), or a purely observationally driven boundary conditions (an empirical model is used for wave energy density), plus waves are not treated in closed field regions. We do not propose that a free parameter of 12 is the correct answer in the solar wind; we only use this value to show how the mechanism behaves in a 3D environment, which will be shown in Section 3.

The computational domain is a Sun-centered 48 × 48 × 48 R_☉ Cartesian grid, decomposed of blocks with 4 × 4 × 4 cells. The smallest cell size is 3/128 R_☉, in a shell of thickness 0.3 R_☉ beginning at the inner boundary. The heliospheric current sheet (HCS) is refined to 3/32 R_☉. The total number of cells in the simulation is 2.4 × 10⁶.

In this section, we present three simulations: one for each dissipation mechanism individually, and one which included both SAW and Hollweg dissipation (described in Sections 2.2 and 2.3). All simulations have the same boundary conditions, as described in Section 2.3. First we present the large-scale simulation results, and then the low corona (R ⩽ 1.225 R_☉). We include the DEMT reconstructions and Ulysses data as evidence that SAW dissipation occurs and is an important effect. Finally, we compare the Alfvén speed profile from the new solar wind model to the previous solar wind model in the SWMF, and discuss the implications for modeling CMEs in the lower corona.

3.1. Corona and Inner Heliosphere

In Figure 3, the dissipation rate for waves (Γ as defined in Equations (15) and (16)) is shown in the X = 0 plane (−12 to 12 R_☉ in each direction). The values for SAW dissipation are shown in panels (a) and (b) and for Hollweg dissipation in panels (c) and (d). Panels (a) and (c) correspond to waves traveling along the magnetic field, and the opposite is shown in panels (b) and (d). The magnetic field projected onto the plane is shown as black lines.

Zoom In Zoom Out Reset image size

Although we parameterized the form of the SAW damping rate, it can been seen in Figure 3 that the behavior of the mechanism is as expected for solar minimum conditions. Radial and latitudinal gradients in density create dissipation along mid-latitude (40°–70°) open magnetic field lines for R ⩽ 2.5 R_☉. We expect dissipation in this region due to the superradial expansion of the magnetic field, which is captured by density gradients. Also as expected, there is reduced dissipation above the poles, where the field line expansion is mostly radial. Latitudinal density gradients produce SAW dissipation near the solar surface at the boundaries of the open and closed magnetic field lines, and extend into the corona surrounding the HCS. Hollweg dissipation occurs most strongly at high latitudes for R ⩽ 4 R_☉, where the magnetic field (and thus the wave energy density) is largest.

For both mechanisms, the dissipation rate differs for the two waves in the simulation (traveling along or opposite to the magnetic field). For the Hollweg dissipation, it is due to the dependence of Γ on the wave energy density. For SAW damping, it occurs due to the dependence on the group velocity of the wave.

Figure 4 shows the solar wind temperature at R = 0.1 AU. The panels correspond to simulations with different wave dissipation choices. The Hollweg dissipation only (panel (a)) produces the most heating directly above the poles. And, because it scales with the magnetic field, Hollweg dissipation produces a temperature minimum at the HCS. SAW dissipation (shown in panel (b)) provides heating in two locations: the subpolar superradially expanding open field lines, and at the boundaries between open and closed magnetic fields. Combining the mechanisms (panel (c)) results in heating at high latitudes (from Hollweg dissipation), subpolar latitude heating (from SAWs), and heating around the HCS (due to SAWs; although in the literature temperature fluctuations near the HCS are attributed to turbulence). This temperature profile is most aligned with observations. Using ion temperature data measured during the first fast latitude scan of Ulysses, we scaled the latitudinal temperature distribution in the fast solar wind to 1 AU (details in Appendix B). The slope of a linear fit applied to the data was positive in the direction from the pole toward the equator for both hemispheres. This result suggests that there could be a heating mechanism at work in the solar wind which is efficient at the boundaries between fast and slow solar wind (or, at the Sun between open and closed magnetic fields).

Zoom In Zoom Out Reset image size

In Figure 5, we show the large-scale solar wind structure from the simulation in which both SAW and Hollweg dissipation were included (corresponding to Figure 4(c)). The meridional slice (−24 to 24 R_☉ in each direction) gives the magnitude of the solar wind speed. The sphere in the center represents the inner boundary, from which selected magnetic field lines are drawn (shown in black). The color contour on the solar surface is the radial magnetic field, with the color bar corresponding to Figure 2(a). The sharp transition between fast and slow solar wind measured by Ulysses is not reproduced in the model; in all simulations the transition occurs over 25° of latitude. The fast solar wind is accelerated along open field lines by the wave pressure gradient, and therefore a bimodal wind is achieved independent of the wave dissipation mechanism. However, the dissipation rate does affect the magnitude of the solar wind speed (see discussion in Appendix A).

Zoom In Zoom Out Reset image size

3.2. Low Corona

In Figure 6, meridional slices are shown from 1.035 R_☉ ⩽ R ⩽ 1.225 R_☉. Panels correspond to simulations with different wave dissipation choices: panel (a): Hollweg dissipation only and panel (b): SAW dissipation only. The contour is the ratio of the electron temperature reconstructed using the DEMT technique applied to STEREO/EUVI data during CR2077 (Vásquez et al. 2010) to the simulation temperature. The white spots in this region indicate locations where the reconstruction method failed (perhaps due to time-dependent phenomena). The magnetic field projected onto the plane is shown as black lines: the outer lines are open field lines, and the inner two lines are a closed streamer line. The thick black boxes emphasize the boundary between open and closed field lines in the northern hemisphere.

Zoom In Zoom Out Reset image size

The DEMT reconstructions show clear regions of temperature enhancement at latitudes of ±60°, which correspond to the boundaries of open and closed magnetic fields in the simulation. In panel (a), the temperature ratio is large in this region, which demonstrates that Hollweg dissipation does not generate appropriate heating at these locations (as can been seen by the red spots in the northern and southern boundaries). In panel (b), it can be seen that SAW damping heats where the DEMT reconstruction shows enhanced temperatures. These enhanced temperature regions are reproduced in location and nearly in magnitude by SAW damping, indicating that the mechanism is likely the dominant dissipation mechanism in this region.

3.3. Profile of Alfvén Speed

Observations of energetic electrons and type II radio bursts indicate that CME-driven shocks can accelerate particles very close to the Sun (Klassen et al. 2002). It is well established that interplanetary CMEs drive shocks at larger distances (Gopalswamy & Kaiser 2002; Claßen & Aurass 2002); however, the radio data show a disconnect between shock formed in the low corona and interplanetary shocks (Cane & Erickson 2005). Gopalswamy et al. (2001) and Mann et al. (2003), among others, explain these observations with the combination of CME acceleration while it propagates through an environment in which the Alfvén speed has a maximum, or hump. The scenario is that CMEs can drive strong shocks close to the Sun which then dissipate around 4 R_☉ due to the maximum in the Alfvén profile. The shock could form again above 6–10 R_☉ and propagate through the heliosphere.

Evans et al. (2008) presented Alfvén speed profiles from several solar wind models in the lower corona. It was found that global MHD models with ad hoc heating generated an Alfvén speed profile which had low maximum values. If the profile on the Sun were low and flat, as in some of the models, then we would expect that slow disturbances would drive strong shocks low in the lower corona, and these shocks would propagate out as interplanetary shocks (in contrast with observations discussed above). Evans et al. (2008) showed that 1D solar wind models which included Alfvén waves and physically motivated Alfvén wave dissipation produced steepened Alfvén profiles, and better agreement with observations. They proposed that the inclusion of Alfvén waves with physically motivated dissipation would steepen the Alfvén speed profile in a global MHD model.

Figure 7 shows Alfvén speed profile with radial distance from an active region. The black line was calculated from the semi-empirical model presented in Mann et al. (2003), specific to the active region in the simulation (see Appendix C for details of the model). The blue line was taken from the simulation in this paper with both SAW and Hollweg dissipation. The red line was taken from a simulation using the model of Cohen et al. (2007), in which the solar wind is heated by a spatially varying polytropic index. Overall, the inclusion of Alfvén waves results in better agreement with the features of the semi-empirical Mann et al. model. It can be seen that the profile of the Alfvén speed is steepened by the inclusion of waves: the slope from the valley to the hump is 139 km s⁻¹ R⁻¹_☉ for the Alfvén wave simulation, and 70 km s⁻¹ R⁻¹_☉ for the polytropic simulation. Also, the value at the hump for the Alfvén speed is over 500 km s⁻¹ for the wave-driven model, 150 km s⁻¹ larger than the hump in the polytropic profile. The speed at the valley for the polytropic model is less than 150 km s⁻¹, whereas the wave-driven Alfvén speed is less than 300 km s⁻¹ over only a short distance. The wave-driven profile inhibits shock formation by slow drivers (v < 300 km s⁻¹), which is more aligned with the study presented in Gopalswamy et al. (2005), which found that only 5 of 72 sampled CMEs that were a source of a large SEP event in the low corona had v < 1000 km s⁻¹.

Zoom In Zoom Out Reset image size

The goal of this paper is to determine in which regions SAW damping could be contributing to the heating of the solar corona. To achieve this goal, we take advantage of the development of a solar wind model driven by Alfvén waves with physical damping mechanisms in the SWMF. Whereas other modeling efforts have focused on wave dissipation mechanisms which dominate only in closed magnetic field regions or along open polar field lines in only 1D, we emphasize the interface region in a 3D model. This region is characterized by gradients in density perpendicular to the magnetic field. These gradients drive the dissipation of SAWs by resonant absorption. In this work, we presented results from the first self-consistent implementation of SAW damping in a 3D global MHD model.

What distinguishes the solar wind model presented in this paper from others in the literature is the combination of (1) data-driven boundary conditions, (2) a momentum source of wave pressure, and (3) the self-consistent calculation of an additional non-MHD energy source through physically motivated wave dissipation, including the first implementation of multiple wave dissipation mechanisms in a 3D environment. Although several global models now apply complex energy equations which treat physical processes such as electron heat conduction and radiative cooling, all include empirical source terms to achieve coronal heating (Lionello et al. 2009; Downs et al. 2010; Usmanov et al. 2011; Riley et al. 2011; Airapetian et al. 2011). Several of these works showed that a combination of heating terms (formulated to operate in quiet Sun, active regions, and open field regions differently) was required to match observations such as emission measurements.

In this work, we show that the same approach must be taken in wave-driven models: a combination of dissipation mechanisms, which focus in different regions, is required in order to reproduce the large-scale solar corona. Additionally, the physically motivated wave dissipation presented in this paper showed the ability to reproduce observational signatures, which removes the need for ad hoc coronal heating terms. However, the SAW dissipation presented in this paper was calculated with a phenomenological dissipation rate, which contained one free parameter. As we use the damping to accelerate and heat the wind only, we constrained this free parameter with the maximum solar wind speed, and therefore this free parameter should not be given any physical meaning. In future studies, we will improve the representation of the dissipation, and the corresponding free parameter will be better constrained, such as with the overall profile of the solar wind, as well as specifics like the wave amplitude, and local plasma properties at the inner boundary of the simulation.

We did not expect a simulation with only SAW damping to reproduce the details of the solar wind, just like a simulation with only Hollweg dissipation (van der Holst et al. 2010; Jin et al. 2012) does not reproduce heating at the boundary of open/closed magnetic field. Including SAW dissipation in the model of van der Holst et al. (2010) improved agreement with observations of coronal temperature both near the Sun and in the inner heliosphere. The electron temperature in the low corona (R < 1.225 R_☉) from DEMT reconstructions has been shown to contain enhanced temperature regions near the boundaries of open and closed magnetic fields (Vásquez et al. 2010). The simulations with SAW produce regions of enhanced temperature which match fairly well, both in location and magnitude. In the heliosphere, we have argued that a fit of the scaled Ulysses temperature data showed enhanced temperature at the boundaries of fast and slow solar wind, which is a signature of SAW dissipation. However, SAW damping alone cannot reproduce the high temperature structure in the open field regions, as they do not heat directly above the poles. It should be combined with other mechanisms to gain a complete picture of the corona.

We demonstrated that the inclusion of SAW dissipation improved the usability of the solar wind model in studying transients, which is crucial for space weather forecasting. As proposed in Evans et al. (2008), the inclusion of Alfvén waves with a physically motivated dissipation mechanism steepened the Alfvén speed profile in the low corona, aligning the Alfvén profile better with observational constraints. This improved Alfvén profile makes this model a more viable solar wind background in which to simulate CME-driven shocks than a polytropic solar wind model.

The next step in developing a physically motivated solar wind model is to treat wave spectral evolution. This aspect has been implemented in an axisymmetric stellar wind model (Airapetian et al. 2010), and is under development in the SWMF (Oran et al. 2010). A spectrum of Alfvén waves is required to treat ion-cyclotron damping, thought to be prevalent in coronal holes (Cranmer et al. 1999; Isenberg et al. 2001; Chandran et al. 2010). Additionally, relaxing the WKB limit (to account for counterstreaming waves) will allow for self-consistent calculation of nonlinear turbulent dissipation of in and outgoing waves (Verdini & Velli 2007; Cranmer et al. 2007; Chandran & Hollweg 2009). Including wave propagation along closed loops in the SWMF will allow us to model SAW damping in this region as well. From CoMP data, Verth et al. (2010) made new estimates of resonant damping along closed loops, and we will compare modeled estimations to this data set. Also, the temperature profiles from the DEMT reconstructions have been used to show that closed loops in which the temperature increased with height were located near the boundaries of open and closed fields lines, which could be due to resonant absorption. These aspects will be pursued in a future study.

We thank NASA Ames for use of the Pleiades supercomputer. R.E. thanks an anonymous referee for the comments and suggestions that greatly improved the manuscript. This research was supported in part by NSF CAREER Grant ATM 0747654, and in part by an appointment to the NASA Postdoctoral Program at Goddard Space Flight Center, administered by Oak Ridge Associated Universities through a contract with NASA.

In the formalism presented in this paper, wave dissipation by SAW damping contains a free parameter (C_SW). This parameter is used to gauge the strength of the dissipation. In this section, we investigate the effect of changing this quantity on the solar wind solution.

In Figure 8, the temperature (solid lines) and solar wind speed (dashed lines) are plotted as a function of heliocentric distance, extracted along a radial line at 60° latitude. We chose this location because SAW dissipation is strong there. The thicker lines correspond to a simulation with stronger wave dissipation (C_SW = 8) and weaker dissipation in the thinner lines (C_SW = 10).

Zoom In Zoom Out Reset image size

The effect of increasing the strength of SAW damping on the solar wind speed is to increase its value near the Sun and reduce its terminal value. This behavior is in agreement with the expected behavior of depositing energy in a stellar wind below the sonic point (Hartmann & MacGregor 1980). We find that the temperature is also larger close to and smaller far from the Sun when the dissipation rate is increased.

Ulysses data in the heliosphere indirectly support the presence of low (subpolar) latitude heating. Following Goldstein et al. (1996), we scaled 1 hr averaged Solar Wind over the Sun (SWOOPS) ion radial temperature data (from the first fast latitude scan, 1994–1995, solar minimum conditions) by performing a power-law fit to data for which v > 700 km s⁻¹ and R = 1.55–3.03 AU, separately in each hemisphere. We found a temperature falloff of r^−0.921 and r^−1.083 in the south and north poles, respectively. (The values reported in Goldstein et al. (1996) were r^−0.81 and r^−1.03.) The Ulysses data (scaled to 1 AU) are plotted as the solid black line in panel (a) of Figure 9. We applied a linear fit (yellow lines) to the scaled Ulysses data to determine whether the latitudinal temperature profile in the fast solar wind, from the poles toward the equator, decreased (as would be expected if heating is strongest at the poles) or increased (as would be expected if most heating occurs near the poles during solar minimum conditions). We found that in both the northern and southern hemispheres, the temperature is almost constant with latitude, except for a slight increase toward the boundary. The slopes of the fit were 227 K per degree in latitude in the northern hemisphere, and 357 K per degree in latitude in the southern hemisphere. This behavior supports heating not being concentrated at the poles, but occurring also at the boundaries of fast and slow solar wind, which we connect to SAW dissipation close to the Sun (as calculated in this paper).

Zoom In Zoom Out Reset image size

The fitting we performed is different than what was originally presented in McComas et al. (2000). The value reported there was the average latitudinal temperature, computed over eight segments of Ulysses's orbit, from 1992 February to 1997 December. For protons, the temperature decreased from the pole toward the equator: the value reported was an average of 2.7 × 10⁵ K at 60°, with a latitudinal dependence of +223 K per degree in latitude (i.e., increasing toward the poles). For alpha particles, however, the temperature increased from the pole to the equator: the value reported was an average of 1.4 × 10⁵ K at 60°, with a latitudinal dependence of −871 K per degree in latitude.

Although the temperature at R = 0.1 AU has not been measured in situ, we can make an estimate of the temperature in the fast solar wind by assuming that the radial behavior we calculated at 1 AU can be used to scale the temperature to 0.1 AU. The result is T = 2.8 MK. In the simulation which includes SAW and Hollweg dissipation, the temperature at 0.1 AU at 60° latitude is 0.6 MK, which suggests that the radial scaling is not valid closer in.

In Section 3.3, we compared the Alfvén speed profile from our model with a semi-empirical model (herein MAN03). MAN03 defines a profile of Alfvén speed near an active region as a function of height r (in units of solar radii) above the solar surface (Mann et al. 2003). The magnetic field strength is expressed analytically as an active region field superimposed on quiet-Sun dipole:

$$\begin{equation} |B| = B_1*(1-r)^{-3} + B_2* r^{-2}, \end{equation} \tag{ C1 }$$

where B₁ = 0.005 G to match the strength of the active region in CR2077 (40 G). The second term was set B₂ = 2.2 G, a typical value for the quiet-Sun field. Above R = 1.8 R_☉, the quiet-Sun field dominates.

The density in MAN03 is modeled as a one-fold Newkirk model (Newkirk 1961) for R < 1.8 R_☉,

$$\begin{equation} N(r) = 1.92*4.2 \times 10^4 (10^{4.32*r^{-1}}), \end{equation} \tag{ C2 }$$

above which they use the empirical form of Mann et al. (1999),

$$\begin{equation} N(r) = 1.92*5.14 \times 10^9 \exp (13.8*(r^{-1}-1)). \end{equation} \tag{ C3 }$$