STELLAR MIXING LENGTH THEORY WITH ENTROPY RAIN

Late-type stars such as our Sun have outer convection zones. The observed solar granulation is a surface manifestation of their existence. Solar granulation and the first few megameters (Mm) of the convection zone have been modeled successfully using mixing length theory (MLT) and numerical simulations with realistic physics included (Stein & Nordlund 1989, 1998; Vögler et al. 2005; Gudiksen et al. 2011; Freytag et al. 2012). As a function of depth, the simulations reproduce some essential features predicted by MLT, in particular the depth dependence of the turbulent root mean square (rms) velocity, ${u}_{\mathrm{rms}}\approx {({F}_{\mathrm{enth}}/\rho )}^{1/3}$, where ${F}_{\mathrm{enth}}$ is the enthalpy flux and ρ is the density. Simulations also seem to confirm an important hypothesis of MLT regarding the gradual increase of the typical convective time and length scales with depth. Given such agreements, there was never any reason to question our basic understanding of convection.

In recent years, local helioseismology has allowed us to determine subsurface flow velocities (Duvall et al. 1997; Gizon & Birch 2005). Helioseismic observations by Hanasoge et al. (2010, 2012) using the deep-focusing time–distance technique have not, however, detected large-scale convection velocities at the expected levels (Miesch et al. 2012); see also Gizon & Birch (2012) for a comparison with both global simulations and radiation hydrodynamics simulations with realistic radiation and ionization physics included. Greer et al. (2015) have suggested that the approach of Hanasoge et al. (2012) may remove too much signal over the time span of the measurements. Using ring-diagram analysis (Gizon & Birch 2005) with appropriately assembled averaging kernels to focus on the deeper layers, Greer et al. (2015)  instead find values of the turbulent rms velocities that are consistent with conventional wisdom. Moreover, they find that at large length scales corresponding to spherical harmonic degrees of 30 or below (corresponding to scales of about $140\,\mathrm{Mm}$ and larger), the rms velocity actually increases with depth. This in itself is remarkable, because velocity perturbations from deeper layers were expected to be transmitted to the surface almost unimpededly (Stix 1981), unlike perturbations in the heat flux, which are efficiently being screened by the convection (Spruit 1977). Thus, an increase in horizontal velocity power with depth is unexpected. However, van Ballegooijen (1986) pointed out that for the scale and depth of giant cells, the screening might be large enough to allow giant cell convection of $100\,{\rm{m}}\,{{\rm{s}}}^{-1}$ to result in only $10\,{\rm{m}}\,{{\rm{s}}}^{-1}$ at the surface, which would be compatible with observations (Hathaway et al. 2013). The presence of the near-surface shear layer of the Sun (Schou et al. 1998) might enhance the screening further.

Given that the deeply focused kernels used by Greer et al. (2015) can still have 5%–10% sensitivity to Doppler shifts arising from flows in the upper layers, such near-surface flows could still leave an imprint on the signal. This effect would be exaggerated further if the near-surface flows were stronger compared to those in deeper layers, as was theoretically expected. On the other hand, if convection in the surface layers is of smaller scale, the signal from those layers will to a large extent be averaged out despite its larger amplitude. It would obviously be important to examine this more thoroughly by applying the kernels of Greer et al. (2015) to realistic simulations. Unfortunately, this has not yet been attempted. Conversely, the deep-focusing technique of Hanasoge et al. (2012) could be extended to allow for imaging of the deeper flow structures and thereby a direct comparison of individual turbulent eddies with those detected by Greer et al. (2015).

It should be mentioned that small flow speeds of giant cell convection have been found by correlation tracking of supergranule proper motions (Hathaway et al. 2013). The typical velocities are of the order of $20\,{\rm{m}}\,{{\rm{s}}}^{-1}$ at spherical degrees of around 10. Those speeds are still an order of magnitude above the helioseismic upper limits of Hanasoge et al. (2012), but Hanasoge & Sreenivasan (2014) argue that these surface measurements translate to lower flow speeds in the deeper and denser layers when considering mass conservation. This argument may be too naive and would obviously be in conflict with the results of Greer et al. (2015), which show an increase with depth at low spherical harmonic degrees. On the other hand, Bogart et al. (2015) have argued that the flows reported by Hathaway et al. (2013) may be non-convective in nature and, in fact, magnetically driven and perhaps related to the torsional oscillations.

Quantitatively, the horizontal flow speed at different scales is characterized by the kinetic energy spectrum, E(k), where k is the wavenumber or inverse length scale. Stein et al. (2007) have shown that the spectra of surface Dopplergrams and correlation tracking collapse with those of Stein & Nordlund onto a single graph such that the spectral velocity ${[{kE}(k)]}^{1/2}$ is proportional to k, i.e., $E(k)\propto k;$ see also Nordlund et al. (2009). This is a remarkable agreement between simulations and observations. While the origin of such a spectrum is theoretically not understood, it should be emphasized that these statements concern the horizontal velocity at the surface and do not address the controversy regarding the spectrum at larger depths, where the flow speeds at large length scales may still be either larger (e.g., Greer et al. 2015) or smaller (e.g., Featherstone & Hindman 2016) than those at the surface. If most of the kinetic energy were to reside on small scales throughout the convection zone even at a depth of several tens of megameters and beneath, E(k) is expected to decrease toward smaller wavenumbers k either like white noise ∝k² or maybe even with a Batchelor spectrum ∝k⁴ (Davidson 2004). However, even in simulations of forced isothermally stratified turbulence, in which there is no larger-scale driving from thermal buoyancy, there is a shallower spectrum, $E(k)\propto {k}^{3/2}$ (Losada et al. 2013). By contrast, if one imagines the flow to be anelastic so that the mass flux $\rho {\boldsymbol{u}}={\boldsymbol{\nabla }}\times {\boldsymbol{\psi }}$ can be written as a vectorial stream function ${\boldsymbol{\psi }}$, and if ${\boldsymbol{\psi }}$ is given by white noise, one should expect a k⁴ Batchelor spectrum. This is qualitatively similar to the results of granule tracking, which reveal an intermediate scaling proportional to k³ (Rieutord et al. 2008; Roudier et al. 2012). While the steeper k⁴ spectrum for $k\lt {k}_{{\rm{f}}}$ might leave some hope that the results of Hanasoge et al. (2012) could be reconciled with these theoretical and observational constraints, this would be virtually impossible for the linear energy spectrum, $E(k)\propto k$.

The issue has recently been examined by Lord et al. (2014), who have shown that simulation results with the MURaM code (Vögler et al. 2005) can be reproduced with a model that is composed of a continuous hierarchy of layers, each with its own driving scale at a scale of four local density scale heights. At wavenumbers below that driving scale, the spectral power falls off in a certain way. The horizontal surface spectrum is a superposition of these contributions, each of which is assumed to decay with height, contributing therefore progressively less with depth. However, both their model and the MURaM results (Lord 2014) show an order of magnitude more power (factor 2.5 in rms velocity) than the Sun at small wavenumbers. To understand this, they show that the observations can be reproduced if below a depth of $10\,\mathrm{Mm}$ the convective energy flux is augmented by an artificially added flux term, which thus significantly reduces the rms velocity of the resolved flow field. They speculate that the flow–temperature correlation entering the enthalpy flux may be larger in the Sun, possibly caused by a magnetic field that may maintain flow correlations and boost the transport of convective flux by smaller scales. This has been partially confirmed by Hotta et al. (2015), who have discussed the possible role of small-scale magnetic fields in suppressing the formation of large-scale flows. Furthermore, Featherstone & Hindman (2016) have found that with increasing Rayleigh numbers, there is more kinetic energy at small scales and less at large scales such that the total kinetic energy is unchanged by this rearrangement of energy.

Given the importance of this subject, it is worthwhile reviewing possible shortcomings in our theoretical understanding and numerical modeling of stellar convection. Both global and local surface simulations of solar-type convection would become numerically unstable with just the physical values of viscosity and radiative diffusivity. Therefore, the radiative flux (in the optically thick layers) is modified in one of two possible ways. (i) The contribution from temperature fluctuations is greatly enhanced, i.e.,

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{rad}}=-K{\boldsymbol{\nabla }}T\quad \to \quad -K{\boldsymbol{\nabla }}\overline{T}-{K}_{\mathrm{SGS}}{\boldsymbol{\nabla }}(T-\overline{T}),\end{eqnarray} \tag{ 1 }$$

where T and $\overline{T}$ are the actual and horizontally averaged temperatures, respectively, K is the radiative conductivity, and ${K}_{\mathrm{SGS}}$ is a subgrid scale (SGS) conductivity. The latter is enhanced by many orders of magnitude relative to K. Furthermore, numerical diffusion operators often do not translate in any obvious way to the physical operators. (ii) Alternatively, in direct numerical simulations (DNS), one uses physical viscosity and diffusivity operators, i.e., the replacement in Equation (1) is not invoked, but the coefficients are enhanced ($K\to {K}_{\mathrm{enh}}$) and exceed the physical ones by many orders of magnitude. Both approaches are problematic.

In many global DNS with enhanced coefficients (e.g., Käpylä et al. 2013), the lower boundary is closed, so at the bottom of the domain all the energy is carried by radiation alone. By choosing an enhanced radiative diffusivity, the radiative flux is increased by a corresponding amount, and therefore also the total flux. The luminosity in those simulations can exceed the solar value by several orders of magnitude. This has a series of consequences. Most notably, the convective velocities are too high in the upper parts where most of the flux is carried by convection (Käpylä et al. 2013). There are two ways to avoid this problem. One is to use simulations with a polytropic hydrostatic reference solution, whose polytropic index is close enough to the adiabatic one so that the convective flux is everywhere a small fraction of the radiative one (Brandenburg et al. 2005), which reduces the convective velocities. Alternatively, one applies the enhanced radiative conductivity only to the temperature fluctuations so as not to disturb the very small radiative energy flux, $-K{\boldsymbol{\nabla }}T$, compared with $-{K}_{\mathrm{enh}}{\boldsymbol{\nabla }}T$, which it would have been in the DNS approach without subtracting $\overline{T}$. The temperature smoothing implied by invoking alternative (i) is also necessary in the simulations with realistic opacities, because the Péclet number (which is similar to the Reynolds number) based on the physical value of K would reach values above 10¹⁰, which cannot be handled by a DNS (Barekat & Brandenburg 2014). In the global simulations, it is common to use the specific entropy gradient instead of ${\boldsymbol{\nabla }}(T-\overline{T})$ (e.g., Käpylä et al. 2013). This is equivalent if $\overline{T}$ is close to the adiabatic value. However, the diffusion coefficient in the SGS term can easily be five orders of magnitude larger than the physical one acting on the mean stratification. This has the consequence of suppressing small-scale turbulent flows and entropy fluctuations, which may prematurely damp the low-entropy fluid parcels that originate within the strong cooling layer at the surface and which will be discussed as entropy rain in the bulk of this paper. If such low-entropy contributions are poorly captured by the simulation, this could be compensated for by a sufficiently strong contribution proportional to the superadiabatic gradient, which in turn would give rise to the excitation of flows on much larger scales than would be compatible with the observations discussed by Lord et al. (2014). The sensitivity of large-scale flow excitation to small changes in the superadiabatic gradient was also discussed by Cossette & Rast (2016), who studied models with different types of surface driving. It is further supported by the work of Featherstone & Hindman (2016), who found that there is more kinetic energy at small scales and less at large scales as the Rayleigh number increases.

There is yet another problem that concerns the global simulations in which a predetermined profile K(z) is used instead of calculating its local value with a physical opacity. Such an approach has been used routinely in studies of compressible convection, especially when stably and unstably stratified layers are combined (Hurlburt et al. 1986; Brandenburg et al. 1996). The simulations of Käpylä et al. (2013) adopted a profile that yields a negative (unstable) radial entropy gradient through most of the layer (see the inset of their Figure 3). However, as will be pointed out below, it is only a tiny surface layer in which the non-convective model is unstable. The rest of the model is a priori stable to convection, but it becomes marginally stable (or perhaps slightly unstable) as a consequence of the resulting turbulence leading to bulk mixing across the deeper layers. This is quite contrary to the models with a predefined K(z) profile, which would be unstable by construction over the entire depth of the convection zone.

We should emphasize that the non-convective solution is mainly of academic interest. It is used to compute, for example, the Rayleigh number, a measure of the degree of instability. However, there is no doubt that the actual stratification of the Sun is close to marginal stability down to a fractional radius of 0.71 before turning decidedly stable, as confirmed by helioseismology (Christensen-Dalsgaard et al. 1991; Basu 1997). Nevertheless, the concept of convection being driven by surface cooling rather than heating from below has been promoted in a number of papers by Stein & Nordlund (1989, 1998) and elaborated upon by Spruit (1997), who introduced the idea that flows in the deeper parts are being driven nonlocally through the entropy rain from the surface. The question is to what extent this affects our understanding of the speed and especially the typical scales of convective motions, how MLT models would need to be modified, and whether this might have any bearing on the interpretation of the observed flow amplitudes at large scales.

To include the nonlocal effects described by Spruit (1997), we must look for a contribution to the enthalpy flux that is not related to the local entropy gradient.⁵ Such a term has been identified in the meteorological context by Deardorff (1966), who describes it as a counter-gradient flux. In Deardorff (1972), he derives an expression for this flux, which depends on the local temperature or entropy fluctuations which, in his case, come from measurements. In the present case, we assume that such fluctuations have their origin in what Spruit (1997, p. 406) refers to as threads, which are thin downdrafts on an almost perfectly isentropic passive fluid upflow between the threads. Spruit already emphasized back then the dynamical importance of the "forest of narrow cool threads that are produced at the solar surface." He further suggested that their absence in simulations with a fixed top boundary at some depth below the actual surface might be "the reason why they would produce large-scale flows with amplitudes that are about two orders of magnitude larger than observed."

In the scenario described by Spruit (1997), convection is driven solely within the surface layers. This is superficially reminiscent of convective overshoot, as modeled in some of the aforementioned papers (Hurlburt et al. 1986; Brandenburg et al. 1996), where convective flows are driven into stable layers. The extent of overshoot depends not only on the stiffness of the stable gradient that the convective plumes are flowing into (Hurlburt et al. 1994), but also on the flow speed. The depth of such a layer can therefore not be determined a priori from purely hydrostatic stability considerations. The lower boundary might then not be sharp. As explained above, in a hydrostatic non-convecting reference model, only a very thin surface layer is a priori unstable to convection. The rest is made unstable purely by bulk mixing. If entropy rain convection is a nonlocal phenomenon, the extent of convection should depend on surface properties and cannot be predicted from the local entropy gradient, similarly to convective overshoot. One might then not be able to understand the relatively sharp demarcation at the bottom of the convection zone, as found in global helioseismology. Of course, once convection has become fully developed, the low-entropy elements descend into buoyantly neutral layers, which is quite different from the usual overshoot. Furthermore, the usual overshoot layer is characterized by negative buoyancy and therefore a downward enthalpy flux (Hurlburt et al. 1986). An important purpose of the present paper is to produce a quantitative model and to demonstrate that, with a hypothetical nonlocal contribution to the flux, it is possible to obtain models that still have a sharply defined lower boundary.

A quantitative model of convection with entropy rain, even with the somewhat hypothetical Deardorff term included, would also be useful to illustrate the qualitative nature of the resulting stratification, and to show whether it is weakly super or subadiabatic. Indeed, as already argued by Spruit (1997), if the local mass fraction of entropy deficient material, descending from the cooling surface, decreases with depth and if the stratification outside the entropy rain were exactly isentropic, the resulting horizontally averaged entropy would increase with depth. This suggests that the entropy rain itself can make an otherwise vanishing radial entropy gradient negative and therefore Schwarzschild unstable, as seen in surface simulations. This raises two important questions. First, to what extent is such a stratification affected by radiative heating—especially toward the bottom of the convection zone, where it would tend to produce a positive (stable) mean entropy gradient in the upflows. Second, would such a negative (unstable) mean entropy gradient always lead to giant cell convection, as has been seen in global convection simulations (Miesch et al. 2008), or could the stratification still be stable to the excitation of large-scale flows?

We postpone addressing the two aforementioned issues connected with the qualitative reasoning of Spruit (1997) to the end of the paper (Section 6), and begin by highlighting the less commonly known fact that in the non-convecting reference model, only the layer in the top $1\,\mathrm{Mm}$ is convectively unstable (Section 2). We then explain the nature of the Deardorff flux (Section 3), present a correspondingly modified mixing length model (Section 4), and give some illustrative numerical solutions (Section 5). We discuss alternative explanations for the lack of giant cell convection in Section 6 and present our conclusions in Section 7.

The main argument for the existence of a highly unstable surface layer comes from the consideration of the associated non-convective reference solution, where the flux is forced to be transported by radiation only, i.e., $F={F}_{\mathrm{rad}}$. This is something that is not normally considered in stellar physics, because we know that such a solution would be unstable to convection and would therefore never be realized. However, as mentioned in the introduction, this solution is of certain academic interest. We postpone the discussion of an explicit numerical solution to Section 5 and present in this section the basic argument only at a qualitative level, making reference to earlier numerical calculations by Barekat & Brandenburg (2014) for a simple model. They considered an opacity law of the form

$$\begin{eqnarray}&&\kappa ={\kappa }_{0}{(\rho /{\rho }_{0})}^{a}{(T/{T}_{0})}^{b},\end{eqnarray} \tag{ 2 }$$

where a and b are adjustable parameters, ${\rho }_{0}$ and T₀ are reference values for density and temperature, respectively, and ${\kappa }_{0}$ gives the overall magnitude of the opacity. The essential point to note here is that the exponents a and b determine the gradient of specific entropy in a purely non-convecting reference model. In thermodynamic equilibrium, the radiative flux must be constant, i.e.,

$$\begin{eqnarray}&&{F}_{\mathrm{rad}}=-K\,{dT}/{dz}=\mathrm{const},\end{eqnarray} \tag{ 3 }$$

where $K=16{\sigma }_{\mathrm{SB}}{T}^{3}/(3\kappa \rho )$ is the radiative conductivity with ${\sigma }_{\mathrm{SB}}$ being the Stefan–Boltzmann constant, and z is the vertical coordinate in a Cartesian coordinate system. For the simple opacity law (2), but with

$$\begin{eqnarray}&&b\lt 4+a,\end{eqnarray} \tag{ 4 }$$

the optically thick regime is characterized by a constant temperature gradient and therefore $K=\mathrm{const}$ (see Appendix A). We have then a polytropic stratification with $\rho \propto {T}^{n}$, where

$$\begin{eqnarray}&&n=(3-b)/(1+a)\end{eqnarray} \tag{ 5 }$$

is the polytropic index. It is larger than −1 when Equation (4) is obeyed, i.e., when the pressure decreases with height, as is required for a physically meaningful solution. For a ratio of specific heats of $\gamma =5/3$, the value of n for marginal Schwarzschild stability is ${n}_{\mathrm{crit}}=1/(\gamma -1)=3/2$. In most of the solar convection zone the dominant opacity is the bound-free absorption owing to the absorption of light during ionization of a bound electron, which is well described by the Kramers-type opacity law with a = 1 and $b=-7/2$, and so n = 3.25, which corresponds to a Schwarzschild-stable solution. Only near the surface, at temperatures typically below 15,000 K, is the dominant opacity the ${{\rm{H}}}^{-}$ opacity. It can no longer be approximated by a simple power law of the type given by Equation (2). However, in limited density and temperature ranges, certain values of a and b can tentatively be specified, e.g., a = 0.5 and b = 7...18. Clearly, for all these values the constraint (4) is violated, so the hydrostatic stratification is no longer polytropic, but solutions can still be constructed numerically (Barekat & Brandenburg 2014) and they demonstrate, not surprisingly, that the stratification is highly unstable. What is more surprising is the fact that even with a combined opacity law of the form

$$\begin{eqnarray}&&{\kappa }^{-1}={\kappa }_{\mathrm{Kr}}^{-1}+{\kappa }_{{{\rm{H}}}^{-}}^{-1},\end{eqnarray} \tag{ 6 }$$

where ${\kappa }_{\mathrm{Kr}}$ and ${\kappa }_{{{\rm{H}}}^{-}}$ are given by Equation (2) with suitable exponents a and b, the solutions of the non-convective reference state is unstable only over a depth of approximately $1\,\mathrm{Mm}$. We return to this at the end of Section 5, where we present numerical solutions.

Of course, as stated earlier, the non-convective reference state is only of academic interest. Already with standard MLT (Biermann 1932; Vitense 1953), which allows for a non-vanishing enthalpy flux, one finds a vastly extended convection zone with a depth of the order of $100\,\mathrm{Mm}$ (Biermann 1938). However, the question now is how this can be affected by the presence of the Deardorff flux. This will be the subject of the rest of this paper. If the Deardorff flux were to become dominant and the stratification subadiabatic, it would become locally stable to the onset of convection. One might further speculate that the typical scale would no longer be controlled by the local pressure scale height, but it might be imprinted from the downdraft pattern just beneath the surface and be therefore comparable to the granulation scale or at least the supergranulation scale (Cossette & Rast 2016). This can have other important consequences that will also be addressed in this paper. It should be pointed out, however, that the surface simulations have so far not produced evidence for subadiabatic stratification. On the other hand, the models presented below predict subadiabatic stratification only a certain distance below the surface, depending on ill-known input parameters.

3.1. Derivation

In the meteorological context, counter-gradient heat flux terms have been noticed for a long time (Ertel 1942; Priestley & Swinbank 1947; Deardorff 1966). They appear naturally when calculating the enthalpy flux ${F}_{\mathrm{enth}}$ using the τ approximation in its minimalistic form (e.g., Blackman & Field 2003). In this approach, one computes the time derivative of ${F}_{\mathrm{enth}}$. In the absence of ionization effects, ${F}_{\mathrm{enth}}$ can be written as

$$\begin{eqnarray}&&{F}_{\mathrm{enth}}=\overline{\rho {u}_{z}{c}_{P}T},\end{eqnarray} \tag{ 7 }$$

where the overbar denotes horizontal averaging. In standard MLT, the enthalpy flux is usually referred to as the convective flux and the kinetic energy flux vanishes because of the assumed perfect symmetry between up- and downflows. In deeper layers especially, however, this is not justified (Stein et al. 2009), so this restriction will later be relaxed.

In the case of strongly stratified layers, it is convenient to use pressure P and specific entropy S as thermodynamic variables, because we later want to ignore pressure fluctuations on the grounds that pressure disturbances are quickly equilibrated by sound waves. Furthermore, we will restrict ourselves to second order correlations in Equation (7) and thus to fluctuations only in the correlation of specific entropy and velocity. Up to some reference value, we have

$$\begin{eqnarray}&&S/{c}_{P}=\mathrm{ln}T-{{\rm{\nabla }}}_{\mathrm{ad}}\mathrm{ln}P,\end{eqnarray} \tag{ 8 }$$

where ${{\rm{\nabla }}}_{\mathrm{ad}}=1-1/\gamma$ and $\gamma ={c}_{P}/{c}_{V}$ is the ratio of specific heats at constant pressure and constant volume, respectively, all constants for a perfect gas with a fixed degree of ionization. Ignoring pressure variations and linearizing $\delta \mathrm{ln}T=\delta T/T$, we can replace ${c}_{P}\delta T$ by $T\delta S$. In the following, we denote fluctuations by lower case characters, i.e., $S=\overline{S}+s$, where $s\equiv \delta S$ are used interchangeably. As argued above, other fluctuations are omitted. We also ignore mean flows ($\overline{{\boldsymbol{U}}}={\bf{0}}$), so there are only velocity fluctuations (${\boldsymbol{U}}={\boldsymbol{u}}$). Focusing thus on the dominant contribution proportional to $\overline{{u}_{z}s}$, we have

$$\begin{eqnarray}&&{F}_{\mathrm{enth}}=\overline{\rho }\overline{T}\,\overline{{u}_{z}s}.\end{eqnarray} \tag{ 9 }$$

Next, we write the time derivative of ${F}_{\mathrm{enth}}$ as

$$\begin{eqnarray}&&\partial {F}_{\mathrm{enth}}/\partial t=\overline{\rho }\overline{T}(\overline{{u}_{z}\dot{s}}+\overline{{\dot{u}}_{z}s}),\end{eqnarray} \tag{ 10 }$$

where dots denote partial time derivatives and changes of the background state have been neglected. Using the governing equations for u_i and s, we have (see Appendix B)

$$\begin{eqnarray}&&\dot{s}=-{u}_{j}{{\rm{\nabla }}}_{j}\overline{S}-s/{\tau }_{\mathrm{cool}}...,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\dot{u}}_{i}=-{g}_{i}s/{c}_{P}+...,\end{eqnarray} \tag{ 12 }$$

where g_i is the ith component of the gravitational acceleration in Cartesian coordinates, namely ${\boldsymbol{g}}=(0,0,-g)$, the ellipses refer to terms that are nonlinear in the fluctuations,

$$\begin{eqnarray}&&{\tau }_{\mathrm{cool}}^{-1}=\iota {c}_{\gamma }{k}_{{\rm{f}}}\ \ {\rm{with}}\ \ \iota ={\ell }{k}_{{\rm{f}}}/(3+{{\ell }}^{2}{k}_{{\rm{f}}}^{2})\end{eqnarray} \tag{ 13 }$$

is the inverse heating and cooling time owing to radiation (Unno & Spiegel 1966; Edwards 1990), ${c}_{\gamma }=16{\sigma }_{\mathrm{SB}}{\overline{T}}^{3}/\overline{\rho }{c}_{P}$ is the photon diffusion speed (Barekat & Brandenburg 2014), ${\ell }=1/\kappa \rho$ is the photon mean-free path, and ${k}_{{\rm{f}}}$ is the typical wavenumber of the fluctuations, which may be associated with the inverse mixing length used routinely in MLT. We mention at this point that the nonlocal nature of the entropy rain must lead to yet another contribution in Equation (11). We expect this to be the result of the nonlinear term (indicated by the ellipses), of which a part later gives rise to the ${{\rm{\nabla }}}_{{\rm{D}}}$ term. This will be motivated further at the end of this section and in Section 3.3, where it will be included in the final expression for ${{\rm{\nabla }}}_{{\rm{D}}}$.

Using Equations (11) and (12) in Equation (10), we have

$$\begin{eqnarray}&&\displaystyle \frac{\partial {F}_{i}^{\mathrm{enth}}}{\partial t}=\overline{\rho }\overline{T}(-\overline{{u}_{i}{u}_{j}}\,{{\rm{\nabla }}}_{j}\overline{S}-{g}_{i}\overline{{s}^{2}}/{c}_{P})-\displaystyle \frac{{F}_{i}^{\mathrm{enth}}}{{\tau }_{\mathrm{cool}}}+{{ \mathcal T }}_{i},\end{eqnarray} \tag{ 14 }$$

where the $\overline{{s}^{2}}$ term is the Deardorff flux and the ${{ \mathcal T }}_{i}$ refer to triple correlations that will be approximated by the quadratic correlation ${F}_{\mathrm{enth}}$ as

$$\begin{eqnarray}&&{{ \mathcal T }}_{i}=-{F}_{i}^{\mathrm{enth}}/\tau ,\end{eqnarray} \tag{ 15 }$$

where τ is a relaxation time due to turbulence (Blackman & Field 2003), which will later be identified with the turnover time. This procedure, in which correlations with pressure fluctuations are also neglected, is called the minimal τ approximation. Important aspects of the closure assumption (15) have been verified numerically for passive scalar transport (Brandenburg et al. 2004). Among other things, they found that the time derivative on the left-hand side of Equation (14) restores causality by turning the otherwise parabolic heat equation into a hyperbolic wave equation. Here, however, we are interested in slow variations such that the time derivative of ${F}_{\mathrm{enth}}$ can be neglected. (Note, however, that this does not imply that the cooling term will be neglected.)

Since the relaxation term in Equation (15) is similar to the cooling term in Equation (14), we can combine the two by introducing the reduced relaxation time ${\tau }_{\mathrm{red}}$ defined through

$$\begin{eqnarray}&&{\tau }_{\mathrm{red}}^{-1}={\tau }_{\mathrm{cool}}^{-1}+{\tau }^{-1}.\end{eqnarray} \tag{ 16 }$$

We can then solve for ${F}_{\mathrm{enth}}$, which appears on the right-hand sides (rhs) of Equations (14) and (15), and find ${{\boldsymbol{F}}}_{\mathrm{enth}}={{\boldsymbol{F}}}_{{\rm{G}}}+{{\boldsymbol{F}}}_{{\rm{D}}}$, where

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{\rm{G}}}=-\tfrac{1}{3}{\tau }_{\mathrm{red}}{u}_{\mathrm{rms}}^{2}\overline{\rho }\,\overline{T}\,{\boldsymbol{\nabla }}\overline{S},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{\rm{D}}}=-{\tau }_{\mathrm{red}}\overline{{s}^{2}}\,{\boldsymbol{g}}\,\overline{\rho }\,\overline{T}/{c}_{P}\end{eqnarray} \tag{ 18 }$$

are the ordinary gradient and the new Deardorff fluxes, respectively, and anisotropies have been ignored for the benefit of simpler notation. Thus, we write $\overline{{u}_{i}{u}_{j}}\approx \tfrac{1}{3}{\delta }_{{ij}}{u}_{\mathrm{rms}}^{2}$, where ${u}_{\mathrm{rms}}$ is the rms velocity of the turbulence. Equations (17) and (18) are written in vectorial forms to highlight the directions of the fluxes: ${{\boldsymbol{F}}}_{{\rm{G}}}$ is counter-gradient and ${{\boldsymbol{F}}}_{{\rm{D}}}$ is countergravity. In Equation (17), the term $\tfrac{1}{3}{\tau }_{\mathrm{red}}{u}_{\mathrm{rms}}^{2}\equiv {\chi }_{{\rm{t}}}$ is the turbulent thermal diffusivity.

In view of astrophysical applications, we replace the specific entropy gradient by the commonly defined superadiabatic gradient, i.e.,

$$\begin{eqnarray}&&-d(\overline{S}/{c}_{P})/{dz}=({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}})/{H}_{P},\end{eqnarray} \tag{ 19 }$$

where ${\rm{\nabla }}=d\mathrm{ln}\overline{T}/d\mathrm{ln}\overline{P}$ is the double-logarithmic temperature gradient and ${H}_{P}=-{(d\mathrm{ln}\overline{P}/{dz})}^{-1}$ is the pressure scale height. Thus, we arrive at

$$\begin{eqnarray}&&{F}_{\mathrm{enth}}=\tfrac{1}{3}\overline{\rho }{c}_{P}\overline{T}\,({\tau }_{\mathrm{red}}{u}_{\mathrm{rms}}^{2}/{H}_{P})({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\nabla }}}_{{\rm{D}}}),\end{eqnarray} \tag{ 20 }$$

where ${{\rm{\nabla }}}_{{\rm{D}}}$ is a new contribution to standard MLT, which results from ${F}_{{\rm{D}}}$. Using the terms on the rhs of Equation (18), we can write it explicitly as

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{{\rm{D}}}=(3/\gamma )(\overline{{s}^{2}}/{c}_{P}^{2})\ {{\rm{Ma}}}^{-2},\end{eqnarray} \tag{ 21 }$$

where ${\rm{Ma}}={u}_{\mathrm{rms}}/{c}_{{\rm{s}}}$ is the Mach number of the turbulence and ${c}_{{\rm{s}}}$ is the sound speed with ${c}_{{\rm{s}}}^{2}=\gamma {{gH}}_{P}$. Note that this Mach number dependence arises in order to cancel the corresponding ${u}_{\mathrm{rms}}^{2}$ factor in the definition (20). The Deardorff term is a contribution to the flux that is always directed outward and results from the transport of fluid elements with entropy fluctuations of either sign, as is illustrated in Figure 1.

Zoom In Zoom Out Reset image size

3.2. Physical Interpretation

3.3. Depth Dependence of the Deardorff Term

To estimate the resulting depth dependence of the Deardorff term in Equation (21), we must know how the $\overline{{s}^{2}}$ associated with the Deardorff term varies with depth. If the entropy rain was purely of the form of Hill vortex-like structures, as discussed above, their filling factor f_s would decrease with increasing depth and density like

$$\begin{eqnarray}&&{f}_{s}\propto {\overline{\rho }}^{-\zeta },\end{eqnarray} \tag{ 22 }$$

where $\zeta =0.8$ has been found for spherical vortex structures descending in an isothermally stratified layer; see Appendix C. However, Equation (23) is not contingent on Hill vortices, but is a consequence of downdrafts along a density gradient. For purely spherical compression one would expect $\zeta =2/3$, while for horizontal compression one has $\zeta =1$.

If we neglect non-ideal (radiative or viscous) effects, as well as entrainment between up- and downflows, the difference ${\rm{\Delta }}S={S}_{\uparrow }-{S}_{\downarrow }$ between the entropy of the slowly rising surroundings, ${S}_{\uparrow }$, and that in the downward propagating vortex, ${S}_{\downarrow }$, would remain constant and equal to the entropy deficit ${\rm{\Delta }}{S}_{0}$ suffered by overturning motions at the surface, which is the only location where radiative losses are significant. On sufficiently short length scales, however, radiative heating from the surroundings would erode this entropy difference and lead to a decrease of the effective ${\rm{\Delta }}S$. We model this by considering ζ an adjustable parameter that is increased relative to the value 0.8 that we expect in the absence of radiation, i.e., f_s decreases more strongly. The resulting fractional entropy difference, ${f}_{s}\,{\rm{\Delta }}{S}_{0}$, therefore decreases faster with depth than for $\zeta =0.8$.

As will become clear from the considerations below, if there is no entrainment from the upflows into the downdrafts and all of the downflows were confined to the small surface area of the vortex structures, the resulting downward-directed kinetic energy flux would increase with depth and eventually exceed the enthalpy flux. This would be unphysical, because convection should lead to outward energy transport. However, if there is entrainment with associated mixing, the actual filling factor (fractional area) of downflows, which we denote by f, would be larger and the mean entropy difference ${\rm{\Delta }}\overline{S}$ would be diluted correspondingly such that

$$\begin{eqnarray}&&f\,{\rm{\Delta }}\overline{S}={f}_{s}\,{\rm{\Delta }}{S}_{0}\quad (\propto {\overline{\rho }}^{-\zeta }\,{\rm{\Delta }}{S}_{0}).\end{eqnarray} \tag{ 23 }$$

We recall that non-ideal effects can be captured by choosing $\zeta \gt 0.8$ in the prescription (23). The physics of entrainment has been modeled in the stellar context by Rieutord & Zahn (1995) and the extent of entrainment has been quantified in realistic surface simulations by Trampedach & Stein (2011), who determined the entrainment length scale as the typical scale height of the mass flux in the up- and downflows separately. They call this scale the mass mixing length and find its value to be comparable to ${H}_{P}$.

To estimate the depth dependence of various horizontal averages, and to compute enthalpy and kinetic energy fluxes, we now compute various averages in the two-stream approximation, in which all horizontal averages are the sum of the fraction $1-f$ of the value in upflows and the fraction f of the value in downflows. Thus, for the mean specific entropy stratification we have

$$\begin{eqnarray}&&\overline{S}=(1-f){\overline{S}}_{\uparrow }+f{\overline{S}}_{\downarrow }={\overline{S}}_{\uparrow }-f\,{\rm{\Delta }}\overline{S},\end{eqnarray} \tag{ 24 }$$

where ${\overline{S}}_{\uparrow }$ and ${\overline{S}}_{\downarrow }$ are the mean specific entropies in up- and downflows, and ${\rm{\Delta }}\overline{S}={\overline{S}}_{\uparrow }-{\overline{S}}_{\downarrow }$ is their difference. We expect that ${\overline{S}}_{\uparrow }\approx {S}_{\uparrow }$ will be constant if there is no heating in the upwellings, while ${\overline{S}}_{\downarrow }\approx (1-{f}_{s}/f){S}_{\uparrow }+({f}_{s}/f){S}_{\downarrow }$ will be dominated by contributions from ${S}_{\uparrow }$ due to entrainment.

To calculate $\overline{{s}^{2}}$, we must first compute the fluctuating quantity $s=S-\overline{S}$ and then average the squared values in up- and downflows, which, using Equation (24), gives

$$\begin{eqnarray}&&\overline{{s}^{2}}=(1-f){({\overline{S}}_{\uparrow }-\overline{S})}^{2}+f{({\overline{S}}_{\downarrow }-\overline{S})}^{2}=\hat{f}{({\rm{\Delta }}\overline{S})}^{2},\end{eqnarray} \tag{ 25 }$$

where $\hat{f}=(1-f)f$ has been introduced as a shorthand. Analogously to the specific entropy, we find for the velocity

$$\begin{eqnarray}&&{\overline{U}}_{z}=(1-f){\overline{U}}_{\uparrow }+f{\overline{U}}_{\downarrow }={\overline{U}}_{\uparrow }-f\,{\rm{\Delta }}\overline{U},\end{eqnarray} \tag{ 26 }$$

where ${\overline{U}}_{\uparrow }\gt 0$ and ${\overline{U}}_{\downarrow }\lt 0$ are the mean up- and downflow velocities with ${\rm{\Delta }}\overline{U}={\overline{U}}_{\uparrow }-{\overline{U}}_{\downarrow }$. Since the densities in up- and downflows are nearly the same, especially in deeper layers, we have ${\overline{U}}_{z}\approx 0$ from mass conservation. With this, we find

$$\begin{eqnarray}&&{u}_{\mathrm{rms}}^{2}\equiv \overline{{{uz}}^{2}}=(1-f){\overline{U}}_{\uparrow }^{2}+f\,{\overline{U}}_{\downarrow }^{2}=\hat{f}{({\rm{\Delta }}\overline{U})}^{2},\end{eqnarray} \tag{ 27 }$$

and thus ${\overline{U}}_{\uparrow }/{u}_{\mathrm{rms}}={[f/(1-f)]}^{1/2}$. Similarly, we calculate

$$\begin{eqnarray}&&\overline{{u}_{z}^{3}}=(1-f){\overline{U}}_{\uparrow }^{3}+f\,{\overline{U}}_{\downarrow }^{3}=-\hat{f}(1-2f){({\rm{\Delta }}\overline{U})}^{3},\end{eqnarray} \tag{ 28 }$$

which is negative for $f\lt 1/2$, and finally

$$\begin{eqnarray}&&\overline{{u}_{z}s}=(1-f){\overline{U}}_{\uparrow }{\overline{S}}_{\uparrow }+f\,{\overline{U}}_{\downarrow }{\overline{S}}_{\downarrow }=\hat{f}\,{\rm{\Delta }}\overline{U}\,{\rm{\Delta }}\overline{S}.\end{eqnarray} \tag{ 29 }$$

Thus, the kinetic energy flux, $\overline{\rho {{\boldsymbol{u}}}^{2}{u}_{z}}/2$, which reduces to $\overline{\rho }\,\overline{{u}_{z}^{3}}/2$ in this two-stream approximation, is given by

$$\begin{eqnarray}&&{F}_{\mathrm{kin}}=-{\phi }_{\mathrm{kin}}\,\overline{\rho }{u}_{\mathrm{rms}}^{3},\end{eqnarray} \tag{ 30 }$$

where ${\phi }_{\mathrm{kin}}=(1/2-f)/{\hat{f}}^{1/2}$ is a positive prefactor (corresponding to downward kinetic energy flux) if $f\lt 1/2$. Stein et al. (2009) find $f\approx 1/3$, nearly independently of depth, which yields ${\phi }_{\mathrm{kin}}\approx \sqrt{2}/4\approx 0.35;$ see Table 1, where we list ${\phi }_{\mathrm{kin}}$ and $-{\overline{U}}_{\downarrow }/{u}_{\mathrm{rms}}={[(1-f)/f]}^{1/2}$ for selected values of f. The enthalpy flux is proportional to $\overline{{u}_{z}s}$ and, using Equations (9) and (29) together with Equations (25) and (27), we find ${F}_{\mathrm{enth}}=\overline{\rho }\overline{T}\,{u}_{\mathrm{rms}}{s}_{\mathrm{rms}}$. In Appendix D we show that $\overline{T}{s}_{\mathrm{rms}}\ ={u}_{\mathrm{rms}}^{2}{k}_{{\rm{f}}}{H}_{P}/({a}_{\mathrm{MLT}}{{\rm{\nabla }}}_{\mathrm{ad}})$, where ${a}_{\mathrm{MLT}}=1/8$ is a geometric factor in standard MLT.⁶ This leads to

$$\begin{eqnarray}&&{F}_{\mathrm{enth}}={\phi }_{\mathrm{enth}}\,\overline{\rho }{u}_{\mathrm{rms}}^{3}\end{eqnarray} \tag{ 31 }$$

with ${\phi }_{\mathrm{enth}}={k}_{{\rm{f}}}{H}_{P}/({a}_{\mathrm{MLT}}{{\rm{\nabla }}}_{\mathrm{ad}})$. This yields ${\phi }_{\mathrm{enth}}\approx 20$, which is rather large. By contrast, Brandenburg et al. (2005) determined a quantity k_u such that ${\phi }_{\mathrm{enth}}={k}_{u}^{-3/2}\approx 4$.

Table 1.  ${\phi }_{\mathrm{kin}}$, $-{\overline{U}}_{\downarrow }/{u}_{\mathrm{rms}}$, and ${\overline{U}}_{\uparrow }/{u}_{\mathrm{rms}}$ for Selected Values of f

f	1/2	1/3	0.14	0.015	0.0006
${\phi }_{\mathrm{kin}}$	0	0.35	1	4	20
$-{\overline{U}}_{\downarrow }/{u}_{\mathrm{rms}}$	1	1.4	2.5	8	40
${\overline{U}}_{\uparrow }/{u}_{\mathrm{rms}}$	1	0.7	0.4	0.12	0.025

Download table as:  ASCIITypeset image

We see that the presence of a kinetic energy flux just modifies the usual expression for the convective flux, which then becomes the sum of enthalpy and kinetic energy fluxes; see Section 3.1. This is compatible with recent simulations of R. F. Stein (2016, private communication), in which the fractional kinetic energy flux increases toward the deeper parts ($\gtrsim 40\,\mathrm{Mm}$) of the domain.

When f becomes small ($\lt 0.14$), ${\phi }_{\mathrm{kin}}$ exceeds unity and for $f\lt 0.015$, ${\phi }_{\mathrm{kin}}$ exceeds the estimate ${\phi }_{\mathrm{enth}}\approx 4$ found by Brandenburg et al. (2005), so the sum of enthalpy and kinetic energy fluxes may become negative, which appears unphysical. In that case, the idea of reconciling the results of Hanasoge et al. (2012) with such convection transporting the solar luminosity could be problematic, unless both up- and downflows were to occur on sufficiently small scales. In the case f = 0.14, the resulting ${\overline{U}}_{\downarrow }$ would still not be particularly fast and only 2.5 times larger than ${u}_{\mathrm{rms}};$ for f = 0.015 we have $-{\overline{U}}_{\downarrow }\approx 8{u}_{\mathrm{rms}};$ see Table 1. Conversely, if the cold entropy blobs were much smaller and were still to contribute significantly to the total energy flux, i.e., if they were much faster, this might not be compatible with an upward directed total energy transport. A possible alternative might be suggested by the compressible simulations of Cattaneo et al. (1991), where the downward-directed kinetic energy flux in the downdrafts was found to balance the upward directed enthalpy flux in the downdrafts. This would imply that convection would transport energy only in the upwellings. This result has however not been confirmed in realistic surface simulations (Stein et al. 1992).

Equation (23) shows how f_s decreases with depth, but its actual value and that of ${{\rm{\nabla }}}_{{\rm{D}}}$ remains undetermined. In the following, we propose a quantitative prescription for ${{\rm{\nabla }}}_{{\rm{D}}}$ by estimating its value within the top few hundred kilometers. In those top layers, we expect ${{\rm{\nabla }}}_{{\rm{D}}}$ to reach its maximum value, ${{\rm{\nabla }}}_{{\rm{D}}}^{\max }$, and it should be a certain fraction of ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$. In Appendix E we show that ${{\rm{\nabla }}}_{{\rm{D}}}^{\max }/{({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}})}_{\max }\approx 3$.

In deeper layers, where the local value of ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$ has become small, ${{\rm{\nabla }}}_{{\rm{D}}}$ should scale with $\overline{{s}^{2}}=\hat{f}{({\rm{\Delta }}\overline{S})}^{2}$, and therefore, using Equation (23), it should be proportional to ${f}_{s}^{2}(z)\propto {\overline{\rho }}^{-2\zeta }$, in addition to an ${{\rm{Ma}}}^{-2}$ factor; see Equation (21). Thus, we have

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{{\rm{D}}}={f}_{s}^{2}\,{{\rm{\nabla }}}_{{\rm{D}}}^{\max },\end{eqnarray} \tag{ 32 }$$

where we have used for f_s the expression

$$\begin{eqnarray}&&{f}_{s}={f}_{s0}{(\overline{\rho }/{\rho }_{* })}^{-\tilde{\zeta }}\quad ({\rm{for}}\,{\rm{}}\overline{\rho }\gt {\rho }_{* }).\end{eqnarray} \tag{ 33 }$$

Here, f_s0 is a prefactor determining the strength of the resulting Deardorff flux and ${\rho }_{* }$ is the density at the point at which ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$ is equal to its maximum value (which is just below the photosphere). The new exponent $\tilde{\zeta }=\zeta -{\rm{\Delta }}\zeta$ takes the scaling of Mach number with density into account, i.e.,

$$\begin{eqnarray}&&{\rm{Ma}}\propto {\overline{\rho }}^{-{\rm{\Delta }}\zeta },\end{eqnarray} \tag{ 34 }$$

where ${\rm{\Delta }}\zeta$ ($\gt 0$) will be computed in Section 4.3.

3.4. Kinetic Energy Flux

The importance of the kinetic energy flux has been stressed for some time as a property of compressible stratified convection (Hurlburt et al. 1984; Cattaneo et al. 1991). This flux is related to the asymmetry of up- and downflows (Stein & Nordlund 1989) and is non-vanishing when $f\ne 1/2;$ see Section 3.3. It is neglected in standard MLT, as has been discussed by Arnett et al. (2015). Just like the Deardorff flux, the kinetic energy flux is also a non-gradient flux, but it is always directed downward for $f\lt 1/2$ and must therefore be overcome by the enthalpy flux so that energy can still be transported outwards. However, since the lowest order correlations in ${F}_{\mathrm{kin}}$ are triple correlations, there is no τ approximation treatment analogous to that of the Deardorff flux. However, by comparing ${F}_{\mathrm{kin}}=-{\phi }_{\mathrm{kin}}\,\overline{\rho }{u}_{\mathrm{rms}}^{3}$ with the enthalpy flux in Equation (20), it is possible to define a corresponding nabla term via ${F}_{\mathrm{kin}}=-\tfrac{1}{3}\overline{\rho }{c}_{P}\overline{T}\,({\tau }_{\mathrm{red}}{u}_{\mathrm{rms}}^{2}/{H}_{P}){{\rm{\nabla }}}_{\mathrm{kin}}$. This yields

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mathrm{kin}}=3\,{\phi }_{\mathrm{kin}}\,{u}_{\mathrm{rms}}{H}_{P}/{\tau }_{\mathrm{red}}{c}_{P}\overline{T},\end{eqnarray} \tag{ 35 }$$

which is obtained analogously to ${{\rm{\nabla }}}_{{\rm{D}}}$ in Equation (21). The prefactor is here defined with a positive sign, so the total non-radiative flux caused by the turbulence is proportional to ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\nabla }}}_{{\rm{D}}}-{{\rm{\nabla }}}_{\mathrm{kin}}$.

Obviously, if ${{\rm{\nabla }}}_{\mathrm{kin}}$ is strictly proportional to ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$, the addition of the ${{\rm{\nabla }}}_{\mathrm{kin}}$ term does not modify standard MLT, provided that ${{\rm{\nabla }}}_{\mathrm{kin}}$ depends just on the local value of ${u}_{\mathrm{rms}}$, which in turn is, again, assumed to depend just on the local value of ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$. This is different for ${{\rm{\nabla }}}_{{\rm{D}}}$, which depends on the entropy deficit produced by cooling near the surface and in this way on the value of ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$ at the position where $\overline{\rho }={\rho }_{* }$. The ${{\rm{\nabla }}}_{{\rm{D}}}$ term is therefore truly nonlocal. In the following, we combine kinetic and enthalpy fluxes into a total contribution ${F}_{\mathrm{conv}}={F}_{\mathrm{enth}}+{F}_{\mathrm{kin}}$ which arises from convection.

4.1. Stratification and Flux Balance

To construct an equilibrium model, we begin by considering first the case without convection, so the flux F is carried by radiation alone. Hydrostatic and thermal equilibrium then imply $d\overline{P}/{dz}=-\overline{\rho }g$ and ${Kd}\overline{T}/{dz}=-F$, or, alternatively for the logarithmic gradients,

$$\begin{eqnarray}&&d\mathrm{ln}\overline{P}/{dz}=-\overline{\rho }g/\overline{P},\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&d\mathrm{ln}\overline{T}/{dz}=-F/(K\overline{T}).\end{eqnarray} \tag{ 37 }$$

The double-logarithmic temperature gradient is obtained by dividing the two equations through each other, i.e.,

$$\begin{eqnarray}&&{\rm{\nabla }}=\displaystyle \frac{d\mathrm{ln}\overline{T}}{d\mathrm{ln}\overline{P}}=\displaystyle \frac{F\overline{P}}{K\overline{T}\,\overline{\rho }g}=\displaystyle \frac{{{Fc}}_{P}{{\rm{\nabla }}}_{\mathrm{ad}}}{{Kg}},\end{eqnarray} \tag{ 38 }$$

where we have used the perfect gas equation of state in the form $\overline{P}/\overline{T}\,\overline{\rho }={c}_{P}-{c}_{V}={c}_{P}(1-1/\gamma )={c}_{P}{{\rm{\nabla }}}_{\mathrm{ad}}$. If the energy is no longer carried by radiation alone, ${\rm{\nabla }}$ cannot be computed from Equation (38), but we have to invoke a suitable theory of convection. In standard MLT, one obtains ${\rm{\nabla }}$ as a solution of a cubic equation (Vitense 1953; Kippenhahn & Weigert 1990). In the following, we consider a modification that accounts for the possibility of a Deardorff flux.

Flux balance implies that the sum of the radiative, enthalpy, and kinetic energy fluxes equals the total flux, i.e.,

$$\begin{eqnarray}&&{F}_{\mathrm{tot}}={F}_{\mathrm{rad}}+{F}_{\mathrm{conv}}.\end{eqnarray} \tag{ 39 }$$

where ${F}_{\mathrm{conv}}={F}_{\mathrm{enth}}+{F}_{\mathrm{kin}}$ is the flux that arises from convection. In MLT, it is customary to express these fluxes in terms of nablas. For a given double-logarithmic temperature gradient ${\rm{\nabla }}$, the radiative flux is evidently

$$\begin{eqnarray}&&{F}_{\mathrm{rad}}=\displaystyle \frac{{Kg}}{{c}_{P}{{\rm{\nabla }}}_{\mathrm{ad}}}\,{\rm{\nabla }},\end{eqnarray} \tag{ 40 }$$

where ${\rm{\nabla }}$ characterizes the actual temperature gradient. We can also define a hypothetical radiative temperature gradient ${{\rm{\nabla }}}_{\mathrm{rad}}$ that would result if all the energy were carried by radiation, so we can write

$$\begin{eqnarray}&&{F}_{\mathrm{tot}}=\displaystyle \frac{{Kg}}{{c}_{P}{{\rm{\nabla }}}_{\mathrm{ad}}}\,{{\rm{\nabla }}}_{\mathrm{rad}},\end{eqnarray} \tag{ 41 }$$

which follows from Equation (38). We note in passing that the ratio between ${F}_{\mathrm{tot}}$ and the radiative flux carried by the adiabatic temperature gradient, ${Kg}/{c}_{P}$, is known in laboratory and theoretical studies of convection as the Nusselt number (Hurlburt et al. 1984), whose local value is thus equal to

$$\begin{eqnarray}&&{\rm{Nu}}={{\rm{\nabla }}}_{\mathrm{rad}}/{{\rm{\nabla }}}_{\mathrm{ad}}.\end{eqnarray} \tag{ 42 }$$

Finally, as explained in Section 3, we have from Equation (20)

$$\begin{eqnarray}&&{F}_{\mathrm{enth}}={F}_{0}({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\nabla }}}_{{\rm{D}}}),\end{eqnarray} \tag{ 43 }$$

with ${F}_{0}=\tfrac{1}{3}\overline{\rho }{c}_{P}\overline{T}{\tau }_{\mathrm{red}}{u}_{\mathrm{rms}}^{2}/{H}_{P}$, and ${F}_{\mathrm{kin}}$ being essentially proportional to ${F}_{\mathrm{enth}};$ see Equations (30) and (31). Flux equilibrium then implies

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mathrm{rad}}={\rm{\nabla }}+\epsilon ({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\nabla }}}_{{\rm{D}}}),\end{eqnarray} \tag{ 44 }$$

where $\epsilon ={F}_{0}{c}_{P}{{\rm{\nabla }}}_{\mathrm{ad}}(1-{\phi }_{\mathrm{kin}}/{\phi }_{\mathrm{enth}})/{Kg}$. However, the expression for F₀ involves the still unknown values of ${u}_{\mathrm{rms}}$ and ${\tau }_{\mathrm{red}}$, which will be discussed in Section 4.2.

4.2. Relaxation Time and Mixing Length

In convection the relevant timescale is the turnover time, which we write as $\tau =1/{u}_{\mathrm{rms}}{k}_{{\rm{f}}}$. We argue that ${k}_{{\rm{f}}}$ should be estimated via the separation between the entropy rain structures and not their thickness. This becomes important in the picture in which the downdraft threads of the entropy rain merge with neighboring ones to form a tree-like structure, as seen in the surface simulations (Stein & Nordlund 1989; Spruit 1997), leading therefore to different scalings of the two length scales (separation and thickness of structures) with depth. It is then not obvious which of these scales are more relevant for determining the ${k}_{{\rm{f}}}$ that is relevant for mixing. In view of the uncertainty regarding the choice of ${k}_{{\rm{f}}}$, as well as for comparison with standard MLT, we consider models with a fixed value of ${k}_{{\rm{f}}}$ as well as the more conventional case in which ${k}_{{\rm{f}}}{H}_{P}\approx \mathrm{const}$. To capture the various cases in one expression, we assume in the following

$$\begin{eqnarray}&&{k}_{{\rm{f}}}={\alpha }_{\mathrm{mix}}^{\beta }{H}_{P}^{-1}{({k}_{{\rm{f}}0}{H}_{P})}^{1-\beta },\end{eqnarray} \tag{ 45 }$$

where $\beta =0$ corresponds to ${k}_{{\rm{f}}}={k}_{{\rm{f}}0}$ with a fixed value ${k}_{{\rm{f}}0}$ of the wavenumber of the entropy rain and $\beta =1$ corresponds to ${k}_{{\rm{f}}}{H}_{P}={\alpha }_{\mathrm{mix}}$, where we have allowed for the possibility of a free mixing length parameter, ${\alpha }_{\mathrm{mix}}$, as is commonly done in standard MLT. It is not to be confused with the parameter ${a}_{\mathrm{MLT}}$ that was introduced in Section 3.3 and Appendix D. Negative values of β would correspond to shrinking scales, but will not be considered here. Likewise, non-integer values of β are conceivable, but will also not be considered here.

Returning now to the discussion of the relaxation time τ in the beginning of this subsection, instead of associating it with the turnover time ${({u}_{\mathrm{rms}}{k}_{{\rm{f}}})}^{-1}$, we will allow for the possibility of an additional dilution factor $\phi (z)$ and write $\tau ={({u}_{\mathrm{rms}}{k}_{{\rm{f}}}\phi )}^{-1}$. Here, $\phi (z)$ increases with depth in a similar fashion as the scale height, so we assume, in analogy with our treatment of ${k}_{{\rm{f}}}$ in Equation (45), the expression $\phi ={({k}_{{\rm{f}}0}{H}_{P}/{\alpha }_{\mathrm{mix}})}^{\beta -\tilde{\beta }}$. This dilution factor only enters in expressions involving the turbulent diffusivity, such as in Equation (17), and hence terms involving τ or ${\tau }_{\mathrm{red}}$. For $\beta =1$, the case $\tilde{\beta }=1$ corresponds to the usual MLT concept, while $\tilde{\beta }=0$ corresponds to a value of ϕ that increases with depth.

In Figure 2 we present illustrative flow structures as well as their depth-dependent horizontal power spectra associated with the three combinations of β and $\tilde{\beta }$ considered here. At the top of the domain, the size and separation of flow structures are the same in all three cases. They remain constant with depth in the case $\beta =\tilde{\beta }=0$ (Case I), while for $\beta =0$ and $\tilde{\beta }=1$ (Case II), the separation increases with depth, but the thickness is still constant. Finally, for $\beta =\tilde{\beta }=1$ (Case III), both thickness and separation increase with depth. The filling factor decreases with depth in the case $\beta =0$ and $\tilde{\beta }=1$, which could potentially make the kinetic energy flux divergent with depth, while for both $\beta =\tilde{\beta }=0$ and $\beta =\tilde{\beta }=1$, the filling factor is independent of height.

Zoom In Zoom Out Reset image size

With these preparations in place, we can write F₀ in Equation (43) in the form

$$\begin{eqnarray}&&{F}_{0}=(\sigma /3\gamma {{\rm{\nabla }}}_{\mathrm{ad}})\,\overline{\rho }{u}_{\mathrm{rms}}{c}_{{\rm{s}}}^{2}{\alpha }_{\mathrm{mix}}^{-\tilde{\beta }}{({k}_{{\rm{f}}0}{H}_{P})}^{-(1-\tilde{\beta })},\end{eqnarray} \tag{ 46 }$$

where

$$\begin{eqnarray}&&\sigma \equiv {\tau }_{\mathrm{red}}{u}_{\mathrm{rms}}{k}_{{\rm{f}}}\phi ={u}_{\mathrm{rms}}/({u}_{\mathrm{rms}}+\iota {c}_{\gamma }/\phi )\end{eqnarray} \tag{ 47 }$$

quantifies the radiative heat exchange between convective elements and the surroundings. The ι term defined in Equation (13) has a maximum at ${\ell }{k}_{{\rm{f}}}=\sqrt{3}$, which typically occurs near the surface (e.g., Barekat & Brandenburg 2014).

In Equation (46), the local value of the turbulent rms velocity always depends on the actual flux transported, and therefore it must also depend on ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\nabla }}}_{{\rm{D}}}$. On dimensional grounds, since ${F}_{\mathrm{enth}}$ is proportional to ${u}_{\mathrm{rms}}^{3}$, we have

$$\begin{eqnarray}&&{u}_{\mathrm{rms}}={c}_{0}{({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\nabla }}}_{{\rm{D}}})}^{1/2},\end{eqnarray} \tag{ 48 }$$

so ${F}_{\mathrm{enth}}$ and ${F}_{\mathrm{kin}}$ are proportional to ${({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\nabla }}}_{{\rm{D}}})}^{3/2}$. Standard mixing length arguments can be used to show that the prefactor in Equation (48) is a fraction of ${c}_{{\rm{s}}}$. As shown in Appendix F, we have

$$\begin{eqnarray}&&{c}_{0}/{c}_{{\rm{s}}}=\sqrt{\sigma {a}_{\mathrm{MLT}}/3\gamma }\,{\alpha }_{\mathrm{mix}}^{-\beta ^{\prime} }{({k}_{{\rm{f}}0}{H}_{P})}^{-(1-\beta ^{\prime} )},\end{eqnarray} \tag{ 49 }$$

where $\beta ^{\prime} =(\beta +\tilde{\beta })/2$. It is then clear that ${F}_{\mathrm{kin}}$ scales with ${H}_{P}$ like ${({k}_{{\rm{f}}0}{H}_{P})}^{-3(1-\beta ^{\prime} )}$. This is not the case for ${F}_{\mathrm{enth}}$, however, because of the factor ${({k}_{{\rm{f}}0}{H}_{P})}^{1-\tilde{\beta }}$ which, together with the ${({k}_{{\rm{f}}0}{H}_{P})}^{1-\beta ^{\prime} }$ factor, gives the scaling proportional to ${({k}_{{\rm{f}}0}{H}_{P})}^{2(1-\beta ^{\prime\prime} )}$, where $\beta ^{\prime\prime} =(\beta ^{\prime} +\tilde{\beta })/2$.

4.3. Equation for the Superadiabatic Gradient

Now that we know F₀, we can solve the equation for the superadiabatic gradient. This leads to an equation that is similar to the cubic equation for ${\rm{\nabla }}$, which is familiar from standard MLT (Kippenhahn & Weigert 1990),

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mathrm{rad}}={\rm{\nabla }}+{\epsilon }_{* }{({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\nabla }}}_{{\rm{D}}})}^{\xi },\end{eqnarray} \tag{ 50 }$$

where $\xi =3/2$ and ${\epsilon }_{* }=(1-{\phi }_{\mathrm{kin}}/{\phi }_{\mathrm{enth}})({\epsilon }_{\mathrm{enth}}+{\epsilon }_{\mathrm{kin}})$ has contributions from the enthalpy and kinetic energy fluxes. These expressions are similar to in Equation (44), except that ${\epsilon }_{\mathrm{enth}}$ is evaluated with c₀ in place of ${u}_{\mathrm{rms}}$, i.e.,

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{enth}}={(\sigma /3\gamma )}^{3/2}\,{a}_{\mathrm{MLT}}^{1/2}\ ({c}_{{\rm{s}}}^{3}/\chi g){\alpha }_{\mathrm{mix}}^{-2\beta ^{\prime\prime} }{({k}_{{\rm{f}}0}{H}_{P})}^{-2(1-\beta ^{\prime\prime} )},\end{eqnarray} \tag{ 51 }$$

where $\chi =K/\overline{\rho }{c}_{P}$ is the radiative diffusivity. This shows that ${\epsilon }_{* }$ is essentially a Péclet number based on ${c}_{{\rm{s}}}$. The contribution from the kinetic energy flux is given by

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{kin}}/{\epsilon }_{\mathrm{enth}}=-{\phi }_{\mathrm{kin}}{a}_{\mathrm{MLT}}{{\rm{\nabla }}}_{\mathrm{ad}}\,{\alpha }_{\mathrm{mix}}^{-\beta }{({k}_{{\rm{f}}0}{H}_{P})}^{-(1-\beta )}.\end{eqnarray} \tag{ 52 }$$

We note in passing that ${\epsilon }_{* }$ is also related to the Rayleigh number, which is commonly defined in laboratory and numerical studies of convection; see Appendix G. Furthermore, because of convection and the resulting bulk mixing, $\overline{S}$ is now approximately constant, and therefore, unlike in the non-convecting reference solution with $K=\mathrm{const}$ (Section 2), K can no longer be constant, but it reaches a minimum at the point where κ is maximum, which turns out to be at a depth of about $1\,\mathrm{Mm}$ in the convection models presented below. Since ${\epsilon }_{* }$ is inversely proportional to K, it reaches a maximum at that depth and falls off both toward the top and the bottom of the convection zone.

An essential difference between Equation (50) and the usual one in MLT is the presence of ${{\rm{\nabla }}}_{{\rm{D}}}$ arising from the Deardorff flux. Within the usual MLT, where ${{\rm{\nabla }}}_{{\rm{D}}}=0$, one finds that ${\rm{\nabla }}$ is slightly above ${{\rm{\nabla }}}_{\mathrm{ad}}$, but now it might instead be slightly above ${{\rm{\nabla }}}_{\mathrm{ad}}-{{\rm{\nabla }}}_{{\rm{D}}}$. There are indeed two possibilities for convecting solutions (${F}_{\mathrm{conv}}\gt 0$), one corresponding to a Schwarzschild-stable solution,

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mathrm{ad}}-{{\rm{\nabla }}}_{{\rm{D}}}\lt {\rm{\nabla }}\lt {{\rm{\nabla }}}_{\mathrm{ad}}\quad {\rm{(stable)}},\end{eqnarray} \tag{ 53 }$$

and one that is Schwarzschild unstable,

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mathrm{ad}}-{{\rm{\nabla }}}_{{\rm{D}}}\lt {{\rm{\nabla }}}_{\mathrm{ad}}\lt {\rm{\nabla }}\quad {\rm{(unstable)}}.\end{eqnarray} \tag{ 54 }$$

Which of the two possibilities is attained depends on the value of ${{\rm{\nabla }}}_{{\rm{D}}}$ and also on details of the solution. As will be discussed in Section 6 below, entropy rain convection may actually still be Schwarzschild unstable without exciting giant cell convection if small-scale turbulent viscosity and diffusivity are strong enough so that the local turbulent Rayleigh number for the deeper layers is subcritical.

In this connection, we note that in standard MLT, one includes the effects of radiative cooling of the convective elements in a different manner than here. Instead of ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$, the effective buoyancy force is written as ${\rm{\nabla }}-{\rm{\nabla }}^{\prime}$, where ${\rm{\nabla }}^{\prime}$ always lies between ${\rm{\nabla }}$ and ${{\rm{\nabla }}}_{\mathrm{ad}}$ (Vitense 1953). Thus, one has ${{\rm{\nabla }}}_{\mathrm{ad}}\lt {\rm{\nabla }}^{\prime} \lt {\rm{\nabla }}$, which resembles Equation (54) with a negative value of ${{\rm{\nabla }}}_{{\rm{D}}}$.

To understand the nature of the solutions of Equation (50), it is instructive to treat ξ as an adjustable parameter. For given values of ${{\rm{\nabla }}}_{\mathrm{rad}}$ and ${\widetilde{{\rm{\nabla }}}}_{\mathrm{ad}}\equiv {{\rm{\nabla }}}_{\mathrm{ad}}-{{\rm{\nabla }}}_{{\rm{D}}}$, the case $\xi =1$ yields

$$\begin{eqnarray}&&{\rm{\nabla }}({\epsilon }_{* })=\displaystyle \frac{{{\rm{\nabla }}}_{\mathrm{rad}}+{\epsilon }_{* }{\widetilde{{\rm{\nabla }}}}_{\mathrm{ad}}}{1+{\epsilon }_{* }}.\end{eqnarray} \tag{ 55 }$$

It shows that ${\rm{\nabla }}\to {{\rm{\nabla }}}_{\mathrm{rad}}$ for ${\epsilon }_{* }\to 0$ (stable surface layers) and ${\rm{\nabla }}\to {\widetilde{{\rm{\nabla }}}}_{\mathrm{ad}}$ for ${\epsilon }_{* }\gg 1$ (deeper layers). Next, to discuss the general case $\xi \ne 1$, we define ${\rm{\Delta }}{\rm{\nabla }}={\rm{\nabla }}-{\widetilde{{\rm{\nabla }}}}_{\mathrm{ad}}$ and ${\rm{\Delta }}{{\rm{\nabla }}}_{\mathrm{rad}}\ ={{\rm{\nabla }}}_{\mathrm{rad}}-{\widetilde{{\rm{\nabla }}}}_{\mathrm{ad}}$. For ${\rm{\Delta }}{{\rm{\nabla }}}_{\mathrm{rad}}\lt 0$ we have ${\rm{\Delta }}{\rm{\nabla }}={\rm{\Delta }}{{\rm{\nabla }}}_{\mathrm{rad}}$, while for ${\rm{\Delta }}{{\rm{\nabla }}}_{\mathrm{rad}}\gt 0$ and ${\rm{\Delta }}{\rm{\nabla }}\ll 1$, a useful approximation is

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\nabla }}\approx {\rm{\Delta }}{{\rm{\nabla }}}_{\mathrm{rad}}^{1/\xi }/(q{\rm{\Delta }}{{\rm{\nabla }}}_{\mathrm{rad}}^{1/\xi -1}+{\epsilon }_{* }^{1/\xi }),\end{eqnarray} \tag{ 56 }$$

where $q={\xi }^{-1}$. It agrees with Equation (55) in the special case $\xi =1$, where ${\rm{\Delta }}{\rm{\nabla }}={\rm{\Delta }}{{\rm{\nabla }}}_{\mathrm{rad}}/(1+{\epsilon }_{* })$. In Figure 3 we plot ${\rm{\nabla }}$ versus ${\epsilon }_{* }$ for ${{\rm{\nabla }}}_{\mathrm{rad}}={10}^{5}$. The approximation yields ${\rm{\nabla }}\gt {{\rm{\nabla }}}_{\mathrm{rad}}$ for ${\epsilon }_{* }\lt {10}^{-3}$, which is unphysical. This can be mitigated by choosing q = 1; see Figure 3.

Zoom In Zoom Out Reset image size

For the relevant case of large values of ${\epsilon }_{* }$, we have

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\nabla }}\approx {({{\rm{\nabla }}}_{\mathrm{rad}}/{\epsilon }_{* })}^{1/\xi }.\end{eqnarray} \tag{ 57 }$$

This relation is useful because, even though both ${{\rm{\nabla }}}_{\mathrm{rad}}$ and ${\epsilon }_{* }$ depend on K, their ratio does not and is given by

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\nabla }}}_{\mathrm{rad}}}{{\epsilon }_{* }}=\displaystyle \frac{3{{\rm{\nabla }}}_{\mathrm{ad}}}{\sigma }\,\displaystyle \frac{{F}_{\mathrm{tot}}/\overline{\rho }}{{c}_{0}g/{k}_{{\rm{f}}}}\propto \displaystyle \frac{{({k}_{{\rm{f}}0}{H}_{P})}^{2-2\beta ^{\prime\prime} }}{\overline{\rho }{c}_{{\rm{s}}}^{3}}\propto {\overline{T}}^{-(1+2\beta ^{\prime\prime} )},\end{eqnarray} \tag{ 58 }$$

where $\beta ^{\prime\prime} =(\beta +3\tilde{\beta })/4$. Thus, ${\rm{\Delta }}{\rm{\nabla }}\propto {\overline{T}}^{-(1+2\beta ^{\prime\prime} )/\xi }$ gives the scaling of ${\rm{\Delta }}{\rm{\nabla }}$ for the bulk of the convection zone. Therefore, looking at Equation (48), we find ${\rm{Ma}}\propto {\overline{T}}^{-m}$ with

$$\begin{eqnarray}&&m=1-\beta ^{\prime} +(1/2+\beta ^{\prime\prime} )/\xi .\end{eqnarray} \tag{ 59 }$$

For $\xi =3/2$, using the relation $3\beta ^{\prime} -2\beta ^{\prime\prime} =\beta$, we find that $m=(4-\beta )/3$ is independent of the value of $\tilde{\beta }$. Thus, we have $m=4/3$ with $\beta =0$ and m = 1 with $\beta =1$. For isentropic stratification, this implies for the Mach number, given by Equation (34), the following scaling: ${\rm{\Delta }}\zeta =8/9$ with $\beta =0$ and ${\rm{\Delta }}\zeta =2/3$ with $\beta =1$.

4.4. Relative Importance of the Deardorff Term

In Section 3.3 we have considered the depth dependence of the Deardorff term via Equations (32) and (33). However, for ${{\rm{\nabla }}}_{{\rm{D}}}$ to be important at increasing depths, it must exceed the subadiabatic gradient ${{\rm{\nabla }}}_{\mathrm{ad}}-{\rm{\nabla }}$, because otherwise it would not be possible for the Deardorff term to make ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\nabla }}}_{{\rm{D}}}$ positive. From Equations (57) and (58) we see that ${\rm{\Delta }}{\rm{\nabla }}$ depends on $\overline{T}$ in a power-law fashion. We are particularly interested in the conditions under which ${\rm{\Delta }}{\rm{\nabla }}$ falls off faster than ${{\rm{\nabla }}}_{{\rm{D}}}$, because that would ensure that the ${{\rm{\nabla }}}_{{\rm{D}}}$ term remains important even at larger depths.

For $\xi =3/2$, and since the convection zone is nearly isentropically stratified ($\overline{T}\propto {\overline{\rho }}^{2/3}$), we have

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\nabla }}\propto {\overline{\rho }}^{-{\zeta }_{{\rm{\Delta }}{\rm{\nabla }}}}\quad {\rm{with}}\quad {\zeta }_{{\rm{\Delta }}{\rm{\nabla }}}=\displaystyle \frac{4}{9}(1+2\beta ^{\prime\prime} ).\end{eqnarray} \tag{ 60 }$$

In Table 2 we compare for various combinations of β and $\tilde{\beta }$ the exponents for ${\rm{\Delta }}{\rm{\nabla }}$ and ${{\rm{\nabla }}}_{{\rm{D}}}$, i.e., $2\tilde{\zeta }=2\zeta$ and $2{\rm{\Delta }}\zeta$, where ${\rm{\Delta }}\zeta =2\,m/3$. In all cases we have chosen $\zeta =1$, i.e., we allow for moderately non-ideal (radiative) effects relative to the ideal case with $\zeta =0.8$. We see that the difference is positive and non-vanishing in all cases, except for $\beta =1$ and $\tilde{\beta }=0$. In the following, we present solutions for all of the remaining three cases. Before doing this, let us recapitulate what led to the threefold dominance of $\tilde{\beta }$ over β. For better illustration, we summarize in Table 3 the various relationships that led to the scaling of ${\rm{\Delta }}{\rm{\nabla }}$ with ${k}_{{\rm{f}}0}{H}_{P}$.

Table 2.  Comparison of ${\zeta }_{{\rm{\Delta }}{\rm{\nabla }}}$ and $2\tilde{\zeta }=2\zeta -2{\rm{\Delta }}\zeta$, and Their Difference for Various Combinations of β and $\tilde{\beta }$ Using $\zeta =1$

	β	$\tilde{\beta }$	$\beta ^{\prime}$	$\beta ^{\prime\prime}$	${\zeta }_{{\rm{\Delta }}{\rm{\nabla }}}$	$2\tilde{\zeta }$	${\zeta }_{{\rm{\Delta }}{\rm{\nabla }}}-2\tilde{\zeta }$
Case I	0	0	0	0	4/9	2/9	2/9
Case II	0	1	1/2	3/4	10/9	2/9	8/9
	1	0	1/2	1/4	2/3	2/3	0
Case III	1	1	1	1	4/3	2/3	2/3

Download table as:  ASCIITypeset image

Table 3.  Illustration of Scaling Relationships with ${k}_{{\rm{f}}0}{H}_{P}$

${F}_{\mathrm{enth}}\propto {u}_{\mathrm{rms}}^{3}{({k}_{{\rm{f}}0}{H}_{P})}^{1-\beta }$	buoyancy force
${F}_{\mathrm{enth}}\propto {u}_{\mathrm{rms}}{\rm{\Delta }}{\rm{\nabla }}\,/{({k}_{{\rm{f}}0}{H}_{P})}^{1-\tilde{\beta }}$	mean-field expression
${u}_{\mathrm{rms}}\ \propto {({\rm{\Delta }}{\rm{\nabla }})}^{1/2}/{({k}_{{\rm{f}}0}{H}_{P})}^{1-\beta ^{\prime} }$	$\beta ^{\prime} \ =(\beta +\tilde{\beta })/2$
${F}_{\mathrm{enth}}\propto {({\rm{\Delta }}{\rm{\nabla }})}^{3/2}/{({k}_{{\rm{f}}0}{H}_{P})}^{2(1-\beta ^{\prime\prime} )}$	$\beta ^{\prime\prime} =(\beta ^{\prime} +\tilde{\beta })/2=(\beta +3\tilde{\beta })/4$
${F}_{\mathrm{kin}}\,\,\propto {({\rm{\Delta }}{\rm{\nabla }})}^{3/2}/{({k}_{{\rm{f}}0}{H}_{P})}^{3(1-\beta ^{\prime} )}$	kinetic energy flux

Download table as:  ASCIITypeset image

In the first expression for ${F}_{\mathrm{enth}}$, β enters because it characterizes the relation between the buoyancy force proportional to $\delta T/\overline{T}$ and advection proportional to ${k}_{{\rm{f}}}{u}_{\mathrm{rms}}^{2};$ see Appendix D. For thin threads, we expect the relevant ${k}_{{\rm{f}}}$ to be large, i.e., $\beta =0$ (Cases I and II in Figure 2). Next, in the second expression for ${F}_{\mathrm{enth}}$, we have used a mean-field expression to relate ${F}_{\mathrm{enth}}$ to ${\rm{\Delta }}{\rm{\nabla }}$ via a turbulent diffusivity proportional to $\tau {u}_{\mathrm{rms}}^{2}\approx {u}_{\mathrm{rms}}/{k}_{{\rm{f}}}\phi$, where, as argued above, the dilution factor ϕ has entered. It is this expression that is closest to conventional MLT, because here we expect $\tilde{\beta }=1$.

The remaining two relationships in Table 3 explain why the $\tilde{\beta }$ term appears three times more dominantly than the β term. By equating the first two expressions for ${F}_{\mathrm{enth}}$ in Table 3, we find first of all the relation between ${u}_{\mathrm{rms}}$ and ${\rm{\Delta }}{\rm{\nabla }}$, where β and $\tilde{\beta }$ contribute equal shares through $\beta ^{\prime} =(\beta +\tilde{\beta })/2$. However, to see the scaling of ${\rm{\Delta }}{\rm{\nabla }}$, we need to go back to the second expression for ${F}_{\mathrm{enth}}$, because it changes only weakly with depth. Now, $\beta ^{\prime}$ and $\tilde{\beta }$ contribute equal shares, and this means that $\tilde{\beta }$ has now become three times more dominant than β through the expression $\beta ^{\prime\prime} =(\beta +3\tilde{\beta })/4$. Thus, for $\beta =0$ and $\tilde{\beta }=1$, ${\rm{\Delta }}{\rm{\nabla }}$ shows nearly the standard scaling with ${H}_{P}$. Furthermore, looking again at the first expression for ${F}_{\mathrm{enth}}$ in Table 3, we see that only β enters, so the scaling of ${u}_{\mathrm{rms}}$ is fully characterized by that of small blobs with negative buoyancy.

In this section we present numerical solutions to demonstrate the effect of the ${{\rm{\nabla }}}_{{\rm{D}}}$ term on the resulting stratification. At the end of this section, we also compare with the non-convecting reference solution mentioned in the introduction. We should emphasize that, although we use solar parameters, our models cannot represent the Sun. because ionization effects have been ignored (we take $\mu =0.6$ for the mean molecular weight in ${c}_{{\rm{p}}}-{c}_{{\rm{v}}}={ \mathcal R }/\mu$, where ${ \mathcal R }$ is the universal gas constant). A rather simple opacity law of the form of Equation (6), with ${\kappa }_{0}={10}^{4}\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, ${\rho }_{0}={10}^{-5}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, ${T}_{0}={\rm{13,000}}\,{\rm{K}}$, and a = 0.5, b = 18, for the exponents in the power-law expression for ${\kappa }_{{{\rm{H}}}^{-}}$ has been used. We also neglect the departure from plane-parallel geometry, so our model can only give qualitative indications.

The system of two differential Equations (36) and (37) decouples by using $\mathrm{ln}P$ as the independent variable. We thus integrate

$$\begin{eqnarray}&&d\mathrm{ln}\overline{T}/d\mathrm{ln}\overline{P}={\rm{\nabla }},\end{eqnarray} \tag{ 61 }$$

using Equations (32), (33), and (50) to compute ${\rm{\nabla }}$. As an initial condition we use ${\overline{T}}_{\mathrm{top}}={T}_{\mathrm{eff}}/{2}^{1/4}$ at a sufficiently low pressure (here ${\overline{P}}_{\mathrm{top}}={10}^{5}\,\mathrm{dyn}\,{\mathrm{cm}}^{-2}$) so as to capture the initially isothermal part of the atmosphere; see, e.g., Böhm-Vitense (1958) or Mihalas (1978). Here, ${T}_{\mathrm{eff}}={({F}_{\mathrm{tot}}/{\sigma }_{\mathrm{SB}})}^{1/4}$ is the effective temperature. We approximate the ${u}_{\mathrm{rms}}$ term in Equation (47) by using the value from the previous step. The Deardorff term is characterized by the assumed value of f_s0 and the value of $\tilde{\zeta }=\zeta -{\rm{\Delta }}\zeta$ (Section 3.3), where $\zeta =1$ and ${\rm{\Delta }}\zeta$ is a function of β and $\tilde{\beta }$ (Section 4.3). Since we integrate from the top downward, no prior knowledge of ${\rho }_{* }$ is needed, because the Deardorff term is invoked only after ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$ has reached its peak value.

Geometrical and optical depths are obtained respectively as

$$\begin{eqnarray}&&-z=\int (\overline{P}/\overline{\rho }g)d\mathrm{ln}\overline{P}\quad {\rm{and}}\ \ \tau =\int (\kappa \overline{P}/g)d\mathrm{ln}\overline{P},\end{eqnarray} \tag{ 62 }$$

where the integration of τ starts at $\mathrm{ln}{\overline{P}}_{\mathrm{top}}$ and that of $-z$ at the position where $\tau =1$, which is referred to as the surface. The factor $\kappa \overline{P}/g$, which is the same as ${H}_{P}$, is retained because of the similarity with that in the expression for τ. The z coordinate is used in some of our plots.

In the following, we present solutions for the three combinations of β and $\tilde{\beta }$ sketched in Figure 2. We recall that only Cases I and II ($\beta =0$ with $\tilde{\beta }=0$ or 1) correspond to small length scales in the deeper layers (see also the power spectra in Figure 2) and are therefore of interest when trying to reconcile the non-detection of convective motions by Hanasoge et al. (2012, 2016) at the theoretically expected levels. We did already emphasize that Case II with $\beta =0$ and $\tilde{\beta }=1$ is likely to lead to large kinetic energy fluxes. However, based on the results presented below, it turns out that in the case $\beta =\tilde{\beta }=0$ (Case I), which was favored by these two requirements (small length scale and non-divergent kinetic energy flux), the Deardorff term is unlikely to have a significant effect, because for $\tilde{\beta }=0$, the gradient term in the enthalpy flux becomes rather inefficient and must therefore be compensated for by a correspondingly larger superadiabatic gradient, and thus, only a rather large Deardorff flux (${f}_{s0}\geqslant 0.5$) can make the resulting stratification sufficiently subadiabatic. This is the case shown in Figure 4, where we present profiles of $\overline{S}/{c}_{P}$, ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$, and ${u}_{\mathrm{rms}}$ for ${f}_{s0}=0$ (no Deardorff term), as well as ${f}_{s0}=0.2$, 0.3, and 0.5. The ${u}_{\mathrm{rms}}$ profiles are basically the same for all values of f_s0 and fall off like ${\overline{P}}^{-1/3}$. This is expected, because ${u}_{\mathrm{rms}}\propto {\overline{T}}^{1/2-m}$ and $m=4/3;$ see Equation (59), where we have used ${c}_{{\rm{s}}}\propto {\overline{T}}^{1/2}$. Thus, for isentropic stratification we have ${u}_{\mathrm{rms}}\propto {\overline{P}}^{(1-2m)/5}={\overline{P}}^{-1/3}$.

Zoom In Zoom Out Reset image size

Next, we show in Figures 5 and 6 Cases II and III with $\tilde{\beta }=1$ and ${f}_{s0}=0$, 0.1, and 0.2. For both $\beta =0$ and $\beta =1$, the Deardorff flux now has a stronger effect and is able to make the deeper parts of the domain Schwarzschild stable even for ${f}_{s0}=0.1$. Those layers would then be Schwarzschild stable and no longer a source of giant cells. Again, the profiles of ${u}_{\mathrm{rms}}$ are similar regardless of the value of f_s0, but fall off more slowly for $\beta =1$ (${u}_{\mathrm{rms}}\propto {\overline{P}}^{-1/5}$), compared to $\beta =0$ (${u}_{\mathrm{rms}}\propto {\overline{P}}^{-1/3}$). This agrees with our theory, because for m = 1 we have ${u}_{\mathrm{rms}}\propto {\overline{P}}^{(1-2m)/5}={\overline{P}}^{-1/5}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given that the stratification is stable in the deeper parts, we can calculate the Brunt–Väisälä frequency of buoyancy oscillations, ${N}_{\mathrm{BV}}$, which is given by ${N}_{\mathrm{BV}}^{2}=-({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}})g/{H}_{P}$. At intermediate and larger depths of the convection zone, we have $g/{H}_{P}\lesssim {10}^{-2}\,{{\rm{s}}}^{-2}$ and ${{\rm{\nabla }}}_{\mathrm{ad}}-{\rm{\nabla }}\lesssim {10}^{-4}$, so the period of buoyancy oscillations would be of the order of days. This is comparable to or less than the turnover time τ. Indeed, one finds that $\tau {N}_{\mathrm{BV}}\approx \sqrt{{{\rm{\nabla }}}_{{\rm{D}}}/{\rm{\Delta }}{\rm{\nabla }}-1}$ exceeds unity in the deeper parts, which is a consequence of ${{\rm{\nabla }}}_{{\rm{D}}}$ falling off with $\overline{\rho }$ with a smaller power than ${\rm{\Delta }}{\rm{\nabla }}$, as discussed in Section 4.4. Nevertheless, the resulting decrease in $\overline{S}/{c}_{P}$ with depth remains always small (10⁻³ or less) compared with the value of ${\rm{\Delta }}{S}_{0}/{c}_{P}$ produced at the surface ($\approx {u}_{\mathrm{rms}}^{2}/{a}_{\mathrm{MLT}}{{\rm{\nabla }}}_{\mathrm{ad}}{c}_{{\rm{s}}}^{2}\approx 0.01$). Thus, based on this, the descending low-entropy blobs would reach the bottom of the convection zone before their negative buoyancy is neutralized by the decreasing average entropy. In other words, they will never perform any actual oscillations before reaching the bottom of the convection zone, i.e., where ${F}_{\mathrm{conv}}=0$.

It may be interesting to note that the depth of the convection zone increases slightly with increasing values of f_s0. Of course, our model is idealized and represents the Sun at best only approximately. Furthermore, as emphasized in Section 1, the depth of the convection zone is well determined seismically, and this should be reproduced by a solar model with realistic atomic physics and appropriately chosen adjustable parameters. However, it is known from the work of Serenelli et al. (2009) that a slight expansion of the solar convection zone would actually be required to compensate for the shrinking that follows from the downward revision of solar abundances (Asplund et al. 2004) which is based on three-dimensional convective atmosphere simulations, compared to previous analysis based on one-dimensional semi-empirical models (Grevesse & Sauval 1998).

Finally, we compare in Figure 7 the standard convective solution ($\beta =\tilde{\beta }=1;$ same as the case with ${f}_{s0}=0$ in Figure 6) with the non-convective radiative reference solution. Not surprisingly, owing to the absence of convection, the same flux can now only be transported with a greatly enhanced negative (unstable) entropy gradient near the surface. However, this layer is now extremely thin ($1.15\,\mathrm{Mm}$) and the peak in ${{\rm{\nabla }}}_{\mathrm{rad}}/{{\rm{\nabla }}}_{\mathrm{ad}}$ is about $100\,\mathrm{km}$ below the $\tau =1$ surface. Note also that its peak value (≈10⁵) is below that in the presence of convection, where it reaches a maximum of ≈4 × 10⁶ at a depth of $\approx 1\,\mathrm{Mm}$. As shown in Equation (42), this value can be interpreted as the local Nusselt number. Note also that the result for $\beta =\tilde{\beta }=0$ (when ${{\rm{\nabla }}}_{{\rm{D}}}$ is weak) is rather similar to that for $\beta =\tilde{\beta }=1;$ see the dashed and dotted lines in Figure 7.

Zoom In Zoom Out Reset image size

Our calculations have demonstrated that for $\tilde{\beta }=1$, regardless of the value of β, bulk mixing changes the non-convecting reference state to a nearly isentropic one. However, whether the mean entropy gradient is slightly stably or slightly unstably stratified depends on the presence of the Deardorff flux. In the model with $\tilde{\beta }=0$, however, bulk mixing is rather inefficient and the stratification would be Schwarzschild unstable unless an unrealistically large Deardorff flux is invoked.

At the end of the introduction, we discussed that the entropy rain itself might create an unstable stratification. Let us now return to this question with more detailed estimates. This will be done in the following section.

Assuming that the upflows are perfectly isentropic, Spruit (1997) argues that the low-entropy material from the top with its decreasing entropy filling factor f_s in deeper layers necessarily leads to a negative mean entropy gradient. Specifically, using Equation (24), one obtains

$$\begin{eqnarray}&&-{H}_{P}\displaystyle \frac{d\overline{S}/{c}_{P}}{{dz}}=\displaystyle \frac{{({\rm{\Delta }}S)}_{0}}{{c}_{P}}\displaystyle \frac{{{df}}_{s}}{d\mathrm{ln}\overline{P}}=\displaystyle \frac{2}{5}\displaystyle \frac{{({\rm{\Delta }}S)}_{0}}{{c}_{P}}{f}_{s}\gt 0,\end{eqnarray} \tag{ 63 }$$

where we have used $d\mathrm{ln}{f}_{s}/d\mathrm{ln}\overline{\rho }=2/3$ and $d\mathrm{ln}\overline{\rho }/d\mathrm{ln}\overline{P}\ =3/5$ for an isentropic layer with $\gamma =5/3$. Thus, the stratification would be Schwarzschild unstable. This is also borne out by the solar simulations with realistic physics, although the computational domains are sufficiently shallow so that the radiative flux is still small in the deeper parts and usually even neglected altogether. Toward the bottom of the convection zone, however, radiation becomes progressively more important and the mean entropy gradient in the upflows may no longer be vanishing.

To estimate the mean entropy gradient in the upflows, we may balance the steady state entropy advection with the negative radiative flux divergence, i.e.,

$$\begin{eqnarray}&&\overline{\rho }\,\overline{T}\,{\overline{U}}_{\uparrow }\,d{\overline{S}}_{\uparrow }/{dz}\approx -{({{dF}}_{\mathrm{rad}}/{dz})}_{\uparrow }.\end{eqnarray} \tag{ 64 }$$

The sign of ${{dF}}_{\mathrm{rad}}/{dz}$ is negative and thus compatible with a positive $d{\overline{S}}_{\uparrow }/{dz}$ in the upflows, but it would only be large enough near the bottom of the convection zone. Higher up, ${{dF}}_{\mathrm{rad}}/{dz}$ becomes smaller and eventually unimportant. On the other hand, ${\overline{U}}_{\uparrow }/{u}_{\mathrm{rms}}$ can be rather small (see Table 1). Furthermore, our considerations neglect the fact that the gas in the upflows expands, so only a fraction of the gas can ascend before it begins to occupy the available surface area. Therefore, the rest of the gas would have to remain stagnant and continue to heat up. In reality, of course, there would be continuous entrainment, resulting in a finite, but still low, effective upward velocity.

We now use selected solutions obtained in Section 5 to estimate the effective fractional upward velocity for radiative heating/cooling to dominate over advection, defined as

$$\begin{eqnarray}&&{\overline{U}}_{\mathrm{eff}}^{\uparrow }/{u}_{\mathrm{rms}}=(-{{dF}}_{\mathrm{rad}}/{dz})/({u}_{\mathrm{rms}}\,\overline{\rho }\,\overline{T}d\overline{S}/{dz}).\end{eqnarray} \tag{ 65 }$$

Here we have used $\overline{T}d\overline{S}/{dz}=g\,{\rm{\Delta }}{\rm{\nabla }}/{{\rm{\nabla }}}_{\mathrm{ad}}$. Figure 8 shows this quantity for Cases I–III, which are shown in Figures 4–6. It turns out that for Case I with ${f}_{s0}=0.5$, the effective fractional upward velocity is around 0.01, while for Case III with ${f}_{s0}=0.1$ it can be even larger. The value 0.01 is compatible with the ${\overline{U}}^{\uparrow }/{u}_{\mathrm{rms}}$ values in Table 1 if we assume f = 0.01 (so ${\overline{U}}^{\uparrow }/{u}_{\mathrm{rms}}=0.1$) and ${\overline{U}}_{\mathrm{eff}}^{\uparrow }/{\overline{U}}^{\uparrow }=0.1$. For Case II with ${f}_{s0}=0.1$, the effective fractional upward velocity is around 10⁻³, which may be unrealistically small.

Zoom In Zoom Out Reset image size

Based on the considerations above, we may conclude that the case for a Schwarzschild-stable entropy gradient depends on model assumptions for the implementation of the Deardorff flux that are potentially in conflict with the mean entropy gradient expected from the two-stream model. Whether or not there is really a potential conflict depends not only on the model parameters, but also on Equation (64) itself. This is because the $\overline{T}$ factor in the entropy equation is outside the derivative, making it impossible to derive the total flux balance assumed in Equation (39); see also Rogachevskii & Kleeorin (2015) for related details.

In view of these complications, it is worthwhile to discuss alternative ways of avoiding giant cell convection. For this purpose, we have to address the global problem of using some coarse-grained form of effective mean-field equations. When considering the equations of mean-field hydrodynamics, in which the small-scale enthalpy and momentum fluxes are parameterized in terms of negative mean entropy and mean velocity gradients, respectively, one finds solar differential rotation as a result of non-diffusive contributions to the Reynolds stress, but in certain parameter regimes a new instability was found to develop (Rüdiger 1989; Tuominen & Rüdiger 1989; Rüdiger & Spahn 1992). This instability was later identified as one that is analogous to Rayleigh–Bénard convection, but now for an already convecting mean state (Tuominen et al. 1994). This instability would lead to giant cell convection.

The existence of giant cell convection is under debate, but assuming that it does not exist in the Sun, one might either hypothesize that the mean-field equations are too simplified or that mean-field convection could be suppressed by sufficiently strong turbulent viscosity and turbulent thermal diffusivity coefficients. These turbulent coefficients define a turbulent Rayleigh number for a layer of thickness d,

$$\begin{eqnarray}&&{{\rm{Ra}}}_{{\rm{t}}}=\displaystyle \frac{{{gd}}^{4}}{{\nu }_{{\rm{t}}}{\chi }_{{\rm{t}}}}{\left(\displaystyle \frac{{\rm{\Delta }}{\rm{\nabla }}}{{H}_{P}}\right)}_{{z}_{* }}.\end{eqnarray} \tag{ 66 }$$

This number would then have to be still below the critical value for convection. Equation (66) differs from the usual one defined through Equation (86) in Appendix G in that, first, $\nu \to {\nu }_{{\rm{t}}}$ and $\chi \to {\chi }_{{\rm{t}}}$ have been substituted, and second, the superadiabatic gradient of the non-convecting reference solution is now replaced by the actual one. This might be a plausible alternative to explaining the absence of giant cell convection if the idea of turning the stratification from Schwarzschild unstable to Schwarzschild stable through the ${{\rm{\nabla }}}_{{\rm{D}}}$ term were to turn out untenable.

We recapitulate that in this alternate explanation, the stratification is Schwarzschild unstable, i.e., ${{\rm{Ra}}}_{{\rm{t}}}\gt 0$, corresponding to Equation (54), but ${{\rm{Ra}}}_{{\rm{t}}}\lt {{\rm{Ra}}}_{{\rm{t}}}^{\mathrm{crit}}$ is still below a certain critical value, so it would be stable by a turbulent version of the Rayleigh–Bénard criterion. Tuominen et al. (1994) found ${{\rm{Ra}}}_{{\rm{t}}}^{\mathrm{crit}}\approx 300$ for a vanishing rotation rate. However, they also estimated that ${{\rm{Ra}}}_{{\rm{t}}}\gt {{\rm{Ra}}}_{{\rm{t}}}^{\mathrm{crit}}$ for plausible solar parameters, so the direct adoption of this idea would be problematic, too. However, one might speculate that a more accurate treatment could lead to stability when allowing, for example, for spatial nonlocality of the turbulent transport (Brandenburg et al. 2008; Rheinhardt & Brandenburg 2012).

In the present work we have suggested that the enthalpy flux in stellar mixing length models should contain an extra nonlocal contribution so that the enthalpy flux is no longer proportional to the local superadiabatic gradient, ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$, but to ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\nabla }}}_{{\rm{D}}}$, where ${{\rm{\nabla }}}_{{\rm{D}}}$ is a new nonlocal contribution that was first identified by Deardorff (1966) in the meteorological context. The significance of this term lies in the fact that it provides an alternative to the usual local entropy gradient term and can transport enthalpy flux outwards—even in a slightly stably stratified layer.

We have presented a modified formulation of stellar MLT that includes the ${{\rm{\nabla }}}_{{\rm{D}}}$ term, in addition to a ${{\rm{\nabla }}}_{\mathrm{kin}}$ term resulting from the kinetic energy flux. The formalism and the final results are similar to those of conventional MLT in that one also arrives here at a cubic equation for ${\rm{\nabla }}$, but the term ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$ is now replaced by ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\nabla }}}_{{\rm{D}}}-{{\rm{\nabla }}}_{\mathrm{kin}}$. This new formulation implies that convection can carry a finite flux while ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$ is still negative and therefore the stratification is Schwarzschild stable, i.e., Equation (53) is obeyed. Consequently, if confirmed, no large length scales are being excited.

The present formulation allows for different treatments of the length scales governing buoyant elements on the one hand and the time and length scales associated with mixing on the other. When both are independent of depth ($\beta =\tilde{\beta }=0$), mixing becomes inefficient at larger depths. Thus, to carry a certain fraction of the enthalpy flux, the superadiabatic gradient needs to be larger than otherwise, making it harder for the Deardorff term to revert the sign of ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$. On the other hand, for a tree-like hierarchy of many downdrafts merging into fewer thin ones at greater depths ($\beta =0$, $\tilde{\beta }=1$), the increasing length scale associated with increasing separation enhances vertical mixing, making the stratification nearly isentropic without the Deardorff term, and slightly subadiabatic with a weak Deardorff term. In that case, however, the filling factor of the downflows decreases with depth as ${\rho }^{-\zeta }$, which may imply an unrealistically large downward kinetic energy flux. This leaves us with the standard flow topology ($\beta =\tilde{\beta }=1$), where both size and separation of structures increase with depth. The Deardorff flux can still cause the stratification to have a subadiabatic gradient, so no giant cell convection would be excited locally, but the flow structures would be large and should be helioseismically detectable, as has been found by Greer et al. (2015) using ring-diagram local helioseismology.

It would be useful to explore the thermodynamic aspects of the present model more thoroughly and to connect with related approaches. An example is the work by Rempel (2004), who studied a semianalytic overshoot model that was driven nonlocally by downdraft plumes, similar to what was suggested by Spruit (1997). Rempel also finds an extended subadiabatic layer in large parts of the model. Furthermore, there are similarities to the nonlocal mixing length model by Xiong & Deng (2001), who, again, find an extended deeper layer that is subadiabatically stratified in the deeper parts of their model. In this connection we emphasize the main difference between nonlocal turbulence owing to the Deardorff flux and usual overshoot: in the latter case the enthalpy flux would go inward as a consequence of the reversed entropy gradient, while in the present model the Deardorff flux goes outward.

Future work could proceed along two separate paths. On the one hand, one must establish the detailed physics leading to the ${{\rm{\nabla }}}_{{\rm{D}}}$ term using models with reduced opacity, in which reliable DNS are still possible, i.e., no SGS terms are added and the primitive equations are solved as stated, without invoking Equation (1). On the other hand, one could study suitably parameterized large eddy simulations that either include a nonlocal Deardorff term of the form given by Equation (32), as discussed in Section 3, or that explicitly release entropy rain at the surface such that the resulting stratification is still slightly stable. This would be particularly useful in global simulations that would otherwise have no entropy rain.

As stimulating as the results of Hanasoge et al. (2012) are, they do require further scrutiny and call for the resolution of the existing conflicts with other helioseismic studies such as those of Greer et al. (2015). Alternatively, global helioseismic techniques for detecting giant cell convection (Lavely & Ritzwoller 1993; Chatterjee & Antia 2009; Woodard 2014) can provide another independent way of detecting deep larger-scale flows (Woodard 2016). Realistic simulations should eventually agree with helioseismic results of flows in deeper layers of the Sun, but at the moment it is still unclear whether a subgrid scale treatment as in Equation (1) adequately captures the small-scale flows that can be responsible for the Deardorff flux and whether they would in principle be able to predict the subtle departures from superadiabatic stratification on subthermal timescales.

I thank the referee for criticism that has led to many improvements, and Evan Anders, Jean-Francois Cossette, Ben Greer, Åke Nordlund, Mark Rast, Matthias Rempel, Matthias Rheinhardt, Bob Stein, Peter Sullivan, Regner Trampedach, and Jörn Warnecke for interesting discussions and comments. I also wish to acknowledge Juri Toomre for mentioning to me the need for modeling entropy rain during my time in Boulder some 24 years ago. This work was supported in part by the Swedish Research Council grant No. 2012-5797, and the Research Council of Norway under the FRINATEK grant 231444. This work utilized the Janus supercomputer, which is supported by the National Science Foundation (award number CNS-0821794), the University of Colorado Boulder, the University of Colorado Denver, and the National Center for Atmospheric Research. The Janus supercomputer is operated by the University of Colorado Boulder.

We show here that, for the non-convecting reference solution, using the Kramers-like opacity law of Equation (2), but not the combined opacity law of Equation (6), we have $K\to \mathrm{const}$ in the deeper, optically thick layers. Dividing Equation (3) by the equation for hydrostatic equilibrium, ${dP}/{dz}=-\rho g$, we have

$$\begin{eqnarray}\displaystyle \frac{{dT}}{{dP}} & = & \displaystyle \frac{{F}_{\mathrm{rad}}}{K\rho g}=\displaystyle \frac{{F}_{\mathrm{rad}}}{{K}_{0}{\rho }_{0}g}\displaystyle \frac{{(\rho /{\rho }_{0})}^{a}}{{(T/{T}_{0})}^{3-b}}\\ & = & \displaystyle \frac{{F}_{\mathrm{rad}}}{{K}_{0}{\rho }_{0}g}\displaystyle \frac{{(P/{P}_{0})}^{a}}{{(T/{T}_{0})}^{3+a-b}},\end{eqnarray} \tag{ 67 }$$

where ${K}_{0}=16{\sigma }_{\mathrm{SB}}{T}_{0}^{3}/(3{\kappa }_{0}{\rho }_{0})$ is a constant and $P/{P}_{0}=(\rho /{\rho }_{0})(T/{T}_{0})$ is the ideal gas equation with a suitably defined constant ${P}_{0}=({c}_{{\rm{p}}}-{c}_{{\rm{v}}}){\rho }_{0}{T}_{0}$. Here, ${\rho }_{0}$ and T₀ are reference values that were defined in Equation (2). Equation (67) can be integrated to give

$$\begin{eqnarray}&&{(T/{T}_{0})}^{4+a-b}=(n+1){{\rm{\nabla }}}_{\mathrm{rad}}^{(0)}{(P/{P}_{0})}^{1+a}+{({T}_{\mathrm{top}}/{T}_{0})}^{4+a-b},\end{eqnarray} \tag{ 68 }$$

where ${{\rm{\nabla }}}_{\mathrm{rad}}^{(0)}={F}_{\mathrm{rad}}{P}_{0}/({K}_{0}{T}_{0}{\rho }_{0}g)$, which is defined analogously to the ${{\rm{\nabla }}}_{\mathrm{rad}}$ without superscript $(0)$ in Equations (38) and (41), and ${T}_{\mathrm{top}}$ is an integration constant that is specified such that $T\to {T}_{\mathrm{top}}$ as $P\to 0$. Note also that $4+a-b=(n+1)(1+a)$, where n was defined in Equation (5) as the polytropic index, so the ratio of $4+a-b$ to $1+a$ is just $n+1$, which enters in front of the ${{\rm{\nabla }}}_{\mathrm{rad}}^{(0)}$ term in Equation (68). Since $K\propto {T}^{3-b}/{\rho }^{1+a}\propto {T}^{4+a-b}/{P}^{1+a}$, we have $K\to \mathrm{const}={K}_{0}$ for $T\gg {T}_{\mathrm{top}}$.

To obtain Equations (11) and (12), which are used to derive the Deardorff flux term in the τ approximation, we start with the equations for specific entropy and velocity in the form (see, e.g., Barekat & Brandenburg 2014)

$$\begin{eqnarray}&&\rho T\displaystyle \frac{{DS}}{{Dt}}=-{\boldsymbol{\nabla }}\cdot {{\boldsymbol{F}}}_{\mathrm{rad}},\end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray}&&\rho \displaystyle \frac{D{\boldsymbol{U}}}{{Dt}}=-{\boldsymbol{\nabla }}P+\rho {\boldsymbol{g}},\end{eqnarray} \tag{ 70 }$$

where $D/{Dt}=\partial /\partial t+{\boldsymbol{U}}\cdot {\boldsymbol{\nabla }}$ is the advective derivative, ${{\boldsymbol{F}}}_{\mathrm{rad}}$ is the radiative flux, and viscosity has been omitted. Subtracting those equations from their averaged ones, we obtain the following set of equations

$$\begin{eqnarray}&&\displaystyle \frac{\partial s}{\partial t}+{\overline{U}}_{j}\displaystyle \frac{\partial s}{\partial {x}_{j}}+{u}_{j}\displaystyle \frac{\partial \overline{S}}{\partial {x}_{j}}=-\displaystyle \frac{1}{\overline{\rho }\,\overline{T}}{\boldsymbol{\nabla }}\cdot \delta {{\boldsymbol{F}}}_{\mathrm{rad}}+{{ \mathcal N }}_{s},\end{eqnarray} \tag{ 71 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {u}_{i}}{\partial t}+{\overline{U}}_{j}\displaystyle \frac{\partial {u}_{i}}{\partial {x}_{j}}+{u}_{j}\displaystyle \frac{\partial {\overline{U}}_{i}}{\partial {x}_{j}}\\ &&\quad =-\displaystyle \frac{1}{\overline{\rho }}{\boldsymbol{\nabla }}p+{g}_{i}(p/\gamma \overline{P}-s/{c}_{P})+{{ \mathcal N }}_{{ui}},\end{eqnarray} \tag{ 72 }$$

where we have used $\delta \rho /\overline{\rho }=p/\gamma \overline{P}-s/{c}_{P}$, which can be obtained from Equation (8) and the perfect gas law by linearization. Pressure fluctuations will again be neglected and $({\boldsymbol{\nabla }}\cdot \delta {{\boldsymbol{F}}}_{\mathrm{rad}})/\overline{\rho }\,\overline{T}$ will be replaced by $-s/{\tau }_{\mathrm{cool}}$, as explained in Section 3. Assuming $\overline{{\boldsymbol{U}}}={\bf{0}}$ and omitting the nonlinear terms ${{ \mathcal N }}_{s}$ and ${{ \mathcal N }}_{{ui}}$ we arrive at Equations (11) and (12).

The solution for a Hill vortex with radius ${a}_{{\rm{H}}}$ and propagation velocity ${u}_{{\rm{H}}}$ is given by a stream function Ψ in spherical coordinates $(r,\theta ,\phi )$ as (e.g., Moffatt & Moore 1978)

$$\begin{eqnarray}{\rm{\Psi }}=\displaystyle \frac{1}{4}\left\{\begin{array}{ll}-3{u}_{{\rm{H}}}(1-{r}^{2}/{a}_{{\rm{H}}}^{2})\varpi & ({\rm{for}}\,{\rm{}}r\lt {a}_{{\rm{H}}}),\\ +2{u}_{{\rm{H}}}(1-{a}_{{\rm{H}}}^{3}/{r}^{3})\varpi & ({\rm{for}}\,{\rm{}}r\gt {a}_{{\rm{H}}}).\end{array}\right.\end{eqnarray} \tag{ 73 }$$

We apply it in Cartesian coordinates as the initial condition for the mass flux as $\rho {\boldsymbol{u}}={\boldsymbol{\nabla }}\times ({\rm{\Psi }}\hat{{\boldsymbol{\phi }}})$, where $\hat{{\boldsymbol{\phi }}}\ =(-y/\varpi ,x/\varpi ,0)$ is the unit vector in the toroidal direction of the vortex, using ${\varpi }^{2}={x}^{2}+{y}^{2}$ and ${r}^{2}={\varpi }^{2}+{z}^{2}$ for cylindrical radius ϖ and spherical radius r. We adopt here an isothermal equation of state, i.e., there is no buoyancy force in this problem. For non-isothermal calculations, but in two dimensions, we refer to the work of Rast (1998). Density and pressure fall off exponentially with height and both the scale height and the sound speed are independent of height. No analytic solution exists in that case, so the Hill vortex solution is at best approximate. We consider a domain of size $L\times L\times 4L$ with $-L/2\lt x,y/{H}_{P}\lt L/2$ and $-7L/2\ \lt z\lt L/2$, where $L=5{H}_{P}$. We choose ${a}_{{\rm{H}}}=0.5{H}_{P}$ and ${u}_{{\rm{H}}}=0.2{c}_{{\rm{s}}}$. The viscosity is $\nu =5\times {10}^{-5}$, so the Reynolds number is ${a}_{{\rm{H}}}{u}_{{\rm{H}}}/\nu =2000$. We use the Pencil Code ⁷ with a resolution of $1152\times 1152\times 4608$ meshpoints.

Figure 9 shows snapshots zoomed into the vortex as it traverses about five scale heights. The filling factor, which is proportional to the radius squared, decreases with depth and is found to scale with the surrounding density like Equation (23) with $\zeta =0.8;$ see the last panel of Figure 9.

Zoom In Zoom Out Reset image size

To find the relation between ${u}_{\mathrm{rms}}$ and the rms values of temperature or entropy fluctuations, it is customary in standard MLT to approximate the steady state momentum equation, ${\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}{\boldsymbol{u}}\approx -{\boldsymbol{g}}\,\delta T/\overline{T}$, by ${u}_{\mathrm{rms}}^{2}{k}_{{\rm{f}}}={a}_{\mathrm{MLT}}\,g\,\delta T/\overline{T}$, where ${a}_{\mathrm{MLT}}\approx 1/8$ is a commonly adopted geometric factor (e.g., Spruit 1974). This leads to

$$\begin{eqnarray}&&{s}_{\mathrm{rms}}/{c}_{P}=\delta T/T={u}_{\mathrm{rms}}^{2}{k}_{{\rm{f}}}/{a}_{\mathrm{MLT}}g=\gamma {{\rm{Ma}}}^{2}{k}_{{\rm{f}}}{H}_{P}/{a}_{\mathrm{MLT}},\end{eqnarray} \tag{ 74 }$$

where we have used ${c}_{{\rm{s}}}^{2}=\gamma {{gH}}_{P}$, and thus $({s}_{\mathrm{rms}}/{c}_{{\rm{p}}}){{\rm{Ma}}}^{-2}\ =\gamma {k}_{{\rm{f}}}{H}_{P}/{a}_{\mathrm{MLT}}$. Furthermore, using ${{\rm{\nabla }}}_{\mathrm{ad}}{c}_{{\rm{p}}}\overline{T}={{gH}}_{P}$, we have $\overline{T}{s}_{\mathrm{rms}}={u}_{\mathrm{rms}}^{2}{k}_{{\rm{f}}}{H}_{P}/({a}_{\mathrm{MLT}}{{\rm{\nabla }}}_{\mathrm{ad}})$, which is used to derive Equation (31).

The purpose of this section is to show that ${{\rm{\nabla }}}_{{\rm{D}}}$ is a certain fraction of ${\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}$ in the top layers, as stated in Equation (32). To have an estimate for $\overline{{s}^{2}}$, we multiply Equation (11) by s and average, so we get

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\displaystyle \frac{\partial \overline{{s}^{2}}}{\partial t}=-\overline{{u}_{j}s}\,{{\rm{\nabla }}}_{j}\overline{S}-\overline{{s}^{2}}/{\tau }_{\mathrm{cool}}+{{ \mathcal T }}_{s}.\end{eqnarray} \tag{ 75 }$$

As for ${F}_{\mathrm{enth}}$ in Equations (14) and (15), we have a triple correlation term, which is here ${{ \mathcal T }}_{s}=-\overline{s{\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}s}$. Again, we adopt the τ approximation and replace ${{ \mathcal T }}_{s}$ by a damping term of the form ${{ \mathcal T }}_{s}=-\overline{{s}^{2}}/\tau$, which, together with the ${\tau }_{\mathrm{cool}}$ term, combines to give ${\tau }_{\mathrm{red}}$ as the relevant timescale. Assuming a statistically steady state, $\partial \overline{{s}^{2}}/\partial t=0$, we derive the following expression for $\overline{{s}^{2}}$:

$$\begin{eqnarray}&&\overline{{s}^{2}}=-{\tau }_{\mathrm{red}}\,\overline{{u}_{j}s}\,{{\rm{\nabla }}}_{j}\overline{S}\quad {\rm{(in}}\,{\rm{the}}\,{\rm{surface}}\,{\rm{layers)}}.\end{eqnarray} \tag{ 76 }$$

This shows that fluctuations of specific entropy are produced when there is an outward flux ($\overline{{u}_{z}s}\gt 0$) and a locally negative (unstable) mean entropy gradient; see also Garaud et al. (2010) for a similar derivation. Inserting this into Equation (18) and using Equation (9), we obtain

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{\rm{D}}}={\tau }_{\mathrm{red}}^{2}\,{\boldsymbol{g}}\,({{\boldsymbol{F}}}_{\mathrm{enth}}\cdot {\boldsymbol{\nabla }}\overline{S}/{c}_{P})\quad {\rm{(in}}\,{\rm{the}}\,{\rm{surface}}\,{\rm{layers)}}.\end{eqnarray} \tag{ 77 }$$

Here, both ${\boldsymbol{g}}$ and ${\boldsymbol{\nabla }}\overline{S}$ point downward, so ${{\boldsymbol{F}}}_{{\rm{D}}}$ points upward and we can write ${F}_{{\rm{D}}}=\lambda {F}_{\mathrm{enth}}$, where $\lambda ={({\tau }_{\mathrm{red}}/{\tau }_{\mathrm{ff}})}^{2}({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}})$ is a coefficient that itself is proportional to the superadiabatic gradient, and ${\tau }_{\mathrm{ff}}={({H}_{P}/g)}^{1/2}$ is the free-fall time (after which a fluid parcel at rest has reached a depth of ${H}_{P}/2$). Since ${F}_{\mathrm{enth}}={F}_{{\rm{G}}}+{F}_{{\rm{D}}}$, this implies ${F}_{{\rm{G}}}=(1-\lambda ){F}_{\mathrm{enth}}$ and therefore $\lambda \lt 1$ in the highly unstable layer at the top, where both ${F}_{{\rm{G}}}$ and ${F}_{\mathrm{enth}}$ are positive. We expect ${F}_{{\rm{D}}}$ to be largest just a few hundred kilometers below the photosphere, so λ should be maximum in the upper parts. To calculate the fraction ${{\rm{\nabla }}}_{{\rm{D}}}/({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}})$, we use the fact that ${F}_{\mathrm{enth}}={F}_{{\rm{G}}}+{F}_{{\rm{D}}}$, together with the part of Equation (20) that relates to ${F}_{{\rm{D}}}$, to write

$$\begin{eqnarray}{F}_{{\rm{D}}} & = & \tfrac{1}{3}\overline{\rho }{c}_{P}\overline{T}\,({\tau }_{\mathrm{red}}{u}_{\mathrm{rms}}^{2}/{H}_{P}){{\rm{\nabla }}}_{{\rm{D}}}\\ & = & {({\tau }_{\mathrm{red}}/{\tau }_{\mathrm{ff}})}^{2}({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}){F}_{\mathrm{enth}}.\end{eqnarray} \tag{ 78 }$$

Starting with Equation (31) and expressing ${\tau }_{\mathrm{red}}$ in terms of σ using Equation (47), we find

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{{\rm{D}}}^{\max }/{({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}})}_{\max }=3\sigma {\phi }_{\mathrm{enth}}{{\rm{\nabla }}}_{\mathrm{ad}}\ {\alpha }_{\mathrm{mix}}^{-\tilde{\beta }}{({k}_{{\rm{f}}0}{H}_{P})}^{-(1-\tilde{\beta })}\approx 3,\end{eqnarray} \tag{ 79 }$$

where we have used $\tilde{\beta }=1$ (appropriate to the near-surface layers), ignored radiative cooling ($\sigma =1$), and assumed ${\alpha }_{\mathrm{mix}}=1.6$ (Stix 2002), as well as ${\phi }_{\mathrm{enth}}\approx 4$ that was found in simulations of Brandenburg et al. (2005), as discussed in Section 3.3.

We emphasize that these equations only characterize the initiation of entropy rain. They cannot be used to compute the Deardorff flux in the deeper layers where we have instead invoked Spruit's concept of nonlocal convection in the form of threads. The $\overline{{s}^{2}}$ associated with those threads is likely to result from a part of the triple correlation term ${{ \mathcal T }}_{s}=-\overline{s{\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}s}$, which gives rise to a negative divergence of the flux of $\overline{{s}^{2}}$ of the form ${\boldsymbol{H}}\equiv \overline{{\boldsymbol{u}}{s}^{2}}$ on the rhs of Equation (75). This flux (which should not be confused with an energy flux) should point downward (s² is largest when ${u}_{z}\lt 0$) and be strongest in the upper layers, so ${\boldsymbol{\nabla }}\cdot {\boldsymbol{H}}\lt 0$, leading to a positive contribution from $-{\boldsymbol{\nabla }}\cdot {\boldsymbol{H}}$ on the rhs of Equation (75). This may explain the $\overline{{s}^{2}}$ associated with the ${{\rm{\nabla }}}_{{\rm{D}}}$ term.

To estimate the ratio ${F}_{{\rm{D}}}/{F}_{{\rm{G}}}$ in the deeper layers, we use Equation (75) in the steady state with $1/{\tau }_{\mathrm{cool}}\to 0$ and ${{ \mathcal T }}_{s}=-{\boldsymbol{\nabla }}\cdot {\boldsymbol{H}}\approx (2\zeta -{\rm{\Delta }}\zeta )\overline{{s}^{2}}{u}_{\mathrm{rms}}/\gamma {H}_{P}$ and obtain

$$\begin{eqnarray}&&0=-\overline{{u}_{j}s}\,{{\rm{\nabla }}}_{j}\overline{S}+(2\zeta -{\rm{\Delta }}\zeta )\overline{{s}^{2}}{u}_{\mathrm{rms}}/\gamma {H}_{P}\,{\rm{(in}}\,{\rm{deeper}}\,{\rm{layers)}}.\end{eqnarray} \tag{ 80 }$$

Multiplying by $\tfrac{1}{3}{\tau }_{\mathrm{red}}{u}_{\mathrm{rms}}^{2}{(\overline{\rho }\overline{T})}^{2}$, using Equations (9) and (17), and expanding the fraction by $g/{c}_{{\rm{p}}}$, we have

$$\begin{eqnarray}&&0={F}_{\mathrm{enth}}\,{F}_{{\rm{G}}}+\tfrac{1}{3}\overline{\rho }{u}_{\mathrm{rms}}^{3}(2\zeta -{\rm{\Delta }}\zeta )({\tau }_{\mathrm{red}}\overline{{s}^{2}}g\,\overline{\rho }\,\overline{T})\,{c}_{{\rm{p}}}\overline{T}/(\gamma {{gH}}_{P}).\end{eqnarray} \tag{ 81 }$$

This implies that ${F}_{{\rm{G}}}\lt 0$ in the deeper layers. The second term in the parentheses is ${F}_{{\rm{D}}}$ (with a plus sign, because $g\gt 0$ is a scalar here); see Equation (18). Furthermore, we use ${c}_{{\rm{p}}}\overline{T}/(\gamma {{gH}}_{P})=1/(\gamma -1)$ and $\overline{\rho }{u}_{\mathrm{rms}}^{3}={F}_{\mathrm{enth}}/{\phi }_{\mathrm{enth}};$ see Equation (31). The ${F}_{\mathrm{enth}}$ terms on both sides cancel, so we have

$$\begin{eqnarray}&&{F}_{{\rm{D}}}/| {F}_{{\rm{G}}}| =3(\gamma -1)/[{\phi }_{\mathrm{enth}}(2\zeta -{\rm{\Delta }}\zeta )]\approx 0.3\,\ldots \,0.5,\end{eqnarray} \tag{ 82 }$$

where $2\zeta -{\rm{\Delta }}\zeta =2\tilde{\zeta }+{\rm{\Delta }}\zeta =10/9$ for Cases I and II, and 14/9 for Case III (Section 4.3), and ${\phi }_{\mathrm{enth}}=4$ has been assumed.

To find the coefficient c₀ given by Equation (49), we use Equation (74) in the form

$$\begin{eqnarray}&&{c}_{P}\delta T={u}_{\mathrm{rms}}^{2}{k}_{{\rm{f}}}{H}_{P}/{{\rm{\nabla }}}_{\mathrm{ad}}{a}_{\mathrm{MLT}}.\end{eqnarray} \tag{ 83 }$$

Inserting this into Equation (7), and using Equation (45), yields

$$\begin{eqnarray}&&{F}_{\mathrm{enth}}=\overline{\rho }{u}_{\mathrm{rms}}^{3}\,{\alpha }_{\mathrm{mix}}^{\beta }{({k}_{{\rm{f}}0}{H}_{P})}^{1-\beta }/{{\rm{\nabla }}}_{\mathrm{ad}}{a}_{\mathrm{MLT}}.\end{eqnarray} \tag{ 84 }$$

Equating this expression with Equation (43), using Equation (46), we can derive the desired expression for ${u}_{\mathrm{rms}}$ in the form

$$\begin{eqnarray}&&{u}_{\mathrm{rms}}^{2}=\displaystyle \frac{\sigma }{3}\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{\gamma }\displaystyle \frac{{a}_{\mathrm{MLT}}}{{({k}_{{\rm{f}}0}{H}_{P})}^{2-2\beta ^{\prime} }}({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\nabla }}}_{{\rm{D}}}),\end{eqnarray} \tag{ 85 }$$

where $\beta ^{\prime} =(\beta +\tilde{\beta })/2$ and ${c}_{{\rm{s}}}^{2}=\gamma {{gH}}_{P}$ have been used. We thus find the coefficient c₀ as stated in Equation (49).

The purpose of this section is to show that ${\epsilon }_{* }$, as defined in Equation (51), is related to the Rayleigh number ${\rm{Ra}}$. In laboratory and numerical studies of convection, it is customary to define ${\rm{Ra}}$ as (e.g., Käpylä et al. 2009)

$$\begin{eqnarray}{\rm{Ra}} & = & \displaystyle \frac{{{gd}}^{4}}{\nu \chi }{\left(-\displaystyle \frac{d\overline{S}/{c}_{P}}{{dz}}\right)}_{{z}_{* }}^{\mathrm{non} \mbox{-} \mathrm{conv}}\approx \displaystyle \frac{{{gd}}^{4}}{\nu \chi }{\left(\displaystyle \frac{{{\rm{\nabla }}}_{\mathrm{rad}}-{{\rm{\nabla }}}_{\mathrm{ad}}}{{H}_{P}}\right)}_{\max }^{\mathrm{non} \mbox{-} \mathrm{conv}}\\ & \approx & \displaystyle \frac{{{gd}}^{4}}{\nu \chi }{\left(\displaystyle \frac{{{\rm{\nabla }}}_{\mathrm{rad}}-{{\rm{\nabla }}}_{\mathrm{ad}}}{{H}_{P}}\right)}_{\max },\end{eqnarray} \tag{ 86 }$$

which is usually evaluated in the middle of the layer at $z={z}_{* }$. In the Sun, however, the maximum value is more relevant. The superscript "non-conv" indicates that the entropy gradient is taken for the non-convecting reference state, d is the thickness of the layer and ν is the viscosity. Introducing the Prandtl number ${\rm{\Pr }}=\nu /\chi$ and using the definition of ${\rm{Nu}}$ in Equation (42), we have

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Pr }}\,{\rm{Ra}}}{{\rm{Nu}}-1}=\displaystyle \frac{{{gd}}^{4}{{\rm{\nabla }}}_{\mathrm{ad}}}{{\chi }^{2}{H}_{P}}.\end{eqnarray} \tag{ 87 }$$

On the other hand, using Equation (51), we find

$$\begin{eqnarray}{\epsilon }_{* }^{2} & = & \displaystyle \frac{{\sigma }^{3}{a}_{\mathrm{MLT}}}{27{\gamma }^{3}}\displaystyle \frac{{({c}_{{\rm{s}}}^{3}/\chi g)}^{2}}{{({k}_{{\rm{f}}0}{H}_{P})}^{4(1-\beta ^{\prime\prime} )}}=\displaystyle \frac{{\sigma }^{3}{a}_{\mathrm{MLT}}}{27{({k}_{{\rm{f}}0}{H}_{P})}^{4(1-\beta ^{\prime\prime} )}}\displaystyle \frac{{H}_{P}^{3}g}{{\chi }^{2}}\\ & = & \displaystyle \frac{{\sigma }^{3}{a}_{\mathrm{MLT}}{({H}_{P}/d)}^{4}}{27{{\rm{\nabla }}}_{\mathrm{ad}}{({k}_{{\rm{f}}0}{H}_{P})}^{4(1-\beta ^{\prime\prime} )}}\displaystyle \frac{{{gd}}^{4}{{\rm{\nabla }}}_{\mathrm{ad}}}{{\chi }^{2}{H}_{P}},\end{eqnarray} \tag{ 88 }$$

and therefore

$$\begin{eqnarray}&&{\epsilon }_{* }^{2}=\displaystyle \frac{{\sigma }^{3}{a}_{\mathrm{MLT}}/27{{\rm{\nabla }}}_{\mathrm{ad}}}{{({k}_{{\rm{f}}0}d)}^{4}{({k}_{{\rm{f}}0}{H}_{P})}^{-4\beta ^{\prime\prime} }}\displaystyle \frac{{\rm{\Pr }}\,{\rm{Ra}}}{{\rm{Nu}}-1},\end{eqnarray} \tag{ 89 }$$

which shows that ${\rm{Ra}}$ is proportional to ${\epsilon }_{* }^{2}$.