Upper bound of heat flux in an anelastic model for Rayleigh–Bénard convection

Abstract

Bounds on heat transfer have been the subject of previous studies concerning convection in the Boussinesq approximation: in the Rayleigh–Bénard configuration, the first result obtained by Howard (J. Fluid Mech., vol. 17, issue 3, 1963, pp. 405–432) states that the dimensionless heat flux $\textit {Nu}$ carried out by convection is such that $\textit {Nu} < (3/64 \ Ra)^{1/2}$ for large values of the Rayleigh number $Ra$, independently of the Prandtl number $Pr$. This is still the best-known upper bound, only with the prefactor improved to $\textit {Nu} -1 < 0.02634 \ Ra^{1/2}$ by Plasting & Kerswell (J. Fluid Mech., vol. 477, 2003, pp. 363–379). In the present paper, this result is extended to compressible convection. An upper bound is obtained for the anelastic liquid approximation, which is similar to an anelastic model used in astrophysics based on a turbulent diffusivity for entropy. The anelastic bound is still scaling as $Ra^{1/2}$, independently of $Pr$, but depends on the dissipation number $\mathcal {D}$ and on the equation of state. For monatomic gases and large Rayleigh numbers, the bound is $\textit {Nu} < 25.8\, Ra^{{1}/{2}} / (1-\mathcal {D}/2 )^{{5}/{2}}$.

1. Introduction

An important landmark of fluid mechanics has been to show that rigorous upper bounds could be obtained from the governing equations on quantities such as energy dissipation (or pressure gradient) in a pipe with a given flow rate, or heat flux between walls maintained at different temperatures (Howard 1963). Concerning Rayleigh–Bénard convection, the method of Howard (1963) consisted of identifying two integral constraints based on energy conservation and on the limit of an entropy balance in the Boussinesq approximation. In the large space of all possible temperature fields, not necessarily satisfying the pointwise governing equations, but satisfying the boundary conditions and the integral constraints, an upper bound on the heat flux was determined. A second method, the so-called ‘background’ method, was developed in the 1990s by Doering & Constantin (1996) and is based on a decomposition of the temperature field into an arbitrary vertical profile (the background profile) satisfying the boundary conditions and an homogeneous three-dimensional, time-dependent field. A spectral condition is said to hold when the ‘dissipation’ contained in the background profile (in fact, the $L^2$ norm of its derivative) is larger than the total possible dissipation of the convective flow. This spectral condition has been shown to be related to the same eigenvalue problem as that involved in the energy stability (Joseph 1976) of the background profile. The problem is finally turned into finding the background profile with the minimum possible dissipation. Kerswell (1999) proved that both methods were dual approaches to the same problem of optimisation. A third method was obtained recently by Seis (2015), with a more intuitive approach. The average heat flux must be constant over height, but cannot be carried by conduction after some distance to the bottom wall and convection must take over. This implies that sufficiently strong vertical velocity components exist there, which is necessarily associated with deformation (hence, viscous dissipation) as those vertical components are zero at the bottom. However, the total viscous dissipation is related to the heat flux and imposes a limit to the convective flux. That constraint leads to the same scaling as that of Howard. Finally, Chernyshenko (2022) makes the connection between what he calls the auxiliary functional method (Howard's method can be considered as a very particular instance of this more general method), the background method (by Doering & Constantin) and the direct method (a formulation of Seis's method).

Upper bounds of the heat flux have not been derived for compressible convection until now. Here, in § 2, we consider a simple model of compressible convection, the anelastic liquid approximation (Anufriev, Jones & Soward 2005). As an anelastic model, acoustic modes have been filtered out. Moreover, entropy is supposed to depend on the superadiabatic temperature only, so that pressure is not a relevant thermodynamic variable. An anelastic model is also used in astrophysics (Lantz & Fan 1999) with nearly the same equations as the anelastic liquid model, although the path to get there has taken a different direction. From a general anelastic model, a subgrid model for turbulence is used to change the conduction term (gradient of temperature) into a gradient of entropy. Again, this anelastic model depends only on a single thermodynamic variable, entropy. We also use a simple equation of state, that of the ideal gases. Finally, we consider a simple geometry, that of a plane layer, in a uniform gravity field perpendicular to the plane layer. The horizontal extent of the layer can be infinite or finite. The vertical depth of the layer is such that compressible effects range from negligible (Boussinesq limit) to extreme values (adiabatic temperature profile reaching zero kelvin at the top).

The maximum principle for parabolic equations plays an important role in our derivation of an upper bound. This also plays a crucial role in the work by Seis (2015). In the Boussinesq model, temperature is bounded below by the cold temperature imposed at the top and bounded above by the hot temperature imposed at the bottom. In a compressible model, adiabatic compression and decompression as well as viscous heating imply that these limits not longer hold for temperature. Instead, we show in § 3 that entropy has a minimum value imposed at the top boundary but no obvious maximum value. That property will be used several times in the paper.

In § 4 we derive an equation for the logarithm of entropy (up to a constant), a quantity that we call $\log$-entropy. We show that, similarly to the entropy flux, the flux of that $\log$-entropy increases with height, or decreases only slightly. Otherwise stated, the sources of $\log$-entropy are either positive (possibly unbounded) or slightly negative (bounded from below). With that equation, we can bound the integral of the gradients of the $\log$-entropy in the layer. The derivation follows then the same principle as that of Seis (2015). Two cases have to be discussed, a case of small Nusselt number and a case of large Nusselt number that is expected to hold when the Rayleigh number is large. In the latter case, a lower bound of vertical velocities (root mean square) is found to be necessary to carry the flux of $\log$-entropy at some finite distance to the lower boundary, leading to a lower bound of viscous dissipation required to generate these velocity components. Coming back to the entropy equation (not $\log$-entropy) in § 5, we obtain an upper bound for dissipation. As shown in § 6, the condition that the lower bound is less than the upper bound of dissipation leads to an upper bound for the heat flux in terms of the governing parameters. Formally, the absolute upper bound is the maximum of the small Nusselt number case and the upper bound obtained under the assumption of a larger Nusselt. As discussed in § 7, the case of large Nusselt numbers should apply to most, if not all, convective flows. In Appendices A and B, we discuss the choice of constants and optimisation of the bound. Appendix C is devoted to the application of our method to an anelastic model with turbulent thermal diffusivity.

2. Governing equations and dimensionless numbers

The fluid (a monatomic ideal gas) is contained in a horizontal layer, between altitudes $0$ and $d$, in a uniform gravity field ${\boldsymbol {g}}$. A superadiabatic temperature difference $\Delta T_{sa}$ is imposed in addition to the adiabatic gradient between the bottom and top boundaries. The governing equations in the anelastic liquid approximation are written in a dimensionless form as follows (see Anufriev et al. 2005):

(2.1)$$\begin{gather} {\boldsymbol{\nabla}} \boldsymbol{\cdot} (\rho_a \boldsymbol{v}) = 0 , \end{gather}$$

(2.2)$$\begin{gather}\frac{\rho_a}{Pr} \frac{\mathrm{D} \boldsymbol{v}}{\mathrm{D} t} ={-} \rho_a {\boldsymbol{\nabla}} \left( \frac{ P}{\rho_a}\right) + Ra \rho_a s \boldsymbol{e}_z + {\boldsymbol{\nabla}} \boldsymbol{\cdot} \tau , \end{gather}$$

(2.3)$$\begin{gather}\rho_a T_a \frac{\mathrm{D} s}{\mathrm{D} t} = \frac{\mathcal{D}}{Ra } \dot{\varepsilon} : \tau + \nabla^2 T , \end{gather}$$

where $\boldsymbol {e}_z$ is the vertical unit vector ($\boldsymbol {e}_x$ and $\boldsymbol {e}_y$ are the horizontal unit vectors) and where the dimensionless governing parameters are the Prandtl $Pr$, Rayleigh $Ra$ and dissipation $\mathcal {D}$ numbers,

(2.4a–c)\begin{equation} Pr = \frac{\eta c_p}{k}, \quad Ra= \frac{\rho_0^2 c_p g \Delta T_{sa} d^3}{T_0 \eta k}, \quad \mathcal{D} = \frac{g d}{c_p T_0}, \end{equation}

where the viscosity $\eta$, the thermal conductivity $k$ and the heat capacity $c_p$ of the gas are uniform and constant. The dimensionless tensors of deformation rate and stress, in the Stokes approximation of zero bulk viscosity (proven correct for monatomic ideal gases; Emanuel 1998), are the following:

(2.5)$$\begin{gather} \dot{\varepsilon}_{ij} = \tfrac{1}{2} \left( \partial_i v_j + \partial_j v_i \right), \end{gather}$$

(2.6)$$\begin{gather}\tau_{ij} = 2 \dot{\varepsilon}_{ij} - \tfrac{2}{3} ( \partial_k v_k) \delta_{ij}. \end{gather}$$

The average temperature $T_0$ and average density $\rho _0$ of the adiabatic profiles are chosen to express dimensionless temperature and density adiabatic profiles (hydrostatic, isentropic) as follows (Curbelo et al. 2019):

(2.7)$$\begin{gather} T_a (z) = 1 - \mathcal{D} \left(z - \tfrac{1}{2} \right) , \end{gather}$$

(2.8)$$\begin{gather}\rho_a (z) = \frac{\mathcal{D}/ \left( 1 - \gamma ^{{-}1} \right) }{\left( 1+\mathcal{D}/2 \right)^{{\gamma}/({\gamma -1})} - \left( 1-\mathcal{D}/2 \right)^{{\gamma}/({\gamma -1})}} \left[ T_a (z) \right] ^{{1}/({\gamma -1})}, \end{gather}$$

where $\gamma = c_p/c_v$ is the ratio of heat capacities (for example, $\gamma = 5/3$ for monatomic gases). The dimensionless gradient of adiabatic temperature is $- \mathcal {D} {\boldsymbol {e}}_z$, uniform and vertical.

Superadiabatic temperature $T$ and entropy $s$ are scaled using $\Delta T_{sa}$ and $c_p \Delta T_{sa}/T_0$, respectively. Space coordinates $(x,y,z)$, time $t$, velocity $\boldsymbol {v}$, rate of deformation tensor $\dot {\varepsilon }$, stress tensor $\tau$ and pressure $P$ are made dimensionless using $d$, $\rho _0 c_p d^2 / k$, $k/(\rho _0 c_p d)$, $k/(\rho _0 c_p d^2)$, $k \eta /(\rho _0 c_p d^2)$ and $k \eta /(\rho _0 c_p d^2)$, respectively.

In the anelastic liquid model, entropy $s$ (or, rather, the superadiabatic entropy in addition to a uniform base value) is assumed to depend on superadiabatic temperature $T$ only,

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

A consequence is that pressure $P$ has no effect on thermodynamic variables. In this model, pressure $P$ is only a Lagrange multiplier associated with the conservation of mass (2.1).

The boundary conditions are given by constant values of temperature or entropy on bottom and top boundaries. In terms of velocity, we impose that the normal component vanishes on both horizontal boundaries and that no work is done by the boundaries. Both non-slip (zero tangential velocity components) or no-stress in the horizontal direction are acceptable.

(2.10)$$\begin{gather} z=0: \quad v_z = 0 \quad \mathrm{and} \quad ( v_x = v_y = 0 \quad \mathrm{or}\quad \partial_z v_x = \partial_z v_y = 0), \end{gather}$$

(2.11)$$\begin{gather}z=0: \quad T = \frac{1}{2} \quad \mathrm{or}\quad s = \frac{1}{2 + \mathcal{D}}, \quad \mathrm{since}\ T_a = 1 + \frac{\mathcal{D}}{2} , \end{gather}$$

(2.12)$$\begin{gather}z=1: \quad v_z = 0 \quad \mathrm{and} \quad (v_x = v_y = 0 \quad \mathrm{or}\quad \partial_z v_x = \partial_z v_y = 0), \end{gather}$$

(2.13)$$\begin{gather}z=1: \quad T ={-}\frac{1}{2} \quad \mathrm{or}\quad s ={-}\frac{1}{2 - \mathcal{D}}, \quad \mathrm{since}\ T_a = 1 - \frac{\mathcal{D}}{2} . \end{gather}$$

In astrophysics, the Fourier law for thermal conduction is replaced by a subgrid model of turbulent diffusion for entropy. This has the consequence that the term $\nabla ^2 T$ in (2.3) is changed for ${\boldsymbol {\nabla }} \boldsymbol {\cdot } ( c \boldsymbol {\nabla } s )$, where the coefficient $c$ accounts for turbulence at small scale. In this paper, we stick to the usual conduction term $\nabla ^2 T$, but the treatment of the ‘turbulent diffusivity’ is discussed in Appendix C.

The parameter $\mathcal {D}$ is the one associated with compressibility, as better seen in the relationship $\mathcal {D} = \alpha g d / c_p = (1-\gamma ^{-1}) \rho g d / K_T$, proportional to the ratio of hydrostatic pressure to incompressibility $K_T$. Its range, $0 < \mathcal {D} < 2$, covers all cases from the Boussinesq limit ($\mathcal {D} \rightarrow 0$) to the most-extreme case of compressibility ($\mathcal {D} \rightarrow 2$) where a temperature of 0 K and a vanishing density are reached at the top of the layer, see (2.7) and (2.8).

We define the Nusselt number as $\textit {Nu} = - \mathrm {d}_z \bar {T} (z=0)$ the average superadiabatic heat flux injected at the bottom and extracted at the top (an overline $\bar {X}$ on any quantity $X$ denotes its horizontal and time average). The additional heat flux conducted along the adiabat does not affect convection and is here uniform and equal to $\mathcal {D} T_0 / \Delta T_{sa}$ in the same dimensionless scale as the superadiabatic heat flux. The proof that the average heat flux conducted at the top is equal to the average heat flux conducted at the bottom can be found in several papers, including Curbelo et al. (2019) where this was also tested in numerical simulations. This reflects the consistency of the set of the governing equations. The left-hand side of (2.3) can be expanded as $\rho _a T_a \,\mathrm {D} s / \mathrm {D} t = \rho _a \,\mathrm {D} (T_a s) / \mathrm {D} t - \rho _a (\mathrm {d} T_a / \mathrm {d} z) v_z s$ and the last term of adiabatic heating (with the expression for the adiabatic gradient $\mathrm {d} T_a / \mathrm {d} z = - \mathcal {D}$) can be seen to balance the average viscous dissipation because it corresponds to the same balance as in the dot product of velocity with the momentum equation (2.2), between the power of buoyancy forces and viscous dissipation.

3. A minimum principle

In the anelastic liquid approximation (2.9), (2.3) can be rewritten using entropy $s$ alone:

(3.1)\begin{equation} \rho_a T_a \frac{\mathrm{D} s}{\mathrm{D} t} = \frac{\mathcal{D}}{Ra } \dot{\varepsilon} : \tau + T_a \nabla^2 s + 2 \boldsymbol{\boldsymbol{\nabla}} T_a \boldsymbol{\cdot}\boldsymbol{\nabla} s , \end{equation}

because $\nabla ^2 T_a = 0$, see (2.7), for our choice of an ideal gas equation of state and uniform gravity. The entropy equation has a suitable form for a maximum principle (Picone 1929; Nirenberg 1953). The second-order operator is elliptic and even uniformly elliptic as the coefficient $T_a$ is above a positive constant in the whole domain, for any choice of the dissipation parameter $0 < \mathcal {D} < 2$. Furthermore, the term of viscous dissipation $(\mathcal {D}/{Ra })\, \dot {\varepsilon } : \tau$ is positive (or zero) everywhere and at all times. It follows from the maximum principle that $s$ cannot take a value smaller than a value it takes at a boundary or at an initial time. As we are interested in statistically stationary solutions, we argue either that the memory of the initial time is lost (the argument developed by Foias, Manley & Temam (1987) can be applied here to prove that initial values smaller than that imposed at the boundaries will eventually be erased by diffusion), or that the initial condition is chosen such that it does not contain values of the entropy $s$ lower than those at the boundaries. We are left with the conclusion that entropy must be larger, everywhere and at all times, than the value assigned at the top boundary,

(3.2)\begin{equation} s \geqslant - \frac{1}{2-\mathcal{D}}. \end{equation}

4. A $\log$-entropy equation

Let us define a constant $s_0 = 1/(2-\mathcal {D}) + 8/(4-\mathcal {D}^2)$, such that $s+s_0$ is always positive from the maximum principle and satisfies

(4.1)\begin{equation} s+s_0 \geqslant \frac{8}{4-\mathcal{D}^2}. \end{equation}

This choice of $s_0$ will simplify subsequent calculations but we show in Appendix A that an optimisation of this choice only leads to a mild improvement. Let us divide (3.1) by $T_a$ and by $s+s_0$, two positive terms. After rearranging some terms, we obtain

(4.2)\begin{equation} \rho_a \frac{\mathrm{D} \mathcal{L} }{\mathrm{D} t} = \frac{\mathcal{D}}{Ra } \frac{ \dot{\varepsilon} : \tau }{T_a (s +s_0 )} + \left| \boldsymbol{\nabla} \mathcal{L} \right| ^2 + 2 \mathrm{d}_z ( \ln T_a ) \partial_z \mathcal{L} + {\nabla}^2 \mathcal{L} , \end{equation}

where $\mathcal {L} = \ln (s +s_0 )$. With our choice of $s_0$, the value of $\mathcal {L}$ on the bottom and top boundaries are $\mathcal {L} (z=0) = \ln (12/(4- \mathcal {D}^2 ))$ and $\mathcal {L} (z=1) = \ln (8/(4- \mathcal {D}^2 ))$. Averaging over horizontal directions and in time, and taking into account (2.1), we obtain an equation for the vertical flux of $\mathcal {L}$,

(4.3)\begin{equation} \mathrm{d}_z \varPhi_{\mathcal{L}} = \frac{\mathcal{D}}{Ra \ T_a} \overline{\frac{ \dot{\varepsilon} : \tau }{s +s_0 } } + \overline{ \left| \boldsymbol{\nabla} \mathcal{L} \right| ^2 } + 2 \mathrm{d}_z ( \ln T_a ) \mathrm{d}_z \bar{\mathcal{L}} , \end{equation}

where $\varPhi _{\mathcal {L}}$ is defined as the average vertical flux of $\mathcal {L}$, at any height $z$, as follows:

(4.4)\begin{equation} \varPhi_{\mathcal{L}} (z) ={-}\mathrm{d}_z \bar{\mathcal{L}} + \rho_a \overline{v_z \mathcal{L}}. \end{equation}

We now average equation (4.2) over the whole layer and in time. This is equivalent to integrating (4.3) between $z=0$ and $z=1$. Our objective here is to obtain an integral bound on $| \boldsymbol {\nabla } \mathcal {L} |^2$. Let us first consider the integral of the last term in (4.3)

(4.5)\begin{align} \int_0^1 2 \mathrm{d}_z ( \ln T_a ) \mathrm{d}_z \bar{\mathcal{L}} \,\mathrm{d} z &={-}2 \mathcal{D} \int_0^1 \frac{\mathrm{d}_z \bar{\mathcal{L}}}{T_a} \mathrm{d} z, \nonumber\\ &={-}2 \mathcal{D} \left( \left[ \frac{\bar{\mathcal{L}}}{T_a} \right]_0^1 - \int_0^1 \mathcal{D} \frac{\bar{\mathcal{L}}}{T_a^2} \mathrm{d} z \right) . \end{align}

The first term is known from the boundary conditions (2.11) and (2.13). We use the minimum principle (4.1), implying $\bar {\mathcal {L}} \geqslant \ln (8/(4-\mathcal {D}^2))$, to bound the last term in (4.5) so that we have

(4.6)\begin{equation} \int_0^1 2 \mathrm{d}_z ( \ln T_a ) \mathrm{d}_z \bar{\mathcal{L}} \,\mathrm{d} z \geqslant \frac{4 \mathcal{D}}{2 + \mathcal{D}} \ln \left( \frac{3}{2} \right) . \end{equation}

The term $\mathrm {d}_z (\rho _a \overline {v_z \mathcal {L}})$ has no integral contribution because $v_z$ vanishes at the bottom and top, we simply have to evaluate the contribution of the diffusion term $-\mathrm {d}_z\bar {\mathcal {L}}$ in the integral of (4.3), given that $\partial _z \mathcal {L} = \partial _z (T/T_a )/ (s+s_0 )$,

(4.7)\begin{equation} - \left[ \mathrm{d}_z \bar{\mathcal{L}} \right]_0^1 = \frac{\textit{Nu}}{12} ( 2 + 5 \mathcal{D} ) + \frac{\mathcal{D}}{4} \frac{2 + \mathcal{D} }{2 - \mathcal{D}} + \frac{\mathcal{D}}{6} \frac{2 - \mathcal{D} }{2 + \mathcal{D}} . \end{equation}

Since the term involving viscous dissipation is positive in (4.3), combining (4.6) and (4.7) leads to the following bound:

(4.8)\begin{equation} \left\langle \left| \boldsymbol{\nabla} \mathcal{L} \right| ^2 \right\rangle \leqslant \frac{\textit{Nu}}{12} ( 2 + 5 \mathcal{D} ) + \frac{\mathcal{D}}{4} \frac{2 + \mathcal{D} }{2 - \mathcal{D}} + \frac{\mathcal{D}}{6} \frac{2 - \mathcal{D} }{2 + \mathcal{D}} - \frac{4 \mathcal{D}}{2 + \mathcal{D}} \ln \left( \frac{3}{2} \right) , \end{equation}

where the bracket denotes time and space average (horizontal and vertical), so that $\langle X \rangle = \int _0^1 \bar {X} \,\mathrm {d} z$ for any variable $X$. The sum of the last two terms is less than zero for all values of $\mathcal {D}$, hence the sum of the last three terms is less than $\mathcal {D}/(2 - \mathcal {D})$ so that we have

(4.9)\begin{equation} \left\langle \left| \boldsymbol{\nabla} \mathcal{L} \right| ^2 \right\rangle \leqslant \frac{\textit{Nu}}{12} ( 2 + 5 \mathcal{D} ) + \frac{\mathcal{D} }{2 - \mathcal{D}} . \end{equation}

This bound is linear in the Nusselt number with a coefficient ranging from $1/6$ at small $\mathcal {D}$ to $1$ when $\mathcal {D}$ reaches its maximum value $2$. In addition, the other term depends on $\mathcal {D}$ only and diverges towards infinity when $\mathcal {D}$ approaches $2$.

At $z=0$, the diffusive part of the $\log$-entropy flux, $- \mathrm {d}_z \bar {\mathcal {L} }$, carries the whole flux:

(4.10)\begin{equation} \varPhi_{\mathcal{L}} (0) ={-} \mathrm{d}_z \bar{\mathcal{L} } (0) = \frac{2 - \mathcal{D}}{6} \left( \textit{Nu} - \frac{\mathcal{D}}{2 + \mathcal{D}} \right) . \end{equation}

Let us first consider the large Nusselt number case, precisely such that

(4.11)\begin{equation} \textit{Nu} > \frac{\mathcal{D}}{2 + \mathcal{D}} + 24 \frac{ \ln \left( \dfrac{3}{2} \right) }{ 2 - \mathcal{D}} . \end{equation}

The flux $\varPhi _{\mathcal {L}} (0)$ is therefore strictly positive and this ensures that $\bar {\mathcal {L}}$ (see (4.10)) is locally decreasing at $z=0$. If it kept decreasing at the same rate as in (4.10), then the flux $\varPhi _{\mathcal {L}}$ would still be carried by diffusion at higher values of $z$. However, $\bar {\mathcal {L}}$ cannot decrease below the minimum value $\ln ( 8 / (4- \mathcal {D}^2) )$, limiting the extension of the diffusive region, and indicating that the convective part of the flux (4.4) must take over. We define the height $\delta$ as the smallest value of $z>0$ where the diffusive component, $-\mathrm {d}_z \bar {\mathcal {L} }$, becomes less than half of $\varPhi _{\mathcal {L} } (0)$,

(4.12)\begin{equation} - \mathrm{d}_z \bar{ \mathcal{L} } (\delta) \leqslant \frac{\varPhi_{\mathcal{L}} (0)}{2} \quad \mathrm{and} \quad - \mathrm{d}_z \bar{ \mathcal{L} } (z) > \frac{\varPhi_{\mathcal{L}} (0)}{2} \quad \mathrm{for} \ 0 \leqslant z < \delta. \end{equation}

It is shown in Appendix B that this choice of half-value is better than any other fraction of the bottom flux. The value of $\bar {\mathcal {L}}$ at $z=\delta$ is

(4.13)\begin{align} \bar{\mathcal{L}} (\delta) &=\bar{\mathcal{L}} (0) + \int_0^\delta {\rm d}_z \bar{\mathcal{L}} \, \mathrm{d} z , \nonumber\\ & < \ln \left( \frac{12}{4 - \mathcal{D}^2} \right) - \frac{\delta}{2} \varPhi_{\mathcal{L}} (0). \end{align}

Let us define $\delta _0$ as

(4.14)\begin{equation} \delta_0 = \frac{12 \ln \left( \dfrac{3}{2} \right)}{(2 - \mathcal{D}) \left( \textit{Nu} - \dfrac{\mathcal{D}}{2 + \mathcal{D}} \right) }, \end{equation}

which is ensured to be less than $1/2$ from (4.11). We prove by contradiction that $\delta$ exists and is less than $\delta _0$. Suppose that this is not the case, then from (4.13) and using (4.10), we would have $\bar {\mathcal {L}} (\delta _0) < \ln ( 8/(4- \mathcal {D}^2 ) )$ which would be less than the minimum value for $\bar {\mathcal {L}}$. This is impossible, hence there exists a value of $\delta$ satisfying (4.12).

The definition of $\delta$ implies that $\bar {\mathcal {L}}$ is decreasing everywhere in the range $0 \leqslant z \leqslant \delta$. This ensures the positivity of the last term of (4.3) in that interval, so that $\varPhi _{\mathcal {L}} (\delta ) \geqslant \varPhi _{\mathcal {L}} (0)$. Because the diffusive part has been divided by a factor two, this means that the convective part of the flux of the $\log$-entropy, at $z=\delta$, must be at least half the boundary flux (4.10)

(4.15)\begin{equation} \rho_a \overline{v_z \mathcal{L}} (\delta) \geqslant \frac{2 - \mathcal{D}}{12} \left( \textit{Nu} - \frac{\mathcal{D}}{2 + \mathcal{D}} \right) . \end{equation}

Using continuity (2.1) which implies that the horizontal average of $v_z$ is zero at all heights, the property $\rho _a \leqslant \rho _{a0}$ and a Cauchy–Schwarz inequality, the convective flux can be bounded as follows:

(4.16)\begin{equation} \rho_a \overline{v_z \mathcal{L}} (\delta) \leqslant \rho_{a0} \overline{v_z (\mathcal{L}-\mathcal{L}_0)} \leqslant \rho_{a0} \sqrt{ \overline{v_z^2} } \sqrt{\overline{\left( \mathcal{L} - \mathcal{L}_0 \right) ^2 } } , \end{equation}

where $\rho _{a0}$ and $\mathcal {L}_0$ are the density of the adiabatic profile and the value of $\mathcal {L}$ at $z=0$, whereas $v_z$ and $\mathcal {L}$ are evaluated at $z=\delta$. Now, we use the gradients of $\mathcal {L}$ to bound $\mathcal {L} (\delta )$:

(4.17)\begin{equation} \mathcal{L} (\delta) = \mathcal{L}_0 + \int_0^\delta \partial_z \mathcal{L} \,\mathrm{d} z . \end{equation}

Using a Cauchy–Schwarz inequality, we have

(4.18)\begin{equation} \left( \mathcal{L}(\delta) - \mathcal{L}_0 \right) ^2 \leqslant \delta \int_0^\delta \left( \partial_z \mathcal{L} \right) ^2 \,\mathrm{d} z . \end{equation}

Taking time and horizontal average, extending the last integral over the whole volume and including all gradient components, we obtain

(4.19)\begin{equation} \overline{\left( \mathcal{L}(\delta) - \mathcal{L}_0 \right) ^2} \leqslant \delta \left\langle \left| \boldsymbol{\nabla} \mathcal{L} \right| ^2 \right\rangle \leqslant \delta_0 \left\langle \left| \boldsymbol{\nabla} \mathcal{L} \right| ^2 \right\rangle , \end{equation}

where $\langle | \boldsymbol {\nabla } \mathcal {L} | ^2 \rangle$ is itself bounded from (4.9). Following the same steps, we can bound $v_z$, at any height $z$, as

(4.20)\begin{equation} \overline{ v_z^2 } (z) \leqslant z \int_0^{z} \overline{(\partial_z v_z )^2} \,\mathrm{d} \zeta . \end{equation}

Now, we need to relate $\overline {(\partial _z v_z )^2}$ to the mean viscous dissipation. From the general expression of viscous dissipation with zero bulk viscosity (Landau & Lifshitz 1966),

(4.21)\begin{equation} \dot{\varepsilon} : \tau = \frac{1}{2} \sum_{i=1}^3 \sum_{j=1}^3 \left[ \partial_i v_j + \partial_j v_i - \frac{2}{3} (\partial_k v_k ) \delta_{ij} \right]^2, \end{equation}

retaining only the three ‘diagonal’ terms $i=j$ among the nine terms, finally dropping $2(\partial _x v_x )^2$ and $2(\partial _y v_y )^2$, we derive

(4.22)\begin{equation} \dot{\varepsilon} : \tau \geqslant 2 (\partial_z v_z )^2 - \tfrac{2}{3} (\partial_k v_k )^2. \end{equation}

Using the anelastic equation of continuity (2.1), this leads to an upper bound of $(\partial _z v_z )^2$,

(4.23)\begin{equation} (\partial_z v_z )^2 \leqslant \frac{1}{2} \dot{\varepsilon} : \tau + \frac{1}{3} \frac{\mathcal{D}^2 v_z^2}{(\gamma -1 )^2 T_a^2} , \end{equation}

which is substituted in (4.20) to obtain

(4.24)\begin{equation} \overline{v_z^2} (z) \leqslant \frac{z}{2} \left\langle \dot{\varepsilon} : \tau \right\rangle + \frac{\mathcal{D}^2}{3 (\gamma -1)^2} z \int_0^{z} \frac{\overline{v_z^2}}{T_a^2} \mathrm{d} \zeta . \end{equation}

We now use this equation for $z \leqslant \delta \leqslant \delta _0$, given that $\delta _0 \leqslant 1/2$ from (4.11) so that $1/T_a^2$ in the integral above is smaller than $1$, and bounding the same integral by extending its range from $\zeta = 0$ to $\zeta = \delta$ for all $z$. Integrating (4.24) from $0$ to $\delta$ leads to the bound

(4.25)\begin{equation} \int_0^\delta \overline{v_z^2} (z) \,\mathrm{d} z \leqslant \frac{\delta ^2}{4} \left\langle \dot{\varepsilon} : \tau \right\rangle + \frac{\mathcal{D}^2}{3 (\gamma -1)^2} \frac{\delta ^2}{2} \int_0^\delta \overline{v_z^2} (z) \,\mathrm{d} z, \end{equation}

which is used to express $\int _0^\delta \overline {v_z^2} (z) \,\mathrm {d} z$ in terms of $\langle \dot {\varepsilon } : \tau \rangle$ so that (4.24) can finally be written at $z=\delta$,

(4.26)\begin{equation} \overline{v_z^2} (\delta ) \leqslant \frac{\delta_0 }{2} \left\langle \dot{\varepsilon} : \tau \right\rangle \frac{1}{1 - \dfrac{\mathcal{D}^2 \delta_0^2}{6 (\gamma -1)^2 }} , \end{equation}

where the denominator of the last fraction can be checked to be positive when $\delta _0 \leqslant 1/2$ which has been shown to follow from (4.11) and the definition of $\delta _0$ (4.14).

With (4.9), (4.14), (4.16), (4.19) and (4.26), we can essentially write (4.15) as a lower bound for the viscous dissipation:

(4.27)\begin{align} \left\langle \dot{\varepsilon} : \tau \right\rangle &\geqslant \frac{2 \,\textit{Nu}^3 (2 - \mathcal{D})^4 }{12^3 \left( \ln \dfrac{3}{2} \right)^2 \rho_{a0}^2 (2 + 5 \mathcal{D}) }\nonumber\\ &\quad \times \frac{\left[ 1 - \dfrac{\mathcal{D}}{\textit{Nu} (2 + \mathcal{D})} \right]^4 \left[ 1 - \dfrac{24 \left( \ln \dfrac{3}{2} \right)^2 \mathcal{D}^2}{(2-\mathcal{D})^2 \left( \textit{Nu} - \dfrac{\mathcal{D}}{2 + \mathcal{D}} \right)^2 (\gamma -1 )^2 } \right]}{1 + \dfrac{12 \mathcal{D}}{\textit{Nu} (2 - \mathcal{D} ) (2 + 5 \mathcal{D})}} . \end{align}

This bound is valid as long as the large-Nusselt-number condition (4.11) is valid. For very large Nusselt numbers (more precisely, $(2- \mathcal {D}) \textit {Nu} \gg 1$), this lower bound becomes

(4.28)\begin{equation} \left\langle \dot{\varepsilon} : \tau \right\rangle \geqslant \frac{2 \, \textit{Nu}^3 (2 - \mathcal{D})^4 }{12^3 \left( \ln \dfrac{3}{2} \right)^2 \rho_{a0}^2 (2 + 5 \mathcal{D})}. \end{equation}

Finally, we come back to the condition imposed on the Nusselt number (4.11) and consider the alternative case of small Nusselt numbers,

(4.29)\begin{equation} \textit{Nu} \leqslant \frac{\mathcal{D}}{2 + \mathcal{D}} + 24 \frac{ \ln \left( \dfrac{3}{2} \right) }{ 2 - \mathcal{D}} . \end{equation}

This expression is itself an upper bound for the Nusselt number. Whenever another upper bound obtained under assumption (4.11) gets lower than (4.29), it must be replaced by (4.29).

5. Upper bound on viscous dissipation for a given heat flux

We have thus obtained (4.27) a lower bound of the viscous dissipation $\langle \dot {\varepsilon } : \tau \rangle$ for a given heat flux $\textit {Nu}$. Now, considering the classical entropy budget, we are going to obtain an upper bound for $\langle \dot {\varepsilon } : \tau \rangle$. Let us divide the governing equation (2.3) by $T_a$, rearrange the diffusion term, and obtain the usual anelastic entropy equation

(5.1)\begin{equation} \rho_a \frac{\mathrm{D} s}{\mathrm{D} t} = \frac{\mathcal{D}}{Ra } \frac{\dot{\varepsilon} : \tau }{T_a} + \frac{\boldsymbol{\nabla} T \boldsymbol{\cdot}\boldsymbol{\nabla} T_a}{T_a^2} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \frac{ \boldsymbol{\nabla} T }{T_a } \right), \end{equation}

The time and space average of the second term on the right-hand side can be evaluated as follows:

(5.2)\begin{align} \left\langle \frac{\boldsymbol{\nabla} T \boldsymbol{\cdot} \boldsymbol{\nabla} T_a}{T_a^2} \right\rangle &={-} \int_0^1 \frac{\mathcal{D}}{T_a^2} \mathrm{d}_z \bar{T} \,\mathrm{d} z = \left[ - \frac{\mathcal{D}}{T_a^2} \bar{T} \right]_0^1 + \int_0^1 \mathrm{d}_z \left( \frac{\mathcal{D}}{T_a^2} \right) \bar{T} \,\mathrm{d} z \nonumber\\ &= \frac{4 \mathcal{D} (4 + \mathcal{D}^2)}{(4 - \mathcal{D}^2)^2} + 2 \mathcal{D}^2 \int_0^1 \frac{1}{T_a^2} \frac{\bar{T} }{T_a} \mathrm{d} z . \end{align}

Once again, the maximum principle applying to the entropy variable $s=T/T_a$ is used to bound the second term.

(5.3)\begin{align} \left\langle \frac{\boldsymbol{\nabla} T \boldsymbol{\cdot} \boldsymbol{\nabla} T_a}{T_a^2} \right\rangle & \geqslant \frac{4 \mathcal{D} (4 + \mathcal{D}^2)}{(4 - \mathcal{D}^2)^2} - \frac{2 \mathcal{D}^2}{2 - \mathcal{D}} \int_0^1 \frac{1}{T_a^2} \mathrm{d} z \nonumber\\ & \geqslant \frac{4 \mathcal{D} \left( 4 - 4 \mathcal{D} - \mathcal{D}^2 \right) }{(4 - \mathcal{D}^2 )^2 } . \end{align}

With that bound (5.3), integrating (5.1) in space and time leads to an upper bound

(5.4)\begin{equation} \frac{\mathcal{D}}{Ra } \left\langle \frac{\dot{\varepsilon} : \tau }{T_a} \right\rangle \leqslant \frac{4 \mathcal{D} }{4 - \mathcal{D}^2 } \textit{Nu} - \frac{4 \mathcal{D} \left( 4 - 4 \mathcal{D} - \mathcal{D}^2 \right) }{(4 - \mathcal{D}^2 )^2 } . \end{equation}

As $T_a \leqslant 1 + \mathcal {D}/2$, this becomes an upper bound on viscous dissipation

(5.5)\begin{equation} \left\langle \dot{\varepsilon} : \tau \right\rangle \leqslant \frac{2 \, Ra }{2 - \mathcal{D} } \textit{Nu} \left[ 1 - \frac{ 4 - 4 \mathcal{D} - \mathcal{D}^2 }{(4 - \mathcal{D}^2 ) \textit{Nu} } \right] . \end{equation}

In the limit of large Nusselt numbers (more precisely $(2- \mathcal {D}) \textit {Nu} \gg 1$), this upper bound becomes

(5.6)\begin{equation} \left\langle \dot{\varepsilon} : \tau \right\rangle \leqslant \frac{2 \, Ra }{2 - \mathcal{D} } \textit{Nu} . \end{equation}

6. Obtaining an upper bound on the heat flux

Combining the upper bound (5.5) and the lower bound (4.27) leads to the following inequality:

(6.1)\begin{align} Ra &\geqslant \frac{\textit{Nu}^2 (2 - \mathcal{D})^5 }{12^3 \left( \ln \dfrac{3}{2} \right)^2 \rho_{a0}^2 (2 + 5 \mathcal{D}) } \nonumber\\ &\quad \times\frac{\left[ 1 - \dfrac{\mathcal{D}}{\textit{Nu} (2 + \mathcal{D})} \right]^4 \left[ 1 - \dfrac{24 \left( \ln \dfrac{3}{2} \right)^2 \mathcal{D}^2}{(2-\mathcal{D})^2 \left( \textit{Nu} - \dfrac{\mathcal{D}}{2 + \mathcal{D}} \right)^2 (\gamma -1 )^2 } \right]}{\left[1 + \dfrac{12 \mathcal{D}}{\textit{Nu} (2 - \mathcal{D} ) (2 + 5 \mathcal{D})}\right] \left[ 1 - \dfrac{ 4 - 4 \mathcal{D} - \mathcal{D}^2 }{(4 - \mathcal{D}^2 ) \textit{Nu} } \right] } . \end{align}

This bound is valid when the condition (4.11) is valid. Rigorously, expression (6.1) provides an implicit Nusselt bound in terms of $Ra$ and $\mathcal {D}$ and we have thus proven that the Nusselt number is bounded by the maximum of this implicit bound and the right-hand side of (4.11). With the condition (4.11), we can bound the last part in the right-hand side of (6.1), the square brackets altogether, to be less than $2^{1/2} (2+5 \mathcal {D} )^{-1/2}$ for all acceptable values of $\textit {Nu}$. Moreover, it can be checked that $\rho _{a0} \geqslant 1 + \frac {3}{4} \mathcal {D}$, so that (6.1) leads to the following bound of Nusselt:

(6.2) \begin{equation} \textit{Nu} \leqslant \frac{12^{{3}/{2}}}{2^{{1}/{4}}} \left(\ln \frac{3}{2} \right) \left(1 + \frac{3}{4} \mathcal{D} \right) \left(2 + 5 \mathcal{D} \right)^{{3}/{4}} \frac{Ra^{{1}/{2}}}{(2 - \mathcal{D})^{{5}/{2}}}, \end{equation}

also valid when (4.11) applies.

In the limit $(2- \mathcal {D}) \textit {Nu} \gg 1$, coming back to (6.1), we just have

(6.3)\begin{equation} Ra \geqslant \frac{\textit{Nu}^2 (2 - \mathcal{D})^5 }{12^3 \left( \ln \dfrac{3}{2} \right)^2 \rho_{a0}^2 (2 + 5 \mathcal{D}) } , \end{equation}

which may be rewritten as

(6.4)\begin{equation} \textit{Nu} \leqslant 12^{3/2} \left( \ln \frac{3}{2} \right) \rho_{a0} (2 + 5 \mathcal{D})^{{1}/{2}} \frac{Ra^{{1}/{2}}}{(2 - \mathcal{D})^{{5}/{2}}} . \end{equation}

When $\mathcal {D}$ varies from $0$ to $2$, the value of $\rho _{a0}$ varies from $1$ to $2.5$ (perfect monatomic gas) and $2+5\mathcal {D}$ is less than $12$, so that we have a simpler bound of the form

(6.5)\begin{equation} \textit{Nu} < 25.8 \, \frac{Ra^{{1}/{2}}}{\left(1 - \dfrac{\mathcal{D}}{2} \right)^{{5}/{2}}} . \end{equation}

We plot the $Ra^{1/2}$ prefactors of these Nusselt laws in figure 1, using a logarithmic scale since the singularity at $\mathcal {D}=2$ leads to large values of $\textit {Nu}$. Keep in mind that (6.2) is a rigorous bound, valid for all Nusselt numbers satisfying (4.11), whereas (6.4) and (6.5) are valid only in the limit of very large Nusselt numbers.

Figure 1. Prefactor of $Ra^{1/2}$ in the Nusselt bound, in the limit of large values of $(2-\mathcal {D}) \textit {Nu}$ for a perfect gas with $\gamma =5/3$, computed from (6.4). The simpler expression computed from (6.5) is also shown as a function of $\mathcal {D}$. The thin line shows the bound (A12) when $s_0$ is optimised (see Appendix A).

7. Conclusions

We have obtained an upper bound for the heat transfer in a compressible model of convection known as the anelastic liquid approximation. The bound has been expressed in an implicit algebraic form (6.1), valid under the condition (4.11) of a sufficiently large Nusselt number. An explicit, less stringent, bound (6.2) has also been obtained, under the same condition. The bound (6.1) takes the simpler expression (6.4) in the limit of infinite Nusselt numbers. This simpler expression can itself be bounded by $\textit {Nu} < 25.8\, Ra^{1/2} / (1-\mathcal {D}/2)^{5/2}$.

We had to treat separately the case of small Nusselt numbers (4.29) and that of large Nusselt numbers (4.11). However, it turns out that the case of small Nusselt numbers is most likely irrelevant. Using the explicit bound (6.2), we can derive the condition in the $(\mathcal {D}, Ra)$ space when it is equal to the small-Nusselt-number upper bound (4.29):

(7.1)\begin{equation} Ra^{{1}/{2}} = \frac{\dfrac{\mathcal{D} \left( 2- \mathcal{D} \right)^{{5}/{2}}}{2 + \mathcal{D}}+ 24 \left( \ln \dfrac{3}{2} \right) \left( 2- \mathcal{D} \right)^{{3}/{2}}}{ \dfrac{12^{{3}/{2}}}{2^{{1}/{4}}} \left( \ln \dfrac{3}{2} \right) \left( 1+ \dfrac{3}{4} \mathcal{D} \right) \left(2 + 5 \mathcal{D} \right)^{{3}/{4}} }. \end{equation}

This limit corresponds to the thick full line in figure 2, with a lower region A where the small-Nusselt-number bound (4.29) applies and an upper region B where the large-Nusselt-number bound (6.2) is valid. The maximum value of the Rayleigh number along the limit between the small- and large-Nusselt-number zones (see figure 2) is equal to $4/3 \simeq 1.333$ when $\mathcal {D}=0$. This value is well below the linear stability limit at $\mathcal {D}=0$ and also for any value of $\mathcal {D}$, as shown by Alboussière & Ricard (2017). At $\mathcal {D}=0$, the results of nonlinear stability by Joseph (1976) ensure that no subcritical flow can develop below the linear stability threshold $Ra = 27 {\rm \pi}^4 / 4$, however we do not have such a result for $\mathcal {D} \neq 0$. We can only remark that it is highly unlikely that a convective flow can be sustained below the limit (7.1) in zone A, so that the actual Nusselt number is most probably equal to $1$ in zone A (heat transport by conduction only).

Figure 2. Isovalues of the bound of the Nusselt number (indicated by colour level and a few isolines), for small Nusselt numbers (4.29) in zone A and large Nusselt numbers (6.2) in zone B, separated by the thick black line corresponding to (7.1).

In order to obtain the large-Nusselt-number bound, we have introduced an unusual quantity, the logarithm of entropy (shifted with a constant so that it is positive everywhere), and have shown that it obeys an equation similar to that of entropy. Its flux has a conduction term and a convective term and its sources are positive (or at least bounded from below). The difference with entropy is that it is possible to derive an $L^2$ upper bound for the gradients of this quantity while we could not do so for the gradients of entropy.

Importantly, obtaining a bound relies heavily on the existence of a minimum principle for entropy. Although entropy is bounded from below only, this limit enables us to bound the diffusive part of the flux close to the hot boundary and also to bound the sources of entropy (or its logarithm) from below. For this reason, it will be difficult to obtain an upper bound for the general anelastic model, or for the complete set of compressible equations as soon as entropy is also a function of pressure, as we do not know the value of pressure at that boundary. We may impose an average pressure, but pressure fluctuations are not known a priori. Hence, we do not know what is the minimum value of the entropy on the boundary.

Chernyshenko (2022) has shown how different existing approaches to obtaining bounds on global quantities can be derived from each other. Here, we have used the ‘direct method’ and we can foresee clearly that our result can be cast into his ‘auxiliary functional’ method. Concerning the ‘background profile’ method, this is less straightforward, because the associated functionals will not be quadratic in the unknown fields (velocity and entropy), nor even in our modified variables (velocity and $\log$-entropy). Thus, even if we can write the problem in terms of a background profile (far from obvious, according to Chernyshenko 2022), this will not provide results simply by solving linear eigenvalue problems, so it seems difficult to imagine that the background profile method will help very much in improving our bound. However, at this stage, this is just the expression of our feelings and should not be taken as a definite statement.

There is a degree of freedom in the choice of the constant $s_0$ added to entropy in order to make it strictly positive. We have chosen it so that the minimum value of $s+s_0$ is $a/(4-\mathcal {D}^2)$ where $a$ is a constant equal to $8$, see (4.1). That choice affects the final upper bound of Nusselt number and we can choose $a(\mathcal {D} )$, for each value of $\mathcal {D}$, to minimise the upper bound. This is done explicitly in Appendix A and the optimal bound is plotted in figure 1, from (A12) for the optimal value $a(\mathcal {D} )$. Although our choice $a=8$ is not exactly the optimal choice, it is very close to the optimal bound. One cannot improve the final upper bound of Nusselt number by more than 10 % with another choice for $0.1 < \mathcal {D} < 2$, whereas a better choice could lower the upper bound by a factor of 1.4 near $\mathcal {D} = 0$. The constant $a=8$ has the advantage of making the algebra simpler.

There is another degree of freedom concerning the fraction of the conduction term needed to define the thickness $\delta _0$. We have decided to consider when $1/2$ of the flux must be carried by convection, but we could have taken any fraction between $0$ and $1$. This is investigated in Appendix B. Actually, this choice of $1/2$ leads to the best final bound in our case. In the Boussinesq case, Seis (2015) found that a fraction $1/3$ was the optimum, but both the governing equations and the nature of the flux are different (heat flux vs flux of $\log$-entropy).

The bound (6.1) and its approximations at large Nusselt numbers (6.4) or (6.5) are not expected to be very tight. In the limit $\mathcal {D} \rightarrow 0$, where the anelastic model should converge towards the Boussinesq model, we can readily see that it is less tight than the original bound by Howard by a factor of nearly $20$, whereas the bound (A12) obtained in Appendix A for the optimal value $a_{min}$ is still a factor of 14 above (see figure 1). Concerning large values of $\mathcal {D}$, our bound is made very large owing to the divergent factor $(2-\mathcal {D})^{-5/2}$. We think this is due to our inability to track the $\log$-entropy flux near the cold boundary and, more fundamentally, this originates from the lack of an upper bound for entropy (we only have a lower bound). However, such a trend is not observed in numerical calculations; on the contrary, an increase of the dissipation number seems to lead to a decrease of the heat flux (Curbelo et al. 2019).

In Appendix C, we consider a model of turbulent diffusivity, rather than thermal conduction, with the gradient of entropy replacing the gradient of temperature. The analysis is pretty similar, except there is a condition on the vertical profile of the diffusivity (or conductivity) which must be verified otherwise a bound on the Nusselt number cannot be obtained rigorously. When the small-scale heat flux is modelled as $-c \boldsymbol {\nabla } s$ that condition is $\mathrm {d}_z ( c / T_a^2 ) > 0$, and when the Fourier law of conduction is considered $-k \boldsymbol {\nabla } T$, the condition is $\mathrm {d}_z (k / T_a ) > 0$. For instance, the original model of uniform diffusivity from Lantz & Fan (1999) does not verify this condition, hence we cannot guarantee an upper bound for the Nusselt number. In contrast, the model of uniform conductivity of Mizerski & Tobias (2011) (or that of uniform conductivity for the Fourier law) does verify the condition and our bound is proven.

Among the possible extensions of this work to other models, it would be interesting to consider fluids of infinite Prandtl number. In that case and in the Boussinesq model, tight upper bounds on the heat flow have been obtained (Doering, Otto & Reznikoff 2006). In addition, in the Boussinesq limit and when Coriolis forces are taken into account, upper bounds have been derived (Tilgner 2022a,b). Extending these results to a compressible model of convection would be relevant to planetary convection. For that purpose, we need to consider different models of equation of state for condensed matter, idealised (Alboussière et al. 2022) or more realistic (Ricard et al. 2022; Ricard & Alboussière 2023) concerning compressible convection in planetary interiors.