THE ROLE OF SUBSURFACE FLOWS IN SOLAR SURFACE CONVECTION: MODELING THE SPECTRUM OF SUPERGRANULAR AND LARGER SCALE FLOWS

Equation (2) can be written as

$$\begin{equation} i{\bf k}_h\cdot {\tilde{\bf u}}_h = -\frac{\partial \tilde{u}_z}{\partial z}-\frac{\tilde{u}_z}{H_{\rho }}, \end{equation} \tag{ 4 }$$

where the overlying tildes indicate the complex Fourier amplitudes resulting from a two-dimensional horizontal Fourier transform at each depth z and k_h is the horizontal mode wave vector. By squaring both sides and taking the two limits of Equation (4), we can define a relationship between the power in horizontal and vertical motions without directly solving for the phases of the modes. For modes smaller than the integral scale, we take the limit $({\partial \tilde{u}_z}/{\partial z}) \gg ({\tilde{u}_z}/{H_{\rho }})$, while for larger scale modes we take $({\partial \tilde{u}_z}/{\partial z}) \ll ({\tilde{u}_z}/{H_{\rho }})$. Even in these limits, defining the relationship between vertical and horizontal power is difficult for two reasons: when squaring Equation (4), the cross terms between horizontal modes $\tilde{u}_x$ and $\tilde{u}_y$ on the left side do not have an a priori known form, and the vertical derivative on the right-hand side cannot be simply related to the wavenumber of the horizontal Fourier modes.

To proceed we make simplifying assumptions that we have empirically verified to hold in stratified (J. W. Lord 2014, in preparation) and incompressible turbulence simulations (Mininni et al. 2006) as appropriate. At small scales, the flow is nearly isotropic and homogeneous, with unstratified homogeneous and isotropic turbulence simulations showing that $k_x^2\tilde{u}_x\tilde{u}_x^* \approx k_y^2\tilde{u}_y\tilde{u}_y^* \approx ({\partial \tilde{u}_z}/{\partial z})({\partial \tilde{u}_z^*}/{\partial z})$, which together with incompressibility, yields $({\partial \tilde{u}_z}/{\partial z})({\partial \tilde{u}_z^*}/{\partial z}) \approx ({1}/{4})k_h^2\tilde{{{\bf u}}}_h\cdot \tilde{{{\bf u}}}_h^*$. The simulations also suggest a relationship between the vertical and horizontal gradients, $({\partial \tilde{u}_z}/{\partial z})({\partial \tilde{u}_z^*}/{\partial z}) \approx ({1}/{4})k_h^2\tilde{u}_z\tilde{u}_z^*$, and together these yield a relationship between horizontal and vertical power:

$$\begin{equation} \tilde{{\bf u}}_h\cdot \tilde{{\bf u}}_h^*=\tilde{u}_z\tilde{u}_z^*, \end{equation} \tag{ 5 }$$

where $\tilde{{\bf u}}_h\cdot \tilde{{\bf u}}_h^*=\tilde{u}_x\tilde{u}_x^*+\tilde{u}_y\tilde{u}_y^*$. At large scales, the cross terms, which result from squaring the left-hand side of Equation (4), are measured in stratified simulations to be small and are set to zero. This implies that $k_x^2\tilde{u}_x\tilde{u}_x^* + k_y^2\tilde{u}_y\tilde{u}_y^* \approx \tilde{u}_z\tilde{u}_z^*/{H_{\rho }^2}$, and the horizontal and vertical power in the modes are related as

$$\begin{equation} \tilde{{\bf u}}_h\cdot \tilde{{\bf u}}_h^*= \frac{2}{k_h^2H_{\rho }^2}\tilde{u}_z\tilde{u}_z^*, \end{equation} \tag{ 6 }$$

where $k_h^2=k_x^2+k_y^2$.

Finally, without approximation, the driving scale that separates these two behaviors (Equation (3)) can be rewritten as

$$\begin{equation} \lambda _h=4\alpha H_{\rho } \frac{u_h}{\vert u_z\vert }, \end{equation} \tag{ 7 }$$

where λ_h = 2π/k_h is the wavelength of the Fourier mode corresponding to a convective cell diameter of d = 2r. By taking α and u_h/u_z ≈ 1 we approximate the driving scale as λ_h ≈ 4H_ρ. It is on the basis of these relationships (Equation (5) at small scales and Equation (6) at large scales with the crossover between them given by the driving scale λ_h = 4H_ρ at each depth) that we calculate the horizontal velocity power spectrum from the vertical.

Analysis of large-scale radiative hydrodynamic simulations of solar convection (for details, see Section 5.1 and J. W. Lord 2014, in preparation) helps validate these relationships. Below 1.3 Mm beneath the photosphere, the mass continuity in the simulations matches the anelastic balance (Equations (2) and (4)) to within a few percent. Near the photosphere, this balance breaks down because of fluid compressibility, particularly at high wavenumbers. We thus restrict our model analysis to depths below 1.3 Mm. At low wavenumbers, p-mode contributions can still be important at the shallowest depths. We remove these when comparing the numerical simulations to the model by averaging the simulation velocities over 30 minutes. This averaging also reduces the amplitude of the high wavenumber convective motions, but preserves the relationships between horizontal and vertical flows of Equations (5) and (6). This is illustrated by Figure 1, in which the horizontal velocity spectra measured at several depths in a solar-like radiative hydrodynamic simulation are plotted. Overplotted are the horizontal velocity spectra deduced from the vertical velocity spectra of the simulation at the same depths using Equations (5) and (6) (dotted and dashed line styles, respectively). The two-component reconstruction of the horizontal velocity spectrum from the vertical reproduces the shape and amplitude of the actual spectrum quite well. Moreover, the driving scale estimate of 4H_ρ is in good agreement with the crossover between the two behaviors. Plotted in Figure 2 is the crossover wavenumber as a function of depth (defined as the smallest wavenumber in the simulations for which the balance in Equation (6) begins to fail, meaning that the difference between the two terms at neighboring larger horizontal wavenumbers is increasingly large). For comparison, 4H_ρ is overplotted in red. They are in good agreement. Note that the discontinuities in the measured values are due to the finite spectral resolution of the simulation; many depths in the simulation appear to have the same crossover scale because there are no modes that can discriminate between them.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We employ a simplified hydrostatic mixing length atmosphere of pure hydrogen and helium in Saha equilibrium, integrating the mixing length equations (Prandtl 1925; Böhm-Vitense 1958) using the observed density and temperature of the solar photosphere as boundary conditions. We take the convective flux to be equal to the photospheric radiative output (6.3 × 10¹⁰ erg cm⁻² s⁻¹) as appropriate for efficient convection and note that this introduces an error in the lower portions of the model where, in the Sun, radiation transports a significant fraction of the energy flux. This error makes a small contribution to the excess model power at the largest scales (Section 4).

Explicitly, we solve the equation for the convective energy flux F_c = (1/2)ρvC_pT(l/H_p)(∇ − ∇') along with that of hydrostatic balance (dP/dz) = −ρg. Here ρ is the fluid density, v is the velocity, C_p is the specific heat at constant pressure, T is the temperature, ∇ is the mean temperature gradient, ∇' is the temperature gradient within the convective cell, P is the pressure, and g = Gm(z)/r(z)² is the gravitational acceleration with r(z) the distance from the Sun's center, m(z) the mass within that radius, and G the gravitational constant. We employ the equation of state P = ρkT/μ, where k is the Boltzmann constant and μ is the mean molecular weight of the plasma, and assume that the convective motions are adiabatic, so that ∇' = ∇_ad = (∂lnT/∂lnP)|_ad, the adiabatic temperature gradient. Finally, the rms convective velocity is given by v² = (1/8)gQ(l²/H_p)(∇ − ∇'), where l is the mixing length measured in units of the pressure scale height H_P. Note that Q = 1 − (∂lnμ/∂lnT)|_P, C_p, ∇_ad, and μ account for the non-ideal effects of hydrogen and helium ionization (where the number density of each ionization state is determined in collisional equilibrium as a Saha balance).

The equations are integrated numerically from the photosphere downward, yielding the convective rms velocity and the local density scale height at each depth.

In all simulations of solar convection, the amplitude of the vertical velocity at low wavenumbers decreases toward the surface where granular scale convection is dominant. The rate of this decrease for the largest scale convective modes is a fundamental uncertainty in our understanding of solar convection. While global simulations predict giant cell convection throughout much of the convection zone (e.g., Miesch et al. 2008), surface observations have only very recently found evidence for weak flows at these scales (Hathaway et al. 2013). Numerical simulations have difficulty directly addressing the photospheric amplitude of large-scale motions. They are either of limited extent in depth (Stein et al. 2009; Ustyugov 2010; J. W. Lord 2014, in preparation) or do not capture the non-ideal and highly compressible nature of the uppermost layers (Miesch et al. 2008).

Because of these uncertainties, we examine two possible vertical velocity amplitude decay profiles. Starting from the depth at which the wavelength of the mode exceeds the integral scale, we decay the modes either by approximating the flow as potential (van Ballegooijen 1986) or by using a cubic polynomial fit to the observed decay of modes with wavelengths between 20 Mm and 50 Mm in the hydrodynamic simulations (shown as diamonds in Figure 3). The potential flow approximation takes the flow to be irrotational, allowing a direct solution to the large-scale continuity balance, written as

$$\begin{equation} \frac{\partial \tilde{\phi }}{\partial z}=-{k}_hH_\rho \tilde{\phi }, \end{equation} \tag{ 8 }$$

where ϕ is the velocity potential with $\tilde{u}_z={\partial \tilde{\phi }}/{\partial z}$. This yields a profile for the velocity amplitude with height

$$\begin{equation} \tilde{u}_z(z)=\tilde{u}_z(z_d) \frac{H_\rho (z)}{H_\rho (z_d)} \exp \left[k_h^2\int _{z_d}^zH_\rho (z^{\prime })dz^{\prime }\right], \end{equation} \tag{ 9 }$$

where z_d is the driving depth. This velocity profile can be integrated numerically for any wavenumber k_h, the results of which are shown with triangles in Figure 3. The polynomial fit, on the other hand, approximates the decay of the modes by a single function determined from the hydrodynamic simulations (solid line in Figure 3). The fit groups the behavior of all modes between 20 and 50 Mm together and is thus inadequate to reproduce the hydrodynamic simulation in detail (see Section 4). It is employed in the model because of its simplicity. The two schemes are quite different in form and together provide a test of the sensitivity of the model to this key unknown function.

Zoom In Zoom Out Reset image size

In summary, we construct the model surface horizontal velocity spectrum as follows. To calculate the spectrum of the vertical velocity we follow these steps.

• 1.  
  Construct a mixing length model of the solar convection zone integrating from the photosphere downward to 200 Mm, the approximate depth of the solar convection zone.
• 2.  
  Determine the wavelength of the largest scale mode allowed at the bottom of the model atmosphere, λ_h = 4H_ρ, and use this as the integral (driving) scale (i.e., the lowest wavenumber mode) in a k^−5/3 turbulent cascade. The highest wavenumber in the spectrum is taken to be the Nyquist frequency of the hydrodynamic simulations on which the model is based (2π/384 km, see Section 5.1), and the spectrum is normalized so that the integrated power is equal to the mixing length velocity squared at the bottom of the model atmosphere.
• 3.  
  Move one step up in the atmosphere (a grid spacing of 64 km is used to again match the hydrodynamic simulations). Decay modes with wavelengths longer than the local integral scale (4H_ρ) using one of the two decay functions discussed in Section 3.2. Compute the integrated power in the decaying modes and normalize the remaining k^−5/3 spectrum by the squared mixing length velocity minus the power in the decaying modes.
• 4.  
  Repeat step 3 until the top of the model atmosphere is reached.

From the vertical velocity spectrum, the horizontal velocity spectrum at any height is computed using Equations (5) and (6).

Both the simplified model and the full three-dimensional radiative hydrodynamic spectra show more power than the Sun at scales larger than supergranulation, with that power increasing monotonically toward lower wavenumbers because large-scale flows are convectively driven in the deep layers of the domains. It is worth noting that if the solar spectrum matched either the simplified model or the radiative hydrodynamic simulation spectrum, giant cell convection would be relatively easy to observe as the power in these large-scale modes would exceed that in supergranulation.

To make a more direct comparison between our numerical simulations and observations, we employed a Coherent Structure Tracking (CST; Roudier et al. 2012) algorithm to infer the horizontal velocities on large scales from measurements of the motions of granules. In Figure 5 we compare the CST horizontal velocity spectra of a large-scale radiative hydrodynamic simulation (solid red curve) using the OPAL equation of state (Rogers & Iglesias 1992) with that of solar observations (dashed black curve) from the Helioseismic and Magnetic Imager (HMI) aboard the Solar Dynamics Observatory. The measurements in both cases employ a 22 hr sequence of continuum intensity images with each HMI image separated by 45 s and each simulation image separated by ∼40 s. We break this sequence into 11 two-hour subsets and use the CST method to compute the velocity for each 2 hr window. The spectrum shown is the average spectrum of those 11 velocity computations. The HMI observations are of a 192 × 192 Mm² region at disk center with a low magnetic activity on 2010 19 June. The simulation solution was computed using the MURaM code (Vögler et al. 2005) in a 196 × 196 × 49 Mm³ domain with 192 × 192 × 64 km³ grid spacing (J. W. Lord 2014, in preparation).

Zoom In Zoom Out Reset image size

We note that the spectra in Figure 5 are truncated at high wavenumbers because the CST method is not reliable for scales smaller than 2.5 Mm. Moreover, low wavenumber modes, those with wavenumbers below k < 0.013 Mm⁻¹ (indicated by the vertical fiducial line and dot-dashed line style in Figure 5) have length scales larger than the integral (driving) scale at the bottom of the simulation domain and consequently have lower photospheric amplitudes than they would likely have in a deeper simulation. Between these extremes are two notable mismatches between the simulation and observation spectra: the simulation shows an excess of power at low wavenumbers and a deficit of power at high wavenumbers when compared to observations. Our very wide and deep simulations resolve supergranular scale motions well but under-resolve granular motions. This leads to an inferred CST velocity with reduced power at high k, a result that is inconsistent with the actual simulation velocities and HMI observations. The excess of low wavenumber power is, on the other hand, a fundamental difference between the resolved motions in the simulation and those in observations and is robust as the CST constrains large-scale motions better than small-scale motions (Roudier et al. 2012). Thus, understanding the observed solar supergranulation spectrum requires understanding the origin of this low wavenumber reduction in power along with any mechanism that may enhance power at supergranular scales.

Our simplified mixing length model suggests that the low wavenumber vertical motions are driven deep in the convection zone and decrease in amplitude toward the surface. The radiative hydrodynamic simulations behave similarly (Section 4), and reducing the vertical flow velocities at depth reduces the low wavenumber horizontal velocity power in the simulated photosphere and improves the match between simulations and observations. We demonstrate this conclusively via simulations in which convective velocities in the deep layers are reduced without changing the mean stratification of the atmosphere (which is also fundamental to the surface spectrum). This was done using an artificial energy transport term. Specifically, we added an artificial flux function to the energy equation that depends only on depth. The artificial flux carries the full solar flux below a specified depth and none of the flux at heights above this. The hyperbolic tangent flux profile employed is 5.12 Mm wide centered at 10 Mm (where 4H_ρ ∼ 20 Mm), effectively supporting radiative losses from the photosphere by depositing the heat where the divergence of the function is nonzero. In Figure 5 (dotted blue curve) we plot the resulting photospheric horizontal velocity spectrum using the same CST method described above. There is substantially reduced power in the photosphere of the artificial flux simulation in those modes that are driven at depths below ∼10 Mm (scales larger than ∼20 Mm). This is the region of the domain for which the artificial energy flux is important and consequently convective (rms) velocities are reduced by a factor of ∼2.5.

The artificial energy flux experiment confirms the hypothesis that low wavenumber modes are driven deep in the simulated convection zone and imprint as horizontal flows in the surface layers. The photospheric power spectrum reflects a hierarchy of driving scales with depth even in fully nonlinear radiative hydrodynamic simulations. It also suggests that neither the radiative hydrodynamic solutions nor the simplified model spectrum capture the true dynamics of the solar convection below ∼10 Mm. In other words, in the Sun, low wavenumber flows carry much less of the convective energy flux or transport the energy at substantially lower velocities than expected based on the simulations or the model. Flow/enthalpy correlations, essential to convective transport, may thus not be correctly captured by hydrodynamic simulations. This may be due to their limited resolution or result from the boundary conditions applied. For example, the open boundary condition commonly employed in radiative hydrodynamic simulations of photospheric convection may smooth perturbations in the inflowing plasma. Alternatively, magnetic fields, not included in the simulations we have discussed in this paper, may maintain flow correlations and allow convective transport on smaller scales or at lower velocities than predicted by purely hydrodynamic models. Preliminary results from magnetized simulations favor this hypothesis, though the underlying mechanisms are still under investigation and the effect so far appears insufficient to explain solar observations (J. W. Lord 2014, in preparation).

The small plateau of power at supergranular scales (solid red curve in Figure 4(a) at k_h ∼ 0.04 Mm⁻¹) reflects the role of helium ionization in determining the convective velocities at depth. Superimposed on the horizontal velocity spectrum in Figure 4(a) we have plotted fiducial vertical lines to highlight the integral (driving) scale at the depths of 50% He i (k/2π = 0.1 Mm⁻¹, 6 Mm depth) and He ii (k/2π = 0.03 Mm⁻¹, 17.5 Mm depth) ionization. The plateau of supergranular power falls between these two fiducial lines. We have also computed the horizontal velocity spectrum for mixing length background atmospheres with equations of state that do not allow He ii or both He i and He ii ionization (Figure 4, dashed (green) and dotted (blue) curves, respectively). These test atmospheres show a continuous power law increase toward low wavenumbers with no feature at supergranular scales. More precisely, these spectra do not show the suppression of power at scales corresponding to the integral scale at the depths of 50% He i or He ii when the ionization processes are disallowed. The differences between the ionizing and non-ionizing spectra at still lower wavenumbers result because the stratification in the deep layers lies along a different adiabat. The velocity differences at depths below helium ionization, reflected in the low wavenumber horizontal velocity spectra in the near surface, are due to differences in the mean stratification as the medium is nearly fully ionized. In the region of partial ionization, convective velocities are also influenced by the availability of ionization energy in heat transport via perturbations about the mean ionization state.

Helium ionization is thus responsible for the small supergranular plateau in the model spectrum, albeit in a curious fashion. The ionization of helium yields a slight reduction in the driving scale mode amplitudes in the partially ionized regions (where the driving scale λ_h ∼ 10 Mm for He i and λ_h ∼ 35 Mm for He ii), producing a small apparent enhancement of power in the upper layers at wavenumbers that lie between these driving scales (where λ_h ∼ 20 Mm). The reduction in mode amplitudes results because ionization energy contributes to the heat transport. In a partially ionized fluid the heat can be transported by ionization state perturbations as well as thermal perturbations (Rast et al. 1993; Rast & Toomre 1993), and convective velocities in the mixing length model are thus reduced in partially ionized regions. (We note that the mixing length model is a local transport model and does not take into account other effects of ionization such as the increased linear (Rast 1991) and nonlinear (Rast 2001) instability of the fluid, though these may play a role in solar convective flows or in our more complete three-dimensional simulations.) The vertical velocity amplitudes of modes that begin their decay at the depths of helium ionization (the integral or driving scale modes at those depths) are thus suppressed, resulting in a reduction in their horizontal velocity power near the surface. Since ionization energy transport depends on ionization state perturbations, with transport in a fully neutral or fully ionized plasma behaving as an ideal gas, modes with peak amplitudes (those with integral scales equal to 4H_ρ) at depths that lie between the partially ionized regions (10 Mm, for example, where λ_h ∼ 20 Mm and k/2π = 0.04 Mm⁻¹) have more power than neighboring modes.

This is illustrated by Figure 4(b), which shows the vertical velocity power at the driving depth (i.e., the depth where the wavelength of the mode is equal to the integral (driving) scale) for mixing length model atmospheres which do or do not allow He i or He ii ionization. The driving scale modes in the regions of partial ionization have lower amplitudes than those outside of it, with minima in the mode amplitudes occurring at the depths of 50% ionization when allowed. The role of hydrogen ionization is difficult to illustrate as its effect is dominant in the surface layers where hydrogen recombination supports radiative losses. The region of H partial ionization is broad in depth and integral to the structure of the radiative boundary layer, and experiments preventing H ionization dramatically alter the mean state of the atmosphere and result in dramatic changes in the velocity spectrum across a wide range of wavenumbers. Thus we do not explicitly consider an atmosphere that disallows hydrogen ionization, but it is clear from Figure 4(b) that the effect of hydrogen ionization on mode amplitudes overlaps that of He i (compare crosses (blue) and triangles (green)).

The simple model we have presented thus suggests that there is an apparent enhancement of photospheric power at ∼20 Mm scales that occurs because helium ionization reduces the flow speeds in the regions of partial ionization which suppresses power at larger (∼35 Mm from He ii) and smaller (∼10 Mm from He i) scales, not because He ii ionization enhances the driving of flows at this depth as has been previously suggested. This apparent enhancement, however, is much smaller than the increased power in observations of solar supergranulation. To investigate the suppression of photospheric power by helium ionization in the context of solar-like convection we use the same Saha equations of state described above in three-dimensional radiative hydrodynamic simulations. We use the MURaM (Vögler et al. 2005) code to run simulations that use 192×192 km² horizontal resolution and 64 km vertical resolution with 1024 × 1024 × 768 grid cells (giving a domain size of 196 × 196 × 49 Mm³). The results presented here are from more than 5 days of solar time after the simulation has reached a relaxed equilibrium (J. W. Lord 2014, in preparation).

Figure 6(a) shows a comparison of the photospheric horizontal velocity spectrum from three such simulations, one that allows H, He i, and He ii ionization (solid red curve), one in which only H and He i ionization are permitted (dashed green curve), and one with only H ionization (dotted blue curve). The resulting horizontal velocity spectra show suppression of photospheric power similar to that seen in the simplified model (Figure 4) when ionization is allowed. The two lowest wavenumber modes (dot-dashed line style) have driving depths outside of the simulation domain and are consequently unreliable and weaker than what would be expected in a deeper simulation. The modes with integral (driving) scales equal to 4H_ρ in the regions of partial helium ionization again have reduced amplitudes. This is particularly apparent at the wavenumbers corresponding to modes that peak in the He ii partial ionization region (k/2π near 0.013 Mm⁻¹ in Figure 6(a)) which is well separated from the effects of hydrogen ionization. Not allowing He ii ionization (dashed green and dotted blue curves) causes a small but significant elevation of photospheric power at those wavenumbers. Disallowing He i ionization (dotted blue curve) induces smaller differences due to the dominant role of hydrogen in the surface layers.

Zoom In Zoom Out Reset image size

The same reduction in the mode amplitude of the vertical velocity spectrum at the driving depths corresponding to partial helium ionization seen in the simplified model is apparent in these hydrodynamic simulations (Figure 6(b)). Modes with scales equal to 4H_ρ at the depths of partial He i and He ii ionization have reduced amplitudes, though this reduction is noisier in the simulation than in the mixing length atmosphere (we note that the rms velocity amplitudes, not shown here, also very clearly increase at the nominal ionization depths when ionization is disallowed). This is due to three primary factors: the spectral resolution of the simulation is limited by the domain width, the three-dimensional simulation is non-local which makes using a single depth a poor representation of the vertical velocity power that reaches the surface, and the intrinsic temporal variation in the power of the modes below ∼10 Mm is long compared to the 5 days of simulation time. Moreover, other nonlinear effects of ionization may play some role, as discussed above. These experiments do, however, confirm, in the context of fully nonlinear three-dimensional radiative hydrodynamic simulations, two important results of the simplified model: the horizontal velocity spectrum in the photosphere reflects the amplitude of the vertical velocity at depth and the reduced amplitude of vertical velocity in the regions of partial helium ionization plays a minor role in shaping the spectrum of supergranular flows at the surface.