2 Governing equations

For the governing model we adopt the Oberbeck–Boussinesq approximation, in which the fluid has constant kinematic viscosity, ${\it\nu}$ , thermal diffusivity, ${\it\kappa}$ , and coefficient of thermal expansion, ${\it\alpha}$ . We non-dimensionalize length by the layer height, $d$ , time by the thermal diffusion timescale, $d^{2}/{\it\kappa}$ , and temperature by $d^{2}H/{\it\kappa}$ , where $H$ is the heating rate in units of temperature per time. The dimensionless Boussinesq equations governing the velocity $\boldsymbol{u}=(u,v,w)$ , temperature $T$ and pressure $p$ are then (Rayleigh Reference Rayleigh1916; Chandrasekhar Reference Chandrasekhar1981)

The dimensionless control parameters are the Rayleigh and Prandtl numbers defined by

where $g$ is the uniform gravitational acceleration in the $-\hat{\boldsymbol{z}}$ direction. The parameter $R$ is standard in the study of IH convection but differs from the Rayleigh number of RB convection, where the temperature scale comes from the boundary conditions instead of the heating rate.

The dimensionless vertical extent is $-1/2\leqslant z\leqslant 1/2$ . At the boundaries we impose no-slip conditions on the velocity, and we fix the temperatures to zero (without loss of generality), so $\boldsymbol{u},T=0$ at $z=\pm 1/2$ . Both horizontal directions are periodic and have the same aspect ratio, ${\it\Gamma}$ , so the horizontal coordinates are bounded by $0\leqslant x,y<{\it\Gamma}$ . As described below, we choose ${\it\Gamma}$ large enough so that the global quantities of interest are not sensitive to ${\it\Gamma}$ .

3 Integral quantities

We are especially interested in the mean fluid temperature and the fractions of internally produced heat escaping across the top and bottom boundaries. These quantities and additional information about vertical structure are captured by the mean temperature profile, $\overline{T}(z)$ , where an overbar denotes an average over the horizontal directions and infinite time. (Such infinite-time averages are approximated in simulations by long finite times.) When the fluid is static, the equilibrium temperature profile is parabolic: $T=(1-4z^{2})/8$ . When $R$ is too small to sustain convection, all initial conditions evolve toward this static state, which is parameter-independent by virtue of the non-dimensionalization. When there is no constraint on horizontal wavenumbers, convection is guaranteed by linear instability of the static state when $R>37\,325$ (Sparrow, Goldstein & Jonsson Reference Sparrow, Goldstein and Jonsson1964), and subcritical convection is possible at smaller $R$ (Tveitereid Reference Tveitereid1978). Here we simulate large $R$ , well above the onset of motion.

Figure 2. Mean vertical temperature profiles, $\overline{T}(z)$ , in the static state (- - - -) and in 3D simulations with $Pr=1$ and $R=10^{6},10^{7},10^{8},10^{9},10^{10}$ (from right to left).

Figure 2 shows the parabolic temperature profile of the static state – that is, the purely conductive state – along with selected mean temperature profiles, $\overline{T}(z)$ , for the 3D simulations described below. The basic qualitative trends in figure 2 are shared by all temperature profiles reported for past simulations (Peckover & Hutchinson Reference Peckover and Hutchinson1974; Mayinger et al. Reference Mayinger, Jahn, Reineke and Steinberner1975; Straus Reference Straus1976; Emara & Kulacki Reference Emara and Kulacki1980; Grötzbach Reference Grötzbach, Iwasa, Tamai and Wada1989; Wörner et al. Reference Wörner, Schmidt and Grötzbach1997; Goluskin & Spiegel Reference Goluskin and Spiegel2012) and experiments (Kulacki & Goldstein Reference Kulacki and Goldstein1972; Mayinger et al. Reference Mayinger, Jahn, Reineke and Steinberner1975; Ralph et al. Reference Ralph, McGreevy, Peckover, Spalding and Afgan1977; Lee et al. Reference Lee, Lee and Suh2007). As $R$ is raised, convection further assists conduction in carrying heat to the boundaries. This decreases the dimensionless fluid temperature and brings the interior closer to isothermal. Thermal boundary layers form at the top and bottom, but the unstable top layer is thinner than the stable bottom layer, reflecting the fact that the majority of internally produced heat escapes across the top boundary.

To quantify the reduction of temperature by convection, we consider the (dimensionless) mean fluid temperature, $\langle T\rangle$ , where angle brackets denote an average over volume and infinite time. Most researchers have instead considered $\overline{T}_{max}$ , the maximum value of $\overline{T}(z)$ . The two quantities behave similarly, but $\langle T\rangle$ is more amenable to mathematical analysis, whereas $\overline{T}_{max}$ is easier to measure in the laboratory.

The mean temperature for any solution to the model (2.1)–(2.3) has been proven to obey the bounds (Lu, Doering & Busse Reference Lu, Doering and Busse2004)

at large $R$ . The quantity $\langle T\rangle$ saturates its upper bound of $1/12$ only in the static state and typically decreases as $R$ is raised. This decrease can be no faster than $O(R^{-1/3})$ as $R\rightarrow \infty$ . Since $T$ has been non-dimensionalized using the heating rate $H$ , this means that the dimensional mean temperature must grow with $H$ no slower than $O(H^{2/3})$ . In the only data reported for $\langle T\rangle$ , which are from 2D DNS, $\langle T\rangle$ decays proportionally to $R^{-1/5}$ (Goluskin & Spiegel Reference Goluskin and Spiegel2012). Other studies report similar $R$ -dependence for $\overline{T}_{max}$ (Goluskin Reference Goluskin2015). Questions motivating our present study include: how different are $\langle T\rangle$ and $\overline{T}_{max}$ , will a novel scaling of $\langle T\rangle$ be found at larger $R$ , and how does $\langle T\rangle$ depend on $Pr$ ?

To quantify the up–down asymmetry that convection induces we examine $\mathscr{F}_{B}$ , the fraction of internally produced heat that flows outward across the bottom boundary. This fraction is related a priori to the mean convective transport, $\langle wT\rangle$ , by $\mathscr{F}_{B}=1/2-\langle wT\rangle$ (Goluskin & Spiegel Reference Goluskin and Spiegel2012). Likewise, the fraction of heat flowing outward across the top boundary must equal $1/2+\langle wT\rangle$ since the fractions crossing the top and bottom boundaries sum to unity. The quantity $\mathscr{F}_{B}$ is not well understood, seemingly having no analogue in convection that is not penetrative, and only the crude bounds $0<\mathscr{F}_{B}\leqslant 1/2$ have been proven mathematically (Goluskin & Spiegel Reference Goluskin and Spiegel2012).

Heat flows outward across the top and bottom boundaries in equal proportion in the static state, meaning that $\mathscr{F}_{B}=1/2$ . In sustained convection $\mathscr{F}_{B}<1/2$ because the mean work exerted by buoyancy forces, which must be positive, is proportional to $\langle wT\rangle =1/2-\mathscr{F}_{B}$ . In all past studies, $\mathscr{F}_{B}$ decreases slowly after the onset of convection as $R$ is raised over several decades (Goluskin Reference Goluskin2015). At larger $R$ , there is stark disagreement between past studies. In 2D simulations with $Pr=1$ , $\mathscr{F}_{B}$ reaches a minimum value of 0.33 and then increases as $R$ is raised further (Goluskin & Spiegel Reference Goluskin and Spiegel2012). In experiments (Kulacki & Goldstein Reference Kulacki and Goldstein1972) and 3D simulations (Wörner et al. Reference Wörner, Schmidt and Grötzbach1997), on the other hand, $\mathscr{F}_{B}$ continues to decrease at the largest $R$ studied. Questions raised by these findings include: are the differing observations of $\mathscr{F}_{B}$ explained by the differences between 2D and 3D flows, and what is the ultimate limit of $\mathscr{F}_{B}$ as $R\rightarrow \infty$ ?

4 Simulation results

We have simulated the governing equations (2.1)–(2.3) using an energy-conserving finite difference method (Verzicco Reference Verzicco1996; Verzicco & Camussi Reference Verzicco and Camussi2003; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015). Time averages are deemed converged when values of $\langle T\rangle$ and $\langle wT\rangle$ over the full post-transient duration agree with their values over half that duration to within 1 %. One way we have checked spatial resolution is from the integral balances $R\langle wT\rangle =\langle |\boldsymbol{{\rm\nabla}}\boldsymbol{u}|^{2}\rangle$ and $\langle T\rangle =\langle |\boldsymbol{{\rm\nabla}}T|^{2}\rangle$ , derived by integrating $\boldsymbol{u}\boldsymbol{\cdot }$ (2.2) and $T\times$ (2.3), respectively. The balances are satisfied to within 1 % for all simulations reported here. We have also checked resolution by repeating some simulations at lower resolution and repeating every 3D simulation with $Pr=1$ and $R\leqslant 2\times 10^{8}$ using the spectral element code nek5000 (Fisher, Lottes & Kerkemeier Reference Fisher, Lottes and Kerkemeier2016). In each case, the nek5000 values of $\langle T\rangle$ and $\langle wT\rangle$ agree with the finite difference values to within 1 %. The Appendix gives further details on convergence, including time spans, meshes, and the resolution of boundary layers and Kolmogorov length scales. We have chosen aspect ratios large enough such that raising ${\it\Gamma}$ by approximately 50 % changes $\langle T\rangle$ and $\langle wT\rangle$ by less than 1 %. As shown in the Appendix, the necessary ${\it\Gamma}$ values decrease as $R$ is raised, from ${\it\Gamma}={\rm\pi}$ when $R\leqslant 2\times 10^{6}$ to ${\it\Gamma}=1$ when $R\geqslant 2\times 10^{9}$ .

Figure 3. Variation of mean quantities with $R$ and $Pr$ in 2D (▪) and 3D (●) simulations. The fraction of internally produced heat flowing outward across the bottom boundary, $\mathscr{F}_{B}$ , is shown (a) for various $R$ with $Pr=1$ and (b) for various $Pr$ with $R=2\times 10^{7}$ . For the same simulations, the dimensionless mean fluid temperature, $\langle T\rangle$ , is shown (c) for various $R$ (compensated by $R^{-1/5}$ ) and (d) for various $Pr$ . The data shown are tabulated in the Appendix, except for the 2D data in panels (a) and (c), which are from Goluskin & Spiegel (Reference Goluskin and Spiegel2012).

Figure 3(a) shows the $R$ -dependence of the down-flowing heat fraction, $\mathscr{F}_{B}$ , in 2D (▪) and 3D (●) simulations with $Pr=1$ . As $R$ is raised through moderate values, $\mathscr{F}_{B}$ falls at similar rates in 2D and 3D. At larger $R$ , the two cases diverge dramatically. In 2D, $\mathscr{F}_{B}$ stops falling around $R=10^{9}$ and then rises slowly. In 3D, $\mathscr{F}_{B}$ continues to decay up to the largest $R$ simulated. This decay is steady but quite slow, being approximated well for $R\geqslant 5\times 10^{7}$ by $\mathscr{F}_{B}\sim 0.80\,R^{0.055}$ . It is hard to anticipate whether this decay will persist, let alone to justify its rate. The value of $\mathscr{F}_{B}$ can be viewed as the result of two competing effects: buoyancy-driven mixing of the cold top boundary layer helps heat escape from the top, which lowers $\mathscr{F}_{B}$ , while shear-driven mixing of the cold bottom boundary layer helps heat escape from the bottom, which raises $\mathscr{F}_{B}$ . Both effects strengthen as $R$ is increased, so the change in $\mathscr{F}_{B}$ depends on which effect strengthens faster.

We expect that the emergence of a large-scale circulation (LSC) in our 2D simulations is partly why the bottom boundary layer mixes more effectively in 2D (and thus why $\mathscr{F}_{B}$ is larger in 2D). Such LSC is seen often in RB convection in both 2D and 3D (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009), but mean flows can grow especially strong in 2D, where the vorticity of the plumes is entirely aligned with that of the LSC. Here we have seen no LSC in 3D, although one may arise in larger domains. Figure 4 shows profiles of root-mean-square velocities in the horizontal and vertical directions for 2D and 3D simulations at three different $R$ . The velocities are stronger in the 2D simulations than in the 3D ones, and this difference increases as $R$ is raised, as in RB convection (van der Poel, Stevens & Lohse Reference van der Poel, Stevens and Lohse2013). In 3D, the profiles are not only weaker than the 2D ones but are also skewed strongly upward at all $R$ , indicating less convective penetration. In 2D, on the other hand, the horizontal velocity profiles become visibly more symmetric as convection strengthens, suggesting a strong LSC that penetrates to the bottom boundary layer. The greater asymmetry of velocity profiles in 3D is consistent with the greater asymmetry of heat fluxes in 3D that is reflected by smaller values of $\mathscr{F}_{B}$ in figure 3(a,b).

Figure 3(c) shows the dependence of $\langle T\rangle$ on $R$ , compensated by $R^{-1/5}$ . When $R$ is large, $\langle T\rangle$ decays like $R^{-1/5}$ in both 2D and 3D, as reflected by nearly flat trends in the compensated plot. (Fits to data for $R\geqslant 10^{8}$ give $\langle T\rangle \sim 1.13\,R^{-0.20}$ in 2D and $\langle T\rangle \sim 1.11\,R^{-0.20}$ in 3D. The 3D fit, while less good than the 2D fit, is robust; excluding any two data points gives an exponent of either $-0.20$ or $-0.21$ .) The decay exponent of $-1/5$ corresponds to the dimensional mean temperature growing with the heating rate like $H^{4/5}$ .

Figure 4. Mean vertical profiles of root-mean-square horizontal velocity (a–c) and vertical velocity (d–f) for 2D (- - - -) and 3D (——) simulations with $Pr=1$ and $R=10^{6}$ (a,d), $10^{8}$ (b,e), $10^{10}$ (c,f). In the 2D cases, $v\equiv 0$ . The 3D profiles are qualitatively consistent with those reported by Wörner et al. (Reference Wörner, Schmidt and Grötzbach1997).

The maximum of $\overline{T}(z)$ on the interior, $\overline{T}_{max}$ , is necessarily larger than $\langle T\rangle$ , but the two values grow closer as $R$ is raised. For instance, in the 3D simulations $\overline{T}_{max}$ is 28 % larger than $\langle T\rangle$ when $R=10^{6}$ but only 8 % larger when $R=10^{10}$ (cf. table 2 in the Appendix). This reflects the flattening of temperature profiles that is evident in figure 2. The decay of $\overline{T}_{max}$ is thus slightly faster than that of $\langle T\rangle$ , being fitted well by $\overline{T}_{max}\sim 1.62\,R^{-0.22}$ in 3D. This exponent is consistent with past experiments (Kulacki & Goldstein Reference Kulacki and Goldstein1972; Jahn & Reineke Reference Jahn and Reineke1974; Mayinger et al. Reference Mayinger, Jahn, Reineke and Steinberner1975; Ralph et al. Reference Ralph, McGreevy, Peckover, Spalding and Afgan1977; Lee et al. Reference Lee, Lee and Suh2007) and 3D simulations (Wörner et al. Reference Wörner, Schmidt and Grötzbach1997), where $\overline{T}_{max}$ falls at rates between $R^{-0.18}$ and $R^{-0.22}$ (Goluskin Reference Goluskin2015). At very large $R$ , we expect $\overline{T}_{max}\approx \langle T\rangle$ once $\overline{T}(z)$ becomes very flat outside of negligibly thin boundary layers.

Figures 3(b) and 3(d) show the $Pr$ -dependence of $\mathscr{F}_{B}$ and $\langle T\rangle$ , respectively, with $R=2\times 10^{7}$ . When $Pr\rightarrow 0$ , deviations from the static temperature profile cannot be maintained, so we expect $\mathscr{F}_{B}$ and $\langle T\rangle$ to approach their static values in both 2D and 3D, much like the Nusselt number in RB convection. When $Pr\rightarrow \infty$ , we expect both quantities to saturate at intermediate values as solutions of the Boussinesq equations approach solutions of the infinite- $Pr$ Boussinesq equations (Wang Reference Wang2004, Reference Wang2008). The findings of Schmalzl, Breuer & Hansen (Reference Schmalzl, Breuer and Hansen2004) on RB convection suggest that the infinite- $Pr$ limits of mean quantities will be similar in 2D and 3D. Behaviour is more complicated at intermediate $Pr$ , where the deviation between 2D and 3D flows is expected to be largest. Although $\mathscr{F}_{B}$ varies monotonically over the $Pr$ range simulated, $\langle T\rangle$ does not. Below we relate this non-monotonicity of $\langle T\rangle$ to that of the RB Nusselt number.

Quantities like $\mathscr{F}_{B}$ that measure vertical asymmetry seem to have no analogues in RB convection. On the other hand, the dimensionless ratio $1/\langle T\rangle$ , which is proportional to the ratio of the heating rate $H$ to the dimensional mean temperature, behaves much like the Nusselt number in RB convection. This motivates us to define a Nusselt number, $N$ , for IH convection and consider its dependence on a diagnostic Rayleigh number, $Ra$ , that differs from the control parameter  $R$ .

The quantities $N$ and $Ra$ should be defined to convey the strengths of convective transport and thermal forcing, respectively. When comparing RB models with different thermal boundary conditions, it is most useful to define $N$ as the ratio of total vertical heat flux to conductive vertical heat flux in the flow, and to define $Ra$ like $R$ but using a temperature scale that is characteristic of the developed flow (Otero et al. Reference Otero, Wittenberg, Worthing and Doering2002; Johnston & Doering Reference Johnston and Doering2009; Wittenberg Reference Wittenberg2010). In the present model, we cannot normalize $N$ using the mean vertical conduction, which is zero, but we can instead consider outward heat flux: upward flux above the height where $\overline{T}_{max}$ occurs, plus downward flux below it (Goluskin & Spiegel Reference Goluskin and Spiegel2012). Total outward flux is fixed on average, but conductive outward flux is proportional to $\overline{T}_{max}$ , so $N$ should be proportional to $1/\overline{T}_{max}$ . The dimensional version of $\overline{T}_{max}$ can be used also as the temperature scale defining $Ra$ , yielding

where $N=1$ and $Ra=R$ in the static state. Alternatively, choosing $\langle T\rangle$ as the temperature scale instead of $\overline{T}_{max}$ gives the similar definitions

Further discussion of how to define Nusselt-number-like quantities in IH convection is given by Goluskin (Reference Goluskin2015).

Nusselt and Rayleigh numbers defined according to (4.1) or (4.2) display a parameter dependence similar to the Nusselt number in RB convection. For instance, the fits to $\overline{T}_{max}$ and $\langle T\rangle$ reported above for 3D simulations with $Pr=1$ , when formulated in terms of Nusselt numbers, become $N\sim 0.038\,Ra^{0.28}$ and $\widetilde{N}\sim 0.039\,\widetilde{Ra}^{0.26}$ , respectively. These exponents are within the ranges seen in RB experiments with similar parameters, although they are slightly smaller than usual (Grossmann & Lohse Reference Grossmann and Lohse2000; Ahlers et al. Reference Ahlers, Grossmann and Lohse2009). The approximate scaling of $\langle T\rangle$ like $R^{-1/5}$ corresponds to $\widetilde{N}$ scaling like $\widetilde{Ra}^{1/4}$ . The closeness of the data to this scaling is shown by figure 5(a), which re-expresses the $\langle T\rangle$ data of figure 3(c) in terms of (compensated) $\widetilde{N}$ and $\widetilde{Ra}$ . Moreover, scaling arguments like those of Malkus (Reference Malkus1954), Kraichnan (Reference Kraichnan1962), Spiegel (Reference Spiegel1971) or Grossmann & Lohse (Reference Grossmann and Lohse2000) predict the same Nusselt number scalings for IH convection as for RB convection, provided that Nusselt and Rayleigh numbers are defined by (4.1) or (4.2). Even mathematical bounds share the same scalings; the bounds (3.1) on $\langle T\rangle$ correspond to $1\leqslant \widetilde{N}<0.025\widetilde{Ra}^{1/2}$ , and upper bounds with this same exponent have been proven for 3D RB convection with various boundary conditions (Constantin & Doering Reference Constantin and Doering1996; Otero et al. Reference Otero, Wittenberg, Worthing and Doering2002; Plasting & Kerswell Reference Plasting and Kerswell2003; Wittenberg Reference Wittenberg2010). Finally, the parallels with the RB Nusselt number extend also to the $Pr$ -dependence of $\widetilde{N}$ , shown in figure 5(b). Both the 2D and 3D trends in this figure resemble the analogous results for RB convection (cf. figure 6b of van der Poel et al. (Reference van der Poel, Stevens and Lohse2013)).

Figure 5. Variation of $\widetilde{N}$ (a) with $\widetilde{Ra}$ (compensated by $\widetilde{Ra}^{1/4}$ ) and (b) with $Pr$ in 2D (▪) and 3D (●) simulations. The $\widetilde{N}$ and $\widetilde{Ra}$ values are calculated according to (4.2a,b ) from the $\langle T\rangle$ values represented in figure 3(c,d).