eurm10 msam10

# A playground for compressible natural convection with a nearly uniform
density

Thierry Alboussièrea Email address for correspondence:
<EMAIL_ADDRESS>Jezabel Curbelob Fabien Dubuffeta Stéphane
Labrossea and Yanick Ricarda aUniv. Lyon, Laboratoire de Géologie de Lyon,
CNRS, ENS-Lyon, Université Lyon 1, France
bDepartament de Matemàtiques, Universitat Politècnica de Catalunya, Barcelona,
Spain

(?; revised ?; accepted ?. - To be entered by editorial office)

###### Abstract

In the quest to understand the basic universal features of compressible
convection, one would like to disentangle genuine consequences of compression
from spatial variations of transport properties. Unfortunately this is
sometimes difficult: for instance, one may choose to consider a fluid with
uniform dynamic viscosity, but then, compressible effects will often generate
a density gradient and consequently the kinematic viscosity will not be
uniform. Similarly, a uniform thermal conductivity leads to a gradient of
thermal diffusivity. In the present work, we consider a very peculiar equation
of state, whereby entropy is solely dependent on density, so that a nearly
isentropic fluid domain is nearly isochoric. Within this class of equations of
state, there is a thermal adiabatic gradient and a key property of
compressible convection is still present, namely its capacity to viscously
dissipate a large fraction of the thermal energy involved, of the order of the
well-named dissipation number. In a series of anelastic approximations, under
the assumption of an infinite Prandtl number, the number of governing
parameters can be brought down to two, the Rayleigh number and the dissipation
number. This framework is proposed as a playground for compressible
convection, an opportunity to extend the vast corpus of theoretical analyses
on the Oberbeck-Boussinesq equations regarding stability, bifurcations or the
determination of upper bounds for the turbulent heat transfer. Here, in a two-
dimensional geometry, we concentrate on the structure of upward and downward
plumes depending on the dissipation number, on the heat flux dependence on the
dissipation number and on the ratio of dissipation to convective heat flux.
For all Rayleigh numbers, a change in the vertical temperature profile is
observed in the range of dissipation number between $0$ and less than $0.4$,
associated with the weakening of ascending plumes. For larger dissipation
numbers, the heat flux dependence on this number is found to be well predicted
by Malkus’s model of critical layers. For dissipation numbers of order unity,
in the limit of large Rayleigh numbers, dissipation becomes related to the
entropy heat flux at each depth, so that the vertical dissipation profile can
be predicted, and consequently so does the total ratio of dissipation to
convective heat flux.

###### keywords:

Rayleigh-Bénard, Equations of state, Compressible convection, Anelastic
Approximations

##  1 Introduction

Natural convection is the response of a fluid with a specific equation of
state subject to a thermal or compositional buoyancy forcing – for instance an
imposed temperature difference in a gravity field – while conservation laws of
mass, momentum and energy apply. Compressibility effects are inevitable, but
in a famous approximation due to Oberbeck [1879] and Boussinesq [1903],
pressure effects are relegated to a secondary role. The Oberbeck-Boussinesq
model is so simple and has become so popular that most theoretical studies of
natural convection are made in its framework. We will concentrate on the
Rayleigh-Bénard configuration (Bénard [1901], Rayleigh [1916]), mostly
relevant to stars and planets. In these large natural objects, where
compressibility plays a large role, fewer theoretical results have been
derived and we think this is essentially due to the absence of a simple set of
equations which could be used as a playground for studies of compressible
convection.

In a simple geometry, the Oberbeck-Boussinesq model has just two dimensionless
parameters, the Rayleigh $Ra$ and Prandtl $Pr$ numbers. The Prandtl number is
only relevant to the inertial effects in the momentum equation. In the limit
of infinite Prandtl numbers as it is the case for convection in the solid
mantle of terrestrial planets, this parameter becomes irrelevant, so that
there is a single governing parameter, the Rayleigh number $Ra$. Since the
stability analysis of Rayleigh [1916], a century of theoretical investigations
were led and thousands of scientific papers have been published using the
Oberbeck-Boussinesq model. As soon as compressibility effects are taken into
account, the number of governing parameters jumps to six [see Curbelo et al.,
2019]: $Pr$, $Ra$, $\alpha T$, $c_{p}/c_{v}$, $T_{h}/T_{c}$, $\alpha
gH/c_{p}$, where the symbols $\alpha$, $c_{p}$, $c_{v}$, $T$, $T_{h}$,
$T_{c}$, $g$, $H$ denote the coefficient of thermal expansion, heat capacity
at constant pressure, heat capacity at constant volume, temperature, hot
imposed temperature, cold imposed temperature, gravity and the height of the
fluid layer, respectively. Depending on the equation of state considered,
there can be fewer parameters ($\alpha T=1$ for ideal gases) or more
parameters needed to describe the fluid. This – and numerical difficulties
mentioned below – explain why there are comparatively few studies devoted to
compressible convection and stresses the need to propose simple approaches
that might enable the community to identify basic features of compressible
effects. Hopefully our work will contribute to this objective.

Carnot [1824] was the first to suggest that the low temperature at high
altitude were due to adiabatic decompression of air in ascending currents,
while descending currents and adiabatic compression would bring air back to
the higher temperature at sea level. This was later generalized by
Schwarzschild [1906] for the temperature profile in convective regions of
stars, while Jeffreys [1930] proved that the stability of compressible
convection was governed by the superadiabatic Rayleigh number with the same
threshold (for moderate compressibility) as obtained by Rayleigh [1916] in the
Boussinesq approximation. Later, stability was studied by a number of authors
[Spiegel, 1965, Giterman and Shteinberg, 1970, Bormann, 2001, Busse, 1967,
Paolucci and Chenoweth, 1987, Fröhlich et al., 1992]. More recently, we
published a model of stability valid for any arbitrary equation of state and
uniform dynamic viscosity and conductivity [Alboussière and Ricard, 2017].

A difficulty with the fully compressible governing equations was soon spotted:
they contain the fast sound wave and the convective timescales. In many
instances those timescales are so different that the numerical task of
computing convection is overwhelming. Anelastic approximations (AA) were
developed for the atmosphere, Earth’s core and stars (Ogura and Phillips
[1961], Braginsky and Roberts [1995], Lantz and Fan [1999]), valid in
convective regions, consisting in an expansion about an isentropic state. The
simplified anelastic liquid approximation ALA was proposed in [Anufriev et
al., 2005] in which the role of pressure fluctuations on other thermodynamic
quantities is neglected. In the stably stratified cases, sound-proof models
have also been developed [Durran, 1989, Lipps, 1990, Vasil et al., 2013] in
the pseudo-incompressible approximation. Lecoanet et al. [2014] note that the
pseudo-incompressible equation of state introduces some inaccuracies in the
thermodynamic variables.

When compressible effects are present, there is usually a significant range of
temperatures in the system, since the adiabatic gradient is a key feature of
compressible convection. The same is true for pressure, density, and so on. A
consequence is that transport coefficients of heat or momentum – thermal
conductivity and (dynamic) viscosity – are usually not uniform. It is then
difficult to distinguish between consequences of compressibility and
consequences of non-uniform transport properties. Even in the classical
Boussinesq model can non-uniform transport coefficients be modelled, they are
called the Non-Oberbeck-Boussinesq (NOB) effects, for instance in Horn et al.
[2013]. In the present paper, we try to minimize the NOB effects. For this
reason, we choose uniform constant thermal conductivity and (dynamic)
viscosity. However, when density varies so do kinematic viscosity and thermal
diffusivity. Hence we make a peculiar choice of equation of state, such that
density is constant when entropy is constant $s(\rho)$. This ensures that a
nearly isentropic convective region is also a region of nearly uniform density
and kinematic viscosity. We will see that the heat capacity $c_{p}$ and
thermal diffusivity are also uniform where entropy is uniform.

In the next section 2, we discuss the general validity of an equation of state
and expand the case $s=s(\rho)$. Using that class of equations of state,
section 3 is devoted to the description of the configuration and to writing
the governing equations and anelastic approximations. In section 4 we first
show results of the initial phase of convection from rest, with a small
superadiabatic Rayleigh number and a large dissipation number, in order to
assess the validity of the different anelastic approximation models. We then
show that a significant change in temperature profiles occurs at small
dissipation number, in section 5.1, namely the disappearance of the top
overshoot on the vertical averaged profile. Top and bottom asymmetry is
further studied in section 5.2 for larger values of the dissipation number.
The basic model of critical boundary layer is applied to the compressible case
in section 6 and provides a model for the change of heat flux when the
dissipation number is increased. In section 7, we introduce the expressions
for the vertical heat flux in the different models (fully compressible and
anelastic approximations), as well as that for the dissipation profile, under
a form that will be suited to understand energy transfers in the final
sections. The numerical results of global dissipation relative to the
convective heat flux are shown for all models and a range of superadiabatic
Rayleigh numbers and dissipation numbers in section 7.1. A definite limit is
observed at large superadiabatic Rayleigh numbers which is further studied in
section 7.2. It is interpreted as a local mesoscale equilibrium state whereby
the entropy flux contribution is found to correspond both to energy
dissipation and to the main part of the heat flux. Finally, in section 8, we
consider the additional effects of inertia, cavity aspect ratio and boundary
conditions, to show that another state of flow can be obtained which does not
correspond to that local equilibrium and exhibits larger dissipation. However,
those last boundary conditions with impermeable vertical walls are less
relevant in the geophysical and astrophysical context. In conclusion (section
9), our study gives support to the mesoscale equilibrium implying that the
vertical profile of dissipation takes the form of the function $\alpha
g/c_{p}$ in compressible convection in the limit of large dissipation and
superadiabatic Rayleigh numbers.

Considering the very specific equation of state considered here, the condition
of infinite Prandtl number, the two-dimensional geometry and the absence of
rotation and magnetic effects, the results of this study should not be applied
to stellar or planetary objects without further investigations. However, they
provide a possible asymptotic behaviour for the large compressibility, large
Rayleigh number, limit. It remains to determine under which conditions that
behaviour will be observed.

##  2 Impossible EoS $\rho(T)$ and possible EoS $\rho(s)$

From gases to solids, a wide range of equations of state (EoS in short) is
possible. From a theoretical point of view, one can wonder what should be a
possible EoS and when a tentative EoS is impossible. One answer is that one
should just start from a fundamental EoS under the form

$e=e(s,\rho).$ (1)

where $s$, $e$ and $\rho$ are the specific entropy, specific energy and
density, respectively. From Gibbs equation
$\mathrm{d}e=T\mathrm{d}s+P/\rho^{2}\mathrm{d}\rho$ ($T$ temperature, $P$
pressure), one just needs $T$ to be positive, if we want to consider real
existing conditions

$\left.\frac{\partial e}{\partial s}\right|_{\rho}>0.$ (2)

However, one rarely starts from a fundamental EoS (1). Usually, one expresses
density $\rho$ as a function of pressure $P$ and temperature $T$. The first
obvious idea when one wishes to get rid of compressible effects – and jump
immediately into the Boussinesq approximation – is to state that density is
independent of pressure

$\rho=\rho(T).$ (3)

We now investigate the consequences of this assumption (3) (see also Grandi
and Passerini [2020]). We derive a general relationship, equation (A7) in
Alboussière and Ricard [2013], on the partial derivative of enthalpy
$h=e+P/\rho$ with respect to pressure at constant temperature

$\left.\frac{\partial h}{\partial P}\right|_{T}=\frac{1-\alpha T}{\rho},$ (4)

which is obtained from Gibbs equation under several forms (using the
differential of $h$ and that of Gibbs free energy $g=h-Ts$) and deriving
Maxwell relations. Note that the right-hand side – assuming (3) and hence
$\alpha=-\rho^{\prime}/\rho$ where the prime denotes derivative with respect
to the single variable $T$ – is a function of $T$ only, that we denote $A$

$A=\frac{1}{\rho}+\frac{\rho^{\prime}T}{\rho^{2}},$ (5)

Equation (4) is integrated to give an expression for the enthalpy

$h=AP+B,$ (6)

where $B$ is another function of temperature. This expression is used to write
$\mathrm{d}h$ which is then substituted in Gibbs equation
$\mathrm{d}h=T\mathrm{d}s+\mathrm{d}P/\rho$ leading to

$\mathrm{d}s=\frac{A^{\prime}P+B^{\prime}}{T}\mathrm{d}T+\frac{\rho^{\prime}}{\rho^{2}}\mathrm{d}P.$
(7)

Considering that $A^{\prime}=T(\rho^{\prime}/\rho^{2})^{\prime}$ from (5),
equation (7) implies the following form for $s$

$s=\frac{\rho^{\prime}}{\rho^{2}}P+C,$ (8)

where $C$ is yet another function of $T$ which turns out to be an integral of
$B^{\prime}/T$. We can now write an expression for the Gibbs free energy
$g=h-Ts$

$g=\frac{P}{\rho}+B-TC.$ (9)

A condition of stability of a material substance [Bazarov, 1989] is that its
Gibbs free energy $g$ should be a concave function of $P$ and $T$. If not, the
substance would split into two different phases that have together a larger
entropy, as for instance in the phase-change region of the Van der Waals
model. Locally, a necessary condition is that the Hessian of $g$ (its matrix
of second partial derivatives) is a negative-definite matrix, i.e. has
alternatively negative and positive leading principal minors according to
Sylvester’s criterion [Gilbert, 1991]. The Hessian matrix is

$\left[\begin{array}[]{l l}\left.\frac{\partial^{2}g}{\partial
T^{2}}\right|_{P}&\frac{\partial^{2}g}{\partial T\partial P}\\\
\frac{\partial^{2}g}{\partial P\partial T}&\left.\frac{\partial^{2}g}{\partial
P^{2}}\right|_{T}\end{array}\right]$ (10)

The first leading principal minor is $\partial^{2}g/\partial T^{2}$ at
constant pressure. From Gibbs equation
$\mathrm{d}g=-s\mathrm{d}T+\mathrm{d}P/\rho$, we have

$\left.\frac{\partial g}{\partial T}\right|_{P}=-s,\hskip
56.9055pt\left.\frac{\partial g}{\partial P}\right|_{T}=\frac{1}{\rho},$ (11)

and hence

$\left.\frac{\partial^{2}g}{\partial T^{2}}\right|_{P}=-\left.\frac{\partial
s}{\partial
T}\right|_{P}=-\left(\frac{\rho^{\prime}}{\rho^{2}}\right)^{\prime}P-C^{\prime}=-\left(\frac{\rho^{\prime}}{\rho^{2}}\right)^{\prime}P-\frac{B^{\prime}}{T},$
(12)

from (8), which can indeed be made negative for an appropriate choice of the
function $\rho(T)$ and $B(T)$. Now, the second and last leading principal
minor (in dimension two) is the determinant of the whole Hessian matrix. From
(11), we can see that the second derivative of $g$ with respect to $P$ will be
zero. The determinant of the Hessian matrix is then just equal to

$\det\left[\begin{array}[]{l l}\left.\frac{\partial^{2}g}{\partial
T^{2}}\right|_{P}&\frac{\partial^{2}g}{\partial T\partial P}\\\
\frac{\partial^{2}g}{\partial P\partial T}&\left.\frac{\partial^{2}g}{\partial
P^{2}}\right|_{T}\end{array}\right]=-\left(\frac{\partial^{2}g}{\partial
P\partial T}\right)^{2}=-\left[\left.\frac{\partial}{\partial
T}\left(\frac{1}{\rho}\right)\right|_{P}\right]^{2}.$ (13)

It is negative, meaning that $g$ is not a concave function of $T$ and $P$. The
only way to save that equation of state would be to make this determinant
zero: since it is the partial derivative of $1/\rho$ with respect to $T$ at
constant $P$, it is zero only when $\rho$ is a constant. Such an equation of
state is not interesting for thermal convection as no buoyancy would exist.
Hence we consider equation (3) as an impossible equation of state. Another
related aspect can be noted from Mayer’s relationship

$c_{p}-c_{v}=-\frac{T}{\rho^{2}}\left.\frac{\partial P}{\partial
T}\right|_{\rho}\left.\frac{\partial\rho}{\partial T}\right|_{P}.$ (14)

That difference is infinite since ${\partial P}/{\partial T}$ is infinite at
constant $\rho$, from Euler’s chain rule $\left.\partial P/\partial
T\right|_{\rho}\left.\partial
T/\partial\rho\right|_{P}\left.\partial\rho/\partial P\right|_{T}=-1$. Hence
the choice of meaningful heat capacities is impossible.

We now investigate another simple form of EoS, such that density is a function
of entropy only and show that it satisfies marginally the criterion of
stability, as noted in Scott [2001]. Let us identify all possible equations of
state such that density is solely a function of entropy, or reciprocally such
that entropy is solely a function of density

$s=s(\rho).$ (15)

The thermodynamics Gibbs equation can be written

$\mathrm{d}e=T\mathrm{d}s-P\mathrm{d}v=\left(Ts^{\prime}+\frac{P}{\rho^{2}}\right)\mathrm{d}\rho,$
(16)

where $v$ is the specific volume ($v=1/\rho$) and the primes now denote the
usual derivative with respect to the single variable $\rho$. It follows from
the previous equation that $e$ is also solely a function of $\rho$,

$e=e(\rho).$ (17)

The next consequence, by definition, is that the heat capacity at constant
volume $c_{v}$ is zero

$c_{v}\equiv\left.\frac{\partial e}{\partial
T}\right|_{v}=\left.\frac{\partial e}{\partial T}\right|_{\rho}=0.$ (18)

This shows that our choice is a limit case of valid equations of state, a
negative $c_{v}$ would not be realistic. Instead of considering that entropy
is a function of density only, had we added a tiny dependence on temperature,
we would probably have been able to obtain a strictly positive and small value
for $c_{v}$ and that equation of state would have been perfectly valid.
Equation (16) can also be written

$Ts^{\prime}=e^{\prime}-\frac{P}{\rho^{2}}.$ (19)

Multiplying (19) by $\rho^{2}$ and deriving with respect to temperature $T$ at
constant pressure $P$ leads to the following expression for the coefficient of
thermal expansion

$\alpha\equiv-\frac{1}{\rho}\left.\frac{\partial\rho}{\partial
T}\right|_{P}=\frac{\rho
s^{\prime}}{-(\rho^{2}e^{\prime})^{\prime}+T(\rho^{2}s^{\prime})^{\prime}},$
(20)

Using Gibbs equation (19) to extract $P/\rho$, we express the specific
enthalpy $h$ as follows

$h\equiv e+\frac{P}{\rho}=e+\rho e^{\prime}-\rho Ts^{\prime}.$ (21)

From (21) and (20), after straightforward but slightly tedious steps, we
derive an expression for the heat capacity at constant pressure

${\color[rgb]{0,0,0}c_{p}\equiv\left.\frac{\partial h}{\partial
T}\right|_{P}=-\rho T\alpha s^{\prime}.}$ (22)

At this point, from (20) and (22), we note that $\alpha T/c_{p}$ – which
multiplied by gravity $g$ expresses the adiabatic gradient – is solely a
function of $\rho$

$\frac{\alpha T}{c_{p}}=-\frac{1}{\rho s^{\prime}},$ (23)

so that, in an isentropic region under a uniform gravity field, one can expect
to observe a uniform adiabatic temperature gradient. The condition
$s^{\prime}<0$ is needed to avoid a negative $\alpha$ or worse a negative
$c_{p}$ according to equation (23). However, $c_{p}$ and $\alpha T$ are
functions of $\rho$ and $T$, hence will not be uniform in an isentropic
region. In order to avoid complexity, we assume – in addition to $s$ being a
function of $\rho$ – that $(\rho^{2}e^{\prime})^{\prime}$ is zero, hence

$e=\frac{K}{\rho},$ (24)

up to an irrelevant additive constant, and where the multiplicative constant
$K$ is a parameter whose value can be freely specified. This eliminates the
temperature dependence of $\alpha T$ and $c_{p}$. With equation (24) we have

$\displaystyle\alpha T$ $\displaystyle=\frac{\rho
s^{\prime}}{\left(\rho^{2}s^{\prime}\right)^{\prime}},$ (25) $\displaystyle
c_{p}$
$\displaystyle=-\frac{\rho^{2}{s^{\prime}}^{2}}{\left(\rho^{2}s^{\prime}\right)^{\prime}}.$
(26)

The condition of stability on the leading principal minors of the Hessian
matrix of $g$ is now examined. With our choice for energy (24) and the form of
$h$ in (21), Gibbs free energy $g\equiv h-Ts$ takes the form $g=-T(\rho
s)^{\prime}$.  Using (11), we obtain

$\displaystyle\left.\frac{\partial^{2}g}{\partial T^{2}}\right|_{P}$
$\displaystyle=-\left.\frac{\partial s}{\partial
T}\right|_{P}=-s^{\prime}\left.\frac{\partial\rho}{\partial T}\right|_{P},$
(27) $\displaystyle\left.\frac{\partial^{2}g}{\partial P^{2}}\right|_{T}$
$\displaystyle=-\frac{1}{\rho^{2}}\left.\frac{\partial\rho}{\partial
P}\right|_{T},$ (28) $\displaystyle\frac{\partial^{2}g}{\partial T\partial P}$
$\displaystyle=-\left.\frac{\partial s}{\partial
P}\right|_{T}=-s^{\prime}\left.\frac{\partial\rho}{\partial P}\right|_{T}.$
(29)

In order to evaluate these second derivatives, we need expressions for the
partial derivatives of density with respect to temperature and density. From
(19), we get

$\left.\frac{\partial
P}{\partial\rho}\right|_{T}=\left(\rho^{2}e^{\prime}\right)^{\prime}-T\left(\rho^{2}s^{\prime}\right)^{\prime}=-T\left(\rho^{2}s^{\prime}\right)^{\prime},$
(30)

owing to our choice $\left(\rho^{2}e^{\prime}\right)^{\prime}=0$. The inverse
of (30) provides ${\partial\rho}/{\partial P}$ while (25) is used to express
${\partial\rho}/{\partial T}$. When substituted in (27), (28) and (29), we
obtain the Hessian matrix

$\left[\begin{array}[]{l l}\left.\frac{\partial^{2}g}{\partial
T^{2}}\right|_{P}&\frac{\partial^{2}g}{\partial T\partial P}\\\
\frac{\partial^{2}g}{\partial T\partial P}&\left.\frac{\partial^{2}g}{\partial
P^{2}}\right|_{T}\end{array}\right]=\left[\begin{array}[]{l
l}\frac{\rho^{2}{s^{\prime}}^{2}}{T\left(\rho^{2}s^{\prime}\right)^{\prime}}&\frac{s^{\prime}}{T\left(\rho^{2}s^{\prime}\right)^{\prime}}\\\
\frac{s^{\prime}}{T\left(\rho^{2}s^{\prime}\right)^{\prime}}&\frac{1}{\rho^{2}T\left(\rho^{2}s^{\prime}\right)^{\prime}}\end{array}\right],$
(31)

The first leading principal minor is negative when
$\left(\rho^{2}s^{\prime}\right)^{\prime}<0$, so that this condition must be
fulfilled. The second minor is the whole determinant of the Hessian matrix and
it is easy to check that it is zero. In that sense the equation of state
$s(\rho)$ is just marginally stable.

There is still a large set of possibilities since we are free to consider any
function $s(\rho)$, provided $\left(\rho^{2}s^{\prime}\right)^{\prime}<0$ and
$s^{\prime}<0$ if one wishes to restrict the analysis to positive values of
$\alpha$. Let us choose a set of such decreasing functions, defined as one of
the following up to an irrelevant additive constant

$\displaystyle s(\rho)$ $\displaystyle=a\ln(\rho),$ (32)
$\displaystyle\mathrm{or}\ \ s(\rho)$ $\displaystyle=a\rho^{n},\hskip
28.45274pt\mathrm{for}\ n>0\ \mathrm{a\ real\ number},$ (33)

where $a<0$ is a negative constant real parameter. With the logarithm function
$s\sim\ln(\rho)$, we have $\alpha T=1$ and a constant $c_{p}=-a$. With a power
law $s\sim\rho^{n}$, we have a constant $\alpha T$ between 1 and 0 whose value
can be tuned by choosing the positive exponent $n$ and $c_{p}$ is a function
of $\rho$

$\alpha T=\frac{1}{n+1},\hskip 28.45274ptc_{p}=-a\frac{n}{n+1}\rho^{n},\hskip
28.45274pt\mathrm{for}\ n>0.$ (34)

Other relations will be needed, namely the expressions of $P$ and $h$

$\displaystyle P$ $\displaystyle=-K-aT\rho,\hskip 28.45274pth=-aT,$ (35)
$\displaystyle\mathrm{or}\ \ P$ $\displaystyle=-K-naT\rho^{n+1},\hskip
28.45274pth=-naT\rho^{n},\hskip 28.45274pt\mathrm{for}\ n>0.$ (36)

We first remark that one of these equations of state is an ideal gas equation:
this is the case of (32) when $K=0$, $a=-c_{p}$ and corresponds to an ideal
gas with $c_{v}=0$. The marginal stability of this equation of state is
reflected by the infinite speed of sound that results from a finite $c_{p}$
and a null $c_{v}$.

Although our equation of state was built from theoretical arguments, one may
try to find real substances with a similar behaviour, at least in some range
of temperature and pressure: a monoatomic gas with large molar mass has a
small $c_{v}$ for instance. Radon gas is a good example. Next, the ratio
$c_{p}/c_{v}$ can be made large (diverging to infinity) near the critical
point, so that radon near the critical point would have the expected behaviour
concerning heat capacities. However, the thermal expansion coefficient also
diverges near the critical point and that does not match our equation of
state.

Among the large class of equations of state such that entropy is a function of
density, driven by a principle of simplicity, we have identified a set of such
equations, with $\alpha T$ constant ranging from $1$ (log function) to zero
asymptotically (power law with $n\rightarrow\infty$). For all of them, $c_{p}$
and the expected adiabatic gradient are functions of density only.

##  3 Rayleigh-Bénard configuration and governing equations

We define the geometric configuration and boundary conditions that will be
investigated in this paper. Different convection models are considered :
complete continuum thermodynamic and dynamic equations (FC for ’fully
compressible’), anelastic approximation (AA), anelastic liquid approximation
(ALA) and a further simplified model (SCA for ’simple compressible
approximation’).

###  3.1 Fully compressible model

For simplicity, we take the infinite Prandtl number approximation which
eliminates inertia. This limit has been studied mathematically [Wang, 2004]
and used for the study of mantle dynamics [Ricard, 2015] for which Prandtl
numbers are estimated around $10^{25}$: the effective kinematic viscosity of
solids is much larger than their thermal diffusivity. Other objects, like the
Earth’s core, the interior of stars and of gaseous planets have low Prandtl
numbers [Schaeffer et al., 2017, Garaud, 2020, Fuentes and Cumming, 2020]. For
infinite Prandtl numbers, the governing equations of thermal convection are
the following

$\displaystyle\frac{\mathrm{D}\rho}{\mathrm{D}t}$
$\displaystyle=-\rho{\boldsymbol{\nabla}}\cdot{\bf u},$ (37)
$\displaystyle{\bf 0}$
$\displaystyle=-{\mathbf{\boldsymbol{\nabla}}}P+\rho{\bf
g}+\eta\left[{\mathbf{\boldsymbol{\nabla}}}^{2}{\bf
u}+\frac{1}{3}{\mathbf{\boldsymbol{\nabla}}}\left(\boldsymbol{\nabla}\cdot{\bf
u}\right)\right],$ (38) $\displaystyle\rho T\frac{\mathrm{D}s}{\mathrm{D}t}$
$\displaystyle=\dot{\epsilon}:\tau+\boldsymbol{\nabla}\cdot\left(k{\mathbf{\boldsymbol{\nabla}}}T\right),$
(39)

where ${\bf u}$ is the velocity field, ${\bf g}$ the gravity field, $\eta$ the
dynamic viscosity of the fluid, $k$ its thermal conductivity,
$\mathrm{D}/\mathrm{D}t=\partial/\partial t+{\bf
u}\cdot{\mathrm{\boldsymbol{\nabla}}}$ is the material derivative,
$\dot{\epsilon}$ denote the tensor of rate of deformation, $\tau$ the
Newtonian stress tensor defined as

$\tau_{ij}=2\eta\left[\dot{\epsilon}_{ij}-\frac{1}{3}\left({\mathbf{\boldsymbol{\nabla}}}\cdot{\bf
u}\right)\delta_{ij}\right],$ (40)

using Stokes’ assumption regarding the bulk viscosity.

We consider a two-dimensional rectangular domain, with horizontal periodic
boundary conditions. In a Cartesian frame $(x,z)$, the horizontal axis is $x$
while the vertical axis is $z$. The height of the cavity is $H$, its length
$L$. The aspect ratio is set to $L/H=4\sqrt{2}$, corresponding to twice the
horizontal period of the stability analysis in the Boussinesq approximation
for an infinite layer. Gravity is uniform ${\bf g}=-g{\bf e}_{z}$ along the
direction of the unit vertical vector ${\bf e}_{z}$. The thermal boundary
conditions are that of a hot temperature $T_{h}$ at the bottom and a cold
temperature $T_{c}$ at the top. At the top and bottom boundaries, the normal
velocity component is zero and so is the tangential viscous stress, i.e.
$\partial u_{x}/\partial z=0$. Because there is no natural constraint on the
horizontal velocity, we impose that the horizontal average of $u_{x}$ is zero
on the top boundary. Finally, instead of imposing some pressure value, we
impose that the average density in the domain is $\rho_{0}$

$\frac{1}{HL}\int\rho\ \mathrm{d}x\mathrm{d}z=\rho_{0}.$ (41)

This condition is an initial condition and that integral cannot change in time
with impermeable or periodic boundaries.

Figure 1: Geometry and boundary conditions. A typical vertical temperature
profile is sketched (solid line) with an adiabatic, isentropic, temperature
profile (dashed line).

The set of equations is complete when an equation of state is specified. In
this paper, as we consider the class of EoS such that entropy is a function of
density (15), we have
$\mathrm{D}s/\mathrm{D}t=s^{\prime}\mathrm{D}\rho/\mathrm{D}t$. Using the
continuity equation (37), equation (39) can be written in the following form

$-\rho^{2}Ts^{\prime}{\mathbf{\boldsymbol{\nabla}}}\cdot{\mathbf{u}}=\dot{\epsilon}:\tau+\boldsymbol{\nabla}\cdot\left(k{\mathbf{\boldsymbol{\nabla}}}T\right),$
(42)

which is now an elliptic equation for temperature. When $P$ is expressed in
terms of $T$ and $\rho$ – see equations (35) or (36) – the Stokes’s equation
(38) also becomes a Poisson equation for velocity (along with the continuity
equation). By the way, it is already clear that the constant $K$ in the
expression for the internal energy (24) and in that for pressure (35) or (36)
is completely irrelevant in the governing equations for convection: internal
energy does not appear explicitly and taking the gradient of pressure
eliminates $K$ from the momentum equation (38).

The next step consists in defining dimensional scales and in writing the
equations in dimensionless form. We have already mentioned a scale for
density, $\rho_{0}$ which is the average density in the domain that remains
constant with the imposed boundary conditions. Next, we define
$T_{0}=(T_{h}+T_{c})/2$ the average temperature of the hot and cold
boundaries. Then, we need to choose either a $\log$ or power law EoS along
with an exponent $n$. We specify $c_{p0}$ the value of $c_{p}$ at the
conditions $T=T_{0}$ and $\rho=\rho_{0}$, which is equivalent to specifying
the constant $a$. From the logarithmic equation $s\sim\log\rho$, we have
$c_{p0}=-a$ while for the power law $s\sim\rho^{n}$ and (34), we have

$c_{p0}=-a\frac{n}{n+1}\rho_{0}^{n},\hskip 28.45274pt\mathrm{for}\ n>0.$ (43)

From (32) and (33), we derive an expression for $s^{\prime}$ which is valid
for both the logarithmic ($n=0$) and power-law ($n>0$) EoS

$s^{\prime}(\rho)=-(n+1)\frac{c_{p0}}{\rho_{0}}\left(\frac{\rho}{\rho_{0}}\right)^{n-1},\hskip
28.45274pt\mathrm{for}\ n\geq 0.$ (44)

Similarly, a generic expression is obtained for the pressure gradient, from
(35) and (36)

${\mathbf{\boldsymbol{\nabla}}}P=(n+1)c_{p0}\left(\frac{\rho}{\rho_{0}}\right)^{n}\left[(n+1)T{\mathbf{\boldsymbol{\nabla}}}\rho+\rho{\mathbf{\boldsymbol{\nabla}}}T\right],\hskip
28.45274pt\mathrm{for}\ \ n\geq 0.$ (45)

We consider a uniform thermal conductivity $k$, so that a scale for thermal
diffusivity is $\kappa=k/(\rho_{0}c_{p0})$. We now make all variables
dimensionless using $H$, $\kappa/H$, $H^{2}/\kappa$, $T_{0}$, $\rho_{0}$,
$c_{p0}$, $\rho_{0}c_{p0}T_{0}$, $\kappa/H^{2}$, $\eta\kappa/H^{2}$, for
length, velocity, time, temperature, density, entropy, pressure, deformation
rate and stress. Using the same symbols for dimensionless variables, the
equations of continuity, momentum and entropy become

$\displaystyle\frac{\mathrm{D}\rho}{\mathrm{D}t}$
$\displaystyle=-\rho{\boldsymbol{\nabla}}\cdot{\bf u},$ (46)
$\displaystyle{\mathbf{0}}$
$\displaystyle=-\frac{Ra_{sa}(n+1)}{\varepsilon}\left(\frac{\rho^{n}}{\mathcal{D}}\left[(n+1)T{\mathbf{\boldsymbol{\nabla}}}\rho+\rho{\mathbf{\boldsymbol{\nabla}}}T\right]+\rho{\mathbf{e}}_{z}\right)+{\mathbf{\boldsymbol{\nabla}}}^{2}{\bf
u}+\frac{1}{3}{\mathbf{\boldsymbol{\nabla}}}\left(\boldsymbol{\nabla}\cdot{\bf
u}\right),$ (47) $\displaystyle 0$
$\displaystyle=-(n+1)\rho^{n+1}T\boldsymbol{\nabla}\cdot{\bf
u}+\frac{\varepsilon\mathcal{D}}{Ra_{sa}}\dot{\epsilon}:\tau+\nabla^{2}T,$
(48)

where the following dimensionless numbers appear, the superadiabatic Rayleigh
number $Ra_{sa}$, the dissipation number $\mathcal{D}$, the ratio of
superadiabatic temperature difference over the average temperature
$\varepsilon$ and implicitly the product $\alpha_{0}T_{0}$, as a function of
$n$:

$\displaystyle Ra_{sa}$ $\displaystyle=\frac{\rho_{0}g\alpha_{0}\Delta
T_{sa}H^{3}}{\eta\kappa},$ (49) $\displaystyle\mathcal{D}$
$\displaystyle=\frac{\alpha_{0}gH}{c_{p0}},$ (50) $\displaystyle\varepsilon$
$\displaystyle=\frac{\Delta T_{sa}}{T_{0}},$ (51)
$\displaystyle\alpha_{0}T_{0}$ $\displaystyle=\frac{1}{n+1}.$ (52)

The dissipation number $\mathcal{D}$ is one of the possible measures for
compressibility, of the same nature as the number of scale heights in
astrophysics [Spiegel and Veronis, 1971]. It was introduced by Gebhart [1962],
motivated by the context of cooling turbine blades by natural convection.
Interestingly, the dissipation number can be defined in the framework of the
Boussinesq approximation, although compressibility is absent and despite the
fact that its value has no effect on the solutions. Moreover, it can be shown
rigorously from the Boussinesq equations that the integral of viscous
dissipation is equal to the product of the dissipation number $\mathcal{D}$
and the convective heat flux in a Rayleigh-Bénard cavity [Howard, 1963].  The
superadiabatic temperature difference $\Delta T_{sa}$ is equal to the
difference between the imposed hot and cold temperatures minus the temperature
difference along the adiabat

$\Delta T_{sa}=T_{h}-T_{c}-\frac{\alpha_{0}gT_{0}H}{c_{p0}}.$ (53)

When writing the dimensionless momentum equation (47), we use (45) to express
the pressure gradient in terms of density and temperature gradient. When
writing the dimensionless entropy equation (48), we use (42) and (44). It can
be checked that the final set of dimensionless equations (46), (47) and (48)
takes a generic form for any real value of $n\geq 0$: the case $n=0$
corresponds to the logarithmic relationship (32) while the cases $n>0$
correspond to the power laws (33). The choice of $n$ amounts to choosing the
product $\alpha_{0}T_{0}$, see (52).

As initial conditions, we set the velocity to zero, and density, pressure and,
temperature fields satisfying the (potentially unstable) hydrostatic
conduction regime, with an additional random temperature field of magnitude
$10^{-6}$. The boundary conditions on the velocity and temperature fields are
the following

$\displaystyle\frac{\partial u_{x}}{\partial z}$ $\displaystyle=0,\hskip
28.45274pt\mathrm{when}\hskip 28.45274ptz=\pm\frac{1}{2},$ (54) $\displaystyle
u_{z}$ $\displaystyle=0,\hskip 28.45274pt\mathrm{when}\hskip
28.45274ptz=\pm\frac{1}{2},$ (55)
$\displaystyle{\color[rgb]{0,0,0}\int_{-L/(2H)}^{L/(2H)}u_{x}\left(x,z=\frac{1}{2}\right)\mathrm{d}x=0,}$
(56) $\displaystyle{T}$ $\displaystyle=\frac{T_{h}}{T_{0}},\hskip
28.45274pt\mathrm{when}\hskip 28.45274ptz=-\frac{1}{2},$ (57)
$\displaystyle{T}$ $\displaystyle=\frac{T_{c}}{T_{0}},\hskip
28.45274pt\mathrm{when}\hskip 28.45274ptz=\frac{1}{2}.$ (58)

The stress-free, non-penetrative conditions (54) and (55) do not constrain the
mean horizontal velocity, hence an arbitrary condition of zero average
horizontal velocity (56) is imposed on the upper boundary. The imposed
temperature ratios are linked to the values of the dissipation number and the
superadiabatic temperature coefficient

$\displaystyle\frac{T_{h}}{T_{0}}$
$\displaystyle=1+\frac{\mathcal{D}+\varepsilon}{2},$ (59)
$\displaystyle\frac{T_{c}}{T_{0}}$
$\displaystyle=1-\frac{\mathcal{D}+\varepsilon}{2}.$ (60)

As shown in Curbelo et al. [2019], the equations of convection with infinite
Prandtl number are subjected to viscous relaxation, and the associated
relaxation time limits the time-steps to $\mathcal{D}/{Ra_{sa}}$ for numerical
calculations. Let us determine here the expression for this relaxation time
scale, for our particular class of EoS. We consider a small planar disturbance
with respect to the steady solution $(\rho=1,\ T=1-\mathcal{D}z,\ {\bf u}={\bf
0})$ with $\epsilon=0$

$\rho^{\prime}=\tilde{\rho}e^{ikx+\omega t},\hskip
28.45274ptT^{\prime}=\tilde{T}e^{ikx+\omega t},\hskip
28.45274ptu_{x}^{\prime}=\tilde{u}_{x}e^{ikx+\omega t},$ (61)

where $\tilde{\rho}$, $\tilde{T}$ and $\tilde{u}_{x}$ are scalars. The
governing equations (46), (47) and (48), are linearized near $z=0$ (the steady
solution is nearly constant $T=1$) and lead to

$\displaystyle\omega\tilde{\rho}$ $\displaystyle=-ik\tilde{u}_{x},$ (62)
$\displaystyle 0$
$\displaystyle=-\frac{Ra_{sa}(n+1)}{\varepsilon\mathcal{D}}\left[(n+1)ik\tilde{\rho}+ik\tilde{T}\right]-\frac{2}{3}k^{2}\tilde{u}_{x},$
(63) $\displaystyle 0$ $\displaystyle=-(n+1)ik\tilde{u}_{x}-k^{2}\tilde{T}.$
(64)

Eliminating $\tilde{u}_{x}$ and $\tilde{T}$, leads to a single equation for
$\tilde{\rho}$

$0=-\frac{Ra_{sa}(n+1)^{2}}{\varepsilon\mathcal{D}}\left[k\tilde{\rho}+\frac{\omega}{k}\tilde{\rho}\right]-\frac{2}{3}k\omega\tilde{\rho},$
(65)

admitting non-trivial solutions when

$\omega=-\frac{3}{2}\frac{Ra_{sa}(n+1)^{2}}{\varepsilon\mathcal{D}}\left[1+\frac{\omega}{k^{2}}\right],$
(66)

implying that the magnitude of the rate of decay $|\omega|$ is bounded from
above as

${\color[rgb]{0,0,0}|\omega|<}\frac{3}{2}\frac{Ra_{sa}(n+1)^{2}}{\varepsilon\mathcal{D}},$
(67)

irrespective of the wavenumber $k$. In practice, we make it slightly safer by
changing the prefactor from $3/2$ to $1$, and our numerical scheme is always
found to be stable with time-steps $\delta t$ smaller than

$\delta
t\leq\frac{\varepsilon\mathcal{D}}{Ra_{sa}(n+1)^{2}}=\frac{\varepsilon\mathcal{D}\left(\alpha_{0}T_{0}\right)^{2}}{Ra_{sa}}.$
(68)

This makes it difficult to calculate flows with large superadiabatic Rayleigh
numbers, small dissipation numbers, small superadiabatic parameters
$\varepsilon$ or small products $\alpha T_{0}$ (large $n$).

We are now going to write a series of anelastic models from the most to the
least faithful approximation of the fully compressible equations.

###  3.2 Anelastic Approximation AA

The first model is called simply the anelastic model and corresponds to the
early models by Ogura and Phillips [1961] for the atmosphere, Lantz and Fan
[1999] for stellar convection and Braginsky and Roberts [1995] for the Earth’s
core. It corresponds to a first order expansion modelling of thermodynamic
variables with respect to a hydrostatic isentropic state. In our case, the
structure of the well-mixed isentropic region is simple, with a uniform
density and uniform temperature gradient: in dimensionless form

$\displaystyle\rho_{a}$ $\displaystyle=1,$ (69) $\displaystyle T_{a}$
$\displaystyle=1-\mathcal{D}z,$ (70)

where $T_{a}$ is the isentropic profile. We have set arbitrarily $T_{a}=1$ at
$z=0$ (mid-height) but we will see later that this does not constrain the
anelastic solution. Let us denote with tildes the two-dimensional and time-
dependent departures of each variable from its isentropic counterpart. From
the standard procedure of linearization of the functions of state about the
adiabatic profile [Anufriev et al., 2005], and with a change in the
dimensional scale for temperature ($\Delta T_{sa}$ instead of $T_{0}$),
pressure ($\rho_{0}c_{p0}\Delta T_{sa}$ instead of $\rho_{0}c_{p0}T_{0}$) and
entropy ($c_{p0}\Delta T_{sa}/T_{0}$ instead of $c_{p0}$), we get the
following dimensionless anelastic equations

$\displaystyle{\mathbf{\boldsymbol{\nabla}}}\cdot{\bf u}$ $\displaystyle=0,$
(71) $\displaystyle{\bf 0}$
$\displaystyle=-\frac{Ra_{sa}}{\mathcal{D}}{\mathbf{\boldsymbol{\nabla}}}\tilde{P}+Ra_{sa}\tilde{s}{\hat{\mathbf{e}}}_{z}+{\mathbf{\boldsymbol{\nabla}}}^{2}{\bf
u},$ (72) $\displaystyle\frac{\mathrm{D}}{\mathrm{D}t}(T_{a}\tilde{s})$
$\displaystyle=-\mathcal{D}u_{z}\tilde{s}+\frac{\mathcal{D}}{Ra_{sa}}\dot{\epsilon}:\tau+\nabla^{2}\tilde{T}.$
(73)

For our equation of state, $\tilde{s}$ is proportional to $\tilde{\rho}$ – see
(44) – and linearizing equations (35) or (36) leads to

$\displaystyle\tilde{s}$ $\displaystyle=-(n+1)\tilde{\rho},$ (74)
$\displaystyle\tilde{\rho}$
$\displaystyle=-\frac{\tilde{T}}{(n+1)T_{a}}+\frac{\tilde{P}}{(n+1)^{2}T_{a}},$
(75)

and therefore

$\tilde{s}=\frac{\tilde{T}}{T_{a}}-\frac{\tilde{P}}{(n+1)T_{a}}.$ (76)

Because of the new temperature scale $\Delta T_{sa}$, the temperature boundary
conditions become

$\tilde{T}\left(z=\pm\frac{1}{2}\right)=\mp\frac{1}{2}.$ (77)

The boundary conditions for pressure are obtained from the condition of mass
conservation. With our choice in (69), the (uniform) adiabatic density profile
corresponds already to the total mass in the fluid layer, the integral of the
departure $\tilde{\rho}$ must be zero at all times. This might be achieved by
imposing an appropriate value of pressure at the top or at the bottom of the
cavity. However, an easier way is to impose that the mean value of $\tilde{P}$
on the top boundary is equal to the mean value at the bottom. This can be seen
on equation (72) by integration along $z$. The integral of density in the
cavity is zero, the viscous term ${\mathbf{\boldsymbol{\nabla}}}^{2}{\bf u}$
in (72) integrates into the difference of averaged viscous traction
$\tau_{zz}=2\partial u_{z}/\partial z$ between top and bottom boundaries. The
continuity equation (71) leads to $\tau_{zz}=-2\partial u_{x}/\partial x$
whose integral of each boundary is zero with periodic conditions on $x$. The
condition on pressure is thus

$\int_{0}^{L/H}\tilde{P}\left(x,z=\frac{1}{2}\right)-\tilde{P}\left(x,z=-\frac{1}{2}\right)\mathrm{d}x=0.$
(78)

An invariance property of these AA equations can be put in evidence. Since the
anelastic equations have been obtained by linearization around the adiabatic
profile, one expects that a shift in the superadiabatic temperature conditions
should leave the solution unchanged, with the same total mass. From a
(possibly time-dependent) solution $({\bf u},\tilde{P},\tilde{T})$ to the
equations above, we just add a constant $c$ to the temperature boundary
conditions, now becoming $\tilde{T}(z=\pm 1/2)=\mp 1/2+c$. We can check that
$({\bf u},\tilde{P}+c/(n+1),\tilde{T}+c)$ is a solution to the AA equations
with the shifted temperature boundary conditions.

###  3.3 Anelastic Liquid Approximation ALA

In that approximation, departures of entropy from the adiabatic profile are
considered to be due only to temperature departures, while departures in
pressure are neglected in equation 76 [Anufriev et al., 2005]. The governing
equations are still (71), (72) and (73) where $\tilde{s}$ is changed in each
instance into

$\tilde{s}=\frac{\tilde{T}}{T_{a}}.$ (79)

When applying the boundary conditions, we now find that it is not possible to
impose the obvious temperature boundary condition (77): if we did so, it is
most likely that the integral of $\tilde{T}/T_{a}$ over the fluid domain would
not be zero. However, because of the anelastic liquid approximation, entropy
departures are linked to $\tilde{T}/T_{a}$ which are themselves proportional
to density departures, because of the equation of state. Hence a non-zero
integral of $\tilde{T}/T_{a}$ implies that the total mass of the fluid is not
conserved at first order. Thus, we keep the condition (78) on pressure, which
ensures total mass conservation and we consider that only the imposed
temperature difference is meaningful between two isothermal boundaries

$\displaystyle\frac{\partial\tilde{T}}{\partial x}=0\hskip
28.45274ptz=\pm\frac{1}{2},$ (80)
$\displaystyle\tilde{T}\left(x,z=-\frac{1}{2}\right)-\tilde{T}\left(x,z=\frac{1}{2}\right)=1.$
(81)

Coming back to the pressure constraint (78), any additive constant to
$\tilde{P}$ is irrelevant since only the gradient of pressure plays a role in
the ALA equations. In conclusion, we may decide to set the mean value of
$\tilde{P}$ to zero (or any other constant) on both hot and cold boundaries.
Equation (78) is changed for

$\int_{0}^{L/H}\tilde{P}\left(x,z=\pm\frac{1}{2}\right)\mathrm{d}x=0.$ (82)

The invariance mentioned for the solutions to the AA equations is no longer
relevant in the ALA equations. Now, the mass balance imposes that the integral
of $T/T_{a}$ must be zero over the whole domain because density fluctuations
are solely functions of entropy fluctuations, which are themselves solely
functions of temperature fluctuations in the ALA approximation. That cannot be
changed by another choice of pressure offset.

###  3.4 Simple Compressible Approximation SCA

We now introduce a new approximation aiming at getting a very simple system of
equations where compressible effects are still present. In the ALA
approximation, the adiabatic temperature profile appears explicitly in the
equations and we consider replacing $T_{a}$ by a constant value equal to $1$,
the value of $T_{a}$ in the mid-plane of the cavity. We certainly lose
connection to thermodynamics with that move, but compressible work is still
present and it will be interesting to investigate which compressible effects
are still well accounted for in this approximation. Note that this SCA model
is equivalent to a version of the ”Extended Boussinesq Approximation”, or EBA
[King et al., 2010], where the background density is assumed to be uniform
(with our EoS, this is the case of all our anelastic models) and where the
background temperature is also assumed to be uniform. The SCA model is still
defined by equations (71), (72) and (73), where the expression for entropy
(79) is changed for

$\tilde{s}=\tilde{T},$ (83)

and $T_{a}$ is also changed for the constant value $1$ is the left-hand side
term of (73). Under this approximation, the equations are very similar to the
classical Boussinesq equations, except for viscous heating and adiabatic
cooling that play a significant role when the dissipation number is of order
one.

The boundary conditions are similar to those for the ALA equations. Now, the
average of the temperature departure $\tilde{T}$ on the whole domain is zero
thanks to the condition on pressure (82). An important invariance is valid
only in the case $\mathcal{D}=0$. This corresponds to the Boussinesq equations
with another change in pressure scale, from $\rho_{0}c_{p0}\Delta T_{sa}$ to
$\alpha_{0}\Delta T_{sa}\rho_{0}gH$. In that case only ($\mathcal{D}=0$), the
equations are invariant by symmetry about the mid-plane. More precisely, if
$(u_{x},u_{z},\tilde{P},\tilde{T})$ is a solution (possibly time-dependent),
then the fields
$(u_{x}(x,-z,t),-u_{z}(x,-z,t),\tilde{P}(x,-z,t),-\tilde{T}(x,-z,t))$
constitute another solution. This implies that the solutions to the
incompressible Rayleigh-Bénard system are bottom-up invariant: ascending and
descending plumes are statistically symmetrical. However, when
$\mathcal{D}\neq 0$, that invariance does no longer hold, for none of the
compressible models presented here, fully compressible (FC), anelastic
approximation (AA), anelastic liquid approximation (ALA) nor the last simple
compressible approximation (SCA). We will have the opportunity to investigate
that non-invariance in the following sections.

All numerical results, FC, AA, ALA and SCA have been obtained with the
software Dedalus [Burns et al., 2020]. For the fully compressible model a
Runge-Kutta model of order 1 was used (RK111 in dedalus) and of order 4 for
the anelastic models (RK443 in dedalus). The number of Fourier modes in the
horizontal direction $n_{x}$ is four times that of Chebyshev modes $n_{z}$ in
the vertical direction. A dealiasing factor of $3/2$ has been used in all
cases. The value of $n_{z}$ we used goes from $32$ at low superadiabatic
Rayleigh numbers to $512$ at $Ra_{sa}=10^{9}$. Time steps have been set by
half the short viscous relaxation time in the fully compressible model and by
a Courant condition for anelastic models with a safety factor set to $0.9$.
Noise on the initial conduction temperature field, of magnitude $10^{-6}$ has
been added in all models to trigger convection.

The complete set of equations FC, AA, ALA and SCA, with boundary conditions,
are written in their explicit forms in appendix A. The parameters of all
simulations are listed in appendix B. The files used for each convection model
FC, AA, ALA and SCA are provided as supplemental materials or available at
github.

##  4 From rest to steady rolls at $Ra_{sa}=10^{4}$, $\mathcal{D}=1.5$

In this section, we just analyze the transition from an unstable
superadiabatic motionless initial state to steady rolls of convection, for a
moderate superadiabatic Rayleigh number of $Ra_{sa}=10^{4}$ and a large
dissipation number of $\mathcal{D}=1.5$. We do that for different values of
the dimensionless parameter $\alpha_{0}T_{0}=1$, $0.5$ and $0.1$ (with $n=0$,
1 and 9) and the different models of convection: FC, AA, ALA and SCA. In
figure 2, we plot the upper and lower heat fluxes (on the top and bottom
boundaries) for all values of $\alpha_{0}T_{0}$ and all models. Only the heat
fluxes above the conduction flux along the adiabatic gradient are represented:
this is straightforward in the anelastic model as the conduction flux along
the adiabat is not computed, while for the fully compressible model we just
remove the contribution of conduction along the adiabatic gradient. Then, the
remaining part of the flux is scaled by the conduction heat flux driven by
$\Delta T_{sa}$, the superadiabatic temperature difference: again, this is
natural in the anelastic models where temperature intervals are already scaled
by $\Delta T_{sa}$, while in the fully compressible model the temperature
scale is $T_{0}$ and the flux has to be rescaled by $\Delta T_{sa}$
corresponding to a division by the superadiabatic fraction $\varepsilon=\Delta
T_{sa}/T_{0}$. In the present FC calculations, the superadiabatic fraction
$\varepsilon$ is set to $0.1$. The initial (unstable) conduction state
corresponds to a heat flux unity, while when a convective steady state of
convection is reached, the heat flux is around 5 or slightly less. This
dimensionless flux is the classical Nusselt number, which will also be used in
the next sections. The blue/green colors correspond to the heat flux at the
upper boundary, while red/purple colors correspond to the heat flux on the
lower boundary (for all values of $\alpha_{0}T_{0}$ and FC, AA, ALA). The
simple compressible model (SCA) is plotted with a black color: it can be shown
easily that upper and lower heat fluxes coincide at all times for this model.
All plots have been shifted in time so that the maximum upper flux is at
$t=0$. The real starting time of the simulation depends on the model
considered and is made visible by a small transient period, as we have started
the simulations from our approximations of the conductive hydrostatic state.
The time needed to develop the convective instability depends on the model of
convection, and weakly on $\alpha_{0}T_{0}$ for FC. For all models FC, AA and
ALA, the curves of upper and lower heat fluxes are rather similar, we will
come back to the small differences below. In all cases, in the beginning of
the convective instability, the heat flux on the upper boundary grows rapidly
to a large value (around 11), while the heat flux on the lower boundary
decreases rapidly to negative values (around -4). Then follows a series of
oscillations of decreasing amplitude, with a phase shift of approximately
$\pi/2$ between upper and lower fluxes, until a steady state is reached with
equal upper and lower fluxes (around 5 or slightly less).

Figure 2: Upper and lower heat flux (Nusselt number) during the initial
transient from rest to steady rolls, at $Ra_{sa}=10^{4}$ and
$\mathcal{D}=1.5$.

Figure 3: Difference between the heat flux (Nusselt number) on the upper
boundary obtained with FC and AA, ALA and SCA respectively, at
$Ra_{sa}=10^{4}$ and $\mathcal{D}=1.5$, for $n=0$ ($\alpha_{0}T_{0}=1$), $n=1$
($\alpha_{0}T_{0}=0.5$) and $n=9$ ($\alpha_{0}T_{0}=0.1$) from left to right.

On Fig. 3, we take the difference on the upper flux between the FC model and
the approximations AA, ALA and SCA. Unsurprisingly, the smallest difference is
obtained with the AA approximation, followed by the ALA and finally the SCA
approximation. We also observe that the difference between AA and ALA
approximations becomes smaller as the product $\alpha_{0}T_{0}$ is decreased.
This was expected from equation (76) as the effect of pressure is divided by
$n+1$, i.e. decreases with $\alpha_{0}T_{0}$.

##  5 Top/bottom asymmetry

The top/bottom symmetry is observed to hold for all models in the limit of
vanishing compressibility effect $\mathcal{D}\longrightarrow 0$. In the case
of the fully compressible model (FC) one must also have a top/bottom
temperature ratio close to one, but that condition has a small effect on the
asymmetry compared to that of the dissipation parameter. However, when
$\mathcal{D}$ is non-zero we will see that a distinct difference appears
between the top and bottom parts of the average temperature profile, or
between raising and descending plumes. Perhaps surprisingly, the asymmetry
becomes very clear from relatively small values of the dissipation number,
$\mathcal{D}\sim 0.1$ (section 5.1 below), and continues to exist when
$\mathcal{D}$ is further increased (section 5.2 below). Also, increasing the
superadiabatic Rayleigh number does not seem to change that asymmetry.  From
this point until the end of the paper, the value of $\alpha_{0}T_{0}$ is set
to $1$ ($n=0$).

###  5.1 Change of temperature profile with moderate compressibility

Figure 4: Time and horizontal averaged superadiabatic temperature profiles
along the vertical direction $z$, for $Ra_{sa}=10^{6}$ (left) and
$Ra_{sa}=10^{8}$ (right), for $\mathcal{D}$ between $0$ and $0.4$,
$\alpha_{0}T_{0}=1$, obtained in the anelastic approximation AA.

We examine now the effect of a small compressibility on the structure of
convection. From this point until the end of the paper, the value of
$\alpha_{0}T_{0}$ is set to $1$ ($n=0$). When the dissipation number
$\mathcal{D}$ is increased from $0$ to a moderate value of $0.1$ to $0.4$, a
change in the averaged vertical temperature profile is observed. In the
absence of compressible effects ($\mathcal{D}\simeq 0$), the temperature
profile is symmetrical with respect to the horizontal mid-plane as a result of
the invariance of the Boussinesq equations under the transformation
$T(x,z)\rightarrow-T(x,-z)$, $u_{x}(x,z)\rightarrow-u_{x}(x,-z)$ and
$u_{z}(x,z)\rightarrow u_{z}(x,-z)$. It is also well-known that overshoots in
the temperature profile occur near the top and bottom thermal boundary layers
[Sotin and Labrosse, 1999]. These overshoots have an amplitude (and extent)
decreasing with increasing Rayleigh numbers, but are always present. We
observe here that the overshoot near the top is nearly eliminated when the
dissipation number $\mathcal{D}$ exceeds 0.2. In Fig. 4, $\mathcal{D}$ is
increased from $0$ to $0.4$ in AA calculations and the time and horizontally
averaged superadiabatic temperature profiles are plotted along the vertical
direction for two values of the superadiabatic Rayleigh number
$Ra_{sa}=10^{6}$ and $Ra_{sa}=10^{8}$. The inset shows a close-up view of the
top overshoot in the temperature profile. A tiny value of $\mathcal{D}=0.05$
already has a noticeable effect, and when $\mathcal{D}=0.2$ most of the change
has been made. Conversely, the bottom overshoot is nearly unchanged. As a
result, the nearly constant mean value of temperature is increased, an
observation that will be related to the behaviour of the Nusselt number in
section 6.

Figure 5: Snapshots of superadiabatic temperature fields for $Ra_{sa}=10^{7}$,
$\alpha_{0}T_{0}=1$, with $\mathcal{D}=0.05$ (top) and $\mathcal{D}=0.2$
(bottom) obtained in the anelastic approximation AA.

The disappearance of the top superadiabatic temperature overshoot under
moderate compressibility can be connected to the fact that ascending plumes
fail to reach the top of the cavity. In Fig. 5, we show snapshots of the
superadiabatic temperature fields, obtained in the anelastic approximation AA,
for $Ra_{sa}=10^{7}$ and two values of the dissipation number
$\mathcal{D}=0.05$ and $\mathcal{D}=0.2$. In the first case, ascending plumes
initiated at the bottom reach the top boundary and spread horizontally,
creating the top temperature overshoot. Similarly descending plumes initiated
at the top reach the bottom and spread. When $\mathcal{D}=0.2$ however, most
ascending plumes get mixed with the surrounding fluid before they reach the
top, hence the absence of top overshoot in the temperature profile. Descending
plumes can still reach the bottom.

One should not think that compressible effects are always associated with
stronger descending plumes. This property seems to be related to the equation
of state. In another paper [Ricard et al., 2022], we consider an equation of
state suitable for planetary mantles and in that case, the opposite is true:
with compressible effects, ascending plumes are stronger.

###  5.2 Asymmetry at larger compressibility

Figure 6: Vertical profiles of conditional entropy, depending on the vertical
velocity scaled by $Ra_{sa}^{2/3}$, for the fully compressible equations FC
and anelastic models AA, ALA and SCA, $Ra_{sa}=3\,\times\,10^{5}$,
$\mathcal{D}=1.2$ and $\alpha_{0}T_{0}=1$. For FC, a value $\epsilon=0.1$ has
been taken for the superadiabatic temperature difference. The dashed profile
is the overall mean entropy profile.

We have seen in the previous section that a relatively moderate dissipation
number of $\mathcal{D}=0.2$ introduces a top-bottom asymmetry in thermal
convection. In this section, we consider a large dissipation number
$\mathcal{D}=1.2$ and investigate the asymmetry of convection depending on the
four models presented earlier: fully compressible FC, anelastic approximation
AA, anelastic liquid approximation ALA and simple compressible approximation
SCA. In Fig. 6, averaged entropy profiles are shown. They are conditional
profiles obtained for selected values of the vertical velocity in bins centred
around the indicated values, i.e. entropy profiles for parcels of a given
vertical velocity. Negative velocity values put the emphasis on descending
plumes, positive values on ascending plumes. We can see that there are more
profiles with negative velocities. This is because there are strong descending
plumes and weak ascending plumes, a tendency already seen in the previous
section for moderate dissipation numbers. The FC results are well recovered in
the anelastic models AA and ALA which show perhaps a slight increase in the
asymmetry with fewer curves corresponding to positive velocity. The entropy
scale of the FC plot is exactly ten times smaller than in the other models as
a result of the choice of the superadiabatic parameter $\epsilon=0.1$ in the
FC calculations. This difference comes from the choice of temperature scale:
it is $T_{0}$ for the FC model and $\Delta T_{sa}$ for the AA, ALA and SCA
models.

Figure 7: Probability density function (pdf) of the vertical velocity $u_{z}$
in logarithmic scale, depending on the altitude $z$, for FC, AA, ALA and SCA
models, $Ra_{sa}=3\,\times\,10^{5}$, $\mathcal{D}=1.2$ and
$\alpha_{0}T_{0}=1$. Velocity is scaled by $Ra_{sa}^{2/3}$.

A striking feature of these profiles, and in particular the overall mean
profiles (dashed lines) is that the top boundary layer has a much larger
entropy drop compared to the bottom layer in the FC, AA and ALA cases. We will
see in the next section that the superadiabatic temperature drop varies slowly
with the dissipation number, and the entropy drop is essentially equal to the
temperature drop divided by temperature. At $\mathcal{D}=1.2$, the adiabatic
temperature ratio is equal to $4$, hence an entropy drop nearly four times
larger at the top compared to that at the bottom. Meanwhile, as can be seen in
the anelastic equation (72), entropy is the driving force for convection. This
explains why descending plumes are stronger than ascending plumes. In that
respect, the situation gets somewhat similar to the case of convection cooled
from the top and thermally insulated bottom, or equivalently to the case of
volumic heating with fixed boundary temperatures [Sotin and Labrosse, 1999].

Note however that other hand-waving arguments can be put forward indicating
that ascending plumes (not descending) should be stronger in compressible
convection. For instance, the increasing value of thermal expansion as
altitude increases is thought be be such an argument leading to the
strengthening of ascending plumes and weakening of descending plumes. This
seems to apply in mantle convection [Ricard et al., 2022]. However, we have
the opposite effect in the present paper, although thermal expansion also
increases with altitude.

In terms of entropy jump, the asymmetry of the thermal boundary layers is very
much reduced in the simple anelastic model SCA (see Fig. 6) because the
adiabatic gradient is ignored in this model and the entropy drop is identical
to the temperature drop. In the SCA case, the top-bottom asymmetry still
exists, but its origin has been shifted in the thermal equation (73). In that
equation, the term $-\mathcal{D}u_{z}\tilde{T}$ is negative in ascending and
descending plumes, and consequently favours descending plumes while impeding
rising plumes. This feature of the thermal equation is really specific to the
SCA model. For the other anelastic models AA and ALA, this is not the case,
because the left-hand side term and the first term on the right-hand side of
equation (73) combine to form a conservative term
$\mathrm{D}(T_{a}\tilde{s})/\mathrm{D}t$.

In Fig. 7, we have plotted the distribution of vertical velocities for the
same simulations as those in Fig. 6. A symmetrical top-bottom configuration
would lead to contour plots symmetrical with respect to the central point
$(z=0,\ u_{z}=0)$. The asymmetry of convection under $\mathcal{D}=1.2$ is
clear, with distributions extending further toward the negative velocities
(descending plumes), while the returning ascending flow is more broadly
distributed on smaller values of positive velocity. Here again, AA and ALA
models capture very well the distribution of velocities obtained in the FC
model, although maybe with a lower probability of large values than for the FC
case. However, the SCA model displays an exaggerated asymmetry and stands
clearly away from the other models.

##  6 Malkus-type model of heat flux

The total heat flux across the fluid layer is split into the heat conducted
along the adiabat and the superadiabatic heat flux. The superadiabatic heat
flux is split itself into the conduction heat flux due to the superadiabatic
temperature difference and the convective heat flux. Finally, as we will see
in section 7, the convective heat flux is split into the convective transport
of enthalpy and the power of viscous forces.

Here, we consider the superadiabatic flux and scale it with the heat conducted
along the superadiabatic temperature difference, which is known as the Nusselt
number $Nu$. In the fully compressible model (FC), we have to subtract the
adiabatic conduction heat flux ($\mathcal{D}$ in this paper) to the total
flux, so that the Nusselt number is defined as

$Nu=\frac{Q_{FC}-\mathcal{D}}{\varepsilon}.$ (84)

One of the simplest models for the heat flux in the Boussinesq approximation
is that proposed by Malkus [1954] (see also Howard [1964]), AKA the “critical
boundary layer” model. In this model, heat flux is independent of the depth of
the fluid layer. In dimensionless terms, this leads to $Nu\sim Ra^{1/3}$. Let
us adapt this model to the case of compressible convection. Let $\delta_{c}$
and $\delta_{h}$ be the dimensionless thicknesses of the top and bottom
boundary layers, while $\Delta T_{c}$ and $\Delta T_{h}$ are the dimensionless
superadiabatic temperature jumps across these layers. The heat flux conducted
through both layers must be equal to the superadiabatic heat flux, so that

$Nu\sim\frac{\Delta T_{c}}{\delta_{c}}=\frac{\Delta T_{h}}{\delta_{h}},$ (85)

while obviously, the sum of both temperature jumps must be equal to the
superadiabatic temperature difference, in dimensionless terms:

$\Delta T_{c}+\Delta T_{h}=1.$ (86)

Concerning local Rayleigh numbers at the scale of each boundary layer, we can
build them from (49) keeping in mind that $\alpha_{0}$ should be changed
accordingly (other parameters are uniform with our equation of state). From
(52), we can relate the local value of $\alpha$ to the local adiabatic
temperature. Finally the local Rayleigh numbers can be written

$\displaystyle Ra_{c}$ $\displaystyle=Ra_{sa}\frac{\delta_{c}^{3}\Delta
T_{c}}{Ta_{c}}=Ra_{sa}\frac{\delta_{c}^{3}\Delta T_{c}}{1-\mathcal{D}/2},$
(87) $\displaystyle Ra_{h}$ $\displaystyle=Ra_{sa}\frac{\delta_{h}^{3}\Delta
T_{h}}{Ta_{h}}=Ra_{sa}\frac{\delta_{h}^{3}\Delta T_{h}}{1+\mathcal{D}/2}.$
(88)

Assuming that the local Rayleigh numbers remain equal to a number $Ra_{BL}$
independent of $\mathcal{D}$, combining (85), (86), (87) and (88) leads to the
following relationship

$Nu\sim\left(\frac{Ra_{sa}}{Ra_{BL}}\right)^{1/3}\left[\left(1-\frac{\mathcal{D}}{2}\right)^{1/4}+\left(1+\frac{\mathcal{D}}{2}\right)^{1/4}\right]^{-4/3}.$
(89)

Assuming that the boundary layer Rayleigh $Ra_{BL}$ does not depend on the
dissipation number, this expression (89) shall provide a prediction on the
effect of the dissipation number on the Nusselt number. The Rayleigh $Ra_{BL}$
defined above should not be confused with the critical Rayleigh number for the
onset of convection (based on the height of the cavity). The latter is close
to Rayleigh’s value $27\pi^{4}/4$ for a nearly uniform density [Alboussière
and Ricard, 2017].

Figure 8: Nusselt number normalized by $Ra_{sa}^{1/3}$ in the anelastic model.
The Nusselt number is defined as the ratio of the superadiabatic heat flux to
the superadiabatic pure conduction flux. The solid lines correspond to the
expression (89) with boundary layer Rayleigh numbers equal to $Ra_{BL}=14.1$,
$Ra_{BL}=17.5$ and $Ra_{BL}=19.4$, such that they are quadratic best fits of
the numerical Nusselt between $\mathcal{D}=0.4$ and $\mathcal{D}=1.8$.

Figure 8 shows the Nusselt number as a function of the dissipation number, for
three values of the superadiabatic Rayleigh number in the higher range
$Ra=10^{7}$, $Ra=10^{8}$ and $Ra=10^{9}$. From $\mathcal{D}=0$ to
approximately $\mathcal{D}=0.4$, a reduction of the Nusselt number is
observed: this should be understood as a consequence of the structure change
discussed in section 5.1 where the ascending plumes are shown to get weaker
with the disappearance of the top overshoot of the superadiabatic temperature
profile. Then for larger values of $\mathcal{D}$, the Nusselt number
increases, as predicted by our model (89) and even a little faster. The three
lines in Fig. 8 correspond to equation (89) each with a boundary Rayleigh
number $Ra_{BL}$ adjusted to fit the numerical Nusselt numbers at different
$Ra$. A general good agreement is obtained from the model of boundary layer
Rayleigh number (89).

It should be noted that a more sophisticated version of heat transfer model
exists for compressible convection. The model of Jones et al. [2022] is
designed as an extension of a model of heat transfer by Grossmann and Lohse
[2000] valid in the Boussinesq approximation. These models provide expressions
for the Nusselt number as a function of the Rayleigh number, Prandtl number
and dissipation number in the compressible case. They apply in a domain of
finite aspect ratio, i.e. bounded by vertical walls, as assumptions are made
on the large-scale velocity field. The Prandtl number can be small or large,
but not infinite. In the compressible case [Jones et al., 2022], the model
takes a different form depending on whether dissipation occurs mostly in the
boundary layers or in the bulk. Unfortunately, these models have been
developed in the no-slip boundary condition and the case of dissipation mostly
in the bulk is associated with small Prandtl numbers, and turbulence cascade.
This makes it rather difficult to apply in our case of infinite Prandtl number
and free-slip boundary conditions, even though dissipation is mainly in the
bulk in the upper range of our superadiabatic Rayleigh numbers. Moreover, the
model by Jones et al. [2022] has been developed for an EoS of a perfect gas.
This model predicts a decrease of the Nusselt number with the dissipation
number for a fixed value of the Rayleigh (and Prandtl) number, which is not
compatible with our results. That discrepancy is however not relevant, given
the mismatch between the conditions of validity of the model and our
configuration.

##  7 Heat flux and Dissipation in the different models

The expression for the vertical heat flux across horizontal planes takes a
specific form for each model of convection. In the fully compressible model
(FC), we take equation (48), integrate by parts the viscous dissipation term
and use the dot-product of equation (47) with velocity ${\bf u}$ to obtain an
expression for the heat flux as a function of height $z$

$Q_{FC}(z)=(n+1)\overline{T\rho^{n+1}u_{z}}-\frac{\varepsilon\mathcal{D}}{Ra_{sa}}\overline{u_{j}\tau_{zj}}-\frac{\mathrm{d}\overline{T}}{\mathrm{d}z},$
(90)

where the overline $\overline{\cdot}$ denotes the average over horizontal
planes, or constant-z surfaces, and over time. This is fully in accordance
with the general expression for the heat flux (see for instance (4.5) in
Curbelo et al. [2019])

$Q_{FC}(z)=\overline{\rho
hu_{z}}-\frac{\varepsilon\mathcal{D}}{Ra_{sa}}\overline{u_{j}\tau_{zj}}-\frac{\mathrm{d}\overline{T}}{\mathrm{d}z},$
(91)

as from (35) and (36), the dimensionless expression of $h$ is

$h=(n+1)T\rho^{n},\hskip 28.45274pt\mathrm{for}\ n\geq 0.$ (92)

In the statistically stationary case considered here, the function $Q_{FC}(z)$
must be independent of $z$, and we may equally denote it $Q_{FC}$.

In the anelastic approximation (AA), equations (73), (72) and (76) lead to the
following expression for the superadiabatic heat flux

$Q_{AA}(z)=\overline{\left(\tilde{T}+\frac{n}{n+1}\tilde{P}\right)u_{z}}-\frac{\mathcal{D}}{Ra_{sa}}\overline{u_{j}\tau_{zj}}-\frac{\mathrm{d}\overline{\tilde{T}}}{\mathrm{d}z}.$
(93)

This is again compatible with the general anelastic expression (see for
instance (4.6) in Curbelo et al. [2019])

$Q_{AA}(z)=\rho_{a}\overline{\tilde{h}u_{z}}-\frac{\mathcal{D}}{Ra_{sa}}\overline{u_{j}\tau_{zj}}-\frac{\mathrm{d}\overline{\tilde{T}}}{\mathrm{d}z},$
(94)

with $\rho_{a}=1$, uniform within the class of equations of state considered
in this paper, and the following linearized expression for enthalpy, from (92)
and using (75)

$\tilde{h}=\tilde{T}+\frac{n}{n+1}\tilde{P},\hskip 28.45274pt\mathrm{for}\
n>0.$ (95)

In statistically stationary cases, the time-averaged heat flux is independent
of $z$ and its value can be obtained by integrating (93) over the height of
the fluid domain, leading to

$Q_{AA}=\left<\tilde{T}u_{z}\right>+\frac{n}{n+1}\left<\tilde{P}u_{z}\right>+1,$
(96)

since the power of viscous forces can be shown to integrate to zero, because
of the continuity equation ${\bf\boldsymbol{\nabla}}\cdot(\rho_{a}{\bf
u})={\bf\boldsymbol{\nabla}}\cdot{\bf u}=0$ in the anelastic approximation, as
$\rho_{a}=1$.

In the anelastic liquid approximation (ALA), the expression for the
superadiabatic flux, obtained in the similar way as for AA, is

$Q_{ALA}(z)=\overline{\left(\tilde{T}+\tilde{P}\right)u_{z}}-\frac{\mathcal{D}}{Ra_{sa}}\overline{u_{j}\tau_{zj}}-\frac{\mathrm{d}\overline{\tilde{T}}}{\mathrm{d}z}.$
(97)

In a statistically stationary situation, the heat flux can be computed as

$Q_{ALA}=\left<\tilde{T}u_{z}\right>+\left<\tilde{P}u_{z}\right>+1.$ (98)

These expressions correspond to the limit $n\rightarrow\infty$ of those
obtained in the general anelastic approximation (AA), which also corresponds
to the limit $\alpha T\rightarrow 0$.

Finally, in the simple compressible approximation (SCA), the expressions for
the heat flux are identical to those in the anelastic liquid approximation
(ALA).

However, the heat flux can be written differently than in the expressions
above. For instance, in the anelastic approximation (AA), instead of (93),
(94) or (96), one can base it on the flux of entropy. From Gibbs equation
$\mathrm{d}h=T\mathrm{d}s+\mathrm{d}P/\rho$, expression (94) can be written

$Q_{AA}(z)=T_{a}\overline{\tilde{s}u_{z}}+\overline{\tilde{P}u_{z}}-\frac{\mathcal{D}}{Ra_{sa}}\overline{u_{j}\tau_{zj}}-\frac{\mathrm{d}\overline{\tilde{T}}}{\mathrm{d}z}.$
(99)

In the same time, from the horizontal and time average of the dot-product of
Navier-Stokes with velocity, one obtains the expression for the viscous
dissipation at each vertical position

$\frac{\mathcal{D}}{Ra_{sa}}\overline{\dot{\epsilon}:\tau}(z)=\mathcal{D}\overline{\tilde{s}u_{z}}-\frac{\mathrm{d}}{\mathrm{d}z}\left[\overline{\tilde{P}u_{z}}-\frac{\mathcal{D}}{Ra_{sa}}\overline{u_{j}\tau_{zj}}\right].$
(100)

Introducing

$G(z)=\overline{\tilde{P}u_{z}}-\frac{\mathcal{D}}{Ra_{sa}}\overline{u_{j}\tau_{zj}},$
(101)

equations (99) and (100) can be rewritten

$\displaystyle Q_{AA}(z)$
$\displaystyle=T_{a}\overline{\tilde{s}u_{z}}+G(z)-\frac{\mathrm{d}\overline{\tilde{T}}}{\mathrm{d}z},$
(102)
$\displaystyle\frac{\mathcal{D}}{Ra_{sa}}\overline{\dot{\epsilon}:\tau}(z)$
$\displaystyle=\mathcal{D}\overline{\tilde{s}u_{z}}-\frac{\mathrm{d}G(z)}{\mathrm{d}z}.$
(103)

It can be noted that these expressions are valid for any general equation of
state and can also be generalized when inertia and possible Lorentz forces are
included: it suffices to add the flux of inertia to the expression of $G(z)$,
as detailed in appendix C. In the anelastic liquid approximation, the
expression for $G(z)$ is unchanged and the flux and dissipation profiles
become

$\displaystyle Q_{ALA}(z)$
$\displaystyle=\overline{\tilde{T}u_{z}}+G(z)-\frac{\mathrm{d}\overline{\tilde{T}}}{\mathrm{d}z},$
(104)
$\displaystyle\frac{\mathcal{D}}{Ra_{sa}}\overline{\dot{\epsilon}:\tau}(z)$
$\displaystyle=\mathcal{D}\frac{\overline{\tilde{T}u_{z}}}{T_{a}}-\frac{\mathrm{d}G(z)}{\mathrm{d}z}.$
(105)

In the simplest SCA approximation, $G(z)$ is still the same and the flux and
dissipation profiles are

$\displaystyle Q_{SCA}(z)$
$\displaystyle=\overline{\tilde{T}u_{z}}+G(z)-\frac{\mathrm{d}\overline{\tilde{T}}}{\mathrm{d}z},$
(106)
$\displaystyle\frac{\mathcal{D}}{Ra_{sa}}\overline{\dot{\epsilon}:\tau}(z)$
$\displaystyle=\mathcal{D}\overline{\tilde{T}u_{z}}-\frac{\mathrm{d}G(z)}{\mathrm{d}z}.$
(107)

###  7.1 Global dissipation

The total Dissipation is expressed as a fraction of the heat flux, more
specifically the superadiabatic convective heat flux. The reference result is
obtained in the Boussinesq approximation, where it has been shown that the
integrated viscous dissipation is equal to the product of the dissipation
number and the convective heat flux [Howard, 1963, Hewitt et al., 1975]. We
express the total viscous dissipation as a fraction of the conduction heat
flux along the superadiabatic gradient and denote it $D_{\nu}$. In the fully
compressible model (FC), the integrated viscous dissipation divided by $\Delta
T_{sa}$ has the following expression

$D_{\nu}=\frac{\mathcal{D}}{Ra_{sa}\,\varepsilon}\left<\dot{\epsilon}:\tau\right>.$
(108)

In the approximated models (AA), (ALA) and (SCA), the expression is

$D_{\nu}=\frac{\mathcal{D}}{Ra_{sa}}\left<\dot{\epsilon}:\tau\right>,$ (109)

because the dimensionless temperature is scaled using $\Delta T_{sa}$ for
these models. Considering the expressions for the Nusselt number in section 6,
the ratio $E$ of viscous dissipation to the convective heat flux takes the
same expression for all models

$E=D_{\nu}/(Nu-1).$ (110)

The classical Boussinesq result on dissipation comes from integrating the dot
product of Navier-Stokes with the velocity field. It follows that viscous
dissipation is exactly equal to the convective heat flux multiplied by the
dissipation number [Howard, 1963], hence ${\color[rgb]{0,0,0}E=}\
D_{\nu}/(Nu-1)=\mathcal{D}$ in the Boussinesq limit. For this reason, we call
the quantity $E-\mathcal{D}$ the ratio of dissipation to convective flux in
excess of $\mathcal{D}$. It is zero in the Boussinesq case and we compute it
for compressible convection FC, AA, ALA and SCA.

The numerical results concerning viscous dissipation are shown in Fig. 9. For
small values of $\mathcal{D}$, all results are close to $0$, indicating that
the Boussinesq results apply: dissipation is equal to the product of the
convective heat flux and the dissipation number. When the dissipation number
is increased, the simple compressible approximation model (SCA) goes to
slightly negative values, whereas all other models go to positive values. So,
for all models except (SCA), viscous dissipation becomes larger than predicted
by the Boussinesq approximation. That departure increases with the dissipation
number $\mathcal{D}$ and also with the superadiabatic Rayleigh number
$Ra_{sa}$ (better seen in Fig. 10).

Figure 9: Excess of the ratio of dissipation to convective heat flux relative
to $\mathcal{D}$, as a function of $\mathcal{D}$, for FC, AA, ALA, SCA models,
$\alpha_{0}T_{0}=1$, $\mathcal{D}$ in
$[0.25,\,0.5,\,0.75,\,1.0,\,1.25,\,1.5]$, $Ra_{sa}$ in $[10^{3},\
10^{3.5},\,10^{4},\ 10^{4.5},\,10^{5},\ 10^{5.5},\,10^{6},\ 10^{6.5}]$. Figure
10: Excess of the ratio of dissipation to convective heat flux relative to
$\mathcal{D}$, as a function of $Ra_{sa}$ (same values as in Fig. 9).

Let us consider some limits to the dissipation results. First, we have
rigorous upper (and lower) bounds, since the papers of Backus [1975] and
Hewitt et al. [1975]. Taking equation (73), dividing by $T_{a}$ and
integrating by parts leads to

$\frac{\mathrm{D}\tilde{s}}{\mathrm{D}t}-{\bf\boldsymbol{\nabla}}\cdot\left(\frac{{\bf\boldsymbol{\nabla}}\tilde{T}}{T_{a}}\right)=\frac{\mathcal{D}}{Ra_{sa}}\frac{\dot{\epsilon}:\tau}{T_{a}}+\frac{{\bf\boldsymbol{\nabla}}\tilde{T}\cdot{\bf\boldsymbol{\nabla}}T_{a}}{T_{a}^{2}},$
(111)

where we recognize the positive sources of entropy due to irreversible viscous
dissipation and conduction on the right-hand side. The upper bound for total
dissipation occurs for negligible conduction (as an entropy source) and in the
case when dissipation takes place at the largest possible temperature
$T_{a}=1+\mathcal{D}/2$, at the bottom of the domain. Integrating (111) over
the whole domain leads to

$\frac{D_{\nu}}{1+\frac{\mathcal{D}}{2}}<Nu\left[\frac{1}{1-\frac{\mathcal{D}}{2}}-\frac{1}{1+\frac{\mathcal{D}}{2}}\right],$
(112)

implying the upper bound

$\frac{D_{\nu}}{Nu-1}<\mathcal{D}+\frac{\mathcal{D}^{2}}{2-\mathcal{D}},$
(113)

in the limit of a large Nusselt number $Nu\gg 1$.

Another possible limit case will be shown to correspond to our numerical
results at large Rayleigh number in the next section 7.2. That limit is that
of a vanishing contribution of the $G(z)$ component of the heat flux defined
in equation (101) in the anelastic approximation. It corresponds to the case
when the heat flux is carried entirely by the flux of entropy
$T_{a}\overline{\tilde{s}u_{z}}$ (outside top and bottom conduction layers)
while the energy dissipation at each height is
$\mathcal{D}\overline{\tilde{s}u_{z}}$, as can be seen from equations (102)
and (103). Under that assumption, the integral of energy dissipation is equal
to

$D_{\nu}=Nu\int_{1/2}^{1/2}\frac{\mathcal{D}}{T_{a}}\mathrm{d}z,$ (114)

leading to the expression

$\frac{D_{\nu}}{Nu-1}=\log\left(\frac{1+\mathcal{D}/2}{1-\mathcal{D}/2}\right),$
(115)

The expressions for the upper bound (113) and for the “entropy flux” model
(115) are plotted in Fig. 9. It can be seen that the upper bound curve (113)
lies far above the numerical results and that the ’entropy flux’ expression
(115) seems to correspond to the limit of the numerical results when the
Rayleigh number is increased, for FC, AA and ALA models. For small
(supercritical) Rayleigh numbers, dissipation is close to the ’Boussinesq’
limit. For the SCA model, the behaviour is different: starting from negative
values at small $Ra_{sa}$ numbers, the ’Boussinesq’ limit is reached for large
$Ra_{sa}$ numbers. However, from a fundamental perspective, the SCA model does
not behave differently from the other models if one considers the consequences
of the ’entropy flux’ assumption, because of the different expressions for the
flux and dissipation. From (106) and (107), it is clear that the assumption of
a negligible $G$ contribution leads to the usual ’Boussinesq’ limit for the
SCA model, just because it is identical to the ’entropy flux’ limit in that
simple compressible approximation. Thus, in the limit of large $Ra_{sa}$, when
the contribution of $G$ is small, it is expected that dissipation becomes
close to the product of the dissipation number and the convective heat flux.
This is what we observe on Figs. 9 and 10.

It is striking that in Fig. 10 at low Rayleigh number ALA, FC and AA results
differ significantly. Dissipation is smallest with AA, then FC and highest
with ALA. Typically AA excess dissipation goes to zero near the critical
Rayleigh number, while ALA excess dissipation shows an increase to some finite
value. FC results are intermediate. When the superadiabatic Rayleigh number
becomes larger, the difference between AA, FC and ALA tends to shrink until
the predicted dissipation is the same above $10^{5}$ or $10^{6}$.

###  7.2 Convergence of heat flux and dissipation profiles

Figure 11: Entropy flux contribution for $\mathcal{D}$ equal to $0.8$ (left)
and $1.6$ (right), for a superadiabatic Rayleigh number up to $10^{9}$ in the
anelastic approximation and $\alpha_{0}T_{0}=1$.

Figure 12: Dissipation profile for $\mathcal{D}$ equal to $0.8$ (left) and
$1.6$ (right), $\alpha_{0}T_{0}=1$, for a superadiabatic Rayleigh number up to
$10^{9}$ in the anelastic approximation.

Figure 13: Relative distance ($L^{1}$-norm) of the entropy to total heat flux
profile (left) and of the dissipation profile to the limit (146).

In this section, we focus on the convergence of the entropy heat flux and of
the dissipation profiles towards universal limit profiles, when the
superadiabatic Rayleigh number becomes large. We have restricted our analysis
to a maximum value of $Ra_{sa}=10^{9}$, so that a spatial resolution of $512$
vertical and $2048$ horizontal modes was able to capture thin plumes and
boundary layers. For simplification and ease of calculation, we consider the
anelastic approximation AA only: from the previous section 7, we have seen
that the global amount of dissipation does not seem to depend on the model FC,
AA or ALA, at large $Ra_{sa}$ numbers, so that we expect the same convergence
for FC and ALA than that from AA (SCA is different, as we have seen).

Figure 11 shows the ratio of the entropy heat flux
$T_{a}\overline{\tilde{s}u_{z}}$ to the uniform total heat flux $Q_{AA}$, see
equation (102). For both values of the dissipation number $\mathcal{D}=0.8$
and $\mathcal{D}=1.6$, the entropy heat flux profile converges (slowly)
towards the uniform value $1$, except in thin boundary layers: their thickness
is of order $Ra_{sa}^{-1/3}$ and conduction is the only way to transfer heat
in the vicinity of the top and bottom boundaries. This implies that the other
flux contributions, called together $G(z)$, see (101) and (102), become more
and more negligible as the superadiabatic Rayleigh number is increased.

Figure 14: Snapshot of the superadiabatic temperature field for a very small
$\mathcal{D}=0.05$ (top), moderate $\mathcal{D}=0.4$ (middle) to a large
dissipation number $\mathcal{D}=1.6$ (bottom), for $Ra_{sa}=10^{9}$, in the
anelastic approximation with $\alpha_{0}T_{0}=1$. In the close-up views, the
distribution of viscous dissipation $\dot{\epsilon}:\tau$ (left) and entropy
flux $\tilde{s}u_{z}$ (right) are shown.

As expected, as $G(z)$ becomes small (and without variations at some small
scale), so does $\mathrm{d}G(z)/\mathrm{d}z$, implying that the profile of
viscous dissipation becomes close to $\mathcal{D}\overline{\tilde{s}u_{z}}$,
see (103). As $T_{a}\overline{\tilde{s}u_{z}}$ is close to $1$ at large
$Ra_{sa}$, we have a dissipation profile close to $1/T_{a}(z)$ (see appendix C
for a general dimensional expression). Figure 12 shows the change in
dissipation profiles as the superadiabatic Rayleigh number increases from
$10^{4}$ to $10^{9}$, for two different values of the dissipation number,
$\mathcal{D}=0.8$ and $\mathcal{D}=1.6$. As expected from (103), the profile
converges towards $1/T_{a}(z)$ in both cases, however the ”entropy” component
of dissipation being multiplied by $\mathcal{D}$, the convergence looks more
obvious for the larger value of $\mathcal{D}$.

We give a quantitative measure of the convergence of the entropy flux
contribution (towards $1$) and dissipation profiles (towards $1/T_{a}(z)$) in
Fig. 13. The measure is defined in each case as the logarithm of the
$\mathrm{L}^{1}$ distance. The larger the dissipation number, the faster the
convergence is.

In an attempt to understand how the heat flux contributions $G(z)$ become
negligible as $Ra_{sa}$ is increased, we plot a snapshot of the superadiabatic
temperature field for three different values of the dissipation number
$\mathcal{D}=0.05,\ 0.4\ \mathrm{and}\ 1.6$ at $Ra_{sa}=10^{9}$ in Fig. 14. At
small $\mathcal{D}$ sparse plumes go from bottom to top or top to bottom and a
velocity field is generated with a length-scale comparable to the height of
the cavity, in the case considered here of infinite Prandtl number. At
moderate $\mathcal{D}$ the top-bottom asymmetry is strong, more plumes are
present and smaller scales are visible. At the largest $\mathcal{D}$, numerous
plumes exist and they cannot cross the whole height of the cell (not even the
descending ones) without their heads detaching from their tails and continue
their course as isolated blobs. The length-scale of typical distance between
adjacent plumes $l$ is reduced significantly compared to the height $H$ of the
cavity. If $U$ is a typical scale for velocity, then we expect the local
viscous dissipation $\overline{\dot{\epsilon}:\tau}$ to scale as
$U^{2}/l^{2}$, while $u_{j}\tau_{zj}$ scales as $U^{2}/l$. Once averaged in
time and horizontally, that quantity depends smoothly on $z$ on the global
scale $H$, so that $\mathrm{d}/\mathrm{d}z(\overline{u_{j}\tau_{zj}})$ scales
as $U^{2}/(lH)$, i.e. $l/H$ smaller than viscous dissipation. Extending this
result on deviatoric stress work to pressure work, this explains that
$\mathrm{d}G(z)/\mathrm{d}z$ becomes negligible compared to dissipation, see
(103).

In Fig. 15, we plot the energy spectrum of one component of the deformation
rate tensor, $\partial u_{z}/\partial z$, for a large Rayleigh number
$Ra=10^{9}$ and different dissipation numbers, in the anelastic approximation
(AA). As the dissipation number is increased, a shift of the whole spectrum
towards larger wavenumbers is observed. At small values of $\mathcal{D}$, the
spectrum has a maximum around $k_{x}=4$, while at large $\mathcal{D}$ the
maximum is close to $k_{x}=10$. This indicates that for a given integral of
viscous dissipation, the fraction of heat flux carried by the work done by
viscous stresses becomes smaller and smaller as the dissipation number
$\mathcal{D}$ is increased.

Figure 15: Energy spectrum of $\partial u_{z}/\partial z$: the Fourier
coefficients of this component of the velocity gradient are computed along the
$x$ direction and the square of their complex magnitude $\mathcal{F}(k_{x})$
is plotted as a function of the wavenumber ($k_{x}=1$ corresponds to a signal
of period $L$ the length of the domain along $x$, and for any value of $k_{x}$
the corresponding mode is $e^{ik_{x}2\pi x/L}$). Those Fourier coefficients
are averaged over a central vertical extent from $z=-0.11$ to $z=0.11$.

This idea of smaller scales for velocity gradients at larger superadiabatic
Rayleigh numbers can also give a hint on why AA, FC and ALA results converge
at large $Ra_{sa}$, as seen in Fig. 10. An estimate of pressure contributions
to entropy fluctuations in Anufriev et al. [2005] – obtained from the analysis
of the order of magnitude of the forces in the momentum equation – is adapted
to a lengthscale of convection $l$: from Stokes’ equation, an order of
magnitude of pressure is $\alpha\rho g(\delta T)l$ and from
$\mathrm{d}s=c_{p}/T\mathrm{d}T-\alpha/\rho\mathrm{d}P$, we evaluate the ratio
of pressure over temperature to be of order $\alpha T\mathcal{D}l/H$. In
Anufriev et al. [2005], the length-scale $l$ is taken to be of order $H$ and
the condition of validity for the ALA approximation is given as $\alpha
T\mathcal{D}\ll 1$. If however, a smaller lengthscale $l$ prevails, that
estimate must be smaller proportionally to $l/H$ making the ALA approximation
more valid.

At large dissipation number and large superadiabatic Rayleigh number –
typically our last case of Fig. 14 – we thus propose that the time and
horizontal average of viscous dissipation is linked to the time and horizontal
average of the product $\tilde{s}u_{z}$ (see equation (103) with negligible
term $\mathrm{d}G/\mathrm{d}z$). However, both quantities
$\dot{\epsilon}:\tau$ and $\tilde{s}u_{z}$ are not pointwise (and timewise)
correlated. This is illustrated in the close-up views in Fig. 14 where
$\dot{\epsilon}:\tau$ and $\tilde{s}u_{z}$ are plotted in a small region of
the fluid domain. The quantity $\tilde{s}u_{z}$ (right close-up) represents
well the plumes while dissipation $\dot{\epsilon}:\tau$ (left close-up) looks
more like a halo around the plumes. This is also a consequence of the
different symmetries concerning each quantity: for a supposedly straight
descending plume, $\tilde{s}u_{z}$ is maximum on the centerline of the plume,
however dissipation must be zero on that centerline (where vertical gradients
are smaller than horizontal gradients, i.e. not ahead of the tip of the
plume).

##  8 Effect of confinement and inertia

Different results were obtained in Currie and Browning [2017], where a larger
dissipation than our limit (114) is obtained. There are actually several
differences in the configuration they have studied: i) they use an EoS of an
ideal gas, ii) they model conduction using the gradient of entropy, iii) they
use a boundary condition of a fixed flux (bottom), iv) they consider inertia
($Pr=1\ \mathrm{or}\ 10$), v) they have a square domain, vi) they have non-
penetrative conditions on lateral walls. In this section, we test changes in
our configuration regarding the last three points iv), v) and vi) that can be
implemented easily in our code. It turns out that we need to make all three
changes to recover the results of Currie and Browning [2017].

In the results presented here, we have kept the same equation of state as in
the beginning of the paper, with $s=s(\rho)$. We have kept the same top and
bottom boundary conditions. We have included inertia with a Prandtl number
equal to ten, $Pr=10$. We have changed the aspect ratio from $4\sqrt{2}$ to
$1$ (square domain). Then, we consider two cases, one with periodic lateral
boundary conditions as before in this paper, one with a non-penetrative
boundary condition. That last case is that considered in Currie and Browning
[2017] and corresponds to a vanishing perpendicular velocity component
$u_{x}=0$ and no shear-stress $\partial u_{z}/\partial x=0$. This is achieved
numerically in dedalus by choosing the so-called SinCos base of functions for
horizontal decomposition, instead of the complete Fourier base for periodic
boundary conditions.

Figure 16, on the left-hand side, shows the averaged vertical profiles of
different components of the heat flux identified in (145): the entropy flux
fraction of the heat flux $\rho_{a}T_{a}\overline{u_{z}\tilde{s}}$, the
kinetic energy flux $1/2\rho_{a}\overline{u^{2}u_{z}}$, the viscous work
$-\overline{u_{i}\tau_{iz}}$ and the pressure work
$\overline{\tilde{P}u_{z}}$. On the right-hand side of Fig. 16, the averaged
profiles of viscous dissipation are plotted for $Pr=10$, $Ra_{sa}=10^{7}$ and
$\mathcal{D}=1.6$. With periodic boundary conditions, we observe some
departure from the results we had previously (without inertia and in a long
cavity) in the top half of the cavity, but this does not change the total
dissipation significantly. On the contrary, with the confined boundary
condition on lateral boundaries (SinCos in dedalus), the fraction of entropy
heat flux exceeds $1$ by 75% in the bottom half, and as a consequence of (102)
and (103), viscous dissipation reaches much larger values there. This brings
the total dissipation in the range of that obtained by Currie and Browning
[2017]. Looking into more details, the flux of kinetic energy is very
significantly negative in the SinCos case, causing a significant increase of
the entropy flux in the lower half of the fluid domain. On Fig. 17, we compare
a snapshot of the superadiabatic temperature field in both cases (periodic or
SinCos). We can see that the vertical ”walls” in the SinCos case – and inertia
– are capable of channeling the descending plume that dissipates its kinetic
energy at the bottom. A large-scale circulation is created. Note that Tilgner
[2011] has led numerical simulations in the fully compressible case with
explicit vertical walls. In the periodic case however, descending plumes still
cannot reach the bottom (even with $Pr=10$ instead of $Pr=\infty$) and finer
scales develop. Another feature of the periodic box is the large scale shear
deformation, induced by Reynolds stresses (with a finite Prandtl number),
which is suppressed in the wall-bounded geometry.

Figure 16: Profiles of heat flux components from (145) except conduction terms
(left) and dissipation profile (right), for an aspect ratio unity and a
Prandtl number equal to 10, $Ra_{sa}=10^{7}$, $\mathcal{D}=1.6$. Comparison
between a periodic boundary conditions (Square periodic, in green) and
horizontal confinement with no shear stress (Square SinCos, in red).

Figure 17: Snapshot of the superadiabatic temperature field for $Pr=10$,
$Ra_{sa}=10^{7}$, $\mathcal{D}=1.6$, $\alpha_{0}T_{0}=1$, in a domain of
aspect ratio unity. On the left-hand side, there are periodic boundary
conditions in the horizontal direction, while there are impermeable no-shear-
stress walls for the calculation on the right-hand side.

So, it seems that there are two possible convection states. One that is
dominated by the entropy flux, where $G(z)$ in the flux profile (102) and in
the dissipation profile (103) is negligible. The length scales of convection
are small. In our work, this state is obtained in the limit of large
superadiabatic Rayleigh numbers, although that limit might be difficult to
reach for small dissipation numbers. A second state, observed by Currie and
Browning [2017], can be seen when inertia is introduced and when vertical
walls are present. Convection is large-scale, descending plumes are stuck to a
wall (the left wall in the particular snapshot in Fig. 17) and viscous
dissipation is enhanced compared to the first state. Our investigations
concerning the effect of inertia and confinement have not been systematic and
complete and they only concern the anelastic approximation AA. We have tested
all combinations of two aspect ratios (our usual $4\sqrt{2}$ and $1$), two
types of boundary conditions in the horizontal direction (periodic or SinCos)
and two values of the Prandtl number ($Pr=\infty$ and $Pr=10$). Yet, that
second state was only seen when we had simultaneously vertical walls (SinCos
boundary condition), inertia ($Pr=10$) and an aspect ratio equal to one.
Incidentally, this configuration is that of most experiments where the
‘ultimate’ regime of thermal convection is investigated, for instance in Roche
et al. [2002].

In their paper, Currie and Browning [2017] have tested the model (that we call
’entropy flux’) leading to the global dissipation (115), but they discard it
on the basis of their numerical results. However, they write ”Often it is
assumed that in the bulk of the convection zone, the total heat flux is just
equal to the convective flux…”, where they call ”convective flux” that part of
flux we call here ”entropy flux”, $T_{a}\overline{\tilde{s}u_{z}}$. This
corresponds precisely to the assumption of negligible contribution of $G(z)$
to the heat flux in (102). So the idea has been expressed already and is
seemingly well accepted in astrophysics, but we could not find a precise
reference for it until now.

At this point, it is important to have in mind the specific nature of the
present study. The equation of state considered here is a peculiar one, a sort
of limit case retaining the thermal features of compressibility (adiabatic
gradient, adiabatic cooling) but minimizing the actual changes in density
(nearly uniform density). Previous works on stellar or gas planet convection
with an ideal gas equation of state have mostly shown that the flux of kinetic
energy corresponds to a significant fraction of the total heat flow [Chan and
Sofia, 1989, Viallet et al., 2013, Featherstone and Hindman, 2016, Käpylä et
al., 2019]. These studies report the existence of deep convective plumes
crossing the whole adiabatic layer. The impact of density stratification in an
adiabatic region is highlighted by Anders et al. [2019] who point out its
effect on descending plumes that can be sometimes compacted and accelerated
downward. However, their authors also insist on the gap between the values of
the stellar Rayleigh number and those actually accessible numerically. In
their three-dimensional calculations, the Rayleigh number based on the heat
flux is restricted to be less than $Ra_{F}=10^{7.5}$ (under a classical
scaling $Nu\sim Ra^{1/3}$, this corresponds to a Rayleigh number based on a
temperature difference of $Ra=Ra_{F}^{3/4}=10^{5.625}$). In these papers, it
is also reported that the numerical models overestimate stellar convective
velocities compared to the observations: this might be related to the relative
small values of the numerically accessible Rayleigh numbers. For instance,
Featherstone and Hindman [2016] report a slow tendency of the spectrum of
convection to shift to larger wavenumbers as the Rayleigh number is increased.
This is also something we observe (see Fig. 15) and we associate this to a
slow convergence toward a regime of heat flow dominated by the entropy flux.
In the context of the Boussinesq approximation, Goluskin et al. [2014] show
how shear flow is generated by convection in a periodic domain, leading to a
reduction of the vertical heat transfer. In their two-dimensional case, they
can reach Rayleigh numbers of $10^{10}$. In addition to the role of the
equations of state, other features of the model can potentially affect the
final structure of the flow: imposed temperatures versus imposed heat flux,
sub-grid-scale models in particular under the form of a Fourier heat flux
proportional to the gradient of entropy… Finally, the very important effects
of rotation and magnetic field are not considered here: for instance, the
dynamics of the Earth’s core is dominated by the influence of the Coriolis
force and one consequence is that kinetic energy is negligible despite the
small value of the Prandtl number [Schaeffer et al., 2017].

##  9 Conclusions

In this paper, we have taken the limit case of a class of equations of state
such that entropy is a function of density. In the assumption of an infinite
Prandtl number, we have written the fully compressible governing equations of
convection as well as anelastic approximations of increasing simplification
AA, ALA and SCA (anelastic, anelastic liquid and simple compressible). Under
that choice, a nearly uniform entropy field implies that density is nearly
uniform. A consequence is that, with a uniform (dynamic) viscosity, we also
have a nearly uniform kinematic viscosity and a uniform thermal diffusivity.
With such a class of equations of state, there is still an adiabatic
temperature gradient and its effect in heat transport is still present. The
idea is to keep features of compressible convection and to discard effects
related to non-uniform fluid properties, sometimes called non-Oberbeck-
Boussinesq (NOB) effects. The anelastic approximation AA is based on a linear
expansion about a state of uniform entropy, in the anelastic liquid
approximation ALA pressure variations are neglected on all thermodynamic
quantities, and in our simple compressible approximation SCA even the
adiabatic gradient of temperature is eliminated. In that last approximation,
thermodynamics is badly treated, however adiabatic heating and cooling is
retained and the mathematical structure of the equations contains the key
ingredients of compressible convection. It has the advantage of simplicity
with just two scalar parameters $Ra_{sa}$ and $\mathcal{D}$. There is no need
to determine a profile of adiabatic temperature. Applied mathematicians might
want to play with that system (see also appendix A.4) and determine
fundamental properties of its solutions, which might then have applications in
the more physical models.

The genuine compressible effects are governed by the dissipation parameter
$\mathcal{D}$. Around $\mathcal{D}=0.1$ (between 0.05 and 0.2), compressible
effects create a top-bottom asymmetry. Ascending plumes starting from the
bottom thermal boundary layer do no longer reach the top. This leads to a
change in the temperature profile with the loss of the overshoot near the top
of the superadiabatic temperature profile. At larger values of the dissipation
number, the change in the global heat flux under a constant superadiabatic
Rayleigh number is compatible with the model of critical boundary layer of
Malkus [1954]. As we have shown in section 6, increasing the dissipation
number leads to an increase of heat transfer owing to the asymmetry of heat
transfer resistance of the top and bottom thermal boundary layers.

With significant compressible effect ($\mathcal{D}$ large enough) and in the
limit of very large superadiabatic Rayleigh numbers, we have shown that a
state of ”local equilibrium” is reached, where the heat flow due to the flux
of entropy fluctuations is accompanied with the corresponding viscous
dissipation at the same height. This is a small-scale process which is
consistent with the observation that heat flux components such as shear-stress
or pressure fluxes are negligible. In that limit, we can predict the vertical
profile of viscous dissipation, as soon as the profile of $\alpha$, $g$ and
$c_{p}$ are known from equation (146). A similar process takes place in the
simplest SCA model: however, due to its extreme simplicity, it ignores the
depth dependence of the expansion coefficient and a constant heat flux with
depth leads to a constant profile of viscous dissipation (see equations (106)
and (107)). Although the result is different – and certainly less relevant to
geophysics and astrophysics – it is still mathematically interesting to
investigate the consequences of the SCA model on dissipation distribution. It
still contains viscous heating and adiabatic cooling in the thermal equation.
These terms – particularly at high superadiabatic Rayleigh numbers – are
capable to drive the system in a state of mesoscale equilibrium between
themselves and lead to a dominant mode of heat transfer due to the flux of
entropy, while other modes (kinetic energy flux, pressure and stress terms)
become small. If this process is understood in the simpler SCA model, then it
could certainly help understand the behaviour of the other models.

Other works, in particular in Currie and Browning [2017], find that a
different type of flows can exist, with large-scale circulation and more
dissipation for the same heat flux. In that case, the vertical profile of
dissipation has larger values in the lower half of the domain where
temperature is large (and less in the upper half): this explains how the
entropy budget is balanced despite an increased global dissipation. We have
observed that this type of convection can be reached only when the Prandtl
number is finite, vertical walls are present and the aspect ratio of the
domain is not large. In geophysical or astrophysical contexts, vertical walls
certainly do not exist and that second type of flow appears unlikely to
develop.

Regarding the different models, we find a good agreement between FC and AA, as
expected, except at large dissipation numbers and small superadiabatic
Rayleigh numbers (above threshold). We also find the ALA approximation to be
good, especially for small values of $\alpha_{0}T_{0}$ (as expected again),
and at large superadiabatic Rayleigh numbers: that last feature may be due to
the fact that we get solutions with small convective scales (mesoscale
equilibrium) for which pressure variations are small. The differences would be
larger in the case of large-scale circulations. The SCA model gives different
results, but there are good reasons for that. As explained above, this model
is not meant to provide realistic results. It should be seen as a simple set
of equations where compressible effects can still be studied.

In future works, it seems important to investigate precisely the conditions of
existence of both types of flow. The following features should be studied:
realistic equations of state, three-dimensional flows, electromagnetic forces,
effect of rotation. If the flux of kinetic energy or the Poynting flux cannot
compete with the entropy flux, then it is likely that the model of “local
equilibrium” applies, embodied by the vertical distribution of dissipation
(146).

##  Appendix A Equations for FC, AA, ALA and SCA models

The velocity boundary conditions are common to all models: stress-free, non-
penetrative conditions, with an additional constraint on the average
horizontal velocity since Galilean invariance does not constrain it.

$\displaystyle\frac{\partial u_{x}}{\partial z}=0,\hskip
28.45274pt\mathrm{when}\hskip 28.45274ptz=\pm\frac{1}{2},$ (116)
$\displaystyle u_{z}=0,\hskip 28.45274pt\mathrm{when}\hskip
28.45274ptz=\pm\frac{1}{2}.$ (117)
$\displaystyle\int_{-L/(2H)}^{L/(2H)}u_{x}\left(x,z=\frac{1}{2}\right)\mathrm{d}x=0,$
(118)

The non-dimensional equations of the different models studied – fully
compressible (FC), anelastic approximation (AA), anelastic liquid
approximation (ALA) and simple compressible approximation model (SCA) – and
other boundary conditions are as follows:

###  A.1 Fully compressible model FC

The governing equations are

$\displaystyle\frac{\mathrm{D}\rho}{\mathrm{D}t}$
$\displaystyle=-\rho{\boldsymbol{\nabla}}\cdot{\bf u},$ (119)
$\displaystyle{\mathbf{0}}$
$\displaystyle=-\frac{Ra_{sa}(n+1)}{\varepsilon}\left(\frac{\rho^{n}}{\mathcal{D}}\left[(n+1)T{\mathbf{\boldsymbol{\nabla}}}\rho+\rho{\mathbf{\boldsymbol{\nabla}}}T\right]+\rho{\mathbf{e}}_{z}\right)+{\mathbf{\boldsymbol{\nabla}}}^{2}{\bf
u}+\frac{1}{3}{\mathbf{\boldsymbol{\nabla}}}\left(\boldsymbol{\nabla}\cdot{\bf
u}\right),$ (120) $\displaystyle 0$
$\displaystyle=-(n+1)\rho^{n+1}T\boldsymbol{\nabla}\cdot{\bf
u}+\frac{\varepsilon\mathcal{D}}{Ra_{sa}}\dot{\epsilon}:\tau+\nabla^{2}T.$
(121)

Temperatures are imposed at the top and bottom

$\displaystyle T$ $\displaystyle=1+\frac{\mathcal{D}+\varepsilon}{2},\hskip
28.45274pt\mathrm{when}\hskip 28.45274ptz=-\frac{1}{2},$ (122) $\displaystyle
T$ $\displaystyle=1-\frac{\mathcal{D}+\varepsilon}{2},\hskip
28.45274pt\mathrm{when}\hskip 28.45274ptz=\frac{1}{2}.$ (123)

###  A.2 Anelastic Approximation AA

The governing equations are

$\displaystyle{\mathbf{\boldsymbol{\nabla}}}\cdot{\bf u}$ $\displaystyle=0,$
(124) $\displaystyle{\bf 0}$
$\displaystyle=-\frac{Ra_{sa}}{\mathcal{D}}{\mathbf{\boldsymbol{\nabla}}}\tilde{P}+Ra_{sa}\left(\frac{\tilde{T}}{T_{a}}-\frac{\tilde{P}}{(n+1)T_{a}}\right){\hat{\mathbf{e}}}_{z}+{\mathbf{\boldsymbol{\nabla}}}^{2}{\bf
u},$ (125)
$\displaystyle\frac{\mathrm{D}}{\mathrm{D}t}\left(\tilde{T}-\frac{\tilde{P}}{n+1}\right)$
$\displaystyle=-\mathcal{D}u_{z}\left(\frac{\tilde{T}}{T_{a}}-\frac{\tilde{P}}{(n+1)T_{a}}\right)+\frac{\mathcal{D}}{Ra_{sa}}\dot{\epsilon}:\tau+\nabla^{2}\tilde{T},$
(126)

and boundary conditions

$\displaystyle\tilde{T}\left(z=\pm\frac{1}{2}\right)=\mp\frac{1}{2}.$ (127)
$\displaystyle\int_{0}^{L/H}\tilde{P}\left(x,z=\frac{1}{2}\right)-\tilde{P}\left(x,z=-\frac{1}{2}\right)\mathrm{d}x=0,$
(128)

###  A.3 Anelastic Liquid Approximation ALA

The governing equations are

$\displaystyle{\mathbf{\boldsymbol{\nabla}}}\cdot{\bf u}$ $\displaystyle=0,$
(129) $\displaystyle{\bf 0}$
$\displaystyle=-\frac{Ra_{sa}}{\mathcal{D}}{\mathbf{\boldsymbol{\nabla}}}\tilde{P}+Ra_{sa}\frac{\tilde{T}}{T_{a}}{\hat{\mathbf{e}}}_{z}+{\mathbf{\boldsymbol{\nabla}}}^{2}{\bf
u},$ (130) $\displaystyle\frac{\mathrm{D}\tilde{T}}{\mathrm{D}t}$
$\displaystyle=-\mathcal{D}u_{z}\frac{\tilde{T}}{T_{a}}+\frac{\mathcal{D}}{Ra_{sa}}\dot{\epsilon}:\tau+\nabla^{2}\tilde{T},$
(131)

and boundary conditions

$\displaystyle\frac{\partial\tilde{T}}{\partial x}=0\hskip
28.45274ptz=\pm\frac{1}{2},$ (132)
$\displaystyle\tilde{T}\left(x,z=-\frac{1}{2}\right)-\tilde{T}\left(x,z=\frac{1}{2}\right)=1,$
(133)
$\displaystyle\int_{0}^{L/H}\tilde{P}\left(x,z=\pm\frac{1}{2}\right)\mathrm{d}x=0.$
(134)

###  A.4 Simple Compressible Approximation SCA

The governing equations are

$\displaystyle{\mathbf{\boldsymbol{\nabla}}}\cdot{\bf u}$ $\displaystyle=0,$
(135) $\displaystyle{\bf 0}$
$\displaystyle=-\frac{Ra_{sa}}{\mathcal{D}}{\mathbf{\boldsymbol{\nabla}}}\tilde{P}+Ra_{sa}\tilde{T}{\hat{\mathbf{e}}}_{z}+{\mathbf{\boldsymbol{\nabla}}}^{2}{\bf
u},$ (136) $\displaystyle\frac{\mathrm{D}\tilde{T}}{\mathrm{D}t}$
$\displaystyle=-\mathcal{D}u_{z}\tilde{T}+\frac{\mathcal{D}}{Ra_{sa}}\dot{\epsilon}:\tau+\nabla^{2}\tilde{T},$
(137)

and boundary conditions

$\displaystyle\frac{\partial\tilde{T}}{\partial x}=0\hskip
28.45274ptz=\pm\frac{1}{2},$ (138)
$\displaystyle\tilde{T}\left(x,z=-\frac{1}{2}\right)-\tilde{T}\left(x,z=\frac{1}{2}\right)=1,$
(139)
$\displaystyle\int_{0}^{L/H}\tilde{P}\left(x,z=\pm\frac{1}{2}\right)\mathrm{d}x=0.$
(140)

##  Appendix B Tables of simulation parameters

The parameters of infinite Prandtl number, $4\sqrt{2}$ aspect ratio,
simulations are the following

| $Ra_{sa}$ | $\mathcal{D}$ | $\alpha_{0}T_{0}$ | $\epsilon$ | $n_{x}$ | $n_{z}$ | model | duration
---|---|---|---|---|---|---|---|---
12 runs | $10^{4}$ | $1.5$ |  $1.0$/$0.5$ | 0.1 | 128 | 32 | FC/AA | 0.3
|  |  | 0.1 |  |  |  | ALA/SCA |
5 runs | $10^{6}$ | 0/0.05/0.1 | 1.0 |  | 256 | 64 | AA | 0.5
|  | 0.2/0.4 |  |  |  |  |  |
5 runs | $10^{8}$ | 0/0.05/0.1 | 1.0 |  | 512 | 128 | AA | $10^{-2}$
|  | 0.2/0.4 |  |  |  |  |  |
2 runs | $10^{7}$ | 0.05/0.2 | 1.0 |  | 256 | 64 | AA | $0.1$
4 runs | $3.0\,\times\,10^{5}$ | 1.2 | 1.0 | 0.1 | 256 | 64 | FC/AA | $0.5$
|  |  |  |  |  |  | ALA/SCA |
192 runs | $10^{3/3.5/4/4.5}$ | 0.25/0.5/0.75 | 1.0 | 0.1 | 256 | 64 | FC/AA | $1.0$
| $10^{5/5.5/6/6.5}$ | 1.0/1.25/1.5 |  |  |  |  | ALA/SCA |
17 runs | $10^{7}$ | 0.001/0.01/0.02 | 1.0 |  | 512 | 128 | AA | 0.1
|  | 0.03/0.04/0.05 |  |  |  |  |  |
|  | 0.07/0.1/0.2/0.3 |  |  |  |  |  |
|  | 0.5/0.7/1.0/1.2 |  |  |  |  |  |
|  | 1.4/1.6/1.8 |  |  |  |  |  |
8 runs | $10^{8}$ | 0.05/0.1/0.2/0.4 | 1.0 |  | 1024 | 256 | AA | 0.03
|  | 0.8/1.2/1.6/1.8 |  |  |  |  |  |
8 runs | $10^{9}$ | 0.05/0.1/0.2/0.4 | 1.0 |  | 2048 | 512 | AA | 0.01
|  | 0.8/1.2/1.6/1.8 |  |  |  |  |  |

Two additional runs have been performed with a Prandlt number equal to $10$ in
a cavity of aspect ratio $1$, with the following parameters

| $Ra_{sa}$ | $\mathcal{D}$ | $\alpha_{0}T_{0}$ | $n_{x}$ | $n_{z}$ | model | duration | lateral
---|---|---|---|---|---|---|---|---
|  |  |  |  |  |  |  | condition
1 run | $10^{7}$ | 1.6 | 1.0 | 512 | 128 | AA | 0.1 | x-periodic
|  |  |  |  |  |  |  | (Fourier)
1 run | $10^{7}$ | 1.6 | 1.0 | 512 | 128 | AA | 0.1 | wall-bounded
|  |  |  |  |  |  |  | (SinCos)

##  Appendix C Dimensional anelastic heat flux and dissipation profiles

From the dimensional anelastic equations for a general equation of state, we
derive expressions for the horizontal and time average of the vertical heat
flux and dissipation profile. The anelastic equations are

$\displaystyle{\bf\boldsymbol{\nabla}}\cdot\left(\rho_{a}{\bf u}\right)$
$\displaystyle=0,$ (141) $\displaystyle\rho_{a}\frac{\mathrm{D}{\bf
u}}{\mathrm{D}t}$
$\displaystyle=-\rho_{a}{\bf\boldsymbol{\nabla}}\left(\frac{\tilde{P}}{\rho_{a}}\right)+\frac{\rho_{a}\alpha_{a}T_{a}g}{c_{pa}}\tilde{s}{\bf\hat{e}}_{z}+{\bf\boldsymbol{\nabla}}\cdot\tau.$
(142)
$\displaystyle\rho_{a}\frac{\mathrm{D}\left(T_{a}\tilde{s}\right)}{\mathrm{D}t}$
$\displaystyle=-\frac{\rho_{a}\alpha_{a}T_{a}g}{c_{pa}}u_{z}\tilde{s}+\dot{\epsilon}:\tau-{\bf\boldsymbol{\nabla}}\cdot\left(\phi_{a}+\tilde{\phi}\right),$
(143)

where $\phi_{a}$ and $\tilde{\phi}$ are the conduction heat flux along the
adiabat and the superadiabatic temperature, respectively. The scalar product
of the Navier-Stokes equation is averaged horizontally and temporally (denoted
by an overline) in the assumption of a statistically stationary flow, to
obtain the dissipation profile (after integrating by parts the last term)

$\overline{\dot{\epsilon}:\tau}(z)=\frac{\rho_{a}\alpha_{a}T_{a}g}{c_{pa}}\overline{u_{z}\tilde{s}}-\frac{\mathrm{d}}{\mathrm{d}z}\left[\rho_{a}\overline{\frac{u^{2}}{2}u_{z}}+\overline{\tilde{P}u_{z}}-\overline{u_{i}\tau_{iz}}\right].$
(144)

Taking the horizontal and temporal average of the energy equation (143),
eliminating $\overline{\dot{\epsilon}:\tau}$ using (144), shows that the
following function is independent of $z$ while it is obviously equal to the
heat flux at the top and bottom

$Q_{AA}(z)=\rho_{a}T_{a}\overline{u_{z}\tilde{s}}+\left[\rho_{a}\overline{\frac{u^{2}}{2}u_{z}}+\overline{\tilde{P}u_{z}}-\overline{u_{i}\tau_{iz}}\right]+\phi_{a}+\overline{\phi}.$
(145)

If the heat flux components in brackets converge toward zero, for instance
when the Rayleigh number increases to large values, then the main part of the
flux is carried by the entropy flux $\rho_{a}T_{a}\overline{u_{z}\tilde{s}}$,
except in small layers at the top and bottom where thermal conduction can
compete. In the statically stationary case, the heat flux is uniform along
$z$. From (144), it can be seen that the dissipation profile converges toward
${\rho_{a}\alpha_{a}T_{a}g}/{c_{pa}}\overline{u_{z}\tilde{s}}$, so that the
dissipation profile becomes a well-defined function of height

$\overline{\dot{\epsilon}:\tau}(z)\simeq\frac{\alpha_{a}g}{c_{pa}}Q_{AA},$
(146)

depending on the vertical profiles of $\alpha_{a}$, $c_{pa}$ and $g$.

###### Acknowledgements.

Acknowledgements: thanks are due to Daniel Lecoanet who provided help on how
to use the Dedalus software. The authors are grateful to the LABEX Lyon
Institute of Origins (ANR-10-LABX-0066) of the Université de Lyon for its
financial support within the program “Investissements d’Avenir”
(ANR-11-IDEX-0007) of the French government operated by the National Research
Agency (ANR). Support was provided by the ICMAT Severo Ochoa Project No.
SEV-2015‐0554. J. Curbelo also acknowledges the support of the RyC project
RYC2018-025169, the Spanish grant [PID2020-114043GB-I00] and the Catalan Grant
[No. 2017SGR1049].

Declaration of Interests. The authors report no conflict of interest.

## References

  * Alboussière and Ricard [2013] T. Alboussière and Y. Ricard.  Reflections on dissipation associated with thermal convection.  _Journal of Fluid Mechanics_ , 725:1469–7645, 2013.
  * Alboussière and Ricard [2017] T. Alboussière and Y. Ricard.  Rayleigh-Bénard stability and the validity of quasi-Boussinesq or quasi-anelastic liquid approximations.  _Journal of Fluid Mechanics_ , 817:264–305, 2017.
  * Anders et al. [2019] E.H. Anders, D. Lecoanet, and B. Brown.  Entropy rain: dilution and compression of thermals in stratified domains.  _The Astrophysical Journal_ , 884(65), 2019.
  * Anufriev et al. [2005] A.P. Anufriev, C.A. Jones, and A.M. Soward.  The Boussinesq and anelastic liquid approximations for convection in the Earth’s core.  _Phys. Earth and Planet. Int._ , 152:163–190, 2005\.
  * Backus [1975] G.E. Backus.  Gross thermodynamics of heat engines in deep interior of Earth.  _Proceedings of the National Academy of Sciences_ , 72(4):1555–1558, 1975.
  * Bazarov [1989] I. P. Bazarov.  _Thermodynamique_.  Ed. MIR, 1989.
  * Bénard [1901] H. Bénard.  _Les Tourbillons cellulaires dans une nappe liquide propageant de la chaleur par convection, en régime permanent_.  Thèse de la Faculté des Sciences de Paris, 1901.
  * Bormann [2001] A.S. Bormann.  The onset of convection in the Rayleigh-Bénard problem for compressible fluids.  _Continuum Mech. Thermodyn._ , 13:9–23, 2001.
  * Boussinesq [1903] J. Boussinesq.  _Théorie analytique de la chaleur, tome 2_ , pages 157–161.  Gauthier-Villars, Paris, 1903.
  * Braginsky and Roberts [1995] S.I. Braginsky and P.H. Roberts.  Equations governing convection in earth’s core and the geodynamo.  _Geophys. Astrophys. Fluid Dynam._ , 79:1–97, 1995.
  * Burns et al. [2020] K.J. Burns, G.M. Vasil, J.S. Oishi, D. Lecoanet, and B.P. Brown.  Dedalus: A flexible framework for numerical simulations with spectral methods.  _Phys. Rev. Research_ , 2:023068, 2020.
  * Busse [1967] F. H. Busse.  The stability of finite amplitude cellular convection and its relation to an extremum principle.  _Journal of Fluid Mechanics_ , 30(4):625–649, 1967.
  * Carnot [1824] S. Carnot.  _Réflexions sur la puissance motrice du feu_.  Bachelier Libraire, 1824.
  * Chan and Sofia [1989] K.L. Chan and S. Sofia.  Turbulent compressible convection in a deep atmosphere. IV. Results of three-dimensional computations.  _The Astrophysical Journal_ , 336:1022–1040, 1989.
  * Curbelo et al. [2019] J. Curbelo, L. Duarte, T. Alboussière, F. Dubuffet, S. Labrosse, and Y. Ricard.  Numerical solutions of compressible convection with an infinite Prandtl number: comparison of the anelastic and anelastic liquid models with the exact equations.  _Journal of Fluid Mechanics_ , 873:646–687, 2019.
  * Currie and Browning [2017] L.K. Currie and M.K. Browning.  The magnitude of viscous dissipation in strongly stratified two-dimensional convection.  _The Astrophysical Journal Letters_ , 845(L17), 2017\.
  * Durran [1989] D.R. Durran.  Improving the anelastic approximation.  _J. of the Atmosph. Sci._ , 46(11):1453–1461, 1989.
  * Featherstone and Hindman [2016] N.A. Featherstone and B.W. Hindman.  The spectral amplitude of stellar convection and its scaling in the high-Rayleigh-number regime.  _The Astrophysical Journal_ , 818(1), 2016.
  * Fröhlich et al. [1992] J. Fröhlich, P. Laure, and R. Peyret.  Large departure from Boussinesq approximation in the Rayleigh-Bénard problem.  _Phys. Fluids_ , 4(7):1355–1372, 1992.
  * Fuentes and Cumming [2020] J.R. Fuentes and A. Cumming.  Penetration of a cooling convective layer into a stably-stratified composition gradient: Entrainment at low prandtl number.  _Physical Review Fluids_ , 5(124501), 2020.
  * Garaud [2020] P. Garaud.  Horizontal shear instabilities at low prandtl number.  _The Astrophysical Journal_ , 901(2), 2020.
  * Gebhart [1962] B. Gebhart.  Effects of viscous dissipation in natural convection.  _Journal of Fluid Mechanics_ , 14(2):225–232, 1962.
  * Gilbert [1991] G. T. Gilbert.  Positive definite matrices and sylvester’s criterion.  _The American Mathematical Monthly_ , 98(1), 1991\.
  * Giterman and Shteinberg [1970] M.Sh. Giterman and V.A. Shteinberg.  Criteria of occurrence of free convection in a compressible viscous heat-conducting fluid.  _J. Appl. Math. Mech._ , 34(2):305–311, 1970\.
  * Goluskin et al. [2014] D. Goluskin, H. Johnston, G.R. Flierl, and E.A. Spiegel.  Convectively driven shear and decreased heat flux.  _Journal of Fluid Mechanics_ , 759:360–385, 2014.
  * Grandi and Passerini [2020] D. Grandi and A. Passerini.  Approximation à la Oberbeck-Boussinesq for fluids with pressure-induced stratified density.  _Geophysical & Astrophysical Fluid Dynamics_, pages 1–24, 2020.
  * Grossmann and Lohse [2000] S. Grossmann and D. Lohse.  Scaling in thermal convection: a unifying theory.  _J. Fluid Mech._ , 407:27–56, 2000.
  * Hewitt et al. [1975] J. M. Hewitt, D. P. McKenzie, and N. O. Weiss.  Dissipative heating in convective flows.  _J. Fluid Mech._ , 68(4):721–738, 1975.
  * Horn et al. [2013] S. Horn, O. Shishkina, and C. Wagner.  On non-Oberbeck-Boussinesq effects in three-dimensional Rayleigh-Bénard convection in glycerol.  _J. Fluid Mech._ , 724:175–202, 2013.
  * Howard [1963] L.N. Howard.  Heat transport by turbulent convection.  _Journal of Fluid Mechanics_ , 17(3):405–432, 1963.
  * Howard [1964] L.N. Howard.  Convection at high Rayleigh number.  _Görtler H., Applied Mechanics_ , 1964.
  * Jeffreys [1930] H. Jeffreys.  The instability of a compressible fluid heated below.  _Proc. of the Cambridge Phil. Soc._ , 26(2):170–172, 1930.
  * Jones et al. [2022] C.A. Jones, K. Mizerski, and M. Kessar.  Fully developed anelastic convection with no-slip boundaries.  _Journal of Fluid Mechanics_ , 930(A13), 2022.
  * Käpylä et al. [2019] P.J. Käpylä, M. Viviani, M.J. Käpylä, A. Brandenburg, and F. Spada.  Effects of a subadiabatic layer on convection and dynamos in spherical wedge simulations.  _Geophysical and Astrophysical Fluid Dynamics_ , 113(1–2), 2019.
  * King et al. [2010] S. D. King, C. Lee, P. E. van Keken, W. Leng, S. Zhong, E. Tan, N. Tosi, and M. C. Kameyama.  A community benchmark for 2-D Cartesian compressible convection in the Earth’s mantle.  _Geophysical Journal International_ , 180:73–87, 2010.
  * Lantz and Fan [1999] S.R. Lantz and Y. Fan.  Anelastic magnetohydrodynamic equations for modeling solar and stellar convection zones.  _Astrophys. Journal_ , 121:247–264, 1999.
  * Lecoanet et al. [2014] D. Lecoanet, B.P. Brown, E.G. Zweibel, K.J. Burns, J.S. Oishi, and G.M. Vasil.  Conduction in low Mach number flows. I. Linear and weakly nonlinear regimes.  _The Astrophysical Journal_ , 797(2), 2014.
  * Lipps [1990] F.B. Lipps.  On the anelastic approximation for deep convection.  _J. of the Atmosph. Sci._ , 47(14):1794–1798, 1990.
  * Malkus [1954] W. V. R. Malkus.  The heat transport and spectrum of thermal turbulence.  _Proceedings of the Royal Society A_ , 225:196–212, 1954\.
  * Oberbeck [1879] A. Oberbeck.  Über die Wärmeleitung des Flüssigkeiten bei Berücksichtigung des Strömungen infolge von Temperaturdifferenzen.  _Ann. Phys. Chem._ , 7:271–292, 1879.
  * Ogura and Phillips [1961] Y. Ogura and N.A. Phillips.  Scale analysis of deep and shallow convection in the atmosphere.  _J. Atm. Sci._ , 19:173–179, 1961.
  * Paolucci and Chenoweth [1987] S. Paolucci and D. R. Chenoweth.  Departures from the Boussinesq approximation in laminar Bénard convection.  _Physics of Fluids_ , 30(5):1561–1564, 1987.
  * Rayleigh [1916] J.W.S. Rayleigh.  On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side.  _Phil. Mag. S._ ,  32(192):529–546, 1916.
  * Ricard [2015] Y. Ricard.  Vol. 7, physics of mantle convection.  In Bercovici D. and Schubert G., editors, _Treatise on Geophysics_. Cambridge University Press, 2015.
  * Ricard et al. [2022] Y. Ricard, T. Alboussière, S. Labrosse, J. Curbelo, and F. Dubuffet.  Fully compressible convection for planetary mantles.  _GJI_ , 2022.
  * Roche et al. [2002] P.-E Roche, B Castaing, B Chabaud, and B Hébral.  Prandtl and rayleigh numbers dependences in rayleigh-bénard convection.  _Europhysics Letters (EPL)_ , 58(5):693–698, jun 2002.
  * Schaeffer et al. [2017] N. Schaeffer, D. Jault, H.-C. Nataf, and A. Fournier.  Turbulent geodynamo simulations: a leap towards earth’s core.  _Geophysical Journal International_ , 211:1–29, 2017.
  * Schwarzschild [1906] K. Schwarzschild.  Über das Gleichgewicht des Sonnenatmosphäre.  _Nachr. Kgl. Ges. d. Wiss. zu Gött. Math. Phys. Klasse_ , 1:41–53, 1906.
  * Scott [2001] N.H. Scott.  Thermoelasticity with thermomechanical constraints.  _International Journal of Non-Linear Mechanics_ , 36:549–564, 2001.
  * Sotin and Labrosse [1999] L. Sotin and S. Labrosse.  Three-dimensional thermal convection in an iso-viscous, infinite Prandtl number fluid heated from within and from below: applications to the transfer of heat through planetary mantles.  _Phys. Earth and Planet. Int._ , 112:171–190, 1999\.
  * Spiegel and Veronis [1971] E. A. Spiegel and G. Veronis.  On the boussinesq approximation for a compressible fluid.  _Astrophys. J._ , 131:442–447, 1971.
  * Spiegel [1965] E.A. Spiegel.  Convective instability in a compressible atmosphere. I.  _Astrophys. Journal_ , 141(3):1068–1090, 1965\.
  * Tilgner [2011] A Tilgner.  Convection in an ideal gas at high Rayleigh numbers.  _Physical Review E_ , 84(2):026323, 2011.
  * Vasil et al. [2013] G.M. Vasil, D. Lecoanet, B.P. Brown, T. Wood, and E.G. Zweibel.  Energy conservation and gravity waves in sound-proof treatments of stellar interiors. II. Lagrangian constrained analysis.  _The Astrophysical Journal_ , 773(169), 2013.
  * Viallet et al. [2013] M. Viallet, C. Meakin, D. Arnett, and M. Mocák.  Turbulent convection in stellar interiors. III. Mean-field analysis and stratification effects.  _The Astrophysical Journal_ , 769(1), 2013.
  * Wang [2004] X. Wang.  Infinite Prandtl number limit of Rayleigh-Bénard convection.  _Communications on pure and applied mathematics_ , 57(10):1265–1282, 2004.