Fully developed anelastic convection with no-slip boundaries

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 measures how stratified the layer is. 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 thin even when the stratification is large, the incompressible boundary layer case. 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.


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. (2008), Jones & Kuzanyan (2009), Jones et al. (2011), Verhoeven et al. (2015), Kessar et al. (2019) and many others. 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 distance d apart, and we assume the turbulent flow is spatially homogeneous in the horizontal direction. We also assume that the convection is in a statistically steady state, so that the time-averages of time-derivative terms can be neglected. 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, N u, 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, P r (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 developed later and updated in Stevens et al. 2013) developed for Boussinesq, i.e. nonstratified, convection. 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 to 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, 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 the Boussinesq convection.
In this paper, we concentrate on the convective regimes which seem to be most relevant to current numerical capabilities, i.e. regimes most easily achieved by numerical experiments. These are the regimes where the thermal dissipation is dominated by the thermal boundary layer contribution. It is those regimes, denoted by I u , I l , II u and II l on the phase diagram, figure 2 of Grossmann & Lohse (2000), that in the Boussinesq case correspond to Rayleigh numbers less than about 1012. 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 focussed on elucidating how well the anelastic approximation performs compared to 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 high Prandtl number are different.
In addition to the dependence on Rayleigh and Prandtl number, our problem depends on how stratified the layer is, which can be estimated by the ratio Γ of the temperature at the bottom of the the layer T B to the temperature at the top T T . When Γ is close to unity the layer is nearly Boussinesq, but when Γ is large there are many temperature and density scale heights within the layer. We aim to derive the functional form of N u(Γ, Ra, P r) and Re(Γ, Ra, P r), 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, which are being developed, but are not yet in a state to determine the prefactors accurately.
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 anrelastic 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 sections § 4 and § 5 we derive the physical arguments used to get 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. § 7 gives the results of our numerical simulations, comparing them with our theoretical results. Our conclusions are in § 8.

Fully developed compressible convection under the anelastic approximation
Consider a layer of compressible perfect gas between two parallel plates, of thickness d, periodic in horizontal directions, the evolution of which is described by the set of the Navier-Stokes, mass conservation and energy equations under the anelastic approximation, 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 µ =ρν, 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 = ∆S at z = 0, s = 0 at z = d.
(2.6) 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,p,ρ andT , are expanded in terms of the small parameter , which is defined precisely in equation (2.10) below, sô (2.7) wherep,ρ andT comprise the adiabatic reference state. Heres is simply a constant, and since s = c v lnp/ρ γ + const., andp/ρ γ is constant, we obtain (2.4b) by choosing the constant appropriately.

The adiabatic reference state
The reference state is the adiabatic static equilibrium governed by dp/dz = −gρ,p = RρT andp = Kρ γ , where R is the gas constant, z = 0 is the bottom of the layer and z = d the top. It 1995; Lantz & Fan, 1999 (2.13)

The Nusselt and Rayleigh numbers in anelastic convection
Next we consider the superadiabatic heat flux. The horizontal average at level z is denoted by h . At the boundaries, all the heat is carried by conduction, and if the total temperatureT = T + T , then the total heat flux at the boundaries is −kd T h /dz = −kdT /dz − k d T h /dz, but the superadiabatic part is obtained by subtracting off the heat flux carried along the adiabat, so we let (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 (2.15) so N u is close to unity near onset, and is large in fully developed convection. For fixed entropy boundary conditions, the Rayleigh number is defined as The anelastic approximation is asymptotically valid in the limit → 0. Note that small superadiabatic temperature gradient does not imply small Rayleigh number Ra, since the diffusion coefficients can be small, in fact to get Ra ∼ O(1) in the limit → 0 we must have 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 , the mass conservation equation by 1/2 and the energy equation by 3/2 . Having taken this limit, the parameter no longer appears in our analysis. However, if anelastic work is compared to fully compressible situations, then a finite value of 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.

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 h . The total entropy drop is the conduction state value ∆S = c p ln Γ/θ. Since the entropy is constant in the bulk interior, we define the entropy drops ∆S B and ∆S T across the bottom and top boundary layers respectively. These will not be equal, with ∆S T normally considerably larger than ∆S B . We must however have ∆S B + ∆S T = ∆S. (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.4) we obtain (2.19a,b) where the ∆ 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, δ th i and δ ν i with i = B, T , which we use to obtain scaling estimates. Numerical simulations indicate that the horizontal velocity 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, ∆U i , so where z = z max,B , z = z max,T are the locations of the local maxima. The thermal boundary layer thickness for the entropy, δ th i , and the viscous boundary layer thickness, δ ν i , can be defined as In the boudary 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 ∼ O(δ ν i ) so the vertical component of inertia is small. The boundary layers are therefore approximately hydrostatic, Inserting (2.23) into (2.19b) leads to Typically the term (θδ th i 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 fluctuation term is constant in both boundary layers, so that in the boundary layers and defining the temperature boundary layer thicknesses similarly to (2.22), 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.25), (2.27) In figure 1(b) we sketch the horizontally-averaged anelastic temperature perturbation T h . This is sometimes referred to as the superadiabatic temperature (e.g. Verhoeven et al. 2015). Note that with our constant entropy boundary conditions, T h is not zero at the boundaries. We show in appendix A that the offsets, T hB at z = 0 and T hT at z = d, are both positive and we show also that in the bulk, turbulent pressure effects make the gradient of T positive as shown in figure 1 (b). This means that the total horizontally averaged temperature gradient in the presence of convection is less negative than the adiabatic reference state, so that on horizontal average the layer is subadiabatically stratified (e.g. Korre et al. 2017). Of course, in some parts of the layer the local temperature gradient must be superadiabatic to drive the convection, but other parts are subadiabatic so that the horizontal average can be subadiabatic.
To obtain the anelastic temperature fluctuation as sketched in figure 1(b), we need to make use of equations (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 = ∆S. A different choice of entropy constant adds an easily found function of z to T , ρ 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 p h is required. Here we choose the natural condition that the anelastic density perturbation vanishes when integrated over the layer, that is 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.

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.

Energy balance
By multiplying the Navier-Stokes equation (2.1) byρu and averaging over the entire volume (recalling that horizontal averages of x and y derivatives are zero) we obtain the following relation g c p ||ρu z s|| = µ||q||, (3.1) stating that the total work 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 equation (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, In deriving this expression, we have made use of (2.8a,d) and since the continuity equation gives and we recall that x or y derivatives vanish on horizontal averaging. As z → d all terms with a factor u z tend to zero, so on using (3.1) we obtain the overall energy balance equation, thus in a stationary state the heat flux entering the fluid volume at z = 0 must be equal to the outgoing heat flux z = d.

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) byT , horizontally averaging, and integrating from 0 to z, where use has been made of equation (3.3). The overall entropy balance equation comes from taking the limit z → d of (3.6), noting u z → 0 in this limit, 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 the net entropy carried out of the region (0, z) by the convective velocity at level z, so that S dif f 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 T 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 1b 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 T h is O(∆T /d) whereas at the boundary it is O(∆T /δ th ), much larger since the boundary layer is thin. The third term is always small compared to the first, because the gradient is O(∆T /d) outside the boundary layers. The integrand is of order O(∆T /δ 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 We now turn to Because of the horizontal averaging, and using equations (2.2), (2.8) and (3.4), the second integral can be written (3.14) We now consider the magnitude of the terms in equation (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 δ ν i , so the velocity gradients appearing in q are of magnitude O(U/δ ν i ). q itself is therefore O(U 2 /(δ ν i ) 2 ), and since the boundary layers are thin their contribution to the first integral in S visc is O(µU 2 /T δ ν i ). In the boundary layers u z is small, but in the bulk we expect u z to be O(U ) and so u 2 z h is O(U 2 ). Because u 2 z h is horizontally averaged, it will vary on a length-scale d with z, so the gradient of u 2 z 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(µU 2 /T d), which is O(δ ν i /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, equation (3.6) is asymptotically equivalent to 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 → d we get (3.16) Note also that in these expressions, the term (∇ · u) 2 /3 in equation (2.5) makes a negligible contribution to (3.16) compared to the gradient terms by using (3.4).

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 Grossman & Lohse (2000) in the Boussinesq case, who start from equation (2.5) of Siggia (1994), in our dimensional units. Here ∆T is the superadiabatic temperature difference between the boundaries. In the Boussinesq limit Γ → 1, the basic state temperature and density tend to constant values T B and ρ B respectively, so the thermal diffusivity κ and kinematic viscosity ν are constants, κ = k/ρ B c p and ν = µ/ρ B . For a perfect gas the coefficient of expansion α = 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) we note that from the z-component of (2.1) p/ρd ∼ gs/c p , so where H is the pressure scale height. In the Boussinesq limit d/H becomes small; the Boussinesq limit Γ → 1 is the thin layer limit (Spiegel & Veronis, 1960). Then (2.4a,b) become (a) Thermal and viscous boundary layers in the case P r > 1. The thermal diffusion is smaller, so the thermal boundary layer is nested inside the viscous boundary layer. (b) The case P r < 1, where the viscous boundary layer is nested inside the thermal boundary layer.
so the entropy variable becomes equivalent to the superdiabatic temperature variable. The entropy jump ∆S across the layer can be written as a superadiabatic temperature jump ∆T = ∆S/αc p . As Γ → 1, from (2.12) ∆S/c p → 1 so ∆T → 1/α = T B . From energy conservation (3.5) the gradient of T h is the same at the top and bottom of the layer, so (3.7) can be written using the constancy ofT in the Boussinesq limit. From (2.14) and (2.15) 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∆T /d. From (2.8d) ∆T can be written gd/c p so (3.20) becomes (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).

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 characterized 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 thus 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 P r = µc p k . (4.1) In figure (2a) 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 δ th i U i /δ ν i , 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 For the entropy boundary layers, from (2.3) where (2.25) has been used and the factor δ th i /δ ν i arises because the horizontal velocity is reduced because the entropy boundary layer is thinner than the viscous boundary layer. Dividing (4.2) by (4.3) we obtain δ ν i δ th 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 (2b). 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 δ th i /δ ν i is no longer required, so (4.3) becomes giving the ratio of the boundary layer thicknesses as in the low Prandtl number case.

The boundary layer ratio problem
In Boussinesq convection, there is a symmetry about the mid-plane 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, and the top and bottom boundary layers can be very different, and the entropy of the bulk interior is significantly different from ∆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 , so ∆S B = ∆S 1 + r s and the entropy in the bulk is r s 1 + r s ∆S, (5.1) and the ratio of the anelastic temperature jumps across the boundary layers as We define the ratio of the velocities at the edge of the boundary layers as 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 δ th T and δ th B with ratio 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, 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, 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 entropy diffusion is small in the bulk interior, the convected entropy flux ρu z s h is approximately constant in the bulk interior. We now multiply the equation of motion (2.1) byρu and average over the bulk interior, ignoring the small viscous term in the bulk, to get Near the boundary layers, the pressure term p will be approximately constant as shown in (2.23), and since u z h = 0, the term involving u z p h will be small there, and we ignore it. Since we expect ρu z s h to be approximately the same just outside the two boundary layers, 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 doesn't matter whether the pressure scale height or the density scale height is used. Adopting the pressure scale height, so that equation (5.9) gives From the incompressible boundary layer equation (2.25) we have r s = Γ r T , so with the other three ratio equations (5.5), (5.7) and (5.11) we have r s = Γ r T , r δ = r T , r U r δ 2 = Γ m , r u = Γ In the case where m = 3/2, appropriate for ideal monatomic gas, (5.14) Since Γ 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 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 Γ , so again rather large values of Ra are required to establish the asymptotic trend. Note that in deriving equation (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 equation (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 Γ 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 Γ we did find a departure from the (5.14) 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 P r (Grossmann & Lohse, 2000) the equations (5.1-5.7) still hold, but the argument for equation (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.

The Nusselt number and Reynolds number scaling laws
When the boundary layers are thin, the overall entropy balance reduced to (3.16), 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 in these layers. So integrating over the boundary layers of thickness δ ν i and assuming T varies little in these layers, 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.25), (5.1) and (2.18) to get (6.4) Inserting this into (6.3) and using the definition (2.16) for the Rayleigh number, the entropy balance equation can be written (6.5)

Now we introduce the Reynolds number near the bottom boundary
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 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 Γ . In the high Prandtl number limit applying (4.4) gives Re B ∼ Ra 1/2 P r −5/6 Γ 1/2 (1 + r s ) −1/2 1 + Γ r 2 u r δ −1/2 , (6.8) while in the low Prandtl number case we get using (4.5) Re B ∼ Ra 1/2 P r −3/4 Γ 1/2 (1 + r s ) −1/2 1 + Γ r 2 u r δ −1/2 . (6.9) We now use the viscous boundary layer balance between advection and diffusion, (4.2), ρ B U B /d ∼ µ/(δ ν B ) 2 , to obtain a balance between N u and Re B . The boundary layer balance becomes using the expression (2.27) for the Nusselt number together with (2.12) and (5.1). Eliminating d/δ th B between these, (6.11) As above, the ratio of the boundary layer thicknesses can be evaluated in terms of the Prandtl number, and (6.11) allows us to eliminate Re B from (6.7) to obtain the Nusselt number as a function of Rayleigh number, At large P r, (4.4) gives the Nusselt number -Rayleigh number scaling in the form (6.13) while at low P r (4.6) gives N u ∼ Ra 1/4 P r 1/8 Γ 5/4 ln Γ Γ − 1 (1 + r s ) −5/4 1 + Γ r 2 u r δ −1/4 . (6.14) If we accept the ratio scalings derived in §5, in the case of a monatomic ideal gas, γ = 5/3, m = 3/2, we can write (1 + r s ) −5/4 1 + Γ r 2 u r δ −1/4 = 1 + Γ 5/3 −5/4 1 + Γ 2/3 −1/4 , (6.15) so as Γ becomes large, N u decreases as ln Γ Γ −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 P r. Analogous results for the case where the dissipation is mainly in the bulk rather than in the boundary layers, as can happen at low P r, are given in Appendix B.
In the Boussinesq limit, Γ is close to unity and θ = (Γ − 1)/Γ is small, so ln Γ/(Γ − 1) → 1 and ρ B → ρ T and T B → T T . Equation (2.1) reduces to the usual Boussinesq equation with s/c p replaced by αT , where for a perfect gas α = 1/T is the coefficient of expansion, consistent with (2.25). The total jump of entropy across the layer, ∆S, is replaced by the total temperature jump ∆T = ∆S/αc p so the Rayleigh number (2.16) reduces to the familiar form Ra = gα∆T d 3 /κν where κ = k/ρc p and ν = µ/ρ are the thermal diffusivity and kinematic viscosity respectively. These are both constant in the Boussinesq limit, and the ratios r u , r δ and r s all go to unity, see §5. Our scaling laws (6.8), (6.9), (6.13) and (6.14) 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 N u ≈ 0.33Ra 1/4 P r −1/12 and Re ≈ 0.039Ra 1/2 P r −5/6 , so (6.13) becomes (6.16) In the low P r limit where (6.14) applies, regime I l of Grossmann & Lohse (2000), they suggest a prefactor corresponding to C N u = 0.76 rather than 0.93. For the Reynolds number, (6.8) becomes Re B ≈ C Re Ra 1/2 P r −5/6 Γ 1/2 (1 + r s ) −1/2 1 + Γ r 2 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.

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, sō 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 −kd T h /dz by −k sT d s h /dz. In the entropy balance equation, entropy diffusion changes S dif f to 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 layers, when 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.13) 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 ρ B d 2 c p /k, the thermal diffusion time at the bottom of the layer. The velocities are in units of k/ρ B dc p , so they correspond to local Peclet numbers, where P e = ReP r.
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 Ra 10 6 3 × 10 6 10 7 3 × 10 6 6 × 10 6 10 6 3 × 10 6 10 6 10 6 10 6 10 6 Pr  192  332  601  295  341  219  260  144  200  260  149   TABLE 1. Data from the numerical runs all corresponding to m = 3/2 polytropes. The first four rows are the input parameters. r δ , rs and ru are the measured boundary layer ratios for each run. The velocities UT and UB are the local maxima at the edge of the boundary layers, measured in velocity units of k/dρBcp, i.e they are Peclet 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 equation (5.14). The N u-theory entries are based on (6.16) with the prefactors CNu as given in the text, and the boundary ratios come from (5.14).
The N u-nblr entries also use (6.16) with the same prefactors, but instead of using (5.14), the numerical boundary layer ratios (nblr) above are used. The P eT -theory and P eB-theory entries come from (6.8) and (6.9). The prefactors used are not those of Grossmann & Lohse (2000), see (6.17), but those given in the text. Again, (5.14) is used to determine the boundary layer ratios. The P eT -nblr and P eB-nblr entries use the numerical boundary layer ratios. 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, Γ = 1, for comparison with the anelastic cases. The density stratification measured by Γ varies over a moderate range only, because for the modest values of Ra that are numerically accessible, large Γ 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 Γ ln Γ/(Γ − 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 Peclet numbers to Reynolds numbers using ReP r = P e. We note that there is no great difference between the top and bottom horizontal velocities, consistent with the weak scaling with Γ 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 P r = 1 and P r = 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 P r 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 P r = 1 with increasing Γ . The most noticeable feature is that at larger Γ 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.14). 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 Γ in (5.14). As expected, the boundary layers become thicker at larger Γ , so that at fixed Ra the thin boundary layer assumption breaks down at large Γ .
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 S h /dz 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 ratios of the top and bottom entropy thicknesses and entropy jumps are denoted by r δ and r s in table 1. The velocity ratios at the top and bottom are denoted by r u . We can compare these with the predicted ratios in (5.14). 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.25), gives r s = Γ r δ , 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 Γ = 1.9438 had errors less than 10% (except the low P r 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 δ and r s large, but not quite as large as predicted by (5.14). 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 Γ 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 Γ is large. 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 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 N u appropriate for each Prandtl number used. This gives C N u = 0.949 for P r = 10, C N u = 0.785 for P r = 1 and C N u = 0.869 for P r = 0.25. Note that in table 1 this means that for the Boussinesq runs D1, D2, and D3, the N u-theory entries are the same as the actual N u entries by construction. None of these prefactors is very far from the values suggested by Grossmann & Lohse (2000), who had C N u = 0.93 at large P r and C N u = 0.76 at small Prandtl number, 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 P r = 1 case is slightly lower than the other values of C N u .
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.16) rather than the theoretically predicted ones gives the result N u-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.16) is inaccurate. For the density ratio 2.71 runs with P r 1, the boundary layer ratios were close to the predicted values, so not surprisingly, N unblr are not significantly better than N u-theory. 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.16). The least squares fit to a straight line in log(N u) 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 case 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 N u is only 61% of the numerical value for run B2, and in the case A5 predicted N u is 71% of the numerical N u. Part of this discrepancy is down to the boundary layer ratios, which become very large at high Γ , and so small inaccuracy can affect the Nusselt number significantly. If we use the numerical boundary layer ratios in (6.16) rather than the theoretical ones for run B2, the predicted N u-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.17), though it is convenient to express this in terms of Peclet number P e = ReP r . If the table 1 parameter values are inserted into (6.17) with the value of C Re quoted there, the values of the Peclet number are consistently a factor of about 5 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 Peclet 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 et al. 2010), use aspect ratios of 1 or less. The prefactor in (6.17) 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 (2002), using water (with P r = 5.5) in a cylinder of aspect ratio unity found Re = 0.085Ra 0.455 . At the run A1 parameters this formula gives P e = 251 consistent with a prefactor of C Re = 0.38 in (6.17), much larger than the Grossmann & Lohse (2002) value. They found very similar Re prefactors for both high and low P r. 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 P r = 10, C Re = 0.400 at P r = 1 and C Re = 0.421 at P r = 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 P e T -theory and P e 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 Peclet numbers, though the theoretical Peclet numbers are generally a few percent lower that the computed Peclet 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.

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 derived in terms of the Rayleigh number, the Prandtl number and the temperature ratio Γ which measures the strength of the stratification. In the Boussinesq limit, Γ → 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 Γ 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 §5 we proposed formulae based on a simple physical picture for these ratios. 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 Γ increases, or whether it becomes smaller at the top boundary than the bottom boundary at large Γ .
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 dependence of Re on the aspect ratio. When the prefactors were determined by normalizing on our Boussinesq runs, the issue was resolved.
We have focussed 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 similar problems which could be addressed which are of great astrophysical interest: 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.
written as q h ∼ U 3 H /2H where U H (z) is the horizontally averaged horizontal velocity and H is the local pressure scale height. We don't know how U H (z) is distributed in z, but we argued in §B1 above that ρU 3 H is approximately the same at the edge of both boundary layers, so a reasonable assumption for the purposes of estimation is that 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 (B 6) for the rest of this section. Equation (B 5) then becomes From ( , (B 9) and combining this with (B 8) and using the definition of the Prandtl number (4.1) we obtain N uRa P r 2 ∼ (m + 1) ln Γ 2 Re 3 B , (B 10) 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 P r, we only use (4.5) for the boundary layer ratio, so (6.11) becomes (Re B P r) 1/2 = N u(Γ − 1)(1 + r s ) Γ ln Γ .