# Dynamics and Scaling of Internally Cooled Convection

Lokahith Agasthya  Caroline Jane Muller

###### Abstract

Our goal is to investigate fundamental properties of the system of internally
cooled convection. The system consists of an upward thermal flux at the lower
boundary and a constant cooling term in the bulk. This simple model represents
idealised dry convection in the atmospheric boundary layer, where the bulk
cooling mimics the radiative cooling to space notably through longwave
radiation. We introduce a lapse-rate dependent dimensionless Rayleigh-number
that determines the behaviour of the system and the critical Rayleigh-number
for transition to convection is obtained numerically. The dynamic behaviour of
the fluid system is described and we also infer the scaling of various
important measured quantities such as the total vertical convective heat flux
and the upward mass flux. We conclude with a short discussion on the relevance
to atmospheric convection and the scope for further investigations of
atmospheric convection using similar simplified approaches.

††journal: Communications in Nonlinear Science and Numerical Simulation

[aff1]organization=Institute of Science and Technology Austria,addressline=Am
Campus 1, city=Klosterneuburg, postcode=3400, country=Austria

##  1 Introduction

The study of thermal, convective systems has a long history going back to the
late 19th century, with the recorded observations of James Thomson ([1]) and
the systematic experiments performed by Henri Bénard [2] setting the stage for
further investigation into the effects of heating a fluid. Lord Rayleigh was
the first to analytically describe the convective instability resulting from
heating a fluid from below [3] using a model system that was named Rayleigh-
Bénard (R-B) convection and one that has become the bedrock of studies on the
behaviour of thermal fluids and phenomena such as pattern formation,
transition to chaos etc.

One of the primary motivations to study R-B convection and other similar model
fluid systems is to understand and characterise the behaviour of the earth’s
atmosphere. The atmospheric circulation of the earth is driven primarily by
the heating of the earth’s surface by the sun, which in turn heats the lowest
level of the atmosphere. In addition to simple thermal effects, atmospheric
convection includes a large number of physical, chemical and biological
processes at different scales. The study of atmospheric convection is usually
performed by solving a large set of coupled non-linear equations which
represent all these processes. These could be General Circulation Models (GCM)
at the global scale or Cloud Resolving Models (CRM) at smaller scales. While
they show realistic behaviour, their complexity, especially the large number
of state variables in the models makes their results hard to interpret and can
even obscure a more fundamental understanding of atmospheric processes.

In the context of atmospheric convection, moist convection (ie., convection in
the presence of moisture) is one of the most important phenomenon in the
tropics. Here, the release of latent heat due to phase changes in water is
very important to the dynamics, while the micro-physical details of the
nature, number and size of hydrometeors (water in liquid or solid form) have
an important feedback on the convective scale dynamics. Recently several
studies have focused on developing simplified models which ignore several
processes at various scales but capture essential features of the real
atmosphere - for some examples, see [4, 5, 6]. These studies considered moist
convection with highly simplified cloud microphysics, moist convection with
varying thermal boundary conditions and moist-convection without considering
the dynamics of liquid water respectively. They found robust evidence of an
up-down asymmetry (where hot, rising updrafts occupy a smaller fraction of the
domain than cold, subsiding downdrafts) and scalings with different parameters
for important measured quantities such as the heat-flux, the mass flux in
clouds, the height of maximum cloud formation, etc. Further exploration of
such simplified convective models with some basic processes of the atmosphere
represented while omitting or greatly simplifying other processes has the
potential to uncover which are the fundamental processes which set the
dynamics of tropical moist convection and which do not greatly affect the
behaviour of moist convection.

In this study, we focus on the simpler case of dry convection which occurs in
the atmosphere when moisture is either absent or present in small enough
quantities that condensation or freezing can be neglected. In other words, we
consider convection of dry air or of moist air (i.e. with water vapor), but
without phase change of water vapor into liquid or ice. For simplicity we will
derive the equations for dry air, but we note that water vapor could easily be
accounted for by replacing temperature with virtual temperature [7]. Dry
convection occurs in the region between the earth’s surface and the cloud-
base, known as the sub-cloud layer. Dry convection is prevalent in the tropics
over dry land and is particularly important in studies of the planetary
boundary layer. Here, the constant cooling of the atmosphere by the emission
of longwave radiation plays a non-negligible role in the dynamics.

To study this convective boundary layer, an idealisation that is most often
made is to study the model fluid system with a fixed temperature at the lower
boundary and (unlike Rayleigh-Bénard convection) a constant rate of diabatic
cooling everywhere in the domain. An analysis of the stability of the system
can be found in Chapter 3.6 of [7], while the results from more recent
simulations are reported in [8].

Berlengiero et. al. in their study [9] (henceforth referred to as B2012)
investigated the convective behaviour of a layer of fluid with a fixed heat-
flux at the lower surface and a uniform bulk-cooling term to understand the
basic features of atmospheric dry convection with a constant radiative cooling
to space. They argued that in some scenarios, particularly over land, fixed-
flux boundary conditions is a better representation of the heating of the
atmosphere by the surface than fixed temperature. The system consists of a
layer of a Boussinesq fluid with a fixed input heat-flux $f_{0}$ at the lower
surface, no vertical heat-flux at the upper boundary, a constant bulk-cooling
$-R$ and a set adiabatic temperature lapse rate $\gamma$. The heat input into
the system and the cooling via the bulk-cooling term must be in global balance
to ensure that a system at thermal equilibrium is achieved.

In B2012, the authors described the dynamics of two different 3D flows, one
with a non-zero, finite temperature stratification and a second flow with no
stratification. The results presented included the vertical temperature
profile for both cases, the volume fraction of the domain occupied by
updraughts and the skewness of the vertical velocity. In particular, it was
found that for both cases, less than half the volume was occupied by upward
velocities, with concentrated intense hot plumes and large regions of
subsiding flow. It was also found that the flow with $\gamma=0$ showed strong
clustering of hot plumes due to the interaction between the rising plumes and
the returning subsiding flow from the upper boundary.

This set-up has also been studied in [10], albeit for a adiabatic lapse-rate
$\gamma$ set to $0$. The focus in their work was on the difference in the
nature of the boundary layer and resulting net heat transfer when the heating
was distributed over a small height near the lower surface and not only the
lower boundary.

In this study, we continue and expand the exploration of the internally cooled
system for a wide-range of parameters, with the goal to derive scalings for
important physical quantities, and identify different behaviors as a function
of key adimensional parameters. We study and characterise how various
important response parameters, particularly the temperature profiles, the up-
down asymmetry and the total convective heat flux as well as the convective
mass flux change with changing the various input parameters of the model.
Further, we also add a discussion on the relevance to the scenario of dry
atmospheric convection.

The article is laid out as follows. Section 2 discusses the basic equations of
our model, a recap of the basic relations derived in B2012 and a description
of the numerical methods. In Section 3 we derive some more analytical results,
define appropriate length, velocity and temperature scales to derive the non-
dimensional parameters characterising the system and end with a discussion on
typical values for various parameters in the atmosphere. Section 4 presents
the main results from our numerical simulations and we conclude with a
discussion on the significance of the results as well as avenues for future
work in Section 5.

##  2 Methodology: Equations and numerical simulations

###  2.1 Equations

As in B2012, we start with the incompressible Boussinesq fluid equations along
with the heat equation with a bulk-cooling term. We thus write the density as
the sum of a reference constant density $\rho_{ref}$ and a perturbation from
this reference density denoted by $\rho^{\prime}$ such that
$\rho=\rho_{ref}+\rho^{\prime}$, with corresponding pressure
$p=p_{ref}(z)+p^{\prime}$ where $p_{ref}(z)$ is in hydrostatic balance with
$\rho_{ref}$ so that $dp_{ref}/dz=-\rho_{ref}g$. The equations for the total
fluid velocity vector $\boldsymbol{u}$ with vertical component $w$ and
temperature fluctuations $T^{\prime}$ about a reference temperature $T_{ref}$
are

$\displaystyle\boldsymbol{\nabla}\cdot\boldsymbol{u}=0,$ (1)
$\displaystyle\partial_{t}\boldsymbol{u}+(\boldsymbol{u}\cdot\boldsymbol{\nabla})\boldsymbol{u}=-(1/\rho_{ref})\boldsymbol{\nabla}p^{\prime}+\nu\nabla^{2}\boldsymbol{u}-\beta
T^{\prime}\boldsymbol{g},$ (2)
$\displaystyle\partial_{t}T^{\prime}+\boldsymbol{u}\cdot\boldsymbol{\nabla}T^{\prime}+\gamma
w=\kappa\nabla^{2}T^{\prime}-R,$ (3)

where $\nu$ is the kinematic viscosity, $\beta$ is the thermal expansion
coefficient, $\boldsymbol{g}=-g\hat{z}$ is the acceleration due to gravity and
$\gamma=g/c_{p}$ is the lapse-rate of the mean stratification (or the dry
adiabatic lapse-rate, see A). $R$ ($>0$) is the bulk-cooling term applied to
the fluid. As discussed by B2012, the term $R$ breaks the up-down symmetry of
the internally cooled system, leading to downward moving fluid ($w<0$)
occupying more than half the domain.

Here we have assumed that the density perturbation $\rho^{\prime}$ is given by
$\rho^{\prime}(T^{\prime})=-\rho_{ref}\beta T^{\prime}$, $\rho_{ref}$ being
the density of the fluid at temperature $T_{ref}$. For the Boussinesq
approximation, it is only important that $T^{\prime}$ remains small enough
that $\rho(T)/\rho_{ref}\ll 1$ or equivalently, $\beta T^{\prime}\ll 1$. For
simplicity, we henceforth drop the primes and investigate equations
$\eqref{eq:incomp}$ \- (3) without the primes.

The velocity and the temperature fields are periodic in the horizontal
directions. At the bottom ($z=0$) and top ($z=L_{z}$) surfaces, they are given
by

$\boldsymbol{u}(z=0)=\boldsymbol{u}(z=L_{z})=0$ (4)

and

$\displaystyle\frac{\partial T}{\partial
z}\biggr{\rvert}_{z=0}=-f_{0},\qquad\frac{\partial T}{\partial
z}\biggr{\rvert}_{z=L_{z}}=-f_{1}.$ (5)

It is noteworthy that this system is invariant under the transformation $T\to
T+\delta T$ where $\delta T$ is some constant temperature since the heat
equation as well as the boundary conditions contain only derivative terms of
the temperature. In moist convection, this is no longer true as the partition
of water into solid, liquid and gaseous phases strongly depends on the precise
thermodynamic temperature.

In addition to the fluid equations, we recall from B2012 that the large-scale
thermal balance for the system, which can easily be derived by taking the
domain average of eqn. 3 in steady state, is given by

$\kappa(f_{0}-f_{1})=RL_{z}.$ (6)

As in B2012, we set $f_{0}$ to be positive and $f_{1}=0$, leading to a balance
between the radiative cooling in the domain and the net heat-flux from the
boundaries. This also ensures that the average temperature of the domain
remains stationary in time.

A steady-state solution for the equations with $\boldsymbol{u}=0$ everywhere
leads to a temperature gradient

$\frac{dT}{dz}=-f_{0}\left(1-\frac{z}{L_{z}}\right)$ (7)

with the corresponding temperature profile given by

$T(z)=T_{bot}-zf_{0}\left(1-\frac{z}{2L_{z}}\right),$ (8)

where $T_{bot}$ is the temperature at $z=0$. The fluid is stably stratified
where $\partial_{z}T>-\gamma$, which is given by $z>L_{z}(1-\gamma/f_{0})$.
Thus we expect that when convection occurs, it occurs only up to a height
$z_{0}$ given by

$z_{0}=L_{z}(1-\gamma/f_{0}).$ (9)

###  2.2 Numerical Methods

In order to understand and investigate the behaviour of the above system, we
solve the system of equations (1)-(5) numerically in a 2D domain box with
fixed aspect ratio $2\pi$ which is periodic in the horizontal direction. The
various input parameters other than the extent of the domain are $\gamma$,
$\kappa$, $\nu$, the product $\beta g$ and $f_{0}$. We $f_{1}$ is set to $0$,
$\rho_{ref}$ to unity while $R=\kappa f_{0}$ in accordance with the condition
of thermal equilibrium (eqn. (6)). We also perform simulations for a
horizontally periodic 3D domain with the same aspect ratio for a few selected
parameters. The equations are solved with a python code employing the Dedalus
spectral solver [11], using Fourier components in the periodic horizontal
directions and Chebyshev polynomial decomposition for the vertical direction.
The code was bench-marked against the results from B2012.

The outputs of the simulations are analysed to calculate various flow
parameters only after the flow has become statistically stationary, wherein
the average kinetic energy and the average temperature of the domain as well
as the average temperature at the top and bottom surfaces of the domain
fluctuate about a constant value. While we report only a few measurements from
3D runs, we have checked that all the below results are qualitatively true for
2D as well as 3D flows, unless otherwise mentioned.

##  3 Model Overview

###  3.1 Exact Relations

In addition to the large-scale thermal balance and the steady-state
temperature profile (eqns. (6) and (8)), we also derive the z-dependent
thermal balance as is standard in the case of Rayleigh-Bénard convection and
other model thermal fluid systems. Rewriting the advective term in flux-form
($\boldsymbol{\nabla}\cdot(\boldsymbol{u}T)$) and integrating the horizontal
statistical average of eqn. (3) from the bottom to a height $z$ gives

$\langle wT\rangle_{A,t}-\kappa\langle\partial_{z}T\rangle_{A,t}=R(L_{z}-z),$
(10)

where $\langle\cdot\rangle_{A,t}$ indicates the statistical average taken at a
fixed height $z$ and we have used the fact that in the steady-state, the time-
derivatives go to $0$ while horizontal derivatives of the flux term average
out to zero due to the periodic boundary conditions. This relation shows that
the total vertical heat transfer by convection decays linearly as a function
of height. Heat is removed uniformly at each height by the bulk radiative
cooling term. In case of $f_{1}$ not being set to $0$, there would be an
addition term $\kappa f_{1}$ on the RHS which represents the outward heat-flux
at the top of the domain.

For our system, we also calculate the global thermal dissipation
$\epsilon_{T}\equiv\kappa\langle|\boldsymbol{\nabla}T|^{2}\rangle$ given by
(see B.1)

$\displaystyle\kappa\langle|\boldsymbol{\nabla}T|^{2}\rangle$
$\displaystyle=\frac{\kappa}{L_{z}}T_{bot}f_{0}-\gamma\langle
wT\rangle-R\langle T\rangle$ $\displaystyle=R\langle
T_{bot}-T\rangle-\gamma\langle wT\rangle.$ (11)

where $\langle\cdot\rangle$ indicates the statistical average taken over the
whole domain and $T_{bot}$ is the measured average temperature at $z=0$. We
can see that that the RHS value still remains invariant under the
transformation $T\to T+\delta T$. The strength of the thermal gradients in the
fluid increases with the average departure of the domain temperature from
$T_{bot}$ and decreases when convective transfer of heat is efficient, that is
$\langle wT\rangle$ is large.

The viscous dissipation $\epsilon\equiv(\nu/2)\sum_{i,j}(\partial
u_{i}/\partial x_{j}+\partial u_{j}/\partial x_{i})^{2}$ is given by (see B.2)

$\frac{\nu}{2}\sum_{i,j}\left(\frac{\partial u_{i}}{\partial
x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)^{2}=\beta g\langle
wT\rangle.$ (12)

We can see that, similar to the Nusselt number in R-B convection, the response
quantity $\langle wT\rangle$ is of fundamental importance and it determines
the convective heat transfer properties of the fluid system as well as the
strength of the thermal and kinetic gradients. In the context of a convecting
atmosphere, this value is the vertical buoyancy-flux, crucial in determining
the thermodynamics of the atmosphere in shallow dry convection. In deep, moist
convection too, the thermal flux in the sub-cloud layer is of crucial
importance - for example it is approximately equal to the pressure-work done
by the fluid and this is closely related to the production of irreversible
entropy in the sub-cloud layer [12].

Analogous to the usual definition for R-B convection we define the Nusselt
number Nus for the current system as

$\text{Nus}=\frac{\langle wT\rangle}{\kappa T_{0}/L_{z}}$ (13)

where $T_{0}$ is a temperature scale of the system to be defined later. The
Nusselt number here compares the vertical heat-flux with the typical
conductive heat flux. When the flow is conductive, we have no fluid motion and
hence $\text{Nus}=0$. In the presence of convection, Nus has a finite value
which increases with the increase in the strength of convection.

###  3.2 Non-dimensional parameters

To quantify the system and non-dimensionalise the equations, we need
appropriately defined temperature, velocity and length scales. Since the
dynamics of the system is set by the bulk-cooling $R$, which has SI unit
dimensions of $\text{Kelvin s}^{-1}$, we can write the temperature scale
$T_{0}$ as

$T_{0}=Rt_{0}$ (14)

where $t_{0}$ is the appropriate time-scale, with length-scale $L_{z}$ giving
the velocity scale $U_{0}=L_{z}/t_{0}$. As is standard for thermal flows, we
write $U_{0}$ as a convective velocity given by

$U_{0}=\sqrt{\beta gT_{0}L_{z}}.$ (15)

Combining the above equations gives finally

$\displaystyle T_{0}=(\beta g)^{-1/3}R^{2/3}L_{z}^{1/3},$ (16) $\displaystyle
U_{0}=(\beta g)^{1/3}R^{1/3}L_{z}^{2/3},$ (17) $\displaystyle t_{0}=(\beta
g)^{-1/3}R^{-1/3}L_{z}^{1/3}.$ (18)

Using these scales, the equations are non-dimensionalised as

$\displaystyle\widehat{\boldsymbol{\nabla}}\cdot\widehat{\boldsymbol{u}}=0,$
(19)
$\displaystyle\partial_{\widehat{t}}\widehat{\boldsymbol{u}}+(\widehat{\boldsymbol{u}}\cdot\widehat{\boldsymbol{\nabla}})\widehat{\boldsymbol{u}}=-\widehat{\boldsymbol{\nabla}}\widehat{p}+\sqrt{\frac{\Pr}{\text{Ra}}}\widehat{\nabla}^{2}\widehat{\boldsymbol{u}}+\widehat{T}\hat{z},$
(20)
$\displaystyle\partial_{\widehat{t}}\widehat{T}+\widehat{\boldsymbol{u}}\cdot\widehat{\boldsymbol{\nabla}}\widehat{T}+(\gamma/f_{0})\sqrt{\Pr\text{Ra}}\widehat{w}=\frac{1}{\sqrt{\Pr\text{Ra}}}\widehat{\nabla}^{2}\widehat{T}-1,$
(21)

where the hat ( $\widehat{}$ ) indicates the non-dimensionalised variable or
operator obtained by dividing by the appropriate dimensional scale. The non-
dimensional parameters are the Rayleigh-Number Ra, the Prandtl number $\Pr$
and the ratio $\gamma/f_{0}$ between the lapse-rate and heat-flux. The
Rayleigh number is given by

$\text{Ra}=(\beta g)^{2/3}\frac{R^{2/3}L_{z}^{10/3}}{\nu\kappa}.$ (22)

This is equivalent to the form derived in B2012, albeit without appealing to a
flux formulation. This equivalence is not surprising as the approach of B2012
(setting the time-scale of the fluid motion equal to the time-scale of
radiative cooling) is identical to our approach in setting the temperature
scale according to the rate of radiative cooling, which gives a non-
dimensional rate of radiative cooling ($RL_{z}/(U_{0}T_{0})$) of unity. Ra can
also be re-written in terms of $f_{0}$ using $R^{2/3}L_{z}^{10/3}=(\kappa
f_{0})^{2/3}L_{z}^{8/3}$ from eqn. (6), which gives $\text{Ra}\propto
f_{0}^{2/3}/\kappa^{1/3}$.

The Prandtl number is the ratio $\nu/\kappa$ between the viscosity and the
thermal conductivity. The ratio $\gamma/f_{0}$ sets the height $z_{0}$ above
which the fluid is stably stratified (see Eqn. (9)). It can indeed be argued
that the height $z_{0}$ rather than $L_{z}$ is the appropriate length scale to
define the convective time and velocity-scales, which leads to a lapse-rate
dependent Rayleigh number $\text{Ra}_{\gamma}$ given by

$\text{Ra}_{\gamma}=(\beta
g)^{2/3}\frac{R^{2/3}z_{0}^{10/3}}{\nu\kappa}=\text{Ra}\left(1-\frac{\gamma}{f_{0}}\right)^{10/3}.$
(23)

In this study, we solve the system of equations in arbitrarily chosen
simulations units - however, all the results presented involve only non-
dimensional quantities and it is the scaling with the non-dimensional
parameters that we are interested in to characterise the system.

While several of our simulations have large differences of temperature between
the lower and upper surface, thus invalidating the assumptions underlying the
Boussinesq approximation, we note that this temperature difference is dictated
by a combination of the temperature gradient at the lower boundary $f_{0}$ and
the lapse-rate $\gamma$. The principle of dynamic similarity for flows with
the same $\text{Ra}_{\gamma}$ can be used to construct an equivalent flow with
$f_{0}$ chosen to be an appropriately small value and also decreasing $\kappa$
and $\nu$ to retain the same $\text{Ra}_{\gamma}$ ($\propto
f_{0}^{2/3}/\kappa^{1/3}$).

###  3.3 Some Orders of Magnitude

In the context of a typical scenario of dry convection in the tropics,
convection occurs in a layer of atmosphere from the surface to cloud base
level, which is $500\text{ m}$ to $3\text{ km}$ high, with the thermal heat
flux (sensible heat flux) from the surface in the order of $10-100\text{ W
m}^{-2}$, radiative cooling to space of $1-2\text{ K day}^{-1}$ and dry air
with a typical $\beta\sim 1/300\text{ K }^{-1}$.

Sensible heat flux in $\text{ W m}^{-2}$ can be converted to an equivalent
thermal gradient by dividing by $(\rho\kappa c_{p})$ where $\rho$ is the
density and $c_{p}$ is the heat capacity of air. When we consider an eddy-
diffusivity like value of thermal conductivity $\kappa\sim 10^{-2}\text{
m}^{2}\text{s}^{-1}$, with $\rho\sim 1\text{ kg m}^{-3}$ and $c_{p}\sim
10^{3}$ $\text{ J kg}^{-1}\text{ K}^{-1}$ we obtain $f_{0}\sim\mathcal{O}(1)$.
Considering the adiabatic lapse-rate to be the typical dry value of
$0.01\text{ K m}^{-1}$ and $\beta g\sim 1/30\text{ m s}^{-2}\text{K}^{-1}$
gives $Ra_{\gamma}$ to be at the least of order $10^{9}$. For the true
molecular value of $\nu$ and $\kappa$, the value of $\text{Ra}_{\gamma}$ would
be closer to $10^{15}$. It is important to note here that even for the lowest
estimates of $f_{0}$, $\gamma$ is less than $1$% of $f_{0}$ iiiIf we consider
real-world observations, it is well-known that ambient temperature
measurements are recorded $2$ metres above the land surface, which is typical
a few Kelvin warmer than the air at the lowest level of the atmosphere on a
sunny day. Thus, the temperature gradient at the boundary here is of the order
of $1$ Kelvin per metre, compared to the dry adiabatic lapse rate of
$0.01\text{ K m}^{-1}$. .

Given that the value of $10\text{ W m}^{-2}$ for the sensible heat-flux is a
global annual average, there do exist scenarios where the sensible heat flux
can be far smaller than this value. However, in these situations it is
unrealistic to model the atmosphere to be locally in equilibrium with
radiative cooling - instead the energy balance would be determined by large-
scale, horizontal heat-fluxes or the cooling term would dominate the heat-
equation. Studies of adjustment in such an out-of-equilibrium state is beyond
the scope of the current study.

In using the current system as a model for atmospheric dry-convection, it is
also important to assess the validity of the Boussinesq-approximation for dry
convection. For a layer of atmosphere around $1$ km thick, temperature
differences $\Delta T$ between the bottom and the top of this layer are of the
order of $\sim 10\text{ K }$, as the atmosphere remains close to a dry
adiabatic profile. This translates to $\beta\Delta T\approx 3\times 10^{-2}$.
Thus, we see that the Boussinesq approximation remains an excellent
approximation in the simulation of dry convection. In the presence of moisture
and the ensuing deep, moist convection that penetrates up to the top of the
troposphere ($\sim 15$ km), the Boussinesq approximation breaks down and the
vertical variation in density is significant.

From the above discussion it is clear that the current model under
investigation is valid for shallow, sub-cloud dry convection where it can be
assumed that the convection and sensible heat-flux is locally in equilibrium
with the cooling of the atmosphere through outgoing long-wave radiation to
space. In Table 1, we show the range of dimensionless parameters simulated in
the current study.

Parameter | Range
---|---
$\text{Ra}_{\gamma}$ | ($50$, $10^{6}$)
$\gamma/f_{0}$ | ($0$, $0.75$)
$\Pr$ | $1$

Table 1: The range of dimensionless parameters studied

$\gamma/f_{0}$ |  Critical $\text{Ra}_{\gamma}$
---|---
$0$ | ($130.162$, $130.9497$)
$0.2$ | ($86.966$, $88.017$)
$0.4$ | ($67.93$, $69.3375$)
$0.75$ | ($77.748$, $78.412$)

Table 2: The range of the critical $Ra_{\gamma}$ for transition to convection
found numerically for different values of $\gamma/f_{0}$ in 2D. Flows with
lower bound of $Ra_{\gamma}$ were found to be stationary while flows with
$Ra_{\gamma}$ equal to the upper bound were found to be convective.

##  4 Results

Figure 1: Instantaneous snapshots of the temperature field $T$ in the
stationary regime for flows with different $\text{Ra}_{\gamma}$ as indicated
in the figure and $\gamma/f_{0}=0.2$. The temperature fields are divided by
the respective temperature scale $T_{0}$ defined in Eqn. (16) and scaled such
their domain average is $0$. Panel (a) shows a stable flow while panels (b),
(c) and (d) show convective flows. The arrows represent velocity vectors with
the arrow length in plot units equal to velocity magnitude divided by $U_{0}$,
$U_{0}$, $2U_{0}$ and $5U_{0}$ respectively where $U_{0}$ is the velocity
scale defined in Eqn. (17).

###  4.1 Transition to convection

Figure 2: (a) Average vertical profiles of $\partial_{z}T$ divided by $f_{0}$.
(b) Average vertical profiles of $wT$ divided by $\kappa T_{0}/L_{z}$. Both
panels share the same legend, with the flows having $\text{Ra}_{\gamma}\sim
10^{4}$ and the averages are statistical, height-wise averages for each flow.

With $\gamma=0$, it is found that for very small values of $f_{0}$ (or
equivalently, small values of $R$), the flow remains motionless
($\boldsymbol{u}=\boldsymbol{0}$) and assumes a temperature profile given by
Eqn. (8). When $f_{0}$ is increased beyond a critical value, the flow becomes
convective with convective rolls that extend up to height $L_{z}$ since
$z_{0}=L_{z}$ when $\gamma=0$. For finite values of $\gamma$, the flow remains
stable for larger values of $f_{0}$ than for the $\gamma=0$ case. On
increasing $f_{0}$ beyond the critical value, the flow becomes convective
where the convective rolls extend up to height $z_{0}$. For the case of
$\gamma=f_{0}$, the flow is always stable as in this case no heat-transfer
occurs from the lower surface to the fluid. We find numerically that the
critical $\text{Ra}_{\gamma}$ for the triggering of convection depends weakly
on the ratio $\gamma/f_{0}$. The numerically estimated values using 2D
simulations are given in Table 2. The goal of the current study is to
characterise the behaviour of the fully convective regime rather than a
detailed study of the fluid instability. We note however that the estimation
of the critical $\text{Ra}_{\gamma}$ for different values of $\gamma$ should
be analytically tractable through the analysis into normal modes of a small
perturbation on the steady-state solution as is usually done for classical
Rayleigh-Bénard convection (for example see [13]) and other instabilities.

Figure 1 shows instantaneous snapshots of the linearly scaled, non
dimensionalised temperature field for 2D flows with $\gamma/f_{0}=0.2$. When
the flow is stable (see panel (a)), the fluid is held motionless by viscosity
and thermal dissipation and the temperature field is perfectly homogeneous in
the horizontal direction. When $f_{0}$ is increased slightly, as in panel (b),
the flow shows a pattern of large, space-filling convective rolls. Larger
values of $f_{0}$ shown in the lower panels show intense and clearly defined
rising convective plumes, where the hot, rising plumes have a large positive
height-wise temperature anomaly, while the gently subsiding regions outside
these plumes have a smaller negative height-wise temperature anomaly. In panel
(d), the hot plume is the narrowest with the highest velocities (seen by
length of the arrows) strongly concentrated in the rising plumes. The reader
should note that the arrows representing the velocity fields are scaled by a
different value in each panel only to ensure clarity in the figure - arrow
lengths must not be used to compare the velocity between different panels.

###  4.2 Flow structures and Profiles

In the convective regime, intense hot plumes arise from the lower surface
along with broad regions of cold subsidence. There exists a small diffusive
boundary layer near the lower surface underneath the intense hot plumes. Here,
the vertical gradient of the temperature field ($\partial_{z}T$) quickly goes
from $-f_{0}$ to $-\gamma$ as seen in Figure 2(a). Above this layer lies a
convective bulk region where the effective mixing of the fluid ensures
$\partial_{z}T\approx-\gamma$ up to a height of approximately $z_{0}$. For
small $\gamma$, this corresponds to the top of the domain and thus there is no
conducting region close to the upper boundary. Above $z_{0}$, $\partial_{z}T$
goes to $0$ at $z=L_{z}$ in a linear fashion, which corresponds to the stable
layer above the convective layer.

Figure 3: Nusselt number of the flow as a function of $\text{Ra}_{\gamma}$ for
$4$ different values of $\gamma/f_{0}$ in 2D and $2$ different values in 3D.

Panel (b) of Figure 2 shows the time-averaged convective flux divided by the
product of a diffusive velocity scale $\kappa/L_{z}$ and the temperature scale
$T_{0}$ \- that is the height-wise Nusselt number. The convective flux attains
a maxima at $z$ roughly corresponding to the end of the lower thermal boundary
layer. Above this layer and up to $z=z_{0}$, ie., in the convective region
where $\partial_{z}T$ remains nearly constant, the flux decreases linearly
consistent with eqn. (10).

For larger values of $\gamma$, there exists a small layer of fluid close to
$z\gtrsim z_{0}$ where the flux is negative - this can be explained by the
formation of cold patches of fluid above the rising thermal plumes. The fast
rising parcels of fluid are cooled rapidly due to the lapse-rate term, leading
to a region where the fluid is colder than the horizontal average while the
fluid is still rising ($w>0$), leading to a negative flux. Above this, the
fluid is no longer convecting, leading to a flux close to $0$. Thus, in the
presence of a large lapse-rate, the domain-averaged convective heat-flux (Nus)
is partially decreased by the strong cooling of fast-rising parcels of fluid.

###  4.3 The Nusselt Number

The convective heat-flux characterises the behaviour of the flow to a large
degree, as it determines the strength of the velocity and temperature
gradients and also plays an import role in understanding the behaviour of the
fluid in the viscous boundary layers. The non dimensionalised heat-flux, the
Nusselt number (defined in Eqn. (13)) is thus particularly important in
understanding the convective behaviour of the system.

The Nusselt number is identically $0$ for conductive flows as the fluid is not
in motion. In the convective regime, Nus increases with increasing
$\text{Ra}_{\gamma}$ as we see in figure 3, which shows the scaling of the
Nusselt number for changing $\text{Ra}_{\gamma}$ and for various values of
$\gamma/f_{0}$. For flows with identical $\text{Ra}_{\gamma}$, Nus is the
largest for the cases with $\gamma=0$. This is due to the fact that for flows
without an imposed lapse rate, the convection is penetrative ($z_{0}=L_{z}$)
and the hot plumes reach the top of the domain. In the case of flows with a
finite lapse rate, the region of the flow above $z_{0}$ contributes negatively
to the total convective heat flux of the system, as already discussed in the
previous section and seen in Figure 2(b). We note that this contribution is
quite small, given that the values are nearly identical for $\gamma/f_{0}$
ranging from $0$ to $0.75$. This also indicates that $\text{Ra}_{\gamma}$ is a
well-chosen parameter to characterise the large-scale behaviour of this model
system.

The magnitude and scaling of the Nusselt number as a function of
$\text{Ra}_{\gamma}$ for 3D flows is identical to the 2D flows. This indicates
that in making estimates of dry convective heat-fluxes with given boundary
conditions, it is sufficient to simulate a 2D domain.

###  4.4 Mass-flux, Updraft Velocities and Up-down Asymmetry

Along with the convective heat-flux, the convective mass-flux also remains an
important parameter to be estimated in thermal flows, particularly in the dry
atmosphere. The mass flux at a given height is calculated as the product of
the volume fraction $V_{f}$ of the domain which is rising ($w>0$) and the
average vertical velocity $U_{u}$ within these rising regions. The total
convective mass flux is the vertical average of the mass-flux at each height
written as $\langle V_{f}U_{u}\rangle$.

Figure 4: (a) Domain-averaged convective mass flux plotted against
$\text{Ra}_{\gamma}$ for $4$ different values of $\gamma/f_{0}$ in 2D and $2$
values of $\gamma/f_{0}$ in 3D. The mass flux is non-dimensionalised by
dividing it by the diffusive velocity scale $\kappa/L_{z}$. (b) Average
velocity $U_{u}$ in regions of the flow where $w>0$ divided by the diffusive
velocity scale as a function of $\text{Ra}_{\gamma}$ for the same parameters
as panel (a). Inset (i) shows the volume fraction $V_{f}$ of the domain
occupied by updrafts ($w>0$).

The variation of the convective mass-flux with $\text{Ra}_{\gamma}$ is shown
in panel (a) of figure 4. We see that the mass flux close to the transition to
convection shows large variation in magnitude depending on the value of
$\gamma/f_{0}$, with the mass flux being greater for the flows with greater
$\gamma/f_{0}$. This is similar to the behaviour of the Nusselt number. For
large enough $Ra_{\gamma}$, all values of mass-flux converge to a single line
on the log-log plot, which corresponds to an exponential dependence on
$\text{Ra}_{\gamma}$. The updraft velocities $U_{u}$ and the volume fraction
$V_{f}$ are shown in panel (b) and inset (i) respectively. For larger values
of $\text{Ra}_{\gamma}$ the values of $U_{u}$ also converge to a single
exponent independent of $\gamma/f_{0}$ while $V_{f}$ converges to a single
value $\sim 0.48$.

The variation of $V_{f}$ with increasing $\text{Ra}_{\gamma}$ with can be
explained from the behavior seen in Figure 1 \- when the flow is convective
and laminar (as is the case for small $\text{Ra}_{\gamma}$ close to the
critical value for transition) the only upward moving regions occur in the
vicinity of the rising plume, and the entire domain is filled by a single,
large convective roll. The rest of the domain is uniformly, weakly subsiding.
This behaviour can be observed for example on closer examination of panel (b)
in Figure 1, where a part of domain extending from approximately $x=3$ to
$x=5.5$ is a clearly delineated upward rising region and the rest of the
domain is subsiding with a much smaller velocity to compensate for the rising
mass in the rising plume. For larger $\text{Ra}_{\gamma}$ flows which are more
turbulent and have greater kinetic energy, the flow contains multiple plumes
and convective rolls with rapid point-wise temporal fluctuations in velocity.
For example, in panel (d) of the same figure, while the region between
$x\approx 1$ to $x\approx 2.5$ is subsiding on average, some grid points
individually have $w>0$ within this region. For the largest values of
$\text{Ra}_{\gamma}$, the flows are highly turbulent and energetic, leading to
a $V_{f}$ value close to half.

The asymmetry between up and down that results from the $-R$ term is important
only in the convective region of the flow. The stable upper layer in flows
with non-zero $\gamma$ are locally in conductive equilibrium and fluid motion
here is only due to the transport of mass and momentum from the convective
region. These regions do not show an up-down asymmetry. Thus, regions with a
thicker stable layer (larger $\gamma/f_{0}$) are overall more symmetric and
the domain average $V_{f}$ is closer to half, independent of
$\text{Ra}_{\gamma}$.

The mass flux and the updraft velocity scales nearly identically in 2D as well
as 3D, though in the 3D case $V_{f}$ converges to a value close to half for
much larger $\text{Ra}_{\gamma}$ (if we assume that the scaling seen in Figure
3(i) holds for larger $\text{Ra}_{\gamma}$). Simulations for 3D flows with
$\text{Ra}_{\gamma}>10^{5}$ were not performed due to the requirement of large
computing resources and are beyond the scope of this study.

While the volume fraction $V_{f}$ gives an insight into the degree of
asymmetry in the flow, it does not capture the relative magnitude of the
velocity of rising and subsiding regions, thus does not fully quantify the
asymmetry. Another measure often used (including in B2012) is the skewness
measure of the vertical velocity field. The skewness $S$ is defined as

$S=\frac{\langle w^{3}\rangle}{(\langle w^{2}\rangle)^{3/2}},$ (24)

where we have used the fact that $\langle w\rangle=0$. A positive skew implies
a left-peaked distribution, where the tail on the positive end is thicker
while the median $w$ is negative. Figure 5 shows $S$ as a function of
$\text{Ra}_{\gamma}$ for various values of $\gamma/f_{0}$. $S$ follows a trend
which is close to the inverse of $V_{f}$ (see figure 4(i)), where the value of
$S$ peaks when the value of $V_{f}$ is at a minima. Similar to $V_{f}$, the
skewness also converges to a value which is between $0$ and $0.5$ for large
$\text{Ra}_{\gamma}$, consistent with the aforementioned small-degree of
asymmetry.

Figure 5: Skewness $S$ of the vertical velocity for flows with varying
$\gamma/f_{0}$ as a function of $\text{Ra}_{\gamma}$ for 2D and 3D
simulations.

It is the $-R$ term that breaks the up-down asymmetry - however we see that
for greater $R$, the flow becomes less asymmetric. Whereas the flows with the
largest degree of asymmetry have $R$ just large enough to overcome the
stabilisation by viscous forces. We note here that the key factor determining
the asymmetry is the relative importance of the bulk-cooling term. Considering
the non-dimensionalised heat-equation, eqn. (21), the LHS value is large when
the typical magnitude of fluctuations in $\widehat{T}$ and
$\widehat{\boldsymbol{u}}$ are large. This is the case when Ra is large. For
cases of smaller and intermediate values of Ra, the magnitude of fluctuations
in $\widehat{T}$ and $\widehat{\boldsymbol{u}}$ are small, so the cooling term
is relatively more important in the dynamics.

##  5 Discussion and Conclusions

We have presented results from 2D as well as 3D simulations of an uniformly
cooled thermal fluid system in the presence of gravity and a lapse rate with a
uniform heat flux at the lower surface and an adiabatic (no heat flux)
boundary condition for the top surface. The fluid is assumed to be
incompressible and we work in the regime of the Boussinesq approximation,
which is a fair approximation for the sub-cloud layer of the earth’s
atmosphere, where the atmosphere is either dry or moist-unsaturated. The
constant, bulk cooling term is set such that it balances the incoming heat
flux at the lower surface, leading to a situation where there is balance
between a constant cooling of the domain and a heating from below, which is
analogous to the atmosphere which is heated by the earth surface and
constantly loses heat to space through longwave radiation. The current work is
a direct extension of previous work by Berlengiero et. al. reported in B2012.

The equations are non-dimensionalised using a temperature and velocity scale
based on the magnitude of the diabatic cooling term. For the length-scale, we
choose the height up to which the domain is unstable to dry convection
$z_{0}$, where $z_{0}$ depends on the ratio between the lapse-rate $\gamma$
and the uniform flux at the lower surface $f_{0}$. This leads to the lapse-
rate dependent Rayleigh number denoted $\text{Ra}_{\gamma}$ which
characterises the system well.

We show that when $\text{Ra}_{\gamma}$ is small, the dissipation of heat and
momentum is sufficient to have a motionless steady-state solution for the
fluid with all the heat transfer occurring through thermal conductivity alone.
Beyond a critical value, the fluid no longer remains motionless and is
unstable to small perturbations which induce a convective motion in the
system. We have further shown that in this convective regime, the scaling of
important measured quantities such as the Nusselt number, the non-
dimensionalised mass-flux and average velocity in updrafts with
$\text{Ra}_{\gamma}$ converge to a single exponent, with some dependence of
the magnitude on the dimensionless ratio $\gamma/f_{0}$. The bulk cooling term
also introduces an up-down asymmetry in the fluid - this asymmetry is
particularly marked when the flow is only weakly convective, ie.
$\text{Ra}_{\gamma}$ is only slightly larger than the critical value. For
strongly convective flows, the degree of the asymmetry as measured by the
volume of the domain occupied by updrafts (regions where $w>0$) and the
skewness of the vertical velocity decreases, where nearly half the domain is
made up of updrafts and the skewness goes close to $0$. While the 3D flows are
identical to the 2D flows in the scaling of heat and mass transfers with
$\text{Ra}_{\gamma}$, they retain a large degree of up-down asymmetry even at
$\text{Ra}_{\gamma}\sim 10^{5}$.

In this study, we have added to the existing body of literature on various
models of thermal convective systems which help us understand atmospheric
convection under different conditions. While the system has been introduced
elsewhere, this is to the best of our knowledge the first study that
systematically varies the parameters for a fixed-flux, internally cooled
convective system, a system with particular relevance to the sub-cloud
atmospheric boundary layer. Our main aim here has been to demonstrate the
simple scaling of the system with respect to various input parameters and
provoke future studies which explore higher Rayleigh numbers as well as mixed
boundary conditions and non-uniform radiative cooling.

Idealised models of convection, including those with moisture and water phase
changes serve to provide a tool to study the complex atmospheric system in a
simplified setting where it is far easier to delineate the dynamic effects and
feedbacks due to individual processes. In the future, we plan to study
idealised models which include moist dynamics as well as global-scale
idealised models which include the effects of rotation.

## Acknowledgments

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Sklodowska-Curie grant
agreement No. 101034413. CM gratefully acknowledges funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation program (Project CLUSTER, Grant Agreement No. 805041).

## Conflict of Interests and Data Availability

The authors declare no conflicts of interest. Datasets and the code used to
generate the data are available from the authors upon reasonable request.

##  Appendix A The heat-equation

We start with the first law of thermodynamics in enthalpy form given by

$c_{p}\frac{DT}{Dt}-\frac{1}{\rho}\frac{Dp}{Dt}=\dot{Q}$ (25)

where $m$ is the mass of the fluid parcel, $c_{p}$ is the specific heat-
capacity at constant pressure and $\dot{Q}$ is the rate of diabatic heating or
cooling per unit mass. We make the assumption ($D_{t}p\approx
w\partial_{z}p_{ref}\approx-\rho_{ref}gw$), which is the condition of
hydrostatic balance. This also assumes that a given parcel always has the same
pressure as its environment, which is a common assumption made in parcel
theory. This gives,

$\partial_{t}T+\boldsymbol{u}\cdot\boldsymbol{\nabla}T+(g/c_{p})w=\dot{Q}/c_{p},$
(26)

where the RHS includes thermal dissipation as well as diabatic cooling, which
is given by $-R$ in this study. $g/c_{p}\equiv\gamma$ is the well known dry
adiabatic lapse-rate.

##  Appendix B Thermal and Viscous dissipation

###  B.1 Thermal Dissipation

The thermal dissipation is defined as

$\epsilon_{T}\equiv\kappa\bigl{\langle}(\partial_{i}T(\boldsymbol{x},t))^{2}\bigr{\rangle}_{V}.$
(27)

Following [siggia1994high]’s analysis for Rayleigh-Bénard convection, we
multiply equation (3) by temperature $T$ and average the product over the
entire domain and time to give

$\frac{1}{2}\frac{d\langle
T^{2}\rangle}{dt}+\frac{1}{2}\bigl{\langle}\boldsymbol{u}\cdot\boldsymbol{\nabla}(T^{2})\bigr{\rangle}+\gamma\bigl{\langle}wT\bigr{\rangle}+\bigl{\langle}RT\bigr{\rangle}\\\
=\kappa\bigl{\langle}T\nabla^{2}T\bigr{\rangle}=\kappa\bigl{\langle}\boldsymbol{\nabla}\cdot(T\boldsymbol{\nabla}T)\bigr{\rangle}-\kappa\bigl{\langle}|\boldsymbol{\nabla}T|^{2}\bigr{\rangle},$
(28)

Since we working in the stationary regime,
($\partial_{t}\langle\cdot\rangle=0$). Further, using the incompressibility
condition (eq. (1)) and the fluid rigid boundary condition ((4)), we get

$\bigl{\langle}\boldsymbol{u}\cdot\boldsymbol{\nabla}(T^{2})\bigr{\rangle}_{V}=\bigl{\langle}\boldsymbol{\nabla}\cdot(\boldsymbol{u}T^{2})\bigr{\rangle}_{V}=0.$
(29)

Then, equation (28) becomes

$\kappa\bigl{\langle}|\boldsymbol{\nabla}T|^{2}\bigr{\rangle}=\kappa\bigl{\langle}\boldsymbol{\nabla}\cdot(T\boldsymbol{\nabla}T)\bigr{\rangle}-\gamma\bigl{\langle}wT\bigr{\rangle}-R\bigl{\langle}T\bigr{\rangle},$
(30)

or

$\epsilon_{T}=\kappa\bigl{\langle}\boldsymbol{\nabla}\cdot(T\boldsymbol{\nabla}T)\bigr{\rangle}-\gamma\bigl{\langle}wT\bigr{\rangle}-R\bigl{\langle}T\bigr{\rangle}.$
(31)

Using the Gauss theorem, the first term of $\epsilon_{T}$ can be written in
terms of a surface integral as

$\kappa\bigl{\langle}\boldsymbol{\nabla}\cdot(T\boldsymbol{\nabla}T)\bigr{\rangle}=\frac{\kappa}{L_{z}}\biggl{[}\Bigl{\langle}T\partial_{z}T\Bigr{\rangle}_{z=L_{z}}-\Bigl{\langle}T\partial_{z}T\Bigr{\rangle}_{z=0}\biggr{]}.$
(32)

Setting $\partial_{z}T|_{z=0}=-f_{0}$ and $\partial_{z}T|_{z=L_{z}}=-f_{1}=0$
gives the expression in eqn. (3.1)

$\kappa\langle|\boldsymbol{\nabla}T|^{2}\rangle=\frac{\kappa}{L_{z}}(T_{0}f_{0})-\gamma\langle
wT\rangle-R\langle T\rangle.$ (33)

###  B.2 Viscous Dissipation

Taking the dot product of eq. (2) with $\boldsymbol{u}$ and taking the
statistical average over the whole volume as above, we get

$\frac{1}{2}\frac{d}{dt}\bigl{\langle}(\boldsymbol{u}\cdot\boldsymbol{u})\bigr{\rangle}+\frac{1}{2}\bigl{\langle}\boldsymbol{u}\cdot\boldsymbol{\nabla}(\boldsymbol{u}\cdot\boldsymbol{u})\bigr{\rangle}\\\
=-\bigl{\langle}\boldsymbol{u}\cdot\boldsymbol{\nabla}p\bigr{\rangle}+\nu\bigl{\langle}\boldsymbol{u}\cdot\nabla^{2}\boldsymbol{u}\bigr{\rangle}+\beta
g\bigl{\langle}wT\bigr{\rangle}.$ (34)

In the stationary state, the terms of the form $d_{t}\langle\cdot\rangle$
vanish. Using the incompressibility condition and the rigid boundary
condition, we have

$\bigl{\langle}\boldsymbol{u}\cdot\boldsymbol{\nabla}(\boldsymbol{u}\cdot\boldsymbol{u})\bigr{\rangle}_{V}=\Bigl{\langle}\boldsymbol{\nabla}\cdot\bigl{[}\boldsymbol{u}(\boldsymbol{u}\cdot\boldsymbol{u})\bigr{]}\Bigr{\rangle}_{V}=0$
(35)
$\langle\boldsymbol{u}\cdot\boldsymbol{\nabla}p\rangle_{V}=\bigl{\langle}\boldsymbol{\nabla}\cdot(\boldsymbol{u}p)\bigr{\rangle}_{V}=0$
(36)
$\begin{split}\langle\boldsymbol{u}\cdot\nabla^{2}\boldsymbol{u}\rangle_{V}&=\frac{1}{2}\bigl{\langle}\boldsymbol{\nabla}\cdot\boldsymbol{\nabla}(\boldsymbol{u}\cdot\boldsymbol{u})\bigr{\rangle}_{V}-\sum_{i,j}\biggl{\langle}\Bigl{(}\frac{\partial
u_{j}}{\partial x_{i}}\Bigr{)}^{2}\biggr{\rangle}_{V}\\\
&=-\sum_{i,j}\biggl{\langle}\Bigl{(}\frac{\partial u_{j}}{\partial
x_{i}}\Bigr{)}^{2}\biggr{\rangle}_{V}\\\
&=-\frac{1}{2}\sum_{i,j}\biggl{\langle}\Bigl{(}\frac{\partial u_{i}}{\partial
x_{j}}+\frac{\partial u_{j}}{\partial
x_{i}}\Bigr{)}^{2}\biggr{\rangle}_{V}.\end{split}$ (37)

So, equation (34) becomes

$\frac{\nu}{2}\sum_{i,j}\biggl{\langle}\Bigl{(}\frac{\partial u_{i}}{\partial
x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\Bigr{)}^{2}\biggr{\rangle}=\beta
g\langle wT\rangle,$ (38)

or

$\epsilon\equiv\frac{\nu}{2}\sum_{i,j}\biggl{\langle}\Bigl{(}\frac{\partial
u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial
x_{i}}\Bigr{)}^{2}\biggr{\rangle}=\beta g\langle wT\rangle.$ (39)

## References

  * Thomson et al. [1882] J. Thomson, et al.,  On a changing tesselated structure in certain liquids,  Proc. Phil. Soc. Glasgow 13 (1882) 464–468.
  * Bénard [1901] H. Bénard,  Les tourbillons cellulaires dans une nappe liquide.-méthodes optiques d’observation et d’enregistrement,  Journal de Physique Théorique et Appliquée 10 (1901) 254–266.
  * Rayleigh [1916] L. Rayleigh,  Lix. on convection currents in a horizontal layer of fluid, when the higher temperature is on the under side,  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 32 (1916) 529–546.
  * Hernandez-Duenas et al. [2013] G. Hernandez-Duenas, A. J. Majda, L. M. Smith, S. N. Stechmann,  Minimal models for precipitating turbulent convection,  Journal of Fluid Mechanics 717 (2013) 576–611.
  * Weidauer and Schumacher [2012] T. Weidauer, J. Schumacher,  Moist turbulent rayleigh-bénard convection with neumann and dirichlet boundary conditions,  Physics of Fluids 24 (2012) 076604.
  * Vallis et al. [2019] G. K. Vallis, D. J. Parker, S. M. Tobias,  A simple system for moist convection: the rainy–bénard model,  Journal of Fluid Mechanics 862 (2019) 162–199.
  * Emanuel [1994] K. A. Emanuel, Atmospheric convection, Oxford University Press, USA, 1994.
  * Hartlep and Busse [2018] T. Hartlep, F. H. Busse,  Convection in an internally cooled fluid layer heated from below,  Geophysical & Astrophysical Fluid Dynamics 112 (2018) 20–35.
  * Berlengiero et al. [2012] M. Berlengiero, K. Emanuel, J. Von Hardenberg, A. Provenzale, E. Spiegel,  Internally cooled convection: a fillip for philip,  Communications in Nonlinear Science and Numerical Simulation 17 (2012) 1998–2007.
  * Lepot et al. [2018] S. Lepot, S. Aumaître, B. Gallet,  Radiative heating achieves the ultimate regime of thermal convection,  Proceedings of the National Academy of Sciences 115 (2018) 8937–8941.
  * Burns et al. [2020] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, B. P. Brown,  Dedalus: A flexible framework for numerical simulations with spectral methods,  Physical Review Research 2 (2020) 023068.
  * Emanuel and Bister [1996] K. A. Emanuel, M. Bister,  Moist convective velocity and buoyancy scales,  Journal of the atmospheric sciences 53 (1996) 3276–3285.
  * Chandrasekhar [2013] S. Chandrasekhar, Hydrodynamic and hydromagnetic stability, Courier Corporation, 2013\.