\section{Introduction}
\label{Intro}

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 \citep{SN89,SN98,VSSCE05,Gudiksen,Freytag12}.
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_{\rm rms}\approx(F_{\rm enth}/\rho)^{1/3}$, where $F_{\rm enth}$
is the enthalpy flux and $\rho$ is the density.
Simulations also seem to confirm an important {\em 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 \citep{Duvall97,GB05}.
Helioseismic observations by \cite{HBBG10,Hanasoge} using the
deep-focusing time-distance technique have not, however, detected
large-scale convection velocities at the expected levels \citep{MFRT12};
see also \cite{GB12} for a comparison with both global simulations
and radiation hydrodynamics simulations with realistic radiation and
ionization physics included.
\cite{Greer} have suggested that the approach of \cite{Hanasoge}
may remove too much signal over the time span of the measurements.
Using instead ring diagram analysis \citep{GB05} with appropriately
assembled averaging kernels to focus on the deeper layers,
\cite{Greer} 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\,{\rm Mm}$ and larger), the rms velocity actually
{\em increases} with depth.
This in itself is remarkable, because velocity perturbations from
deeper layers were expected to be transmitted to the surface almost
unimpededly \citep{Sti81}, unlike perturbations in the heat flux,
which are efficiently being screened by the convection \citep{Spr77}.
Thus, an increase in horizontal velocity power with depth
is unexpected.
However, \cite{vanBall86} 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 \citep{Hathaway}.
The presence of the near-surface shear layer of the Sun \citep{Schou98},
might enhance the screening further.

Given that the deeply focused kernels used by \cite{Greer} 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 \cite{Greer} to realistic simulations.
Unfortunately, this has not yet been attempted.
Conversely, the deep-focusing technique of \cite{Hanasoge} 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 \cite{Greer}.

It should be mentioned that small flow speeds of giant cell convection
have been found by correlation tracking of supergranule proper
motions \citep{Hathaway}.
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 \cite{Hanasoge}, but \cite{HS14} 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 \cite{Greer}, which show an {\em increase} with depth at
low spherical harmonic degrees.
On the other hand, \cite{BBB15} have argued that the flows reported by
\cite{Hathaway} 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.
\cite{SBN07} 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 \cite{NSA09}.
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 \citep[e.g.][]{Greer}
or smaller \citep[e.g.][]{FH16} 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 $\propto k^2$ or maybe even with a Batchelor
spectrum $\propto k^4$ \citep{Dav04}.
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}$ \citep{Los13}.
By contrast, if one imagines the flow to be anelastic so that the mass
flux $\rho\mbox{\boldmath $u$} {}=\mbox{\boldmath $\nabla$} {}\times\bm{\psi}$ can be written as a vectorial stream
function $\bm{\psi}$, and if $\bm{\psi}$ is given by white noise, one should
expect a $k^4$ Batchelor spectrum.
This is qualitatively similar to the results of granule tracking,
which reveal an intermediate scaling proportional to $k^3$
\citep{Rieutord08,Roudier12}.
While the steeper $k^4$ spectrum for $k<k_{\rm f}$ might leave some hope that
the results of \cite{Hanasoge} 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 \cite{Lord}, who have shown that
simulation results with the MURaM code \citep{VSSCE05} 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 \citep{Lord14} 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\,{\rm 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 being 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 \cite{Hotta}, who have discussed
the possible role of small-scale magnetic fields in suppressing the
formation of large-scale flows.
Furthermore, \cite{FH16} 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{equation}
\mbox{\boldmath $F$} {}_{\rm rad}=-K\mbox{\boldmath $\nabla$} {} T \quad \rightarrow \quad
-K\mbox{\boldmath $\nabla$} {}\overline{T}-K_{\rm SGS}\mbox{\boldmath $\nabla$} {}(T-\overline{T}),
\label{FradSplit}
\end{equation}
where $T$ and $\overline{T}$ are the actual and horizontally averaged
temperatures, respectively, $K$ is the radiative conductivity,
and $K_{\rm 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 {\em direct} numerical simulations (DNS), one uses
physical viscosity and diffusivity operators, i.e., the replacement
in \Eq{FradSplit} is not invoked, but the coefficients are enhanced
($K\to K_{\rm enh}$) and exceed the physical
ones by many orders of magnitude.
Both approaches are problematic.

In many global DNS with enhanced coefficients \citep[e.g.,][]{KMCWB13},
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 \citep{KMCWB13}.
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 \citep{BCNS05}, 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\mbox{\boldmath $\nabla$} {} T$, compared with $-K_{\rm enh}\mbox{\boldmath $\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\'eclet number (which is similar to the Reynolds number)
based on the physical value of $K$ would reach values above
$10^{10}$, which cannot be handled by a DNS \citep{BB14}.
In the global simulations, it is common to use the specific entropy gradient
instead of $\mbox{\boldmath $\nabla$} {}(T-\overline{T})$ \citep[e.g.][]{KMCWB13}.
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 {\em 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 what would
be compatible with the observations discussed by \cite{Lord}.
The sensitivity of large-scale flow excitation to small changes
in the superadiabatic gradient was also discussed by \cite{CR16},
who studied models with different types of surface driving.
It is further supported by the work of \cite{FH16}, 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 \citep{HTM86,BJNRST96}.
The simulations of \cite{KMCWB13} adopted a profile that yields
a negative (unstable) radial entropy gradient through most of
the layer (see the inset of their Fig.~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 \citep{CDGT91,Basu}.
Nevertheless, the concept of convection being driven by surface cooling
rather than heating from below has been promoted in a number of papers
by \cite{SN89,SN98} and elaborated upon by \cite{Spr97}, 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 \cite{Spr97}, we
must look for a contribution to the enthalpy flux that is not related
to the local entropy gradient.\footnote{As will be discussed in
\Sec{UnstablyStratified}, \cite{Spr97} argues that the low-entropy
material from the top always leads to a negative mean entropy gradient.
However, here we argue that the mean entropy gradient is only one
contribution to a mean-field (here one-dimensional)
parameterization of the enthalpy flux, and that
there is another one that is not locally connected with it.}
Such a term has been identified in the meteorological context by \cite{Dea66},
who describes it as a counter-gradient flux.
In \cite{Dea72}, 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 \cite{Spr97} refers to as {\it threads}, which are thin downdrafts on
an almost perfectly isentropic passive fluid upflow between the threads.
Spruit already emphasized back then (p.\ 406) the dynamical importance of the
``forest of narrow cool threads that are produced at the solar surface.''
He further suggested (p.\ 406) 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 \cite{Spr97}, convection is driven solely
within the surface layers.
This is superficially reminiscent of convective overshoot, as modeled
in some of the aforementioned papers \citep{HTM86,BJNRST96},
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 \citep{HTMZ94}, 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 \citep{HTM86}.
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 \cite{Spr97}, 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 {\em 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 \citep{Miesch08}, 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 \cite{Spr97} to the end of the paper
(\Sec{UnstablyStratified}), and begin by highlighting the less
commonly known fact that in the non-convecting reference model, only the
layer in the top $1\,{\rm Mm}$ is convectively unstable (\Sec{HighlyUnstable}).
We then explain the nature of the Deardorff flux (\Sec{Deardorff}),
present a correspondingly modified mixing length model (\Sec{Modified})
and give some illustrative numerical solutions (\Sec{NumSol}).
We discuss alternative explanations for the lack of giant cell convection in
\Sec{UnstablyStratified} and present our conclusions in \Sec{Conclusions}.

\section{A highly unstable surface layer}
\label{HighlyUnstable}

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_{\rm 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 \Sec{NumSol}
and present in this section the basic argument only at a qualitative
level, making reference to earlier numerical calculations by \cite{BB14}
for a simple model.
They considered an opacity law of the form
\begin{equation}
\kappa=\kappa_0 (\rho/\rho_0)^a (T/T_0)^b,
\label{Kramers}
\end{equation}
where $a$ and $b$ are adjustable parameters, $\rho_0$ and $T_0$ 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{equation}
F_{\rm rad}=-K\,{\rm d} {} T/{\rm d} {} z={\rm const}  {},
\label{FradEqn}
\end{equation}
where $K=16\sigma_{\rm SB} T^3/(3\kappa\rho)$ is the radiative conductivity
with $\sigma_{\rm SB}$ being the Stefan--Boltzmann constant, and $z$ is the
vertical coordinate in a Cartesian coordinate system.
For the simple opacity law \eq{Kramers}, but with
\begin{equation}
b<4+a,
\label{bpoly}
\end{equation}
the optically thick regime is characterized by a constant
temperature gradient and therefore $K={\rm const}  {}$ (see \App{Kconst}).
We have then a polytropic stratification with $\rho\propto T^n$, where
\begin{equation}
n=(3-b)/(1+a)
\label{n_Def}
\end{equation}
is the polytropic index.
It is larger than $-1$ when \Eq{bpoly} 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_{\rm 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$, so $n=3.25$, which corresponds to a
Schwarzschild-stable solution.
Only near the surface, at temperatures typically below $15,000\,{\rm K}$, the
dominant opacity is the ${\rm H}^{-}$ opacity.
It can no longer be approximated by a simple power law of the
type given by \Eq{Kramers}.
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 \eq{bpoly} is violated, so the
hydrostatic stratification is no longer polytropic,
but solutions can still be constructed numerically \citep{BB14} 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{equation}
\kappa^{-1}=\kappa^{-1}_{\rm Kr}+\kappa^{-1}_{{\rm H}^{-}},
\label{doubleKramers}
\end{equation}
where $\kappa_{\rm Kr}$ and $\kappa_{{\rm H}^{-}}$ are given by
\Eq{Kramers} with suitable exponents $a$ and $b$,
the solutions of the non-convective reference state is unstable only
over a depth of approximately $1\,{\rm Mm}$.
We return to this at the end of \Sec{NumSol},
where we present numerical solutions.

Of course, as stated earlier, the non-convective reference state
is only of academic interest.
Already with standard MLT \citep{Bie32,Vit53}, which allows for a
non-vanishing enthalpy flux, one finds a vastly extended convection zone
with a depth of the order of $100\,{\rm Mm}$ \citep{Bie38}.
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 \citep{CR16}.
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.

\section{The Deardorff flux}
\label{Deardorff}

\subsection{Derivation}
\label{Derivation}

In the meteorological context, counter-gradient heat flux terms have been
noticed for a long time \citep{Ert42,PS47,Dea66}.
They appear naturally when calculating the enthalpy flux $F_{\rm enth}$
using the $\tau$ approximation in its minimalistic form \citep[e.g.][]{BF03}.
In this approach, one computes the time derivative of $F_{\rm enth}$.
In the absence of ionization effects, $F_{\rm enth}$ can be written as
\begin{equation}
F_{\rm enth}=\overline{\rho u_z c_{P} T},
\label{Fconv_orig}
\end{equation}
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
\citep{SNGBS09}, 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 \Eq{Fconv_orig} and thus to fluctuations only in the correlation of
specific entropy and velocity.
Up to some reference value, we have
\begin{equation}
S/c_{P}=\ln T-\nabla_{\rm ad}\ln P,
\label{Seqn}
\end{equation}
where $\nabla_{\rm 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\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{\mbox{\boldmath $U$}}{}}{=\bm{0}$), so there are only velocity
fluctuations ($\mbox{\boldmath $U$} {}=\mbox{\boldmath $u$} {}$).
Focusing thus on the dominant contribution proportional to
$\overline{u_z s}$, we have
\begin{equation}
F_{\rm enth}=\overline{\rho}\, \overline{T}\,\overline{u_z s}.
\label{Fenth}
\end{equation}
Next, we write the time derivative of $F_{\rm enth}$ as
\begin{equation}
\partial F_{\rm enth}/\partial t=\overline{\rho}\, \overline{T}
\left(\,\overline{u_z \dot{s}}+\overline{\dot{u}_z s}\,\right),
\label{dFconv}
\end{equation}
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 \App{deriv10and11})
\begin{eqnarray}
&&\dot{s}=-u_j\nabla_j\overline{S}-s/\tau_{\rm cool}...,
\label{sdot}
\\
&&\dot{u}_i=-g_i s/c_{P}+...,
\label{udot}
\end{eqnarray}
where $g_i$ is the $i$th component of the gravitational acceleration
in Cartesian coordinates, namely $\mbox{\boldmath $g$} {}=(0,0,-g)$, the ellipses
refer to terms that are nonlinear in the fluctuations,
\begin{equation}
\tau_{\rm cool}^{-1}=\iota c_\gammak_{\rm f}
\;\;\mbox{with}\;\; \iota=\ellk_{\rm f}/(3+\ell^2k_{\rm f}^2)
\label{taucool}
\end{equation}
is the inverse heating and cooling time owing to radiation
\citep{US66,Edw90}, $c_\gamma=16\sigma_{\rm SB} \overline{T}^3/\overline{\rho} c_{P}$
is the photon diffusion speed \citep{BB14},
$\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 \Eq{sdot}.
We expect this to be the result of the nonlinear term
(indicated by the ellipses),
of which a part later gives rise to the $\nabla_{\rm D}$ term.
This will be motivated further at the end of this section
and in \Sec{DepthDependence}, where it will be included
in the final expression for $\nabla_{\rm D}$.

Using \Eqs{sdot}{udot} in \Eq{dFconv}, we have
\begin{equation}
{\partial F_i^{\rm enth}\over\partial t}\!=\!
\overline{\rho}\, \overline{T} \left(-\overline{u_i u_j}\,\nabla_j\overline{S}
-g_i \overline{s^2}/c_{P} \right)-{F_i^{\rm enth}\over\tau_{\rm cool}}
+{\cal T}_i,
\label{dotFconv}
\end{equation}
where the $\overline{s^2}$ term is the Deardorff flux and
the ${\cal T}_i$ refer to triple correlations that will be approximated
by the quadratic correlation $F_{\rm enth}$ as
\begin{equation}
{\cal T}_i=-F_i^{\rm enth}/\tau,
\label{MTAclosure}
\end{equation}
where $\tau$ is a relaxation time due to turbulence \citep{BF03},
which will later be identified with the turnover time.
This procedure, in which correlations with pressure fluctuations
are also neglected, is called the minimal $\tau$ approximation.
Important aspects of the closure assumption \eq{MTAclosure} have been
verified numerically for passive scalar transport \citep{BKM04}.
Among other things, they found that the time derivative on the
left-hand side of \Eq{dotFconv} 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_{\rm enth}$ can be neglected.
(Note, however, that this does not imply that the cooling term
will be neglected.)

Since the relaxation term in \Eq{MTAclosure} is similar to the cooling
term in \Eq{dotFconv}, we can combine the two by introducing 
the reduced relaxation time $\tau_{\rm red}$ defined through
\begin{equation}
\tau_{\rm red}^{-1}=\tau_{\rm cool}^{-1}+\tau^{-1}.
\label{tau}
\end{equation}
We can then solve for $F_{\rm enth}$, which appears on the
right-hand sides (rhs) of \Eqs{dotFconv}{MTAclosure}, and find
$\mbox{\boldmath $F$} {}_{\rm enth}=\mbox{\boldmath $F$} {}_{\rm G}+\mbox{\boldmath $F$} {}_{\rm D}$, where
\begin{eqnarray}
\label{FG}
\mbox{\boldmath $F$} {}_{\rm G}
&=&
-{\textstyle{1\over3}}\tau_{\rm red}u_{\rm rms}^2
\,\overline{\rho}\,\overline{T}\,\mbox{\boldmath $\nabla$} {}\overline{S},
\\
\label{FD}
\mbox{\boldmath $F$} {}_{\rm D}
&=&
-\tau_{\rm red}\overline{s^2}\,\mbox{\boldmath $g$} {}\,\overline{\rho}\,\overline{T}/c_{P}
\end{eqnarray}
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{\textstyle{1\over3}}\delta_{ij}u_{\rm rms}^2$,
where $u_{\rm rms}$ is the rms velocity of the turbulence.
\EEqs{FG}{FD} are written in vectorial forms to highlight the directions
of the fluxes: $\mbox{\boldmath $F$} {}_{\rm G}$ is counter-gradient and $\mbox{\boldmath $F$} {}_{\rm D}$
is countergravity.
In \Eq{FG}, the term ${\textstyle{1\over3}}\tau_{\rm red}u_{\rm 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{equation}
-{\rm d} {}(\overline{S}/c_{P})/{\rm d} {} z=(\nabla-\nabla_{\rm ad})/H_{\!P},
\end{equation}
where $\nabla={\rm d} {}\ln \overline{T}/{\rm d} {}\ln \overline{P}$ is the double-logarithmic
temperature gradient and $H_{\!P}=-({\rm d} {}\ln\overline{P}/{\rm d} {} z)^{-1}$ is the
pressure scale height.
Thus, we arrive at
\begin{equation}
F_{\rm enth}={\textstyle{1\over3}}\overline{\rho} c_{P} \overline{T} \, (\tau_{\rm red}u_{\rm rms}^2/H_{\!P})
\left(\nabla-\nabla_{\rm ad}+\nabla_{\rm D}\right),
\label{Fconv_urms2}
\end{equation}
where $\nabla_{\rm D}$ is a new contribution to standard MLT, which
results from $F_{\rm D}$.
Using the terms on the rhs of \Eq{FD}, we can write it explicitly as
\begin{equation}
\nabla_{\rm D}=(3/\gamma)\, (\overline{s^2}/c_{P}^2)\; \mbox{\rm Ma}^{-2},
\label{nabD}
\end{equation}
where $\mbox{\rm Ma}=u_{\rm 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_{\rm rms}^2$ factor in the definition \eq{Fconv_urms2}.
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 \Fig{sketch}.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{sketch}
\end{center}\caption[]{
Sketch illustrating overshoot in a stably stratified layer (a),
growth of perturbations in an unstably stratified layer (b),
buoyant rise of a blob with positive entropy perturbation (c), and
descent of a blob with negative entropy perturbation and hence
negative buoyancy (d).
In cases (b)--(d), the turbulent flux is upward (increasing $z$).
Case (d) is relevant to entropy rain.
}\label{sketch}\end{figure}

\subsection{Physical interpretation}
\label{Interpretation}

We recall that in a stably (unstably) stratified layer, shown in
\Fig{sketch}(a) and (b), and under the assumption of pressure equilibrium,
the entropy perturbation of a blob with respect to its current surroundings
decreases (increases) after a small ascent.
By contrast, in the marginally or nearly marginally stratified cases,
the perturbation remains constant if one ignores mixing
with the surroundings; see panels (c) and (d).
Therefore, a positive entropy
perturbation ($s>0$) always corresponds to a positive temperature
perturbation and a negative density perturbation.
Thus, the fluid parcel is buoyant and moves upward ($u_z>0$), so $u_zs>0$,
giving a positive contribution to $F_{\rm enth}$.
Likewise, for $s<0$, the parcel is cooler and heavier and sinks ($u_z<0$),
and again, $u_zs>0$, giving a positive contribution to $F_{\rm enth}$.

In laboratory convection, strong positive entropy perturbations can
be driven at the lower boundary, and strong negative ones at the top
boundary.
In the Sun, however, only the strong radiative cooling in the
photosphere provides a significant source of entropy perturbations.
This became clear with the emergence of the realistic solar convection
simulations of \cite{SN89}.
These simulations showed for the first time the strong vertical asymmetry
of solar convection resulting from the physics of ionization and strong
cooling via the ${\rm H}^{-}$ opacity.
This led to the notion of {\em entropy rain} and Spruit's description
of solar convection as a nonlocal phenomenon that is driven solely by
surface cooling.
Our considerations suggest that the resulting
(one-dimensional) mean-field energy flux is
parameterized not just in terms of the local superadiabatic gradient.

In analogy with laboratory convection \citep{HCL87,Castaing}, Spruit refers
to the downdrafts as threads.
Interestingly, the laboratory experiments show that these threads persist
even at large Rayleigh numbers.
For such threads to persist, one could imagine them to be stabilized
by their intrinsic vorticity, akin to a Hill vortex \citep{Hil94}.
These are vortex rings that can stay concentrated over large distances.
Vortex rings have also been reported by \cite{SNGBS09}, but many of
those are associated with mushroom shapes, indicating that they widen
and soon break up and mix with their surroundings.
Numerical studies (see \App{HillVortices}) suggest that Hill vortices
also survive in a highly stratified isothermal layer, and that their
persistence increases significantly with resolution and decreasing
viscosity.
On the other hand, the presence of background turbulence provides an
effective turbulent viscosity, which contributes to a premature
break up of small-scale vortices.
However, those studies were done with an isothermal equation of state
and thus ignore the effects of continued driving by a persistent negative
entropy perturbation between the vortex and its surroundings.

Studies of the Deardorff flux in the meteorological context suggest that
the flux-carrying plumes are associated with so-called coherent structures
\citep{DeRoode04} and those may not simply be the vortex-like structures
envisaged above.
It is clear, however, that the structures seen in the Earth's
atmosphere are driven by the boundary layer on the ground, and sometimes
also by the upper inversion layer \citep{DeRoode04}.
\cite{Ple07} discusses the Deardorff flux in connection with other
parameterizations of nonlocal fluxes such as the
transilient matrix approach \citep{Stu84,Stu93}.
The close connection between counter-gradient fluxes and nonlocal
transport has been elaborated upon by \cite{vDop01} and \cite{Buske07}.
Thus, while the approach presented in \Sec{Derivation} may capture the
physical phenomenon of the Deardorff flux qualitatively correctly, it
is quite possible that the nature of the underlying coherent structures
may require additional refinements.

\subsection{Depth dependence of the Deardorff term}
\label{DepthDependence}

To estimate the resulting depth dependence of the Deardorff term in
\Eq{nabD}, 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{equation}
f_s\propto\overline{\rho}^{\,-\zeta},
\label{fsDef}
\end{equation}
where $\zeta=0.8$ has been found for spherical vortex structures
descending in an isothermally stratified layer; see \App{HillVortices}.
However, \Eq{fsDef} 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
$\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 $\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 $\Delta S$.
We model this by considering $\zeta$ an adjustable parameter that is
{\em 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\,\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 {\em 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 {\em mean} entropy difference
$\Delta\overline{S}$ would be diluted correspondingly such that
\begin{equation}
f\,\Delta\overline{S}=f_s\,\Delta S_0\quad
\left(\propto\overline{\rho}^{\,-\zeta}\,\Delta S_0\right).
\label{fsDef}
\end{equation}
We recall that non-ideal effects can be captured by choosing
$\zeta>0.8$ in the prescription \eq{fsDef}.
The physics of entrainment has been modeled in the stellar context by
\cite{RZ95} and the extent of entrainment has been quantified in
the realistic surface simulations by \cite{TS11}, 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{equation}
\overline{S}=(1-f)\overline{S}_\uparrow+f\overline{S}_\downarrow=\overline{S}_\uparrow-f\,\Delta\overline{S},
\label{meanSfill}
\end{equation}
where $\overline{S}_\uparrow$ and $\overline{S}_\downarrow$
are the mean specific entropies in up- and downflows, and
$\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 \Eq{meanSfill}, gives
\begin{equation}
\overline{s^2}=(1-f)\,(\overline{S}_\uparrow-\overline{S})^2
+f\,(\overline{S}_\downarrow-\overline{S})^2=\hat{f}\,(\Delta\overline{S})^2,
\label{s2m}
\end{equation}
where $\hat{f}=(1-f)f$ has been introduced as a shorthand.
Analogously to the specific entropy, we find for the velocity
\begin{equation}
\overline{U}_z=(1-f)\overline{U}_\uparrow+f\overline{U}_\downarrow
=\overline{U}_\uparrow-f\,\Delta\overline{U},
\end{equation}
where $\overline{U}_\uparrow>0$ and $\overline{U}_\downarrow<0$ are the mean up- and
downflow velocities with $\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\approx0$ from mass conservation.
With this, we find
\begin{equation}
u_{\rm rms}^2\equiv\overline{u_z^2}=(1-f)\,\overline{U}_\uparrow^2
+f\,\overline{U}_\downarrow^2=\hat{f}\,(\Delta\overline{U})^2,
\label{u2m}
\end{equation}
and thus $\overline{U}_\uparrow/u_{\rm rms}=[f/(1-f)]^{1/2}$.
Similarly, we calculate
\begin{equation}
\overline{u_z^3}=(1-f)\,\overline{U}_\uparrow^3+f\,\overline{U}_\downarrow^3
=-\hat{f}\,(1-2f)\,(\Delta\overline{U})^3,
\end{equation}
which is negative for $f<1/2$, and finally
\begin{equation}
\overline{u_z s}=(1-f)\,\overline{U}_\uparrow\overline{S}_\uparrow
+f\,\overline{U}_\downarrow\overline{S}_\downarrow=\hat{f}\,\Delta\overline{U}\,\Delta\overline{S}.
\label{uzs_eqn}
\end{equation}
Thus, the kinetic energy flux, $\overline{\rho\mbox{\boldmath $u$} {}^2 u_z}/2$, which
reduces to $\overline{\rho}\,\overline{u_z^3}/2$ in this two-stream approximation,
is given by
\begin{equation}
F_{\rm kin}=-\phi_{\rm kin}\,\overline{\rho}u_{\rm rms}^3,
\label{Fkin}
\end{equation}
where $\phi_{\rm kin}=(1/2-f)/\hat{f}^{1/2}$ is a positive prefactor
(corresponding to downward kinetic energy flux) if $f<1/2$.
\cite{SNGBS09} find $f\approx1/3$, nearly independently of depth,
which yields $\phi_{\rm kin}\approx\sqrt{2}/4\approx0.35$; see
\Tab{SummaryPHIKIN}, where we list $\phi_{\rm kin}$ and
$-\overline{U}_\downarrow/u_{\rm rms}=[(1-f)/f]^{{1/2}}$ for selected values of $f$.
The enthalpy flux is proportional to $\overline{u_z s}$ and, using
\Eqs{Fenth}{uzs_eqn} together with \Eqs{s2m}{u2m}, we find
$F_{\rm enth}=\overline{\rho}\, \overline{T}\,u_{\rm rms}s_{\rm rms}$.
In \App{DeltaS_Ma2} we show that
$\overline{T}s_{\rm rms}=u_{\rm rms}^2k_{\rm f} H_{\!P}/(a_{\rm MLT}\nabla_{\rm ad})$, where $a_{\rm MLT}=1/8$
is a geometric factor in standard MLT.\footnote{Not to be confused with
the parameter $\alpha_{\rm mix}$ of \Sec{RelaxationTime} below.}
This leads to
\begin{equation}
F_{\rm enth}=\phi_{\rm enth}\,\overline{\rho}u_{\rm rms}^3
\label{Fenth2}
\end{equation}
with $\phi_{\rm enth}=\kfH_{\!P}/(a_{\rm MLT}\nabla_{\rm ad})$.
This yields $\phi_{\rm enth}\approx20$, which is rather large.
By contrast, \cite{BCNS05} determined a quantity $k_u$ such that
$\phi_{\rm enth}=k_u^{-3/2}\approx4$.

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 \Sec{Derivation}.
This is compatible with recent simulations of R.\ F.\ Stein
(2016, private communication), in which the fractional kinetic energy flux
increaes toward the deeper parts ($\ga40\,{\rm Mm}$) of the domain.

\begin{table}[t!]\caption{
$\phi_{\rm kin}$, $-\overline{U}_\downarrow/u_{\rm rms}$, and $\overline{U}_\uparrow/u_{\rm rms}$
for selected values of $f$.
}\centerline{\begin{tabular}{ccccccc}
$f$                        & 1/2 & 1/3  & 0.14  & 0.015 & 0.0006\\
\hline
$\phi_{\rm kin}$           &  0  & 0.35 &   1   &   4   &   20 \\
$-\overline{U}_\downarrow/u_{\rm rms}$ &  1  & 1.4  &  2.5  &   8   &   40 \\
$\overline{U}_\uparrow/u_{\rm rms}$    &  1  & 0.7  &  0.4  &  0.12 &  0.025
\label{SummaryPHIKIN}\end{tabular}}\end{table}

When $f$ becomes small ($<0.14$), $\phi_{\rm kin}$ exceeds unity and for
$f<0.015$, $\phi_{\rm kin}$ exceeds the estimate $\phi_{\rm enth}\approx4$
found by \cite{BCNS05}, 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
\cite{Hanasoge} 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_{\rm rms}$; for $f=0.015$ we have $-\overline{U}_\downarrow\approx8u_{\rm rms}$;
see \Tab{SummaryPHIKIN}.
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 \cite{Cat91}, 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 \citep{SBN92}.

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

In deeper layers, where the local value of $\nabla-\nabla_{\rm ad}$ has become
small, $\nabla_{\rm D}$ should scale with $\overline{s^2}=\hat{f}\,(\Delta\overline{S})^2$,
and therefore, using \Eq{fsDef}, it should be proportional to
$f_s^2(z)\propto\overline{\rho}^{\,-2\zeta}$, in addition to an $\mbox{\rm Ma}^{-2}$ factor;
see \Eq{nabD}.
Thus, we have
\begin{equation}
\nabla_{\rm D}=f_s^2\,\nabla_{\rm D}^{\max},
\label{nabD_ansatz}
\end{equation}
where we have used for $f_s$ the expression
\begin{equation}
f_s=f_{s0} (\overline{\rho}/\rho_\ast)^{-\tilde\zeta}
\quad\mbox{(for $\overline{\rho}>\rho_\ast)$}.
\label{fzeta}
\end{equation}
Here, $f_{s0}$ is a prefactor determining the strength of the
resulting Deardorff flux and $\rho_\ast$ is the density at the point
at which $\nabla-\nabla_{\rm ad}$ is equal to its maximum value
(which is just below the photosphere).
The new exponent $\tilde\zeta=\zeta-\Delta\zeta$ takes the scaling of
Mach number with density into account, i.e.,
\begin{equation}
\mbox{\rm Ma}\propto\overline{\rho}^{\,-\Delta\zeta},
\label{Madep}
\end{equation}
where $\Delta\zeta$ ($>0$) will be computed in \Sec{SuperadiabaticGradient}.

\subsection{Kinetic energy flux}

The importance of the kinetic energy flux has been stressed for some
time as a property of compressible stratified convection \citep{HTM84,Cat91}.
This flux is related to the asymmetry of up- and downflows \citep{SN89}
and is non-vanishing when $f\neq1/2$; see \Sec{DepthDependence}.
It is neglected in standard MLT, as has been discussed by \cite{Arnett15}.
Just like the Deardorff flux, the kinetic energy flux is also a
non-gradient flux, but it is always directed downward for $f<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_{\rm kin}$
are triple correlations, there is no $\tau$ approximation treatment
analogous to that of the Deardorff flux.
However, by comparing $F_{\rm kin}=-\phi_{\rm kin}\,\overline{\rho}u_{\rm rms}^3$
with the enthalpy flux in \Eq{Fconv_urms2}, it is possible to define
a corresponding nabla term via $F_{\rm kin}=-{\textstyle{1\over3}}\overline{\rho}
c_{P} \overline{T} \, (\tau_{\rm red}u_{\rm rms}^2/H_{\!P})\nabla_{\rm kin}$.
This yields
\begin{equation}
\nabla_{\rm kin}=3\,\phi_{\rm kin}\,u_{\rm rms} H_{\!P}/\tau_{\rm red}c_{P}\overline{T},
\end{equation}
which is obtained analogously to $\nabla_{\rm D}$ in \Eq{nabD}.
The prefactor is here defined with a positive sign, so
the total non-radiative flux caused by the turbulence is
proportional to $\nabla-\nabla_{\rm ad}+\nabla_{\rm D}-\nabla_{\rm kin}$.

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

\section{Modified mixing length model}
\label{Modified}

\subsection{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 ${\rm d} {}\overline{P}/{\rm d} {} z=-\overline{\rho} g$
and $K{\rm d} {}\overline{T}/{\rm d} {} z=-F$, or, alternatively for the logarithmic gradients,
\begin{eqnarray}
{\rm d} {}\ln \overline{P}/{\rm d} {} z&=&-\overline{\rho} g/\overline{P},
\label{dlnP}
\\
{\rm d} {}\ln \overline{T}/{\rm d} {} z&=&-F/(K\overline{T}).
\label{dlnT}
\end{eqnarray}
The double-logarithmic temperature gradient is obtained by dividing
the two equations through each other, i.e.,
\begin{equation}
\nabla={{\rm d} {}\ln \overline{T}\over{\rm d} {}\ln \overline{P}}
={F\overline{P}\over K\overline{T}\,\overline{\rho} g}={Fc_{P}\nabla_{\rm ad} \over Kg},
\label{nabla}
\end{equation}
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}\nabla_{\rm ad}$.
If the energy is no longer carried by radiation alone, $\nabla$ cannot be
computed from \Eq{nabla}, but we have to invoke a suitable theory of
convection.
In standard MLT, one obtains $\nabla$ as a solution of a cubic equation
\citep{Vit53,KW90}.
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{equation}
F_{\rm tot}=F_{\rm rad}+F_{\rm conv}.
\label{FluxBalance}
\end{equation}
where $F_{\rm conv}=F_{\rm enth}+F_{\rm 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 $\nabla$,
the radiative flux is evidently
\begin{equation}
F_{\rm rad}={Kg\over c_{P}\nabla_{\rm ad}}\,\nabla,
\end{equation}
where $\nabla$ characterizes the {\em actual} temperature gradient.
We can also define a {\em hypothetical} radiative temperature gradient
$\nabla_{\rm rad}$ that would result if all the energy were carried by radiation,
so we can write 
\begin{equation}
F_{\rm tot}={Kg\over c_{P}\nabla_{\rm ad}}\,\nabla_{\rm rad},
\label{nabrad_Def}
\end{equation}
which follows from \Eq{nabla}.
We note in passing that the ratio between $F_{\rm 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
\citep{HTM84}, whose local value is thus equal to
\begin{equation}
\mbox{Nu}=\nabla_{\rm rad}/\nabla_{\rm ad}.
\label{Nu}
\end{equation}
Finally, as explained in \Sec{Deardorff}, we have from \Eq{Fconv_urms2}
\begin{equation}
F_{\rm enth}=F_0\left(\nabla-\nabla_{\rm ad}+\nabla_{\rm D}\right),
\label{Fconv_urms1}
\end{equation}
with $F_0={\textstyle{1\over3}}\overline{\rho} c_{P} \overline{T} \tau_{\rm red}u_{\rm rms}^2/H_{\!P}$,
and $F_{\rm kin}$ being essentially proportional to $F_{\rm enth}$;
cf.\ \Eqs{Fkin}{Fenth2}.
Flux equilibrium then implies
\begin{equation}
\nabla_{\rm rad}=\nabla+\epsilon\left(\nabla-\nabla_{\rm ad}+\nabla_{\rm D}\right),
\label{epsilon_eqn}
\end{equation}
where $\epsilon=F_0 c_{P}\nabla_{\rm ad}(1-\phi_{\rm kin}/\phi_{\rm enth})/Kg$.
However, the expression for $F_0$ involves the still unknown
values of $u_{\rm rms}$ and $\tau_{\rm red}$, which will be discussed in
\Sec{RelaxationTime}.

\begin{figure*}[t!]\begin{center}
\includegraphics[width=\textwidth]{pbetacases2}
\end{center}\caption[]{
Illustrative flow structures (upper row) and corresponding
horizontal power spectra (lower row) associated with the three
combinations of $\beta$ and ${\tilde{\beta}}$ considered in this paper.
Light shades correspond to large logarithmic power, which is seen to
extend over large values of $k$ for $\beta=0$ (Cases~I and II)
and is confined to progressively smaller $k$ at larger depths
when $\beta=1$ (Case~III).
}\label{pbetacases2}\end{figure*}

\subsection{Relaxation time and mixing length}
\label{RelaxationTime}

In convection the relevant time scale is the turnover time,
which we write as $\tau=1/\urmsk_{\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 \citep{SN89,Spr97},
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{\rm const}  {}$.
To capture the various cases in one expression, we assume in the following
\begin{equation}
k_{\rm f}=\alpha_{\rm mix}^{\beta}H_{\!P}^{-1} (\kfzH_{\!P})^{1-\beta},
\label{kfbeta}
\end{equation}
where $\beta=0$ corresponds to $k_{\rm f}=k_{\rm f0}$ with a fixed value $k_{\rm f0}$
of the wavenumber of the entropy rain and $\beta=1$ corresponds to
$\kfH_{\!P}=\alpha_{\rm mix}$, where we have allowed for the possibility
of a free mixing length parameter, $\alpha_{\rm mix}$, as is commonly
done in standard MLT.
It is not to be confused with the parameter $a_{\rm MLT}$ that
was introduced in \Sec{DepthDependence} and \App{DeltaS_Ma2}.
Negative values of $\beta$ would correspond to shrinking scales,
but will not be considered here.
Likewise, non-integer values of $\beta$ are conceivable,
but will also not be considered here.

Returning now to the discussion of the relaxation time $\tau$ in the beginning
of this subsection, instead of associating it with the turnover time
$(\urmsk_{\rm f})^{-1}$, we will allow for the possibility of an additional
dilution factor $\phi(z)$ and write $\tau=(\urmsk_{\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
\Eq{kfbeta}, the expression $\phi=(\kfzH_{\!P}/\alpha_{\rm mix})^{\beta-{\tilde{\beta}}}$.
This dilution factor only enters in expressions involving the turbulent
diffusivity, such as in \Eq{FG}, and hence terms involving $\tau$ or $\tau_{\rm red}$.
For $\beta=1$, the case ${\tilde{\beta}}=1$ corresponds to the usual MLT concept,
while ${\tilde{\beta}}=0$ corresponds to a value of $\phi$ that increases with depth.

In \Fig{pbetacases2} we present illustrative flow structures
as well as their depth-dependent horizontal power spectra associated
with the three combinations of $\beta$ 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.

With these preparations in place, we can write $F_0$ in \Eq{Fconv_urms1}
in the form
\begin{equation}
F_0=(\sigma/3\gamma\nabla_{\rm ad}) \, \overline{\rho} u_{\rm rms} c_{\rm s}^2
\alpha_{\rm mix}^{-{\tilde{\beta}}} (k_{\rm f0} H_{\!P})^{-(1-{\tilde{\beta}})},
\label{F0}
\end{equation}
where
\begin{equation}
\sigma\equiv\tau_{\rm red}\urmsk_{\rm f}\phi=u_{\rm rms}/(u_{\rm rms}+\iota c_\gamma/\phi)
\label{urms}
\end{equation}
quantifies the radiative heat exchange between convective elements
and the surroundings.
The $\iota$ term defined in \Eq{taucool}, has a maximum at
$\ellk_{\rm f}=\sqrt{3}$, which typically occurs near the surface
\citep[e.g.][]{BB14}.

In \Eq{F0}, the local value of the turbulent rms velocity always depends
on the actual flux transported, and therefore it must also depend on
$\nabla-\nabla_{\rm ad}+\nabla_{\rm D}$.
On dimensional grounds, since $F_{\rm enth}$ is proportional to $u_{\rm rms}^3$,
we have
\begin{equation}
u_{\rm rms}=c_0\left(\nabla-\nabla_{\rm ad}+\nabla_{\rm D}\right)^{1/2},
\label{prefactor}
\end{equation}
so $F_{\rm enth}$ and $F_{\rm kin}$ are proportional to
$(\nabla-\nabla_{\rm ad}+\nabla_{\rm D})^{3/2}$.
Standard mixing length arguments can be used to show that
the prefactor in \Eq{prefactor} is a fraction of $c_{\rm s}$.
As shown in \App{App}, we have
\begin{equation}
c_0/c_{\rm s}=\sqrt{\sigma a_{\rm MLT}/3\gamma}\,
\alpha_{\rm mix}^{-\beta'} (k_{\rm f0} H_{\!P})^{-(1-\beta')},
\label{urms0}
\end{equation}
where $\beta'=(\beta+{\tilde{\beta}})/2$.
It is then clear that $F_{\rm kin}$ scales with $H_{\!P}$ like
$(k_{\rm f0} H_{\!P})^{-3(1-\beta')}$.
This is not the case for $F_{\rm enth}$, however, because of
the factor $(k_{\rm f0} H_{\!P})^{1-{\tilde{\beta}}}$ which, together with the
$(k_{\rm f0} H_{\!P})^{1-\beta'}$ factor, gives the scaling proportional
to $(k_{\rm f0} H_{\!P})^{2(1-\beta'')}$, where $\beta''=(\beta'+{\tilde{\beta}})/2$.

\subsection{Equation for the superadiabatic gradient}
\label{SuperadiabaticGradient}

Now that we know $F_0$, we can solve the equation for the superadiabatic
gradient.
This leads to an equation that is similar to the cubic equation for $\nabla$,
which is familiar from standard MLT \citep{KW90},
\begin{equation}
\nabla_{\rm rad}=\nabla+\epsilon_\ast\left(\nabla-\nabla_{\rm ad}+\nabla_{\rm D}\right)^\xi,
\label{nabla_final}
\end{equation}
where $\xi=3/2$, and $\epsilon_\ast=(1-\phi_{\rm kin}/\phi_{\rm enth})
(\epsilon_{\rm enth}+\epsilon_{\rm kin})$
has contributions from the enthalpy and kinetic energy fluxes.
These expressions are similar to $\epsilon$ in \Eq{epsilon_eqn}, except that
$\epsilon_{\rm enth}$ is evaluated with $c_0$ in place of $u_{\rm rms}$, i.e.,
\begin{equation}
\epsilon_{\rm enth}=(\sigma/3\gamma)^{3/2}\,a_{\rm MLT}^{1/2}\;(c_{\rm s}^3/\chi g)
\,\alpha_{\rm mix}^{-2\beta''} (\kfzH_{\!P})^{-2(1-\beta'')},
\label{eps2}
\end{equation}
where $\chi=K/\overline{\rho} c_{P}$ is the radiative diffusivity.
This shows that $\epsilon_\ast$ is essentially a P\'eclet number
based on $c_{\rm s}$.
The contribution from the kinetic energy flux is given by
\begin{equation}
\epsilon_{\rm kin}/\epsilon_{\rm enth}=-\phi_{\rm kin}
a_{\rm MLT}\nabla_{\rm ad} \, \alpha_{\rm mix}^{-\beta} (\kfzH_{\!P})^{-(1-\beta)}.
\label{eps3}
\end{equation}

We note in passing that $\epsilon_\ast$ is also related to the
Rayleigh number, which is commonly defined in laboratory and numerical
studies of convection; see \App{Rayleigh}.
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={\rm const}  {}$
(\Sec{HighlyUnstable}), $K$ can no longer be constant,
but it reaches a minimum at the point where
$\kappa$ is maximum, which turns out to be at a depth of about $1\,{\rm Mm}$
in the convection models presented below.
Since $\epsilon_\ast$ 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 \Eq{nabla_final} and the usual one in
MLT is the presence of $\nabla_{\rm D}$ arising from the Deardorff flux.
Within the usual MLT, where $\nabla_{\rm D}=0$, one finds that $\nabla$ is slightly
above $\nabla_{\rm ad}$, but now it might instead be slightly above
$\nabla_{\rm ad}-\nabla_{\rm D}$.
There are indeed two possibilities for convecting solutions
($F_{\rm conv}>0$), one corresponding to a Schwarzschild-stable solution,
\begin{equation}
\nabla_{\rm ad}-\nabla_{\rm D} < \nabla < \nabla_{\rm ad} \quad\mbox{(stable)},
\label{stable}
\end{equation}
and one that is Schwarzschild unstable,
\begin{equation}
\nabla_{\rm ad}-\nabla_{\rm D} < \nabla_{\rm ad} < \nabla \quad\mbox{(unstable)}.
\label{unstable}
\end{equation}
Which of the two possibilities is attained depends on the value of
$\nabla_{\rm D}$ and also on details of the solution.
As will be discussed in \Sec{UnstablyStratified} 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 $\nabla-\nabla_{\rm ad}$, the effective buoyancy force is written as
$\nabla-\nabla'$, where $\nabla'$ always lies between $\nabla$
and $\nabla_{\rm ad}$ \citep{Vit53}.
Thus, one has $\nabla_{\rm ad} < \nabla' < \nabla$, which resembles
\Eq{unstable} with a negative value of $\nabla_{\rm D}$.

To understand the nature of the solutions of \Eq{nabla_final},
it is instructive to treat $\xi$ as an adjustable parameter.
For given values of $\nabla_{\rm rad}$ and $\widetilde\nabla_{\rm ad}\equiv\nabla_{\rm ad}-\nabla_{\rm D}$,
the case $\xi=1$ yields
\begin{equation}
\nabla(\epsilon_\ast)={\nabla_{\rm rad}+\epsilon_\ast\widetilde\nabla_{\rm ad}\over1+\epsilon_\ast}.
\label{xi1}
\end{equation}
It shows that $\nabla\to\nabla_{\rm rad}$ for $\epsilon_\ast\to0$ (stable surface
layers) and $\nabla\to\widetilde\nabla_{\rm ad}$ for $\epsilon_\ast\gg1$ (deeper layers).
Next, to discuss the general case $\xi\neq1$, we define
$\Delta\nabla=\nabla-\widetilde\nabla_{\rm ad}$ and $\Delta\nabla_{\rm rad}=\nabla_{\rm rad}-\widetilde\nabla_{\rm ad}$.
For $\Delta\nabla_{\rm rad}<0$ we have $\Delta\nabla=\Delta\nabla_{\rm rad}$, while for
$\Delta\nabla_{\rm rad}>0$ and $\Delta\nabla\ll1$, a useful approximation is
\begin{equation}
\Delta\nabla\approx\Delta\nabla_{\rm rad}^{1/\xi}\left/\left(q
\Delta\nabla_{\rm rad}^{1/\xi-1}+\epsilon_\ast^{1/\xi}\right)\right.,
\label{DeltaNabla}
\end{equation}
where $q=\xi^{-1}$.
It agrees with \Eq{xi1} in the special case $\xi=1$, where
$\Delta\nabla=\Delta\nabla_{\rm rad}/(1+\epsilon_\ast)$.
In \Fig{pcubic_single} we plot $\nabla$ versus $\epsilon_\ast$
for $\nabla_{\rm rad}=10^5$.
The approximation yields $\nabla>\nabla_{\rm rad}$ for $\epsilon_\ast<10^{-3}$,
which is unphysical.
This can be mitigated by choosing $q=1$; see \Fig{pcubic_single}.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pcubic_single}
\end{center}\caption[]{
Solution of \Eq{nabla_final} for $\xi=3/2$ and $\nabla_{\rm rad}=10^5$
(solid black), compared with the approximation \eq{DeltaNabla}
for $q=\xi^{-1}$ (dashed blue) and $q=1$ (dash--dotted red),
and $\xi=1$ (thin orange).
The dotted line gives an additional unphysical solution
of \Eq{nabla_final} for $\xi=3/2$.
The limiting cases $\nabla=\nabla_{\rm ad}$ and $\nabla_{\rm rad}$ are shown as
thin horizontal green lines.
}\label{pcubic_single}\end{figure}

For the relevant case of large values of $\epsilon_\ast$, we have
\begin{equation}
\Delta\nabla\approx(\nabla_{\rm rad}/\epsilon_\ast)^{1/\xi}.
\label{DelNabApprox}
\end{equation}
This relation is useful because, even though both
$\nabla_{\rm rad}$ and $\epsilon_\ast$ depend on $K$, their ratio does not
and is given by
\begin{equation}
{\nabla_{\rm rad}\over\epsilon_\ast}
={3 \nabla_{\rm ad} \over \sigma} \, {F_{\rm tot}/\overline{\rho} \over c_0 g/k_{\rm f}}
\propto {(k_{\rm f0} H_{\!P})^{2-2\beta''} \over \overline{\rho} c_{\rm s}^3 }
\propto\overline{T}^{\,-(1+2\beta'')},
\label{epsratio}
\end{equation}
where $\beta''=(\beta+3{\tilde{\beta}})/4$.
Thus, $\Delta\nabla\propto\overline{T}^{\,-(1+2\beta'')/\xi}$ gives the
scaling of $\Delta\nabla$ for the bulk of the convection zone.
Therefore, looking at \Eq{prefactor}, we find
$\mbox{\rm Ma}\propto\overline{T}^{\,-m}$ with
\begin{equation}
m=1-\beta'+(1/2+\beta'')/\xi.
\label{mexpression}
\end{equation}
For $\xi=3/2$, using the relation $3\beta'-2\beta''=\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 \Eq{Madep}, the following scaling: $\Delta\zeta=8/9$ with $\beta=0$
and $\Delta\zeta=2/3$ with $\beta=1$.

\subsection{Relative importance of the Deardorff term}
\label{RelativeImportance}

In \Sec{DepthDependence} we have considered the depth dependence
of the Deardorff term via \Eqs{nabD_ansatz}{fzeta}.
However, for $\nabla_{\rm D}$ to be important at increasing depths, it must
exceed the {\em sub}adiabatic gradient $\nabla_{\rm ad}-\nabla$, because otherwise
it would not be possible for the Deardorff term to make
$\nabla-\nabla_{\rm ad}+\nabla_{\rm D}$ positive.
From \Eqs{DelNabApprox}{epsratio} we see that $\Delta\nabla$ depends
on $\overline{T}$ in a power-law fashion.
We are particularly interested in the conditions under which
$\Delta\nabla$ falls off faster than $\nabla_{\rm D}$, because that would
ensure that the $\nabla_{\rm D}$ term remains important even at larger depths.

\begin{table}[t!]\caption{
Comparison of $\zeta_{\Delta\nabla}$, $2\tilde\zeta=2\zeta-2\Delta\zeta$,
and their difference for various combinations of $\beta$ and ${\tilde{\beta}}$ 
using $\zeta=1$.
}\centerline{\begin{tabular}{lccccccc}
& $\beta$ & ${\tilde{\beta}}$ & $\beta'$ & $\beta''$ &
$\zeta_{\Delta\nabla}$
& $2\tilde\zeta$ & $\!\!\zeta_{\Delta\nabla}-2\tilde\zeta$ \\
\hline
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 \\
\label{SummaryBETA}\end{tabular}}\end{table}

For $\xi=3/2$, and since the convection zone is nearly isentropically
stratified ($\overline{T}\propto\overline{\rho}^{\,2/3}$), we have
\begin{equation}
\Delta\nabla\propto\overline{\rho}^{\,-\zeta_{\Delta\nabla}}
\quad\mbox{with}\quad \zeta_{\Delta\nabla}={4\over9}(1+2\beta'').
\end{equation}
In \Tab{SummaryBETA} we compare for various combinations of $\beta$
and ${\tilde{\beta}}$ the exponents for $\Delta\nabla$ and $\nabla_{\rm D}$, i.e.,
$2\tilde\zeta=2\zeta$ and $2\Delta\zeta$, where $\Delta\zeta=2m/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 $\beta$.
For better illustration, we summarize in \Tab{Summary} the various
relationships that led to the scaling of $\Delta\nabla$ with $\kfzH_{\!P}$.

In the first expression for $F_{\rm enth}$, $\beta$ enters because it
characterizes the relation between the buoyancy force proportional to
$\delta T/\overline{T}$ and advection proportional to $k_{\rm f}u_{\rm rms}^2$; see
\App{DeltaS_Ma2}.
For thin threads, we expect the relevant $k_{\rm f}$ to be large, i.e., $\beta=0$
(Cases~I and II in \Fig{pbetacases2}).
Next, in the second expression for $F_{\rm enth}$, we have used a mean-field
expression to relate $F_{\rm enth}$ to $\Delta\nabla$ via a turbulent
diffusivity proportional to $\tauu_{\rm rms}^2\approxu_{\rm rms}/k_{\rm f}\phi$, where,
as argued above, the dilution factor $\phi$ has entered.
It is this expression that is closest to conventional MLT, because here
we expect ${\tilde{\beta}}=1$.

\begin{table}[t!]\caption{
Illustration of scaling relationships with $\kfzH_{\!P}$.
}\centerline{\begin{tabular}{ll}
\hline
\hline
$F_{\rm enth}\proptou_{\rm rms}^3\,(\kfzH_{\!P})^{1-\beta}$ &
  $\!\!$buoyancy force\\
$F_{\rm enth}\proptou_{\rm rms}\Delta\nabla\,/\,(\kfzH_{\!P})^{1-{\tilde{\beta}}}$   &
  $\!\!$mean-field expression \\
$u_{\rm rms}\;   \propto(\Delta\nabla)^{1/2}/(\kfzH_{\!P})^{1-\beta'}$ &
  $\!\!$$\beta'\;=(\beta\!+\!{\tilde{\beta}})/2$ \\
$F_{\rm enth}\propto(\Delta\nabla)^{3/2}/(\kfzH_{\!P})^{2(1-\beta'')}$ &
  $\!\!$$\beta''=(\beta'\!+\!{\tilde{\beta}})/2=(\beta\!+\!3{\tilde{\beta}})/4 $ \\
$F_{\rm kin}\;\;\propto(\Delta\nabla)^{3/2}/(\kfzH_{\!P})^{3(1-\beta')}$ &
  $\!\!$kinetic energy flux \\
\hline
\label{Summary}\end{tabular}}\end{table}

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

\section{Numerical solutions}
\label{NumSol}

In this section we present numerical solutions to demonstrate the effect
of the $\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}={\cal R}/\mu$, where ${\cal R}$ is the universal gas constant).
A rather simple opacity law of the form of \Eq{doubleKramers} with
$\kappa_0=10^4\,{\rm cm}^2\,{\rm g}^{-1}$, $\rho_0=10^{-5}\,{\rm g}\,{\rm cm}^{-3}$, $T_0=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 \eq{dlnP} and \eq{dlnT} decouples
by using $\ln P$ as the independent variable.
We thus integrate
\begin{equation}
{\rm d} {}\ln \overline{T}/{\rm d} {}\ln\overline{P}=\nabla,
\label{dlnT_final}
\end{equation}
using Equations \eq{nabD_ansatz}, \eq{fzeta}, and \eq{nabla_final}
to compute $\nabla$.
As an initial condition we use $\overline{T}_{\rm top}=T_{\rm eff}/2^{1/4}$
at a sufficiently low pressure (here $\overline{P}_{\rm top}=10^5\,{\rm dyn}\,{\rm cm}^{-2}$),
so as to capture the initially isothermal part of the atmosphere; see,
e.g., \cite{Vit58} or \cite{Mih78}.
Here, $T_{\rm eff}=(F_{\rm tot}/\sigma_{\rm SB})^{1/4}$ is the effective temperature.
We approximate the $u_{\rm rms}$ term in \Eq{urms} 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-\Delta\zeta$ (\Sec{DepthDependence}),
where $\zeta=1$ and $\Delta\zeta$ is a function of $\beta$ and ${\tilde{\beta}}$
(\Sec{SuperadiabaticGradient}).
Since we integrate from the top downward, no prior knowledge of
$\rho_\ast$ is needed, because the Deardorff term is invoked only
after $\nabla-\nabla_{\rm ad}$ has reached its peak value.

Geometrical and optical depths are obtained respectively as
\begin{equation}
-z=\int(\overline{P}/\overline{\rho} g)\,{\rm d} {}\ln \overline{P}
\quad\mbox{and}\;\;
\tau=\int(\kappa \overline{P}/g)\,{\rm d} {}\ln \overline{P},
\end{equation}
where the integration of $\tau$ starts at $\ln\overline{P}_{\rm 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 $\tau$.
The $z$ coordinate is used in some of our plots.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pcomp5_lnppp_kfHp00}
\end{center}\caption[]{
Profiles of $\overline{S}/c_{P}$, $\nabla-\nabla_{\rm ad}$, and $u_{\rm rms}$ for $f_{s0}=0$
($\nabla_{\rm D}=0$), as well as $f_{s0}=0.2$, $0.3$, and $0.5$ for
$\beta={\tilde{\beta}}=0$ (Case~I) with $\tilde\zeta=1/9$.
The location of the surface ($\tau=1$) is indicated by vertical
dash-dotted lines and geometric depths below the surface are indicated
in the middle panel, starting with tick marks at $100$, $200$, and
$500\,{\rm km}$, and continuing with $1$, $2$, and $5\,{\rm Mm}$, etc.
The inset in the middle panel shows $\nabla-\nabla_{\rm ad}$ over a narrower
range as a function of $-z$.
}\label{pcomp5_lnppp_kfHp00}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pcomp5_lnppp_kfHp01}
\end{center}\caption[]{
Same as \Fig{pcomp5_lnppp_kfHp00}, but for $\beta=0$, ${\tilde{\beta}}=1$
(Case~II) with $\tilde\zeta=1/9$, and $f_{s0}=0$, $0.1$, and $0.2$.
}\label{pcomp5_lnppp_kfHp01}\end{figure}

In the following, we present solutions for the three combinations
of $\beta$ and ${\tilde{\beta}}$ sketched in \Fig{pbetacases2}.
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 \Fig{pbetacases2}) and are therefore
of interest when trying to reconcile the non-detection of convective
motions by \cite{Hanasoge,HGS16} 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}\ge0.5$)
can make the resulting stratification sufficiently subadiabatic.
This is the case shown in \Fig{pcomp5_lnppp_kfHp00}, where we present
profiles of $\overline{S}/c_{P}$, $\nabla-\nabla_{\rm ad}$, and $u_{\rm rms}$ for $f_{s0}=0$
(no Deardorff term), as well as $f_{s0}=0.2$, $0.3$, and $0.5$.
The $u_{\rm 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_{\rm rms}\propto\overline{T}^{1/2-m}$ and $m=4/3$;
see \Eq{mexpression}, where we have used $c_{\rm s}\propto\overline{T}^{1/2}$.
Thus, for isentropic stratification we have
$u_{\rm rms}\propto\overline{P}^{(1-2m)/5}=\overline{P}^{\,-1/3}$.

Next, we show in \Figs{pcomp5_lnppp_kfHp01}{pcomp5_lnppp_kfHp} 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_{\rm rms}$ are similar regardless of the
value of $f_{s0}$, but fall off more slowly for $\beta=1$
($u_{\rm rms}\propto\overline{P}^{\,-1/5}$), compared to $\beta=0$
($u_{\rm rms}\propto\overline{P}^{\,-1/3}$).
This agrees with our theory, because for $m=1$ we have
$u_{\rm rms}\propto\overline{P}^{(1-2m)/5}=\overline{P}^{\,-1/5}$.

Given that the stratification is stable in the deeper parts, we can
calculate the Brunt-V\"ais\"al\"a frequency of buoyancy oscillations,
$N_{\rm BV}$, which is given by $N_{\rm BV}^2=-(\nabla-\nabla_{\rm ad})g/H_{\!P}$.
At intermediate and larger depths of the convection zone, we have
$g/H_{\!P}\la10^{-2}\,{\rm s}^{-2}$ and $\nabla_{\rm ad}-\nabla\la10^{-4}$, so the
period of buoyancy oscillations would be of the order of days.
This is comparable to or less than the turnover time $\tau$.
Indeed, one finds that
$\tau N_{\rm BV}\approx\sqrt{\nabla_{\rm D}/\Delta\nabla-1}$
exceeds unity in the deeper parts, which is a consequence of $\nabla_{\rm D}$
falling off with $\overline{\rho}$ with a smaller power than $\Delta\nabla$,
as discussed in \Sec{RelativeImportance}.
Nevertheless, the resulting decrease in $\overline{S}/c_{P}$ with
depth remains always small ($10^{-3}$ or less) compared
with the value of $\Delta S_0/c_{P}$ produced at the surface
($\approxu_{\rm rms}^2/a_{\rm MLT}\nabla_{\rm ad}c_{\rm s}^2\approx0.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_{\rm conv}=0$.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pcomp5_lnppp_kfHp}
\end{center}\caption[]{
Same as \Fig{pcomp5_lnppp_kfHp00}, but for $\beta={\tilde{\beta}}=1$ (Case~III)
$\tilde\zeta=1/3$, and $f_{s0}=0$, $0.1$, and $0.2$.
}\label{pcomp5_lnppp_kfHp}\end{figure}

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 \Sec{Intro}, 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 \cite{Ser09} 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 \citep{Asp04} which is based on three-dimensional
convective atmosphere simulations, compared to previous analysis
based on one-dimensional semi-empirical models \citep{GS98}.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pcomp5_rad_ppp}
\end{center}\caption[]{
Comparison of $\overline{S}/c_{P}$ (top) and $\nabla_{\rm rad}/\nabla_{\rm ad}$ (bottom)
between the non-convective radiative reference solution ($\nabla=\nabla_{\rm rad}$)
and standard convective solutions ($\nabla_{\rm D}=0$) with
$\beta={\tilde{\beta}}=0$ (red, dashed) and $\beta={\tilde{\beta}}=1$ (blue, dotted).
In both panels, the location of the $\tau=1$ surface is indicated by
vertical dash-dotted lines and geometric depths below the surface are
indicated for the non-convective solution, starting with tick marks at
$50$, $100$, and $200\,{\rm km}$, etc.
The location of the a priori unstable layer ($\nabla_{\rm rad}>\nabla_{\rm ad}$) is marked
in the lower panel by a small gray strip.
For the other solutions, the depths are different; see those
in \Fig{pcomp5_lnppp_kfHp} for the solution with $f_{s0}=0$, which is
the same as here for $\beta={\tilde{\beta}}=1$.
}\label{pcomp5_rad_ppp}\end{figure}

Finally, we compare in \Fig{pcomp5_rad_ppp} the standard convective
solution ($\beta={\tilde{\beta}}=1$; same as the case with $f_{s0}=0$ in
\Fig{pcomp5_lnppp_kfHp}) 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\,{\rm Mm}$) and the peak in 
$\nabla_{\rm rad}/\nabla_{\rm ad}$ is about $100\,{\rm km}$ below the $\tau=1$ surface.
Note also that its peak value ($\approx10^5$) is below that in the
presence of convection, where it reaches a maximum of $\approx4\times10^6$
at a depth of $\approx1\,{\rm Mm}$.
As shown in \Eq{Nu}, this value can be interpreted as the local Nusselt number.
Note also that the result for $\beta={\tilde{\beta}}=0$ (when $\nabla_{\rm D}$ is weak)
is rather similar to that for $\beta={\tilde{\beta}}=1$; see the dashed and
dotted lines in \Fig{pcomp5_rad_ppp}.

Our calculations have demonstrated that for ${\tilde{\beta}}=1$, regardless of the value
of $\beta$, 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.

\section{Alternative considerations}
\label{UnstablyStratified}

Assuming that the upflows are perfectly isentropic, \cite{Spr97} 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 \Eq{meanSfill}, one obtains
\begin{equation}
-H_{\!P}{{\rm d} {}\overline{S}/c_{P}\over{\rm d} {} z}
={(\Delta S)_0\overc_{P}}{{\rm d} {} f_s\over{\rm d} {}\ln\overline{P}}
={2\over5}{(\Delta S)_0\overc_{P}}f_s > 0,
\end{equation}
where we have used ${\rm d} {}\ln f_s/{\rm d} {}\ln\overline{\rho}=2/3$
and ${\rm d} {}\ln\overline{\rho}/{\rm d} {}\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{equation}
\overline{\rho}\,\overline{T}\,\overline{U}_{\uparrow}\,{\rm d} {}\overline{S}_\uparrow/{\rm d} {} z
\approx-({\rm d} {} F_{\rm rad}/{\rm d} {} z)_\uparrow.
\label{dSdz_up}
\end{equation}
The sign of ${\rm d} {} F_{\rm rad}/{\rm d} {} z$ is negative and thus compatible
with a positive ${\rm d} {}\overline{S}_\uparrow/{\rm d} {} z$ in the upflows, but it
would only be large enough near the bottom of the convection zone.
Higher up, ${\rm d} {} F_{\rm rad}/{\rm d} {} z$ becomes smaller and eventually unimportant.
On the other hand, $\overline{U}_\uparrow/u_{\rm rms}$ can be rather small
(see \Tab{SummaryPHIKIN}).
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.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{p2}
\end{center}\caption[]{
Dependence of $\overline{U}^\uparrow_{\rm eff}/u_{\rm rms}$ on depths $-z$
for Case~I with $f_{s0}=0.5$ (solid red line), as well as Cases~II and III
with $f_{s0}=0.1$ (dashed blue and dotted green lines, respectively).
}\label{p2}\end{figure}

We now use selected solutions obtained in \Sec{NumSol} to estimate
the effective fractional upward velocity for radiative heating/cooling
to dominate over advection, defined as
\begin{equation}
\overline{U}^\uparrow_{\rm eff}/u_{\rm rms}=
(-{\rm d} {} F_{\rm rad}/{\rm d} {} z)\,/\,(u_{\rm rms}\,\overline{\rho}\,\overline{T}{\rm d} {}\overline{S}/{\rm d} {} z).
\end{equation}
Here we have used $\overline{T}{\rm d} {}\overline{S}/{\rm d} {} z=g\,\Delta\nabla/\nabla_{\rm ad}$.
\Fig{p2} shows this quantity for Cases~I--III, which are shown in
\Figss{pcomp5_lnppp_kfHp00}{pcomp5_lnppp_kfHp}.
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_{\rm rms}$
values in \Tab{SummaryPHIKIN} if we assume $f=0.01$
(so $\overline{U}^\uparrow\!/u_{\rm rms}=0.1$) and
$\overline{U}^\uparrow_{\rm eff}/\overline{U}^\uparrow=0.1$.
For Case~II with $f_{s0}=0.1$, the effective fractional upward velocity is
around $10^{-3}$, which may be unrealistically small.

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 \Eq{dSdz_up} 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 \Eq{FluxBalance}; see also \cite{RK15} 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 \citep{TR89,R89,RS92}.
This instability was later identified as one that is analogous to
Rayleigh-B\'enard convection, but now for an already convecting mean
state \citep{TBMR94}.
This instability would lead to giant cell convection.

The existence of giant cell convection is under debate, but assuming that it
does {\em 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{equation}
\mbox{\rm Ra}_{\rm t}={gd^4\over\nu_{\rm t}\chi_{\rm t}}\,
\left({\Delta\nabla\over H_{\!P}}\right)_{z_*}.
\label{Raturb}
\end{equation}
This number would then have to be still below the critical value for
convection.
\EEq{Raturb} differs from the usual one defined through \Eq{RaDef} in
\App{Rayleigh} 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 $\nabla_{\rm D}$ term
were to turn out untenable.

We recapitulate that in this alternate explanation, the stratification
is Schwarzschild unstable, i.e.\ $\mbox{\rm Ra}_{\rm t}>0$, corresponding to \Eq{unstable},
but $\mbox{\rm Ra}_{\rm t}<\mbox{\rm Ra}_{\rm t}^{\rm crit}$ is still below a certain critical value, so it would
be stable by a turbulent version of the Rayleigh-B\'enard criterion.
\cite{TBMR94} found $\mbox{\rm Ra}_{\rm t}^{\rm crit}\approx300$ for a vanishing rotation rate.
However, they also estimated that $\mbox{\rm Ra}_{\rm t}>\mbox{\rm Ra}_{\rm t}^{\rm 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 \citep{BRS08,RB12}.

\section{Conclusions}
\label{Conclusions}

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, $\nabla-\nabla_{\rm ad}$, but to $\nabla-\nabla_{\rm ad}+\nabla_{\rm D}$,
where $\nabla_{\rm D}$ is a new nonlocal contribution that was first identified
by \cite{Dea66} 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 $\nabla_{\rm D}$ term, in addition to a $\nabla_{\rm 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 $\nabla$, but the term
$\nabla-\nabla_{\rm ad}$ is now replaced by $\nabla-\nabla_{\rm ad}+\nabla_{\rm D}-\nabla_{\rm kin}$.
This new formulation implies that convection can carry a finite flux while
$\nabla-\nabla_{\rm ad}$ is still negative and therefore the stratification
is Schwarzschild stable, i.e., \Eq{stable} 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 $\nabla-\nabla_{\rm 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 \cite{Greer} 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 \cite{Rem04}, who studied a semianalytic
overshoot model that was driven nonlocally by downdraft plumes,
similar to what was suggested by \cite{Spr97}.
Rempel also finds an extended subadiabatic layer in large parts
of the model.
Furthermore, there are similarities to the nonlocal mixing length model by
\cite{XD01}, 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
$\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 \Eq{FradSplit}.
On the other hand, one could study suitably parameterized large eddy
simulations that either include a nonlocal Deardorff term of the form
given by \Eq{nabD_ansatz}, as discussed in \Sec{Deardorff},
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 \cite{Hanasoge} are, they do require
further scrutiny and call for the resolution of the existing conflicts
with other helioseismic studies such as those of \cite{Greer}.
Alternatively, global helioseismic techniques for detecting giant cell
convection \citep{LR93,CA09,Woo14} can provide another independent way
of detecting deep larger-scale flows \citep{Woo16}.
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
\Eq{FradSplit} 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 time scales.

\acknowledgments
I thank the referee for his/her criticism that has led to many
improvements, and Evan Anders, Jean-Francois Cossette, Ben Greer,
{\AA}ke Nordlund, Mark Rast, Matthias Rempel, Matthias Rheinhardt,
Bob Stein, Peter Sullivan, Regner Trampedach, and J\"orn 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.

