Fully developed anelastic convection with no-slip boundaries

Abstract

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^7$ have been performed and our theoretical predictions are compared with the numerical results.

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_u$, $I_l$, $II_u$ and $II_l$ on the phase diagram, figure 2 of Grossmann & Lohse (2000), correspond to Rayleigh numbers less than about $10^{12}$ 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 $\varGamma$ of the temperature at the bottom of the layer $T_B$ to the temperature at the top $T_T$. When $\varGamma$ is close to unity the layer is nearly Boussinesq, but when $\varGamma$ is large there are many temperature and density scale heights within the layer. We aim to derive the functional form of $Nu(\varGamma,Ra,Pr)$ and $Re(\varGamma,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,

(2.1)$$\begin{gather} \frac{\partial\boldsymbol{u}}{\partial t}+\left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\boldsymbol{u}={-} \boldsymbol{\nabla}\left(\frac{p}{\bar{\rho}}\right)+\frac{g}{c_{p}}s \hat{\boldsymbol{e}}_{z}+\frac{\mu}{\bar{\rho}}\left[\nabla^{2}\boldsymbol{u}+ \frac{1}{3}\boldsymbol{\nabla}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}\right)\right], \end{gather}$$

(2.2)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\bar{\rho}\boldsymbol{u}\right)=0, \end{gather}$$

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

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

where here and below suffix notation is used,

(2.5)\begin{equation} q={\frac{\partial u_i}{\partial x_j} \frac{\partial u_i}{\partial x_j}}+\frac{1}{3}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}\right)^{2},\end{equation}

$\boldsymbol {u}$ being the fluid velocity, $p$ the pressure, $\rho$ 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 $\mu =\bar {\rho }\nu$, 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 $\boldsymbol {u} = 0$ there. We consider the constant entropy boundary conditions

(2.6a,b)\begin{equation} s= \Delta S \quad \textrm{at} \ z=0, \quad s = 0 \quad \textrm{at} \ z=d. \end{equation}

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

(2.7)\begin{equation} \left.\begin{gathered} \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{gathered}\right\} \end{equation}

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 } + \textrm {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 $\mathrm {d}\bar {p}/\mathrm {d}z=-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

(2.8a–c)$$\begin{gather} \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\left(m+1\right)}\left(1-\theta\frac{z}{d}\right)^{m+1}, \end{gather}$$

(2.8d–g)$$\begin{gather}\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 \varGamma= \frac{T_B}{T_T} = \frac{1}{1-\theta}, \end{gather}$$

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, $\varGamma >1$, is a convenient measure of the compressibility. Here $\varGamma \to 1$ is the Boussinesq limit, and highly compressible layers have $\varGamma$ large. Note that $\varGamma ^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(\varGamma -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 with a slightly more negative temperature gradient, so ${\tilde {T}}_B=T_B$, ${\tilde {T}}_T < T_T$. The conduction state is

(2.9a–c)$$\begin{gather} \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}}\left({\tilde m}+1\right)}\left(1-{\tilde \theta}\frac{z}{d}\right)^{{\tilde m}+1}, \end{gather}$$

(2.9d–f)$$\begin{gather}{\widetilde{\Delta{T}}}=T_{B}-{\tilde{T}}_{T}>0, \quad {\tilde{\theta}}=\frac{{\widetilde{\Delta T}}}{T_{B}}, \quad {\tilde{m}}=\frac{gd}{R {\widetilde{\Delta T}}} -1. \end{gather}$$

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

(2.10)\begin{equation} \epsilon = {\tilde \theta} \frac{{\tilde m}+1-\gamma {\tilde m}}{\gamma} ={-}\frac{d}{T_{B}}\left[\frac{\mathrm{d}\tilde{T}}{\mathrm{d}z}+\frac{g}{c_{p}}\right] \ll 1 , \end{equation}

and the entropy in the conduction state is

(2.11)\begin{align} {\tilde{s}}=c_{v} \ln\frac{\tilde{p}}{\tilde{\rho}^{\gamma}} + \mathrm{const.} =\frac{\epsilon c_p}{\theta} \ln\left(1-\theta\frac{z}{d}\right)+\mathrm{const.}, \ \textrm{so} \ s = \frac{c_p}{\theta} \ln\left(1-\theta\frac{z}{d}\right)+ \mathrm{const.}, \end{align}

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

(2.12)\begin{equation} \Delta S = \frac{ c_p}{\theta} \ln \varGamma = \frac{c_p \varGamma \ln \varGamma}{\varGamma-1}, \end{equation}

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)

(2.13)\begin{equation} -\frac{\boldsymbol{\nabla} p}{\bar{\rho}}-\frac{\rho}{\bar{\rho}}g\hat{\boldsymbol{e}}_{z} ={-}\boldsymbol{\nabla}\left(\frac{p}{\bar{\rho}}\right)+\frac{g}{c_{p}}s\hat{\boldsymbol{e}}_{z} + O(\epsilon) .\end{equation}

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 $\hat {T} = \bar {T} + \epsilon T$, then the total heat flux at the boundaries is $-k \,\textrm {d}{\langle \hat {T} \rangle _h} /\textrm {d}z = -k \,\textrm {d} \bar {T} /\textrm {d}z - k \epsilon \,\textrm {d}\langle T \rangle _h /\textrm {d}z$, but the superadiabatic part is obtained by subtracting off the heat flux carried along the adiabat, so

(2.14)\begin{equation} F^{{\scriptsize{super}}} = \left.- k \frac{\textrm{d} \left\langle T \right\rangle_h }{\textrm{d}z}\right|_{z=0}. \end{equation}

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

(2.15)\begin{equation} Nu = \frac{F^{{\scriptsize{super}}} d}{k T_B}, \end{equation}

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

(2.16)\begin{equation} 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}.\end{equation}

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

(2.17a,b) \begin{equation} \frac{k}{\rho_B c_p} \sim ( g d^3 \epsilon )^{1/2}, \quad \frac{\mu}{\rho_B } \sim (g d^3 \epsilon)^{1/2}{,} \end{equation}

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 $\vert \vert \ \vert \vert$ 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 \varGamma / \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

(2.18)\begin{equation} \Delta S_B + \Delta S_T = \Delta S. \end{equation}

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).

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.

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

(2.19a,b)\begin{equation} \frac{\left(\Delta\rho\right)_{i}}{\rho_{i}}\approx\frac{1}{\gamma-1}\left[\frac{\left(\Delta T\right)_{i}}{T_{i}}-\gamma\frac{\left(\Delta s\right)_{i}}{c_{p}}\right],\quad\frac{\left(\Delta p\right)_{i}}{p_{i}}\approx\frac{\gamma}{\gamma-1}\left[\frac{\left(\Delta T\right)_{i}}{T_{i}}-\frac{\left(\Delta s\right)_{i}}{c_{p}}\right], \end{equation}

where the $\varDelta$ 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 and viscous boundary layer thicknesses, $\delta _i^{th}$ and $\delta _i^{\nu }$ with $i=B,T$, which we use to obtain scaling estimates. Numerical simulations indicate that the horizontal velocity

(2.20)\begin{equation} U_H = ( \langle u_x^2 \rangle_h + \langle u_y^2 \rangle_h )^{1/2} \end{equation}

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

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

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^{\nu }$, can be defined as

(2.22)\begin{equation} \delta_i^{th} = \left\{\left. - \frac{1}{\Delta S_i} \frac{\textrm{d} \left\langle S \right\rangle_h}{\textrm{d}z} \right\vert_{z=z_i} \right\}^{{-}1}, \quad \delta_i^{\nu} = \left\{\left. \pm \frac{1}{\Delta U_i} \frac{\textrm{d}U_H }{\textrm{d}z} \right \vert_{z=z_i} \right\}^{{-}1}, {\quad z_i = 0, d}. \end{equation}

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^{\nu } U_i / d)$ so the vertical component of inertia is small. The boundary layers are therefore approximately hydrostatic,

(2.23)\begin{equation} \left(\Delta p\right)_{i}\approx\frac{g}{c_{p}}\rho_{i}\Delta s_{i}\delta_{i}^{th}.\end{equation}

Inserting (2.23) into (2.19b) leads to

(2.24)\begin{equation} \frac{\left(\Delta T\right)_{i}}{T_{i}}\approx\frac{\left(\Delta s\right)_{i}}{c_{p}}\left(1+\theta\frac{\delta_{i}^{th}}{d}\frac{T_{B}}{T_{i}}\right).\end{equation}

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

(2.25a,b)\begin{equation} \frac{\left(\Delta T\right)_{i}}{T_{i}}\approx\frac{\left(\Delta s\right)_{i}}{c_{p}} \quad \textrm{and} \quad \frac{\left(\boldsymbol{\nabla} T\right)_{i}}{T_{i}} \approx\frac{\left(\boldsymbol{\nabla} s \right)_{i}}{c_{p}}. \end{equation}

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

(2.26)\begin{equation} \delta_i^{T} = \left\{ \left.-\frac{1}{\Delta T_i} \frac{\textrm{d} \left\langle T \right\rangle_h}{\textrm{d}z} \right\vert_{z=z_i} \right\}^{{-}1}, \end{equation}

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),

(2.27)\begin{equation} Nu = \frac{d}{\delta^{th}_T} \frac{\Delta S_T}{\varGamma c_p} = \frac{d}{\delta^{th}_B} \frac{\Delta S_B}{c_p} . \end{equation}

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 $\textrm {d} \langle T \rangle _h /\textrm {d}z >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})} = s T_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

(2.28)\begin{equation} \vert \vert \rho \vert \vert = 0 \quad \Rightarrow \quad \langle p \rangle_{h,T} = \langle p \rangle_{h,B}. \end{equation}

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 } \boldsymbol{u}$ and averaging over the entire volume (recalling that horizontal averages of $x$ and $y$ derivatives are zero), we obtain the relation

(3.1)\begin{equation} \frac{g}{c_{p}} \vert \vert \bar{\rho}u_{z}s \vert \vert =\mu\vert \vert q\vert \vert ,\end{equation}

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$,

(3.2)\begin{align} {F^{{\scriptsize{super}}}}&= {\left\langle \bar{\rho}\bar{T}u_{z}s\right\rangle _{h} - k \frac{\textrm{d}\left\langle T \right\rangle_h }{\textrm{d}z} +\frac{g}{c_p}\int_{0}^{z}\left\langle \bar{\rho}u_{z}s\right\rangle _{h}\,\mathrm{d}z}\nonumber\\ &\quad -{\mu\int_{0}^{z}\left\langle q\right\rangle_h \,\mathrm{d}z -2\mu\left[ \left\langle u_{z} \frac{\textrm{d}u_{z}}{\mathrm{d}z} \right\rangle_{h}-\frac{m\Delta {\bar{T}}}{\bar{T} d}\left\langle u_{z}^{2}\right\rangle _{h}\right].} \end{align}

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

(3.3)\begin{align} \frac{\partial u_i}{\partial x_j} \frac{\partial u_j}{\partial x_i} - (\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}})^2 &= {\frac{\partial^2}{\partial x_i \partial x_j}(u_i u_j) - 2 \frac{\partial }{\partial x_j} (u_j \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}})} \nonumber\\ &={\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{align}

since the continuity equation gives

(3.4)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{u}} = \frac{ u_z m \Delta {\bar{T}}}{ {\bar{T}} d}, \end{equation}

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,

(3.5)\begin{equation} F^{{\scriptsize{super}}}= \left.- k \frac{\textrm{d}\left\langle T \right\rangle_h}{\textrm{d}z}\right|_{z=0}= \left.- k \frac{\textrm{d}\left\langle T \right\rangle_h}{\textrm{d}z}\right|_{z=d};\end{equation}

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$,

(3.6)\begin{align} \left\langle \bar{\rho}u_{z}s\right\rangle_h &= \left.-\frac{k}{T_B} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \right|_{z=0} + \left.\frac{k}{\bar{T}} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \right|_z - \int_0^z \frac{k \Delta {\bar{T}}}{{\bar{T}}^2 d} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \, \textrm{d}z + \int_0^z \frac{\mu}{\bar{T}} \left\langle q \right\rangle_h \, \textrm{d}z \nonumber\\ &\quad +\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 \, \textrm{d}z, \end{align}

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,

(3.7)\begin{equation} \left.\frac{k}{T_B} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \right|_{z=0} - \left.\frac{k}{T_T} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \right|_{z=d} ={-} \int_0^d \frac{k \Delta {\bar{T}}}{{\bar{T}}^2 d} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \, \textrm{d}z + \int_0^d \frac{\mu}{\bar{T}} \left\langle q \right\rangle_h \, \textrm{d}z. \end{equation}

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

(3.8)\begin{equation} S_{{\scriptsize{conv}}} = S_{{\scriptsize{diff}}} + S_{{\scriptsize{visc}}}, \end{equation}

where

(3.9)\begin{equation} S_{{\scriptsize{conv}}} = \left\langle \bar{\rho}u_{z}s\right\rangle_h, \end{equation}

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

(3.10)\begin{equation} S_{{\scriptsize{diff}}} =\left.- \frac{k}{T_B} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \right|_{z=0} + \left.\frac{k}{\bar{T}} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \right|_z - \int_0^z \frac{k \Delta {\bar{T}}}{{\bar{T}}^2 d} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \, \textrm{d}z \end{equation}

so that $S_{{\scriptsize {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

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

(3.12) $$\begin{gather} S_{{\scriptsize{diff}}} \approx - \left.\frac{k}{T_B} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \right|_{z=0} + \left.\frac{k}{T_T} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \right|_{z=d} \quad \textrm{if $z=d.$ } \end{gather}$$

We now turn to

(3.13)\begin{equation} S_{{\scriptsize{visc}}} = \int_0^z \frac{\mu}{\bar{T}} \left\langle q \right\rangle_h \,\textrm{d}z + {\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 \, \textrm{d}z.} \end{equation}

Because of the horizontal averaging, and using (2.2), (2.8a–c) and (3.4), the second integral can be written as

(3.14)\begin{align} &{\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 \, \textrm{d}z } \nonumber\\ &\quad = \int_{0}^{z} \frac{\mu}{\bar{T}} \left[ \frac{\partial^2}{\partial z^2} \left\langle u_{z}^{2}\right\rangle_h - \frac{m \Delta {\bar{T}}}{{\bar{T}} d} \frac{\partial}{\partial z} \left\langle u_{z}^{2}\right\rangle_h - \frac{m (\Delta {\bar{T}})^2}{{\bar{T}}^2 d^2} \left\langle u_{z}^{2}\right\rangle_h \right] \textrm{d}z . \end{align}

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^{\nu }$, so the velocity gradients appearing in $q$ are of magnitude $O(U/\delta _i^{\nu })$. Therefore, $q$ itself is $O(U^2/{(\delta _i^{\nu }})^2)$, and since the boundary layers are thin their contribution to the first integral in $S_{{\scriptsize {visc}}}$ is $O(\mu U^2/ {\bar {T}} \delta _i^{\nu })$. 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/ T d)$, which is $O( \delta _i^{\nu }/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

(3.15)\begin{equation} \left\langle \bar{\rho}u_{z}s\right\rangle_h \sim{-} \left.\frac{k}{T_B} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \right|_{z=0} + \int_0^z \frac{\mu}{\bar{T}} \left\langle q \right\rangle_h \, \textrm{d}z . \end{equation}

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

(3.16)\begin{equation} \left.\frac{k}{T_B} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \right|_{z=0} -\left. \frac{k}{T_T} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle T \right\rangle_h \right|_{z=d} = \frac{F^{{\scriptsize{super}}} \Delta {\bar{T}}}{T_B T_T} \sim \int_0^d \frac{\mu}{\bar{T}} \left\langle q \right\rangle_h \, \textrm{d}z. \end{equation}

Note also that in these expressions, the term $(\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {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),

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

in our dimensional units. Here $\Delta T$ is the superadiabatic temperature difference between the boundaries. In the Boussinesq limit $\varGamma \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 g s /c_p$, so

(3.18)\begin{equation} \frac{p}{\bar{p}} \sim \frac {g {\bar{\rho}} d}{{\bar{p}}} \frac{s}{c_p} \sim \frac{d}{H} \frac{s}{c_p}, \end{equation}

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

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

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 $\varGamma \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

(3.20)\begin{equation} -\left.\frac{\Delta {\bar{T}}}{T_B^2} \frac{\textrm{d}}{\textrm{d}z} \langle T \rangle_h \right\vert_{z=0} ={-} \frac{k \Delta {\bar{T}} }{T_B^2 d} \int_0^d \frac{\textrm{d}}{\textrm{d}z} \langle T \rangle_h \, \textrm{d}z + \frac{\mu}{T_B} \int_0^d \langle q \rangle_h \, \textrm{d}z, \end{equation}

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

(3.21)\begin{equation} Nu ={-}\left.\frac{d}{T_B} \frac{\textrm{d}}{\textrm{d}z} \langle T \rangle_h \right\vert_{z=0} \quad \textrm{or} \quad Nu ={-}\left.\frac{d}{\Delta T} \frac{\textrm{d}}{\textrm{d}z} \langle T \rangle_h \right\vert_{z=0}, \end{equation}

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 $g d / c_p$, so (3.20) becomes

(3.22)\begin{align} \frac{k}{c_p T_B} Nu g \alpha \Delta T = \frac{k}{c_p T_B} g \alpha \Delta T + \frac{\mu}{T_B} \int_0^d \langle q \rangle_h \,\textrm{d}z \quad \textrm{or} \quad (Nu - 1) g a \Delta T \kappa = \nu \int_0^d \langle q \rangle_h \,\textrm{d}z , \end{align}

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

(4.1)\begin{equation} Pr = \frac{\mu c_p}{k}. \end{equation}

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^{\nu }$, 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

(4.2)\begin{equation} \frac{\rho_i U_i^2}{d} \sim \frac{\mu U_i }{(\delta_i^{\nu} )^2 } \quad \textrm{so} \quad U_i \sim \frac{\mu d}{\rho_i (\delta_i^{\nu})^2 }. \end{equation}

For the entropy boundary layers, from (2.3),

(4.3)\begin{equation} \frac{\rho_i T_i U_i s}{d} \frac{\delta_i^{th}}{\delta_i^{\nu}} \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 \textrm{so} \quad U_i = \frac{k d\delta_i^{\nu}}{\rho_i c_p (\delta_i^{th})^3}, \end{equation}

where (2.25a,b) has been used and the factor $\delta _i^{th}/{\delta _i^{\nu }}$ 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

(4.4)\begin{equation} \frac{\delta_i^{\nu}}{\delta_i^{th}} \sim Pr^{1/3}, \end{equation}

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.

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.

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^{\nu }$ is no longer required, and (4.3) becomes

(4.5)\begin{equation} \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 \textrm{so} \quad U_i = \frac{k d }{\rho_i c_p (\delta_i^{th})^2}, \end{equation}

giving the ratio of the boundary layer thicknesses as

(4.6)\begin{equation} \frac{\delta_i^{\nu}}{\delta_i^{th}} \sim Pr^{1/2} \end{equation}

in the low-Prandtl-number case.

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

(5.1)\begin{equation} r_s = \frac{\Delta S_T}{\Delta S_B}, \quad \textrm{so} \quad \Delta S_B=\frac{\Delta S}{1+r_s} \quad \textrm{and the entropy in the bulk is} \ \frac{r_s}{1+r_s} \Delta S, \end{equation}

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

(5.2)\begin{equation} r_T = \frac{\Delta T_T}{\Delta T_B}. \end{equation}

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

(5.3)\begin{equation} r_u = \frac{U_T}{U_B}. \end{equation}

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 ^{th}_T$ and $\delta ^{th}_B$ with ratio

(5.4)\begin{equation} r_\delta = \frac{\delta^{th}_T}{\delta^{th}_B}. \end{equation}

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,

(5.5)\begin{equation} \left.-k \frac{\textrm{d}T}{\textrm{d}z}\right\vert_i \sim k \frac{\Delta T_i}{\delta^{th}_i} \Rightarrow r_\delta = r_T. \end{equation}

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,

(5.6)\begin{equation} \frac{\delta^{th}_T}{\delta^{th}_B} = \frac{\delta^{\nu}_T}{\delta^{\nu}_B} = r_\delta, \end{equation}

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

(5.7)\begin{equation} \frac{\rho_B U_B^2 }{d} \sim \frac{\mu U_B}{{\delta^{\nu}_B}^2}, \quad \frac{\rho_T U_T^2 }{d} \sim \frac{\mu U_T}{{\delta^{\nu}_T}^2} \quad \Rightarrow\quad r_u {r_\delta}^2 \sim \frac{\rho_B}{\rho_T} = \varGamma^{m}. \end{equation}

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_z s \rangle _h$ is approximately constant in the bulk interior. We now multiply the equation of motion (2.1) by ${\bar {\rho }} {\boldsymbol {u}}$ and average over the bulk interior, ignoring the small viscous term in the bulk, to get

(5.8) \begin{equation} \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}}\left\langle {\bar \rho} u_{z}s\right\rangle _{h} . \end{equation}

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_z s \rangle _h$ to be approximately the same just outside the two boundary layers,

(5.9) \begin{equation} \frac{\partial}{\partial z}(\bar{\rho}\langle u_{z}u^{2}\rangle _{h}) \vert_T \approx \frac{\partial}{\partial z}(\bar{\rho}\langle u_{z}u^{2}\rangle _{h})\vert_B, \end{equation}

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,

(5.10)\begin{equation} H_T = \frac{ d}{(m+1) ( \varGamma-1)}, \quad H_B = \frac{ \varGamma d}{ (m+1) (\varGamma-1)} \quad \Rightarrow \quad \frac{H_T}{H_B} = \varGamma^{{-}1}, \end{equation}

so that (5.9) gives

(5.11)\begin{equation} \frac{\rho_T u_T^3}{H_T} \sim \frac{\rho_B u_B^3}{H_B} \Rightarrow r_u^3 \sim \varGamma^{m-1} \Rightarrow r_u \sim \varGamma^{({m-1})/{3}}. \end{equation}

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

(5.12a–d)\begin{equation} r_s = \varGamma r_T, \quad r_\delta= r_T, \quad r_U {r_\delta}^2 = \varGamma^{m}, \quad {r_u} = \varGamma^{({m-1})/{3}}, \end{equation}

with solution

(5.13a–d)\begin{equation} r_u = \varGamma^{({m-1})/{3}}, \quad r_\delta = \varGamma^{({2m+1})/{6}}, \quad r_s = \varGamma^{({2m+7})/{6}}, \quad r_T = \varGamma^{({2m+1})/{6}}. \end{equation}

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

(5.14a–d)\begin{equation} r_u = \varGamma^{{1}/{6}}, \quad r_\delta = \varGamma^{{2}/{3}}, \quad r_s = \varGamma^{{5}/{3}}, \quad r_T = \varGamma^{{2}/{3}}. \end{equation}

Since $\varGamma$ 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 $\varGamma$, 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 $\varGamma$ 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 $\varGamma$ 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),

(6.1)\begin{equation} \frac{F^{{\scriptsize{super}}} \Delta {\bar{T}}}{T_B T_T} \sim \int_0^d \frac{\mu}{\bar{T}} \left\langle q \right\rangle_h \, \textrm{d}z. \end{equation}

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

(6.2)\begin{equation} q={\frac{\partial u_i}{\partial x_j} \frac{\partial u_i}{\partial x_j}}+\frac{1}{3}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}\right)^{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},\end{equation}

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

(6.3)\begin{equation} \frac{F^{{\scriptsize{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{\varGamma r_u^2}{r_\delta}\right){,} \end{equation}

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

(6.4)\begin{equation} F^{{\scriptsize{super}}} = \frac{k \Delta T_B}{\delta_B^{th}} = \frac{k T_B \Delta S}{c_p (1+r_s) \delta_B^{th}}. \end{equation}

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

(6.5)\begin{equation} \frac{k^2 Ra \varGamma}{c_p^2 (1+r_s) d^2 \rho_B^2} \frac{\delta_B^{\nu}}{\delta_B^{th}} \sim U_B^2 \left( 1+ \frac{\varGamma r_u^2}{r_\delta}\right). \end{equation}

Now we introduce the Reynolds number near the bottom boundary

(6.6)\begin{equation} Re_B = \frac{\rho_B U_B d}{\mu}, \end{equation}

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

(6.7)\begin{equation} \frac{Ra \varGamma}{Pr^2 (1+r_s) } \frac{\delta_B^{\nu}}{\delta_B^{th}} \sim Re_B^2 \left( 1+ \frac{\varGamma r_u^2}{r_\delta}\right). \end{equation}

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 $\varGamma$. In the high-Prandtl-number limit, applying (4.4) gives

(6.8)\begin{equation} Re_B \sim Ra^{1/2} Pr^{{-}5/6} \varGamma^{1/2} (1+r_s)^{{-}1/2} \left( 1+ \frac{\varGamma r_u^2}{r_\delta}\right)^{{-}1/2}, \end{equation}

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

(6.9)\begin{equation} Re_B \sim Ra^{1/2} Pr^{{-}3/4} \varGamma^{1/2} (1+r_s)^{{-}1/2} \left( 1+ \frac{\varGamma r_u^2}{r_\delta}\right)^{{-}1/2}. \end{equation}

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

(6.10)\begin{equation} Re_B \sim \left( \frac{d}{\delta_B^{\nu}} \right)^2, \end{equation}

while

(6.11)\begin{equation} Nu = \frac{d}{\delta_B^{th}} \frac{\varGamma \ln \varGamma}{(1+r_s)(\varGamma-1)}, \end{equation}

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

(6.12)\begin{equation} Re_B^{1/2} = \frac{\delta_B^{th}}{\delta_B^{\nu}} \frac{Nu (\varGamma - 1)(1+r_s)}{\varGamma \ln \varGamma}. \end{equation}

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,

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

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

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

while at low $Pr$, (4.6) gives

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

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

(6.16)\begin{equation} (1+r_s)^{{-}5/4} \left( 1+ \frac{\varGamma r_u^2}{r_\delta}\right)^{{-}1/4} = (1+\varGamma^{5/3})^{{-}5/4} (1+\varGamma^{2/3})^{{-}1/4}, \end{equation}

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

(6.17)\begin{equation} Nu \sim Ra^{1/4} Pr^{{-}1/12} \frac{ \varGamma^{5/4} \ln \varGamma}{ \varGamma-1} (1+\varGamma^{5/3})^{{-}5/4} (1+ \varGamma^{2/3})^{{-}1/4}\quad \textrm{for} \ Pr \gg 1 \end{equation}

and

(6.18)\begin{equation} Nu \sim Ra^{1/4} Pr^{1/8} \frac{ \varGamma^{5/4} \ln \varGamma}{ \varGamma-1} (1+\varGamma^{5/3})^{{-}5/4} (1+ \varGamma^{2/3})^{{-}1/4}\quad \textrm{for} \ Pr\ll 1. \end{equation}

The equivalent formulae for the Reynolds number are

(6.19)\begin{equation} Re_B \sim Ra^{1/2} Pr^{{-}5/6} \varGamma^{1/2} (1+\varGamma^{5/3})^{{-}1/2} (1+ \varGamma^{2/3})^{{-}1/2}\quad \textrm{for} \ Pr \gg 1 \end{equation}

and

(6.20)\begin{equation} Re_B \sim Ra^{1/2} Pr^{{-}3/4} \varGamma^{1/2} (1+\varGamma^{5/3})^{{-}1/2} (1+ \varGamma^{2/3})^{{-}1/2} \quad \textrm{for} \ Pr \ll 1. \end{equation}

As $\varGamma$ becomes large, $Nu$ decreases as $\ln \varGamma \varGamma ^{-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 $\varGamma$ 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, $\varGamma$ is close to unity and $\theta =(\varGamma -1)/\varGamma$ is small, so $\ln \varGamma / (\varGamma -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.33 Ra^{1/4} Pr^{-1/12}$ and $Re \approx 0.039 Ra^{1/2} Pr^{-5/6}$, so (6.14) becomes

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

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

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

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

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

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\,\textrm {d}{ \langle T \rangle }_h /\textrm {d}z$ by $-k_s {\bar {T}} \,\textrm {d}{\langle s \rangle }_h /\textrm {d}z$. In the entropy balance equation, entropy diffusion $S_{{\scriptsize {diff}}}$ changes to

(7.2)\begin{equation} S_{{\scriptsize{diff}}} =\left.- {k_s} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle s \right\rangle_h \right|_{z=0} + \left.{k_s} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle s \right\rangle_h \right|_z - \int_0^z \frac{k_s \Delta {\bar{T}}}{{\bar{T}} d} \frac{\mathrm{d}}{\mathrm{d}z} \left\langle s \right\rangle_h \, \textrm{d}z . \end{equation}

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=Re Pr$.

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, $\varGamma =1$, for comparison with the anelastic cases. The density stratification measured by $\varGamma$ varies over a moderate range only, because for the modest values of $Ra$ that are numerically accessible, large $\varGamma$ 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 $\varGamma \ln \varGamma / (\varGamma -1)$ and the constant is chosen so that it is zero at the bottom boundary $z=0$.

Figure 3. Horizontally averaged entropy (in units of $c_p$) and horizontal mean velocity profiles (Péclet number units) from the numerical simulations for $\varGamma =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$.

Figure 4. Horizontally averaged entropy $\langle S \rangle _h$ and horizontal mean velocity $U_H$ profiles for (a,b) $\varGamma =2.924$, $Ra=3 \times 10^6$, $Pr=1$, run A4: (c,d) $\varGamma =4.6416$, $Ra=6 \times 10^6$, $Pr=1$, run A5.

Table 1. Data from the numerical runs all corresponding to $m=3/2$ polytropes. The first four rows are the input parameters. Here $r_\delta$, $r_s$ and $r_u$ are the measured boundary layer ratios for each run. The velocities $U_T$ and $U_B$ are the local maxima at the edge of the boundary layers, measured in velocity units of $k/d \rho _B c_p$, 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_{Nu}$ 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_T$-theory and $Pe_B$-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_T$-nblr and $Pe_B$-nblr entries use the numerical boundary layer ratios.

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 $Re Pr = Pe$. We note that there is no great difference between the top and bottom horizontal velocities, consistent with the weak scaling with $\varGamma$ 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 sequence at constant $Pr=1$ with increasing $\varGamma$. The most noticeable feature is that at larger $\varGamma$ 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 $\varGamma$ in (5.14a–d). As expected, the boundary layers become thicker at larger $\varGamma$, so that at fixed $Ra$ the thin boundary layer assumption breaks down at large $\varGamma$.

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 $\textrm {d} \langle S \rangle _h /\textrm {d}z$ at $z=0$ and $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= \varGamma 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 $\varGamma =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 $\varGamma$ 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 $\varGamma$ 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

(7.3)\begin{equation} Nu = \left.- \frac{\textrm{d}}{c_p} \frac{\textrm{d} \langle S \rangle_h}{\textrm{d}z} \right\vert_{z=0} = \left.-\frac{\textrm{d}}{\varGamma c_p} \frac{\textrm{d} \langle S \rangle_h}{\textrm{d}z} \right\vert_{z=d}, \end{equation}

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 \geqslant 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^{1/4}$ 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 $\varGamma$ 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 $\varGamma$, 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 = Re Pr$. If the table 1 parameter values are inserted into (6.22a,b) with the value of $C_{Re}$ quoted there, the values of the Péclet number are consistently a factor of about five too small compared with the numerical values of $U_T$ and $U_B$ 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.085 Ra^{0.455}$. At the run A1 parameters this formula gives $Pe=251$ consistent with a prefactor of $C_{Re}=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 $C_{Re}=0.354$ for $Pr=10$, $C_{Re}=0.400$ for $Pr=1$ and $C_{Re}=0.421$ for $Pr=0.25$. Reassuringly, these are all quite close to the Qiu & Tong (2001) value of $C_{Re}=0.38$. We therefore use these three values of $C_{Re}$ at the appropriate Prandtl number in all our theory calculations. With these prefactors, the numerical results for $U_T$ and $U_B$ agree reasonably well with the predicted $Pe_T$-theory and $Pe_B$-theory results. The results have some scatter, which seems to reflect the scatter in our computed $r_U$. 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 Péclet numbers, though the theoretical Péclet numbers are generally a few percent lower that 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 $\varGamma$ which measures the strength of the stratification. In the Boussinesq limit, $\varGamma \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 $II_u$ and $II_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 $\varGamma$ 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 $\varGamma > 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 $\varGamma$ increases, or whether it becomes smaller at the top boundary than the bottom boundary at large $\varGamma$.

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.