1. Introduction
In this paper, we consider a layer of fluid trapped between two parallel horizontal plates held at different temperatures. The dimensionless equations of motions for the Boussinesq approximation are
\begin{gather}
\frac{1}{{\textit{Pr}}}\,(\partial_t \boldsymbol{u} +
\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}
\boldsymbol{u}) - \Delta \boldsymbol{u}
+\boldsymbol{\nabla} p = {\textit{Ra}}\,T\,\boldsymbol{e}_2,
\end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}
\begin{gather}\partial_t T + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} T -\Delta T =0, \end{gather}
where the Rayleigh number 
${\textit {Ra}}$ is defined as
\begin{equation} {\textit{Ra}}=\frac{g\alpha\,\delta T\,h^3}{\kappa \nu}, \end{equation}
and the Prandtl number 
${\textit {Pr}}$ is
\begin{equation} {\textit{Pr}}=\frac{\nu}{\kappa}. \end{equation}
In these definitions, 
$g$ is the gravitational constant, 
$\alpha$ is the thermal expansion coefficient, 
$\nu$ is the kinematic viscosity, 
$\kappa$ is the thermal diffusivity, 
$h$ is the distance between the plates, and 
$\delta T =T_{{bottom}}-T_{{top}}$ is the temperature gap.
Lengths are measured in units of 
$h$, time in units of 
${h^2}/{\kappa }$, and temperature in units of 
$\delta T$. In the rectangular domain 
$\varOmega =[0,\varGamma ]\times [0,1]$ the velocity 
$\boldsymbol {u}=u_1(\boldsymbol {x},t)\,\boldsymbol {e}_1+u_2(\boldsymbol {x},t)\,\boldsymbol {e}_2$ and temperature 
$T=T(\boldsymbol {x},t)$ are initialized at 
$t=0$, where
\begin{gather} \boldsymbol{u}(\boldsymbol{x},0) = \boldsymbol{u}_0(\boldsymbol{x}), \end{gather}
\begin{gather}T(\boldsymbol{x},0) = T_0(\boldsymbol{x}). \end{gather}The boundary conditions for the temperature are
\begin{gather} T=0 \quad \text{at } x_2=1, \end{gather}
\begin{gather}T=1 \quad \text{at } x_2=0, \end{gather}while we assume Navier-slip boundary conditions for the velocity field, i.e.
\begin{gather} u_2=0, \partial_2 u_1={-}\frac{1}{L_s}\,u_1\quad \text{at } x_2=1, \end{gather}
\begin{gather}u_2=0, \partial_2 u_1=\frac{1}{L_s}\,u_1\quad \text{at } x_2=0, \end{gather}
where 
$L_s$ is the constant slip length. Here, 
$x_2=\boldsymbol {x}\boldsymbol {\cdot } \boldsymbol {e}_2$. In the horizontal variable 
$x_1=\boldsymbol {x}\boldsymbol {\cdot } \boldsymbol {e}_1$, all variables, including the pressure 
$p=p(\boldsymbol {x},t)$, are periodic.
We are interested in quantifying the heat transport in the upward direction as measured by the non-dimensional Nusselt number
\begin{equation} {\textit{Nu}}=\left\langle \int_0^1(u_2T-\partial_2T)\,{\rm d}\kern 0.06em x_2\right\rangle, \end{equation}where
\begin{equation} \langle \cdot\rangle=\limsup_{t\rightarrow \infty}\frac 1t\int_0^t\frac{1}{\varGamma}\int_0^{\varGamma}(\cdot)\, {\rm d}\kern 0.06em x_1\, {\rm d}s. \end{equation}This number, of utmost relevance in geophysics and industrial applications (Plumley & Julien Reference Plumley and Julien2019), is predicted to obey a power-law scaling of the type
\begin{equation} {\textit{Nu}}\sim {\textit{Ra}}^{\alpha}\,{\textit{Pr}}^{\beta}. \end{equation}
Although physical arguments suggest certain scaling exponents 
$\alpha$ and 
$\beta$ (Malkus Reference Malkus1954; Kraichnan Reference Kraichnan1962; Spiegel Reference Spiegel1971; Siggia Reference Siggia1994), where transitions between scalings could occur (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009), these theories need to be validated. While experiments are expensive and difficult (Ahlers Reference Ahlers2006), numerical studies are limited by the lack of computational power in reaching high-Rayleigh-number regimes (Plumley & Julien Reference Plumley and Julien2019). Even in two spatial dimensions, there have been recent debates about the presence/absence of evidence of the ‘ultimate scaling’ (scaling that holds in the regime of very large Rayleigh numbers) for the Nusselt number (Zhu et al. Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018, Reference Zhu, Mathai, Stevens, Verzicco and Lohse2019; Doering, Toppaladoddi & Wettlaufer Reference Doering, Toppaladoddi and Wettlaufer2019; Doering Reference Doering2020). We use mathematical analysis in order to derive universal upper bounds for the Nusselt number, which serve as a rigorous indication for the scaling exponents holding in the turbulent regime 
${\textit {Ra}}\rightarrow \infty$. The properties of boundary layers and their thicknesses play a central role in the scaling laws for the Nusselt number in Rayleigh–Bénard convection (see Nobili (Reference Nobili2023) and references therein). For this reason, it is interesting to study how heat transport properties change when varying the boundary conditions. In particular, we ask the following question: are there boundary conditions inhibiting or enhancing heat transport compared to the classical no-slip boundary conditions? While many theoretical studies focus on no-slip boundary conditions (Doering & Constantin Reference Doering and Constantin1996, Reference Doering and Constantin2001; Doering, Otto & Reznikoff Reference Doering, Otto and Reznikoff2006; Otto & Seis Reference Otto and Seis2011; Goluskin & Doering Reference Goluskin and Doering2016; Tobasco & Doering Reference Tobasco and Doering2017), other reasonable boundary conditions have been far less explored. In the early 2000s, Ierley, Kerswell & Plasting (Reference Ierley, Kerswell and Plasting2006) considered the Rayleigh–Bénard convection problem at infinite Prandtl number and free-slip boundary conditions; their computational result 
${\textit {Nu}}\lesssim {\textit {Ra}}^{{5}/{12}}$ is obtained by combining the Busse asymptotic expansion in multiple boundary layer solutions (multi-
$\alpha$ solutions) and the Constantin and Doering background field method approach (Doering & Constantin Reference Doering and Constantin1994, Reference Doering and Constantin1996). Inspired by this result and a numerical study in the thesis of Otero (Reference Otero2002), Doering and Whitehead rigorously proved
\begin{equation} {\textit{Nu}}\lesssim {\textit{Ra}}^{{5}/{12}} \end{equation}for the two-dimensional finite Prandtl model (Whitehead & Doering Reference Whitehead and Doering2011) and for the three-dimensional infinite Prandtl number model (Whitehead & Doering Reference Whitehead and Doering2012) using an elaborate application of the background field method. By a perturbation argument around Stokes equations, Wang & Whitehead (Reference Wang and Whitehead2013) proved
\begin{equation} {\textit{Nu}}\lesssim {\textit{Ra}}^{{5}/{12}}+{Gr}^2\,{\textit{Ra}}^{1/4} \end{equation}
in three dimensions for small Grashof number (
${Gr}={{\textit {Ra}}}/{{\textit {Pr}}}$).
In many physical situations, the Navier-slip boundary conditions are used to describe the presence of slip at the solid–liquid interface (Neto et al. Reference Neto, Evans, Bonaccurso, Butt and Craig2005; Uthe, Sader & Pelton Reference Uthe, Sader and Pelton2022). We notice that for any finite 
$L_s> 0$, these conditions imply vorticity production at the boundary. In the limit of infinite slip length, the Navier-slip boundary conditions reduce to free-slip boundary conditions, while in the limit 
$L_s\rightarrow 0$, they converge to the no-slip boundary conditions.
Inspired by the seminal paper of Whitehead & Doering (Reference Whitehead and Doering2011), Drivas, Nguyen and the second author of this paper considered the two-dimensional Rayleigh–Bénard convection model with Navier-slip boundary conditions, and rigorously proved the upper bound
\begin{equation} {\textit{Nu}}\lesssim {\textit{Ra}}^{{5}/{12}}+L_s^{{-}2}\,{\textit{Ra}}^{1/2} \end{equation}
when 
${\textit {Pr}}\geq {\textit {Ra}}^{3/4}\,L_s^{-1}$ and 
$L_s\gtrsim 1$ (Drivas, Nguyen & Nobili Reference Drivas, Nguyen and Nobili2022). The authors refer to (1.15) as an interpolation bound: assuming 
$L_s\sim {\textit {Ra}}^{\alpha }$ with 
$\alpha \geq 0$ (the justification of this choice can be found in the appendix of Bleitner & Nobili Reference Bleitner and Nobili2024), it can be deduced easily that
\begin{equation} {\textit{Nu}}\lesssim \begin{cases} {\textit{Ra}}^{{5}/{12}}, & \text{if } \alpha\geq \frac{1}{24},\\ {\textit{Ra}}^{{1}/{2}-2\alpha}, & \text{if } 0\leq\alpha\leq \frac{1}{24}. \end{cases} \end{equation}In Bleitner & Nobili (Reference Bleitner and Nobili2024), the authors generalized the result in Drivas et al. (Reference Drivas, Nguyen and Nobili2022) to the case of rough walls and Navier-slip boundary conditions, proving interpolation bounds exhibiting explicit dependency on the spatially varying friction coefficient and curvature.
Given the model (1.1), (1.2), (1.3) with (1.8) and (1.9), the objective of this paper is twofold: on one hand we want to apply the so-called direct method (as outlined in Otto & Seis Reference Otto and Seis2011) to derive upper bounds on the Nusselt number exhibiting a transparent relation with the thermal boundary layer's thickness. For discussions about different approaches to derive bounds on the Nusselt number, we refer the reader to Chernyshenko (Reference Chernyshenko2022) and references therein. On the other hand, we aim to improve the bound in Drivas et al. (Reference Drivas, Nguyen and Nobili2022) by refining the estimates on the pressure. Our new result is stated in the following theorem.
Theorem 1.1 Suppose 
$\boldsymbol {u}_0\in W^{1,4}$ and 
$0\leq T_0\leq 1$.
If 
$L_s=\infty$ (i.e. 
$\boldsymbol {u}$ satisfies free-slip boundary conditions), then
\begin{equation} {\textit{Nu}} \lesssim {\textit{Ra}}^{{5}/{12}}. \end{equation} If 
$1\leq L_s<\infty$, then
\begin{equation} {\textit{Nu}} \lesssim {\textit{Ra}}^{{5}/{12}}+L_s^{-{1}/{6}}\,{\textit{Pr}}^{-{1}/{6}}\,{\textit{Ra}}^{{1}/{2}}. \end{equation} If 
$0< L_s<1$, then
\begin{equation} {\textit{Nu}} \lesssim L_s^{-{1}/{3}}\,{\textit{Ra}}^{{1}/{3}} + L_s^{-{2}/{3}}\,{\textit{Pr}}^{-{1}/{6}}\,{\textit{Ra}}^{{1}/{2}}+ L_s^{-{2}/{13}}\,{\textit{Ra}}^{{5}/{13}}+{\textit{Ra}}^{{5}/{12}}. \end{equation} We reserve the discussion on physical implications of this result for the conclusions, in § 4. We first remark that for 
$L_s\geq 1$, this result improves the upper bound in (1.15): in fact, if 
$L_s\sim {\textit {Ra}}^{\alpha }$, with 
$\alpha \geq 0$, then our new result yields
\begin{equation} {\textit{Nu}}\lesssim \begin{cases} {\textit{Ra}}^{{5}/{12}}, & \text{if } {\textit{Pr}}\geq {\textit{Ra}}^{1/2-\alpha},\\ {\textit{Pr}}^{{-}1/6}\,{\textit{Ra}}^{{1}/{2}-{\alpha}/{6}}, & \text{if } {\textit{Pr}}\leq{\textit{Ra}}^{1/2-\alpha} . \end{cases} \end{equation}
In particular, when 
$\alpha =0$, we notice a crossover at 
${\textit {Pr}}\sim {\textit {Ra}}^{1/2}$ between the 
$Ra^{1/2}$ and 
${\textit {Pr}}^{-1/6}\,{\textit {Ra}}^{1/2}$ scaling regimes. This is reminiscent of the upper bound in Choffrut, Nobili & Otto (Reference Choffrut, Nobili and Otto2016) for no-slip boundary conditions. Our result also covers the case 
$0< L_s<1$, and in this region, we can detect four scaling regimes depending on the magnitude of the Prandtl number and 
$L_s<1$. The dominating terms in (1.18) and (1.19) are summarized in table 1. Observe that for 
${\textit {Pr}}\rightarrow \infty$, the term 
$L_s^{-1/3}\,{\textit {Ra}}^{1/3}$ is dominating in the region 
$0< L_s<{\textit {Ra}}^{-{2}/{7}}$. On one hand, when 
${\textit {Pr}}\rightarrow \infty$, this seems to indicate the (expected) transition from Navier-slip to no-slip boundary conditions in the bounds. In order to contextualize this remark, we recall that for no-slip boundary conditions (
$L_s=0$), the upper bound 
${\textit {Nu}}\lesssim {\textit {Ra}}^{1/3}$ was proven when 
${\textit {Pr}}=\infty$ (Otto & Seis Reference Otto and Seis2011) and when 
${\textit {Pr}}\gtrsim {\textit {Ra}}^{1/3}$ (Choffrut et al. Reference Choffrut, Nobili and Otto2016). On the other hand, we stress that our bounding method breaks down in the limit 
$L_s\rightarrow 0$, and consequently all but the last scaling prefactors in (1.19) blow up.
Table 1. Overview of the results in Theorem 1.1. The colouring corresponds to the cases 
$L_s\leq 1$ and 
$1\leq L_s$. In all the other (uncoloured) cases, 
$L_s$ may be smaller or larger than 1.
Differently from the result in Whitehead & Doering (Reference Whitehead and Doering2011), the proof of our theorem does not rely on the background field method but rather exploits the regularity properties of the flow through a localization principle. In fact, the Nusselt number can be localized in the vertical variable:
\begin{equation} {\textit{Nu}}=\frac{1}{\delta}\left\langle \int_0^{\delta} (u_2T-\partial_2 T)\, {\rm d}\kern 0.06em x_2\right\rangle\leq \frac{1}{\delta}\left\langle \int_0^{\delta} u_2T\, {\rm d}\kern 0.06em x_2\right\rangle+\frac{1}{\delta}. \end{equation}
Here, we used the boundary conditions for the temperature and the maximum principle 
$\sup _x|T(x,t)|\leq 1$ for all 
$t$. Notice that this localization principle comes from the fact that the (long-time and horizontal average of the) heat flux is the same for each 
$x_2\in (0,1)$, and this can be deduced by merely using the non-penetration boundary conditions 
$u_2=0$ at 
$x_2=0$, 
$x_2=1$ (see the proof of Lemma 2.1 below).
The second crucial point in our proof is the interpolation inequality
\begin{gather} \frac{1}{\delta}\left\langle \int_{0}^{\delta} u_2T\,{\rm d}\kern 0.06em x_2\right\rangle\leq \frac{1}{2}\left\langle \int_0^{1}|\partial_2 T|^2\, {\rm d}\kern 0.06em x_2\right\rangle+C\delta^3\left\langle \int_0^{1}|\omega|^2\, {\rm d}\kern 0.06em x_2\right\rangle^{1/2}\left\langle \int_0^{1}|\partial_1 \omega|^2\, {\rm d}\kern 0.06em x_2\right\rangle^{1/2}, \end{gather}
where 
$C$ is some positive constant, which was first derived (in a slightly different form) in Drivas et al. (Reference Drivas, Nguyen and Nobili2022); it is proved in Lemma 3.1 below. We again observe that this bound holds only relying on the assumption 
$u_2=0$ at the boundaries 
$x_2=0$, 
$x_2=1$.
In Lemma 2.5, we prove the new pressure estimate
\begin{equation} \|p\|_{H^1}\leq C\left( \frac{1}{L_s}\,\|\partial_2 \boldsymbol{u}\|_{L^2}+\frac{1}{{\textit{Pr}}}\,\|\omega\|_{L^2}\,\|\omega\|_{L^r}+{\textit{Ra}}\,\|T\|_{L^2}\right), \end{equation}
where 
$L^2$ is the space of square integrable functions, while 
$H^1$ is the space of functions in 
$L^2$ with square integrable gradients. The space 
$L^r$ instead consists of functions whose 
$r$th power is integrable. See definitions in (1.25). Notice that the improvement (compared to the pressure estimate in Proposition 2.7 in Drivas et al. Reference Drivas, Nguyen and Nobili2022) lies in the first term on the right-hand side, stemming from the new trace-type estimate
\begin{equation} \left|\int_0^{\varGamma}(p\,\partial_1u_1|_{x_2=1}+p\,\partial_1u_1|_{x_2=0})\, {\rm d}\kern 0.06em x_1\right|\leq 3\,\|p\|_{H^1}\,\|\partial_2 \boldsymbol{u}\|_{L^2}, \end{equation}
used to control the boundary terms in the 
$H^1$ pressure identity (2.49). We observe that the new pressure estimate enables us to treat the case of small slip length, i.e. 
$0< L_s\leq 1$. This small slip length regime was not treatable in Drivas et al. (Reference Drivas, Nguyen and Nobili2022); see Remark 2.6.
1.1. Organization and notations
The paper is divided in two sections: in § 2, we prove all the a priori estimates that we will need in order to prove the main theorem in § 3. The crucial localization and interpolation lemmas are proven in Lemmas 2.1 and 3.1, respectively. The improvement of the upper bounds on the Nusselt number stems from the new pressure estimates in Lemma 2.5. In § 4, we contextualize our result and give a physical interpretation of our bounds.
Throughout the paper, we will use the following Lebesgue and Sobolev norms:
\begin{equation} \|\,f\|_{L^p}^p = \int_{\varOmega} |\,f|^p\ {{\rm d} x}, \quad \|\,f\|_{W^{1,p}}^p = \int_{\varOmega} |\,f|^p\ {{\rm d} x} + \int_{\varOmega} |\boldsymbol{\nabla} f|^p\ {{\rm d} x}, \quad \|\,f\|_{H^1} = \|\,f\|_{W^{1,2}}, \end{equation}
for any 
$1\leq p <\infty$.
2. Identities and a priori bounds
In this section, we derive a priori bounds for the Rayleigh–Bénard convection problem with Navier-slip boundary conditions. These will be used in the proof of Theorem 1.1 in § 3.
Recall that the temperature equation enjoys a maximum principle: if 
$0\leq T_0(\boldsymbol {x})\leq 1$, then
\begin{equation} 0\leq T(\boldsymbol{x},t)\leq 1,\quad \text{for all } \boldsymbol{x},t. \end{equation} The following lemma will allow us to localize the Nusselt number in a strip of height 
$\delta >0$, indicating the thermal boundary layer. This is the key ingredient of the direct method, and will be used later to bound the heat transfer.
Lemma 2.1 (Localization of the Nusselt number)
The Nusselt number 
${\textit {Nu}}$ defined in (1.10) is independent of 
$x_2$, i.e.
\begin{equation} {\textit{Nu}}=\langle u_2T-\partial_2 T\rangle,\quad \text{for all } x_2\in [0,1]. \end{equation}
In particular, for any 
$\delta \in (0,1)$, we have
\begin{equation} {\textit{Nu}}=\frac{1}{\delta}\left\langle \int_0^{\delta}(u_2T-\partial_2 T)\, {\rm d}\kern 0.06em x_2\right\rangle. \end{equation}Proof. Taking the long-time and horizontal average of the equation for 
$T$, we have
\begin{equation} 0=\langle \partial_t T\rangle+\langle \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u} T-\boldsymbol{\nabla} T)\rangle= \partial_{2}\langle u_2T-\partial_2 T\rangle. \end{equation}
Hence 
$\langle u_2T-\partial _2 T\rangle$ is constant in 
$x_2$, proving (2.2). The identity (2.3) is a direct consequence of (2.2).
Notice that from the boundary condition for 
$T$ at 
$x_2=0$, i.e. (1.8b), and the maximum principle (2.1), it follows that
\begin{equation} -\int_0^\delta \partial_2 T \, {\rm d}\kern 0.06em x_2={-}T(x_1,\delta)+1\leq 1. \end{equation}As a consequence, we obtain
\begin{equation} {\textit{Nu}} \leq \frac{1}{\delta}\left\langle \int_0^{\delta}u_2T\, {\rm d}\kern 0.06em x_2\right\rangle + \frac{1}{\delta}. \end{equation}Thanks to (2.2), we can now derive another useful identity relating temperature gradients (naturally emerging in (2.46) and (3.1)) to the Nusselt number.
Lemma 2.2 (Representation of the Nusselt number)
The Nusselt number, defined in (1.10), has the following alternative representation:
\begin{equation} {\textit{Nu}}=\left\langle \int_0^{1}|\boldsymbol{\nabla} T|^2\, {\rm d}\kern 0.06em x_2\right\rangle. \end{equation}Proof. Testing (1.3), the temperature equation with 
$T$, and integrating by parts, we obtain
\begin{equation} \frac{1}{2}\,\frac{{\rm d}}{{\rm d}t}\|T\|_{L^2}^2={-}\|\boldsymbol{\nabla} T\|_{L^2}^2-\int_0^{\varGamma}\partial_2 T|_{x_2=0}\, {\rm d}\kern 0.06em x_1, \end{equation}
where we used the incompressibility condition (1.2), the non-penetration condition 
$\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {e}_2=u_2=0$ at 
$x_2=\{0,1\}$, and the boundary conditions (1.8) for 
$T$. The statement follows from taking the long-time averages, observing that 
$\limsup _{t\rightarrow \infty }\int _0^t ({{\rm d}}/{{\rm d}s})({1}/{\varGamma })$ 
$\int _0^{\varGamma }\int _0^1|T|^2\, {\rm d}\kern 0.06em x_2\,{\rm d}\kern 0.06em x_1\, {\rm d}s=0$ thanks to the maximum principle for 
$T$ (2.1), and using 
${\textit {Nu}}=\langle u_2T-\partial _2 T\rangle |_{x_2=0}$ by (2.2).
The subsequent Lemma 2.3 provides a bound on the long-time (and spatial) average of the velocity gradient, naturally arising from the interpolation estimate in Lemma 3.1. Moreover this bound will be used to control the vorticity gradient in Lemma 2.4.
Lemma 2.3 (Energy)
Let 
$0< L_s<\infty$, and suppose that 
$\boldsymbol {u}_0\in L^2$. Then there exists a constant 
$C=C(\varGamma )>0$ such that
\begin{equation} \|\boldsymbol{u}(t)\|_{L^2}\leq \|\boldsymbol{u}_0\|_{L^2} + C\max\{1,L_s\}\,{\textit{Ra}} \end{equation}
for all times 
$t\in [0,\infty )$, and
\begin{equation} \left\langle\int_0^1 |\boldsymbol{\nabla} \boldsymbol{u}|^2\,{\rm d}\kern 0.06em x_2\right\rangle\leq {\textit{Nu}}\,{\textit{Ra}}. \end{equation}
If 
$L_s=\infty$, then the bound (2.9) simplifies to
\begin{equation} \|\boldsymbol{u}(t)\|_{L^2}\leq \|\boldsymbol{u}_0\|_{L^2} + C\,{\textit{Ra}}, \end{equation}and (2.10) remains valid.
Proof. At first, we assume 
$0< L_s<\infty$ and test the velocity equation (1.1) with 
$\boldsymbol {u}$, integrate by parts, and use the boundary conditions (1.9) for 
$u$ and the incompressibility condition (1.2) to find
\begin{align} \frac{1}{2\,{\textit{Pr}}}\,\frac{{\rm d}}{{\rm d}t}\int_{\varOmega} |\boldsymbol{u}|^2\, {{\rm d} x}+ \int_{\varOmega} |\boldsymbol{\nabla} \boldsymbol{u}|^2\, {{\rm d} x}+\frac{1}{L_s}\int_{0}^{\varGamma} (u_1^2|_{x_2=0}+u_1^2|_{x_2=1})\, {\rm d}\kern 0.06em x_1={\textit{Ra}} \int_{\varOmega} Tu_2\, {{\rm d} x}.\end{align}The fundamental theorem of calculus implies
\begin{equation} u_1(\boldsymbol{x}) = u_1(x_1,0) + \int_0^{x_2} \partial_2 u_1(x_1,z)\, {\rm d}z, \end{equation}and using Young's and Hölder's inequality gives
\begin{align} u_1^2(\boldsymbol{x}) &= u_1^2(x_1,0) + 2 u_1(x_1,0)\int_0^{x_2} \partial_2 u_1(x_1,z)\, {\rm d}z + \left(\int_0^{x_2} \partial_2 u_1(x_1,z)\, {\rm d}z\right)^2 \end{align}
\begin{align} &\leq 2 u_1^2(x_1,0) + 2\left(\int_0^{x_2} |\partial_2 u_1(x_1,z)|\, {\rm d}z\right)^2 \end{align}
\begin{align} &\leq 2 u_1^2(x_1,0) + 2\int_0^{x_2} |\partial_2 u_1(x_1,z)|^2\, {\rm d}z, \end{align}which after integration implies
\begin{equation} \|u_1\|_{L^2}^2 \leq 2 \int_0^\varGamma u_1^2\vert_{x_2=0}\, {\rm d}\kern 0.06em x_1 + 2 \|\partial_2 u_1\|_{L^2}^2.\end{equation}
By (1.9), one has 
$u_2=0$ on the boundaries, therefore the analogous estimate shows that
\begin{equation} \|u_2\|_{L^2}^2 \leq \|\partial_2 u_2\|_{L^2}^2. \end{equation}The full vector norm of the velocity can now be split into the norms of its components, and by (2.17) and (2.18), we have
\begin{align} \|\boldsymbol{u}\|_{L^2}^2 &= \|u_1\|_{L^2}^2 + \|u_2\|_{L^2}^2 \end{align}
\begin{align} &\leq 2 \left(\|\partial_2 u_1\|_{L^2}^2 + \int_{0}^{\varGamma} (u_1^2|_{x_2=0}+u_1^2|_{x_2=1})\, {\rm d}\kern 0.06em x_1 + \|\partial_2 u_2\|_{L^2}^2\right) \end{align}
\begin{align} &\leq 2 \left(\|\boldsymbol{\nabla} \boldsymbol{u}\|_{L^2}^2 + \int_{0}^{\varGamma} (u_1^2|_{x_2=0}+u_1^2|_{x_2=1})\, {\rm d}\kern 0.06em x_1\right), \end{align}which applied to (2.12) yields
\begin{align} &\frac
{1}{2\,{\textit{Pr}}}\,\frac{{\rm d}}{{\rm
d}t}\int_{\varOmega} |\boldsymbol{u}|^2\, {{\rm d} x}+
\frac{1}{C}\min\{1,L_s^{{-}1}\}\,\|\boldsymbol{u}\|_{L^2}^2\nonumber\\
&\quad \leq{\textit{Ra}} \int_{\varOmega} Tu_2\, {{\rm
d} x}
\leq{\textit{Ra}}\,\|T\|_{L^2}\,\|u_2\|_{L^2}\nonumber\\ &\quad \leq
\frac{1}{4\epsilon}\,\|T\|_{L^2}^2\,{\textit{Ra}}^2 +
\epsilon\,\|u_2\|_{L^2}^2
\leq
\frac{1}{4\epsilon}\,\varGamma\,{\textit{Ra}}^2 +
\epsilon\,\|u_2\|_{L^2}^2,
\end{align} for some constant 
$C>0$, where we used Hölder's inequality, Young's inequality and 
$\|T\|_{L^{\infty }}\leq 1$ because of the maximum principle (2.1). Setting 
$\epsilon = ({1}/{2C})\min \{1,L_s^{-1}\}$ implies
\begin{equation} \frac {1}{{\textit{Pr}}}\,\frac{{\rm d}}{{\rm d}t}\|\boldsymbol{u}\|_{L^2}^2+ \frac{1}{C}\min\{1,L_s^{{-}1}\}\,\|\boldsymbol{u}\|_{L^2}^2\leq C\max\{1,L_s\}\,\varGamma\,{\textit{Ra}}^2, \end{equation}and Grönwall's inequality now yields (2.9). Taking the long-time average of (2.12), using (2.9), one has
\begin{equation} \left\langle\int_0^1 |\boldsymbol{\nabla} \boldsymbol{u}|^2\, {\rm d}\kern 0.06em x_2\right\rangle+\frac{1}{L_s}\left\langle (u_1^2|_{x_2=0} + u_1^2|_{x_2=1})\right\rangle={\textit{Ra}} \left\langle \int_0^1 Tu_2\, {\rm d}\kern 0.06em x_2\right\rangle. \end{equation}
The claim follows by observing that, due to the boundary conditions for 
$T$, we have
\begin{equation} {\textit{Nu}}=\left\langle \int_0^1 Tu_2\, {\rm d}\kern 0.06em x_2\right\rangle+1. \end{equation}
If 
$L_s=\infty$, then integrating the first component of (1.1) in space yields
\begin{equation} \frac{1}{{\textit{Pr}}}\,\frac{{\rm d}}{{\rm d}t} \int_{\varOmega} u_1 \, {{\rm d} x} ={-} \frac{1}{{\textit{Pr}}}\int_{\varOmega} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}u_1 \, {{\rm d} x} + \int_{\varOmega} \Delta u_1 \, {{\rm d} x} - \int_{\varOmega} \partial_1 p. \end{equation}
The first term on the right-hand side of (2.26) vanishes after integration by parts due to the incompressibility condition (1.2) and the boundary conditions (1.9). Similarly, the second and third terms on the right-hand side vanish due to Stokes’ theorem, (1.9) and the periodicity in the horizontal direction, showing that the spatial average of 
$u_1$ is conserved. Therefore, due to the Galilean symmetry of the system, we can assume 
$u_1$ to be average free. Consequently, the Poincaré inequality (Evans Reference Evans1998, § 5.8.1) implies that there exists a constant 
$C=C(\varGamma )>0$ such that
\begin{equation} \|u_1\|_{L^2}\leq C\,\|\boldsymbol{\nabla} u_1\|_{L^2}, \end{equation}and combined with (2.18), we find
\begin{equation} \|\boldsymbol{u}\|_{L^2} \leq C\,\|\boldsymbol{\nabla} \boldsymbol{u}\|_{L^2}, \end{equation}
the analogue of (2.21). Using (2.28) instead of (2.21), the arguments corresponding to (2.22)–(2.25) yield the bounds for 
$L_s=\infty$.
In order to bound the second derivatives of 
$\boldsymbol {u}$, we will exploit the equation for the vorticity 
$\omega =\partial _1 u_2-\partial _2 u_1$:
$$\begin{gather} {\textit{Pr}}^{{-}1}\,(\partial_t \omega + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \omega) - \Delta \omega = {\textit{Ra}}\,\partial_1T \quad \text{in } \varOmega, \end{gather}$$
$$\begin{gather}\omega = \frac{1}{L_s}\,u_1 \quad \text{at } x_2=1, \end{gather}$$
$$\begin{gather}\omega ={-}\frac{1}{L_s}\,u_1 \quad \text{at } x_2=0. \end{gather}$$
Note that in the two-dimensional setting, the vorticity is a scalar function, and for any 
$0< L_s\leq \infty$, we have
\begin{equation} \|\boldsymbol{\nabla} \boldsymbol{u}\|_{L^2}=\|\omega\|_{L^2}, \quad \|\boldsymbol{\nabla} \boldsymbol{u}\|_{L^p}\leq C\,\|\omega\|_{L^p}. \end{equation} While the identity in 
$L^2$ follows from a direct computation, the inequality in 
$L^p$ follows by elliptic regularity; in fact, let 
$\psi$ be the the stream function for 
$\boldsymbol {u}$, i.e. 
$\boldsymbol {u}=\boldsymbol {\nabla }^{\perp }\psi =(-\partial _2 \psi, \partial _1 \psi )$. Since 
$\partial _1 \psi = u_2 = 0$ at 
$x_2=1$ and 
$x_2 = 0$, and 
$\psi$ is defined only up to a constant, we can choose it such that 
$\psi = 0$ on 
$x_2=0$. Therefore, using the fundamental theorem of calculus,
\begin{equation} \psi(x_1,1)= \psi(x_1,0) + \int_0^1 \partial_2 \psi(x_1,z) \ {\rm d}z ={-}\int_0^1 u_1(x_1,z) \, {\rm d}z, \end{equation}
and since 
$\psi$ is constant at 
$x_2=1$, averaging (2.33) in 
$x_1$ yields 
$\psi (x_1,1)=-({1}/{\varGamma })\int _{\varOmega } u_1\, {{\rm d} x}$. Combining these observations with the direct computation 
$\Delta \psi =$ 
$\boldsymbol {\nabla }^\perp \boldsymbol {\cdot } \boldsymbol {\nabla }^\perp \psi = \boldsymbol {\nabla }^\perp \boldsymbol {\cdot } \boldsymbol {u} = \omega$ shows that 
$\psi$ is a solution of
\begin{gather} \Delta \psi=\omega \quad \text{in } \varOmega, \end{gather}
\begin{gather}\psi ={-}\frac{1}{\varGamma}\int_\varOmega u_1 \, {{\rm d} x}\quad \text{at } x_2=1, \end{gather}
\begin{gather}\psi = 0\quad \text{at } x_2=0, \end{gather}
and 
$\tilde \psi = \psi + x_2 \int _\varOmega u_1\, {\rm d}\boldsymbol {x}$ solves
\begin{gather} \Delta \tilde\psi=\omega \quad \text{in } \varOmega, \end{gather}
\begin{gather}\psi = 0\quad \text{at } x_2 \in \{0,1\}. \end{gather}One has
\begin{align} \|\boldsymbol{\nabla} \boldsymbol{u}\|_{L^p}^p &=\|\boldsymbol{\nabla} u_1\|_{L^p}^p+\|\boldsymbol{\nabla} u_2\|_{L^p}^p =\|-\boldsymbol{\nabla} \partial_2 \psi\|_{L^p}^p+\|\boldsymbol{\nabla} \partial_1 \psi\|_{L^p}^p \end{align}
\begin{align} &\leq \|\boldsymbol{\nabla}^2 \psi\|_{L^p}^p+\|\boldsymbol{\nabla}^2 \psi\|_{L^p}^p = \|\boldsymbol{\nabla}^2 \tilde \psi\|_{L^p}^p+\|\boldsymbol{\nabla}^2 \tilde\psi\|_{L^p}^p \leq C\|\omega\|_{L^p}^p, \end{align}where we used the Calderon–Zygmund estimate (Gilbarg & Trudinger Reference Gilbarg and Trudinger1977, § 9.4) in the last inequality.
The following lemma provides a higher-order version of Lemma 2.3. The uniform in time bound will later be used to control the pressure terms arising due to the vorticity production on the boundary, while the long-time average estimate will be used in the proof of the main theorem to estimate the thickness of the thermal boundary layer.
Lemma 2.4 (Enstrophy)
Suppose 
$\boldsymbol {u}_0\in W^{1,4}$ and 
$0< L_s \leq \infty$. Then there exists a constant 
$C=C(\varGamma )>0$ such that for all 
$t>0$,
\begin{equation} \|\omega(t)\|_{L^4}\leq C\max \left\{1, L_s^{{-}3}\right\}(\|\boldsymbol{u}_0\|_{W^{1,4}}+{\textit{Ra}}) \end{equation}and
\begin{equation} \left\langle\int_0^1 |\boldsymbol{\nabla} \omega|^2\,{\rm d}\kern 0.06em x_2\right\rangle\leq\frac{1}{L_s}\left|\left\langle p\,\partial_1 u_1|_{x_2=1}\right\rangle+\left\langle p\,\partial_1 u_1|_{x_2=0}\right\rangle\right|+{\textit{Nu}}\,{\textit{Ra}}^{3/2}. \end{equation}Proof. For the proof of (2.41), we refer the reader to Drivas et al. (Reference Drivas, Nguyen and Nobili2022, Lemma 2.11), where the case 
$L_s\geq 1$ is covered. The same argument also yields the bound in the case 
$L_s<1$. Testing the equation with 
$\omega$ and integrating by parts, we obtain
\begin{align}
&\frac{1}{2\,{\textit{Pr}}}\,\frac{{\rm d}}{{\rm
d}t}\left(\|\omega\|_{L^2}^2+\frac{1}{L_s}\,\|u_1\|_{L^2(x_2=0)}^2+\frac{1}{L_s}\,\|u_1\|_{L^2(x_2=1)}^2\right)
+ \|\boldsymbol{\nabla} \omega\|_{L^2}^2\, {{\rm d}
x}\nonumber\\
&\quad +\frac{1}{L_s}\left[\int_0^{\varGamma}\partial_1u_1 p|_{x_2=0}\, {\rm d}\kern 0.06em
x_1 + \int_0^\Gamma \partial_1 u_1 p \vert_{x_2=1} {\rm d}\kern0.7pt x_1\right] = {\textit{Ra}}\int_{\varOmega}\omega\, \partial_1 T\,
{{\rm d} x},
\end{align}
where we used 
$-\partial _2\omega =\Delta u_1=({1}/{{\textit {Pr}}})(\partial _tu_1+(\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla } )u_1)-\partial _1 p$ and the fact that 
$\omega (\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla } )u_1 = \pm ({1}/{L_s}) u_1^2\,\partial _1 u_1$ at 
$x_2=\{0,1\}$ vanishes after integration (in 
$x_1$) due to periodicity. We integrate (2.43) in time, and notice that thanks to (2.32) and the trace estimate (Evans Reference Evans1998, § 5.5), the first bracket on the left-hand side is bounded by the 
$H^1$-norm of 
$\boldsymbol {u}$:
\begin{equation} \|\omega\|_{L^2}^2 + \frac{1}{L_s}\,\|u_1\|_{L^2(x_2=0)}^2+\frac{1}{L_s}\,\|u_1\|_{L^2(x_2=1)}^2 \leq C(L_s)\,\|\boldsymbol{u}\|_{H^1}^2. \end{equation}Note also that due to Hölder's inequality, (2.41) additionally yields
\begin{equation} \|\omega(t)\|_{L^2}\leq C\max \left\{1, L_s^{{-}3}\right\}(\|\boldsymbol{u}_0\|_{W^{1,4}}+{\textit{Ra}}), \end{equation}
which, combined with (2.9) and (2.32), implies that 
$\|\boldsymbol {u}(t)\|_{H^1}$ is universally bounded in time. Therefore, claim (2.42) follows from taking the space and long-time average of (2.43), using the fact that the long-time average of the first term in (2.43) vanishes due to the argument above, and observing
\begin{equation} \left\langle\int_{\varOmega}\omega \partial_1 T\, {\rm d}\kern 0.06em x_2\right\rangle \leq \left\langle\int_{\varOmega}|\omega|^2 \,{\rm d}\kern 0.06em x_2\right\rangle^{1/2} \left\langle\int_{\varOmega}|\boldsymbol{\nabla} T|^2 \,{\rm d}\kern 0.06em x_2\right\rangle^{1/2} \leq ({\textit{Nu}}\,{\textit{Ra}})^{1/2}\,{\textit{Nu}}^{1/2}, \end{equation} Notice that the pressure term appears at the boundary in (2.42), and for this reason, we need to control its 
$H^1$-norm. The next lemma will provide control for this term. Taking the divergence of (1.1), it is easy to see that the pressure 
$p$ satisfies
\begin{gather} \Delta p ={-}\frac{1}{{\textit{Pr}}}\,\boldsymbol{\nabla} \boldsymbol{u}^{\textrm{T}}\boldsymbol{\colon}\boldsymbol{\nabla} \boldsymbol{u}+{\textit{Ra}}\, \partial_2 T \quad \text{in } \varOmega, \end{gather}
\begin{gather}\partial_2 p = \frac{1}{L_s}\,\partial_1u_1\quad \text{at } x_2=1, \end{gather}
\begin{gather}-\partial_2 p = \frac{1}{L_s}\,\partial_1u_1- {\textit{Ra}} \quad \text{at } x_2=0, \end{gather}where the boundary conditions are derived by tracing the second component of (1.1) on the boundary.
Lemma 2.5 (Pressure bound)
Let 
$r>2$. Then there exists a constant 
$C=C(r,\varGamma )>0$ such that
\begin{equation} \|p\|_{H^1}\leq C\left( \frac{1}{L_s}\,\|\partial_2 \boldsymbol{u}\|_{L^2}+\frac{1}{{\textit{Pr}}}\,\|\omega\|_{L^2}\,\|\omega\|_{L^r}+{\textit{Ra}}\,\|T\|_{L^2}\right). \end{equation}Proof. The proof is a consequence of the following two claims:
\begin{align} \|\boldsymbol{\nabla}
p\|_{L^2}^2&=\frac{1}{L_s}\int_0^{\varGamma}(p\,\partial_1u_1|_{x_2=1}+p\,\partial_1u_1|_{x_2=0})\,
{\rm d}\kern 0.06em
x_1+\frac{1}{{\textit{Pr}}}\int_{\varOmega}p\,\boldsymbol{\nabla}
\boldsymbol{u}^{\textrm{T}}\boldsymbol{\colon}\boldsymbol{\nabla}
\boldsymbol{u}\, {{\rm d} x}\nonumber\\
&\quad +{\textit{Ra}}\int_{\varOmega}\partial_2 p T\, {{\rm d} x},
\end{align}
\begin{gather}\left|\int_0^{\varGamma}(p\,\partial_1u_1|_{x_2=1}+p\,\partial_1u_1|_{x_2=0})\, {\rm d}\kern 0.06em x_1\right|\leq 3\,\|p\|_{H^1}\,\|\partial_2 \boldsymbol{u}\|_{L^2}. \end{gather} In fact, applying the Poincaré inequality (
$\kern 1.5pt p$ can be assumed to have zero mean, since 
$p-\langle p\rangle$ satisfies (2.47)) and combining (2.49) and (2.50), we obtain
\begin{align} \| p\|_{H^1}^2 &\lesssim\frac{1}{L_s}\int_0^{\varGamma}(p\,\partial_1u_1|_{x_2=1}+p\,\partial_1u_1|_{x_2=0})\, {\rm d}\kern 0.06em x_1+\frac{1}{{\textit{Pr}}}\int_{\varOmega}p\,\boldsymbol{\nabla} \boldsymbol{u}^T\boldsymbol{\colon}\boldsymbol{\nabla} \boldsymbol{u}\, {{\rm d} x}\nonumber\\ &\quad\ +{\textit{Ra}}\int_{\varOmega}\partial_2 p T\, {{\rm d} x} \end{align}
\begin{align} &\lesssim L_s^{{-}1}\,\|p\|_{H^1}\,\|\partial_2 \boldsymbol{u}\|_{L^2}+{\textit{Pr}}^{{-}1}\,\|p\|_{L^q}\,\|\boldsymbol{\nabla} \boldsymbol{u}\|_{L^2}\,\|\boldsymbol{\nabla} \boldsymbol{u}\|_{L^r}+{\textit{Ra}}\,\|p\|_{H^1}\,\|T\|_{L^2} \end{align}
\begin{align} &\lesssim L_s^{{-}1}\,\|p\|_{H^1}\,\|\partial_2 \boldsymbol{u}\|_{L^2}+{\textit{Pr}}^{{-}1}\,\|p\|_{H^1}\,\|\boldsymbol{\nabla} \boldsymbol{u}\|_{L^2}\,\|\boldsymbol{\nabla} \boldsymbol{u}\|_{L^r}+{\textit{Ra}}\,\|p\|_{H^1}\,\|T\|_{L^2}, \end{align}
where 
$f\lesssim g$ indicates that there exists 
$C>0$ such that 
$f\leq C g$, and 
$q$ and 
$r$ are related by 
$1/p+1/r=1/2$. In (2.52) we used the trace estimate, and in (2.53) the fact that 
$\|p\|_{L^q}\leq C\,\|p\|_{H^1}$ for any 
$q\in (2,\infty )$ by Sobolev embedding in bounded domains. It is left to prove the claims.
Argument for (2.49). Integrating by parts and using (2.47a) gives
\begin{align} \|\boldsymbol{\nabla} p\|_{L^2}^2 &= \int_0^\varGamma p\,\partial_2 p \vert_{x_2=1} \, {\rm d}\kern 0.06em x_1 - \int_0^\varGamma p\,\partial_2 p \vert_{x_2=0} \, {\rm d}\kern 0.06em x_1 - \int_{\varOmega} p\,\Delta p \, {{\rm d} x} \end{align}
\begin{align} &= \frac{1}{L_s} \int_0^\varGamma p (\partial_1 u_1 \vert_{x_2=1} + \partial_1 u_1\vert_{x_2=0}) \, {\rm d}\kern 0.06em x_1 - {\textit{Ra}} \int_0^\varGamma p \vert_{x_2=0} \, {\rm d}\kern 0.06em x_1\nonumber\\ &\quad +\frac{1}{{\textit{Pr}}} \int_{\varOmega} p\,\boldsymbol{\nabla} \boldsymbol{u}^{\textrm{T}}\boldsymbol{\colon}\boldsymbol{\nabla} \boldsymbol{u} \, {{\rm d} x} - {\textit{Ra}} \int_{\varOmega}\partial_2 T p \, {{\rm d} x} \end{align}
\begin{align} &= \frac{1}{L_s} \int_0^\varGamma p (\partial_1 u_1 \vert_{x_2=1} + \partial_1 u_1\vert_{x_2=0}) \, {\rm d}\kern 0.06em x_1\nonumber\\ &\quad +\frac{1}{{\textit{Pr}}} \int_{\varOmega} p\,\boldsymbol{\nabla} \boldsymbol{u}^{\textrm{T}}\boldsymbol{\colon}\boldsymbol{\nabla} \boldsymbol{u} \, {{\rm d} x} + {\textit{Ra}} \int_{\varOmega} T\,\partial_2p \, {{\rm d} x}, \end{align}
where in the last identity we used 
$- {\textit {Ra}} \int _{\varOmega }\partial _2 T p- {\textit {Ra}} \int _0^\varGamma p \vert _{x_2=0} \, {\rm d}\kern 0.06em x_1= {\textit {Ra}} \int _{\varOmega } T\,\partial _2 p \, {{\rm d} x}$ thanks to the boundary conditions for 
$T$.
Argument for (2.50). Since 
$\boldsymbol {u}$ is divergence-free, we have
\begin{equation} 0=\int_{\varOmega}p({-}1+2x_2)\,\partial_2(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u})\, {{\rm d} x}=\int_{\varOmega}p({-}1+2x_2)\,\boldsymbol{\nabla} \boldsymbol{\cdot} (\partial_2\boldsymbol{u})\, {{\rm d}\kern0.06em x}, \end{equation}and integration by parts yields
\begin{align} \int_{\varOmega}p({-}1+2x_2)\,\boldsymbol{\nabla} \boldsymbol{\cdot} (\partial_2\boldsymbol{u})\, {{\rm d} x} &={-}\int_{\varOmega}\boldsymbol{\nabla} p\boldsymbol{\cdot} ({-}1+2x_2)\,\partial_2 \boldsymbol{u}\, {{\rm d} x} - 2\int_{\varOmega}p\,\partial_2 u_2\, {{\rm d} x}\nonumber\\ &\quad -\left(\int_0^{\varGamma}p\,\partial_1 u_1|_{x_2=1}+p\,\partial_1 u_1|_{x_2=0}\right)\, {\rm d}\kern 0.06em x_1, \end{align}
where we used 
$\partial _1 u_1=-\partial _2 u_2$ by incompressibility. Combining the two identities, we obtain
\begin{align} \left|\int_0^{\varGamma}(p\,\partial_1 u_1|_{x_2=1}+p\,\partial_1 u_1|_{x_2=0})\, {\rm d}\kern 0.06em x_1\right| &\leq 2\,\|p\|_{L^{2}}\,\|\partial_2 u_2\|_{L^2}+\|\partial_2\boldsymbol{u}\|_{L^2}\,\|\boldsymbol{\nabla} p\|_{L^2} \end{align}
\begin{align} &\leq 3\,\|p\|_{H^1}\,\|\partial_2\boldsymbol{u}\|_{L^2}. \end{align}Remark 2.6 The pressure bound in Proposition 2.7 of Drivas et al. (Reference Drivas, Nguyen and Nobili2022) is given by
\begin{equation} \|p\|_{H^1}\leq C\left( \frac{1}{L_s}\,\|\partial_1 \omega\|_{L^2}+\frac{1}{{\textit{Pr}}}\,\|\omega\|_{L^2}\,\|\omega\|_{L^r}+{\textit{Ra}}\,\|T\|_{L^2}\right). \end{equation}
In the subsequent analysis, the authors bound the term 
$({1}/{L_s})\,\|\partial _1 \omega \|_{L^2}$ from above with 
$\|\boldsymbol {\nabla } \omega \|_{L^2}$, imposing the condition that 
$L_s\geq 1$. In contrast, using the refined estimate (2.48), we are able to treat any slip length 
$L_s>0$, improving the scaling of 
${\textit {Nu}}$ with respect to 
${\textit {Ra}}$.
We conclude this section with bounds on second derivatives of the velocity field. First, we relate the 
$L^2$-norm of 
$\boldsymbol {\nabla } ^2 \boldsymbol {u}$ to the 
$L^2$-norm of 
$\boldsymbol {\nabla } \omega$.
Lemma 2.7
\begin{equation} \|\boldsymbol{\nabla} ^2 \boldsymbol{u}\|_{L^2}\leq\|\boldsymbol{\nabla} \omega\|_{L^2}. \end{equation}Proof. First, we show that 
$\|\boldsymbol {\nabla } ^2 \boldsymbol {u}\|_{L^2}\leq \|\Delta \boldsymbol {u}\|_{L^2}$: integrating by parts twice yields
\begin{align} \|\boldsymbol{\nabla}^2 \boldsymbol{u}\|_{L^2(\varOmega)}^2&=\int_{\varOmega}\partial_i\partial_j u_k\,\partial_i\partial_j u_k\, {{\rm d} x} \end{align}
\begin{align} &=\int_{\varOmega}\partial_i^2 u_k\,\partial_j^2 u_k \, {{\rm d} x}-\int_{\partial\varOmega}\partial_i^2 u_k\,\partial_j u_k n_j \, {\rm d}S+\int_{\partial\varOmega}\partial_i\partial_j u_k\,\partial_j u_k n_i \, {\rm d}S \end{align}
\begin{align} &=\|\Delta \boldsymbol{u}\|_{L^2(\varOmega)}^2-\int_0^\varGamma\partial_1^2 u_k\,\partial_2 u_k n_2 \, {\rm d}\kern 0.06em x_1+\int_0^\varGamma\partial_2\partial_1 u_k\,\partial_1 u_k n_2 \, {\rm d}\kern 0.06em x_1, \end{align}
where we used periodicity in the horizontal direction and the fact that the terms with 
$i=j$ cancel. Note that due to (1.9), the boundary terms have a sign:
\begin{align} &-\int_0^\varGamma\partial_1^2 u_k\,\partial_2 u_k n_2 \, {\rm d}\kern 0.06em x_1+\int_0^\varGamma\partial_2\partial_1 u_k\,\partial_1 u_k n_2 \, {\rm d}\kern 0.06em x_1\nonumber\\ &\qquad ={-}\int_0^\varGamma\partial_1^2 u_1\,\partial_2 u_1 n_2 \, {\rm d}\kern 0.06em x_1+\int_0^\varGamma\partial_2\partial_1 u_1\,\partial_1 u_1 n_2 \, {\rm d}\kern 0.06em x_1 \end{align}
\begin{align} &\qquad = \frac{1}{L_s}\int_0^\varGamma\partial_1^2 u_1 u_1 \, {\rm d}\kern 0.06em x_1-\int_0^\varGamma(\partial_1 u_1)^2\, {\rm d}\kern 0.06em x_1 ={-}\frac{2}{L_s}\int_0^\varGamma(\partial_1 u_1)^2\, {\rm d}\kern 0.06em x_1 \leq 0, \end{align}where in the last identity we used the periodicity in the horizontal direction. This proves the first claim.
Now, a direct computation yields 
$\Delta \boldsymbol {u}=\boldsymbol {\nabla }^{\perp }\omega$, and from this follows 
$\|\Delta \boldsymbol {u}\|_{L^2(\varOmega )}=\|\boldsymbol {\nabla }\omega \|_{L^2(\varOmega )}$. We conclude that
\begin{equation} \|\boldsymbol{\nabla}^2 \boldsymbol{u}\|_{L^2(\varOmega)} \leq \|\Delta \boldsymbol{u}\|_{L^2(\varOmega)}=\|\boldsymbol{\nabla}\omega\|_{L^2(\varOmega)}, \end{equation}yielding (2.62).
Next, we relate 
$\|\boldsymbol {\nabla}^2 \boldsymbol {u}\|_{L^2}$ to 
${\textit {Nu}},{\textit {Ra}}$ and 
$L_s$, via an upper bound for 
$\|\boldsymbol {\nabla } \omega \|_{L^2}$. This is done by combining (2.42), the pressure bound (2.48) and the boundary integral estimate (2.50). The resulting bound on the long-time average of the velocity Hessian together with the corresponding bound on the velocity gradient of Lemma 2.10 are key ingredients in order to estimate the boundary layer thickness in the proof of the main theorem.
Lemma 2.8 (Hessian bound)
Let 
$L_s>0$ and 
$\boldsymbol {u}_0\in W^{1,4}$. Then there exists a constant 
$C=C(\varGamma )>0$ such that
\begin{align} \left\langle\|\boldsymbol{\nabla}^2
\boldsymbol{u}\|_{L^2}^2 \right\rangle &\leq C \left(
L_s^{{-}2}+\frac{\max\{1,L_s^{{-}3}\}\,(\|\boldsymbol{u}_0\|_{W^{1,4}}+{\textit{Ra}})}{{\textit{Pr}}\,
L_s}\right.\nonumber\\ &\quad \left. +\,L_s^{{-}1}\,{\textit{Nu}}^{{-}1/2}\,{\textit{Ra}}^{1/2}+{\textit{Ra}}^{1/2}\vphantom{\left(
L_s^{{-}2}+\frac{\max\{1,L_s^{{-}3}\}\,(\|\boldsymbol{u}_0\|_{W^{1,4}}+{\textit{Ra}})}{{\textit{Pr}}\,
L_s}\right.}\right){\textit{Nu}}\,{\textit{Ra}}.
\end{align}Proof. From (2.62) and (2.42), we have
\begin{equation} \left\langle \int_0^1 |\boldsymbol{\nabla}^2 \boldsymbol{u}|^2\, {{\rm d} x}\right\rangle\leq L_s^{{-}1}\left(\left\langle p\,\partial_1 u_1|_{x_2=1}\right\rangle+\left\langle p\,\partial_1 u_1|_{x_2=0}\right\rangle\right)+{\textit{Nu}}\,{\textit{Ra}}^{3/2}. \end{equation}Using (2.50) and (2.48), we obtain
\begin{align} \left\langle \int_0^1 |\boldsymbol{\nabla}^2 \boldsymbol{u}|^2\, {{\rm d} x}\right\rangle &\leq L_s^{{-}1}\langle\|p\|_{H^1}\|\partial_2 \boldsymbol{u}\|_{L^2}\rangle+ {\textit{Nu}}\,{\textit{Ra}}^{3/2} \end{align}
\begin{align} &\leq L_s^{{-}2}\langle\|\partial_2 \boldsymbol{u}\|_{L^2}^2\rangle+L_s^{{-}1}\,{\textit{Pr}}^{{-}1}\langle\|\omega\|_{L^2}\,\|\omega\|_{L^4}\,\|\partial_2 \boldsymbol{u}\|_{L^2}\rangle\nonumber\\ &\quad +L_s^{{-}1}\,{\textit{Ra}}\,\langle\|T\|_{L^2}\,\|\partial_2 \boldsymbol{u}\|_{L^2}\rangle+{\textit{Nu}}\,{\textit{Ra}}^{3/2}. \end{align}
For the 
$L^4$-norm of the vorticity, the pointwise (in time) bound (2.41) states
\begin{equation} \|\omega(t)\|_{L^4}\leq C\max \left\{1, L_s^{{-}3}\right\}(\|\boldsymbol{u}_0\|_{W^{1,4}}+{\textit{Ra}}), \end{equation}
for some constant 
$C>0$ depending only on 
$\varGamma$. Therefore, additionally using the identity (2.32), (2.72) can be estimated as follows:
\begin{align} \left\langle \int_0^1 |\boldsymbol{\nabla}^2 \boldsymbol{u}|^2\, {\rm d}\kern 0.06em x_2\right\rangle &\leq L_s^{{-}2}\langle\|\partial_2 \boldsymbol{u}\|_{L^2}^2\rangle +\frac{C\max\{1,L_s^{{-}3}\}}{{\textit{Pr}}\,L_s}\,\langle\|\boldsymbol{\nabla} \boldsymbol{u}\|_{L^2}(\|\boldsymbol{u}_0\|_{W^{1,4}}+{\textit{Ra}})\,\|\partial_2 \boldsymbol{u}\|_{L^2}\rangle \nonumber\\ &\quad +L_s^{{-}1}\,{\textit{Ra}}\,\langle\|T\|_{L^2}\,\|\partial_2 \boldsymbol{u}\|_{L^2}\rangle+{\textit{Nu}}\,{\textit{Ra}}^{3/2}. \end{align}
Finally, using 
$\langle \|\partial _2 \boldsymbol {u}\|_{L^2}^2\rangle \leq \langle \|\boldsymbol {\nabla } \boldsymbol {u}\|_{L^2}^2\rangle$, the upper bound in (2.10) and the maximum principle for the temperature (2.1), we deduce
\begin{align} \left\langle \int_0^1 |\boldsymbol{\nabla}^2 \boldsymbol{u}|^2\, {\rm d}\kern 0.06em x_2\right\rangle &\leq L_s^{{-}2}\,{\textit{Nu}}\,{\textit{Ra}} +\frac{C\max\{1,L_s^{{-}3}\}}{{\textit{Pr}}\,L_s}\,(\|\boldsymbol{u}_0\|_{W^{1,4}}+{\textit{Ra}})\,{\textit{Nu}}\,{\textit{Ra}} \nonumber\\ &\quad +L_s^{{-}1}\varGamma^{1/2}\,{\textit{Nu}}^{1/2}\,{\textit{Ra}}^{3/2}+{\textit{Nu}}\,{\textit{Ra}}^{3/2}. \end{align}3. Proof of Theorem 1.1
A crucial ingredient of the proof of Theorem 1.1 is the following interpolation bound, relating the integral of the product 
$u_2T$ with 
$\|\boldsymbol {\nabla } T\|_{L^2}$, 
$\|\boldsymbol {\nabla } \boldsymbol {u}\|_{L^2}$ and 
$\|\boldsymbol{\nabla}^2 \boldsymbol {u}\|_{L^2}$. Notice that, in turn, these quantities are estimated in terms of Nusselt, Rayleigh and Prandtl numbers in (2.7), (2.10) and (2.69).
Lemma 3.1 The interpolation bound
\begin{equation} \frac{1}{\delta}\left\langle \int_{0}^{\delta} u_2T\,{\rm d}\kern 0.06em x_2\right\rangle\leq \frac{1}{2}\left\langle \int_0^1|\partial_2 T|^2\, {\rm d}\kern 0.06em x_2\right\rangle+C\delta^3\left\langle \int_0^1|\boldsymbol{\nabla} \boldsymbol{u}|^2\, {\rm d}\kern 0.06em x_2\right\rangle^{1/2}\left\langle \int_0^1|\boldsymbol{\nabla}^2 \boldsymbol{u}|^2\, {\rm d}\kern 0.06em x_2\right\rangle^{1/2} \end{equation}
holds for some constant 
$C>0$.
Proof. First, notice that due to incompressibility 
$\partial _2u_2=-\partial _1u_1$ and horizontal periodicity,
\begin{equation} \partial_2\,\frac{1}{\varGamma}\int_0^\varGamma u_2(x_1,x_2)\, {\rm d}\kern 0.06em x_1={-}\frac{1}{\varGamma}\int_0^\varGamma \partial_1 u_1(x_1,x_2)\, {\rm d}\kern 0.06em x_1 = 0, \end{equation}implying
\begin{equation} \frac{1}{\varGamma}\int_0^\varGamma u_2(x_1,x_2)\, {\rm d}\kern 0.06em x_1=0 \end{equation}
thanks to the boundary conditions 
$u_2=0$ at 
$x_2=\{0,1\}$. Let 
$\theta (x_1,x_2):=T(x_1,x_2)-1$; then 
$\theta (x_1, x_2)=0$ at 
$x_2=0$, and by (3.3),
\begin{equation} \frac{1}{\delta}\left\langle \int_{0}^{\delta} u_2T\,{\rm d}\kern 0.06em x_2\right\rangle=\frac{1}{\delta}\left\langle \int_{0}^{\delta} u_2\theta\,{\rm d}\kern 0.06em x_2\right\rangle. \end{equation}
By the fundamental theorem of calculus and the homogeneous Dirichlet boundary conditions for 
$\theta$ and 
$u_2$, we have
\begin{gather} |\theta(x_1,x_2)|\leq x_2^{1/2}\,\|\partial_2\theta\|_{L^2(0,1)}, \end{gather}
\begin{gather}|u_2(x_1,x_2)|\leq x_2\,\|\partial_2u_2\|_{L^{\infty}(0,1)}. \end{gather}
Furthermore, since 
$\int _0^1\partial _2u_2(x_1,x_2)\, {\rm d}\kern 0.06em x_2=0$, there exists 
$\xi =\xi (x_1)\in (0,1)$ such that 
$\partial _2u_2(x_1,\xi )=0$. Then
\begin{equation} |\partial_2u_2(x_1,x_2)|^2=\left|\int_{\xi}^{x_2}\partial_2(\partial_2 u_2(x_1,z))^2\, {\rm d}z\right|=\left|2\int_{\xi}^{x_2}\partial_2 u_2\,\partial_2^2 u_2 \,{\rm d}z\right|, \end{equation}which implies
\begin{equation} |u_2(x_1,x_2)|\leq x_2\,\|\partial_2u_2(x_1,\cdot)\|_{L^2(0,1)}^{1/2}\,\|\partial_2^2u_2(x_1,\cdot)\|_{L^2(0,1)}^{1/2}. \end{equation}Combining these estimates, we obtain
\begin{align} & \frac{1}{\delta}\left\langle \int_0^{\delta}u_2\theta\, {\rm d}\kern 0.06em x_2\right\rangle \nonumber\\ &\qquad\leq \frac{1}{\delta}\left\langle \int_0^{\delta}x_2^{3/2}\, \|\partial_2u_2(x_1,\cdot)\|_{L^2(0,1)}\,\|\partial_2^2u_2(x_1,\cdot)\|_{L^2(0,1)}\,\|\partial_2\theta(x_1,\cdot)\|_{L^2(0,1)}^2\, {\rm d}\kern 0.06em x_2\right\rangle \end{align}
\begin{align} &\qquad\leq \delta^{3/2}\left\langle \|\partial_2u_2(x_1,\cdot)\|_{L^2(0,1)}\,\|\partial_2^2u_2(x_1,\cdot)\|_{L^2(0,1)}\right\rangle^{1/2}\left\langle \|\partial_2\theta(x_1,\cdot)\|_{L^2(0,1)}^2 \right\rangle^{1/2} \end{align}
\begin{align} &\qquad\leq C\delta^{3}\left\langle \|\partial_2u_2(x_1,\cdot)\|_{L^2(0,1)}^2\right\rangle^{1/2}\left\langle \|\partial_2^2u_2(x_1,\cdot)\|_{L^2(0,1)}^2\right\rangle^{1/2} \!+\!\frac 12\left\langle \|\partial_2\theta(x_1,\cdot)\|_{L^2(0,1)}^2 \right\rangle. \end{align}3.1. Case 
$L_s=\infty$
In this subsection, we prove the upper bound (1.17) when 
$L_s=\infty$ in (1.9), re-deriving the seminal result of Whitehead & Doering (Reference Whitehead and Doering2011) with a different technique. While the proof in Whitehead & Doering (Reference Whitehead and Doering2011) is a sophisticated application of the background field method, our new proof for the 
${\textit {Ra}}^{{5}/{12}}$ scaling is a pure partial differential equation argument based on the combination of the localization principle together with the interpolation bound (1.22).
We first notice that setting 
$L_s=\infty$ in (2.10), (2.32) and (2.42), we obtain
\begin{gather} \left\langle \int_0^{1}|\omega|^2\, {\rm d}\kern 0.06em x_2\right\rangle=\left\langle \int_0^{1}|\boldsymbol{\nabla} \boldsymbol{u}|^2\, {\rm d}\kern 0.06em x_2\right\rangle\leq {\textit{Nu}}\,{\textit{Ra}}, \end{gather}
\begin{gather}\left\langle \int_0^{1}|\partial_1\omega|^2\, {\rm d}\kern 0.06em x_2\right\rangle\leq\left\langle \int_0^{1}|\boldsymbol{\nabla} \omega|^2\, {\rm d}\kern 0.06em x_2\right\rangle \leq {\textit{Nu}}\,{\textit{Ra}}^{3/2}. \end{gather}Combining these upper bounds and the identity (2.7) in the localization estimate (1.21), we find
\begin{equation} {\textit{Nu}}\leq \frac{1}{2}\,{\textit{Nu}} +C\delta^3({\textit{Nu}}\,{\textit{Ra}})^{1/2}({\textit{Ra}}^{3/2}\,{\textit{Nu}})^{1/2}+\frac{1}{\delta}, \end{equation}which yields
\begin{equation} \frac 12\,{\textit{Nu}}\leq C\delta^3\,{\textit{Nu}}\,{\textit{Ra}}^{5/4}+\frac{1}{\delta}. \end{equation}
Optimizing in 
$\delta$,
\begin{equation} \delta\sim\frac{1}{{\textit{Ra}}^{{5}/{16}}\,{\textit{Nu}}^{{1}/{4}}}, \end{equation}we deduce
\begin{equation} {\textit{Nu}}\lesssim{\textit{Ra}}^{{5}/{12}}. \end{equation}3.2. Case $0< L_s<\infty$
Using the localization of the Nusselt number (2.6) and (3.1), we have
\begin{align} {\textit{Nu}}&\leq \frac{1}{\delta}\left\langle \int_{0}^{\delta} u_2T\,{\rm d}\kern 0.06em x_2\right\rangle+\frac{1}{\delta} \end{align}
\begin{align} &\leq \frac{1}{2}\left\langle \int_0^{1}|\boldsymbol{\nabla} T|^2\, {\rm d}\kern 0.06em x_2\right\rangle+C\delta^3\left\langle \int_0^{1}|\boldsymbol{\nabla} \boldsymbol{u}|^2\, {\rm d}\kern 0.06em x_2\right\rangle^{1/2}\left\langle \int_0^{1}|\boldsymbol{\nabla}^2 \boldsymbol{u}|^2\, {\rm d}\kern 0.06em x_2\right\rangle^{1/2}+\frac{1}{\delta}, \end{align}yielding
\begin{equation} \frac 12\,{\textit{Nu}}\leq C\delta^3\left\langle \int_0^{\delta}|\boldsymbol{\nabla} \boldsymbol{u}|^2\, {\rm d}\kern 0.06em x_2\right\rangle^{1/2}\left\langle \int_0^{1}|\boldsymbol{\nabla}^2 \boldsymbol{u}|^2\, {\rm d}\kern 0.06em x_2\right\rangle^{1/2}+\frac{1}{\delta}, \end{equation}
where we used (2.7). Finally, we insert the bounds (2.10) and (2.69) in the last inequality, and if 
${\textit {Ra}}$ is large enough such that 
$\|\boldsymbol {u}_0\|_{W^{1,4}}\lesssim {\textit {Ra}}$, then up to redefining constants, we find
\begin{align} {\textit{Nu}}&\leq
C\delta^3\left(L_s^{{-}1}\,{\textit{Nu}}\,{\textit{Ra}}+\max\{1,L_s^{-{3}/{2}}\}\,L_s^{-{1}/{2}}\,{\textit{Pr}}^{-{1}/{2}}\,{\textit{Nu}}\,{\textit{Ra}}^{3/2}\right.\nonumber\\
&\left.\quad +\,L_s^{-{1}/{2}}\,{\textit{Nu}}^{3/4}\,{\textit{Ra}}^{5/4}+{\textit{Nu}}\,{\textit{Ra}}^{5/4}\right)+\frac{2}{\delta}.
\end{align}
We first cover the case 
$L_s\geq 1$. Then, as 
$L_s,{\textit {Nu}}\geq 1$, the first and third terms on the right-hand side of (3.21) are dominated by 
${\textit {Nu}}\,{\textit {Ra}}^{5/4}$, hence
\begin{equation} {\textit{Nu}}\leq C\delta^{3}\,{\textit{Nu}}\left(L_s^{-{1}/{2}}\,{\textit{Pr}}^{-{1}/{2}}\,{\textit{Ra}}^{{3}/{2}}+{\textit{Ra}}^{{5}/{4}}\right)+ 2\delta^{{-}1}, \end{equation}
and optimizing by setting 
$\delta = {\textit {Nu}}^{-{1}/{4}}(L_s^{-{1}/{2}}\,{\textit {Pr}}^{-{1}/{2}}\,{\textit {Ra}}^{{3}/{2}}+{\textit {Ra}}^{{5}/{4}})^{-{1}/{4}}$ implies
\begin{equation} {\textit{Nu}} \lesssim L_s^{-{1}/{6}}\,{\textit{Pr}}^{-{1}/{6}}\,{\textit{Ra}}^{{1}/{2}}+{\textit{Ra}}^{{5}/{12}}. \end{equation}
If instead 
$L_s<1$, then (3.21) is given by
\begin{equation} {\textit{Nu}}\leq C\delta^3\left(L_s^{{-}1}\,{\textit{Nu}}\,{\textit{Ra}}+L_s^{{-}2}\,{\textit{Pr}}^{-{1}/{2}}\,{\textit{Nu}}\,{\textit{Ra}}^{3/2}+ L_s^{-{1}/{2}}\,{\textit{Nu}}^{3/4}\,{\textit{Ra}}^{5/4}+{\textit{Nu}}\,{\textit{Ra}}^{5/4}\right)+2 \delta^{{-}1}. \end{equation}Optimizing by setting
\begin{equation} \delta= \begin{cases}
{\textit{Nu}}^{-{1}/{4}}\left(L_s^{{-}1}\,{\textit{Ra}}\right.\\ \quad\left. +\, L_s^{{-}2}\,{\textit{Pr}}^{{-}1/2}\,{\textit{Ra}}^{3/2}+{\textit{Ra}}^{5/4}\right)^{-{1}/{4}},
& \text{if }
L_s^{-{1}/{2}}\,{\textit{Ra}}^{-{1}/{4}}\\ &\quad +\,L_s^{-{3}/{2}}\,{\textit{Pr}}^{{-}1/2}\,{\textit{Ra}}^{1/4}+L_s^{{1}/{2}}
\geq {\textit{Nu}}^{-{1}/{4}},\\
L_s^{{1}/{8}}\,{\textit{Nu}}^{-{3}/{16}}\,{\textit{Ra}}^{-{5}/{16}},
& \text{if }
L_s^{-{1}/{2}}\,{\textit{Ra}}^{-{1}/{4}}+L_s^{-{3}/{2}}\,{\textit{Pr}}^{{-}1/2}\,{\textit{Ra}}^{1/4}\\ &\quad +\,L_s^{{1}/{2}}
\leq {\textit{Nu}}^{-{1}/{4}}, \end{cases}
\end{equation}yields
\begin{equation} {\textit{Nu}} \lesssim L_s^{-{1}/{3}}\,{\textit{Ra}}^{{1}/{3}} + L_s^{-{2}/{3}}\,{\textit{Pr}}^{-{1}/{6}}\,{\textit{Ra}}^{{1}/{2}}+ L_s^{-{2}/{13}}\,{\textit{Ra}}^{{5}/{13}}+{\textit{Ra}}^{{5}/{12}}. \end{equation}4. Conclusion
Our work is motivated by recent investigations in Whitehead & Doering (Reference Whitehead and Doering2011) and Drivas et al. (Reference Drivas, Nguyen and Nobili2022) indicating that variations in the boundary conditions for the velocity affect heat transport properties in Rayleigh–Bénard convection. The Navier-slip boundary conditions are prescribed in the situation in which some slip occurs on the surface, and at the same time stress is exerted on the fluid. In particular, the use of these boundary conditions is justified when the scale of interest goes down to microns or below, when the no-slip boundary conditions cease to be valid (see Choi & Kim Reference Choi and Kim2006). Mathematically, these boundary conditions are referred to as ‘interpolation’ boundary conditions since they represent a situation in between the physically relevant, but very difficult to treat, no-slip boundary conditions and the easier, but unphysical, free-slip boundary conditions. When free-slip boundary conditions are considered, the key tool in order to control the growth of the vertical velocity 
$u_2$ near the boundaries is the enstrophy balance (Whitehead & Doering Reference Whitehead and Doering2011, Reference Whitehead and Doering2012; Wang & Whitehead Reference Wang and Whitehead2013). When Navier-slip boundary conditions are considered instead, the control of enstrophy production becomes difficult since the pressure term and the vertical derivative of 
$u_2$, 
$\partial _2u_2=-\partial _1 u_1$, appear at the boundary; see (2.42). This problem was already tackled in Drivas et al. (Reference Drivas, Nguyen and Nobili2022), and in this paper, we improve this control by refining the trace estimates; see (2.50). In Bleitner & Nobili (Reference Bleitner and Nobili2024), the analysis in Drivas et al. (Reference Drivas, Nguyen and Nobili2022) was extended to rough boundary conditions, capturing the dependency of the scaling from the spatially varying friction coefficient and curvature. To our knowledge, no other rigorous results are available for the Nusselt number in this set-up.
In the same setting as considered in the present paper, direct numerical simulations (DNS) were performed in Huang & He (Reference Huang and He2022) to study the scaling of the Nusselt number in the convective roll state and in the zonal flow (turbulent state). The authors found that, when 
${L_s}/{\lambda _0}\lesssim 10$ (here the slip length is normalized by 
$\lambda _0$, the thermal boundary layer thickness for the no-slip plates), 
${\textit {Nu}}\sim {\textit {Ra}}^{0.31}$ in the convection roll state. In particular, they observe that in the convection roll state, for a fixed 
${\textit {Ra}}$, the heat transfer 
${\textit {Nu}}$ increases with increasing 
${L_s}/{\lambda _0}$. On the other hand, they observe that in the zonal flow, the heat transport 
${\textit {Nu}}$ decreases with increasing 
${L_s}/{\lambda _0}$ at fixed 
${\textit {Ra}}$. They explain this phenomenon as being due to the reduction of the vertical Reynolds number: as 
${L_s}/{\lambda _0}$ increases, the zonal flow becomes stronger, and it suppresses the vertical velocity, leading to a decrease in 
${\textit {Nu}}$. In this regime, their DNS indicate the scaling 
${\textit {Nu}}\sim {\textit {Ra}}^{0.16}$. As can be noticed in table 1, our rigorous upper bounds indicate a behaviour in agreement with the results in Huang & He (Reference Huang and He2022): in the turbulent regime, as the slip length 
$L_s$ increases, the heat transport decreases. In fact, observe that in each 
${\textit {Pr}}$ region, the heat transport increases from 
${\textit {Nu}}\lesssim {\textit {Ra}}^{{5}/{12}}$ when 
$L_s=\infty$, to 
${\textit {Nu}}\lesssim L_s^{-2/3}\,{\textit {Pr}}^{-1/6}\,{\textit {Ra}}^{{1}/{2}}$ when 
$L_s<1$ (see table 1, where the upper bounds are ordered from the smallest to the largest in each 
$Pr$ region). Let us recall that for no-slip boundary conditions (
$L_s=0$) Otto, Choffrut and the second author of this paper find 
${\textit {Nu}}\lesssim ({\textit {Ra}}\ln ({\textit {Ra}}))^{1/3}$ when 
${\textit {Pr}}\gtrsim ({\textit {Ra}}\ln ({\textit {Ra}}))^{1/3}$ and 
${\textit {Nu}}\lesssim {\textit {Pr}}^{-1/2}({\textit {Ra}}\ln ({\textit {Ra}}))^{1/2}$ when 
${\textit {Pr}}\lesssim ({\textit {Ra}}\ln ({\textit {Ra}}))^{1/3}$ (Choffrut et al. Reference Choffrut, Nobili and Otto2016), while Doering & Constantin (Reference Doering and Constantin1996) prove 
${\textit {Nu}}\lesssim {\textit {Ra}}^{{1}/{2}}$ uniformly in 
${\textit {Pr}}$. In the regime of small 
$L_s$, our upper bound is 
${\textit {Nu}}\lesssim L_s^{-2/3}\,{\textit {Pr}}^{-1/6}\,{\textit {Ra}}^{1/2}$ for any Prandtl number. We cannot directly compare the results in the present paper with the results available for no-slip boundary conditions in Choffrut et al. (Reference Choffrut, Nobili and Otto2016) and Doering & Constantin (Reference Doering and Constantin1996) since our analysis breaks down in the limit 
$L_s\rightarrow 0$. With respect to the long-standing open problem regarding the scaling laws for the Nusselt number in the ‘ultimate state’ (Zhu et al. Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018, Reference Zhu, Mathai, Stevens, Verzicco and Lohse2019; Doering et al. Reference Doering, Toppaladoddi and Wettlaufer2019; Doering Reference Doering2020), we can say the following. On the one hand, the DNS in Huang & He (Reference Huang and He2022) and our analysis seem to suggest that, for any Prandtl number region and any 
$L_s\neq 0$, the Malkus scaling 
${\textit {Nu}}\sim {\textit {Ra}}^{{1}/{3}}$ is not achievable in the turbulent regime. This is quite evident when looking at the results in table 1, since (as mentioned before) the upper bounds are ordered and increase from 
${\textit {Ra}}^{{5}/{12}}$ as 
$L_s$ decreases. However, we warn that the reality of heat transport might be much more complicated than this nice picture. With our arguments we are able to derive only upper bounds, and cannot exclude the possibility of smaller scaling exponents. Indeed, one of the biggest theoretical challenges for this problem is to find (non-trivial) lower bounds for the Nusselt number. Even just exhibiting solutions that attain certain upper bounds in some simplified situation (in the spirit of Souza & Doering Reference Souza and Doering2015) would be a big step ahead in the theory. Interesting ideas in this direction have been developed by Tobasco & Doering (Reference Tobasco and Doering2017) and Tobasco (Reference Tobasco2022): inspired by problems arising in the study of energy-driven pattern formation in materials science, the authors design a two-dimensional ‘branching’ flow that transports at rate 
${\textit {Nu}}\sim {{Pe}}^{2/3}$ (up to a logarithmic correction) for a passive tracer that diffuses and is advected by a divergence-free velocity field with a fixed enstrophy budget 
$\langle \|\nabla \boldsymbol {u}\|_{L^2}^2\rangle \sim {{Pe}}^2$. Here, 
${{Pe}}$ is the Péclet number. The remarkable lower bound that they derive using ‘branching’ techniques is responsible for the logarithmic corrections. Although their result in particular shows that the 
${\textit {Ra}}^{1/2}$ ‘ultimate’ scaling is attained by flows that do not solve (1.1)–(1.3), these arguments can be used to detect mechanisms that make the scalings 
${Ra}^{1/2}$ and 
${Ra}^{1/3}$ realizable for buoyancy-driven solutions.
Finally, we want to conclude with a remark. Our upper bounds on the Nusselt number quantify heat transport properties of the (thermal) boundary layer without characterizing it. The only information that we can extrapolate from our analysis is about the (thermal) boundary layer thickness (see figure 1). On the other hand, in the excellent work Gie & Whitehead (Reference Gie and Whitehead2019), the authors are able to describe the flow in the boundary layer. In a three-dimensional periodic channel, they study the Boussinesq system at very small viscosities and with Navier-slip boundary conditions. Through the explicit construction of a corrector, they show that the behaviour of the dynamics in the boundary layers is governed by the Prandtl equations. In particular, these equations can be linearized for Navier-slip boundary conditions, implying that in this specific case, the boundary layers are non-turbulent. From this, the authors deduce that the ‘ultimate’ state, which is based on the hypothesis of a turbulent boundary layer (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009), cannot exist in this particular set-up.
Figure 1. Thermal boundary layers arising at 
${\textit {Pr}}=1$, 
$L_s=1$ and 
${\textit {Ra}}=10^8$. This snapshot is generated by a simulation.
Acknowledgements
The authors acknowledge support by the Deutsche Forschungsgemeinschaft (DFG) within the Research Training Group GRK 2583 ‘Modeling, Simulation and Optimization of Fluid Dynamic Applications.’
Declaration of interests
The authors report no conflict of interest.
Open access
