Chaotic behavior and controlling chaos in a fast-slow plankton-fish model

1.   Introduction

The Rayleigh-Bénard convection model describes a buoyancy driven flow of a fluid heated from below and cooled from above. In a finite box with height $h$, the impermeable horizontal bottom an top plates are held at temperature $T_b$ and $T_t$ respectively, with $T_t < T_b$. Temperature differences trigger density differences within the fluid layers, which, in turn, generate convective motions. As the intensity of the buoyancy forces grows, the dynamics undergo a series of bifurcations: convection rolls are destabilized by plume-shaped structures forming at the boundaries and the motion becomes unpredictable, chaotic, turbulent. This model is a paradigm of nonlinear dynamics and pattern formation and has a myriad of applications in engineering of heat transfer ^[1]. The motion within the finite box is governed by an advection-diffusion equation for the temperature field $T$ coupled with the Navier-Stokes equations for the velocity field $\bf{u}$. Unless otherwise stated, we will focus on three-dimensional dynamics and denote by $\bf{x} = (\bf{x}', z) = (x, y, z)$ a vector in $\mathbb{R}^3$ and by $\bf{e_z} = (0, 0, 1)^t$ the vertical normal vector. In dimensionless variables^*, the velocity field $\bf{u} = (\bf{u}', u^z) = (u^x, u^y, u^z): \mathbb{R}^3\times (0, \infty)\rightarrow \mathbb{R}^3$ and the scalar temperature $T: \mathbb{R}^3\times(0, \infty)\rightarrow \mathbb{R}$ evolve in the domain $\Omega = [0, L_x]\times [0, L_y]\times [0, 1]$ according to

^*see the Boussinesq equations in [3,Section 2.1] and their non-dimensionalization choosing length scales determined by the vertical gap height $h$, time scales defined by the thermal diffusion $\kappa$ and temperature scales determined by the temperature gap $T_b-T_t$.

This is also called Boussinesq system as, in its derivation, density variations generating the vertical buoyancy force appearing in the right-hand side of the velocity equation, are neglected elsewhere in the velocity and temperature equation (Boussinesq approximation). Two dimensionless parameters appear in the system after the non-dimensionalization:

is the (control parameter) Rayleigh number and

is the (material parameter) Prandtl number. Here $\alpha$ is the thermal expansion coefficient, $g$ is the acceleration of gravity, $\nu$ is the kinematic viscosity and $\kappa$ is the thermal diffusivity. The pressure field $p(\bf{x}, t)$ appears as a Lagrange multiplier enforcing the incompressibility condition $(\nabla \cdot \bf{u} = 0)$. In the non-dimensional varibles, the temperature boundary conditions are

For technical convenience we will suppose all functions to be periodic in the horizontal variables $\bf{x}'$ with period $L = (L_x, L_y)$. We refer to ^[2] for details about the derivation of the system.

One of the most important open problems in fluid dynamics is to determine the rate at which heat is transferred in turbulent convection. This theoretical problem finds uncountable many applications in meteorology, oceanography and industry ^[1]. The natural nondimensional quantity for gauging heat transfer in the upward direction is the Nusselt number ${\rm{Nu}}$ defined as the ratio between the convective to the conductive heat flux

where

By similarity law, the Nusselt number ${\rm{Nu}}$ must be a function of the non-dimensional parameters appearing in the system, and thus obeys a functional relation of the type

Theoretical predictions (Malkus ^[4], Kraichnan ^[5], Spiegel ^[6], Siggia ^[7]) and experiments (see ^[8] and references therein) suggest a power-law scaling of the form

where $\alpha$ and $\beta$ depend on the characteristic of the flow. In Section 2 we will illustrate two physically motivated heuristic arguments predicting different scaling behaviors. The challenge for physicists, mathematicians, and engineers is to identify the value of the powers $\alpha$ and $\beta$, as ${\rm{Ra}}$ and ${\rm{Pr}}$ varies. As described by the stability theory, there exists a critical Rayleigh number, ${\rm{Ra}}_c$, under which the pure conduction profile is stable. As Rayleigh number grows over this critical value, convection rolls appears and are eventually destabilized by plumes-like structures arsing from the boundary layers. At very high Rayleigh numbers the dynamics become turbulent. These three phases can be appreciated in Figure 1. As heat transport is enhanced in turbulent regimes, in the whole manuscript we will assume ${\rm{Ra}}\gg 1$. The Prandtl number, instead will be mostly assumed to be large, because of limitation in our analysis. In fact, the smaller the Prandtl number, the stronger is the inertial term in the Navier-Stokes equations.

In this paper we will review some of the major contributions concerning bounds for the Nusselt number in the Rayleigh-Bénard convection problem (1.1)–(1.2) in the last thirty years, distinguishing between no-slip, free-slip and Navier-slip type boundary conditions for the velocity field. Most numerical simulations and rigorous results assume no-slip boundary conditions in the horizontal plates and periodic sidewalls. Whether or not the no-slip conditions are the more suitable for this problem, remains a matter of debate. This is a general and very important question that has been posed by the founders of fluid dynamics and is still subject of ongoing research. Regarding the issue of detecting suitable boundary conditions for the Navier-Stokes equations, Serrin ^[10] writes "The conditions to be satisfied at a solid boundary are more difficult to state, and subject to some controversy. Stokes argued that the fluid must adhere to the solid, since the contrary assumption implies an infinitely greater resistance to the sliding of one portion of fluid past another than to the sliding of fluid over a solid. Another and perhaps stronger argument in favor of adherence, at least for the case of liquids at ordinary conditions, is found in experiments with tube viscometers where the fourth power of the diameter law is conclusively verified. Although these facts are quite convincing when moderate pressures and low surface stresses are involved, they do not apply in all cases; indeed in high altitude aerodynamics an adherence condition is no longer true". Many works took inspiration from Serrin's observation. Here we mention ^[11], where the authors addressed theoretical questions related to the choice of boundary conditions.

In the last decade, free-slip boundary conditions for the Rayleigh-Bénard convection problem have been considered and rigorous upper bounds on the Nusselt number indicate the possibility of enhanced heat transport compared to the no-slip case. But which one of those boundary conditions better represent the physical picture? On the one hand, as remarked in ^[12], (especially) in "turbulent regimes" the boundary conditions might not be perfectly no-slip. On the other hand, the mathematically interesting free-slip boundary conditions are non-physical because of the absence of vorticity production at the boundary, as we will see later in Section 4. H. Navier in 1824 ^[13] suggested the boundary conditions

where $D$ is the rate of strain tensor, $\bar t$ and $\bar n$ are the tangential and normal vector to the boundary and $\alpha > 0$ is the friction coefficient. These conditions go by the name of Navier-slip boundary conditions and can be interpreted as an interpolation between no-slip and free-slip boundary conditions. Since, as observed in ^[14], very few surfaces are truly no-slip or free-slip, the Navier-slip boundary conditions appear reasonable and have been studied in a multitude of applied and theoretical problems, among which we want to mention the works on vanishing viscosity limits for the two-dimensional Navier-Stokes equations ^[15,16].

In this manuscript we will focus on

$1)$ No-slip boundary conditions:

$2)$ Free-slip boundary conditions:

$3)$ Navier-slip boundary conditions:

where $\rm{\mathcal{L}_s}$ is the slip length.

For these three types of boundary conditions we will critically review results concerning rigorous bounds on the Nusselt number and highlight differences in the analyses. As we will see later, while the model equipped with no-slip boundary conditions for $\bf{u}$ is most commonly investigated, the other two types of boundary conditions and their role have been less explored. The goal of this review is to show how the change of boundary conditions affects heat transport properties of the fluid.

The review is organized as follows: In Section 2 we will go through two significant physical arguments predicting different scaling behavior for the Nusselt number. Moreover, starting from definition (1.3) and using a-priori estimates, we will obtain other identities involving the Nusselt number. Depending from the strategy adopted to produce bounds on the Nusselt number (introduced in the third subsection), each identity will later be used singularly or combined. Section 3 is divided in two parts: in Subsection 3.1 we will review some important contributions for the infinite Prandtl number system. In particular we will introduce the background field method and discuss its success and limitation. Then we will present an alternative argument proposed by Constantin and Doering in 1999 and show how this led Otto and Seis in 2011 to produce the optimal (so far) upper bound on the Nusselt number for no-slip Rayleigh-Bénard convection at infinite Prandtl number. In Subsection 3.2 we will go through three arguments to derive bounds for the Nusselt number in the finite Prandtl number case. Rigorous results for the system complemented with the other types of boundary conditions considered in this paper, namely free-slip and Navier-slip boundary conditions, are reviewed in Section 4. In all sections we will examine results in a (partial) chronological order and try to guide the reader through the reasoning behind development of new strategies and, eventually, optimal results.

2.   Preliminaries

2.1.   Physical scalings

In the fifties Malkus ^[4], considering fluids with very high viscosity, performed experiments in which he noticed sharp transitions in the slope of the ${\rm{Nu}}- {\rm{Ra}}$ relation and suggested the scaling

for very high Rayleigh numbers by a marginal stability argument. This is based on the assumption that the heat flux is limited by transport across the (emergent) thermal boundary layer. We illustrate Malkus' argument here. Since the main temperature drop happens near the boundary, we can assume that in a bottom boundary layer of thickness $\delta$ (to be determined) the temperature drops from $1$ to its average $\frac 12$. Thanks to the average $\langle\cdot\rangle$, we may extract the Nusselt number from the boundary layer, where by the no-slip boundary condition we have $(\bf{u} T-\nabla T)e_z\approx \partial_z T$ so that ${\rm{Nu}}\sim \frac{1}{\delta}$. We think of the boundary layer as a pure conduction state in the interval $0\leq z\leq \delta$. Marginal stability refers to the assumption that this state is borderline stable, meaning that its Rayleigh number is critical, which in view of the definition of the Rayleigh number, means

from which, because ${\rm{Ra}}_c\sim 1$, we infer $\delta\sim {\rm{Ra}}^{-\frac 13}\, .$ Inserting this in the scaling of Nusselt number above one finds

The same conclusion can be achieved by rescaling the equations according to

In this way, neglecting the inertial term in the Navier-Stokes equations, we end up with the parameter-free system

Since for the latter system, it is natural to expect that the heat flux is universal, i.e., $\widetilde{ {\rm{Nu}}}\sim 1$, we also obtain ${\rm{Nu}}\sim {\rm{Ra}}^{\frac{1}{3}}.$ Most recently Malkus's theory was refereed to as classical theory ^[1].

Kraichnan in ^[5] and Spiegel ^[6], instead, proposed the scaling

based on the assumption that the heat flux is limited by the transport of the fluid across the bulk. Recalling that, by definition, the convective and conductive heat flux are $q_{\rm{conv}}\sim \rho vc\Delta T$ and $q_{\rm{cond}}\sim \frac{\rho c \kappa \Delta T}{h}$ respectively and assuming

then, the Nusselt number is

We can obtain the same scaling by neglecting the diffusivity and the viscosity term in the model

and rescaling according to

Then, imitating the previous argument, the Nusselt number for the parameter-free system

is

The Kreichnan-Spiegel theory was recently referred to as ultimate theory ^[1].

In ^[6], Spiegel writes "What is suggested by these arguments is that at sufficiently high ${\rm{Ra}}$, turbulent breakdown of the thermal boundary layer occurs and causes a transition from (2.1) to (2.3). No such transition has been detected experimentally, but this is presumably explained by the limitation of the experiments to what in stellar terms are modest Rayleigh numbers. The need to confirm (or deny) (2.3), which is intimately connected with basic ideas of stellar structure theory, poses a great challenge to the experimentalists.". Malkus's scaling ${\rm{Nu}}\sim {\rm{Ra}}^{\frac 13}$ has been been confirmed by experiments at (relatively) high Prandtl numbers (cf. ^[8] for a list of experimental results) while (up to today) it remains no indication of the Spiegel-Kraichnan scaling ${\rm{Nu}}\sim {\rm{Pr}}^{\frac12} {\rm{Ra}}^{\frac 12}$.

2.2.   Other representations of the Nusselt number

The Nusselt number defined in (1.3) admits other useful representations. As we will see later in the preliminaries, their employment depends on the method chosen to bound the Nusselt number. First, let us notice that if the initial temperature $T(\bf{x}, 0) = T_0(\bf{x})$ is chosen such that $0\leq T_0(\bf{x})\leq 1$, then, by the standard maximum principle for parabolic equations, we have

The first representation can be obtained by taking the long-time and horizontal average of the temperature equation and using $\frac 12 \langle\partial_tT^2 \rangle = 0$ as a consequence of the maximum principle, we have

which implies that for all $z\in (0, 1)$

Recalling the definition of Nusselt number in (1.3), this identity yields

Therefore the Nusselt number is independent of $z$, and it is the same for each horizontal layer of fluid. Another representation is derived by testing the temperature equation with $T$ and integrating by parts. In fact, using the incompressibility together with the temperature boundary conditions and (2.5) we obtain

Instead, formally testing the Navier-Stokes equations with $\bf{u}$ and integrating by parts we have

Note that this energy balance (formally) holds for all boundary conditions considered in this article (i.e., (1.5)–(1.7)) as it uses only that $u^z = 0$ at $z = \{0, 1\}$ and the incompressibility condition.

If the energy of the velocity remains bounded in time, averaging in (2.8) we obtain

This identity produces yet another useful representation of the Nusselt number

Clearly this computation is rigorous in two dimensions, but only formal in three dimensions. In this last case, defining Leray-solutions we will obtain an energy inequality in (2.9).

2.3.   Methods to bound the Nusselt number

In this section we briefly introduce the reader to the most common techniques used to bound the Nusselt number. Some of these methods will be carefully explained and used in the next sections.

The first rigorous result was proven by Howards ^[17] by transforming the problem of finding upper bounds on ${\rm{Nu}}$ into a variational problem. Howard found

under the hypothesis of "statistical stationarity", meaning that horizontal averages are time-independent^†. Later Busse extended Howard's result to solutions with multiple boundary-layer structure ^[18] obtaining the same bound with an improved constant prefactor.

^†this holds for stationary flows.

The background field method, proposed (in this context) by Doering and Constantin in ^[19], consists in decomposing the temperature field $T$ as

where $\tau$ is a steady background profile satisfying the driven boundary conditions

and $\theta$ represents the temperature fluctuations satisfying

The fluctuations evolve according to

and the using the representation (2.7), the Nusselt number can be rewritten as

where, for fixed $\tau$, $\mathcal{Q}^{\tau}$ is a quadratic functional in $\theta$ and $\bf{u}$. In this way we can transform the problem of finding the optimal upper bound on the Nusselt number into a variational problem: find the background profile $\tau$ for which $\mathcal{Q}^{\tau}\{\theta, \bf{u}\}\geq 0$ for all $\theta$ and minimizes the integral $\int_0^1|\tau'|^2\, dz$. The background field method will be carefully described in Section 3.1.1 and we refer to ^[20] for a recent review on it.

A generalization of the background field method was proposed in ^[21] to improve the available bounds. The method is based on the following observation: if $V(t) = V\{\theta(\cdot, t), \bf{u}(\cdot, t)\}$ is a uniformly in time bounded and differentiable function with $\lim_{t\rightarrow \infty}\int_0^t\frac{d}{ds}V(s)\, ds = 0$ then, in order to prove

it suffices to show

where

In fact, time-averaging this inequality we recover (2.14). It is easy to see that the auxiliary functional method and the background field method coincide when the auxiliary functional $V$ is quadratic in $\bf{u}$ and $\theta$ ^[22]. The power of this method is to allow the functional $V$ to be more general (than the quadratic form $\mathcal{Q}^{\tau}$), to the expenses of an increasing analytical complexity which eventually requires computer-assisted investigation ^[23,24]. The auxiliary functional method has proven to be successful to produce sharp bounds for (system) of ordinary differential equations ^[25], but the application of this method for producing optimal bounds on the Nusselt number in the Rayleigh-Bénard convection can only be speculated.

Another method to produce upper bounds based on the vertical localization of the Nusselt number was first used in ^[26] and then proposed by ^[27] as an alternative method to the background field method for a variety of fluid dynamics problems. This method is based on the observation that (2.6) can be averaged over any horizontal layer of fluid; in particular, averaging near the bottom boundary layer $0\leq z\leq\delta$ we are led to the definition

^‡In particular, because of the boundary conditions on the temperature and the fact that $T\geq 0$ for every $\bf{x}, t$, by the maximum principle, we have

^‡Notice here that, if $|\int_0^t\iint u^zT\, d \bf{x}'\, ds| < g(z)\in L^1$, then the switch in order of the $\langle\cdot \rangle$ averages and the integral in $z$ is justified by the Lebesgue dominated convergence theorem. The condition will be generally satisfied thanks to the a-priori (energy) estimates for $u$ and $T$.

The goal is now to exploit the equations of motion to derive regularity bounds on $u^z$ and $T$ in the strip $[0, L_x]\times[0, L_y]\times[0, \delta]$. As we will see later in Section 3.2, identity (2.9) will be particularly useful, when adopting this method. Compared to the background field method, this approach has the advantage to be free from the rigid constraint imposed by the variational method. Nevertheless in Section 3.1.1 we explain that these two methods coincide for a particular choice of background profile.

3.   No-slip boundary conditions

3.1.   Infinite Prandtl number

For some fluids the Prandtl number is so high (see ^[7], table 1) that the effect of inertia in the Navier-Stokes equations is almost negligible. This motivates the study of the infinite-Prandtl number model, obtained by setting ${\rm{Pr}} = \infty$ in the momentum equation in (1.1)

Using the incompressibility condition it is easy to see that (3.1) can be rewritten as the fourth-order elliptic equation

where $\Delta_{{{\rm{H}}}} = \partial_x^2+\partial_y^2\,$ and the pressure is eliminated. Notice that by the incompressibility condition we infer that $\partial_z u^z = -\nabla_H\cdot \bf{u}' = 0$ at $z = \{0, 1\}$ since $\bf{u}' = 0$ at $z = \{0, 1\}$. From relation (3.2), it is now clear that the temperature is "instantaneously" slaved to the velocity field. Xiaoming Wang in ^[28,29] reasonably asked whether this reduction is justified and whether the statistical properties of the flow in the ${\rm{Pr}} = \infty$ model and in the finite (but large) ${\rm{Pr}}$ number model can be close in some sense. The author proved that, indeed, the global attractor of the Boussinesq system at large Prandtl number converges to that of the infinite-Prandtl number model, giving some indication on the possible closedness of the respective statistical properties.

3.1.1.   Story of the background field method in five acts

Introduction of the background field method: The first application of the background field method to find upper bounds on the Nusselt number for the Rayleigh-Bénard convection problem appears in the Doering and Constantin 1996 paper ^[19]. We reproduce their arguments, adapting them to the infinite Prandtl number case. We start by splitting $T$ according to (2.10). Testing equation (2.12) with $\theta$ and integrating by parts in $\Omega$ we obtain

where we used incompressibility together with the condition $u^z = 0$ at $z = \{0, 1\}$ to eliminate the term $\frac{1}{2}\int_{\Omega}u\cdot\nabla\theta^2$; then, averaging we get

The combination of this identity with

yields

Thanks to this new representation of the Nusselt number we can formulate the following

Criterion 3.1. If $\tau:[0, 1]\rightarrow \mathbb{R}$ satisfies

for all fluctuations $\theta(\bf{x}', z)$ satisfying (2.11) and related to $u^z$ through (3.2), then

Condition (3.6) is often referred to as "stability constraint" or "spectral constraint" of the background profile. The strategy just described goes by the name of background field (or background flow) method. Already introduced in 1992 by Constantin and Doering to bound the energy dissipation in shear driven turbulence ^[30], this method has been extensively used for estimating the Nusselt number. In fact, an upper bound is immediately obtained by computing the Dirichlet energy of a stable (in the sense specified above) background profile. Thus, each profile satisfying (3.6), produces an upper bound for the Nusselt number.

The first example of a marginally stable profile is (a smooth approximation of)

plotted in Figure 2, for which one easily shows ${\rm{Nu}}\lesssim {\rm{Ra}}^{\frac 12}$.

We use this example to explain the adjective "marginal" in the stability condition. In analogy to the Malkus' argument reported in Section 2, this refers to the fact that the boundary layer $[0, \delta]$ adjusts itself to be "borderline" stable. We are going to see that the thickness of $\delta$ determines the validity of (3.6) and eventually the bound on the Nusselt number. The bigger $\delta$ is, the smaller the upper bound on the Nusselt number. Obviously there will be artificially thin boundary layers that will ensure the stability condition, but for optimality of the upper bound, we want to detect the thickest stable boundary layer. Any boundary layers thicker than this, will be unstable.

Now we illustrate the idea behind the proof of the marginal stability of profile (3.8), as this will set the basis for future computations. Applying the Cauchy-Schwarz inequality we have

Thanks to the boundary conditions for $u^z$ and $\theta$, we can use the Poincaré estimates

together with (2.7) and (2.9) to get

Optimizing on the boundary layers thickness, we select

yielding

Therefore we proved

Theorem 3.2. For solutions of (1.1)–(1.2) with no-slip boundary conditions (1.5) the following upper bound holds:

We remark that, in order to obtain this bound we could equivalently use the localization principle in (2.15), thus avoiding going through the background field method. In fact by the Cauchy-Schwarz inequality and estimates (3.9) and (3.10), we obtain

Balancing, the optimal $\delta$ satisfies

leading to

Therefore the localization method is equivalent to the background field method with the choice (3.8). More discussions and interesting observations about the relation between the background field method (and auxiliary functional method) and the localization method may be found in ^[22].

Blooming of the background field method: In ^[31], Doering and Constantin choose (a smooth approximation of) the profile

plotted in Figure 3, and improve Theorem 3.2 by employing finer estimates for the vertical velocity $u^z$. Here we describe their argument.

Thanks to periodicity in the horizontal directions, passing to Fourier-series offers a big advantage. In fact, the unknown functions can be written as Fourier series with $z-$ dependent coefficients^§

^§Here we assume, for simplicity $L = L_x = L_y$.

Then the equation for (3.2) decomposes in

where $k = \frac{2\pi}{L}|\bf{j}|$ and the boundary conditions for $u^z$ are

Decomposing $\mathcal{Q}^{\tau}$ mode by mode in the translation invariant horizontal directions, it is clear that $\mathcal{Q}^{\tau}$ will be non-negative, if for each horizontal wavenumber $k$

is non-negative. Here the symbol $^{*}$ denotes the complex conjugate. Because of the choice of the background profile in (3.11) we have

The crucial element leading to the new bound is contained in the way the authors estimate the mixed term. To illustrate the argument, let us focus on the lower boundary layer. By Poincaré's inequality applied to $u^z$, $\theta$ (see (3.9) and (3.10)) and $\partial_z u^z$ (thanks to the fact that $\partial_z u^z = 0$ at $z = \{0, 1\}$) we have

Again, due to the boundary conditions $u^z = 0$ and $\partial_z u^z = 0$ at $z = \{0, 1\}$, the $L^{\infty}$ norm of $u^z$ can be bounded by using the following interpolation estimate

Lemma 3.3. Let $v(z):[0, 1]\rightarrow \mathbb{R}$ be a smooth function satisfying $v = 0$ and $\partial_z v = 0$ at $z = \{0, 1\}$, then

Applying this lemma to $v = u^z$ and using Young's inequality we obtain

In order to conclude the argument, we now reconnect to the equation. In fact, thanks to the instantaneous slaving in (3.12) one can easily deduce

for some numerical constant $C > 0$. In fact, squaring (3.12) and integrating by parts we have

Choosing $a = \frac{7}{10}$, we obtain the desired bound. Finally going back to (3.15) we have

Then, from (3.14)

The choice $\delta\sim {\rm{Ra}}^{-\frac 25}$ in (3.7) yields

Theorem 3.4. For solutions of (1.1)–(1.2) with no-slip boundary conditions (1.5) the following upper bound holds:

From the previous argument we see that the choice of bounding the mixed term with the $L^{\infty}-$norm of the second derivative of $u^z$ "maximizes" the power of $\delta$. In fact this last method gives us a $\delta^{\frac 52}$ prefactor against the $\delta^2$ prefactor of the previous section. Moreover notice that in this proof, the instantaneous slaving of $u^z$ to $\theta$ given by (3.2) was crucial. It is now natural to ask whether the profile chosen in (3.11) is optimal. Before jumping to the next two seminal results, let us comment on the background profile chosen in (3.11): The temperature field is supposed to be laminar in very thin boundary layers and this justifies the linear profile, decreasing from $z = 0$ to $z = \delta$ and from $z = 1-\delta$ to $z = 1$. In the bulk instead the temperature is not expected to vary wildly.

Although the profile in (3.11) turned out to be marginally stable for some choice of very small $\delta$, it does not reflect the idea that the fluid layers must be "stably stratified". This means that the lighter, warmer parcels of fluids are expected to be on the top of the cold, denser parcels. This observation leads us to the next result.

Almost crowning of the background field method: In ^[32], Doering, Otto and Westdickenberg (neé Reznikoff) show stability of the logarithmic background profile

with

In this case, the derivative of the background profile $\tau$ produces weighted terms which are subtler to control, as we are going to describe now. With the choice profileserrin plotted in Figure 4, the quadratic form becomes

Adding and subtracting $(\frac{\lambda}{z}+\frac{\lambda}{1-z}) \langle u^z\theta \rangle$ in the boundary layers we rewrite

where

and

Smuggling in the weight $z^2$ and applying Hölder and Cauchy-Schwarz inequalities, the last term in $\mathcal{Q}^{\tau}_{\rm{B}}$ can be estimated as follows:

By the fundamental theorem of calculus

and, by a direct computation

The crucial estimate in this paper is contained in the following

Lemma 3.5. For $\theta\in L^2(\Omega)$ and $w$ satisfying (3.12) and (3.13), the bound

holds.

The combination of this lemma applied to $w = u^z$ with the estimates above yields

Arguing in the same way for $\mathcal{Q}^{\tau}_{\rm{T}}$ and imposing $\delta^3\ln((1-\delta)/\delta) = \frac{1}{2 {\rm{Ra}}}$, one finds that $\delta\sim \frac{1}{ {\rm{Ra}}\ln {\rm{Ra}}}$ is optimal, and, since ${\rm{Nu}}\leq \int_0^1(\tau')^2\, dz\sim \frac{1}{\delta}$, this proves the following

Theorem 3.6. For solutions of (1.1)–(1.2) with no-slip boundary conditions (1.5) the following upper bound holds:

This result improves Theorem 3.4 and was the first one capturing (with an upper bound) Malkus' scaling ${\rm{Ra}}^{\frac 13}$ up to a logarithmic correction only using the background field method. The choice of the logarithmic background reveled itself to be "optimal" in a sense that we will specify at the end of this section. Roughly speaking, the logarithmic profile respects the stable stratification since it (monotonically) grows in the bulk. And the slow growth in the bulk maintains stability.

Improving the (logarithmic) failure: The logarithmic correction in (3.22) was later improved by Otto and Seis in ^[26] using the following idea: going back to (3.19) and focusing on (3.20), the authors proved that, for $(\lambda {\rm{Ra}})^{-\frac 13}\leq1$ the upper bound

holds. Estimate (3.23) is crucial for reducing the power of $\lambda$ in the bound: in fact, applying it in the estimate

where we used the Poincaré-type inequality

the authors then conclude

Choosing $\varepsilon = \frac 14$ and using $\int_0^1\frac{ \langle\theta u^z \rangle}{z}\, dz\geq 0$ by (3.21), we have

The optimal $\delta$ is

and, as a result, we showed

Theorem 3.7. For solutions of (1.1)–(1.2) with no-slip boundary conditions (1.5) the following upper bound holds:

The new estimate for the velocity in (3.23) improves the logarithmic correction substantially, but the bound still fails to reproduce the classical scaling ${\rm{Ra}}^{\frac 13}$. It is now natural to ask

● can we find a "better" profile which captures the Malkus' scaling without corrections? or

● is it possible to refine the regularity estimates to further reduce the logarithmic correction?

We will (partially) answer these questions in the next sections.

Limitations of the background field method: As we have seen in the previous section, the background field method transforms the problem of finding bounds on the Nusselt number into a variational problem. If we would find the minimizer of this problem, we would get the optimal upper bound that this method is able to produce. In particular, the shape of the optimal $\tau$ would give us precise information about the (statistical) behavior of the temperature.

Let us rescale the infinite Prandtl number system according to the rescaling (2.2) suggested by Malkus' marginal stability argument, and set $H = {\rm{Ra}}^{\frac 13}$ to obtain

Let us define the Nusselt number associated to the background flow method

with $\mathcal{Q}^{\tau}$ being

where $\theta$ is instantaneously "slaved" to $u^z$ through

Clearly,

and the challenge is to understand how far the Nusselt number ${\rm{Nu}}$ is from the number ${\rm{Nu_{bf}}}$ produced by the background field method.

In ^[33] the authors prove

Theorem 3.8. For solutions of (1.1)–(1.2) with no-slip boundary conditions (1.5), ${\rm{Pr}} = \infty$ and ${\rm{Ra}}\gg 1$ the following upper bound holds:

This ansatz-free lower bound is proven by extracting local information on $\tau$ from the non-local stability condition (3.25), which written in Fourier-series reads:

with ${\bf{k}}\in \frac{2\pi}{L}\mathbb{Z}^{d-1}\setminus\{0\}$ and

Notice that in the quadratic form only the modulus of $k$ appears, so that the stability condition is independent of the dimension. The authors initially show the result in a simplified setting. First, considering profiles $\tau$ that are stable even under perturbation that have horizontal wave-length much larger than $H$, letting the side-length $L\rightarrow \infty$, implying $k\in \mathbb{R}$. Second, reducing the stability condition to the second term in (3.25), which is the leading order term in the stability condition. Under these assumptions the authors characterize all stable profiles:

Proposition 3.9. If $\frac{d\tau}{dz}$ satisfies

for all $k\in \mathbb{R}$ and all $\hat{ u^z}$ satisfying (3.27), then

and

This means that, if $\tau$ satisfies the reduced stability condition (3.28), then it must be increasing and "logarithmically" growing in the bulk. We remark that, in ^[34] it is proved that the result above still holds true if the lateral size $L$ is of order $H$. Extending this result to the full stability condition requires some careful analysis: by subtle averaging of the stability condition, the authors in ^[33] construct a non-negative convolution kernel with the help of which they can express the positivity on average approximately in the bulk. Also the logarithmic growth can be recovered at least approximately in the bulk. To finally achieve the lower bound, further estimates connecting the bulk with the boundary layers are provided. This result proves the optimality of the logarithmic profile.

Interestingly, by a mixture of analytical and numerical analysis, Ierely, Kerswell and Plasting in 2006 ^[35] anticipated (although slightly underestimating) the result in ^[33], showing ${\rm{Nu}}\sim {\rm{Ra}}^{0.33173}(\ln {\rm{Ra}})^{0.0325}$.

3.1.2.   The Constantin & Doering '99 argument and the best upper bound at infinite Prandtl number

In 1999 Constantin & Doering ^[36] caught the Malkus scaling ${\rm{Ra}}^{\frac 13}$ (up to a logarithm) without optimizing on the background profiles, but extracting finer information from the equation for the velocity. This paper marks the beginning of a new and powerful way of thinking about this problem. The authors choose $\tau$ to be (a smooth approximation of) the profile

in

An application of the maximum principle $\|\theta\|_{L^{\infty}(\Omega)}\leq 1$ to estimate the mixed term

yields

where the bound $|\tau'(z)|\lesssim \frac{1}{\delta}$ was used.

As we have seen in Section 3.1.1, the interpolation estimate (3.16) was the crucial ingredient for proving the bound ${\rm{Nu}}\lesssim {\rm{Ra}}^{\frac 25}$, but turned out to be sub-optimal. In this paper the authors proved the much finer estimate

Proposition 3.10. Suppose $w$ solves the problem

with horizontally periodic boundary conditions. Then, for any $\alpha\in (0, 1)$ there exists a positive constant $C_{\alpha}$such that the bound

holds, where

This maximal regularity estimate was obtained by decomposing the operator $B = \partial_z^2(\Delta^2)^{-1} \Delta_{{{\rm{H}}}}$ into the sum of non translationally invariant operators, whose kernel is not explicit. Inserting the long-time average of (3.31) into (3.30) and using

to estimate the $C^{0, \alpha}$ norm in the logarithmic correction further^¶, the authors obtain

^¶This estimate is simply obtained by testing (2.12) with $\Delta \theta$, integrating by parts, and using the Cauchy-Schwartz and Youngs inequalities together with the interpolation inequality $\|\nabla \theta\|_{L^4(\Omega)}^2\leq C\|\theta\|_{L^{\infty}(\Omega)}\|\Delta \theta\|_{L^2(\Omega)}\,$ and $\|\nabla \bf{u}\|_{L^2(\Omega)}^2\leq C {\rm{Ra}}^2$, where the last inequality is derived by testing the Navier-Stokes equations, integrating by parts and using the maximum principle for the temperature.

and the choice $\delta\sim {\rm{Ra}}^{-\frac 13}(1+\log_+ {\rm{Ra}})^{-\frac 23}$ yields

We notice that the logarithmic correction in this bound is bigger than the one produced in ^[26] using the logarithmic profile in the background field method. Nevertheless, we are now going to discuss that, by a refinement of this argument, Otto and Seis in ^[26] proved the best upper bound (up to today). We first observe that, by the maximum principle and Poincaré's inequality (applied twice in $z$) the bound

holds and is the same as (3.30) (which was obtained by (2.15)). The key of the result in ^[26] is the estimate on the hessian of $u^z$

obtained by a clever combination of the maximal regularity estimate

and the weighted estimates in (3.21). Inserting (3.33) into the bound for the (localized) Nusselt number, one gets

Using Young's inequality we are led to

Choosing

then

yielding^||

^||Here one uses the fact that $XY\ln^{-\frac 13}(XY)\leq Y$ implies $XY\leq Y\ln^{\frac 13}(Y)$.

Theorem 3.11. For solutions of (1.1)–(1.2) with no-slip boundary conditions (1.5) the following bound holds:

Besides the improvement of the logarithmic correction, the interesting aspect of this result is that the background field method has not been used. That is, this bound relies only on the maximum principle for the temperature and on regularity estimates derived from Stokes equation. We notice that the authors argue that it is sufficient to establish (3.34) for (3.2) with the free-slip boundary conditions

for details see Proof of Lemma 4 in ^[26]. Passing to this boundary conditions is fundamental in order to extend the problem to the whole space by periodicity and therefore being able to use Fourier series techniques in both (horizontal and vertical) variables^**. In Fourier variables, the hessian of $u^z$ can be written as

^**this time the symbol $\hat{\cdot}$ denotes the Fourier transform in $\bf{x} = (\bf{x}', z)$. With $k$ and $\xi$ we denote the horizontal and vertical Fourier variable, respectively.

and, from this representation, mimicking a Littlewood-Paley decomposition, one can analyze the small-intermediate and large length scales separately. This "separation of scales" is crucial for obtaining the desired bound, since the $L^{\infty}-$norm is critical for Calderón-Zygmund type estimates.

3.2.   Finite Prandtl number

In this section we review some results concerning upper bounds on the Nusselt number at finite Prandtl number. Differently from the previous section, in this case the temperature $T$ and the velocity $u^z$ are no longer instantaneously slaved and the equation for $u^z$ can no longer be written as a fourth order boundary value problem as in the case ${\rm{Pr}} = \infty$ (recall (3.2)). As a consequence, we will see that the application of the background field method becomes somewhat unnatural and new approaches have to be considered. In order to present the results concerning the three-dimensional case, we will sometimes implicitly assume regularity of the Navier-Stokes equations. Nevertheless the arguments can be made rigorous by working with Leray weak solutions, i.e., assuming

3.2.1.   Background field method approach

Assuming regularity for the solution of the Navier-Stokes equations, adding the balances (2.8) and (3.3) we have

Since $\bf{u}$, $T$ (and therefore $\theta$) stay bounded in $L^2$ for all time,

we may take the long-time limit obtaining

for some $a > 0$. Here we used that

since $\langle u^z \rangle_{\rm{H}} = 0$ by the incompressibility condition. Combining this identity with (3.4) we obtain the new representation of the Nusselt number

Clearly this coincides with (3.5) when $a = 0$. If we can select a background profile $\tau$ with $\tau(0) = 1, \tau(1) = 0$ such that the quadratic form

is positive definite for all $\theta$ satisfying (2.12), then an upper bound is given by

In ^[19], Doering and Constantin choose the following background profile

Inserting this choice of $\tau$ (see Figure 5) in the quadratic form with $a = 1$, and focusing on the mixed term we get

By Cauchy-Schwarz, Poincaré and Young's inequality we obtain

and the same argument applies to the top boundary layer. Then the quadratic form is bounded from below by

So, clearly, if $\delta$ is chosen to be

then

In particular, this shows

Theorem 3.12. For solutions of (1.1)–(1.2) with no-slip boundary conditions (1.5) it holds

We notice that in this argument neither the incompressibility, nor the full set of boundary conditions was exploited. Therefore this result holds for all boundary conditions such that $u^z = 0$. In order to go beyond this result, the system have to be exploited extensively.

3.2.2.   A dynamical systems approach

In ^[29], Wang shows that the global attractors of the infinite-Prandtl and the large-Prandtl number convection model remain close, by viewing the Boussinesq system as a small perturbation of the infinite Prandtl number model. Subsequently he derived bounds on the Nusselt number at finite (but large) Prandtl number starting from the Constantin and Doering '99 result in ^[36]. Following up on his previous works ^[28,29,37], the central idea in Wang's result in ^[38] is to view the Navier-Stokes equations as a perturbation around the stationary Stokes equation,

where $f = {\rm{Ra}} T \bf{e_z}-\frac{1}{ {\rm{Pr}}}(\partial_t \bf{u}+ \bf{u}\cdot \nabla \bf{u})$ and apply the bound in ^[36] (reviewed in Section 3.1.2).

Using (3.5) and inverting the Stokes operator we have

Following the argument in Section 3.1.2 (applied to the Stokes equation), choosing the background profile (3.29) then

provided

A new argument is needed to estimate the last two terms in (3.38). In fact, for the first term, using the profile (3.29) together with the Cauchy-Schwarz inequality, Poincaré's estimate (twice, in $z$) and Stokes regularity one finds

for some constant $c = c(\Omega) > 0$. Similarly, using the maximum principle, the bound

holds, where the spatial $L^2-$norm of the nonlinear term was estimated by the Poincaré estimate and Agmon inequality in three dimensions, i.e.,

The combination of these estimates in (3.38) yields

where we used again the Poincaré inequality. Finally, inserting the choice of $\delta$ in (3.39) and the (a-priori) upper bounds

in (3.40), Wang obtains

Theorem 3.13. For solutions of (1.1)–(1.2) with no-slip boundary conditions (1.5) the following upper bound holds:

We remark that the crucial (a-priori) upper bounds used to derive the final estimate on Nusselt number, were obtained by energy and enstrophy inequalities. In particular a Gronwall-type argument is needed to ensure the existence of an absorbing ball such that

under the large-Prandtl number condition $\Pr\gtrsim {\rm{Ra}}$.

3.2.3.   Maximal regularity approach

The result in ^[39] improves Wang's result by perturbing around the nonstationary Stokes equation, i.e., considering

and proving a maximal regularity estimate of the type

where the norm $\|\cdot\|$ is strong enough to control the nonlinear term

and sufficiently weak so that

The interpolation norm, which is able to grant this controls, is

In fact, the $L^1$-weighted norm appears natural to bound the nonlinearity:

where we used the Cauchy-Schwarz inequality together with Hardy's inequality and the dissipation bound in (2.9). Since the interpolation norm is critical for maximal regularity estimates (since the $L^{\infty}$ and the $L^1-$weighted norm are both critical with this respect), the authors proved the following

Proposition 3.14. Suppose $\bf{u}, p$ and $\bf{f}$ satisfy

and $\bf{f}$ is horizontally banded, i.e.,

Then

holds.

As we can see, the authors do not obtain the full maximal regularity estimate due to cancellations that occurs as an effect of using the bandedness assumption. Nevertheless the crucial estimate

is not affected and allows to conclude the argument. In fact, if only horizontal intermediate wavelengths could be taken into account, then, inserting the maximal regularity estimate in the bound

and choosing $\delta\sim ((\frac{ {\rm{Nu}}}{ {\rm{Pr}}}+1) {\rm{Ra}})^{-\frac 13}$, we would be able to conclude

This estimate is unfortunately not achieved as the high and low wavelength will weight in the final estimate. Nevertheless, by the incompressibility condition and the a-priori energy inequality (2.9), the authors show that the region of intermediate wavelengths (where estimate (3.41) holds) is large and that the correction from the desired bound is relatively small. In conclusion, we have

Theorem 3.15. For solutions of (1.1)–(1.2) with no-slip boundary conditions (1.5), the following upper bound holds:

4.   Free-slip and Navier-slip boundary conditions

In the previous sections, the no-slip boundary conditions

^†have been crucial to obtain upper bounds of the type ${\rm{Nu}}\lesssim (\log {\rm{Ra}})^{\gamma} {\rm{Ra}}^{\frac 13}$. In fact by the application of Poincaré's estimate in $z$ we bounded

^††Recall from Section 3.1 that this is an equivalent way of writing the no-slip boundary conditions using incompressibility.

With free-slip and Navier-slip boundary conditions this estimate is no longer available and it is therefore more difficult to appeal to maximal regularity estimates, which turned out to be crucial for producing optimal bounds in the no-slip setting. In the next section we see how to circumvent this problem and get sharp quantitative estimates on the flow. We want to notice here that, in terms of regularity estimates, the free-slip boundary conditions present a big advantage compared to the no-slip boundary conditions. In fact the free-slip boundary conditions (1.6) are equivalent to

since $\partial_z^2 u^z = -\nabla_H\cdot\partial_z \bf{u}' = 0$ thanks to the incompressibility condition. Now one can easily see that the equations can be naturally extended in the vertical direction by periodicity. This observation was used in ^[26] to derive a maximal regularity estimate and in ^[40] to show the existence of a global attractor in each affine space where the velocity has fixed spatial average. One major advantage is that the extension by periodicity allows to use Fourier-series techniques also in the vertical direction. In the no-slip case, this extension was not possible, and, in order to derive maximal regularity bounds we were forced to decompose the operator and exploit Fourier-series techniques only in the horizontal variables.

Due to limitations in the analysis at finite Prandtl number, in the next sections we will often consider the two dimensional Rayleigh-Bénard system (1.1). In this case we will denote with $\bf{x} = (x, z)$ a vector in $\Omega = [0, L]\times[0, 1]$ and with $\bf{u} = (u^x, u^z)$ the velocity field.

4.1.   Free-slip boundary conditions

4.1.1.   Bounds in two dimensions and at finite Prandtl number

In his PhD thesis in 2002 ^[41], Jesse Otero reported numerical indication of the scaling ${\rm{Nu}}\sim {\rm{Ra}}^{\frac{5}{12}}$ for the two-dimensional Rayleigh-Bénard convection problem with free slip boundary conditions. Later in 2011, Doering and Whitehead in ^[42] proved rigorously the following result

Theorem 4.1. For solutions of (1.1)–(1.2) and (1.6) in two-dimensions the following bound holds:

The argument takes advantage of information on the scalar vorticity $\omega = \partial_x u^z-\partial_z u^x$ satisfying

The enstrophy balance is

from which one can show ^[40] that

and obtain

Combining this, with the long-time averages of (3.3) and (2.8), the authors obtain the representation

where

with $b\in (0, 1)$ and $a > 0$. By standard arguments it is not difficult to show that

where $\mathcal{R}$ denotes the real part of a complex number. Since the first two terms of the right-hand side are already non-negative (for some choice of $a$ and $b$), the authors need to argue for the (potentially negative) mixed term. We start with an observation: since

there exists a $z_0\in (0, 1)$ (no matters where it lies in $(0, 1)$!) such that

Then the "missing" derivative (for applying Poincaré's estimate) is partially recovered by the estimate

where we used the Youngs inequality. Then

Combining this estimate with

which holds thanks to the assumptions on $\theta$, we showed

Lemma 4.2. For $w\in H^2(0, 1)$ with $w = 0$ at $z = \{0, 1\}$ and $\zeta\in H^1(0, 1)$ with $\zeta = 0$ at $z = \{0, 1\}$, the bound

holds for some positive constants $c_1$.

We observe now that by using the (2d) relation $\omega = \partial_z u^x-\partial_x u^z$ and the divergence-free condition in Fourier variables, one obtains

and easily find the bound

In particular, (4.5) and (4.6) together with the choice $c_1 = \frac{\sqrt{3}}{6}$ yield

Combining this a-priori bound with Lemma 4.2 applied to $w = \hat{ u^z}$ and $\zeta = \hat{\theta}$, we find

Inserting (4.8) in the lower bound for $\hat{ \mathcal{Q}^{\tau}}$, we immediately see that $\hat{ \mathcal{Q}^{\tau}}\geq 0$ if

Minimizing on $\delta$ (over $k$), choosing $\delta\sim {\rm{Ra}}^{-\frac{5}{12}}$ where $k\sim {\rm{Ra}}^{\frac 14}$ is the minimal wavenumber, we obtain the desired bound ${\rm{Nu}}\lesssim {\rm{Ra}}^{\frac{5}{12}}$.

Recent numerical simulations ^[43] indicate the clear asymptotic scaling

for steady $2d$ convection rolls between free-slip boundary conditions. It is an interesting problem to bridge this gap and see whether the result in ^[42] is optimal or steady $2d$ convection rolls (slightly) inhibit heat transport.

4.1.2.   Bounds in three dimensions

In three dimensions, the vorticity is no longer a scalar and the enstrophy balance (4.3) does not hold. Then the variational problem to study is: find the optimal $\tau$ that minimizes the right-hand side of

constrained to

This is the same as (4.4) with $a = 0$. In ^[44] and ^[35] the authors apply the background field method to the free-slip problem and their numerical simulations indicate

for stably stratified profiles in the bulk, with $\tau' = p > 0$. Choosing this family of background profiles, Doering and Whitehead in ^[12] rigorously proved

Theorem 4.3. For solutions of (1.1)–(1.2) and (1.6) in three dimensions the following upper bound holds:

Confirming the optimality of the family of stably stratified profiles and the numerical result in ^[44]. The argument is almost identical to the one for Theorem 4.1, with the crucial bound (4.7) replaced by

where $\zeta$ imitates the role of vorticity in two dimensions and is referred to as "pseudo-vorticity" (or "pseudo-enstrophy") and it is defined by^‡‡

^‡‡In real space the pseudo vorticity is defined through the relation

Recalling the elliptic fourth-order boundary value problem written in Fourier variables (3.12), the pseudo-vorticity satisfies

An interesting point discussed in ^[12] concerns the relation between the slope $p$ of the profile $\tau$ and the optimal constant $c$ in (4.9). Fist, the authors clarify that, a stratified profile is the bulk or a balance parameter $b$ with a non-stratified profile, produce the same scaling. In fact the authors demonstrate how setting $p = 0$ and optimizing in $b$ produce a prefactor $c$ that is almost indistinguishable from the case in which $b = 0$ and $p$ is optimized. In particular, setting $b = 0$ and $p = \frac{3}{29}$, the authors find ${\rm{Nu}}\lesssim0.2982 {\rm{Ra}}^{\frac{5}{12}}$. As the slope parameter is nearly the same to $p\sim 0.103$ computed by Plasting and Ierley in 2005, the analysis of Whitehead and Doering is rather sharp. Second, the author showed that the best constant ($0.28764$) is achieved for $p = -\frac{1}{17}$ and $b = \frac{5}{12}$. This in particular tells that, if $b > 0$, the (optimal) bound is produced with a non-stably stratified profile (since the slope parameter is negative).

Combining the methods in ^[38] and ^[12], Wang and Whitehead in ^[14] succeeded in treating the three-dimensional finite Prandtl number case as a perturbation of the infinite Prandtl number problem problem, i.e., studying

where $\bf{f} = \partial_t \bf{u}+ \bf{u}\cdot \nabla \bf{u}$ and ${\rm{Pr}}$ is supposed to be (very) large, allowing them to apply the argument in ^[12]. In particular, the authors eliminate the pressure and study the relation

and benefit from rewriting the equation in Fourier space

where $\zeta$ is the pseudo-vorticity defined in (4.10). Modulo the $O(\frac{1}{ {\rm{Pr}}})$ inertial terms, this relationship is identical to (4.12).

Also here, as in Section 3.2.2, the idea is to use the background field method, splitting the quadratic form

into a part (containing $O(1)$ contributions already present at infinite Prandtl number) which is positive definite and a rest, of order $o(\frac{1}{\Pr})$ (and potentially negative) that needs to be controlled with other methods. They proved

Theorem 4.4. Suppose that $\bf{u}, T$ solve (1.1)–(1.2) and (1.6). Then there exists a non-dimensional constant $c_0$ such that if the Grashof number ${\rm{Gr}}$ is sufficiently small, i.e.,

then

Here $c_0$ is a small number determined by ensuring the simultaneous validity of

These are the fundamental bounds necessary to control the $O(\frac{1}{\Pr})$ negative contributions in the quadratic form.

4.2.   Navier-slip boundary conditions

As observed before, while free-slip conditions are considered to be "non-physical" due to the absence of vorticity production at the boundary (i.e., $\omega = 0$ at $z = \{0, 1\}$), the Navier-slip conditions are adapt to describe surfaces where some slip occurs, but there is some stress on the fluid. In fact, in realty most of the surfaces are neither perfectly free slip not no slip ^[12]. If we think of this boundary conditions as interpolating between no-slip and free-slip (as the slip length $\rm{\mathcal{L}_s}$ increases from $0$ to $\infty$), we expect that the bounds on Nusselt number will interpolate as well, between the ${\rm{Ra}}^{\frac 12}$ and the ${\rm{Ra}}^{\frac{5}{12}}$ bounds. In two dimensions, the vorticity satisfies

and the enstrophy balance is

Notice that as $\rm{\mathcal{L}_s}\rightarrow \infty$, this balance reduces to (4.3). Following an idea in ^[16] applied to the two-dimensional unforced Navier-Stokes equations the authors in ^[45] prove

Lemma 4.5. Let $\omega$ satisfy (4.14), $\rm{\mathcal{L}_s} \geq 1$ and $p \in [1, \infty)$. Then

Thanks to this lemma, combining the averaged version on (4.15) with the balance

we can write a new representation of the Nusselt number:

where

with parameters $a > 0$, $b\in(0, 1)$ and $M > 0$ to be chosen at the end. The combination of the arguments in ^[42] (to bound the mixed term $\langle\int_0^1\tau' u^z\theta\, dz \rangle$), Sobolev trace inequalities together with the regularity estimate on the pressure contained in

Proposition 4.6. The pressure in (1.1) satisfies

and for any $r\in (2, \infty)$ the estimate

holds.

yield the lower bound for $\mathcal{Q}^{\tau}$

Optimizing in $a, b, M$ and choosing $\delta\sim {\rm{Ra}}^{-\frac{5}{12}}$, Drivas, Nguyen and Nobili in ^[45] achieve the following result:

Theorem 4.7. Let $u, T$ solve (1.1)–(1.2) with Navier-slip boundary conditions (1.7). Then

Moreover, for the case ${\rm{Pr}} = \infty$ in dimension $d\geq 2$, the authors show that the bound

follows from the same arguments. Using similar techniques to the one in ^[12], this upper bound was independently (and previously) noticed by Whitehead in unpublished notes.

5.   Conclusions and open problems

In this paper we reviewed a selection of contributions in the field of quantitative bounds on the Nusselt number for the classical Rayleigh-Bénard convection problem between fixed impermeable horizontal boundaries. In the following chart we summarize the best upper bounds available, distinguishing between no-slip, Navier-slip and free-slip boundary conditions.

There are various questions and open problems that remain open:

Ultimate scaling: ${\rm{Ra}}^{\frac 12}$ or ${\rm{Ra}}^{\frac 13}$? Focusing on the no-slip boundary conditions, we see that, on the one hand, upper bounds of the type ${\rm{Nu}}\lesssim {\rm{Ra}}^{\frac 13}(\log {\rm{Ra}})^{\gamma}$ hold only for ${\rm{Pr}} = \infty$ or for very large Prandtl number (i.e., ${\rm{Pr}}\gtrsim ({\rm{Ra}}\ln {\rm{Ra}})^{\frac 13}$). On the other hand, the upper bound ${\rm{Nu}}\lesssim {\rm{Ra}}^{\frac 12}$ holds uniformly in ${\rm{Pr}}$ and it improves the bound ${\rm{Nu}}\lesssim ({\rm{Ra}}\ln {\rm{Ra}} {\rm{Pr}}^{-1})^{\frac 12}$ for low Prandtl numbers. At the same time, this bound is the same no matter what the boundary conditions are, as long as $u^z = 0$ at $z = \{0, 1\}$. Thus the analysis is not yet able to confirm or rule out any of the scaling theories. Also, direct numerical simulations see their limitations in power of prediction when ${\rm{Ra}}$ becomes very big and often the data points reported are not sufficient to conclude a scaling behavior as argued in ^[46,47]. Nevertheless, in ^[48] the authors show that the three dimensional simulations-data fit ${\rm{Nu}}\sim {\rm{Ra}}^{0.331}$ in a convection cell with small aspect ratio, $\Pr = 1$ and ${\rm{Ra}}\in[10^{10}, 10^{15}]$.

Do boundary conditions matter in the ultimate regime? Whether the boundary conditions matter in the ultimate–turbulent–state, remains to be understood. Nevertheless the results in the chart seem to indicate that free-slip enhance vertical heat transport compared to no-slip boundary conditions. On the other hand, for the choice $\rm{\mathcal{L}_s}\sim {\rm{Ra}}^{\alpha}$ and $\Pr\geq {\rm{Ra}}^{\frac 34-\alpha}$ the bound in (4.16) reads

If we compare this bound for small slip-lengths with the upper bound for no-slip boundary conditions at relatively small Prandtl numbers (${\rm{Nu}}\lesssim {\rm{Ra}}^{\frac 12}$) we observe that the Navier-slip conditions seem to slightly inhibit heat transport.

Dimensions: can we say something more in $2d$? Again looking at the chart we can see that little difference can be appreciated between the results in two and three dimensions, especially for no-slip boundary conditions. While the analysis in three dimensions is limited by the complexity of the Navier-Stokes equations, the two dimensional dynamics are accessible and might reveal significant mechanism in heat transport. In fact, the two dimensional Rayleigh-Bénard problem displays many of the physical and turbulent transport feature of three dimensional convection and the structure and regularity of the flow allows finer estimates. In fact, as seen in ^[12,45] balances involving the vorticity (and its gradient) can be supplemented to the energy balance to discover important relations. The next challenge in mathematical analysis would be to improve the upper bounds in ^[39], as the arguments in this paper are blind towards dimensions. In particular, in two dimensions, can one extend the regime of validity of the bound ${\rm{Nu}}\lesssim {\rm{Ra}}^{\frac 13}$ at finite Prandtl number? or can one prove a bound of the type ${\rm{Nu}}\lesssim {\rm{Ra}}^{\beta}$ with $\beta < \frac 13$ in some (possibly low) Prandtl number regime?

Are these results optimal? Let us notice that all rigorous results are in the form of upper bounds, which therefore give us indications about the "worst case" scenario. It remains to be seen whether these estimates are sharp. This is very challenging from a mathematical standpoint and, at the current time, no lower bounds for the Nusselt number have been derived. Constructing exact solutions that saturate the bounds is probably impossible, given the complexity of the system, even at infinite Prandtl number. Nevertheless, in this direction, we want to mention the works ^[49,50] in which Tobasco and Doering consider a passive tracer advected by an incompressible velocity field $\bf{u}$, satisfying $\langle|\nabla \bf{u}|^2 \rangle\leq \rm{Pe}^2$, where $\rm{Pe}$ is the Peclet number. By adopting optimal designing techniques in energy-driven pattern formation in material science, the authors succeeded in designing an incompressible flow satisfying the enstrophy constraint such that ${\rm{Nu}}\gtrsim \rm{Pe}^{\frac 23}(\log\rm{Pe})^{-\frac 43}$. Together with the upper bound ${\rm{Nu}}\lesssim \rm{Pe}^{\frac 23}$ ^[51], this lower bound demonstrate ${\rm{Nu}}\sim \rm{Pe}^{\frac 23}$ possibly up to a logarithm. In connection with the Rayleigh-Bénard convection problem, this result shows the existence of a flow which is not buoyancy driven but nevertheless realizes ${\rm{Nu}}\sim {\rm{Ra}}^{\frac 12}$.

What are the bounds in case of rough boundaries? Another important problem to study is the role of rough boundaries in the bounds. Numerical studies ^[52] indicate that roughness can enhance heat transport. To the author's knowledge, the only rigorous result available is given by Goluskin and Doering in 2016. In their paper ^[53], they assume no-slip boundary conditions on a rough vertical walls. More precisely, the (vertical) boundaries are continuous (and piecewise differentiable functions) of the horizontal coordinates (that is $z = h^T(x, y)$ and $z = h^B(x, y)$) with finite $L^2$ gradients. By a (not standard) application of the background field method, the authors show

The scaling exponent $\frac 12$ was observed in the experiments reported in ^[54] where all boundaries (including sidewalls) were rough and, most recently, in numerical simulations using sinusoidally rough upper and lower surfaces in two dimensions ^[55]. Nevertheless the bound derived by Goluskin and Doering depends on a constant of the $L^2$-norm of the gradient of $h$. It is expected that this condition can be relaxed, allowing much weaker regularity. As remarked in the introduction, it is not clear whether the no-slip boundary conditions are the most suitable to describe the physical problem. Citing ^[56] "For many years it has been observed that there is no compelling argument to justify the standard "no-slip" boundary condition of textbook continuum hydrodynamics, which states that fluid at a solid surface has no relative velocity to it ^[57]. However, this assumption successfully describes much everyday experience. Indeed, Feynman ^[57] noted in his Lectures that the no-slip condition explains why large particles are easy to remove by blowing past a surface, but small particles are not. Similarly, readers who wash dishes have noticed that it is difficult to remove all the soap just by running water — a dishcloth is needed for effective cleaning. Why?". It is then reasonable to consider surfaces where some slip occurs but there is still some stress on the fluid ^[12,58]. This suggests to consider (in two dimensions) the Navier-slip conditions (1.4), where $\alpha = \alpha(\bf{x})$ is a given positive twice continuously differentiable function defined on $\partial \Omega$, which indicates the roughness of the surface. Future research with the aim of extending the results in ^[45] to this rough-setting is planned.

Acknowledgments

The author thanks Charlie Doering for the uncountable explanations and discussions, for sharing his vision and passion on this fundamental problem. The author is grateful to Philippe Vieweg for providing the pictures in Figure 1, to Steffen Pottel for the many suggestions regarding the manuscript and to the anonymous reviewer for useful feedback. This work was partially supported by the DFG-grant GrK2583 and TRR181.

Conflict of interest

The author declares no conflict of interest.