Improved bounds on horizontal convection

For the problem of horizontal convection the Nusselt number based on entropy production is bounded from above by $C\,Ra^{1/3}$ as the horizontal convective Rayleigh number $Ra\rightarrow \infty$ for some constant $C$ (Siggers et al., J. Fluid Mech., vol. 517, 2004, pp. 55–70). We re-examine the variational arguments leading to this ‘ultimate regime’ by using the Wentzel–Kramers–Brillouin method to solve the variational problem in the $Ra\rightarrow \infty$ limit and exhibiting solutions that achieve the ultimate $Ra^{1/3}$ scaling. As expected, the optimizing flows have a boundary layer of thickness ${\sim}Ra^{-1/3}$ pressed against the non-uniformly heated surface; but the variational solutions also have rapid oscillatory variation with wavelength ${\sim}Ra^{-1/3}$ along the wall. As a result of the exact solution of the variational problem, the constant $C$ is smaller than the previous estimate by a factor of $2.5$ for no-slip and $1.6$ for no-stress boundary conditions. This modest reduction in $C$ indicates that the inequalities used by Siggers et al. (J. Fluid Mech., vol. 517, 2004, pp. 55–70) are surprisingly accurate.

883 A41-2 C. B. Rocha and others layer is heated at the bottom and cooled at the top. In that context HC is interesting because there is a restrictive bound on the rate of dissipation of kinetic energy and on net vertical buoyancy flux (Paparella & Young 2002;Scotti & White 2011;Gayen et al. 2013).
These strands of HC research are entwined because buoyancy (or heat) transport is a prime index of the strength of the HC, and also of the strength of ocean circulation. From this perspective, the strong bound on kinetic energy dissipation is of secondary importance except, perhaps, for its role in limiting the buoyancy transport of HC.
The strength of HC is measured by the Nusselt number Nu. (Notation is defined in § 2.) The problem of developing bounds on HC Nu was first addressed by Siggers, Kerswell & Balmforth (2004), SKB hereafter, using the background method of Doering & Constantin (1996). The main result of SKB is that the HC Nusselt number, based on entropy production, has the Ra → ∞ bound Nu < C Nu Ra 1/3 , (1.1) where C Nu is a constant and Ra is the HC Rayleigh number. Siggers et al. (2004) avoided detailed solution of the relevant variational problem and instead expeditiously obtained (1.1) via relatively simple inequalities, such as those developed here in appendix A. Winters & Young (2009) obtained the upper bound (1.1) using a different set of inequalities due to Howard (1972). The approach of Winters & Young (2009) results in a substantially larger constant C Nu than that of SKB. Siggers et al. (2004) refer to the scaling in (1.1) as the 'ultimate regime' of HC; see also Shishkina, Grossmann & Lohse (2016). The expectation is that at sufficiently high Ra, HC will achieve the scaling Nu ∼ Ra 1/3 , and that further increases in Ra will not change the exponent from 1/3. To better understand the underpinnings of the hypothetical ultimate regime, in § 4 we exhibit the Ra → ∞ solution of a relevant variational problem and thus find smaller values of the constant C Nu . While the exponent is unchanged from 1/3, the variational solution is of interest because it may contain clues as to the structure of the Boussinesq flows that might achieve the ultimate scaling.

Formulation of the horizontal convection problem
Consider a Boussinesq fluid with density ρ = ρ 0 (1 − g −1 b), where ρ 0 is a constant reference density and b is the 'buoyancy'. If the fluid is stratified by temperature variations then b = gα(T − T 0 ) where T 0 is a reference temperature and α is the thermal expansion coefficient. The Boussinesq equations of motion are The kinematic viscosity is ν, the thermal diffusivity is κ and ∆ = ∂ 2 x + ∂ 2 y + ∂ 2 z is the Laplacian.
We suppose the fluid occupies a rectangular domain with depth h, length x and width y ; we assume periodicity in the x and y directions. At the bottom surface (z = 0) and top surface (z = h) the boundary conditions on the velocity u = (u, v, w) are w = 0 and for the viscous boundary condition either no-slip, u = v = 0, or no-stress, u z = v z = 0. At z = 0 the buoyancy boundary condition is no-flux, κb z = 0 and at the top, z = h, the boundary condition is b = b s (x), where the top surface buoyancy b s is a prescribed function of x. As an illustrative surface buoyancy field we use where k def = 2π/ x . This choice is convenient for the analytic work in later sections. Figure 1 shows a two-dimensional solution ( y = 0) illustrating the formation of the buoyancy boundary layer adjacent to the non-uniform upper surface. The solution in figure 1 is unsteady due to vacillation of the plume.
The problem is characterized by four non-dimensional parameters: the Rayleigh and Prandtl numbers and the aspect ratios x /h and y /h. In (2.4), the coordinate x is distinguished by the assumption that the imposed surface buoyancy varies only along the x-axis. Hence the definition of Ra in (2.5) employs x . Rosevear, Gayen & Griffiths (2017) have recently considered surface buoyancy distributions with two-dimensional variation, such as b s ∝ sin 2 kx sin 2 ky. Here, however, we confine attention to the most basic version of HC in which variation of the surface buoyancy b s is only along the x-axis.
We use an overbar to denote an average over x, y and t, taken at any fixed z and to denote a total volume average over x, y, z and t. Using this notation, we recall some results from Paparella & Young (2002) that are used below.

C. B. Rocha and others
Taking u · (2.1) , we obtain the kinetic energy power integral where ε def = ν |∇u| 2 is the rate of dissipation of kinetic energy and wb is rate of conversion between potential and kinetic energy.
In § 3 the inequality in (2.9) is employed to obtain Ra → ∞ bounds on Nu. We have not been able to obtain an analytic estimate of the difference between −b(0) and the largest value allowed by the extremum principle, namely b . However, in the example shown in figure 1, with Ra = 6.4 × 10 9 , the bottom buoyancy isb(0) ≈ −0.755b and thus the right of inequality (2.9) is approximately 33 % larger than hε. In figure 2 we show results from a suite of calculations with Pr = 1 and varying Ra. These numerical results suggest, but do not prove, that as Ra → ∞ the ratiob(0)/b approaches a limiting value that is greater than −1. The solution in figure 1 is in this hypothetical regime: figure 2 shows that increasing Ra by a factor of ten does not substantially changeb(0)/b . Chiu-Webster, Hinch & Lister (2008) consider very viscous HC, that is Pr = ∞; their figure 25 showsb(0)/b extrapolated to Ra = ∞. Based on this numerical result, Chiu-Webster et al. (2008) also conclude that the bottom buoyancy remains greater than the minimum value found on the upper surface. The main result from these numerical investigations is that inequality (2.9) entails an overestimate of hε that in some cases is as large as 30 %. Young (2010) shows that if the Mach number is zero, then the small divergence of the exact non-Boussinesq velocity, U, due, for instance, to thermal expansion, can be diagnosed from within the Boussinesq approximation as ∇ · U ≈ κ b/g. Using this result, the energy source on the right of (2.8) can alternatively be written as h −zκ b = P∇ · U , where the total pressure P is well approximated by the hydrostatic background −ρ 0 gz; see also the appendix of Wang & Huang (2005). In other words, the energy source on the right of (2.8) is the conversion of internal energy into kinetic energy by 'piston work' −P dV.
3. Bounds on the horizontal-convective Nusselt number Siggers et al. (2004) discuss the difficulty of defining the Nusselt number of HC. We follow Paparella & Young (2002) and SKB by first introducing the diffusive dissipation of buoyancy variance, (3.1) Then the Nusselt number is Nu 3.2 ÷ 10 9 6.4 ÷ 10 9 6.4 ÷ 10 10 0.5 FIGURE 2. The 'instantaneous bottom buoyancy', defined as an (x, y)-average of b(x, y, 0, t). This is a suite of five no-stress solutions with the sinusoidal b s in (2.4); all solutions started with initial buoyancy −0.7b . Control parameters are Pr = 1, x /h = 4 and y = 0 (two-dimensional solutions); the Rayleigh number Ra is indicated in the legend. These computations were performed with tools developed by the Dedalus project: a spectral framework for solving partial differential equations (Burns et al. 2019, www.dedalus-project.org).
where χ diff def = κ |∇b diff | 2 is the buoyancy dissipation of the diffusive solution i.e. κ b diff = 0 with b diff satisfying the same boundary conditions as b.
An advantage in bounding χ for HC is that the elementary bound (2.9) takes care of the kinetic energy dissipation ε. Bounds on HC are simpler than those of Rayleigh-Bénard convection because in the Rayleigh-Bénard problem it is necessary to deal with ε and χ simultaneously (Doering & Constantin 1996;Kerswell 2001). In the HC problem the strategy is to first ignore the momentum equations (2.1) and obtain a bound on χ considering only the buoyancy equation (2.2) and incompressibility (2.3). The resulting bound in (3.17) below applies equally to a passive scalar and contains |∇u| 2 (but not ν). Dynamical information is then finally injected by replacing |∇u| 2 with ε/ν and using (2.9) to bound ε. This approach identifies the ε-bound in (2.9) as the only dynamical information incorporated into the χ -bound.
We begin by introducing a 'comparison function', c(x), which is any time-independent function that satisfies the same boundary conditions as b(x, t): c z (x, y, 0) = 0 and c(x, y, h) = b s (x). Forming c(2.2) , we find χ = κ ∇b · ∇c − bu · ∇c .
where the function τ is dimensionless. Without loss of generality we choose the Boussinesq reference density ρ 0 so that the horizontal average of τ is zero:τ = 0. As µh → ∞ the comparison function (3.4) forms a boundary layer adjacent to the non-uniform surface at z = h. Thus the integrands on the right-hand side of the identity (3.3) are localized in the boundary layer at z = h. This observation shows that χ and Nu can be determined solely by the statistical properties of b and u within the boundary layer. Thus, at least as far as Nu is concerned, conditions in the bulk of the flow can be ignored.
We decompose the buoyancy as so that a(x, t) has homogeneous boundary conditions: a(x, y, h, t) = 0 and a z (x, y, 0, t) = 0. Substituting the decomposition (3.7) into (3.3) we find An upper bound on χ follows from (3.8) by noting that there is an η such that au · ∇c η |∇a| 2 |∇u| 2 . (3.9) With the family of comparison functions in (3.4), η is a function of µ that is most precisely determined by solving the variational problem sequestered in § 4. Following SKB and using Young's inequality, equation (3.9) implies that where ω is any positive number. Substituting inequality (3.10) into the identity (3.8) one obtains Judiciously choosing ω = η/κ makes the middle term on the right of (3.12) equal to χ /2, so that (3.12) collapses to The inequality (3.13) is the main tool used to upper-bound Nu in (3.2). The next step is to determine the dependence of η on the boundary-layer parameter µ in the limit µ → ∞. The result from § 4 is that as µ → ∞ where the non-dimensional factor K depends on the surface-buoyancy shape function τ (x), the viscous boundary condition, but not the aspect ratio of the domain e.g. for τ = cos kx see (4.51) and (4.52). We proceed by substituting the µ → ∞ asymptotic result (3.14) into (3.13). Evaluating the first term on the right of (3.13) in the µ → ∞ limit we then have an asymptotic inequality The best upper bound on χ is obtained by finding the µ that minimizes the right-hand side of (3.15). This is The best upper bound is therefore The bound in (3.17) follows from the buoyancy equation (2.2) and the incompressibility condition (2.3) -the momentum equation has not been used. Dynamical information is supplied using inequality (2.9) to bound |∇u| 2 in (3.17) by κb /νh, resulting in After normalization by χ diff , (3.18) results in the Ra 1/3 -bound in (1.1). Using b s (x) = b cos kx as an illustration, we express (3.18) in non-dimensional variables. In this case τ 2 = 1/2 and χ diff = κkb 2 tanh kh/(2h), so that The aspect ratio, h/ x = kh/(2π), enters the bound only through the normalization factor χ diff : the constant C χ in (3.18) does not depend on aspect ratio. We have four different estimates of the constant K, leading to four different values for the prefactor C Nu in (3.19). The simplest results for C Nu follow from the asymptotic inequalities of SKB, (3.20) The analogous result from appendix A is K = 0.3458 in (A 14), so that

C. B. Rocha and others
The values of C Nu above are independent of boundary conditions. More precise results for K are provided by the solution of the variational problem in § 4, culminating in the exact asymptotic expressions for K in (4.51) and (4.52). These result in C no slip Nu = 0.0525 and C no stress It is remarkable that the relatively crude estimates in (3.20) and (3.21) are so close to (3.22). The bounds in (3.22) are based on the cosh µz comparison function introduced in (3.6). Further improvements might be possible with other comparison functions, or perhaps by strengthening the bound on kinetic energy dissipation in (2.9) so that it agrees more closely with clues provided by the numerical solution of figure 2. The results from § 4, however, show that such improvements will not alter the exponent 1/3 in (3.19): the 1/3 is a consequence of η ∼ µ −1 in (3.14), and we now proceed to show constructively that there are trial functions that achieve η ∼ µ −1 .

The variational problem for η
The quantity η, introduced in (3.9), is defined as where the incompressible vector field v satisfies the same homogeneous boundary conditions as u and the scalar field θ satisfies the same homogeneous boundary conditions as a in (3.7). The comparison function c(x, z; µ) is defined in (3.4) and we are interested in the boundary-layer limit µ → ∞ corresponding to Ra → ∞.
The SKB-type inequalities in appendix A result in Appendix A, however, does not exhibit (v, θ ) that actually achieve the η ∼ µ −1 scaling suggested by the asymptotic inequality in (4.2), and nor do SKB. Thus it is possible that η µ −n with an exponent n that is larger than 1: this would be consistent with inequality (4.2). One might expect that it would be easy to settle this question by guessing a 'trial' (v, θ) that achieves the scaling η ∼ µ −1 suggested by (4.2). We found, however, that 'obvious' guesses at (θ , v) resulted in µ −2 for no-slip boundary conditions and µ −3/2 for no-stress. Via (3.17), these alternative asymptotic inequalities would dramatically lower the 1/3 exponent in (1.1) to 1/5 and 1/4, respectively. Unfortunately the obvious guesses do not at all resemble the optimal (v, θ ) determined in this section.
(Another way to view the inequalities in appendix A is to say that the optimal solution determined there is exhibited as a velocity field with u = v = 0 and with b(z) and w(z) as the functions θ(z) and ω(z) that optimize the functionals in (A 6). The 'one-dimensional' velocity field u = ω(z)ẑ does not satisfy the continuity equation (2.3). In this sense, the inequalities in appendix A do not take full advantage of ∇ · u = 0.) The main result of this section is to establish the asymptotic equality in (3.14) by showing that there are (v, θ), with ∇ · v = 0, that actually produce the µ −1 scaling suggested by (4.2). We also calculate the constant K in (3.14); unlike the constant 0.3485 in (4.2), this K depends on the choice between no-stress and no-slip boundary conditions and provides the strongest version of the upper bound on χ in (3.13).

A variational problem
The straightforward approach to obtaining η in (4.1) is to calculate the maximum of the right-hand side using variational methods. With this hope, one can reformulate the problem posed in (4.1) by introducing the constrained functional In (4.3) the Lagrange multipliers ξ and η ensure normalization of θ and v, and φ(x) enforces incompressibility of v.
The Euler-Lagrange equations corresponding to extremization of (4.3) are Forming θ(4.4) and v · (4.5) one has The identities (4.7) imply that ξ = η, (4.8) and so we can regard the system in (4.4) through to (4.6) as an eigenproblem with eigenvalue η. We take ξ = η > 0 and it follows from (4.7) that we are seeking the most negative value of θv · ∇c , corresponding to the largest positive value of η. In other words, η in (4.1) is the most positive eigenvalue of the Euler-Lagrange equations (4.4) through to (4.6), and from the identities in (4.7) Thus the Lagrange multiplier η in (4.3) is the same as the η in (3.9) and (4.1). We begin by considering two-dimensional solutions with a stream function ψ such as v = (ψ z , 0, −ψ x ). (4.10) With this simplification, the Euler-Lagrange equations in (4.4) through to (4.6) reduce to 11) c x θ z − c z θ x = η∆ 2 ψ. (4.12) In terms of ψ and θ, the functional in (4.3) is (4.13) 883 A41-10

C. B. Rocha and others
We attempted a numerical assault on (4.11) and (4.12). This failed miserably: we show below in § 4.5 that as µ → ∞ the numerical solution develops small-scale oscillations in x with a wavelength ∼µ −1 . The solution is also boundary layered in z. Given the boundary-layer structure of the comparison function c(x, z; µ) in (3.3), the boundary layer in (ψ, θ) was expected and was built into our numerical solution. The fast oscillation in x, however, was not anticipated. Consequently the available x-resolution became inadequate at only moderately large values of µ and it was not possible to make a reliable estimate of η in the limit µ → ∞. But once one is alerted to the existence of rapid oscillations, the Wentzel-Kramers-Brillouin (WKB) method is suggested as a means of attacking (4.11) and (4.12). Moreover, with this insight, it is also easy to manufacture trial functions that produce the asymptotic lower bounds of the form K low b µ < η. (4.14) These trial-function lower bounds are complementary to the upper bound in (4.2) and show definitively that η is proportional to µ −1 as µ → ∞. Moreover, one can systematically judge the quality of trial functions: the best trial function produces the largest value of K low . We found that systematic examination of three trial functions developed the intuition that is required to successfully apply the WKB method to (4.11) and (4.12). Thus we begin by estimating K low with trial functions.
(As for the miserable failure mentioned above: we applied a Fourier series in x to (4.11) and (4.12). With the sinusoidal b s in (2.4), symmetry considerations show that the series has the form ψ =ψ 1 (z) sin kx +ψ 3 (z) sin 3kx + · · · (4.15) θ =θ 0 (z) +θ 2 (z) cos 2kx +θ 4 (z) cos 4kx + · · · (4.16) The most efficient approach is to keep N terms in the ψ-series and N + 1 in the θ-series. The Galerkin method was used to obtain a system of 2N + 1 ordinary differential equations in z, with analytically determined coefficients. We then solved the resulting eigenvalue problem for η using Matlab's bvp5c boundary value solver (a collocation method). For given µ, increasing N eventually results in satisfactory convergence of the eigenvalue η. For example, with µ = 10, N = 5 was enough. But as µ is increased, ever larger values of N are required for convergence, indicating that the solution is developing ever smaller scales in x. With hindsight, we see that (4.15) and (4.16) are attempting to represent a large-µ solution such as the example shown figure 3; this requires N of order 1000. As µ → ∞, η ∼ µ −χ and a main goal is to determine the exponent χ. But the values of µ that were accessible with truncations of (4.15) and (4.16) were too small to reliably indicate the exponent χ. The result χ = 1 follows from numerically assisted application of the WKB approximation in § 4.5.) 4.2. The first trial function The first trial function we consider is ψ 1 = a ψ e αZ Z 2 sin µS and θ 1 = a θ e αZ Z cos µS, Horizontal convection 883 A41-11 In (4.18), τ is the dimensionless surface buoyancy introduced in (3.6) and the disposable constants α and γ are selected to maximize E[v, θ]. The amplitudes a ψ and a θ in (4.17) are determined by enforcing the normalization conditions that ( ψ) 2 = |∇θ| 2 = 1. The trial function ψ 1 in (4.17) satisfies no-slip conditions at the top surface (corresponding to Z = 0). In the limit hµ → ∞ one can safely ignore the bottom boundary conditions where both the comparison function c(x, z; µ) and the trial functions are exponentially small.
The non-obvious aspect of the trial function in (4.17) is the rapid oscillation of cos µS and sin µS in the limit µ → ∞. This rapid oscillation of the variational solution is key to all results in this section. Now we evaluate the functional E[v 1 , θ 1 ] using (4.17). We freely use µ → ∞ to evaluate all x-integrals; typical µ → ∞ simplifications are τ (x) 2 sin 2 µS → 1 2 τ 2 , τ (x) 4 sin 2 µS → 1 2 τ 4 , (4.19a,b) where the overbar denotes an x-average. After some calculation we find that the normalization conditions are where r def = γ /α. We drop the final O(µ) term in (4.21) relative to the much larger µ 3 term.
After the normalization conditions are applied, the functional is (4.23) Thus, as in the upper bound in (A 5), only the vertical advection by ψ 1x is asymptotically important. The leading-order approximation to the functional is then The parameters α and r = γ /α are determined so as to maximize E[v 1 , θ 1 ]. Thus α = 1/2 for all τ . With τ (x) = cos kx we have τ 2 = 1/2 and τ 4 = 3/8. In this case, the optimal value of r = γ /α is 1.01819 and (4.26) The constant K low 1 is smaller by a factor of approximately 5 than the upper-bound constant 0.3458 in (4.2).

The second trial function
To obtain a larger value of K low we consider a second trial function with the same horizontal structure as (ψ 1 , θ 1 ), but with general z-structure ψ 2 = µ −3/2 P(Z) sin µS and θ 2 = µ −1/2 T(Z) cos µS. (4.27a,b) Here S(x) is defined in (4.18). The factors µ −3/2 and µ −1/2 are included so that P and T are order unity as µ → ∞. Using simplifications such as (4.19) and (4.20) one can evaluate the x-integrals in the functional E[v 2 , θ 2 ] to obtain (4.28) Setting the variational derivatives with respect to T and P to zero we obtain the Euler-Lagrange equations. Again one can show that ξ = η and thus we obtain the eigenproblem where the eigenvalue is λ def = b /ηµ = 1/K low 2 . Solving (4.29) and (4.30) numerically with τ = cos kx and no-slip boundary conditions, P(0) = P Z (0) = 0, and optimizing using the parameter γ , we find (4.31) Comparing (4.26) with (4.31), we see that the more elaborate trial function in (4.27) has increased K low by only 3.4 %.

The third trial function
The small difference between K low 1 and K low 2 indicates that the simple vertical structure assumed in (4.17) is already close to optimal for no-slip boundary conditions. Comparison of P(Z) and T(z) with Z 2 e Z/2 and Ze Z/2 confirms this expectation (not shown). This motivates a third trial function in which we use the vertical structure in (4.17) but admit general structure in x ψ 3 = µ −3/2 Z 2 e αZ Q(x) and θ 3 = µ −1/2 Ze αZ J(x). (4.32a,b) Evaluating the z-integrals in (4.3), the functional is Horizontal convection 883 A41-13 FIGURE 3. Lowest no-slip eigenfunction of the third trial function in (4.32); αµ = 200 and α = 1/2 and the solution uses 1024 points on the interval (0, 2π). The magnitude of the eigenfunction is arbitrary. The function τ = cos kx is also shown; the rapidly oscillatory eigenfunction is exponentially small away from the maxima of cos 2 kx.
where the overbar denotes an x-average. Setting the variational derivatives with respect to J and Q to zero produces the Euler-Lagrange equations. Again one can show that ξ = η so that the Euler-Lagrange equations are where the eigenvalue is (4.36)
The recherché eigenfunction at µ = 400 is shown in figure 3: the eigenfunction is concentrated in regions of width ∼µ −1/2 centred on kx = 0 and kx = π, and has rapid oscillations with wavelength scaling as µ −1 µ −1/2 . Over the rest of the range, the eigenfunction is exponentially small. As µ → ∞ the amplitude of the increasingly concentrated eigenfunction must be also be increased in order to maintain the normalization. In the limit the oscillations become increasingly violent, even as the eigenvalue Λ tends to a limiting value, Λ ∞ . Ultimately the pathological eigenfunction is loaded only at the points kx = 0 and kx = π where cos 2 kx = 1. This peculiar limiting structure is also characteristic of the WKB solution of (4.11) and (4.12) (and for those partial differential equations it seems impossible to obtain a numerical solution with µ sufficiently large to glimpse the limit).

The WKB solution
With τ = cos kx and µ → ∞ we can attack (4.34) and (4.35) using the WKB approximation. We make some preliminary simplifications by introducing The WKB ansatz is where E(X) denotes the eikonal function. Substituting (4.41) into (4.39) and (4.40), and retaining only the leading-order terms, one obtains a 2 × 2 set of linear equations where the local wavenumber is q def = E X and Λ ∞ is the µ → ∞ limiting value of Λ(µ). The condition for a nontrivial solution (J,Q) of (4.42) and (4.43) is that the determinant is zero, or 3q 6 + 5q 4 + (5 − 9Λ 2 ∞ cos 2 X)q 2 + 3 = 0. (4.44) In principle, q(X) is determined by solving this bicubic equation. Fortunately, however, the numerical solution has taught us that Λ ∞ is obtained by requiring that (4.44) has a double root where cos 2 X = 1. For Λ < Λ ∞ , (4.44) has two pure imaginary roots and four other complex roots: see figure 4. Now any root q with a non-zero imaginary part produces exponential growth of the WKB solution in (4.41). It is not possible Horizontal convection 883 A41-15 to construct a solution of the form (4.41) that decays away from x = 0 and x = π if all six roots have non-zero imaginary parts. But in figure 4(b), where Λ = Λ ∞ , there is a double root of the bicubic, yielding four purely real roots for q. These four oscillatory solutions can, in principle, be superposed to construct WKB solutions with the appropriate behaviour. We do not attempt this difficult construction: the main point is that the µ → ∞ eigenfunction is increasingly concentrated in small regions surrounding the points where cos 2 X = 1 and Λ ∞ is determined by requiring that (4.44) has a double root for q 2 at those same points. Thus, with cos 2 X → 1 in (4.44), we set the discriminant of the bicubic equation to zero, and so obtain (4.45) with real solutions Λ ∞ = ±1.25073. To maximize K low 3 we again take α = 1/2 in (4.36) and thus the lower-bound constant produced by the third trial function is (4.46) K low 3 is larger than K low 1 and K low 2 in (4.26) and (4.31) by approximately 30 %, and is smaller than the upper-bound constant 0.3458 in (4.2) by a factor of 3.5. 4.5. Solution of the Euler-Lagrange equations in the limit µ → ∞ Guided by the WKB solution in the previous section, we stop playing with trial functions and now aim to solve the full Euler-Lagrange equations in (4.11) and (4.12) in the limit µ → ∞. The WKB ansatz is is the boundary layer coordinate. The eikonal function E depends only on x, and therefore, for example, ψ x = µ −1/2 q cos(µE) Ψ + µ −3/2 sin(µE) Ψ x and ψ z = µ −1/2 sin(µE) Ψ Z , (4.48a,b) where q def = E x ; the local wavenumber of the WKB solution is µq. Substituting into the Euler-Lagrange equations (4.11) and (4.12), and retaining only the dominant terms, we obtain (4.50) where the eigenvalue is λ def = b /(ηµ). The x-dependence in (4.49) and (4.50) is parametric through the function τ (x) and in principle we can solve the system at each value of x to determine q(x). But the eigenvalue λ cannot depend on x and must be obtained using an approach analogous to that in the discussion surrounding (4.42) through to (4.46).
To illustrate this solution we solve (4.49) and (4.50) by shooting from Z = 0 with appropriate boundary conditions (either no-slip or no-stress) and requiring that Ψ , Ψ and Θ decay as Z → −∞. This numerical solution produces a relation between 4.6. The three-dimensional variational problem The results in (4.51) and (4.52) are obtained by considering the two-dimensional variational problem. In this section we show that in the limit µ → ∞, consideration of the three-dimensional (3-D) variational problem does not alter (4.51) and (4.52). The 3-D problem in (4.4) through to (4.6) can also be solved using the WKB ansatz, which in this case is (4.53) Variables denoted by capitals, U, V, etc., are functions of x and Z = µ(z − h). Again the eikonal function E depends only on x. The 'spanwise' or y-wavenumber is µm where m is a constant; the factor of µ is included so that the y-wavenumber has the same scaling as the x-wavenumber q(x) = E x . To improve the appearance of subsequent equations, the definition of Φ in (4.53) includes a factor of b . Substituting into (4.4) through to (4.6), and retaining only the leading-order terms, one finds where λ def = b /(µη) is the eigenvalue and τ (x) is the non-dimensional function introduced in (3.6). Now we observe that there is a transformation that isotropizes (4.54) through to (4.58) in the (x, y)-plane: introducing we see that (4.54) through to (4.58) is equivalent to the previously solved twodimensional problem with an incompressible velocity field corresponding to (Ũ, W). Thus the construction in figure 5 applies, except that the ordinate is nowq. Thus the µ → ∞ solution of the η-variational problem is not unique: any wavevector (q, m) in the (x, y) plane, withq determined via figure 5, produces a solution of the variational problem with the K values in (4.51) and (4.52). This includes the special flow with q = E x = 0 i.e. a flow in the (y, z)-plane, with no velocity along the x-axis.
We have checked this conclusion by constructing a trial function analogous to (ψ 1 , θ 1 ) in (4.17), except that trial stream function ψ is a function of (y, z). Provided that the spanwise wavenumber is ∼µ -as assumed in (4.53) -this flow in the (y, z)-plane achieves the scaling η ∼ µ −1 .
In terms of horizontal convection this conclusion is, perhaps, puzzling: a flow with u = 0 cannot produce enhanced horizontal transport of buoyancy in the x-direction. But in considering this puzzle one falls into the trap of ascribing too much physical significance to the solutions of the variational problem. Hypothetical solutions of the Boussinesq equations that achieve the Nu ∼ Ra 1/3 scaling might share some structural properties with the variational solution without having the degeneracy of the 3-D variational solutions.

Discussion and conclusion
The inequalities developed in appendix A show that as µ → ∞ the variational constant η in (4.1) is less than the comparison-function boundary-layer thickness, µ −1 . But this asymptotic inequality, η µ −1 , does not eliminate the possibility that η scales with a higher power of µ. For example, η ∼ µ −2 is consistent with η µ −1 ; the exponent −2 reduces the exponent in the Nusselt number bound (1.1) from 1/3 to 1/5. The main achievement of this paper is to eliminate that possibility: we show that the asymptotic inequality η µ −1 of appendix A is sharp in the sense that η ∼ µ −1 . We do this by: (i) exhibiting trial functions that achieve η ∼ µ −1 , and (ii) solving the Euler-Lagrange equations in the limit µ → ∞. Both approaches show that incompressible flows that achieve η ∼ µ −1 do so by oscillating rapidly in the horizontal directions with wavelength ∼µ −1 .
The solution of the η-variational problem reduces the constant C Nu in (1.1) by approximately a factor of two from the result of SKB. For example, with the sinusoidal surface buoyancy in (2.4), and no-slip boundary conditions, our C no slip Nu in (3.22) is smaller by a factor of 2.5 than the C SKB Nu in (3.20). Siggers et al. (2004) discussed the physical significance of their bound by using ballpark oceanographic numbers and comparing their bound with the planetary-scale ocean heat flux, estimated by Munk & Wunsch (1998) to be of order 2 × 10 15 W. The upper bound is greater by a factor of approximately 10 3 and SKB speculate 'an improvement of the prefactor . . . might lead to a more telling result'. Thus it is disappointing that 883 A41-18 C. B. Rocha and others the arduous solution of the variational problem in § 4 produces an underwhelming improvement in the bound. The structure of flows that achieve the bound is rather surprising. As anticipated by SKB and Winters & Young (2009) these flows have a boundary-layer thickness scaling as µ −1 ∼ (νκ/b ) 1/3 . The additional feature resulting from the analysis in § 4 is that the boundary-layer stream function must also develop rolls with a horizontal wavelength ∼(νκ/b ) 1/3 : the optimal flow is boundary-layered in z and rapidly oscillatory in x with the same length scale (νκ/b ) 1/3 in both directions. Although the optimal solution has the singular structure indicated in figure 3, the scaling in (1.1) can also achieved by less pathological flows, such as the trial function in (4.17).
We are not aware of any numerical or laboratory studies of HC that come close to achieving the ultimate HC scaling Nu ∼ Ra 1/3 . Instead, as emphasized by Gayen, Griffiths & Hughes (2014), all evidence is consistent with the scaling Nu ∼ Ra 1/5 , first obtained by Rossby (1965) with the assumption of a visco-diffusive boundary layer. At high Ra the numerical solutions in Gayen et al. (2014) show the development of turbulent three-dimensional flow in the boundary layer and clear violation of the viscodiffusive balance assumed by Rossby. The onset of boundary-layer turbulence does produce an enhancement in buoyancy transport. Yet with further increase in Ra, the turbulent solutions return to the Nu ∼ Ra 1/5 scaling, though with a constant C Nu larger than that in the visco-diffusive Rossby regime.
Of course the Ra's accessed by the direct numerical simulations of Gayen et al. (2014) are at most 6 × 10 11 : further regime changes may occur at higher Ra. But the tenacity of Nu ∼ Ra 1/5 , apparently even in the turbulence regime, suggests that a new approach to bounding HC Nusselt number is required. Recent progress on Rayleigh-Bénard bounds suggests approaches that might also work for HC. For example, Whitehead & Doering (2011) show that for two-dimensional Rayleigh-Bénard convection, with no-stress boundary conditions, Nu is bounded from above by Ra 5/12 , which is a significant improvement on well known Ra 1/2 -bounds (see also Wen et al. (2015)). Another possibility is that instead of merely using a projection of the buoyancy equation onto a comparison function -the identity in (3.8) -better bounds might be achieved by constructing optimal flows that satisfy the buoyancy equation and transport buoyancy along the wall i.e. the HC analogue of the optimal wall-to-wall flows of Hassanzadeh, Chini & Doering (2014). In passing from (A 12) to (A 13) we have used the inequality 4 ω 2 z |∇v| 2 from Doering & Constantin (1996). With η defined in (4.1) we have It turns out that (A 14) leads to a slightly tighter upper bound than the corresponding result of SKB -see the discussion surrounding (3.20) and (3.21). Thus the cosh µz comparison function is slightly superior to the piecewise linear comparison function of SKB.