# Fully developed anelastic convection with no-slip boundaries

Chris Jones, Krzysztof Mizerski, Mouloud Kessar

<figure><img src="image_1.png" /></figure>

## To cite this version:

Chris Jones, Krzysztof Mizerski, Mouloud Kessar. Fully developed anelastic convection with no-slip boundaries. Journal of Fluid Mechanics, 2022, 930, pp.A13. 10.1017/jfm.2021.905. hal-03658272

HAL Id: hal-03658272

https://u-paris.hal.science/hal-03658272

Submitted on 5 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Fully developed anelastic convection with no-slip boundaries

Chris A. Jones<fnref n="1" />,†, Krzysztof A. Mizerski<fnref n="2" /> and Mouloud Kessar<fnref n="3" />

(Received 27 January 2021; revised 11 August 2021; accepted 9 October 2021)

Anelastic convection at high Rayleigh number in a plane parallel layer with no slip boundaries is considered. Energy and entropy balance equations are derived, and they are used to develop scaling laws for the heat transport and the Reynolds number. The appearance of an entropy structure consisting of a well-mixed uniform interior, bounded by thin layers with entropy jumps across them, makes it possible to derive explicit forms for these scaling laws. These are given in terms of the Rayleigh number, the Prandtl number and the bottom to top temperature ratio, which also measures how much the density varies across the layer. The top and bottom boundary layers are examined and they are found to be very different, unlike in the Boussinesq case. Elucidating the structure of these boundary layers plays a crucial part in determining the scaling laws. Physical arguments governing these boundary layers are presented, concentrating on the case in which the boundary layers are so thin that temperature and density vary little across them, even though there may be substantial temperature and density variations across the whole layer. Different scaling laws are found, depending on whether the viscous dissipation is primarily in the boundary layers or in the bulk. The cases of both high and low Prandtl number are considered. Numerical simulations of no-slip anelastic convection up to a Rayleigh number of 10⁷ have been performed and our theoretical predictions are compared with the numerical results.

**Key words:** Bénard convection, boundary layer structure, stratified turbulence

†Email address for correspondence: c.a.jones@maths.leeds.ac.uk

<footnote n="1">¹Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, UK</footnote>

<footnote n="2">²Department of Magnetism, Institute of Geophysics, Polish Academy of Sciences, ul. Ksiecia Janusza 64, 01-452 Warsaw, Poland</footnote>

<footnote n="3">³Université de Paris, Institut de physique du globe de Paris, CNRS, IGN, F-75005 Paris, France</footnote>
# 1. Introduction

The problem of the influence of density stratification on developed convection is of great importance from the astrophysical point of view. Giant planets and stars are typically strongly stratified and the anelastic approximation, see, e.g. Ogura & Phillips (1962), Gough (1969) and Lantz & Fan (1999), is commonly used to describe convection in their interiors, e.g. Toomre et al. (1976), Glatzmaier & Roberts (1995), Brun & Toomre (2002), Browning et al. (2006), Miesch et al. (2000), Jones & Kuzanyan (2009), Jones et al. (2011), Verhoeven, Wiesehöfer & Stellmach (2015), Kessar et al. (2019) and many others. Note that ‘stratified’ here means that the density is continuously varying and can span several density scale heights, but the layer is convecting rather than stably stratified. The anelastic approximation is based on the convective system being a small departure from the adiabatic state, which is appropriate for large-scale systems with developed, turbulent and, thus, strongly mixing convective regions. The small departure from adiabaticity induces convective velocities much smaller than the speed of sound, so sound waves are neglected in the dynamics. We consider a plane layer of fluid between two parallel plates a distance *d* apart, and we assume that the convection is in a statistically steady state, so that the time averages of time-derivative terms can be neglected. In numerical simulations there are sidewalls or periodic boundary conditions, but here we assume these are chosen to have a negligible influence on the horizontally averaged flow quantities. Most astrophysical applications are in spherical geometry, but the simpler plane layer problem is a natural place to start our investigation of high-Rayleigh-number anelastic convection.

In convection theory, we try to determine the dependencies of the superadiabatic heat flux and the convective velocities (measured by the Nusselt, *Nu*, and Reynolds *Re* numbers) on the driving force measured by the imposed temperature difference between the top and bottom plates, i.e. on the Rayleigh number, *Ra*, and on the Prandtl number, *Pr* (the ratio of the fluid kinematic viscosity to its thermal diffusivity). Here we aim to discover how these dependencies are affected by the stratification. We rely strongly on the theory of Grossmann & Lohse (2000) (further advanced and updated in Stevens et al. 2013) developed for Boussinesq convection, where the density scale height has to be much larger than the fluid depth. However, compressible convection differs strongly from the Boussinesq case, with the latter mostly corresponding to experimental situations. There are two crucial differences, which have very important consequences for the dynamics of convection. Firstly, in the compressible case the viscous heating and the work of the buoyancy force are no longer negligible compared with the heat transport. Secondly, in stratified convection the boundary layers and the flow velocities are different at the top of the layer and the bottom of the layer (Verhoeven et al. 2015), unlike the Boussinesq case where the top and bottom boundary layers have the same structure and the temperature of the well-mixed interior is exactly half-way between the temperature of the top and bottom plates. So although our approach is based on that of Grossmann & Lohse (2000), there are additional novel features required in the compressible convection case. We develop the theory of fully developed convection in stratified systems and study the dependence of the total superadiabatic heat flux and the amplitude of convective flow on the number of density scale heights in the layer. The scaling laws, i.e. the dependencies of the Nusselt and Reynolds numbers on the Rayleigh and Prandtl numbers, are the same as in Boussinesq convection.

In this paper we concentrate on the convective regimes which seem to be most relevant to current numerical capabilities, i.e. the regimes most easily achieved by numerical experiments. These are the regimes where the thermal dissipation is dominated by the thermal boundary layer contribution. These regimes, denoted by *I*<sub>u</sub>, *I*<sub>l</sub>, *II*<sub>u</sub> and *II*<sub>l</sub> on the
phase diagram, figure 2 of Grossmann & Lohse (2000), correspond to Rayleigh numbers less than about 10¹² in the Boussinesq case. It must be noted, however, that contrary to the Boussinesq case, which is well established by numerous experimental and numerical investigations, there are to date no experiments on fully turbulent stratified convection, due to the difficulties of achieving significant density stratification in laboratory settings. Some experiments are being developed using the centrifugal acceleration in rapid rotating systems to enhance the effective gravity (Menaut et al. 2019). There have also been some numerical investigations of anelastic convection in a plane layer, mostly focused on elucidating how well the anelastic approximation performs compared with fully compressible convection, Verhoeven et al. (2015) and Curbelo et al. (2019). This latter paper notes that the top and bottom boundary layer structures that occur in the case of a high Prandtl number are different.

In addition to the dependence on the Rayleigh and Prandtl numbers, our problem depends on how stratified the layer is, which can be estimated by the ratio $\Gamma$ of the temperature at the bottom of the layer $T_B$ to the temperature at the top $T_T$. When $\Gamma$ is close to unity the layer is nearly Boussinesq, but when $\Gamma$ is large there are many temperature and density scale heights within the layer. We aim to derive the functional form of $Nu(\Gamma, Ra, Pr)$ and $Re(\Gamma, Ra, Pr)$, but we cannot derive reliable numerical values for the prefactors in anelastic convection. Since experiments are not available, this will require high resolution high Ra simulations.

In § 2 we derive the anelastic equations and the reference states we use, and outline the structure of high-Rayleigh-number anelastic convection. Further details of the form of the anelastic temperature perturbation are given in Appendix A. In § 3 we derive the energy and entropy production integrals, which are the fundamental building blocks for developing convective scaling laws. In §§ 4 and 5 we derive the physical arguments used to obtain the key properties of the top and bottom boundary layers. In § 6 we derive the scaling laws for the case where the viscous dissipation is primarily in the boundary layers. The case where the dissipation is mainly in the bulk is dealt with in Appendix B. Section 7 gives the results of our numerical simulations, comparing them with our theoretical results. Our conclusions are presented in § 8.

## 2. Fully developed compressible convection under the anelastic approximation

Consider a layer of compressible perfect gas with two parallel boundaries a distance $d$ apart, periodic in the horizontal directions, the evolution of which is described by the set of the Navier-Stokes, mass conservation and energy equations under the anelastic approximation,

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} = -\nabla \left(\frac{p}{\bar{\rho}}\right) + \frac{g}{c_p} s \hat{e}_z + \frac{\mu}{\bar{\rho}} \left[\nabla^2 \mathbf{u} + \frac{1}{3} \nabla (\nabla \cdot \mathbf{u})\right], \quad (2.1) $$

$$ \nabla \cdot (\bar{\rho} \mathbf{u}) = 0, \quad (2.2) $$

$$ \bar{\rho} \bar{T} \left[ \frac{\partial s}{\partial t} + \mathbf{u} \cdot \nabla s \right] = k \nabla^2 T + \mu \left[ q + \frac{\partial u_i}{\partial x_j} \frac{\partial u_j}{\partial x_i} - (\nabla \cdot \mathbf{u})^2 \right], \quad (2.3) $$

$$ \frac{p}{\bar{\rho}} = \frac{\rho}{\bar{\rho}} + \frac{T}{\bar{T}}, \quad s = c_v \frac{p}{\bar{\rho}} - c_p \frac{\rho}{\bar{\rho}}, \quad \gamma = \frac{c_p}{c_v}, \quad c_p - c_v = R, \quad (2.4a-d) $$

where here and below suffix notation is used,

$$ q = \frac{\partial u_i}{\partial x_j} \frac{\partial u_j}{\partial x_i} + \frac{1}{3} (\nabla \cdot \mathbf{u})^2, \quad (2.5) $$
u being the fluid velocity, p the pressure, ρ the density, T the temperature and s the entropy. Barred variables are adiabatic reference state variables, unbarred variables denote the perturbation from the reference state due to convection. The dynamic viscosity μ = ρ̄v, the thermal conductivity k, gravity g and the specific heats at constant pressure, c_p, and constant volume, c_v, are all assumed constant. The bounding plates are no-slip and impenetrable, so u = 0 there. We consider the constant entropy boundary conditions

$$s = \Delta S \quad \text{at } z = 0, \quad s = 0 \quad \text{at } z = d. \qquad (2.6a,b)$$

Note that we do not replace the thermal diffusion term in (2.3) by an entropy diffusion term, as is often done in anelastic approaches. With our no-slip boundaries, there will be boundary layers which may be laminar even at very high Rayleigh number, and entropy diffusion is not appropriate if laminar boundary layers are present. We discuss the additional issues raised by adopting constant temperature boundary conditions rather than constant entropy conditions in Appendix A.

In the anelastic approximation, the full variables, $\hat{p}$, $\hat{\rho}$ and $\hat{T}$, are expanded in terms of the small parameter $\epsilon$, which is defined precisely in (2.10) below, so

$$\left. \begin{aligned} \hat{p} = \bar{p} + \epsilon p, \quad \hat{\rho} = \bar{\rho} + \epsilon \rho, \quad \hat{T} = \bar{T} + \epsilon T, \quad u \sim (\epsilon g d)^{1/2}, \\ t \sim \left(\frac{\epsilon g}{d}\right)^{-1/2}, \quad \hat{s} = \bar{s} + \epsilon s, \end{aligned} \right\} \qquad (2.7)$$

where $\bar{p}$, $\bar{\rho}$ and $\bar{T}$ comprise the adiabatic reference state. Here $\bar{s}$ is simply a constant, and since $s = c_v \ln \hat{p}/\hat{\rho}^\gamma + \text{const.}$, and $\bar{p}/\bar{\rho}^\gamma$ is constant, we obtain (2.4b) by choosing the constant appropriately.

## 2.1. The adiabatic reference state

The reference state is the adiabatic static equilibrium governed by $d\bar{p}/dz = -g\bar{\rho}$, $\bar{p} = R\bar{\rho}\bar{T}$, together with the polytropic law $\bar{p}/p_B = (\bar{\rho}/\rho_B)^\gamma$. Here R is the gas constant, $z = 0$ is the bottom of the layer and $z = d$ the top, and $p_B, \rho_B$ are the pressure and density at the bottom of the layer. It follows that

$$\bar{T} = T_B \left(1 - \theta \frac{z}{d}\right), \quad \bar{\rho} = \rho_B \left(1 - \theta \frac{z}{d}\right)^m, \quad \bar{p} = \frac{gd\rho_B}{\theta(m+1)} \left(1 - \theta \frac{z}{d}\right)^{m+1}, \quad (2.8a-c)$$

$$\frac{gd}{c_p} = \Delta \bar{T} = T_B - T_T > 0, \quad \theta = \frac{\Delta \bar{T}}{T_B}, \quad m = \frac{1}{\gamma - 1}, \quad \Gamma = \frac{T_B}{T_T} = \frac{1}{1 - \theta}, \quad (2.8d-g)$$

which defines $\theta$, and the polytropic index $m$. We use subscripts T and B to denote conditions at the top and bottom boundary, respectively. The temperature ratio, $\Gamma > 1$, is a convenient measure of the compressibility. Here $\Gamma \to 1$ is the Boussinesq limit, and highly compressible layers have $\Gamma$ large. Note that $\Gamma^m = \rho_B/\rho_T$ is the ratio of the highest to lowest density in the layer. The density ratio can be very large in astrophysical applications, the density of the bottom of the solar convection zone being $\sim 10^6$ times the density at the top. Sometimes the number of density scale heights, $N_\rho$ (the scale height being defined at the top of the layer), that fit into the layer is used to measure compressibility, $N_\rho = m(\Gamma - 1)$.

## 2.2. The conduction state

The adiabatic reference state satisfies $\nabla^2 \bar{T} = 0$, but since it is isentropic, it does not satisfy the entropy boundary conditions. The anelastic conduction state is also a polytrope, but
Fully developed anelastic convection with no-slip boundaries

with a slightly more negative temperature gradient, so $\tilde{T}_B = T_B$, $\tilde{T}_T < T_T$. The conduction state is

$$ \tilde{T} = T_B \left(1 - \tilde{\theta} \frac{z}{d}\right), \quad \tilde{\rho} = \rho_B \left(1 - \tilde{\theta} \frac{z}{d}\right)^{\tilde{m}}, \quad \tilde{p} = \frac{gd\rho_B}{\tilde{\theta}(\tilde{m}+1)} \left(1 - \tilde{\theta} \frac{z}{d}\right)^{\tilde{m}+1}, \quad (2.9a-c) $$

$$ \Delta \tilde{T} = T_B - \tilde{T}_T > 0, \quad \tilde{\theta} = \frac{\Delta \tilde{T}}{T_B}, \quad \tilde{m} = \frac{gd}{R \Delta \tilde{T}} - 1. \quad (2.9d-f) $$

The small anelastic parameter $\epsilon$ is now defined as

$$ \epsilon = \tilde{\theta} \frac{\tilde{m} + 1 - \gamma \tilde{m}}{\gamma} = - \frac{d}{T_B} \left[ \frac{d\tilde{T}}{dz} + \frac{g}{c_p} \right] \ll 1, \quad (2.10) $$

and the entropy in the conduction state is

$$ \tilde{s} = c_v \ln \frac{\tilde{p}}{\tilde{\rho}^\gamma} + \text{const.} = \frac{\epsilon c_p}{\theta} \ln \left(1 - \theta \frac{z}{d}\right) + \text{const.}, \quad \text{so } s = \frac{c_p}{\theta} \ln \left(1 - \theta \frac{z}{d}\right) + \text{const.}, \quad (2.11) $$

which is the scaled entropy, see (2.7), correct to $O(\epsilon)$ since $\tilde{\theta}$ and $\theta$ differ by only $O(\epsilon)$. Since the boundaries have fixed entropy, the entropy at the boundaries in the conduction state defines the entropy drop across the layer for all Rayleigh numbers, so

$$ \Delta S = \frac{c_p}{\theta} \ln \Gamma = \frac{c_p \Gamma \ln \Gamma}{\Gamma - 1}, \quad (2.12) $$

relating to the entropy boundary conditions (2.6a,b). Note that as our entropy variable $s$ is scaled by $\epsilon$, the entropy drop is $O(\epsilon)$. Some anelastic papers take the conduction state as the reference state, and some take the adiabatic state as the reference state. Taking the conduction state as the reference state is appropriate when convection near the critical Rayleigh number is considered, but at the large Rayleigh numbers considered here, the conduction state is less relevant. Although the conduction state tends to the adiabatic reference state as $\epsilon \to 0$, the thermodynamic variables are not the same in the two systems: $T=0$ with respect to the adiabatic state corresponds to $T = T_B z/d$ if the conduction state is chosen as the reference state.

In (2.1) we have made use of (2.4b), (2.8a-c) and (2.10) to write (Braginsky & Roberts 1995; Lantz & Fan 1999)

$$ -\frac{\nabla p}{\bar{\rho}} - \frac{\rho}{\bar{\rho}} g \hat{e}_z = -\nabla \left(\frac{p}{\bar{\rho}}\right) + \frac{g}{c_p} s \hat{e}_z + O(\epsilon). \quad (2.13) $$

## 2.3. The Nusselt and Rayleigh numbers in anelastic convection

Next we consider the superadiabatic heat flux. The horizontal average at level $z$ is denoted by $\langle \rangle_h$. At the boundaries, all the heat is carried by conduction, and if the total temperature $\tilde{T} = \tilde{T} + \epsilon T$, then the total heat flux at the boundaries is $-k \frac{d\langle \tilde{T} \rangle_h}{dz} = -k \frac{d\tilde{T}}{dz} - k \epsilon \frac{d\langle T \rangle_h}{dz}$, but the superadiabatic part is obtained by subtracting off the heat flux carried along the adiabat, so

$$ F^{\text{super}} = -k \frac{d \langle T \rangle_h}{dz} \Big|_{z=0}. \quad (2.14) $$

The Nusselt number in anelastic convection is defined as the ratio of the superadiabatic heat flux divided by the heat conducted down the conduction state superadiabatic gradient.
Note that the flux conducted down the adiabatic gradient is ignored in this definition, so

$$Nu = \frac{F^{super}d}{kT_B}, \qquad (2.15)$$

so Nu is close to unity near onset, and is large in fully developed convection. For fixed entropy boundary conditions, the Rayleigh number is defined as

$$Ra = \frac{g\Delta S d^3 \rho_B^2}{\mu k} \approx \frac{c_p \Delta S \Delta \bar{T} d^2 \rho_B^2}{\mu k}. \qquad (2.16)$$

The anelastic approximation is asymptotically valid in the limit $\epsilon \to 0$. Note that a small superadiabatic temperature gradient does not imply a small Rayleigh number Ra, since the diffusion coefficients can be small, in fact to get $Ra \sim O(1)$ in the limit $\epsilon \to 0$ we must have

$$\frac{k}{\rho_B c_p} \sim (gd^3\epsilon)^{1/2}, \quad \frac{\mu}{\rho_B} \sim (gd^3\epsilon)^{1/2}, \qquad (2.17a,b)$$

allowing large but finite Ra even when the superadiabatic gradient is small. To derive the anelastic equations (2.1) to (2.5), we insert (2.7) into the full compressible equations and divide the momentum equation by $\epsilon$, the mass conservation equation by $\epsilon^{1/2}$ and the energy equation by $\epsilon^{3/2}$. Having taken this limit, the parameter $\epsilon$ no longer appears in our analysis. However, if anelastic work is compared with fully compressible situations, then a finite value of $\epsilon$ must be chosen, and the anelastic results are only approximate, though there is a growing body of evidence that the anelastic approximation does capture the main features of subsonic compressible convection.

## 2.4. High-Rayleigh-number convection

We have the following physical picture in mind. In strongly turbulent convection we expect the entropy $s$ to be well-mixed away from boundary layers near $z = 0$, $d$. We denote the global spatial average over the convecting layer by $|| ||$ and the horizontal average at level $z$ by $\langle \rangle_h$. The total entropy drop is the conduction state value $\Delta S = c_p \ln \Gamma/\theta$. Since the entropy is constant in the bulk interior, we define the entropy drops $\Delta S_B$ and $\Delta S_T$ across the bottom and top boundary layers, respectively. These will not be equal, with $\Delta S_T$ normally considerably larger than $\Delta S_B$. We must however have

$$\Delta S_B + \Delta S_T = \Delta S. \qquad (2.18)$$

We consider only the case where both the top and bottom boundary layers are laminar. At extremely high Ra these layers may become turbulent as can happen in the Boussinesq case. The laminar boundary layer case is simpler, and gives predictions which can be broadly compared with numerical simulations, though it is difficult for numerical simulations to get into the fully developed large Rayleigh and Nusselt number regime we are aiming at here. A schematic picture of the horizontally averaged entropy profile expected in highly supercritical anelastic convection is shown in figure 1(a).

Since the heat flux through the boundary layers is determined by thermal diffusion rather than entropy diffusion, we need to express the temperature jumps across the thermal boundary layers in terms of the entropy jumps. From (2.4a–d) we obtain

$$\frac{(\Delta \rho)_i}{\rho_i} \approx \frac{1}{\gamma - 1} \left[ \frac{(\Delta T)_i}{T_i} - \gamma \frac{(\Delta s)_i}{c_p} \right], \quad \frac{(\Delta p)_i}{p_i} \approx \frac{\gamma}{\gamma - 1} \left[ \frac{(\Delta T)_i}{T_i} - \frac{(\Delta s)_i}{c_p} \right], \qquad (2.19a,b)$$

where the $\Delta$ quantities refer to the jump in the horizontally averaged value across the boundary layer and the subscript $i$ stands either for $B$ or $T$. We also define the thermal
<figure><img src="image_4.png" /><figcaption>Figure 1. (a) A schematic picture of the entropy profile in developed convection. (b) A schematic picture of the anelastic temperature perturbation in developed convection.</figcaption></figure>

and viscous boundary layer thicknesses, $\delta_i^{th}$ and $\delta_i^{v}$ with $i = B, T$, which we use to obtain scaling estimates. Numerical simulations indicate that the horizontal velocity

$$U_H = (\langle u_x^2 \rangle_h + \langle u_y^2 \rangle_h)^{1/2} \qquad (2.20)$$

has local maxima close to both boundaries (see, e.g. figure 3(b) below), so these maxima are a convenient way to define the velocity jumps across the viscous boundary layers, $\Delta U_i$, so

$$\Delta U_B = U_H(z = z_{max,B}), \quad \Delta U_T = U_H(z = z_{max,T}), \qquad (2.21a,b)$$

where $z = z_{max,B}$, $z = z_{max,T}$ are the locations of the local maxima. The thermal boundary layer thickness for the entropy, $\delta_i^{th}$, and the viscous boundary layer thickness, $\delta_i^{v}$, can be defined as

$$\delta_i^{th} = \left\{ -\frac{1}{\Delta S_i} \frac{d \langle S \rangle_h}{dz} \bigg|_{z=z_i} \right\}^{-1}, \quad \delta_i^{v} = \left\{ \pm \frac{1}{\Delta U_i} \frac{dU_H}{dz} \bigg|_{z=z_i} \right\}^{-1}, \quad z_i = 0, d. \qquad (2.22)$$

In the boundary layers, the dominant balance in the z-component of the Navier-Stokes equation occurs between the pressure gradient and the buoyancy force. Mass conservation in the boundary layers means $u_{z,i} \sim O(\delta_i^{v} U_i / d)$ so the vertical component of inertia is small. The boundary layers are therefore approximately hydrostatic,

$$(\Delta p)_i \approx \frac{g}{c_p} \rho_i \Delta s_i \delta_i^{th}. \qquad (2.23)$$

Inserting (2.23) into (2.19b) leads to

$$\frac{(\Delta T)_i}{T_i} \approx \frac{(\Delta s)_i}{c_p} \left( 1 + \theta \frac{\delta_i^{th} T_B}{d T_i} \right). \qquad (2.24)$$

Typically the term $(\theta \delta_i^{th} T_B) / (T_i d)$ resulting from the pressure jump across the boundary layers is expected to be small because the boundary layer is thin. However, in simulations where the Rayleigh number is bounded above by numerical constraints, the top boundary layer may not be as thin as desired for accurate asymptotics to apply, and the factor $T_B / T_T$
can be large in layers containing many scale heights. We refer to the case where the pressure term is not negligible as the compressible boundary layer case. However, in this work we assume the boundary layers are incompressible, which is valid provided $T_B/T_T$ remains finite as the Rayleigh number increases and the boundary layers become very thin. Then the pressure term makes a negligible contribution in both boundary layers, so that in these boundary layers

$$ \frac{(\Delta T)_i}{T_i} \approx \frac{(\Delta s)_i}{c_p} \quad \text{and} \quad \frac{(\nabla T)_i}{T_i} \approx \frac{(\nabla s)_i}{c_p}. \qquad (2.25a,b) $$

Defining the temperature boundary layer thicknesses similarly to (2.22),

$$ \delta_i^T = \left\{ - \frac{1}{\Delta T_i} \frac{d \langle T \rangle_h}{dz} \Big|_{z=z_i} \right\}^{-1}, \qquad (2.26) $$

the temperature boundary layer thicknesses are the same as the entropy boundary layer thicknesses. Note that in the compressible boundary layer case the entropy and temperature boundary layer thicknesses will be different. For incompressible boundary layers and high Rayleigh number, the Nusselt number can be written in terms of the boundary layer thicknesses, using (2.15), (2.14), (2.26) and (2.25a,b),

$$ Nu = \frac{d \Delta S_T}{\delta_T^{th} \Gamma c_p} = \frac{d \Delta S_B}{\delta_B^{th} c_p}. \qquad (2.27) $$

In figure 1(b) we sketch the horizontally averaged anelastic temperature perturbation $\langle T \rangle_h$. This is sometimes referred to as the superadiabatic temperature (e.g. Verhoeven et al. 2015). Note that with our constant entropy boundary conditions, $\langle T \rangle_h$ is not zero at the boundaries. We show in Appendix A that the offsets, $\langle T \rangle_{hB}$ at $z=0$ and $\langle T \rangle_{hT}$ at $z=d$, are both positive and we show also that in the bulk, turbulent pressure effects make the gradient of $\langle T \rangle_h$ positive as shown in figure 1(b). Since $T$ is the difference between the actual temperature and the adiabatic reference state temperature, it might be assumed that if $d\langle T \rangle_h/dz > 0$ the layer is stably stratified. However, we must be careful not to confuse the anelastic temperature with the potential temperature, which in our problem is $T_B(T/\bar{T} - (\gamma - 1)p/\gamma\bar{p}) = sT_B/c_p$ where the reference temperature of potential temperature is taken as $T_B$. It is the potential temperature gradient which must be negative for convective instability, not the anelastic temperature gradient.

To obtain the anelastic temperature fluctuation as sketched in figure 1(b), we need to make use of (2.4a) and (2.4b), so we need to make a specific choice for entropy at the boundaries. Here we have chosen to take the entropy at the top boundary as zero, so the entropy at the bottom boundary is $s = \Delta S$. A different choice of entropy constant adds an easily found function of $z$ to $T$, $\rho$ and $p$ but this does not affect the velocity field obtained. One further point is that if (2.1) is horizontally averaged, the horizontal average satisfies a first-order differential equation in $z$ (see Appendix A for details), so a boundary condition on $\langle p \rangle_h$ is required. Here we choose the natural condition that mass is conserved over the layer, so that

$$ ||\rho|| = 0 \Rightarrow \langle p \rangle_{h,T} = \langle p \rangle_{h,B}. \qquad (2.28) $$

This means that the total mass in the layer is the same as in the adiabatic reference state. As we see in Appendix A, this means the horizontally averaged anelastic pressure perturbations at the top and bottom of the layer must be equal.
### 3. Energy and entropy production integrals

Understanding the energy transfer and entropy production in convective flow is the key to understanding the physics of compressible convection. Therefore, we derive now a few exact relations which will allow us to study some general aspects of the dynamics of developed compressible convection. We assume any initial transients in the convection have been eliminated, and we are in a statistically steady state. Formally, this means we consider time-averaged quantities throughout the paper.

#### 3.1. Energy balance

By multiplying the Navier–Stokes equation (2.1) by $\bar{\rho}u$ and averaging over the entire volume (recalling that horizontal averages of $x$ and $y$ derivatives are zero), we obtain the relation

$$ \frac{g}{c_p} ||\bar{\rho} u_z s|| = \mu ||q||, \qquad (3.1) $$

stating that the total rate of working per unit volume of the buoyancy force is equal to the total viscous dissipation in the fluid per unit volume. In deriving (3.1) use has been made of the no-slip boundary conditions and (2.2).

We derive the superadiabatic heat flux in the system at every $z$ by averaging over a horizontal plane and integrating the heat equation (2.3) from 0 to $z$,

$$ F^{\text{super}} = \langle \bar{\rho} \bar{T} u_z s \rangle_h - k \frac{d}{dz} \langle T \rangle_h + \frac{g}{c_p} \int_0^z \langle \bar{\rho} u_z s \rangle_h dz \\ - \mu \int_0^z \langle q \rangle_h dz - 2\mu \left[ \left\langle u_z \frac{du_z}{dz} \right\rangle_h - \frac{m \Delta \bar{T}}{\bar{T} d} \left\langle u_z^2 \right\rangle_h \right]. \qquad (3.2) $$

In deriving this expression, we have made use of (2.8a,d) and

$$ \begin{aligned} \frac{\partial u_i}{\partial x_j} \frac{\partial u_j}{\partial x_i} - (\nabla \cdot \mathbf{u})^2 &= \frac{\partial^2}{\partial x_i \partial x_j} (u_i u_j) - 2 \frac{\partial}{\partial x_j} (u_j \nabla \cdot \mathbf{u}) \\ &= \frac{\partial^2}{\partial x_i \partial x_j} (u_i u_j) + 2 \frac{\partial}{\partial x_j} \left( u_j u_z \frac{m \Delta \bar{T}}{\bar{T}} \right) \end{aligned} \qquad (3.3) $$

since the continuity equation gives

$$ \nabla \cdot \mathbf{u} = \frac{u_z m \Delta \bar{T}}{\bar{T} d}, \qquad (3.4) $$

and we recall that $x$ or $y$ derivatives vanish on horizontal averaging. As $z \to d$, all terms with a factor $u_z$ tend to zero, so on using (3.1) we obtain the overall energy balance equation,

$$ F^{\text{super}} = -k \left. \frac{d \langle T \rangle_h}{dz} \right|_{z=0} = -k \left. \frac{d \langle T \rangle_h}{dz} \right|_{z=d}; \qquad (3.5) $$

thus, in a stationary state the heat flux entering the fluid volume at $z = 0$ must be equal to the outgoing heat flux at $z = d$.
## 3.2. Entropy balance

This energy balance equation alone is not very helpful for evaluating the Nusselt number. We need the entropy balance equation, obtained by dividing the energy equation (2.3) by $\bar{T}$, horizontally averaging, and integrating from 0 to $z$,

$$ \langle \bar{\rho} u_z s \rangle_h = - \frac{k}{T_B} \left. \frac{d}{dz} \langle T \rangle_h \right|_{z=0} + \frac{k}{\bar{T}} \left. \frac{d}{dz} \langle T \rangle_h \right|_z - \int_0^z \frac{k \Delta \bar{T}}{\bar{T}^2 d} \frac{d}{dz} \langle T \rangle_h dz + \int_0^z \frac{\mu}{\bar{T}} \langle q \rangle_h dz \\ + \int_0^z \frac{\mu}{\bar{T}} \left\langle \frac{\partial^2}{\partial x_i \partial x_j} (u_i u_j) - 2 \frac{\partial}{\partial x_j} \left( u_j \frac{\partial u_i}{\partial x_i} \right) \right\rangle_h dz, \quad (3.6) $$

where use has been made of (3.3). The overall entropy balance equation comes from taking the limit $z \to d$ of (3.6), noting $u_z \to 0$ in this limit,

$$ \left. \frac{k}{T_B} \frac{d}{dz} \langle T \rangle_h \right|_{z=0} - \left. \frac{k}{T_T} \frac{d}{dz} \langle T \rangle_h \right|_{z=d} = - \int_0^d \frac{k \Delta \bar{T}}{\bar{T}^2 d} \frac{d}{dz} \langle T \rangle_h dz + \int_0^d \frac{\mu}{\bar{T}} \langle q \rangle_h dz. \quad (3.7) $$

Equations (3.6), (3.7) look complicated, but they simplify considerably when the top and bottom boundary layers are thin. We start with (3.6), which we write as

$$ S_{\text{conv}} = S_{\text{diff}} + S_{\text{visc}}, \quad (3.8) $$

where

$$ S_{\text{conv}} = \langle \bar{\rho} u_z s \rangle_h, \quad (3.9) $$

the net entropy carried out of the region (0, z) by the convective velocity at level z,

$$ S_{\text{diff}} = - \frac{k}{T_B} \left. \frac{d}{dz} \langle T \rangle_h \right|_{z=0} + \frac{k}{\bar{T}} \left. \frac{d}{dz} \langle T \rangle_h \right|_z - \int_0^z \frac{k \Delta \bar{T}}{\bar{T}^2 d} \frac{d}{dz} \langle T \rangle_h dz \quad (3.10) $$

so that $S_{\text{diff}}$ is the entropy balance in region (0, z) of the entropy change due to thermal diffusion. This is divided into the first term, which represents the positive entropy being conducted into our region at the bottom boundary (the gradient of $\langle T \rangle_h$ is negative there; see figure 1b), the second term is the entropy conducted across level z, and the third term is the entropy production by internal diffusion in our given region. By looking at figure 1(b) it is apparent that when the boundary layers are thin, the first term is much larger than the other two except when z is in the boundary layers. If z is in the bulk, the gradient of $\langle T \rangle_h$ is $O(\Delta \bar{T}/d)$ whereas at the boundary it is $O(\Delta \bar{T}/\delta^{th})$, much larger since the boundary layer is thin. The third term is always small compared with the first, because the gradient is $O(\Delta \bar{T}/d)$ outside the boundary layers. The integrand is of order $O(\Delta \bar{T}/\delta^{th})$ in the boundary layers, but because they are thin this only contributes a small amount to the integral. So when the boundary layers are thin

$$ S_{\text{diff}} \approx - \frac{k}{T_B} \left. \frac{d}{dz} \langle T \rangle_h \right|_{z=0} \quad \text{if } z \text{ is in the bulk,} \quad (3.11) $$

$$ S_{\text{diff}} \approx - \frac{k}{T_B} \left. \frac{d}{dz} \langle T \rangle_h \right|_{z=0} + \frac{k}{T_T} \left. \frac{d}{dz} \langle T \rangle_h \right|_{z=d} \quad \text{if } z=d. \quad (3.12) $$

We now turn to

$$ S_{\text{visc}} = \int_0^z \frac{\mu}{\bar{T}} \langle q \rangle_h dz + \int_0^z \frac{\mu}{\bar{T}} \left\langle \frac{\partial^2}{\partial x_i \partial x_j} (u_i u_j) - 2 \frac{\partial}{\partial x_j} \left( u_j \frac{\partial u_i}{\partial x_i} \right) \right\rangle_h dz. \quad (3.13) $$
Because of the horizontal averaging, and using (2.2), (2.8a-c) and (3.4), the second integral can be written as

$$ \int_0^z \frac{\mu}{\bar{T}} \left\langle \frac{\partial^2}{\partial x_i \partial x_j} (u_i u_j) - 2 \frac{\partial}{\partial x_j} \left( u_j \frac{\partial u_i}{\partial x_i} \right) \right\rangle_h dz \\ = \int_0^z \frac{\mu}{\bar{T}} \left[ \frac{\partial^2}{\partial z^2} \langle u_z^2 \rangle_h - \frac{m \Delta \bar{T}}{\bar{T} d} \frac{\partial}{\partial z} \langle u_z^2 \rangle_h - \frac{m (\Delta \bar{T})^2}{\bar{T}^2 d^2} \langle u_z^2 \rangle_h \right] dz. \quad (3.14) $$

We now consider the magnitude of the terms in (3.13). If the root-mean-square velocity in the bulk is $U$, we expect the horizontal velocity to vary from 0 to $O(U)$ across the boundary layers of thickness $\delta_i^v$, so the velocity gradients appearing in $q$ are of magnitude $O(U/\delta_i^v)$. Therefore, $q$ itself is $O(U^2/(\delta_i^v)^2)$, and since the boundary layers are thin their contribution to the first integral in $S_{visc}$ is $O(\mu U^2/\bar{T}\delta_i^v)$. In the boundary layers $u_z$ is small, but in the bulk we expect $u_z$ to be $O(U)$ and so $\langle u_z^2 \rangle_h$ is $O(U^2)$. Because $\langle u_z^2 \rangle_h$ is horizontally averaged, it will vary on a length scale $d$ with $z$, so the gradient of $\langle u_z^2 \rangle_h$ is $O(U^2/d)$ and the second derivative is $O(U^2/d^2)$. From (3.14), the order of magnitude of the second term in (3.13) is then $O(\mu U^2/Td)$, which is $O(\delta_i^v/d)$ smaller than the contribution from the term due to $q$. Therefore, provided the Rayleigh number is high enough for the boundary layers to be thin, (3.6) is asymptotically equivalent to

$$ \langle \bar{\rho} u_z s \rangle_h \sim - \frac{k}{T_B} \frac{d}{dz} \langle T \rangle_h \Big|_{z=0} + \int_0^z \frac{\mu}{\bar{T}} \langle q \rangle_h dz. \quad (3.15) $$

when $z$ is in the bulk. Note that this simplification still holds if the dissipation in the bulk is larger than the dissipation in the boundary layers, which can happen, as noted by Grossmann & Lohse (2000). When the dissipation is primarily in the boundary layers, the left-hand side of (3.15) is constant in the bulk, which we exploit later. In either case, as $z \to d$, we get

$$ \frac{k}{T_B} \frac{d}{dz} \langle T \rangle_h \Big|_{z=0} - \frac{k}{T_T} \frac{d}{dz} \langle T \rangle_h \Big|_{z=d} = \frac{F^{\text{super}} \Delta \bar{T}}{T_B T_T} \sim \int_0^d \frac{\mu}{\bar{T}} \langle q \rangle_h dz. \quad (3.16) $$

Note also that in these expressions, the term $(\nabla \cdot \mathbf{u})^2/3$ in (2.5) makes a negligible contribution to (3.16) compared with the gradient terms, by using (3.4).

### 3.3. The Boussinesq limit

At first sight, it appears that our entropy balance formulation of the equation for dissipation (3.7) is fundamentally different from that used by Grossmann & Lohse (2000) in the Boussinesq case, who start from (2.3) of Siggia (1994),

$$ (Nu - 1)Ra = \left\langle \left( \frac{\partial u_i}{\partial x_j} \right)^2 \right\rangle, \quad \text{or} \quad g\alpha\kappa \Delta T (Nu - 1) = \nu \int_0^d \langle q \rangle_h dz \quad (3.17) $$

in our dimensional units. Here $\Delta T$ is the superadiabatic temperature difference between the boundaries. In the Boussinesq limit $\Gamma \to 1$, the basic state temperature and density tend to constant values $T_B$ and $\rho_B$, respectively, so the thermal diffusivity $\kappa$ and kinematic viscosity $\nu$ are constants, $\kappa = k/\rho_B c_p$ and $\nu = \mu/\rho_B$. For a perfect gas, the coefficient of expansion $\alpha = 1/T_B$. In this subsection we show that (3.17) is in fact the Boussinesq
limit of the entropy balance equation (3.7), which we use to derive the scaling laws in § 6 below. Our entropy balance formulation is a generalization of the Grossmann & Lohse (2000) formulation, which is now seen to be a limiting case of the more general entropy balance approach.

Following Spiegel & Veronis (1960), from the z-component of (2.1), $p/\bar{\rho}d \sim gs/c_p$, so

$$ \frac{p}{\bar{p}} \sim \frac{g\bar{\rho}d}{\bar{p}} \frac{s}{c_p} \sim \frac{d}{H} \frac{s}{c_p}, \qquad (3.18) $$

where $H$ is the pressure scale height. In the Boussinesq limit $d/H$ becomes small; the Boussinesq limit $\Gamma \to 1$ is the thin layer limit (Spiegel & Veronis 1960). Then (2.4a,b) become

$$ \frac{T}{T_B} \sim -\frac{\rho}{\rho_B}, \quad \frac{s}{c_p} \sim \alpha T, \qquad (3.19a,b) $$

so the entropy variable becomes equivalent to the superadiabatic temperature variable. The entropy jump $\Delta S$ across the layer can be written as a superadiabatic temperature jump $\Delta T = \Delta S/\alpha c_p$. As $\Gamma \to 1$, from (2.12), $\Delta S/c_p \to 1$, so $\Delta T \to 1/\alpha = T_B$. From energy conservation, (3.5), the gradient of $\langle T \rangle_h$ is the same at the top and bottom of the layer, so (3.7) can be written as

$$ -\frac{\Delta \bar{T}}{T_B^2} \left. \frac{d}{dz} \langle T \rangle_h \right|_{z=0} = -\frac{k \Delta \bar{T}}{T_B^2 d} \int_0^d \frac{d}{dz} \langle T \rangle_h dz + \frac{\mu}{T_B} \int_0^d \langle q \rangle_h dz, \qquad (3.20) $$

using the constancy of $\bar{T}$ in the Boussinesq limit. From (2.14) and (2.15),

$$ Nu = - \left. \frac{d}{T_B} \frac{d}{dz} \langle T \rangle_h \right|_{z=0} \quad \text{or} \quad Nu = - \left. \frac{d}{\Delta T} \frac{d}{dz} \langle T \rangle_h \right|_{z=0}, \qquad (3.21) $$

which is the familiar form of the Boussinesq Nusselt number, the ratio of the total heat flux at the bottom to the conducted heat flux $-k\Delta T/d$. From (2.8d), $\Delta \bar{T}$ can be written as $gd/c_p$, so (3.20) becomes

$$ \frac{k}{c_p T_B} \mathrm{Nug}\alpha \Delta T = \frac{k}{c_p T_B} g\alpha \Delta T + \frac{\mu}{T_B} \int_0^d \langle q \rangle_h dz \quad \text{or} \quad (Nu - 1)ga\Delta T\kappa = \nu \int_0^d \langle q \rangle_h dz, \qquad (3.22) $$

which is (3.17), showing that the dissipation integral, which plays a key role in Grossmann & Lohse's (2000) analysis, is indeed the Boussinesq limit of the entropy balance equation (3.7).

## 4. The boundary layers and Prandtl number effects

As in the Boussinesq case, the thermal and viscous boundary layers can be nested inside each other when the Prandtl number is different from unity.

A central idea in the theory of fully developed Boussinesq convection is based on the assumption that the structure of turbulent convective flow is always characterised by the presence of a large-scale convective roll called the wind of turbulence, Grossmann & Lohse (2000). This idea, which in the non-stratified case stems from vast numerical and experimental evidence, is retained in the case of anelastic convection. However, the significant stratification in the anelastic case breaks the Boussinesq up-down symmetry
and so we must distinguish between the magnitude of the wind of turbulence near the bottom of the bulk and its magnitude near the top of the bulk, denoted by $U_B$ and $U_T$, respectively, which can now significantly differ. So the label $U_i$ can denote either the horizontal velocity just outside the top or bottom viscous boundary layers. We also assume that this wind of turbulence forms boundary layers with a horizontal length scale comparable to the layer depth $d$. It is of course an assumption that such layers form in anelastic convection, but they are observed to occur in incompressible flow, and the limited simulations we have available gives this idea some support. Whereas the results in §§ 2 and 3 are asymptotically valid in the anelastic framework in the limit of large $Ra$, what follows is dependent on the Grossmann & Lohse (2000) approach being valid for anelastic convection.

The Prandtl number is a constant in this problem, given by

$$ Pr = \frac{\mu c_p}{k}. \qquad (4.1) $$

In figure 2(a) the high-Prandtl-number case is shown, with the thinner thermal boundary layer nested inside the viscous boundary layer. The velocity at the edge of the thermal boundary layer is then $\delta_i^{th}U_i/\delta_i^v$, assuming the velocity falls off linearly inside the viscous boundary layer over the thinner thermal boundary layer. In the boundary layers, advection balances diffusion, so from the momentum equation (2.1) we estimate that

$$ \frac{\rho_i U_i^2}{d} \sim \frac{\mu U_i}{(\delta_i^v)^2} \quad \text{so} \quad U_i \sim \frac{\mu d}{\rho_i (\delta_i^v)^2}. \qquad (4.2) $$

For the entropy boundary layers, from (2.3),

$$ \frac{\rho_i T_i U_i s}{d} \frac{\delta_i^{th}}{\delta_i^v} \sim k \nabla^2 T \approx \frac{k T_i}{c_p} \nabla^2 s \sim \frac{k T_i s}{c_p (\delta_i^{th})^2}, \quad \text{so} \quad U_i = \frac{k d \delta_i^v}{\rho_i c_p (\delta_i^{th})^3}, \qquad (4.3) $$

where (2.25a,b) has been used and the factor $\delta_i^{th}/\delta_i^v$ arises because the horizontal velocity is reduced since the entropy boundary layer is thinner than the viscous boundary layer. Dividing (4.2) by (4.3) we obtain

$$ \frac{\delta_i^v}{\delta_i^{th}} \sim Pr^{1/3}, \qquad (4.4) $$

giving the ratio of the viscous to thermal boundary layer thickness for the high-Prandtl-number case. Note that although the top and bottom boundary layers have different thicknesses, this ratio is the same for both layers.

For the low-Prandtl-number case, the viscous boundary layer lies inside the thermal boundary layer; see figure 2(b). Now the velocity at the edge of the thermal boundary layer is the same as that at the edge of the viscous boundary layer, so the velocity reduction factor $\delta_i^{th}/\delta_i^v$ is no longer required, and (4.3) becomes

$$ \frac{\rho_i T_i U_i s}{d} \sim k \nabla^2 T \approx \frac{k T_i}{c_p} \nabla^2 s \sim \frac{k T_i s}{c_p (\delta_i^{th})^2}, \quad \text{so} \quad U_i = \frac{k d}{\rho_i c_p (\delta_i^{th})^2}, \qquad (4.5) $$

giving the ratio of the boundary layer thicknesses as

$$ \frac{\delta_i^v}{\delta_i^{th}} \sim Pr^{1/2} \qquad (4.6) $$

in the low-Prandtl-number case.
<figure><img src="image_5.png" /><figcaption>Figure 2. (a) Thermal and viscous boundary layers in the case $Pr > 1$. The thermal diffusion is smaller, so the thermal boundary layer is nested inside the viscous boundary layer. (b) The case $Pr < 1$, where the viscous boundary layer is nested inside the thermal boundary layer.</figcaption></figure>

## 5. The boundary layer ratio problem

In Boussinesq convection there is a symmetry about the midplane, which means that the top and bottom boundary layers have the same thickness and structure, and the temperature of the bulk interior is midway between that of the boundaries. In anelastic convection this symmetry no longer holds, so the top and bottom boundary layers can be very different, and the entropy of the bulk interior is significantly different from $\Delta S/2$. This raises the question of how the ratios of the thicknesses of the top and bottom boundary layers, the ratio of the bulk horizontal velocities just outside the boundary layers and the ratio of the entropy jumps across the layers are actually determined. We assume the incompressible boundary layer case holds throughout this section.

We write the ratio of the entropy jumps across the boundary layers as

$$r_s = \frac{\Delta S_T}{\Delta S_B}, \quad \text{so} \quad \Delta S_B = \frac{\Delta S}{1 + r_s} \quad \text{and the entropy in the bulk is} \quad \frac{r_s}{1 + r_s} \Delta S, \qquad (5.1)$$

and the ratio of the anelastic temperature jumps across the boundary layers is

$$r_T = \frac{\Delta T_T}{\Delta T_B}. \qquad (5.2)$$

We define the ratio of the velocities at the edge of the boundary layers as

$$r_u = \frac{U_T}{U_B}. \qquad (5.3)$$

The last important ratio is the ratio of the thicknesses of the boundary layers. In general the viscous and thermal boundary layers will be of different thickness, but here we start with the thermal boundary layers which have thicknesses at the top and bottom of $\delta_T^{th}$ and $\delta_B^{th}$ with ratio

$$r_{\delta} = \frac{\delta_{T}^{th}}{\delta_{B}^{th}}. \qquad (5.4)$$

We have four unknown ratios, so we need four equations to determine them. Our first equation uses the fact that the heat flux passing through the bottom boundary must
equal the heat flux passing through the top boundary. Since this heat flux is entirely by conduction close to the boundary,

$$-k \left. \frac{dT}{dz} \right|_i \sim k \frac{\Delta T_i}{\delta_i^{th}} \Rightarrow r_\delta = r_T. \qquad (5.5)$$

We can use the balance of advection and diffusion in the boundary layers exploited in § 4 to obtain another equation relating the boundary layer ratios. In § 4 we saw that the ratio of the thermal boundary layer thicknesses was the same as the ratio of the viscous boundary layer thicknesses,

$$\frac{\delta_T^{th}}{\delta_B^{th}} = \frac{\delta_T^\nu}{\delta_B^\nu} = r_\delta, \qquad (5.6)$$

so we use the viscous boundary balance equation (4.2) to estimate

$$\frac{\rho_B U_B^2}{d} \sim \frac{\mu U_B}{\delta_B^{\nu 2}}, \quad \frac{\rho_T U_T^2}{d} \sim \frac{\mu U_T}{\delta_T^{\nu 2}} \quad \Rightarrow \quad r_u r_\delta^2 \sim \frac{\rho_B}{\rho_T} = \Gamma^m. \qquad (5.7)$$

We now need an equation for the ratio of bulk large-scale flow velocities at the top and bottom of the layer, $r_U$. We consider first the case where the viscous dissipation occurs primarily in the boundary layers, which is likely to be true in numerical simulations with no-slip boundaries. Since the entropy production occurs in the boundary layers and is relatively small in the interior, and since entropy diffusion is small in the bulk interior, the convected entropy flux $\langle \bar{\rho}u_zs \rangle_h$ is approximately constant in the bulk interior. We now multiply the equation of motion (2.1) by $\bar{\rho}u$ and average over the bulk interior, ignoring the small viscous term in the bulk, to get

$$\frac{1}{2} \frac{\partial}{\partial z} (\bar{\rho} \langle u_z u^2 \rangle_h) \approx -\frac{\partial}{\partial z} \langle u_z p \rangle_h + \frac{g}{c_p} \langle \bar{\rho} u_z s \rangle_h. \qquad (5.8)$$

Near the boundary layers, the pressure term $p$ will be approximately constant as shown in (2.23), and since $\langle u_z \rangle_h = 0$, the term involving $\langle u_z p \rangle_h$ will be small there, and we ignore it. Since we expect $\langle \bar{\rho}u_zs \rangle_h$ to be approximately the same just outside the two boundary layers,

$$\frac{\partial}{\partial z} (\bar{\rho} \langle u_z u^2 \rangle_h)|_T \approx \frac{\partial}{\partial z} (\bar{\rho} \langle u_z u^2 \rangle_h)|_B, \qquad (5.9)$$

where here $T$ and $B$ refer to conditions at the top of the bulk and the bottom of the bulk, respectively. In the turbulent bulk interior (unlike the boundary layers), we expect the three components of velocity to have similar magnitudes. It remains to estimate the length scale associated with the $z$-derivative, and this is perhaps the most uncertain part of the analysis. Astrophysical mixing length theory uses the pressure scale height, or a multiple of the pressure scale height, as the mixing length for the vertical length scale. Since we are only interested in the top and bottom ratios here, our results are independent of what multiple of the scale height is chosen. Some support for the mixing length idea can be derived from Kessar et al. (2019), which shows that convective length scales decrease in the bulk towards the top of the layer. We also note that because we are only concerned with ratios, it does not matter whether the pressure scale height or the density scale height
is used. Adopting the pressure scale height,

$$H_T = \frac{d}{(m+1)(\Gamma - 1)}, \quad H_B = \frac{\Gamma d}{(m+1)(\Gamma - 1)} \quad \Rightarrow \quad \frac{H_T}{H_B} = \Gamma^{-1}, \quad (5.10)$$

so that (5.9) gives

$$\frac{\rho_T u_T^3}{H_T} \sim \frac{\rho_B u_B^3}{H_B} \Rightarrow r_u^3 \sim \Gamma^{m-1} \Rightarrow r_u \sim \Gamma^{(m-1)/3}. \quad (5.11)$$

From the incompressible boundary layer equation (2.25a,b) we have $r_s = \Gamma r_T$, so with the other three ratio equations (5.5), (5.7) and (5.11) we have

$$r_s = \Gamma r_T, \quad r_\delta = r_T, \quad r_U r_\delta^2 = \Gamma^m, \quad r_u = \Gamma^{(m-1)/3}, \quad (5.12a-d)$$

with solution

$$r_u = \Gamma^{(m-1)/3}, \quad r_\delta = \Gamma^{(2m+1)/6}, \quad r_s = \Gamma^{(2m+7)/6}, \quad r_T = \Gamma^{(2m+1)/6}. \quad (5.13a-d)$$

In the case where $m = 3/2$, appropriate for ideal monatomic gas,

$$r_u = \Gamma^{1/6}, \quad r_\delta = \Gamma^{2/3}, \quad r_s = \Gamma^{5/3}, \quad r_T = \Gamma^{2/3}. \quad (5.14a-d)$$

Since $\Gamma$ is always greater than 1 and can be large, we see that the entropy in the bulk is much closer to the entropy of the bottom boundary than to the entropy of the top boundary. The top boundary layer is thicker than the bottom boundary layer. The challenge for numerical simulations is to get to a sufficiently high Rayleigh number that the top boundary layer is truly thin, as required by our asymptotic analysis. The ratio of the bulk velocities at the top and bottom is only weakly dependent on $\Gamma$, so again rather large values of $Ra$ are required to establish the asymptotic trend.

Note that in deriving (5.7) we assumed that the horizontal length scale for advection along the boundary layer was $d$, as did Grossmann & Lohse (2000). We found that choosing the vertical length scales in the bulk to be based on the pressure scale height in (5.10) agreed reasonably with the numerics, see § 7 below, so a natural question is whether the horizontal length scale near the top boundary becomes less than $d$ when the layer is strongly stratified. The picture from our numerics is somewhat mixed, and is discussed further in § 7 below. For moderate stratification, the numerics suggest the boundary layers do appear to extend to $d$ at both the top and bottom boundary; including a factor $\Gamma$ in the horizontal length scales in the boundary layers gives poorer agreement with numerical estimates of the boundary layer thickness ratio. However, at the largest values of $\Gamma$ we did find a departure from the (5.14a-d) scalings which could be accounted for by some reduction in the horizontal length scale near the top boundary.

In the case where the viscous dissipation is mainly in the bulk, which happens at low $Pr$ (Grossmann & Lohse 2000), (5.1)-(5.7) still hold, but the argument for (5.11) breaks down because the entropy flux is no longer approximately constant in the bulk because viscous dissipation in the bulk is no longer negligible. This case is discussed in Appendix B.
## 6. The Nusselt number and Reynolds number scaling laws

When the boundary layers are thin, the overall entropy balance reduces to (3.16),

$$ \frac{F^{\text{super}} \Delta \bar{T}}{T_B T_T} \sim \int_0^d \frac{\mu}{\bar{T}} \langle q \rangle_h \, dz. \qquad (6.1) $$

If the dissipation is mainly in laminar boundary layers, the z-derivatives of the horizontal velocities will dominate terms in the expression for $q$, so

$$ q = \frac{\partial u_i}{\partial x_j} \frac{\partial u_i}{\partial x_j} + \frac{1}{3} (\nabla \cdot \mathbf{u})^2 \approx \left(\frac{\partial u_x}{\partial z}\right)^2 + \left(\frac{\partial u_y}{\partial z}\right)^2 \sim \frac{U_i^2}{(\delta_i^\nu)^2}, \qquad (6.2) $$

in these layers. So integrating over the boundary layers of thickness $\delta_i^\nu$ and assuming $\bar{T}$ varies little in these layers,

$$ \frac{F^{\text{super}} \Delta \bar{T}}{T_B T_T} \sim \frac{\mu U_B^2}{T_B \delta_B^\nu} + \frac{\mu U_T^2}{T_T \delta_T^\nu} = \frac{\mu U_B^2}{T_B \delta_B^\nu} \left( 1 + \frac{\Gamma r_u^2}{r_\delta} \right), \qquad (6.3) $$

where we have used the ratios (5.3), (5.4) and (2.8g). Now we can write the superadiabatic flux in terms of the thermal boundary layer thicknesses, using (2.14), (2.26), (2.25a,b), (5.1) and (2.18) to get

$$ F^{\text{super}} = \frac{k \Delta T_B}{\delta_B^{th}} = \frac{k T_B \Delta S}{c_p (1 + r_s) \delta_B^{th}}. \qquad (6.4) $$

Inserting this into (6.3) and using the definition (2.16) for the Rayleigh number, the entropy balance equation can be written as

$$ \frac{k^2 R a \Gamma}{c_p^2 (1 + r_s) d^2 \rho_B^2} \delta_B^\nu \sim U_B^2 \left( 1 + \frac{\Gamma r_u^2}{r_\delta} \right). \qquad (6.5) $$

Now we introduce the Reynolds number near the bottom boundary

$$ Re_B = \frac{\rho_B U_B d}{\mu}, \qquad (6.6) $$

noting that the Reynolds number near the top, $Re_T$, is given by $r_u Re_B$. We also use the definition of the Prandtl number, (4.1), to write (6.5) as

$$ \frac{R a \Gamma}{P r^2 (1 + r_s)} \frac{\delta_B^\nu}{\delta_B^{th}} \sim Re_B^2 \left( 1 + \frac{\Gamma r_u^2}{r_\delta} \right). \qquad (6.7) $$

The entropy balance equation has thus given us a relation between the Reynolds number and the Rayleigh number, which is similar to that of regime I of Grossmann & Lohse (2000) but with additional factors of $\Gamma$. In the high-Prandtl-number limit, applying (4.4) gives

$$ Re_B \sim R a^{1/2} P r^{-5/6} \Gamma^{1/2} (1 + r_s)^{-1/2} \left( 1 + \frac{\Gamma r_u^2}{r_\delta} \right)^{-1/2}, \qquad (6.8) $$

while in the low-Prandtl-number case we get, using (4.6),

$$ Re_B \sim R a^{1/2} P r^{-3/4} \Gamma^{1/2} (1 + r_s)^{-1/2} \left( 1 + \frac{\Gamma r_u^2}{r_\delta} \right)^{-1/2}. \qquad (6.9) $$

We now use the viscous boundary layer balance between advection and diffusion, (4.2), $\rho_B U_B/d \sim \mu/(\delta_B^\nu)^2$, to obtain a balance between $Nu$ and $Re_B$. The boundary layer balance
becomes

$$Re_B \sim \left(\frac{d}{\delta_B^\nu}\right)^2, \qquad (6.10)$$

while

$$Nu = \frac{d}{\delta_B^{th}} \frac{\Gamma \ln \Gamma}{(1 + r_s)(\Gamma - 1)}, \qquad (6.11)$$

using the expression (2.27) for the Nusselt number together with (2.12) and (5.1). Eliminating $d/\delta_B^{th}$ between these two equations yields

$$Re_B^{1/2} = \frac{\delta_B^{th} Nu(\Gamma - 1)(1 + r_s)}{\delta_B^\nu \Gamma \ln \Gamma}. \qquad (6.12)$$

As above, the ratio of the boundary layer thicknesses can be evaluated in terms of the Prandtl number, and (6.12) allows us to eliminate $Re_B$ from (6.7) to obtain the Nusselt number as a function of the Rayleigh number,

$$\frac{Ra\Gamma}{Pr^2(1+r_s)} \left(\frac{\delta_B^\nu}{\delta_B^{th}}\right)^5 \sim \left(\frac{Nu(\Gamma - 1)(1+r_s)}{\Gamma \ln \Gamma}\right)^4 \left(1 + \frac{\Gamma r_u^2}{r_\delta}\right). \qquad (6.13)$$

At large $Pr$, (4.4) gives the Nusselt number vs. Rayleigh number scaling in the form

$$Nu \sim Ra^{1/4} Pr^{-1/12} \frac{\Gamma^{5/4} \ln \Gamma}{\Gamma - 1} (1 + r_s)^{-5/4} \left(1 + \frac{\Gamma r_u^2}{r_\delta}\right)^{-1/4}, \qquad (6.14)$$

while at low $Pr$, (4.6) gives

$$Nu \sim Ra^{1/4} Pr^{1/8} \frac{\Gamma^{5/4} \ln \Gamma}{\Gamma - 1} (1 + r_s)^{-5/4} \left(1 + \frac{\Gamma r_u^2}{r_\delta}\right)^{-1/4}. \qquad (6.15)$$

If we accept the ratio scalings derived in § 5, in the case of a monatomic ideal gas, $\gamma = 5/3$, $m = 3/2$, we can write

$$(1 + r_s)^{-5/4} \left(1 + \frac{\Gamma r_u^2}{r_\delta}\right)^{-1/4} = (1 + \Gamma^{5/3})^{-5/4} (1 + \Gamma^{2/3})^{-1/4}, \qquad (6.16)$$

so we can write the main results (6.14) and (6.15) in terms of $Ra$, $Pr$ and $\Gamma$ only,

$$Nu \sim Ra^{1/4} Pr^{-1/12} \frac{\Gamma^{5/4} \ln \Gamma}{\Gamma - 1} (1 + \Gamma^{5/3})^{-5/4} (1 + \Gamma^{2/3})^{-1/4} \quad \text{for } Pr \gg 1. \qquad (6.17)$$

and

$$Nu \sim Ra^{1/4} Pr^{1/8} \frac{\Gamma^{5/4} \ln \Gamma}{\Gamma - 1} (1 + \Gamma^{5/3})^{-5/4} (1 + \Gamma^{2/3})^{-1/4} \quad \text{for } Pr \ll 1. \qquad (6.18)$$

The equivalent formulae for the Reynolds number are

$$Re_B \sim Ra^{1/2} Pr^{-5/6} \Gamma^{1/2} (1 + \Gamma^{5/3})^{-1/2} (1 + \Gamma^{2/3})^{-1/2} \quad \text{for } Pr \gg 1. \qquad (6.19)$$

and

$$Re_B \sim Ra^{1/2} Pr^{-3/4} \Gamma^{1/2} (1 + \Gamma^{5/3})^{-1/2} (1 + \Gamma^{2/3})^{-1/2} \quad \text{for } Pr \ll 1. \qquad (6.20)$$

As $\Gamma$ becomes large, $Nu$ decreases as $\ln \Gamma \Gamma^{-2}$. So we expect the Nusselt number, the dimensionless measure of the heat transport, to be considerably smaller when the layer
is strongly compressible, i.e. when Γ is large and there are many density scale heights in the layer at a given Ra and Pr. Analogous results for the case where the dissipation is mainly in the bulk rather than in the boundary layers, as can happen at low Pr, are given in Appendix B.

In the Boussinesq limit, Γ is close to unity and $\theta = (\Gamma - 1)/\Gamma$ is small, so $\ln \Gamma / (\Gamma - 1) \to 1$ and $\rho_B \to \rho_T$ and $T_B \to T_T$. Equation (2.1) reduces to the usual Boussinesq equation with $s/c_p$ replaced by $\alpha T$, where for a perfect gas, $\alpha = 1/\bar{T}$ is the coefficient of expansion, consistent with (2.25a,b). The total jump of entropy across the layer, $\Delta S$, is replaced by the total temperature jump $\Delta T = \Delta S/\alpha c_p$, so the Rayleigh number (2.16) reduces to the familiar form $Ra = g\alpha\Delta T d^3/\kappa\nu$, where $\kappa = k/\bar{\rho}c_p$ and $\nu = \mu/\bar{\rho}$ are the thermal diffusivity and kinematic viscosity, respectively. These are both constant in the Boussinesq limit, and the ratios $r_u$, $r_\delta$ and $r_s$ all go to unity; see § 5. Our scaling laws (6.8), (6.9), (6.14) and (6.15) all reduce to those of regimes $I_u$ and $I_l$ of Grossmann & Lohse (2000). Grossmann & Lohse also give suggested prefactors for their scaling laws in table 2 of their paper, and since our anelastic scaling laws reduce to theirs in the Boussinesq limit, our prefactors should in theory be consistent with theirs. In practice, the prefactors depend on the aspect ratio of the experiments (or numerical experiments) used to determine them (see, e.g. Chong et al. 2018). For the case of the high-Prandtl-number regime $I_u$, their values were $Nu \approx 0.33Ra^{1/4}Pr^{-1/12}$ and $Re \approx 0.039Ra^{1/2}Pr^{-5/6}$, so (6.14) becomes

$$Nu \approx C_{Nu} Ra^{1/4} Pr^{-1/12} \frac{\Gamma^{5/4} \ln \Gamma}{\Gamma - 1} (1 + r_s)^{-5/4} \left(1 + \frac{\Gamma r_u^2}{r_\delta}\right)^{-1/4}, \quad C_{Nu} = 0.93. \quad (6.21a,b)$$

In the low Pr limit where (6.15) applies, regime $I_l$ of Grossmann & Lohse (2000), they suggest a prefactor corresponding to $C_{Nu} = 0.76$ rather than 0.93. For the Reynolds number, (6.8) becomes

$$Re_B \approx C_{Re} Ra^{1/2} Pr^{-5/6} \Gamma^{1/2} (1 + r_s)^{-1/2} \left(1 + \frac{\Gamma r_u^2}{r_\delta}\right)^{-1/2}, \quad C_{Re} = 0.078. \quad (6.22a,b)$$

There is some uncertainty about the prefactor $C_{Re}$, discussed in § 7 below. Prefactors in the case $I_l$ and in the case where dissipation is mainly in the bulk, their case $II_l$, (see Appendix B) can also be found.

## 7. The numerical results and discussion

We have tested the theoretical predictions of our asymptotic theory using numerical simulations of high-Rayleigh-number plane layer anelastic convection. The numerical code differs from the theory in one respect, as it uses entropy diffusion $k_s$ rather than temperature diffusion in the heat equation, so

$$\bar{\rho} \bar{T} \left[ \frac{\partial s}{\partial t} + \mathbf{u} \cdot \nabla s \right] = k_s \nabla \cdot \bar{T} \nabla s + \mu \left[ q + \frac{\partial u_i}{\partial x_j} \frac{\partial u_j}{\partial x_i} - (\nabla \cdot \mathbf{u})^2 \right], \quad (7.1)$$

where $k_s$ is constant, replaces (2.3). This simplifies the code because entropy is the only anelastic thermodynamic variable computed, and it can be justified in circumstances where turbulent diffusion dominates laminar diffusion (Braginsky & Roberts 1995). Constant entropy boundary conditions were used in the code. The energy balance equation (3.5) is
only changed by replacing $-k \frac{d\langle T \rangle_h}{dz}$ by $-k_s \frac{\tilde{T}}{\tilde{T}} \frac{d\langle s \rangle_h}{dz}$. In the entropy balance equation, entropy diffusion $S_{diff}$ changes to

$$S_{diff} = -k_s \left. \frac{d}{dz} \langle s \rangle_h \right|_{z=0} + k_s \left. \frac{d}{dz} \langle s \rangle_h \right|_z - \int_0^z \frac{k_s \Delta \tilde{T}}{\tilde{T} d} \frac{d}{dz} \langle s \rangle_h dz. \quad (7.2)$$

Just as in the temperature diffusion case, the last two terms are negligible when the entropy boundary layers are thin, except when $z$ is in the boundary layer, where the second term is significant. In the case when the viscous dissipation is primarily in the boundary layers, the argument leading to (5.9) still holds, so the ratios satisfy (5.13a-d) in the entropy diffusion case as well as in the temperature diffusion case. It is therefore reasonable to compare with the numerical results for the entropy diffusion case in the expectation that they will be reasonably close to the temperature diffusion case.

The numerical code is described in Kessar et al. (2019), though for that paper, stress-free boundary conditions were applied, whereas no-slip boundaries where imposed in the runs described here. The unit of length is taken as $d$, the unit of time is $\rho_B d^2 c_p / k$, the thermal diffusion time at the bottom of the layer. The velocities are in units of $k / \rho_B d c_p$, so they correspond to local Péclet numbers, where $Pe = RePr$.

All the runs have polytropic index $m = 3/2$. The code assumes periodic boundary conditions in the horizontal $x$ and $y$ directions, with aspect ratio 2, that is, the period in the horizontal directions is $2d$.

Table 1 gives the parameters used in the eleven runs, which span a range of Prandtl numbers and are at Rayleigh numbers which are as large as is numerically feasible, bearing in mind the need to resolve the small-scale structures that develop. The last three runs are for Boussinesq cases, $\Gamma = 1$, for comparison with the anelastic cases. The density stratification measured by $\Gamma$ varies over a moderate range only, because for the modest values of $Ra$ that are numerically accessible, large $\Gamma$ leads to a top boundary layer which is no longer thin, so our theory will no longer be valid. In figure 3 the entropy profiles (in units of $c_p$) and the horizontal velocity profiles are shown for the three runs A1, B1 and C1, and the profiles for A4 and A5 are shown in figure 4. The entropy profiles are constructed by horizontal averaging and time averaging the vertical profiles. From (2.12) the entropy difference between the boundaries is $\Gamma \ln \Gamma / (\Gamma - 1)$ and the constant is chosen so that it is zero at the bottom boundary $z = 0$.

From figure 3 it is immediately apparent that the entropy is indeed rather constant in the well-mixed bulk interior, consistent with a key assumption of the theory. It is also clear that the jump in entropy across the top boundary layer is greater than that across the bottom boundary layer, and that the top entropy boundary layer is thicker than the bottom entropy boundary layer. This is consistent with the boundary layer ratios found in § 5. The velocity profiles have local maxima near the boundaries, which gives a convenient definition for the top and bottom Reynolds numbers, $Re_B$ and $Re_T$, after converting Péclet numbers to Reynolds numbers using $RePr = Pe$. We note that there is no great difference between the top and bottom horizontal velocities, consistent with the weak scaling with $\Gamma$ found in (5.11). This result is slightly surprising, because astrophysical mixing length theory predicts faster velocities where the fluid is less dense, but in our problem the boundaries play an important role. The low-Prandtl-number case, figures 3(e) and 3(f), has a slightly different entropy boundary layer profile from those of the $Pr = 1$ and $Pr = 10$ cases, with a more gradual matching on to the uniform entropy interior, which is particularly noticeable in the upper boundary layer. This suggests it may be necessary to go to higher $Ra$ at low $Pr$ before accurate agreement with a theory that assumes thin boundary layers can be obtained. The cases shown in figure 4, together with figures 3(a) and 3(b), form a
Fully developed anelastic convection with no-slip boundaries

<table><caption>Table 1. Data from the numerical runs all corresponding to *m* = 3/2 polytropes. The first four rows are the input parameters. Here *r*<sub>δ</sub>, *r*<sub>s</sub> and *r*<sub>u</sub> are the measured boundary layer ratios for each run. The velocities *U*<sub>T</sub> and *U*<sub>B</sub> are the local maxima at the edge of the boundary layers, measured in velocity units of *k*/dρ<sub>Bc</sub><sub>p</sub>, i.e they are Péclet numbers based on the diffusivity at the base of the layer. The theoretical predictions for the boundary layer ratios are given in the next three rows; see (5.14a–d). The Nu-theory entries are based on (6.21a,b) with the prefactors C<sub>Nu</sub> as given in the text, and the boundary ratios come from (5.14a–d). The Nu-nblr entries also use (6.21a,b) with the same prefactors, but instead of using (5.14a–d), the numerical boundary layer ratios (nblr) above are used. The Pe<sub>T</sub>-theory and Pe<sub>B</sub>-theory entries come from (6.8) and (6.9). The prefactors used are not those of Grossmann & Lohse (2000), see (6.22a,b), but those given in the text. Again, (5.14a–d) is used to determine the boundary layer ratios. The Pe<sub>T</sub>-nblr and Pe<sub>B</sub>-nblr entries use the numerical boundary layer ratios.</caption><thead><tr><th>Run</th><th>A1<br>B1</th><th>A2<br>B2</th><th>A3<br>C1</th><th>A4<br>D1</th><th>A5<br>D2</th><th>D3</th></tr></thead><tbody><tr><td>Ra</td><td>10<sup>6</sup></td><td>3 × 10<sup>6</sup></td><td>10<sup>7</sup></td><td>3 × 10<sup>6</sup></td><td>6 × 10<sup>6</sup></td><td></td></tr><tr><td></td><td>10<sup>6</sup></td><td>3 × 10<sup>6</sup></td><td>10<sup>6</sup></td><td>10<sup>6</sup></td><td>10<sup>6</sup></td><td>10<sup>6</sup></td></tr><tr><td>Pr</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td></td></tr><tr><td></td><td>10</td><td>10</td><td>0.25</td><td>1</td><td>10</td><td>0.25</td></tr><tr><td>Γ</td><td>1.9438</td><td>1.9438</td><td>1.9438</td><td>2.924</td><td>4.6416</td><td></td></tr><tr><td></td><td>1.9438</td><td>4.6416</td><td>1.9438</td><td>1</td><td>1</td><td>1</td></tr><tr><td>ρ<sub>b</sub>/ρ<sub>T</sub></td><td>2.71</td><td>2.71</td><td>2.71</td><td>5</td><td>10</td><td></td></tr><tr><td></td><td>2.71</td><td>10</td><td>2.71</td><td>1</td><td>1</td><td>1</td></tr><tr><td>r<sub>δ</sub></td><td>1.65 ± 0.01</td><td>1.70 ± 0.01</td><td>1.69 ± 0.01</td><td>2.05 ± 0.01</td><td>2.27 ± 0.02</td><td></td></tr><tr><td></td><td>1.48 ± 0.01</td><td>2.17 ± 0.01</td><td>2.05 ± 0.04</td><td>1</td><td>1</td><td>1</td></tr><tr><td>r<sub>s</sub></td><td>3.21 ± 0.02</td><td>3.29 ± 0.01</td><td>3.31 ± 0.01</td><td>5.99 ± 0.02</td><td>10.62 ± 0.03</td><td></td></tr><tr><td></td><td>2.91 ± 0.02</td><td>10.0 ± 0.02</td><td>4.04 ± 0.04</td><td>1</td><td>1</td><td>1</td></tr><tr><td>r<sub>u</sub></td><td>0.92 ± 0.01</td><td>0.92 ± 0.01</td><td>0.93 ± 0.01</td><td>0.96 ± 0.01</td><td>1.06 ± 0.02</td><td></td></tr><tr><td></td><td>1.17 ± 0.01</td><td>1.38 ± 0.01</td><td>0.82 ± 0.02</td><td>1</td><td>1</td><td>1</td></tr><tr><td>Nu</td><td>5.90 ± 0.02</td><td>7.75 ± 0.03</td><td>10.69 ± 0.03</td><td>5.35 ± 0.02</td><td>4.26 ± 0.02</td><td></td></tr><tr><td></td><td>6.18 ± 0.02</td><td>4.09 ± 0.02</td><td>5.00 ± 0.04</td><td>8.78 ± 0.03</td><td>8.76 ± 0.04</td><td>8.17 ± 0.04</td></tr><tr><td>U<sub>T</sub></td><td>186 ± 3</td><td>321 ± 3</td><td>570 ± 3</td><td>298 ± 3</td><td>368 ± 3</td><td></td></tr><tr><td></td><td>287 ± 2</td><td>410 ± 3</td><td>125 ± 4</td><td>200 ± 3</td><td>260 ± 5</td><td>149 ± 3</td></tr><tr><td>U<sub>B</sub></td><td>202 ± 3</td><td>348 ± 3</td><td>613 ± 3</td><td>311 ± 3</td><td>348 ± 3</td><td></td></tr><tr><td></td><td>244 ± 2</td><td>298 ± 3</td><td>154 ± 4</td><td>200 ± 3</td><td>260 ± 5</td><td>149 ± 3</td></tr><tr><td>Γ<sup>2/3</sup></td><td>1.557</td><td>1.557</td><td>1.557</td><td>2.045</td><td>2.783</td><td></td></tr><tr><td></td><td>1.557</td><td>2.783</td><td>1.557</td><td>1</td><td>1</td><td>1</td></tr><tr><td>Γ<sup>5/3</sup></td><td>3.028</td><td>3.028</td><td>3.028</td><td>5.979</td><td>12.915</td><td></td></tr><tr><td></td><td>3.028</td><td>12.915</td><td>3.028</td><td>1</td><td>1</td><td>1</td></tr><tr><td>Γ<sup>1/6</sup></td><td>1.117</td><td>1.117</td><td>1.117</td><td>1.196</td><td>1.292</td><td></td></tr><tr><td></td><td>1.117</td><td>1.292</td><td>1.117</td><td>1</td><td>1</td><td>1</td></tr><tr><td>Nu-theory</td><td>5.56</td><td>7.32</td><td>9.89</td><td>4.65</td><td>2.98</td><td></td></tr><tr><td></td><td>5.55</td><td>2.50</td><td>5.18</td><td>8.78</td><td>8.76</td><td>8.17</td></tr><tr><td>Nu-nblr</td><td>5.60</td><td>7.22</td><td>9.67</td><td>4.97</td><td>3.86</td><td></td></tr><tr><td></td><td>5.63</td><td>3.11</td><td>4.37</td><td>8.78</td><td>8.76</td><td>8.17</td></tr><tr><td>Pe<sub>T</sub>-theory</td><td>194</td><td>336</td><td>614</td><td>307</td><td>376</td><td></td></tr><tr><td></td><td>252</td><td>345</td><td>144</td><td>200</td><td>260</td><td>149</td></tr><tr><td>Pe<sub>B</sub>-theory</td><td>174</td><td>301</td><td>549</td><td>257</td><td>291</td><td></td></tr><tr><td></td><td>226</td><td>267</td><td>129</td><td>200</td><td>260</td><td>149</td></tr><tr><td>Pe<sub>T</sub>-nblr</td><td>177</td><td>306</td><td>559</td><td>283</td><td>361</td><td></td></tr><tr><td></td><td>256</td><td>358</td><td>118</td><td>200</td><td>260</td><td>149</td></tr><tr><td>Pe<sub>B</sub>-nblr</td><td>192</td><td>332</td><td>601</td><td>295</td><td>341</td><td></td></tr><tr><td></td><td>219</td><td>260</td><td>144</td><td>200</td><td>260</td><td>149</td></tr></tbody></table>
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 3. Horizontally averaged entropy (in units of $c_p$) and horizontal mean velocity profiles (Péclet number units) from the numerical simulations for $\Gamma = 1.9438$, $Ra = 10^6$, runs A1, B1 and C1. (a) Entropy profile at $Pr = 1$. (b) Horizontal velocity profile at $Pr = 1$. (c) Entropy profile at $Pr = 10$. (d) Horizontal velocity profile at $Pr = 10$. (e) Entropy profile at $Pr = 0.25$. (f) Horizontal velocity profile at $Pr = 0.25$.</figcaption></figure>

sequence at constant $Pr = 1$ with increasing $\Gamma$. The most noticeable feature is that at larger $\Gamma$ the entropy of the mixed interior becomes close to that of the bottom boundary, so the boundary layer ratio $r_S$ increases rapidly, consistent with the prediction of (5.14a–d). Also notable is that the velocity ratio of the maximum horizontal velocities, $r_U$, is never far from unity, again consistent with the weak scaling with $\Gamma$ in (5.14a–d). As expected, the boundary layers become thicker at larger $\Gamma$, so that at fixed $Ra$ the thin boundary layer assumption breaks down at large $\Gamma$.

In table 1 we compare various results of the simulations against the theoretical predictions of §§5 and 6. To evaluate the entropy boundary layer thicknesses for our numerical data, we use the definition (2.22), so the gradients $d\langle S \rangle_h/dz$ at $z=0$ and
Fully developed anelastic convection with no-slip boundaries

<figure><img src="image_12.png" /><figcaption>Figure 4. Horizontally averaged entropy $\langle S \rangle_h$ and horizontal mean velocity $U_H$ profiles for (a,b) $\Gamma = 2.924$, $Ra = 3 \times 10^6$, $Pr = 1$, run A4: (c,d) $\Gamma = 4.6416$, $Ra = 6 \times 10^6$, $Pr = 1$, run A5.</figcaption></figure>

$z=1$ were obtained by differentiating a cubic spline representation of the entropy, and the entropy jumps were obtained by averaging the entropy over the well-mixed region, assuming constant entropy there. The ratio of the top and bottom entropy boundary layer thicknesses and the ratio of the entropy jumps are denoted by $r_\delta$ and $r_s$ in table 1. The velocity ratios at the top and bottom are denoted by $r_u$. We used 256 mesh points in the $z$-direction for all runs, with an additional 512 mesh point run for case A3, the highest Ra case, to check that our resolution was sufficient for the table 1 data. A 256 point mesh in the $z$-direction gives about 8 mesh points in the thinner lower boundary layer, which was found to be adequate in the high Ra Rayleigh-Bénard studies of Shishkina et al. (2010). The code is spectral in the horizontal directions, and we used $256 \times 256$ Fourier modes in our runs, giving a resolution comparable to that used by Gastine, Wicht & Aurnou (2015) for similar Ra in spherical shell Rayleigh-Bénard convection, and we did one $512 \times 512$ run to check that this was adequate. The energy spectrum of the $256 \times 256 \times 256$ runs was checked to ensure that the resolution is sufficient to resolve the Kolmogorov length scale. A typical run on the computational cluster used 256 processors and settled down to a nonlinear saturated state after a few hours. However, the typical run time was around a hundred hours, the long run being needed to get reliable time-averaged quantities.

We can compare our numerically obtained ratios with the predicted ratios in (5.14a–d). We see that there is some variation in the numerical results, but they are roughly in agreement with the predicted results. Considering that the top boundary layer is not that thin, as can be seen in figures 3 and 4, these results are as good as can be expected. We also tested how well the equations leading to our boundary layer ratios compare individually with the numerical results. Equation (5.5) expresses the fact that the heat flux is the same
at the top and bottom, and together with the incompressible boundary layer equation (2.25a,b), gives $r_s = \Gamma r_\delta$, which agrees with our numerics rather well, to the 1% level. The viscous boundary layer balance, (5.7), has less accurate agreement with the numerics. All the runs at $\Gamma = 1.9438$ had errors less than 10% (except the low $Pr$ run where, as we saw in figure 3, the boundary layer structure is slightly different), but the runs with density ratio 10, A5 and B2, have $r_\delta$ and $r_s$ large, but not quite as large as predicted by (5.14a–d). The likely explanation is that the horizontal length scale at the top boundary layer is getting smaller than at the bottom boundary, though not by as much as the factor $\Gamma$ predicted for the vertical length scale ratio. There is therefore some doubt as to whether it is correct to have $d$ as the horizontal length scale in the boundary layers when $\Gamma$ is large. Some support for a different horizontal length scale at the two boundaries can be seen in figure 2(b) of Kessar et al. (2019). The convection pattern near the top boundary has broad upwellings surrounded by a network of thin rapidly descending regions, whereas close to the bottom boundary there is no well-defined network, but just a general turbulent behaviour. Even though our boundaries are no-slip rather than stress-free, we also found a network of thin descending regions near the top, but a more turbulent regime near the bottom. The persistent network of descending regions may limit the ‘wind of turbulence’ near the top to a length scale corresponding to that of the upwelling regions rather than that of the layer depth, whereas no such restriction applies at the bottom. The effective $d$ in the top boundary layer balance (5.7) might therefore be somewhat reduced, which would reduce $r_\delta$ consistent with our numerics. Further research is needed to elucidate this issue. The velocity ratio equation (5.11) correctly predicts that the velocity ratio is always close to unity, but is less reliable at predicting whether it is above or below unity. However, the higher density ratio runs do have $r_u > 1$, as predicted by (5.11).

We also evaluated the Nusselt number from the data, using

$$Nu = - \frac{d}{c_p} \left. \frac{d\langle S \rangle_h}{dz} \right|_{z=0} = - \frac{d}{\Gamma c_p} \left. \frac{d\langle S \rangle_h}{dz} \right|_{z=d}, \quad (7.3)$$

again determining the gradients from our spline representation. When the run has been integrated long enough, initial transients in the numerical run are eliminated and these two estimates of the Nusselt number become close, and the average value is used in table 1. Because the flow is turbulent, the Nusselt number fluctuates continuously at about the 10% level, so a long time-average is used. The finite length of the run means there is a small uncertainty due to the fluctuations not exactly cancelling, which we estimate as error bars in table 1.

In table 1 we also give the value of the Nusselt number calculated by our theory. We used the Boussinesq runs D1, D2 and D3 to determine the prefactors $C_{Nu}$ appropriate for each Prandtl number used. This gives $C_{Nu} = 0.949$ for $Pr = 10$, $C_{Nu} = 0.785$ for $Pr = 1$ and $C_{Nu} = 0.869$ for $Pr = 0.25$. Note that in table 1 this means that for the Boussinesq runs D1, D2, and D3, the *Nu*-theory entries are the same as the actual *Nu* entries by construction. None of these prefactors is very far from the values suggested by Grossmann & Lohse (2000), who had $C_{Nu} = 0.93$ at large $Pr$ and $C_{Nu} = 0.76$ at small $Pr$, suggesting that the differences in aspect ratio and geometry only make relatively small changes to the Nusselt number. For consistency, we use these same prefactors in all runs. Since our main interest is in the compressible cases, we have not explored why the prefactor for the $Pr = 1$ case is slightly lower than the other values of $C_{Nu}$.

We first consider the Nusselt number for the cases where the density ratio is only 2.71, runs A1, A2, A3, B1 and C1. We note that all the predicted values are not too far off the numerical values, though the predicted values are generally a little lower than
the actual numerical values. Using the numerically calculated boundary layer ratios in formula (6.21a,b) rather than the theoretically predicted ones gives the result *Nu*-nblr in table 1. These are only available after the simulation is run, but they are helpful for testing whether small errors are due to slightly inaccurate boundary layer ratios, or whether the formula (6.21a,b) is inaccurate. For the density ratio 2.71 runs with *Pr* ≥ 1, the boundary layer ratios were close to the predicted values, so not surprisingly, the *Nu*-nblr values are not significantly more accurate than the *Nu*-theory values. We conclude that the under-prediction of the Nusselt number in these cases, which is less than about the 10% level, is due to the viscous dissipation not being completely in the boundary layers, as required by the theory. We believe that at higher *Ra*, where the dissipation progressively goes into the boundary layers, and using longer runs to average out the fluctuations, the small discrepancy will disappear. Using runs A1, A2 and A3, where only *Ra* varies, we can test the *Ra*¹/₄ power law predicted in (6.21a,b). The least squares fit to a straight line in log(*Nu*) vs log(*Ra*) space has a slope of 0.257, rather close to the predicted slope.

We now look at the larger Γ cases, runs A4, A5 and B2, corresponding to the more
compressible cases. In the run A4, at density ratio 5, the predicted and numerical Nusselt
numbers are reasonably close, but in the most extreme cases of density ratio 10 the
predicted *Nu* is only 61% of the numerical value for run B2, and in the run A5 the predicted
*Nu* is 71% of the numerical *Nu*. Part of this discrepancy is down to the boundary layer
ratios, which become very large at high Γ, and so any small inaccuracy can affect the
Nusselt number significantly. If we use the numerical boundary layer ratios in (6.21a,b)
rather than the theoretical ones for run B2, the predicted *Nu*-nblr rises to 3.12, but even
this is only 76% of the numerical value, and A5 similarly improves but still is too low.
We again conclude that in runs B2 and A5 the assumption that the dissipation occurs
dominantly in the boundary layers is suspect (particularly near the top boundary) and that
higher *Ra* is needed before it becomes robustly valid.

We now consider the Reynolds number formula, (6.22a,b), though it is convenient
to express this in terms of Péclet number Pe = RePr. If the table 1 parameter values
are inserted into (6.22a,b) with the value of CRe quoted there, the values of the Péclet
number are consistently a factor of about five too small compared with the numerical
values of UT and UB in table 1. There are a number of reasons for this, but the two most
important are (i) the power law dependence of the Péclet number with Rayleigh number in
Boussinesq convection is slightly less than the predicted 0.5 (Grossmann & Lohse 2002),
and (ii) our runs are for aspect ratio 2, whereas experiments, and the numerical simulations
that simulate them (e.g. Silano, Sreenivasan & Verzicco 2010), use aspect ratios of 1 or
less. The prefactor in (6.22a,b) is based on experiments at large Ra which mostly used
aspect ratios less than unity. The experiments of Qiu & Tong (2001), see also figure 1 of
Grossmann & Lohse (2000), using water (with Pr = 5.5) in a cylinder of aspect ratio unity
found Re = 0.085Ra0.455. At the run A1 parameters this formula gives Pe = 251 consistent
with a prefactor of CRe = 0.38 in (6.22a,b), much larger than the Grossmann & Lohse
(2000) value. They found very similar Re prefactors for both high and low Pr. If we adopt
the same procedure as we did to get the Nusselt number prefactors, and normalise using
the Boussinesq runs D1, D2 and D3 we obtain CRe = 0.354 for Pr = 10, CRe = 0.400
for Pr = 1 and CRe = 0.421 for Pr = 0.25. Reassuringly, these are all quite close to the
Qiu & Tong (2001) value of CRe = 0.38. We therefore use these three values of CRe at
the appropriate Prandtl number in all our theory calculations. With these prefactors, the
numerical results for UT and UB agree reasonably well with the predicted PeT-theory and
PeB-theory results. The results have some scatter, which seems to reflect the scatter in our
computed ru. If we use the computed boundary layer ratios, rather than the asymptotically
predicted ratios, there is less scatter in the comparison between computed and theoretical
C.A. Jones, K.A. Mizerski and M. Kessar

Péclet numbers, though the theoretical Péclet numbers are generally a few percent lower than the computed Péclet numbers. Given that the boundary layers are not very thin, these small discrepancies are not unexpected, and overall the predicted Reynolds numbers are in reasonable agreement with those of our § 6 asymptotic theory.

## 8. Conclusions

The scaling laws for heat flux and Reynolds number at high-Rayleigh-number convection have been derived from the energy balance and entropy balance equations derived in § 3. These scaling laws are formulated in terms of the Rayleigh number, the Prandtl number and the temperature ratio $\Gamma$ which measures the strength of the stratification. In the Boussinesq limit, $\Gamma \to 1$, they reduce to the scaling laws of Grossmann & Lohse (2000). The existence of the well-mixed entropy state, with the entropy changes being mainly confined to thin boundary layers, makes it possible to estimate the terms in the entropy balance equation, so allowing Nusselt number and Reynolds number relationships to be established. The cases treated are those where the viscous dissipation occurs in the boundary layers, the cases labelled as $I_u$ and $I_l$ by Grossmann & Lohse (2000), the subscripts referring to the high- and low-Prandtl-number regimes, and the cases where the viscous dissipation is primarily in the bulk, the cases $I_u$ and $I_l$. A limitation of the theory is that both the entropy boundary layers do have to be thin for the theory to be valid. For the top boundary layer to be thin when the stratification is strong, the Rayleigh number has to be very large, which is numerically difficult, so the range of $\Gamma$ which can be tested both numerically and asymptotically is quite limited. This condition that the top boundary layer is thin is equivalent to the condition that the boundary layers are incompressible, so that a rather simple relationship holds between temperature and entropy within the boundary layers. The more difficult case where the boundary layers are compressible has not yet been solved in closed form, but it is likely to be significantly different from our solutions.

A feature of this high-Rayleigh-number anelastic problem is that the top and bottom boundary layers have a different structure, so to determine the scaling laws, boundary layer ratios for the top and bottom boundary layers have to be established. The three key ratios are those for the boundary layer widths, the boundary layer entropy jumps and the horizontal velocities just outside the boundary layers. In § 4 we found that although the top and bottom boundary layers are different, the ratios of the viscous and thermal (entropy) boundary layer thicknesses at each boundary depend only on the Prandtl number. In § 5 we proposed formulae based on a simple physical picture for these ratios. The higher density at the bottom means that the viscous shear must be larger there (so a thinner boundary layer) to bring the fluid to rest at the boundary. The boundary layer ratios will be sensitive to the assumptions made about the viscosity, e.g. whether there is constant dynamic viscosity (the assumption made here) or constant kinematic viscosity. The ratios will also change if the viscosity is temperature dependent, as the top and bottom boundaries have significantly different temperatures if $\Gamma > 1$.

We have performed some numerical simulations to test these proposed boundary layer ratios, and within the constraints imposed by the numerics, namely not very high $Ra$, we find broad agreement between the theory and the numerics. Another important assumption for Grossmann-Lohse theory to be valid is the existence of a wind of turbulence. Our numerics suggest that this feature persists in our simulations. There is, however, still some uncertainty about whether the horizontal length scale of that wind, which controls the boundary layers, remains at the vertical length scale $d$ as the stratification $\Gamma$ increases, or whether it becomes smaller at the top boundary than the bottom boundary at large $\Gamma$.
Fully developed anelastic convection with no-slip boundaries

We have also tested the theoretically derived Nusselt number and Reynolds number relationships against the numerics, in the case where the viscous dissipation is mainly in the boundary layers, the only numerically accessible case. Using the prefactors determined in the Boussinesq case, which are the only free parameters in the theory, the Nusselt numbers obtained are in reasonable agreement with the theory, again noting the numerical limitations preventing accurate agreement. A problem was encountered when comparing with the theoretical Reynolds numbers, in that the theory using the original Grossmann–Lohse prefactors gave smaller *Re* than did the numerics. However, the disagreement seems to be due more to issues with the Boussinesq problem rather than to its extension to the anelastic case, in particular to the difficulty of establishing a single prefactor over a huge range of *Ra* and to the dependence of *Re* on the aspect ratio. When the prefactors were determined by normalizing on our Boussinesq runs, the issue was resolved.

We have focused here on the case of no-slip boundaries, as this seems the simplest case in which scaling laws can be derived from first principles without introducing arbitrary constants into the formulae. There are, however, many problems of great astrophysical interest that could be similarly addressed: the case of stress-free boundaries is thought to be particularly relevant to stellar convection zones. Even within the context of our simplified no-slip problem, the case of compressible boundary layers would be of interest. We found it most convenient to consider fixed entropy boundary conditions, but other boundary conditions, such as fixed temperature or fixed flux are of interest too. Another issue that could be explored are the differences between temperature diffusion and entropy diffusion cases. In our particular problem, with incompressible boundary layers, the differences appear to be quite minor, but this is not necessarily the case if more challenging cases are addressed. Given the growing importance of the anelastic approximation in exploring a very wide range of exciting astrophysical problems, a firmer understanding of the fundamental behaviour of high-Rayleigh-number anelastic convection would be very valuable.

**Acknowledgements.** We thank three anonymous reviewers for their helpful comments. The computational work was performed on the ARC clusters, part of the high performance computing facilities at the University of Leeds, and on the COSMA Data Centre system at Durham University, operated on behalf of the STFC DiRAC HPC.

**Funding.** This work was partially funded by the STFC grant ST/S00047X/1 held at the University of Leeds. The partial funding from the subvention of the Ministry of Science and Higher Education in Poland as a part of the statutory activity and the support of the National Science Centre of Poland (grant no. 2017/26/E/ST3/00554) is gratefully acknowledged.

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

**Author ORCIDs.**

Chris A. Jones https://orcid.org/0000-0002-1105-3920;

Krzysztof A. Mizerski https://orcid.org/0000-0003-0106-675X.

**Appendix A. Form of the anelastic temperature perturbation**

Taking the horizontal average of the anelastic continuity equation (2.2), and using the $u_z = 0$ boundary conditions gives

$$ \langle u_z \rangle_h = 0. \qquad (A1) $$
Using (2.1) and (2.13), the z-component of the anelastic equation of motion can be written as

$$ \bar{\rho} \frac{\partial u_z}{\partial t} + \nabla \cdot (\bar{\rho} u_z \mathbf{u}) + \nabla p = -g\rho + \mu \left( \nabla^2 u_z + \frac{1}{3} \frac{\partial}{\partial z} \nabla \cdot \mathbf{u} \right). \quad (A2) $$

Taking the horizontal average of (A2) we see that, using (A1), the viscous term vanishes to leave

$$ \frac{d}{dz} \langle \bar{\rho} u_z^2 \rangle_h + \frac{d \langle p \rangle_h}{dz} = -g \langle \rho \rangle_h = -g \bar{\rho} \left( \frac{\langle p \rangle_h}{\bar{p}} - \frac{\langle T \rangle_h}{\bar{T}} \right) \quad (A3) $$

using (2.4a). In the bulk, entropy is well mixed, so it is constant there, so differentiating (2.4b) and using (2.4a),

$$ R \frac{d}{dz} \left( \frac{\langle p \rangle_h}{\bar{p}} \right) = c_p \frac{d}{dz} \left( \frac{\langle T \rangle_h}{\bar{T}} \right) \quad (A4) $$

in the bulk. Using the adiabatic reference state hydrostatic and perfect gas equations, with (2.8a-c), this can be written as

$$ \frac{d\langle p \rangle_h}{dz} = c_p \bar{\rho} \frac{d\langle T \rangle_h}{dz} - \frac{g\bar{\rho}}{\bar{p}} \langle p \rangle_h + \frac{g\bar{\rho}}{\bar{T}} \langle T \rangle_h, \quad (A5) $$

and on substituting this into (A3) we obtain

$$ \frac{d\langle T \rangle_h}{dz} = - \frac{1}{c_p \bar{\rho}} \frac{d}{dz} \langle \bar{\rho} u_z^2 \rangle_h, \quad (A6) $$

which is valid in the bulk. Integrating this across the bulk from $z = \delta_B^{th}$ to $z = d - \delta_T^{th}$, and assuming $u_z$ is negligible close to the boundaries,

$$ \langle T_{\text{bulk}}(d - \delta_T^{th}) \rangle_h - \langle T_{\text{bulk}}(\delta_B^{th}) \rangle_h = \int_{\delta_B^{th}}^{d-\delta_T^{th}} \langle \bar{\rho} u_z^2 \rangle_h \frac{d}{dz} \left( \frac{1}{c_p \bar{\rho}} \right) dz = \Delta T_{\text{vel}} > 0, \quad (A7) $$

since $\bar{\rho}$ is monotonic decreasing with $z$. This establishes that in the bulk the gradient $d\langle T \rangle_h/dz$ is positive on average, corresponding to a subadiabatic horizontally averaged temperature gradient. We denote this jump in $T$ across the bulk by $\Delta T_{vel}$ because it is physically connected to the pressure changes induced by the fluid velocity. A natural question is how large $\Delta T_{vel}$ is compared with the jumps in $\langle T \rangle_h$ across the boundary layers, $\Delta T_B$ and $\Delta T_T$? Formally they are both of the same order of magnitude in the anelastic approximation, but $\Delta T_{vel}$ will be small if we are close to Boussinesq or if the Rayleigh number is small. Numerical evidence is sparse, but figure 4 from Verhoeven et al. (2015) suggests that for their parameters, $Ra = 10^6$, $\rho_B/\rho_T = 2.72$ and $Pr = 0.7$, their $\Delta T_{vel}$ was small.

### A.1. Positivity of the temperature offsets

We now consider the temperature offsets at the bottom and top boundaries, $\langle T \rangle_{h,T}$ and $\langle T \rangle_{h,B}$. Without numerical simulations, we cannot determine their magnitude, but we can show that they must both be positive, a useful check on future simulations.
Fully developed anelastic convection with no-slip boundaries

By examining the sum of the temperature jumps across the layer in figure 1(b) we can
see that

$$
\langle T \rangle_{h,T} - \langle T \rangle_{h,B} + \Delta T_B + \Delta T_T = \Delta T_{vel.} \quad (A8)
$$

Using the incompressible boundary layer forms for the temperature jumps across the
boundary layers, (2.25a,b), and the formulae for the ratios of these jumps,

$$
r_s = \frac{\Delta S_T}{\Delta S_B}, \quad r_T = \frac{\Delta T_T}{\Delta T_B}, \quad r_s = \Gamma r_T, \qquad (A9a-c)
$$

(A8) becomes

$$
\langle T \rangle_{h,T} - \langle T \rangle_{h,B} + \frac{T_B \Delta S}{c_p} \frac{(1+r_T)}{(1+\Gamma r_T)} = \Delta T_{vel}. \quad (A10)
$$

A second equation for the temperature offsets can be derived from the boundary
conditions. From (2.4a) and (2.4b) we can deduce that

$$
s = \frac{c_p T}{\bar{T}} - \frac{p}{\bar{\rho} \bar{T}}. \tag{A11}
$$

At $z=0$ and $z=d$ this gives

$$
\Delta S = \frac{c_p \langle T \rangle_{h,B}}{T_B} - \frac{\langle p \rangle_{h,B}}{\rho_B T_B}, \quad 0 = \frac{c_p \langle T \rangle_{h,T}}{T_T} - \frac{\langle p \rangle_{h,T}}{\rho_T T_T}. \quad (\text{A12a,b})
$$

We now use the mass conservation equation (2.28) to set the pressure perturbations on the
top and bottom boundary equal, giving

$$
\frac{T_B \Delta S}{c_p} = \langle T \rangle_{h,B} - \langle T \rangle_{h,T} \frac{\rho_T}{\rho_B}. \quad (A13)
$$

Equations (A10) and (A13) are two equations for the temperature offsets $\langle T \rangle_{h,B}$ and $\langle T \rangle_{h,T}$,
and using $\rho_B/\rho_T = \Gamma^m$ the solutions are

$$
\langle T \rangle_{h,B} = \frac{\Delta T_{vel}}{\Gamma^m - 1} + \frac{T_B \Delta S}{c_p} \left\{ \frac{\Gamma^m (1 + \Gamma r_T) - (1 + r_T)}{(\Gamma^m - 1)(1 + \Gamma r_T)} \right\}, \quad (A14)
$$

$$
\langle T \rangle_{h,T} = \frac{\Gamma^m \Delta T_{vel}}{\Gamma^m - 1} + \frac{T_B \Delta S}{c_p} \left\{ \frac{(\Gamma - 1) \Gamma^m r_T}{(\Gamma^m - 1)(1 + \Gamma r_T)} \right\}. \quad (A15)
$$

Since $\Gamma > 1$ and $\Delta T_{vel} > 0$, it follows that both quantities are positive whatever $\Delta T_{vel}$ is.
It is not possible to decide which offset is larger without having more information about
$\Delta T_{vel}$, but these results confirm that figure 1(b) is a plausible sketch of the temperature
perturbation, and will be helpful in testing numerical simulations.

**Appendix B. The case when the dissipation in the bulk dominates the dissipation in the boundary layers**

Grossmann & Lohse (2000) point out that at low *Pr* and large *Ra* it is possible for the viscous dissipation in the bulk to be larger than the viscous dissipation in the boundary layers. When this occurs, our arguments about the boundary layer ratios in § 5 and the scaling laws in § 6 need revising. We now consider this scenario.
C.A. Jones, K.A. Mizerski and M. Kessar

## B.1. The boundary layer ratios

When the viscous dissipation is mainly in the bulk, (5.1)–(5.7) still hold, but the argument for (5.11) breaks down because the entropy flux is no longer approximately constant in the bulk, since viscous dissipation in the bulk is no longer negligible. We can, however, use the energy flux equation (3.2) because when the viscous dissipation is mainly in the bulk, the work done by buoyancy must balance the viscous dissipation in the bulk, since now the viscous dissipation in the boundary layers is negligible.

So (3.1) becomes

$$ \frac{g}{c_p} \int_{bulk} \langle \bar{\rho} u_z s \rangle_h dz = \mu \int_{bulk} \langle q \rangle_h dz, \qquad (B1) $$

and since thermal diffusion and the last two terms in (3.2) are negligible in the bulk when the boundary layers are thin, it follows that $\langle \bar{\rho} \bar{T} u_z s \rangle_h$ will be approximately the same just outside the two boundary layers at $z = \delta_B^\nu$ and $z = d - \delta_T^\nu$, so

$$ \rho_B T_B \langle u_z s \rangle_h |_{z=\delta_B^\nu} \approx \rho_T T_T \langle u_z s \rangle_h |_{z=d-\delta_T^\nu}. \qquad (B2) $$

Note this is different from the case where the dissipation was mainly in the boundary layers, when $\langle \bar{\rho} u_z s \rangle_h$ is approximately constant. As we did in § 5, we horizontally average the dot product of $u$ and (2.1), and apply it just outside the boundary layers, at $z = \delta_B^\nu$ and $z = d - \delta_T^\nu$. Here we are justified in neglecting the pressure term as we did in § 5, and we also neglect the viscous term. This is not obvious when most of the viscous dissipation is in the bulk, but following Grossmann & Lohse (2000), we envisage a turbulent cascade, where the dissipation at larger scales is dominated by the inertial term. We therefore adopt

$$ \frac{1}{2} \frac{\partial}{\partial z} (\bar{\rho} \langle u_z u^2 \rangle_h) \approx -\frac{\partial}{\partial z} \langle u_z p \rangle_h + \frac{g}{c_p} \langle \bar{\rho} u_z s \rangle_h \approx \frac{g}{c_p} \langle \bar{\rho} u_z s \rangle_h \qquad (B3) $$

at $z = \delta_B^\nu$ and $z = d - \delta_T^\nu$. Then since we expect all velocity components to be of similar magnitude in the bulk, using (B2),

$$ \rho_B T_B \frac{U_B^3}{H_B} \approx \rho_T T_T \frac{U_T^3}{H_T} \Rightarrow r_u \sim \Gamma^{m/3}, \qquad (B4) $$

where the pressure scale heights are defined in (5.10). This result differs from (5.11), where the dissipation is in the boundary layers, so that now the horizontal velocity at the top is expected to be considerably faster than the velocity at the bottom, whereas (5.13a–d) predicts only a weak dependence on $\Gamma$.

## B.2. The scaling laws when dissipation is in the bulk

We still expect thin boundary layers even when the dissipation is mainly in the bulk, so (6.1),

$$ \frac{F^{\text{super}} \Delta \bar{T}}{T_B T_T} \sim \int_0^d \mu \langle q \rangle_h dz, \qquad (B5) $$

still applies, but unlike the boundary layer dissipation case, we do not know how the dissipation is distributed over the interior. We therefore assume that the dissipation in the interior can be written as $\langle q \rangle_h \sim U_H^3/2H$, where $U_H(z)$ is the horizontally averaged horizontal velocity and $H$ is the local pressure scale height. We do not know how $U_H(z)$ is distributed in $z$, but we argued in § B.1 above that $\rho U_H^3$ is approximately the same at the
edge of both boundary layers, so a reasonable assumption for the purposes of estimation is that

$$ \rho U_H^3 \sim \text{const.} \approx \rho_B U_B^3 \approx \rho_T U_T^3. \qquad (B6) $$

The form of $U_H(z)$ from our numerical results suggests this might overestimate the dissipation integrated over the whole layer, but nevertheless we adopt (B6) for the rest of this section. Equation (B5) then becomes

$$ \frac{F^{\text{super}} \Delta \bar{T}}{T_B T_T} \sim \int_0^d \frac{\bar{\rho}_B U_H^3}{2 \bar{T} H} dz = \int_0^d \frac{\rho_B U_B^3 (m+1) \Delta \bar{T}}{2 d \bar{T}^2} dz = \frac{\rho_B U_B^3 (m+1) (\Gamma - 1)}{2 T_b}. \quad (B7) $$

From (2.15), $F^{\text{super}} = \text{Nuk}T_B/d$, and writing $U_B$ in terms of the bottom Reynolds number using (6.6),

$$ \frac{\text{Nuk}\Delta\bar{T}}{dT_T} \sim \frac{\mu^3 \text{Re}_B^3 (m+1)(\Gamma-1)}{2T_b \rho_B^2 d^2}. \qquad (B8) $$

Combining (2.12) and (2.16), we can write the Rayleigh number as

$$ Ra = \frac{c_p^2 \Delta \bar{T} d^2 \rho_B^2 \Gamma \ln \Gamma}{\mu k (\Gamma - 1)}, \qquad (B9) $$

and combining this with (B8) and using the definition of the Prandtl number (4.1) we obtain

$$ \frac{NuRa}{Pr^2} \sim \frac{(m+1) \ln \Gamma}{2} Re_B^3, \qquad (B10) $$

which is the entropy balance equation in the case where the dissipation is mainly in the bulk rather than the boundary layers. We now use the same boundary layer balance equation as before, but since we expect bulk dissipation only to dominate at low $Pr$, we only use (4.5) for the boundary layer ratio, so (6.12) becomes

$$ (Re_B Pr)^{1/2} = \frac{Nu(\Gamma - 1)(1 + r_s)}{\Gamma \ln \Gamma}. \qquad (B11) $$

Combining (B10) and (B11) we get the Nusselt number in terms of the Rayleigh number in this case,

$$ Nu \sim Ra^{1/5} Pr^{1/5} \left( \frac{2}{m+1} \right)^{1/5} \frac{\Gamma^{6/5} \ln \Gamma}{(\Gamma-1)^{6/5}} (1+r_s)^{-6/5}. \qquad (B12) $$

REFERENCES

BRAGINSKY, S.I. & ROBERTS, P.H. 1995 Equations governing convection in Earth's core and the geodynamo. *Geophys. Astrophys. Fluid Dyn.* **79**, 1–97.

BROWNING, M.K., MIESCH, M.S., BRUN, A.S. & TOOMRE, J. 2006 Dynamo action in the Solar convection zone and tachocline: pumping and organization of toroidal fields. *Astrophys. J.* **648**, L157–L160.

BRUN, A.S. & TOOMRE, J. 2002 Turbulent convection under the influence of rotation: sustaining a strong differential rotation. *Astrophys. J.* **570**, 865–885.

CHONG, K.L., WAGNER, S., KACZOROWSKI, M., SHISHKINA, O. & XIA, K.-Q. 2018 Effect of Prandtl number on heat transport enhancement in Rayleigh–Bénard convection under geometrical confinement. *Phys. Rev. Fluids* **3**, 013501.

CURBELO, J., DUARTE, L., ALBOUSSIÈRE, T., DUBUFFET, F., LABROSSE, S. & RICARD, Y. 2019 Numerical solutions of compressible convection with an infinite Prandtl number: comparison of the anelastic and anelastic liquid models with the exact equations. *J. Fluid Mech.* **873**, 646–687.
GASTINE, T., WICHT, J. & AURNOU, J.M. 2015 Turbulent Rayleigh-Bénard convection in spherical shells. *J. Fluid Mech.* **778**, 721–764.

GLATZMAIER, G.A. & ROBERTS, P.H. 1995 A three-dimensional convective dynamo solution with rotating and finitely conducting inner core and mantle. *Phys. Earth Planet. Inter.* **91**, 63–75.

GOUGH, D.O. 1969 The anelastic approximation for thermal convection. *J. Atmos. Sci.* **26**, 448–456.

GROSSMANN, S. & LOHSE, D. 2000 Scaling in thermal convection: a unifying theory. *J. Fluid Mech.* **407**, 27–56.

GROSSMANN, S. & LOHSE, D. 2002 Prandtl and Rayleigh number dependence of the Reynolds number in turbulent thermal convection. *Phys. Rev. E* **66** (1), 016305.

JONES, C.A. & KUZANYAN, K. 2009 Compressible convection in the deep atmospheres of giant planets. *Icarus* **204**, 227–238.

JONES, C.A., BORONSKI, P., BRUN, A.S., GLATZMAIER, G.A., GASTINE, T., MIESCH, M.S. & WICHT, J. 2011 Anelastic convection-driven dynamo benchmarks. *Icarus* **216**, 120–135.

KESSAR, M., HUGHES, D.W., KERSALÉ, E., MIZERSKI, K.A. & TOBIAS, S.M. 2019 Scale selection in the stratified convection of the solar photosphere. *Astrophys. J.* **874**, 103–117.

LANTZ, S.R. & FAN, Y. 1999 Anelastic magnetohydrodynamic equations for modelling solar and stellar convection zones. *Astrophys. J. Suppl.* **121**, 247–264.

MENAUT, R., CORRE, Y., HUGUET, L., LE REUN, T., ALBOUSSIÈRE, T., BERGMAN, M., DEGUEN, R., LABROSSE, S. & MOULIN, M. 2019 Experimental study of convection in the compressible regime. *Phys. Rev. Fluids* **4**, 033502.

MIESCH, M.S., ELIOTT, J.R., TOOMRE, J., CLUNE, T.L., GLATZMAIER, G.A. & GILMAN, P.A. 2000 Three-dimensional spherical simulations of solar convection. I. Differential rotation and pattern evolution achieved with laminar and turbulent states. *Astrophys. J.* **532**, 593–615.

OGURA, Y. & PHILLIPS, N.A. 1962 Scale analysis of deep and shallow convection in the atmosphere. *J. Atmos. Sci.* **19**, 173–179.

QIU, X.-L. & TONG, P. 2001 Large-scale velocity structures in turbulent thermal convection. *Phys. Rev. E* **64** (3), 036304.

SIGGIA, E.D. 1994 High Rayleigh number convection. *Annu. Rev. Fluid Mech.* **26**, 137–168.

SHISHKINA, O., STEVENS, R.J.A.M., GROSSMANN, S. & LOHSE, D. 2010 Boundary layer structure in turbulent thermal convection and its consequences for the required numerical resolution. *New J. Phys.* **12** (7), 075022.

SILANO, G., SREENIVASAN, K.R. & VERZICCO, R. 2010 Numerical simulations of Rayleigh-Bénard convection for Prandtl numbers between $10^{-1}$ and $10^4$ and Rayleigh numbers between $10^5$ and $10^9$. *J. Fluid Mech.* **662**, 409–446.

SPIEGEL, E.A. & VERONIS, G. 1960 On the Boussinesq approximation for a compressible fluid. *Astrophys. J.* **131**, 442–447.

STEVENS, R.J.A.M., VAN DER POEL, E.P., GROSSMANN, S. & LOHSE, D. 2013 The unifying theory of scaling in thermal convection: the updated prefactors. *J. Fluid Mech.* **730**, 295–308.

TOOMRE, J., ZAHN, J.-P., LATOUR, J. & SPIEGEL, E.A. 1976 Stellar convection theory. II – single-mode study of the second convection zone in an A-type star. *Astrophys. J.* **207**, 545–563.

VERHOEVEN, J., WIESEHÖFER, T. & STELLMACH, S. 2015 Anelastic versus fully compressible turbulent Rayleigh-Bénard convection. *Astrophys. J.* **805**, 62.