Steady but generally unstable solutions of the 2D Boussinesq equations are obtained for no-slip boundary conditions and Prandtl number 7. The primary solution that bifurcates from the conduction state at Rayleigh number *Ra* ≈ 1708 has been calculated up to *Ra* ≈ 5.10^{6} and its Nusselt number is *Nu* ∼ 0.143 *Ra*^{0.28} with a delicate spiral structure in the temperature field. Another solution that maximizes *Nu* over the horizontal wavenumber has been calculated up to *Ra* = 10^{9} and scales as *Nu* ∼ 0.115 *Ra*^{0.31} for 10^{7} < *Ra* ≤ 10^{9}, quite similar to 3D turbulent data that show *Nu* ∼ 0.105 *Ra*^{0.31} in that range. The optimum solution is a simple yet multi-scale coherent solution whose horizontal wavenumber scales as 0.133 *Ra*^{0.217}. That solution is unstable to larger scale perturbations and in particular to mean shear flows, yet it appears to be relevant as a backbone for turbulent solutions, possibly setting the scale, strength, and spacing of elemental plumes.

Rayleigh-Bénard convection is the buoyancy-driven motion of a fluid contained between two horizontal plates and heated from below. It is a paradigmatic problem with a rich array of nonlinear physics from deterministic chaos^{1} and period doubling^{2} to pattern formation^{3} and turbulence.^{4} A fundamental issue is to determine how the heat flux $H$ scales with the temperature difference Δ*T* between the bottom and top plates. There is a critical Δ*T _{c}* such that heat transport is by conduction with $H\u223c\Delta T$ for Δ

*T*< Δ

*T*but macroscopic fluid motion develops for Δ

_{c}*T*> Δ

*T*and enhances heat transport. Bifurcations take place as Δ

_{c}*T*is increased, leading to turbulent flow and $H\u223c ( \Delta T ) 4 / 3 $ and perhaps even as high as $H\u223c ( \Delta T ) 3 / 2 $. Understanding these bifurcations and the actual asymptotic scaling of $H$ have been the subjects of many experimental, theoretical, and numerical studies focusing on

*turbulent*convection. Here, we consider

*coherent*convection—steady flows that may be the permanent form of the coherent structures (

*“plumes,” “thermals,”*and large scale

*“winds”*) permeating turbulent convection.

In Rayleigh-Bénard convection, the fluid is contained between two infinite horizontal plates and the Boussinesq approximation is made so the governing equations are

for the velocity ** v** and the temperature

*T*. Incompressibility

**∇**⋅

**= 0 is maintained by the kinematic pressure**

*v**p*and $\u2212g y \u02c6 $ is the constant acceleration of gravity. The fluid density is

*ρ*≈

*ρ*

_{0}(1 −

*α*

_{V}

*T*), with

*α*

_{V}≥ 0 the volumetric thermal expansion coefficient and

*ρ*

_{0}a constant reference density,

*ν*> 0 is the kinematic viscosity, and

*κ*> 0 is the thermal diffusivity. The bottom plate temperature is fixed at Δ

*T*/2 and the upper plate is at −Δ

*T*/2. The conduction state is the solution

**= 0,**

*v**T*= −

*y*Δ

*T*/

*H*to (1) and (2), with −

*h*≤

*y*≤

*h*=

*H*/2. This is a steady solution for all values of the parameters, but Rayleigh showed that it is linearly unstable when the

*Rayleigh*number

*Ra*=

*gα*

_{V}Δ

*T*

*H*

^{3}/(

*νκ*) is larger than a critical value

*Ra*, independent of the

_{c}*Prandtl*number

*Pr*=

*ν*/

*κ*.

We consider no-slip boundary conditions for which *Ra _{c}* ≈ 1708 for horizontal wavenumber

^{5}

*α*= 3.116/

_{c}*H*. Convection develops and enhances heat transport for

*Ra*>

*Ra*. For statistically steady solutions, integration of (2) yields the heat flux,

_{c}^{6}

where $ H 0 =\kappa \Delta T/H$ is the heat flux in the absence of fluid motion and $v= y \u02c6 \u22c5v$ is the vertical velocity with $ \u2009 \u2217 \u2009 \xaf $ and 〈∗〉 denoting horizontal and domain averages, respectively.

A classic scaling argument is that heat transport is determined by marginal stability of the thermal boundary layers and of the mean temperature gradient in the interior.^{6} A simple version of this argument assumes an isothermal *T* = 0 interior with boundary layers of thickness *δ* such that^{7} *Ra*^{(δ)} ≡ *gα*_{V}(Δ*T*/2) *δ*^{3}/(*νκ*) ≈ 1708/16, yielding heat flux

independent of *Pr*, where the *Nusselt* number $Nu=H/ H 0 $. This scaling also follows from assuming^{8} that the heat flux becomes independent of the height *H* as *Ra* → ∞.

Another simple argument imagines a fluid at temperature *T* = 0, with cold plumes at temperature −Δ*T*/2 free-falling a distance *H* = 2*h* from top to bottom plates (and hot Δ*T*/2 plumes “free-rising” from bottom to top) at average speed $V= g \u2032 h $, where *g*′ = *gα*_{V}Δ*T*/2 is the reduced gravity. This yields the heat flux

This *inertial* scaling also follows from the standard turbulence assumption^{8} that heat transport becomes independent of *ν* and *κ* in the limit *Vh*/*ν* → ∞, *Vh*/*κ* → ∞.

Kraichnan’s mixing length theory^{9} yields various scalings in distinct regions of the (*Ra*, *Pr*) parameter space. He obtains *Nu* ∼ *Ra*^{1/3}, for fixed *Pr*, but predicts a transition to *Nu* ∼ *Ra*^{1/2}(ln*Ra*)^{−3/2} for very large *Ra* ≳ 10^{12} when the shear boundary layers would be turbulent. The comprehensive scaling theory of Grossmann and Lohse^{10} incorporates classic Kolmogorov scaling of velocity and temperature dissipation rates in the bulk together with thermal and viscous boundary layers. Their theory fits many experimental data sets.^{11} They predict a transition to an “ultimate” regime similar to Kraichnan’s for *Ra* ≳ 10^{14} and argue that their log correction would yield an effective scaling^{12} of about *Nu* ∼ *Ra*^{0.38} in the range 10^{12} < *Ra* < 10^{15}. Some experimental results show transition to an ultimate regime^{13–16} while others do not.^{17,18} However, all the data are well-fitted by *Nu* ∼ 0.105 *Ra*^{0.31} for *Ra* < 10^{11} and that is the regime for which we report on unstable coherent states, in particular optimum transport solutions with *Nu* ∼ 0.115 *Ra*^{0.31}.

Rigorous upper bounds on heat transport for no-slip boundary conditions^{19–23} yield *Nu* − 1 ≲ 0.026 *Ra*^{1/2} as *Ra* → ∞, for any *Pr*, showing that free-fall scaling (5) cannot hold for large *Pr*. Rigorous bounds for free-slip^{24} yield *Nu* − 1 ≲ 0.2295 *Ra*^{5/12} that conflicts with (5) for any *Pr* (the bound has been improved to *Nu* − 1 ≲ 0.106 *Ra*^{5/12}).^{29}

We consider 2D flow with $v=u ( x , y , t ) x \u02c6 +v ( x , y , t ) y \u02c6 $ and use *h* = *H*/2, $V= g \u2032 h $, *τ* = *h*/*V*, and Δ*T*/2 as our characteristic length, velocity, time, and temperature scales, respectively. Eliminating *p* by taking the $ y \u02c6 $ component of the curl of (1) yields

where *ν* and *κ* are now non-dimensionalized by $Vh= g \u2032 h 3 $ and relate to the Rayleigh and Prandtl numbers as

All results in this paper are for *Pr* = 7 (water). The horizontal velocity *u* is obtained from *v* and ∂_{x}*u* + ∂_{y}*v* = 0 except for its horizontal average $ u \xaf ( y , t ) $ that is determined by (8), where the overbar denotes a horizontal average. Equations (6)–(8) are considered with no-slip boundary conditions,

together with periodicity of period *L* = 2*π*/*α* in the *x* direction.

We look for *steady* solutions, ∂_{t} = 0, that bifurcate from the conduction state at *Ra* ≈ 1708 for wavenumber *α* ≈ 3.116/2 corresponding to aspect ratio *L*/*H* ≈ 2.016. Those convective solutions obey mirror symmetry,

as well as shift-reflect symmetry,

Equations (6) and (7) are discretized using a Fourier expansion in *x* and Chebyshev integration^{25} in *y* with time-marching to steady states. Two distinct codes have been written; code 1 uses a 3rd order time-accurate scheme^{26} and code 2 uses a semi-implicit forward-backward Euler time discretization together with Chebyshev tau of the 2nd kind^{27} for the *v* equation. Code 2 typically uses inconsistent time integration, with smaller time steps for smaller wavenumbers and for the *v* equation than for the *T* equation, to speed up or enable convergence to steady state. Both codes are dealiased in *x* and *y* using the 2/3 rule. The results of both codes overlap or connect very well as seen in Figs. 1–3, discussed below. The results also match a 3rd code based on finite differences and a Newton-Krylov iteration to steady state (Sondak, Smith, Waleffe 2015, not shown here). Mirror symmetry (11) is imposed and eliminates mean flow (8). Shift-reflect symmetry (12) is not imposed explicitly but is satisfied by the steady solutions presented here.

We show two branches of nonlinear steady solutions that bifurcate from the conduction state at *Ra* ≈ 1708, *α* ≈ 1.558. The *primary branch* (PB) has fixed horizontal wavenumber *α* (rounded here to *α* = *π*/2⇔*L*/*H* = 2). The *optimum branch* (OB) adjusts *α* to maximize heat flux *Nu*. The (*Ra*, *Nu*) curves for both solutions are shown in Figures 1 and 2. The steady primary solution is stable up to *Ra* ≈ 53 000 where it spawns a time-periodic solution. The maximum and minimum *Nu* achieved by that periodic solution are shown in Fig. 1; that time periodic solution was computed with the time-accurate code 1. The unstable steady state was continued with code 2 up to *Ra* = 5. 10^{6} (symbols ∘ in Fig. 2). Mirror symmetry (11) suffices to stabilize the optimum solution in its fundamental periodic domain. That optimum solution was calculated up to *Ra* = 1.4 10^{6} with code 1 and *Ra* = 10^{9} with code 2 (symbols * in Fig. 2). The optimum solution at *Ra* = 10^{9} is well-resolved with 200 Chebyshev polynomials in *y* and 200 Fourier modes in *x*, after dealiasing.

The heat flux for the unstable primary solution scales as *Nu* − 1 ≈ 0.143 *Ra*^{0.28} according to a least square fit in 5.10^{5} ≤ *Ra* ≤ 5.10^{6}. That solution develops a delicate spiral structure in the temperature field but not in the velocity (Fig. 4). The winding of these temperature spirals continuously increases with Rayleigh number. Convergence of our algorithm to that unstable steady solution becomes quite slow for increasing *Ra*, apparently because of the center region of these spiral structures where *v*, *T* ≈ 0.

The optimum solution increases wavenumber *α* with *Ra* as *α* ≈ 0.133 *Ra*^{0.217} (Fig. 3) to achieve an optimum heat flux *Nu* − 1 ≈ 0.115 *Ra*^{0.31} (Fig. 2), from a least square fit in 10^{7} < *Ra* ≤ 10^{9}. This optimum heat flux scaling is quite similar to the heat flux observed in 3D turbulent convection experiments. Indeed, the optimum branch in Fig. 2 is only slightly above the cryogenic helium gas data fit^{18} *Nu* ≈ 0.088 *Ra*^{0.32} and the pressurized *SF*_{6} gas data fit^{16} *Nu* ≈ 0.105 *Ra*^{0.312} (both in a cylinder of aspect ratio 1/2 and both for *Ra* < 10^{11}). The aspect ratio 4 data^{18} (Table 1) are well-fitted by *Nu* − 1 ≈ 0.102 *Ra*^{0.31} in 10^{8} ≤ *Ra* ≤ 10^{10} and lie just below the optimum transport data with *Nu* − 1 ≈ 0.115 *Ra*^{0.31}. This close agreement is remarkable given the differences in dimension, 2D optimum vs. 3D data, and Prandtl number, *Pr* = 7 for the 2D optimum vs. *Pr* ≈ 0.7 in the helium gas experiments.

The optimizing wavenumber *α* = *α ^{opt}*(

*Ra*) (Fig. 3) shows an undulation between 1708 <

*Ra*≲ 10

^{6}that is linked to a temperature spiral structure in the optimum transport solution as well (Fig. 5). A temperature updraft of hot fluid (and downdraft of cold fluid) develops between

*Ra*= 1708 and

*Ra*≈ 10

^{4}, but a spiral structure with slight downdraft of warm fluid (and updraft of cool) develops in 10

^{4}≲

*Ra*≲ 10

^{5}, corresponding to the bump in the curve

*α*=

*α*(

^{opt}*Ra*) in Fig. 3. The warm spiral begins winding back up at

*Ra*≈ 10

^{5}. The spiral does not appear to ever wind back down for higher

*Ra*, developing instead an increasingly jagged structure (Fig. 5, right), and

*α*approaches the power law scaling

^{opt}*α*≃ 0.133

^{opt}*Ra*

^{0.217}. The structure of optimum

*Nu*solutions depends on Prandtl number

*Pr*. This is investigated in forthcoming work which also shows that the locally optimum steady solution discussed here is in fact the global optimum over

*α*.

The close agreement between the optimum transport 2D steady state solutions and 3D turbulent data in Fig. 2 is intriguing. Although this agreement could be fortuitous, it strongly suggests that a single unstable steady solution may capture key statistical features of fully developed turbulent flows, such as the net heat flux and the mean temperature profile as well as the strength and scale of elemental plumes, as in shear flows.^{25,28} A search for such maximum momentum transport solutions in shear flows was initiated over 10 yr ago but not completed because of the higher computational complexity of those 3D solutions.^{30}

Our results are connected with Malkus’ theory of turbulent convection^{6} and subsequent work on upper bounds.^{19–24} Two key ingredients of Malkus’ theory are maximum heat transport and marginal stability of the mean temperature profile and the smallest scales of motion. Both ingredients are included in our calculations that maximize heat transport over horizontal wavenumber *α* and track steady state solutions that bifurcate from the marginal stability critical point at *Ra* ≈ 1708, *α* ≈ 1.558. We conjecture that our 2D results are in fact the 3D optimum transport solutions (for infinite or periodic horizontal directions). The upper bound results also assume 2D optimizers.

Why the optimum transport solution should capture gross 3D turbulence characteristics might be understood as a “winner-take-all” effect, where the optimum solution consumes all available potential energy so no other flow can be sustained. The optimum solution appears to be stable when mirror symmetry (11) is imposed and length scales are restricted to be less or equal to the optimum wavelength. The optimum solution is unstable to larger scale perturbations, in particular to subharmonics where plumes merge and form bigger plumes in a cyclic or quasi-cyclic fashion. It is also unstable to mean shear flow (8) when mirror symmetry is allowed to be broken. These instabilities would be the source of the “turbulence” but the underlying unstable coherent solutions control the heat transport.

The authors would like to thank Charles Doering (University of Michigan) and David Sondak (University of Wisconsin-Madison) for helpful comments and discussions and Detlef Lohse (University of Twente) for providing many recent references on the ultimate regime and the Grossmann-Lohse scaling theory. This research was partially supported by NSF Grant Nos. DMS-0807349 (FW, AB) and DMS-1008396 (AB, LMS).

## REFERENCES

Estimated from linear stability about the conduction state but with free “outflow” boundary conditions at the top edge of the boundary layer, that is, $\u2202yv=\u2202y3v=\u2202y(T+y)=0$.