3.1. FC model

For simplicity, we take the infinite Prandtl number approximation which eliminates inertia. This limit has been studied mathematically (Wang Reference Wang2004) and used for the study of mantle dynamics (Ricard Reference Ricard2015) for which Prandtl numbers are estimated around $10^{25}$: the effective kinematic viscosity of solids is much larger than their thermal diffusivity. Other objects, such as the Earth's core, the interior of stars and of gaseous planets have low Prandtl numbers (Schaeffer et al. Reference Schaeffer, Jault, Nataf and Fournier2017; Fuentes & Cumming Reference Fuentes and Cumming2020; Garaud Reference Garaud2020). For infinite Prandtl numbers, the governing equations of thermal convection are the following

(3.1)\begin{gather} \frac{\mathrm{D} \rho}{\mathrm{D} t} ={-} \rho {\boldsymbol{\nabla}} \boldsymbol{\cdot} \boldsymbol{u}, \end{gather}

(3.2)\begin{gather}\boldsymbol{0} ={-} {\boldsymbol{\nabla}} P + \rho \boldsymbol{g} + \eta \left[ {\boldsymbol{\nabla}}^2 \boldsymbol{u} +\tfrac{1}{3} {\boldsymbol{\nabla}} \left(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} \right) \right], \end{gather}

(3.3)\begin{gather}\rho T \frac{\mathrm{D} s}{\mathrm{D} t} = \dot{{\textsf{$\boldsymbol{\epsilon}$}}} : {\textsf{$\boldsymbol{\tau}$}} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( k {\boldsymbol{\nabla}} T \right) , \end{gather}

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

(3.4)\begin{equation} \tau _{ij} = 2 \eta \left[ \dot{\epsilon}_{ij} - \tfrac{1}{3} \left( {\boldsymbol{\nabla}} \boldsymbol{\cdot} \boldsymbol{u} \right) \delta_{ij} \right], \end{equation}

using Stokes’ assumption regarding the bulk viscosity.

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

(3.5)\begin{equation} \frac{1}{H L} \int \rho\,\mathrm{d} x \,\mathrm{d} z = \rho _0 . \end{equation}

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

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

(3.6)\begin{equation} - \rho ^2 T s' {\boldsymbol{\nabla}} \boldsymbol{\cdot} \boldsymbol{u} = \dot{{\textsf{$\boldsymbol{\epsilon}$}}} : {\textsf{$\boldsymbol{\tau}$}} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( k {\boldsymbol{\nabla}} T \right) , \end{equation}

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

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

(3.7)\begin{equation} c_{p0} ={-}a \frac{n}{n+1} \rho _0^n,\quad \mathrm{for}\ n>0 . \end{equation}

From (2.32) and (2.33), we derive an expression for $s'$ which is valid for both the logarithmic ($n=0$) and power-law ($n>0$) EoSs

(3.8)\begin{equation} s' (\rho ) ={-} (n+1) \frac{c_{p0}}{\rho _0} \left( \frac{\rho}{\rho _0} \right) ^{n-1},\quad \mathrm{for}\ n \geq 0 . \end{equation}

Similarly, a generic expression is obtained for the pressure gradient, from (2.35a,b) and (2.36a,b),

(3.9)\begin{equation} {\boldsymbol{\nabla}} P = (n+1) c_{p0} \left( \frac{\rho}{\rho _0} \right)^n \left[ (n+1)T {\boldsymbol{\nabla}} \rho + \rho {\boldsymbol{\nabla}} T \right],\quad \mathrm{for}\ n \geq 0 . \end{equation}

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

(3.10)\begin{gather} \frac{\mathrm{D} \rho}{\mathrm{D} t} ={-} \rho {\boldsymbol{\nabla}} \boldsymbol{\cdot} \boldsymbol{u}, \end{gather}

(3.11)\begin{gather}\boldsymbol{0} ={-} \frac{Ra_{sa} (n+1)}{\varepsilon } \left( \frac{ \rho ^n }{ \mathcal{D}} \left[ (n+1) T {\boldsymbol{\nabla}} \rho + \rho {\boldsymbol{\nabla}} T \right] + \rho \boldsymbol{e}_z \right) + {\boldsymbol{\nabla}}^2 \boldsymbol{u} + \frac{1}{3} {\boldsymbol{\nabla}} \left(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} \right), \end{gather}

(3.12)\begin{gather}0 ={-} (n+1) \rho ^{n+1} T \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} + \frac{\varepsilon\mathcal{D}}{Ra_{sa}} \dot{{\textsf{$\boldsymbol{\epsilon}$}}}: {\textsf{$\boldsymbol{\tau}$}} + \nabla ^2 T, \end{gather}

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

(3.13)\begin{gather} Ra_{sa} = \frac{\rho _0 g \alpha _0 {\rm \Delta} T_{sa} H^3}{\eta \kappa}, \end{gather}

(3.14)\begin{gather}\mathcal{D} = \frac{\alpha _0 g H}{c_{p0}}, \end{gather}

(3.15)\begin{gather}\varepsilon = \frac{{\rm \Delta} T_{sa}}{T_0}, \end{gather}

(3.16)\begin{gather}\alpha _0 T_0 = \frac{1}{n+1}. \end{gather}

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

(3.17)\begin{equation} {\rm \Delta} T_{sa} = T_h - T_c - \frac{\alpha _0 g T_0 H}{c_{p0}}. \end{equation}

When writing the dimensionless momentum equation (3.11), we use (3.9) to express the pressure gradient in terms of density and temperature gradient. When writing the dimensionless entropy equation (3.12), we use (3.6) and (3.8). It can be checked that the final set of dimensionless equations (3.10), (3.11) and (3.12) takes a generic form for any real value of $n \geq 0$: the case $n=0$ corresponds to the logarithmic relationship (2.32) whereas the cases $n>0$ correspond to the power laws (2.33). The choice of $n$ amounts to choosing the product $\alpha _0 T_0$, see (3.16).

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

(3.18)\begin{gather} \frac{\partial u_x}{\partial z} = 0, \quad \mathrm{when}\ z={\pm} \frac{1}{2}, \end{gather}

(3.19)\begin{gather}u_z = 0, \quad \mathrm{when}\ z={\pm} \tfrac{1}{2}, \end{gather}

(3.20)\begin{gather}{\int_{{-}L/(2H)}^{L/(2H)} u_x \left( x, z=\tfrac{1}{2} \right) \mathrm{d} x = 0, } \end{gather}

(3.21)\begin{gather}{T} = \frac{T_h}{T_0}, \quad \mathrm{when}\ z={-}\frac{1}{2}, \end{gather}

(3.22)\begin{gather}{T} = \frac{T_c}{T_0}, \quad \mathrm{when}\ z=\frac{1}{2}. \end{gather}

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

(3.23)\begin{gather} \frac{T_h}{T_0} = 1 + \frac{\mathcal{D} + \varepsilon }{2} , \end{gather}

(3.24)\begin{gather}\frac{T_c}{T_0} = 1 - \frac{\mathcal{D} + \varepsilon }{2} . \end{gather}

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

(3.25a–c)\begin{equation} \rho ' = \tilde{\rho} \,{\rm e}^{{\rm i} k x + \omega t},\quad T ' = \tilde{T} \,{\rm e}^{{\rm i} k x + \omega t},\quad u_x' = \tilde{u}_x \,{\rm e}^{{\rm i} k x + \omega t}, \end{equation}

where $\tilde {\rho }$, $\tilde {T}$ and $\tilde {u}_x$ are scalars. The governing equations (3.10), (3.11) and (3.12) are linearized near $z=0$ (the steady solution is nearly constant $T=1$) and lead to

(3.26)\begin{gather} \omega \tilde{\rho} ={-} {\rm i} k \tilde{u}_x , \end{gather}

(3.27)\begin{gather}0 ={-} \frac{Ra_{sa} (n+1)}{\varepsilon \mathcal{D}} \left[(n+1) {\rm i} k \tilde{\rho} + {\rm i} k \tilde{T} \right] - \frac{2}{3} k^2 \tilde{u}_x, \end{gather}

(3.28)\begin{gather}0 ={-} (n+1) {\rm i} k \tilde{u}_x - k^2 \tilde{T} . \end{gather}

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

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

admitting non-trivial solutions when

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

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

(3.31)\begin{equation} { | \omega | < } \frac{3}{2} \frac{Ra_{sa} (n+1)^2}{\varepsilon \mathcal{D}}, \end{equation}

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

(3.32)\begin{equation} \delta t \leq \frac{\varepsilon \mathcal{D}}{Ra_{sa} (n+1)^2} = \frac{\varepsilon \mathcal{D} \left( \alpha _0 T_0 \right) ^2}{Ra_{sa}}. \end{equation}

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

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

3.2. Anelastic approximation

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

(3.33)\begin{gather} \rho _a = 1, \end{gather}

(3.34)\begin{gather}T_a = 1 - \mathcal{D} z , \end{gather}

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

(3.35)\begin{gather} {\boldsymbol{\nabla}} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

(3.36)\begin{gather}\boldsymbol{0} ={-}\frac{Ra_{sa}}{\mathcal{D}} {\boldsymbol{\nabla}} \tilde{P} + Ra_{sa} \tilde{s} {\hat{\boldsymbol{e}}}_z + {\boldsymbol{\nabla}}^2 \boldsymbol{u} , \end{gather}

(3.37)\begin{gather}\frac{ \mathrm{D} }{\mathrm{D} t} (T_a \tilde{s} ) ={-} \mathcal{D} u_z \tilde{s} + \frac{\mathcal{D}}{Ra_{sa}} \dot{\epsilon} : \tau + \nabla ^2 \tilde{T}. \end{gather}

For our EoS, $\tilde {s}$ is proportional to $\tilde {\rho }$ (see (3.8)) and linearizing (2.35a,b) or (2.36a,b) leads to

(3.38)\begin{gather} \tilde{s} ={-} (n+1) \tilde{\rho} , \end{gather}

(3.39)\begin{gather}\tilde{\rho}={-} \frac{\tilde{T}}{(n+1) T_a} + \frac{\tilde{P}}{(n+1)^2 T_a}, \end{gather}

and, therefore,

(3.40)\begin{equation} \tilde{s} = \frac{\tilde{T}}{T_a} - \frac{\tilde{P}}{(n+1)T_a}. \end{equation}

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

(3.41)\begin{equation} \tilde{T} \left( z ={\pm} \tfrac{1}{2} \right) ={\mp} \tfrac{1}{2} . \end{equation}

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

(3.42)\begin{equation} \int _0^{L/H} \tilde{P} \left( x, z=\tfrac{1}{2} \right) - \tilde{P} \left( x, z={-} \tfrac{1}{2} \right) \mathrm{d} x = 0 . \end{equation}

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

3.3. Anelastic liquid approximation

In that approximation, departures of entropy from the adiabatic profile are considered to be due only to temperature departures, whereas departures in pressure are neglected in (3.40) (Anufriev et al. Reference Anufriev, Jones and Soward2005). The governing equations are still (3.35), (3.36) and (3.37) where $\tilde {s}$ is changed in each instance into

(3.43)\begin{equation} \tilde{s} = \frac{\tilde{T}}{T_a} . \end{equation}

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

(3.44)\begin{gather} \frac{\partial \tilde{T}}{\partial x} = 0, \quad z ={\pm} \frac{1}{2}, \end{gather}

(3.45)\begin{gather}\tilde{T} \left( x,z={-} \tfrac{1}{2} \right) - \tilde{T} \left( x,z= \tfrac{1}{2} \right) = 1 . \end{gather}

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

(3.46)\begin{equation} \int _0^{L/H} \tilde{P} \left( x, z={\pm} \tfrac{1}{2} \right) \mathrm{d} x = 0 . \end{equation}

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

3.4. Simple compressible approximation

We now introduce a new approximation aiming at getting a very simple system of equations where compressible effects are still present. In the ALA approximation, the adiabatic temperature profile appears explicitly in the equations and we consider replacing $T_a$ by a constant value equal to $1$, the value of $T_a$ in the mid-plane of the cavity. We certainly lose connection to thermodynamics with that move, but compressible work is still present and it will be interesting to investigate which compressible effects are still well accounted for in this approximation. Note that this SCA model is equivalent to a version of the ‘extended Boussinesq approximation’ (EBA) (King et al. Reference King, Lee, van Keken, Leng, Zhong, Tan, Tosi and Kameyama2010), where the background density is assumed to be uniform (with our EoS, this is the case of all our anelastic models) and where the background temperature is also assumed to be uniform. The SCA model is still defined by (3.35), (3.36) and (3.37), where the expression for entropy (3.43) is changed for

(3.47)\begin{equation} \tilde{s} = \tilde{T}, \end{equation}

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

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

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

The complete set of equations FC, AA, ALA and SCA, with boundary conditions, are written in their explicit forms in Appendix A. The parameters of all simulations are listed in Appendix B. The files used for each convection model FC, AA, ALA and SCA are provided as supplemental materials available at https://doi.org/10.1017/jfm.2022.216 or available at github https://github.com/RayleighBenardModels/Compressible_PLG.

5.1. Change of temperature profile with moderate compressibility

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

Figure 4.

Time and horizontal averaged superadiabatic temperature profiles along the vertical direction $z$, for (a) $Ra_{sa}=10^6$ and (b) $Ra_{sa}=10^8$, for $\mathcal {D}$ between $0$ and $0.4$, $\alpha _0 T_0 = 1$, obtained in the AA.

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

Figure 5.

Snapshots of superadiabatic temperature fields for $Ra_{sa}=10^7$, $\alpha _0 T_0 = 1$, with (a) $\mathcal {D}=0.05$ and (b) $\mathcal {D}=0.2$ obtained in the AA.

One should not think that compressible effects are always associated with stronger descending plumes. This property seems to be related to the EoS. In another paper (Ricard et al. Reference Ricard, Alboussière, Labrosse, Curbelo and Dubuffet2022), we consider an EoS suitable for planetary mantles and in that case, the opposite is true: with compressible effects, ascending plumes are stronger.

5.2. Asymmetry at larger compressibility

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

Figure 6.

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

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

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

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

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

Figure 7.

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

7.1. Global dissipation

The total dissipation is expressed as a fraction of the heat flux, more specifically the superadiabatic convective heat flux. The reference result is obtained in the Boussinesq approximation, where it has been shown that the integrated viscous dissipation is equal to the product of the dissipation number and the convective heat flux (Howard Reference Howard1963; Hewitt, McKenzie & Weiss Reference Hewitt, McKenzie and Weiss1975). We express the total viscous dissipation as a fraction of the conduction heat flux along the superadiabatic gradient and denote it $D_\nu$. In the FC model, the integrated viscous dissipation divided by ${\rm \Delta} T_{sa}$ has the following expression

(7.19)\begin{equation} D_\nu = \frac{\mathcal{D}}{Ra_{sa} \varepsilon} \left\langle \dot{{\textsf{$\boldsymbol{\epsilon}$}}} : {\textsf{$\boldsymbol{\tau}$}} \right\rangle . \end{equation}

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

(7.20)\begin{equation} D_\nu = \frac{\mathcal{D}}{Ra_{sa}} \left\langle \dot{{\textsf{$\boldsymbol{\epsilon}$}}} : {\textsf{$\boldsymbol{\tau}$}} \right\rangle, \end{equation}

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

(7.21)\begin{equation} E = D_\nu / (Nu -1). \end{equation}

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

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

Figure 9.

Excess of the ratio of dissipation to convective heat flux relative to $\mathcal {D}$, as a function of $\mathcal {D}$, for FC, AA, ALA and SCA models, $\alpha _0 T_0 = 1$, $\mathcal {D}$ in $[0.25, 0.5,\ 0.75,\ 1.0,\ 1.25,\ 1.5]$, $Ra_{sa}$ in $[10^3,\ 10^{3.5}, 10^4,\ 10^{4.5},\ 10^5,\ 10^{5.5},\ 10^6,\ 10^{6.5}]$.

Figure 10.

Excess of the ratio of dissipation to convective heat flux relative to $\mathcal {D}$, as a function of $Ra_{sa}$ (same values as in figure 9).

Let us consider some limits to the dissipation results. First, we have rigorous upper (and lower) bounds, since the papers of Backus (Reference Backus1975) and Hewitt et al. (Reference Hewitt, McKenzie and Weiss1975). Taking (3.37), dividing by $T_a$ and integrating by parts leads to

(7.22)\begin{equation} \frac{\mathrm{D} \tilde{s} }{\mathrm{D} t} - \boldsymbol{\boldsymbol{\nabla}} \boldsymbol{\cdot} \left( \frac{\boldsymbol{\boldsymbol{\nabla}} \tilde{T}}{T_a} \right) = \frac{ \mathcal{D}}{Ra_{sa}} \frac{ \dot{\epsilon} : \tau}{T_a} + \frac{\boldsymbol{\boldsymbol{\nabla}} \tilde{T} \boldsymbol{\cdot} \boldsymbol{\boldsymbol{\nabla}} T_a }{T_a^2} , \end{equation}

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

(7.23)\begin{equation} \frac{ D_\nu }{1 + \dfrac{\mathcal{D}}{2}} < Nu \left[\frac{1}{ 1 - \dfrac{\mathcal{D}}{2}} - \frac{1}{ 1 + \dfrac{\mathcal{D}}{2}} \right] , \end{equation}

implying the upper bound

(7.24)\begin{equation} \frac{D_\nu}{Nu -1} < \mathcal{D} + \frac{\mathcal{D}^2}{ 2 - \mathcal{D} }, \end{equation}

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

Another possible limit case is shown to correspond to our numerical results at large Rayleigh number in § 7.2. That limit is that of a vanishing contribution of the $G(z)$ component of the heat flux defined in (7.12) in the AA. It corresponds to the case when the heat flux is carried entirely by the flux of entropy $T_a \overline {\tilde {s} u_z}$ (outside top and bottom conduction layers) whereas the energy dissipation at each height is $\mathcal {D} \overline {\tilde {s} u_z}$, as can be seen from (7.13) and (7.14). Under that assumption, the integral of energy dissipation is equal to

(7.25)\begin{equation} D_\nu = Nu \int _{1/2}^{1/2} \frac{\mathcal{D}}{T_a} \, \mathrm{d} z , \end{equation}

leading to the expression

(7.26)\begin{equation} \frac{D_\nu}{Nu -1} = \log \left(\frac{1+\mathcal{D}/2}{1-\mathcal{D}/2} \right) .\end{equation}

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

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

7.2. Convergence of heat flux and dissipation profiles

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

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

Figure 11.

Entropy flux contribution for $\mathcal {D}$ equal to (a) $0.8$ and (b) $1.6$, for a superadiabatic Rayleigh number up to $10^9$ in the AA and $\alpha _0 T_0 = 1$.

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

Figure 12.

Dissipation profile for $\mathcal {D}$ equal to (a) $0.8$ and (b) $1.6$, $\alpha _0 T_0 = 1$, for a superadiabatic Rayleigh number up to $10^9$ in the AA.

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

Figure 13.

Relative distance ($L^1$-norm) of the entropy to total heat flux profile (a) and of the dissipation profile to the limit (C6).

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

Figure 14.

Snapshot of the superadiabatic temperature field for (a) a very small $\mathcal {D}=0.05$, (b) moderate $\mathcal {D}=0.4$ to (c) a large dissipation number $\mathcal {D}=1.6$, for $Ra_{sa} = 10^9$, in the AA with $\alpha _0 T_0 = 1$. In the enlarged views, the distribution of viscous dissipation $\dot {\epsilon }:\tau$ (left) and entropy flux $\tilde {s} u_z$ (right) are shown.

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

Figure 15.

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

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

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