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 Reference Howard1963). Concerning Rayleigh–Bénard convection, the method of Howard (Reference Howard1963) 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 (Reference Doering and Constantin1996) 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 Reference Joseph1976) of the background profile. The problem is finally turned into finding the background profile with the minimum possible dissipation. Kerswell (Reference Kerswell1999) proved that both methods were dual approaches to the same problem of optimisation. A third method was obtained recently by Seis (Reference Seis2015), 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 (Reference Chernyshenko2022) 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 Reference Anufriev, Jones and Soward2005). 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 Reference Lantz and Fan1999) 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 (Reference Seis2015). 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 (Reference Seis2015). 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. Reference Anufriev, Jones and Soward2005):
$$\begin{gather} {\boldsymbol{\nabla}} \boldsymbol{\cdot} (\rho_a \boldsymbol{v}) = 0 , \end{gather}$$
$$\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}$$
$$\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,
\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 Reference Emanuel1998), are the following:
$$\begin{gather} \dot{\varepsilon}_{ij} = \tfrac{1}{2} \left( \partial_i v_j + \partial_j v_i \right), \end{gather}$$
$$\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. Reference Curbelo, Duarte, Alboussière, Dubuffet, Labrosse and Ricard2019):
$$\begin{gather} T_a (z) = 1 - \mathcal{D} \left(z - \tfrac{1}{2} \right) , \end{gather}$$
$$\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,
\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.
$$\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}$$
$$\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}$$
$$\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}$$
$$\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. (Reference Curbelo, Duarte, Alboussière, Dubuffet, Labrosse and Ricard2019) 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:
\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 Reference Picone1929; Nirenberg Reference Nirenberg1953). 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 (Reference Foias, Manley and Temam1987) 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,
\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
\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
\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}$,
\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:
\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)
\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
\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 )$,
\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:
\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
\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:
\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
\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)$,
\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
\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
\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)
\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:
\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 )$:
\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
\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
\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
\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 Reference Landau and Lifshitz1966),
\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
\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$,
\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
\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
\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$,
\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:
\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
\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,
\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
\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:
\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.
\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
\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
\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
\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:
\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:
\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
\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
\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
\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.
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):
\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 (Reference Alboussière and Ricard2017). At 
$\mathcal {D}=0$, the results of nonlinear stability by Joseph (Reference Joseph1976) 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).
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 (Reference Chernyshenko2022) 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 Reference Chernyshenko2022), 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 (Reference Seis2015) 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. Reference Curbelo, Duarte, Alboussière, Dubuffet, Labrosse and Ricard2019).
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 (Reference Lantz and Fan1999) 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 (Reference Mizerski and Tobias2011) (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 Reference Doering, Otto and Reznikoff2006). In addition, in the Boussinesq limit and when Coriolis forces are taken into account, upper bounds have been derived (Tilgner Reference Tilgner2022a,Reference Tilgnerb). 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. Reference Alboussière, Curbelo, Dubuffet, Labrosse and Ricard2022) or more realistic (Ricard et al. Reference Ricard, Alboussière, Labrosse, Curbelo and Dubuffet2022; Ricard & Alboussière Reference Ricard and Alboussière2023) concerning compressible convection in planetary interiors.
Acknowledgements
The authors thank the referees for useful suggestions that improved the paper. They also thank the ANR for consistently not funding their experimental project of compressible convection, leaving them more time to work on theoretical aspects of the subject.
Declaration of interests
The authors report no conflict of interest.
Appendix A. On the choice of the entropy offset $s_0$
In § 4, we have chosen to add the constant 
$s_0$ to entropy
\begin{equation} s_0 = \frac{1}{2- \mathcal{D}} + \frac{8}{4 - \mathcal{D}^2}. \end{equation}
The choice of the first part 
$1/(2- \mathcal {D})$ is obvious as it ensures that 
$s+s_0$ is positive everywhere and always, by virtue of the minimum principle. The additional positive term 
$8/(4 - \mathcal {D}^2)$ is arbitrary and we investigate in this appendix whether another choice could improve the bound on the Nusselt number. To this end, let us consider a constant 
$s_0$ of the form
\begin{equation} s_0 = \frac{1}{2- \mathcal{D}} + \frac{a}{4 - \mathcal{D}^2}, \end{equation}
where 
$a$ is a strictly positive real number, which we shall make a function of 
$\mathcal {D}$ when needed. Equation (4.6) now becomes
\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{4+a}{a} \right) . \end{equation}Equation (4.7) becomes
\begin{equation} {-}\left[ \mathrm{d}_z \bar{\mathcal{L}} \right]_0^1 = 2 \,\textit{Nu} \left[ \frac{2 + \mathcal{D}}{a} - \frac{2 - \mathcal{D}}{4+a} \right] + 2 \mathcal{D} \left[ \frac{2 + \mathcal{D} }{(2 - \mathcal{D}) a} + \frac{2 - \mathcal{D}}{(2 + \mathcal{D}) (4+a)} \right] . \end{equation}
The corresponding bound for the square of gradients of 
$\mathcal {L}$ obtained in (4.8) now becomes
\begin{align} \left\langle \left| {\boldsymbol{\nabla}} \mathcal{L} \right| ^2 \right\rangle & \leqslant 2 \textit{Nu} \left[ \frac{2 + \mathcal{D}}{a} - \frac{2 - \mathcal{D}}{4+a} \right] + 2 \mathcal{D} \left[ \frac{2 + \mathcal{D} }{(2 - \mathcal{D}) a} + \frac{2 - \mathcal{D}}{(2 + \mathcal{D}) (4+a)} \right] \nonumber\\ &\quad - \frac{4 \mathcal{D}}{2 + \mathcal{D}} \ln \left( \frac{4+a}{a} \right) , \end{align}and the corresponding equation to (4.9) becomes
\begin{equation} \left\langle \left| {\boldsymbol{\nabla}} \mathcal{L} \right| ^2 \right\rangle \leqslant 2 \textit{Nu} \left[ \frac{2 + \mathcal{D}}{a} - \frac{2 - \mathcal{D}}{4+a} \right] + \frac{8 \mathcal{D}}{(2-\mathcal{D}) a}. \end{equation}
Let us now consider the flux of 
$\mathcal {L}$ at the bottom
\begin{equation} \varPhi_{\mathcal{L}} (0) ={-} \mathrm{d}_z \bar{ \mathcal{L} } (0) = \frac{2 (2 - \mathcal{D}) }{4+a} \left( \textit{Nu} - \frac{\mathcal{D}}{2 + \mathcal{D}} \right) . \end{equation}
This provides a maximal value of 
$\bar {\mathcal {L}}$ at 
$z=\delta$
\begin{equation} \bar{\mathcal{L}} (\delta) \leqslant \ln \left( \frac{a+4}{4 - \mathcal{D}^2} \right) - \frac{\delta}{2} \varPhi_{\mathcal{L} } (0). \end{equation}
In order to avoid values below the minimum entropy, 
$\delta$ must be less than 
$\delta _0$, with expression (4.14) changed for
\begin{equation} \delta_0 = \frac{4+a }{(2 - \mathcal{D}) \left( \textit{Nu} - \dfrac{\mathcal{D}}{2 + \mathcal{D}} \right) } \ln \left( \frac{a+4}{a} \right) . \end{equation}
At 
$z=\delta$, the convective flux must be at least half the flux of 
$\mathcal {L}$ at 
$z=0$, so (4.15) becomes
\begin{equation} \rho_a \overline{v_z \mathcal{L}} (\delta) \geqslant \frac{2 - \mathcal{D} }{4+a} \left( \textit{Nu} - \frac{\mathcal{D}}{2 + \mathcal{D}} \right). \end{equation}
Using (4.16), (4.19) and (4.26) which are valid for all values of 
$a$, and using (A6) and (A9), (A10) can be written in the limit of large Nusselt numbers as
\begin{equation} \left\langle \dot{\varepsilon} : \tau \right\rangle \geqslant \frac{a \, \textit{Nu}^3 (2 - \mathcal{D})^4 }{(a+4)^3 \left( \ln \dfrac{a+4}{a} \right)^2 \rho_{a0}^2 (8 + 4 \mathcal{D} + 2 \mathcal{D} a) }, \end{equation}
which is the equivalent of (4.28) for a general value of 
$a$ instead of 
$a=8$.
The upper bound (5.6) from the classical entropy equation is still valid. Combined with (A11), this leads to the following bound for the Nusselt number:
\begin{equation} \textit{Nu} \leqslant (a+4)^{{3}/{2}} Ra^{{1}/{2}} \frac{\rho_{a0}}{(2-\mathcal{D})^{{5}/{2}}} \ln \left( \frac{a+4}{a} \right) \left( \frac{16}{a} + \frac{8 \mathcal{D}}{a} + 4 \mathcal{D} \right)^{{1}/{2}}, \end{equation}
which is identical to (6.4) when 
$a=8$. Now we can optimise 
$a$ to lower the bound (A12). This is done by finding the value of 
$a$ providing the minimum of the following prefactor function 
$f(a,\mathcal {D})$ for each 
$\mathcal {D}$:
\begin{equation} f(a,\mathcal{D}) = (a+4)^{{3}/{2}} \ln \left( \frac{a+4}{a} \right) \left( \frac{4}{a} + \frac{2 \mathcal{D}}{a} + \mathcal{D} \right)^{{1}/{2}}. \end{equation}
The numerical solution to this problem is shown in figure 3. This value is close to 
$8$ for large values of 
$\mathcal {D}$ (actually reaching 7.75 near 
$\mathcal {D}=2$) and diverges to infinity near the Boussinesq limit 
$\mathcal {D}=0$. The gain in term of the bound of Nusselt is defined as the ratio of the prefactor with optimal 
$a$ compared with the prefactor with 
$a=8$. It is shown in figure 3(b). It is close to 1 except for small values of 
$\mathcal {D}$, falling just below 
$0.7$ near 
$\mathcal {D}=0$. The bound (A12) with the optimal value 
$a_{min}$ is shown in figure 1 in the main text.
(a) Value 
$a_{min}$ such that the prefactor of the Nusselt bound (A12) is minimal. (b) Reduction of the prefactor of Nusselt's bound by taking the value 
$a_{min}$ compared with 
$a=8$.
Appendix B. On the definition of the boundary layer thickness $\delta _0$
Another arbitrary choice we have made has been to consider where the convective flux of 
$\bar {\mathcal {L}}$ must be at least one-half of the flux 
$\varPhi _{\mathcal {L}} (0)$ at the bottom (
$z=0$). Instead, let us consider the point where the conduction flux of 
$\bar {\mathcal {L}}$ falls below a factor 
$b$ of the bottom flux (with 
$0 < b < 1$) and the convective flux must be more than 
$(1-b)$ times the flux at the bottom. This changes (4.13) into
\begin{equation} \bar{\mathcal{L}} (\delta) \leqslant \ln \left( \frac{12}{4 - \mathcal{D}^2} \right) - b \delta \varPhi_{\mathcal{L} } (0), \end{equation}and the expression for (4.14) becomes
\begin{equation} \delta_0 = \frac{6 \ln \left( \dfrac{3}{2} \right)}{b \ (2 - \mathcal{D}) \left( \textit{Nu} - \dfrac{\mathcal{D}}{2 + \mathcal{D}} \right) }. \end{equation}Then (4.15) is changed for
\begin{equation} \rho_a \overline{v_z \mathcal{L}} (\delta) \geqslant (1-b) \frac{2 - \mathcal{D}}{6} \left( \textit{Nu} - \frac{\mathcal{D}}{2 + \mathcal{D}} \right) . \end{equation}Following the same steps as in the main text, in the limit of large Nusselt numbers, the lower bound for dissipation (4.28) becomes
\begin{equation} \left\langle \dot{\varepsilon} : \tau \right\rangle \geqslant \frac{4 (1-b)^2b^2 \, \textit{Nu}^3 (2 - \mathcal{D})^4 }{6^3 \left( \ln \dfrac{3}{2} \right)^2 \rho_{a0}^2 (2 + 5 \mathcal{D}) }. \end{equation}
The lower bound is now proportional to 
$b^2 (1-b)^2$. The largest lower bound is then obtained for 
$b=1/2$. The choice we made is indeed optimal. Because the upper bound (5.6) is independent of any choice of 
$\delta _0$, the final bound on the Nusselt number (6.4) is also optimal with respect to the parameter 
$b$.
Appendix C. Application to models of turbulent diffusivity
The idea stems from the model of Lantz & Fan (Reference Lantz and Fan1999), where small scale turbulence transports entropy with a turbulent diffusivity and superadiabatic conduction heat transport 
$-{\boldsymbol {\nabla } T}$ is changed for 
$-\rho _a T_a {\boldsymbol {\nabla }} s$. In this original paper, turbulent diffusivity was taken uniform and used implicitly to determine a Rayleigh number much smaller than it would be evaluated with a molecular diffusivity. However, we can make it more general and consider a superadiabatic heat transfer of the form 
$- c (z) {\boldsymbol {\nabla }} s$. When 
$c (z) = \rho _a T_a$ this corresponds to the most widely used model of Lantz & Fan (Reference Lantz and Fan1999) mimicking a uniform diffusivity and when 
$c (z) = T_a$ this corresponds to a uniform turbulent thermal conductivity (rather than diffusivity; see Mizerski & Tobias Reference Mizerski and Tobias2011). Turbulent diffusivity or conductivity (at small scale) are certainly functions of height, in stars and gas giants, and can be tailored to match detailed calculations or a scaling law. Any choice is possible for 
$c (z)$, as long as it is strictly positive, but we show that we can obtain a Nusselt upper bound for some profiles of 
$c$ (e.g. uniform conductivity), whereas we cannot for other profiles (including the case of uniform diffusivity). The condition when a bound is obtained is uncovered in the following.
Of course we need to change the definition of the Nusselt number 
$\textit {Nu}=- c (z=0) \mathrm {d}_z \bar {s} (z=0)$. Equation (2.3) is changed for
\begin{equation} \rho_a T_a \frac{\mathrm{D} s}{\mathrm{D} t} = \frac{\mathcal{D}}{Ra } \dot{\varepsilon} : \tau + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( c {\boldsymbol{\nabla}} s \right) . \end{equation}
None of the results concerning the minimum principle are affected, the same variable 
$\mathcal {L}$ is considered, with the same constant 
$s_0$ as in § 4. The change in the diffusion term leads to the following equation for 
$\mathcal {L}$ instead of (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 )} + \frac{c}{T_a} \left| {\boldsymbol{\nabla}} \mathcal{L} \right| ^2 + \frac{c}{T_a} \mathrm{d}_z ( \ln T_a ) \partial_z \mathcal{L} + {\boldsymbol{\nabla}} \boldsymbol{\cdot} \left( \frac{c}{T_a} {\boldsymbol{\nabla}} \mathcal{L} \right) . \end{equation}
We also need to change the definition of the flux of 
$\mathcal {L}$,
\begin{equation} \varPhi_{\mathcal{L}} (z) ={-} \frac{c}{T_a} \mathrm{d}_z \bar{\mathcal{L}} + \rho_a \overline{v_z \mathcal{L}}, \end{equation}instead of (4.4) and the equation for the averaged vertical flux (4.3) becomes
\begin{equation} \mathrm{d}_z \varPhi_{\mathcal{L}} = \frac{\mathcal{D}}{Ra \,T_a} \overline{\frac{ \dot{\varepsilon} : \tau }{s +s_0 } } + \frac{c}{T_a} \overline{ \left| {\boldsymbol{\nabla}} \mathcal{L} \right| ^2 } + \frac{c}{T_a} \mathrm{d}_z ( \ln T_a ) \mathrm{d}_z \bar{\mathcal{L}} . \end{equation}
Integrating from 
$z=0$ to 
$z=1$ will provide a bound on 
$\langle | {\boldsymbol {\nabla }} \mathcal {L} | ^2 \rangle$. The integral of the last term in (C4) is
\begin{align} \int_0^1 \frac{c}{T_a} \mathrm{d}_z ( \ln T_a ) \mathrm{d}_z \bar{\mathcal{L}} \,\mathrm{d} z &={-} \int_0^1 \frac{c \mathcal{D}}{T_a^2} \mathrm{d}_z \bar{\mathcal{L}} \,\mathrm{d} z, \nonumber\\ &={-} \left[ \frac{c \mathcal{D} \bar{\mathcal{L}}}{T_a^2} \right]_0^1 + \mathcal{D} \int_0^1 \mathrm{d}_z\left( \frac{c}{T_a^2} \right) \bar{\mathcal{L}} \,\mathrm{d} z . \end{align}
Now comes the condition on 
$c (z)$. If 
$\mathrm {d}_z( c / T_a^2 ) < 0$ for some value of 
$z$, then we cannot have a lower bound of this term since 
$\bar {\mathcal {L}}$ has no upper bound. In the other case, if 
$\mathrm {d}_z( c / T_a^2 ) > 0$, we have
\begin{equation} \int_0^1 \frac{c}{T_a} \mathrm{d}_z ( \ln T_a ) \,\mathrm{d}_z \bar{\mathcal{L}} \mathrm{d} z \geqslant \frac{\mathcal{D}}{\left( 1+ \dfrac{\mathcal{D}}{2} \right)^2} c (0) \ln \left(\frac{3}{2} \right). \end{equation}The integral of the left-hand side of (C4) is evaluated exactly,
\begin{equation} \left[ \varPhi_{\mathcal{L}} \right]_0^1 = \frac{\textit{Nu}}{12} \left( 2 + 5 \mathcal{D} \right), \end{equation}instead of (4.7). Integrating (C4) leads finally to
\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}}{\left( 1+ \dfrac{\mathcal{D}}{2} \right)^2} c (0) \ln \left(\frac{3}{2} \right) . \end{equation}The reasoning is then similar to the anelastic liquid case, with (4.10) becoming
\begin{equation} \varPhi_{\mathcal{L}} (0) ={-} \frac{c (0)}{T_a (0)} \mathrm{d}_z \overline{ \mathcal{L} } (0) = \frac{2 - \mathcal{D}}{6} \textit{Nu} . \end{equation}Expression (4.13) is still valid and (4.14) becomes
\begin{equation} \delta_0 = \frac{12 c \left(\dfrac{1}{2} \right) \ln \left( \dfrac{3}{2} \right) \left(1+ \dfrac{\mathcal{D}}{2} \right)}{(2 - \mathcal{D}) \textit{Nu} } , \end{equation}
where the condition 
$\mathrm {d}_z (c / T_a^2) > 0$ has been used, along with the restriction 
$\delta _0 < 1/2$. The condition on the convective flux (4.15) now becomes
\begin{equation} \rho_a \overline{v_z \mathcal{L}} (\delta) \geqslant \frac{2 - \mathcal{D}}{12} \textit{Nu} . \end{equation}All equations from (4.16) to (4.26) are equally valid for the model with turbulent diffusivity. Equation (C11) leads to the following mirror equation to (4.27), when using (C8), (C10), (4.16), (4.19) and (4.26):
\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 c^2 \left(\dfrac{1}{2}\right) \left(1\,{+}\, \dfrac{\mathcal{D}}{2}
\right)^2 (2 \,{+}\, 5 \mathcal{D}) } \frac{\left[ 1 \,{-}\, \dfrac{24
c^2\left(\dfrac{1}{2}\right) \left( \ln \dfrac{3}{2} \right)^2 \mathcal{D}^2}{(2\,{-}\,\mathcal{D})^2 \textit{Nu}^2 (\gamma -1
)^2 } \right]}{\left[ 1 \,{-}\, \dfrac{12 \mathcal{D} c (0) \ln
\left( \dfrac{3}{2} \right) }{(1+\mathcal{D}/2)^2 \textit{Nu} (2+5 \mathcal{D}) } \right] }. \end{align}
We now derive a lower bound for dissipation, similarly as in § 5. Dividing (C1) by 
$T_a$ and integrating by parts leads to
\begin{equation} \rho_a \frac{\mathrm{D} s}{\mathrm{D} t} = \frac{\mathcal{D}}{Ra } \frac{\dot{\varepsilon} : \tau }{T_a} + c \frac{{\boldsymbol{\nabla}} T_a \boldsymbol{\cdot} {\boldsymbol{\nabla}} s}{T_a^2} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \frac{c}{T_a} {\boldsymbol{\nabla}} s \right) . \end{equation}We aim to bound the integral of the middle term on the right-hand side:
\begin{align} \left\langle \frac{{\boldsymbol{\nabla}} T_a \boldsymbol{\cdot} {\boldsymbol{\nabla}} s}{T_a^2} \right\rangle &= \int_0^1 - \frac{ c \mathcal{D}}{T_a^2} \mathrm{d}_z \bar{s} \,\mathrm{d} z , \end{align}
\begin{align} &={-} \mathcal{D} \int_0^1 \mathrm{d}_z \left( \frac{ c }{T_a^2} \bar{s} \right) \mathrm{d} z + \mathcal{D} \int_0^1 \mathrm{d}_z \left( \frac{ c }{T_a^2} \right) \bar{s} \,\mathrm{d} z . \end{align}
The first term on the right-hand side is evaluated exactly from the boundary conditions, and the last term can only be bounded if the same condition as before 
$\mathrm {d}_z ( c / T_a^2 ) > 0$ is met, since we have no upper bound for 
$\bar {s}$ but only a lower bound 
$-1/ (2 - \mathcal {D})$. Assuming this condition holds, we have
\begin{align} \left\langle \frac{{\boldsymbol{\nabla}} T_a \boldsymbol{\cdot} {\boldsymbol{\nabla}} s}{T_a^2} \right\rangle & \geqslant - \mathcal{D} \left[ \frac{ c }{T_a^2} \bar{s} \right]_0^1 -\frac{\mathcal{D} }{2 - \mathcal{D} } \left[ \frac{ c }{T_a^2} \right]_0^1 , \end{align}
\begin{align} & \geqslant \frac{4 \mathcal{D} c (0)}{\left( 4 - \mathcal{D}^2 \right) \left( 1 + \dfrac{\mathcal{D}}{2} \right)^2 } . \end{align}The integral of (C13) provides then an upper bound for dissipation
\begin{equation} \left\langle \dot{\varepsilon} : \tau \right\rangle \leqslant \frac{ 2 Ra }{2- \mathcal{D}} \textit{Nu} \left[ 1 - \frac{ c (0) }{\left( 1 + \dfrac{\mathcal{D}}{2} \right)^2 \textit{Nu} } \right] . \end{equation}Combining (C13) and (C18) lead to an upper bound on the Nusselt number,
\begin{align} \textit{Nu} & \leqslant
\frac{24 \sqrt{3} Ra^{1/2} \rho_{a0} \ln\left( \dfrac{3}{2}
\right) c \left(\dfrac{1}{2} \right) \left( 1 +
\dfrac{\mathcal{D}}{2} \right) \sqrt{2 + 5 \mathcal{D}}
}{(2 -\mathcal{D})^{5/2} } \nonumber\\ & \quad\times
\frac{\left[1 - \dfrac{ c (0) }{\textit{Nu} \left( 1 +
\dfrac{\mathcal{D}}{2} \right)^2 } \right]^{1/2} \left[ 1 -
\dfrac{12 \mathcal{D} c (0) \ln \left( \dfrac{3}{2} \right)
}{(1+\mathcal{D}/2)^2 \textit{Nu} (2+5 \mathcal{D}) }
\right]^{1/2} }{\left[ 1 - \dfrac{24 c^2\left(\dfrac{1}{2}\right) \left( \ln
\dfrac{3}{2} \right)^2 \mathcal{D}^2}{(2-\mathcal{D})^2
\textit{Nu}^2 (\gamma -1 )^2 } \right]^{1/2}},
\end{align}which for very large values of the Nusselt number becomes
\begin{equation} \textit{Nu} \leqslant \frac{24 \sqrt{3} Ra^{{1}/{2}} \rho_{a0} \ln\left( \dfrac{3}{2} \right) c\left(\dfrac{1}{2} \right) \left( 1 + \dfrac{\mathcal{D}}{2} \right) \sqrt{2 + 5 \mathcal{D}} }{(2 -\mathcal{D})^{{5}/{2}} } . \end{equation}
This expression is rather similar to that obtained with classical uniform thermal conduction (6.3). Let us recall that the important message is that we have to make this assumption 
$\mathrm {d}_z (c / T_a^2 )> 0$ in order to ascertain the validity of this bound. This is not the case of the original model of Lantz & Fan (Reference Lantz and Fan1999), with 
$c = \rho _a T_a$: in that case, 
$c/T_a^2 = \rho _a / T_a$, which from (2.7) and (2.8) is proportional to 
$T_a^{({2-\gamma })/({\gamma -1})}$, hence decreasing with 
$z$ increasing, like 
$T_a$, since 
$1 < \gamma < 2$. In contrast, a model mimicking a uniform thermal conductivity is 
$c = T_a$, see Mizerski & Tobias (Reference Mizerski and Tobias2011). In that case, 
$c / T_a^2 = 1/T_a$ is a function increasing with 
$z$ increasing, and the bound (C20) is valid.
Coming back to the main text, we have considered the usual Fourier law of conduction heat transfer 
$-k \boldsymbol {\nabla } T$ with a uniform conductivity 
$k=1$. If we make that conductivity a function of height 
$z$, we can see that a similar condition holds on 
$k(z)$. In (4.5), when integrating by parts, the function 
$\mathrm {d}_z (k / T_a)$ will appear and its sign is crucial when using the minimum principle for entropy, in order to obtain the lower bound (4.6). The same condition must hold to derive (5.3) from (5.2).
Open access
