Steady Rayleigh–Bénard convection between stress-free boundaries

Natural convection is the buoyancy-driven flow resulting from unstable density variations, typically due to thermal or compositional inhomogeneities, in the presence of a gravitational field. It remains the focus of experimental, computational, and theoretical research worldwide, in large part because buoyancy-driven flows are central to engineering heat transport, atmosphere and ocean dynamics, climate science, geodynamics, and stellar physics. Rayleigh–Bénard convection, in which a layer of fluid is confined between isothermal horizontal boundaries with the higher temperature on the underside (Rayleigh, 1916) is studied extensively as a relatively simple system displaying the essential phenomena. Beyond the importance of buoyancy-driven flow in applications, Rayleigh’s model has served for more than a century as a primary paradigm of nonlinear physics (Malkus & Veronis, 1958), complex dynamics (Lorenz, 1963), pattern formation (Newell & Whitehead, 1969), and turbulence (Kadanoff, 2001).

A central feature of Rayleigh–Bénard convection is the Nusselt number Nu, the factor by which convection enhances heat transport relative to conduction alone. A fundamental challenge for the field is to understand how Nu depends on the dimensionless control parameters: the Rayleigh number Ra, which is proportional to the imposed temperature difference across the layer, the fluid’s Prandtl number $\Pran$, and geometric parameters like the domain’s width-to-height aspect ratio $\Gamma$. Rayleigh (1916) studied the bifurcation from the static conduction state (where $\mbox{{Nu}}=1$) to convection (where $\mbox{{Nu}}>1$) when Ra exceeds a $\Pran$-independent finite value. In the strongly nonlinear large-Ra regime relevant to many applications, convective turbulence is characterized by chaotic plumes that emerge from thin thermal boundary layers and stir a statistically well-mixed bulk. Power-law behavior where Nu scales like $\Pran^{\beta}\mbox{{Ra}}^{\gamma}$ is often presumed for heat transport in the turbulent regime, but heuristic theories—i.e., physical arguments relying on uncontrolled approximations—yield various predictions for the scaling exponents. Rigorous upper bounds on Nu derived from the equations of motion place restrictions on possible asymptotic exponents but do not imply unique values. Meanwhile, direct numerical simulations (DNS) and laboratory experiments designed to respect the approximations employed in Rayleigh’s model have produced extensive data on Nu over wide ranges of Ra, $\Pran$, and $\Gamma$. Even so, consensus regarding the asymptotic large-Ra behavior of Nu remains to be achieved (Chillà & Schumacher, 2012; Doering, 2020).

In addition to the turbulent convection generally observed at large Ra, there are much simpler steady solutions to the equations of motion such as the pair of steady counter-rotating rolls shown in figure 1. Steady coherent flows are not typically seen in large-Ra simulations or experiments because they are dynamically unstable. Nonetheless, they are part of the global attractor for the infinite-dimensional dynamical system defined by Rayleigh’s model, and recent results suggest that steady rolls may be one of the key coherent states comprising the ‘backbone’ of turbulent convection. In the case of no-slip top and bottom boundaries Waleffe et al. (2015) and Sondak et al. (2015) found that, over the range of Rayleigh numbers they explored, two-dimensional (2D) steady rolls display Nu values very close to those of three-dimensional (3D) convective turbulence, provided that the horizontal period of the rolls is tuned to maximize Nu at each value of Ra.

Here we report computations of steady 2D convective rolls in the case of stress-free top and bottom boundaries. We have carried out computations using spectral methods over eight orders of magnitude in Ra, four orders of magnitude in $\Pran$, and more than an order of magnitude in the aspect ratio $\Gamma$, defined as the width-to-height ratio of a pair of rolls. As in the no-slip case, our steady states share many features with time-dependent simulations between stress-free boundaries (Paul et al., 2012; Wang et al., 2020). Moreover, the results verify predictions about the $Ra\rightarrow\infty$ limit made by Chini & Cox (2009) who extended an approach initiated by Robinson (1967) to construct matched asymptotic approximations of steady rolls between stress-free boundaries. In particular, our computations agree quantitatively with the asymptotic prediction that $\mbox{{Nu}}=\mathcal{O}(\mbox{{Ra}}^{1/3})$ uniformly in $\Pran$ with a $\Gamma$-dependent prefactor that assumes its maximum value at $\Gamma\approx 1.9$, and with corresponding asymptotic predictions about the Reynolds number that are derived in the Appendix. The rest of this paper is organized as follows. The equations governing Rayleigh–Bénard convection and our numerical scheme for computing steady solutions are outlined in § 2. The computational results are presented in § 3, followed by further discussion in § 4.

The Boussinesq approximation to the Navier–Stokes equations used by Rayleigh (1916) to model convection in a 2D fluid layer are, in dimensionless variables,

where $\mathbf{u}=u\mathbf{\hat{x}}+w\mathbf{\hat{z}}$ is the velocity, $p$ is the pressure, and $T$ is the temperature. The system has been nondimensionalized using the layer thickness $h$, the thermal diffusion time $h^{2}/\kappa$, where $\kappa$ is the thermal diffusivity, and the temperature drop $\Delta$ from the bottom boundary to the top one.

The dimensionless spatial domain is $(x,z)\in[0,\Gamma]\times[0,1]$, and all dependent variables are taken to be $\Gamma$-periodic in $x$. At the lower ($z=0$) and upper ($z=1$) boundaries, the temperature satisfies isothermal conditions while the velocity field satisfies no-penetration and stress-free boundary conditions:

The three dimensionless parameters of the problem are the aspect ratio $\Gamma$, the Prandtl number $\Pran=\nu/\kappa$, where $\nu$ is the kinematic viscosity, and the Rayleigh number $Ra=g\alpha\Delta h^{3}/\nu\kappa$, where $-g\mathbf{\hat{z}}$ is the gravitational acceleration vector and $\alpha$ is the thermal expansion coefficient. A single pair of the steady rolls computed here fits in the domain, meaning the aspect ratio of the pair is $\Gamma$ while that of each individual roll is $\Gamma/2$.

The static conduction state, for which $\mathbf{u}=\mathbf{0}$ and $T=1-z$, solves 1 and 2 at all parameter values. Rayleigh (1916) showed that rolls vertically spanning the layer with aspect ratio $\Gamma$ bifurcate supercritically from the conduction state as Ra increases past

where $k=2\pi/\Gamma$ is the wavenumber of the fundamental period of the domain. The conduction state is absolutely stable if $\mbox{{Ra}}<\mbox{{Ra}}_{c}(k)$ for all $k$ admitted by the domain; see, e.g., Goluskin (2015).

The Nusselt number is defined as the ratio of total mean heat flux in the vertical direction to the flux from conduction alone:

where $w$ and $T$ are dimensionless solutions of 1 and $\langle\cdot\rangle$ indicates an average over space and infinite time. (For steady states the time average is not needed.) The governing equations imply the equivalent expressions

the latter of which self-evidently ensures $\mbox{{Nu}}>1$ for all sustained convection. Another emergent measure of the intensity of convection is the bulk Reynolds number defined using the dimensional root-mean-squared velocity $U_{rms}$, which in terms of dimensional quantities is $\Rey={U_{rms}h}/{\nu}$. We choose our reference frame such that $\left\langle u\right\rangle=0$, so in dimensionless terms

We compute steady ($\partial_{t}=0$) solutions of 1 using a vorticity–stream function formulation,

where the stream function $\psi$ is defined by $\mathbf{u}=\mathbf{\hat{x}}\partial_{z}\psi-\mathbf{\hat{z}}\partial_{x}\psi$, the (negative) scalar vorticity is $\omega=\partial_{x}w-\partial_{z}u=-\nabla^{2}\psi$, and $\theta$ is the deviation of the temperature field $T$ from the conduction profile $1-z$. The boundary conditions used in our computations are that $\psi$, $\nabla^{2}\psi$, and $\theta$ vanish on both boundaries. The latter two conditions follow from the stress-free and fixed-temperature conditions, respectively. Impenetrability of the boundaries implies that $\psi$ is constant on each boundary, and choosing the reference frame where $\left\langle u\right\rangle=0$ requires these constants to be identical. Their value can be fixed to zero since translating $\psi$ by a constant does not affect the dynamics.

We solve the time-independent equations 7 numerically using a Newton–GMRES (generalised minimal residual) iterative scheme. Starting with an initial iterate $(\omega^{0},\psi^{0},\theta^{0})$ that does not exactly solve 7, each iteration of the Newton’s method applies a correction until the resulting iterates have converged to a solution of 7. Following Wen et al. (2015b) and Wen & Chini (2018), the linear partial differential equations for the corrections are

where the superscript $i$ denotes the $i^{th}$ Newton iterate, the corrections are defined as

and vanish on the boundaries, and

are the residuals of the nonlinear steady equations 7. We simplify the implementation by setting $F^{\psi}_{res}=0$, in which case $\triangle^{\psi}$ can be obtained by solving $\nabla^{2}\triangle^{\psi}=-\triangle^{\omega}$ for a given $\triangle^{\omega}$. After this simplification, the pair 8a and 8c can be solved simultaneously for $\triangle^{\omega}$ and $\triangle^{\theta}$.

For each iteration of Newton’s method, we solve 8a and 8c iteratively using the GMRES method (Trefethen & Bau III, 1997). The spatial discretization is spectral, using a Fourier series in $x$ and a Chebyshev collocation method in $z$ (Trefethen, 2000). The $\nabla^{2}$ operator is used as a preconditioner to accelerate convergence of the GMRES iterations. The roll states of interest have centro-reflection symmetries (cf. figure 1),

which allow the full fields to be recovered from their values on one quarter of the domain, so we encode these symmetries to reduce the number of unknowns. The GMRES iterations are stopped once the $L^{2}$-norm of the relative residual of (8$a,c$) is less than $10^{-2}$, and the Newton iterations are stopped once the $L^{2}$-norm of the relative residual of (7$a,c$) is less than $10^{-10}$. For Ra not far above the critical value $\mbox{{Ra}}_{c}(k)$, convergence to rolls of period $\Gamma=2\pi/k$ is accomplished by choosing the initial iterate with $\omega^{0}=-\sqrt{\mbox{{Ra}}\Pran}\sin(\pi z)\sin(kx)$ and $\theta^{0}=-0.1[\sin(2\pi z)+\sin(\pi z)\cos(kx)]$. For each $Pr$, results from smaller Ra (or $\Gamma$) are used as the initial iterate for larger Ra (or $\Gamma$).

We computed steady rolls over a wide range of Ra starting just above the value $\mbox{{Ra}}_{c}(k)$ at which the rolls bifurcate from the conduction state and ranging up to $10^{9}$ or higher depending on the other parameters. Computations were carried out for $\Pran=10^{-2},10^{-1},1,10,10^{2}$ and a range of values of $\Gamma$ such that the fundamental wavenumber $k=2\pi/\Gamma$ lies in $1/2\leq k\leq 10$. Data for all the $\Gamma=2$ cases are tabulated in the Supplementary Material.

Our computations reach sufficiently large Ra to show clear asymptotic scalings of bulk quantities:

Both of these scalings are predicted by the asymptotic analysis of Chini & Cox (2009), although only the Nu scaling was stated explicitly there. Chini & Cox gave an asymptotic prediction for the prefactor $c_{n}(k)$ but not for $c_{r}(k)$. Using their asymptotic approximations for the stream function and vorticity within each convection roll, we derived an expression for $c_{r}(k)$ in terms of $c_{n}(k)$ that is presented in the Appendix.

Figure 2 shows the Ra-dependence of the compensated quantities $\mbox{{Nu}}/\mbox{{Ra}}^{1/3}$ and $\Rey\Pran/\mbox{{Ra}}^{2/3}$ for rolls of aspect ratio $\Gamma=2$ ($k=\pi$) at various $\Pran$. Rolls of this aspect ratio bifurcate from the conduction state at the Rayleigh number $\mbox{{Ra}}_{c}(\pi)=8\pi^{4}\approx 779$. It is clear from figure 2 that both Nu and the Péclet number $\Rey\Pran$ become independent of $\Pran$ as $\mbox{{Ra}}\to\infty$, as predicted by the asymptotics of Chini & Cox (2009), and also as Ra decreases towards the $\Pran$-independent value $\mbox{{Ra}}_{c}$. Convergence to the large-Ra asymptotic scaling is slower when $\Pran$ is larger, at least over the four decades of $\Pran$ considered here. Numerical values of Nu and $\Rey$ at large Ra suggest scaling prefactors that are indistinguishable from the values $c_{n}(\pi)\approx 0.2723$ and $c_{r}(\pi)\approx 0.0978$ predicted by asymptotic analysis.

Nusselt and Reynolds numbers of steady rolls converge to the asymptotic scalings 12 over the full range $1/2\leq k\leq 10$ for which we have computed steady rolls. This is evident in figure 3 where the $k$-dependence of the compensated quantities $\mbox{{Nu}}/\mbox{{Ra}}^{1/3}$ and $\Rey\Pran/\mbox{{Ra}}^{2/3}$ is shown for various Ra in the $\Pran=1$ case. As Ra increases these quantities converge to asymptotic curves that we have called $c_{n}(k)$ and $c_{r}(k)$. It is clear from the figure that this convergence is slower when $k$ is larger, and that $\Rey$ reaches its asymptotic scaling sooner than Nu does. Since rolls with $k=\pi/\sqrt{2}$ are the first to bifurcate from the conduction state, at $\mbox{{Ra}}_{c}(\pi/\sqrt{2})=27\pi^{4}/4$, this is the $k$ that initially maximizes both Nu and $\Rey$. As $\mbox{{Ra}}\to\infty$, the $k$ values that maximize Nu and $\Rey$ approach the asymptotic values $k\approx 3.31$ ($\Gamma\approx 1.9$) and $k\approx 1.4$ ($\Gamma\approx 4.5$), respectively, where the corresponding maximal prefactors are $c_{n}\approx 0.273$ and $c_{r}\approx 0.117$.

Both Nu and the Péclet number $\Rey\Pran$ of steady rolls become nearly independent of $\Pran$ as Ra grows large. The large-Ra coalescence of data for different $\Pran$ is evident for the $k=\pi$ case in figure 2, as is the fact that $\Pran$ can have a substantial effect in the pre-asymptotic regime. To show that $\Pran$-independence at large Ra occurs over the full range $1/2\leq k\leq 10$ of our computations, figure 4 depicts the $k$-dependence of compensated Nu and $\Rey$ at large Ra for various $\Pran$. All $\mbox{{Nu}}/\mbox{{Ra}}^{1/3}$ and $\Rey\Pran/\mbox{{Ra}}^{2/3}$ values plotted in figure 4 fall close to the asymptotic predictions for $c_{n}(k)$ and $c_{r}(k)$ that, at leading order in the asymptotic small parameter $\mbox{{Ra}}^{-1/3}$, are independent of $\Pran$.

The asymptotic scaling of steady rolls at large Ra is reflected not only in the collapse of rescaled bulk quantities like $\mbox{{Nu}}/\mbox{{Ra}}^{1/3}$ and $\Rey\Pran/\mbox{{Ra}}^{2/3}$ but also in the collapse of the boundary and internal layer profiles when the appropriate spatial variable is stretched by $\mbox{{Ra}}^{1/3}$. Figure 5 shows this collapse of the temperature and vorticity profiles at the bottom boundary and at the left edge of the periodic domain for the case where $\Pran=1$ and $\Gamma=2$ ($k=\pi$). Coincidence of these scaled profiles at large Ra confirms that the thickness of both the thermal and vorticity layers scale like $\mbox{{Ra}}^{-1/3}$ on all four edges of a single convection roll, while both fields are strongly homogenized in the interior with $T\sim 1/2$ and $\omega\sim 0.522\mbox{{Ra}}^{2/3}$. Profiles at other $\Pran$ are not shown but collapse similarly. These findings confirm the deduction of a homogenized interior by Chini & Cox (2009), as well as their prediction that the core vorticity magnitude is asymptotic to $\sqrt{c_{n}(\pi)}\mbox{{Ra}}^{2/3}\approx 0.5218\mbox{{Ra}}^{2/3}$ uniformly in $\Pran$.

Another quantity of interest is the kinetic energy dissipation rate per unit mass,

where ^∗ denotes dimensional variables. The corresponding bulk viscous dissipation coefficient $C=\langle\varepsilon\rangle h/U^{3}_{rms}$ can be expressed in dimensionless variables as

Identity 5 gives $C=\Rey^{-3}\Pran^{-2}\mbox{{Ra}}(\mbox{{Nu}}-1)$, so the asymptotic scalings 12 imply

That is, $C$ depends asymptotically on Ra and $\Pran$ via the distinguished combination $\Pran\mbox{{Ra}}^{-2/3}$ that is asymptotic to $\Rey^{-1}$. This scaling of the dissipation coefficient is characteristic of flows without viscous boundary layers, such as laminar Couette or Poiseuille flow, consistent with the steady velocity fields computed here (cf. figure 1). Indeed, for stress-free steady convection, viscous dissipation is dominated by that in the homogenized core since the vorticity is of the same asymptotic magnitude in the core as in the thin vorticity layers.

The average of dissipation over time and horizontal directions, denoted $\overline{\varepsilon}(z)$, has been used to compare convection between the cases of stress-free and no-slip boundaries. In 3D simulations of the stress-free case at $\mbox{{Ra}}=5\times 10^{6}$, Petschel et al. (2013) found that the normalized profile $\overline{\varepsilon}(z)/\langle\varepsilon\rangle$ exhibits ‘dissipation layers’ near the boundaries that depend strongly on $\Pran$. In steady rolls, on the other hand, we find that $\overline{\varepsilon}(z)/\langle\varepsilon\rangle$ is independent of $\Pran$ at asymptotically large Ra, as shown in figure S1 of the Supplementary Material.

The steady rolls we have computed share many features with unsteady flows from DNS of Rayleigh–Bénard convection with isothermal stress-free boundary conditions. In recent simulations, Wang et al. (2020) found multistability between unsteady states exhibiting various numbers of roll pairs in wide 2D domains. Each of the coexisting states suggested scalings approaching the $\mbox{{Nu}}=\mathcal{O}(Ra^{1/3})$ and $\Rey=\mathcal{O}(\mbox{{Ra}}^{2/3})$ asymptotic behaviour of steady rolls. As in the steady case, the prefactors of these scalings depended on the mean aspect ratios of the unsteady rolls. The highest Nusselt numbers among Wang et al.’s data occur in five-roll-pair states in a $\Gamma=16$ domain—meaning each roll pair has $\Gamma\approx 3.2$ on average—but steady $\Gamma=3.2$ rolls have still larger Nu. At $\mbox{{Ra}}=10^{9}$ and $\Pran=10$, for example, the DNS exhibit $\mbox{{Nu}}=198.01$ and $\Rey=10135$ while steady $\Gamma=3.2$ rolls at the same parameters yield the larger values of $\mbox{{Nu}}=253.61$ and $\Rey=11333$, and comparisons at other Ra are similar (cf. Table S4 of the Supplementary Material). Figure 6 shows the Nu of these five-roll-pair DNS states along with the larger Nu of the steady rolls computed here for various $\Pran$ and $\Gamma=2$. The steady rolls also achieve larger Nu values than have been attained in other unsteady simulations with stress-free boundaries in 2D (van der Poel et al., 2014; Goluskin et al., 2014) and in 3D (Petschel et al., 2013; Pandey et al., 2014; Pandey & Verma, 2016; Pandey et al., 2016).

Comparing steady rolls in the stress-free case with those previously computed in the no-slip case, there are significant differences in their dependence on the aspect ratio $\Gamma$. With stress-free boundaries, $\mbox{{Nu}}=\mathcal{O}(Ra^{1/3})$ for each $\Gamma$ as $\mbox{{Ra}}\to\infty$, with maximal asymptotic heat transport attained by rolls of optimal aspect ratio $\Gamma\approx 1.9$. In the no-slip computations of Waleffe et al. (2015) and Sondak et al. (2015), on the other hand, the $\Gamma$ values that maximize Nu decrease towards zero proportionally to $\mbox{{Ra}}^{-0.22}$ at large Ra. (A similar phenomenon occurs in porous medium Rayleigh–Bénard convection; see Wen et al. (2015b).) The no-slip steady rolls display Nu scaling like $\mbox{{Ra}}^{0.28}$ when $\Gamma$ is fixed but scaling like $\mbox{{Ra}}^{0.31}$ when the optimal $\Gamma$ is chosen to maximize Nu at each Ra. The measured exponent 0.31 is unlikely to be exact, so it remains possible that the asymptotic scaling of optimal-$\Gamma$ steady rolls is $\mbox{{Nu}}=\mathcal{O}(\mbox{{Ra}}^{1/3})$ in the no-slip case, as in the stress-free case.

Steady rolls in the stress-free and no-slip cases differ also in their dependence on the Prandtl number. Only with stress-free boundaries do the Nusselt number Nu and Péclet number $\Rey\Pran$ apparently become uniform in $\Pran$ as $\mbox{{Ra}}\to\infty$. To see how this $Pr$-independence emerges in the stress-free scenario, first note that area-integrated work by the buoyancy forces must balance area-integrated viscous dissipation—i.e., $\mbox{{Ra}}\left\langle wT\right\rangle=\left\langle|\bnabla u|^{2}+|\bnabla w|^{2}\right\rangle$ in the dimensionless formulation of 4 and 5. The former integral is dominated by plumes since the roll’s core is isothermal, so it scales proportionally to the dimensional quantity $\alpha\Delta g\delta hU_{rms}$, where $\delta$ is the dimensional plume thickness. The latter integral is dominated by the core since the vorticity is $\mathcal{O}(U_{rms}/h)$ everywhere in the stress-free case, so this integral scales proportionally to the dimensional quantity $\nu(U_{rms}/h)^{2}h^{2}$. Balancing advection with diffusion of temperature anomalies in the thermal boundary layers, which also have thickness $\delta$, requires that $U_{rms}$ scale in proportion to $\kappa h/\delta^{2}$. Combining this scaling relationship with that from the integral balance gives the dimensionless thermal boundary layer thickness as $\delta/h=\mathcal{O}(\mbox{{Ra}}^{-1/3})$—and so $\mbox{{Nu}}=\mathcal{O}(Ra^{-1/3})$—uniformly in $\Pran$. These relationships also imply that $U_{rms}$ is proportional to $(\kappa/h)\mbox{{Ra}}^{2/3}$, and so the Péclet number $\Rey\Pran=U_{rms}h/\kappa$ scales as $\mbox{{Ra}}^{2/3}$ uniformly in $\Pran$. Ultimately, it is the passivity of the vorticity boundary layers that results in the $Pr$-independence of these emergent bulk quantities. The vorticity layers not only make no contribution to the total dissipation at leading order, they also have no leading-order effect on the stream function that is responsible for the convective flux $\langle wT\rangle$.

The $\Rey=\mathcal{O}(\mbox{{Ra}}^{2/3})$ scaling found at large Ra for steady rolls and approximately evidenced in the DNS of Wang et al. (2020) means that buoyancy forces can sustain substantially faster-than-free-fall velocities. Indeed, if flow speeds were limited by the maximum buoyancy acceleration acting across the layer height then dimensional characteristic velocities could not be of larger order than $\sqrt{g\alpha\Delta h}$, and $\Rey$ could not be larger than $\mathcal{O}(\mbox{{Ra}}^{1/2})$. Such $\Rey$ may be expected if the bulk flow were dominated by effectively independent rising and falling plumes. Significantly higher speeds apparently persist within coherent convection rolls, whether steady or unsteady.

Although steady rolls cannot give heat transport larger than $\mbox{{Nu}}=\mathcal{O}(\mbox{{Ra}}^{1/3})$ as $\mbox{{Ra}}\to\infty$ in the stress-free case (Chini & Cox, 2009), it is an open question whether larger Nu can result from time-dependent flows or other types of steady states in either 2D or 3D. Rigorous upper bounds on Nu derived from the governing equations—bounds depending on Ra that apply to all flows regardless of whether they are steady or unsteady and stable or unstable—do not rule out Nu growing faster than $\mbox{{Ra}}^{1/3}$. Specifically, for the 2D stress-free case Whitehead & Doering (2011) proved analytically that $\mbox{{Nu}}\leq 0.289\,\mbox{{Ra}}^{5/12}$ uniformly in both $\Pran$ and domain aspect ratio. Wen et al. (2015a) improved the prefactor of this bound by solving the relevant variational problem numerically, computing bounds up to large finite Ra with a prefactor depending weakly on $\Gamma$. The numerical upper bound they computed for $\Gamma=2\sqrt{2}$ is shown in figure 6; its scaling at large Ra is $\mbox{{Nu}}\leq 0.106\,\mbox{{Ra}}^{5/12}$. For 3D flows in the stress-free case only the larger upper bound $\mbox{{Nu}}\leq\mathcal{O}(\mbox{{Ra}}^{1/2})$ has been proved (Doering & Constantin, 1996). It remains to be seen whether upper bounds smaller than $\mathcal{O}(\mbox{{Ra}}^{5/12})$ or $\mathcal{O}(\mbox{{Ra}}^{1/2})$ can be proved for 2D or 3D flows, respectively, and whether there exists any sequence of solutions for which Nu grows faster than $\mathcal{O}{(\mbox{{Ra}}^{1/3})}$. In view of available evidence it is possible that, at each Ra and $\Pran$, the steady 2D roll with the largest value of Nu—i.e., with Nu maximized over $\Gamma$—transports more heat than any other 2D or 3D solution. We are aware of no counterexamples to this possibility, either in the stress-free case studied here or in the no-slip case.

We thank L.M. Smith, D. Sondak and F. Waleffe for helpful discussions. This research was supported in part by US National Science Foundation awards DMS-1515161 and DMS-1813003, Canadian NSERC Discovery Grants Program awards RGPIN-2018-04263, RGPAS-2018-522657 and DGECR-2018-00371, and computational resources and services provided by Advanced Research Computing at the University of Michigan.

The authors report no conflict of interest.