RETRACTED – Reynolds-number scaling and convergence time scale in two-dimensional Rayleigh–Bénard convection

Abstract

An equation for the evolution of mean kinetic energy, $E$, in a two-dimensional (2-D) or 3-D Rayleigh–Bénard system with domain height $L$ is derived. Assuming classical Nusselt-number scaling, $Nu \sim Ra^{1/3}$, and that mean enstrophy, in the absence of a downscale energy cascade, scales as $\sim E/L^2$, we find that the Reynolds number scales as $Re \sim Pr^{-1}Ra^{2/3}$ in the 2-D system, where $Ra$ is the Rayleigh number and $Pr$ the Prandtl number. Using the evolution equation and the Reynolds-number scaling, it is shown that $\tilde {\tau } \gtrsim Pr^{-1/2}Ra^{1/2}$, where $\tilde {\tau }$ is the non-dimensional convergence time scale. For the 3-D system, we make the estimate $\tilde {\tau } \gtrsim Ra^{1/6}$ for $Pr = 1$. It is estimated that the total computational cost of reaching the high $Ra$ limit in a simulation is comparable between two and three dimensions. The predictions are compared with data from direct numerical simulations.

1. Introduction

The problem of scaling in Rayleigh–Bénard convection has a long history represented by a huge body of literature. In this introduction, we will not make any attempt to review this literature but concentrate on an issue which is of great relevance for the current debate, namely the differences and similarities between the two-dimensional (2-D) and 3-D Rayleigh–Bénard systems. For general reviews on the problem the reader is referred to Siggia (1994) and Chillà & Schumacher (2012), and for a review on the current debate the reader is referred to Lindborg (2023). The essay by Doering (2020) is also highly recommended.

Despite the fundamental difference between 2-D and 3-D turbulence, a clear distinction is rarely made between two and three dimensions in the theoretical discussion of Rayleigh–Bénard convection. The arguments for and against classical Nusselt-number$/$Rayleigh-number scaling, $Nu \sim Ra^{1/3}$ (Malkus 1954), and non-classical scaling, $Nu \sim Ra^{1/2}$ (Spiegel 1963), $Nu \sim Ra^{1/2} [ \ln (Ra)]^ {-3/2}$ (Kraichnan 1962) or $Nu \sim Ra^{2/7}$ (Castaing et al. 1989), are most often discussed as if they were equally relevant for the 2-D and 3-D systems. Making a series of 2-D direct numerical simulations (DNS), Zhu et al. (2018) claim that they found evidence of a transition to the so-called ‘ultimate regime’ of heat transfer predicted by Kraichnan (1962) and in modified form by Grossmann & Lohse (2011). The theory of Kraichnan (1962) is based on the assumption that the boundary layers in a convection cell in the limit of high $Ra$ are of the same type as classical turbulent boundary layers in a 3-D shear flow and that the friction law of such a boundary layer can be approximated as

(1.1)\begin{equation} u_{\tau}= {\bar{\kappa}} {U} [\ln(Re)]^{{-}1} , \end{equation}

where $u_\tau$ is the friction velocity, $U$ is the free-stream velocity and $\bar {\kappa }$ is the Kármán constant. The modified theory of Grossmann & Lohse (2011) is based on the same general assumption but uses a slightly modified form of (1.1). In order to apply either the Kraichnan (1962) or the Grossmann & Lohse (2011) theory to the 2-D case it must be assumed that a friction law of this type is valid also in two dimensions. However, there is no reason to believe this, because 2-D and 3-D turbulence are fundamentally different. Falkovich & Vladimirova (2018) investigated 2-D Couette and Poiseuille flow and presented strong numerical and theoretical evidence showing that 2-D Couette flow will never become turbulent no matter how large the Reynolds number is, while the high-Reynolds-number Poiseuille flow indeed is turbulent but exhibits a completely different boundary layer structure than the corresponding 3-D flow, with a friction law of the form $u_{\tau } \sim U Re^{-1/4}$. This friction law was also found in 2-D turbulence experiments in soap films by Tran et al. (2010). Zhu et al. (2018) claim that their simulations, indeed, show that the boundary layers in 2-D Rayleigh–Bénard convection are of the same form as in a 3-D shear flow, including a logarithmic dependence of an appropriately defined mean velocity, $U$. However, instead of the classical logarithmic law, they claim that they have discovered a law of the form

(1.2) \begin{equation} U(z) = u_{\tau} \Big ( \frac{1}{\bar{\kappa}} \ln \left (\frac{u_{\tau} z}{\nu} \right ) + B(Re) \Big), \end{equation}

where $z$ is the distance from the wall, $\nu$ is the kinematic viscosity and $B(Re)$ is an increasing function of $Re$, instead of being a constant, as in a standard 3-D boundary layer. It is clear that such a Reynolds-number dependence is not consistent with (1.1) or the corresponding friction law used by Grossmann & Lohse (2011), which are both derived under the assumption that $B$ is a constant. This is not discussed by Zhu et al. (2018). Doering, Toppoladoddi & Wettlaufer (2019) criticise the claim of Zhu et al. (2018) and conclude that their data are more consistent with classical Nusselt-number scaling in accordance with the theory of Malkus (1954) than with the prediction of Kraichnan (1962). However, instead of questioning the consistency of the claim, they question the curve fitting procedure of Zhu et al. (2018). In a reply, Zhu et al. (2019) answer that they have carried out two more simulations at even higher $Ra$ and that their curve fit now shows ‘overwhelming evidence’ of a transition to the ultimate regime. The discussion does not address the underlying theoretical assumption that the 2-D and 3-D cases can be analysed in the same way.

As for the Reynolds-number scaling, it is often observed in numerical simulations that $Re$ has a stronger scaling with $Ra$ in two dimensions as compared with three dimensions, with $Re \sim Ra^{\beta }$, where $\beta > 1/2$ instead of $\beta < 1/2$, as is most often observed in three dimensions. van der Poel, Stevens & Lohse (2013) note this but do not analyse the reason for the difference or discuss if it will prevail in the limit of high $Ra$. Exact steady but unstable solutions to the 2-D problem have been analysed by Chini & Cox (2009), Wen et al. (2020) for stress-free boundary conditions and by Waleffe, Boonkasame & Smith (2015) and Wen, Goluskin & Doering (2022) for no-slip boundary conditions. The solutions show classical Nusselt-number scaling, while the Reynolds-number scaling is quite different between the two cases. For stress-free boundary conditions a clean scaling of the form $Re \sim Pr^{-1} Ra^{2/3}$ is obtained where Pr is the Prandtl number. For no-slip boundary conditions, on the other hand, Wen et al. (2022) find that the scaling is much weaker in $Ra$, with $\beta = 0.47$ for rolls with $Nu$-maximising aspect ratios. Without addressing the issue whether the same type of Reynolds-number scaling is to be expected in two and three dimensions, Wen et al. (2022) point out that $\beta = 0.47$ is close to the value $\beta = 0.46$ which was observed in 3-D DNS (Iyer et al. 2020) up to $Ra = 10^{15}$. These are a few examples of studies in which the issue whether 2-D and 3-D Rayleigh–Bénard convection are fundamentally different is left unaddressed.

In a previous paper (Lindborg 2023), it was argued that the velocity and thermal boundary layer widths in the 3-D system are proportional to the smallest length scales that can develop in 3-D turbulence, that are the Kolmogorov and Batchelor scales. Using this assumption, it was shown that classical Nusselt-number scaling is recovered in the limit of high $Ra$. Moreover, the fundamental scaling relation of 3-D turbulence was used to deduce the Reynolds-number scaling. This relation states that

(1.3)\begin{equation} \frac{\epsilon}{E^{3/2}/L} = {\mbox{const.}} \quad {\mbox{as}} \quad Re = \frac{E^{1/2} L}{\nu} \rightarrow \infty , \end{equation}

where $\epsilon$ is mean kinetic energy dissipation, $E$ is mean turbulent kinetic energy and $L$ is the turbulent integral length scale. The relation (1.3) is most often written using $u^{3}$ instead of $E^{3/2}$, where $u$ is a characteristic turbulent velocity. It has been experimentally verified in a wide range of turbulent flows (Sreenivasan 1998) and numerically verified in 3-D Rayleigh–Bénard convection (Pandey et al. 2022). Using $Nu \sim Ra^{1/3}$ together with (1.3), Lindborg (2023) derived the Reynolds-number scaling

(1.4)\begin{equation} Re \sim Pr^{{-}2/3} Ra^{4/9}, \end{equation}

for the 3-D system. The relation (1.4) was previously derived by a number of other investigators (e.g. Kraichnan 1962; Siggia 1994) under different assumptions. The scaling (1.4) with respect to Rayleigh number, which is our main focus, is in very good agreement with experimental and numerical data. Ashkenazi & Steinberg (1999) report $Re \sim Ra^{0.43}$ from high-Rayleigh-number experiments and Iyer et al. (2020) report $Re \sim Ra^{0.46}$ from DNS of convection in a low aspect ratio cylindrical cell at $Ra \in [10^{9}, 10^{15}]$. It is also interesting to note that the exact steady 3-D solution which recently was found by Motoki, Kawahara & Shimizu (2021) exhibits $Ra^{4/9}$-scaling, while the 2-D solution derived by the same authors exhibits a considerably faster increase of $Re$ with $Ra$.

For the 2-D system with stress-free boundary conditions Whitehead & Doering (2011) rigorously proved that the Nusselt number is limited by $0.2891 Ra^{5/12}$, ruling out $Ra^{1/2}$-scaling, including possible logarithmic corrections. For the 2-D system with no-slip boundary conditions, there is, as far as the author can see, no reason to believe that the theory of Kraichnan (1962) or Grossmann & Lohse (2011) is applicable, even though no proof of this has yet been presented. Therefore, we will assume that classical Nusselt-number scaling holds in two dimensions and focus on the Reynolds-number scaling and the convergence time scale. There is reason to believe that dissipation is much weaker in two than in three dimensions and that (1.3) does not hold in two dimensions. As will be argued, the weaker dissipation implies that the Reynolds-number scaling will be different in two than in three dimensions and that the 2-D system will converge on a much longer time scale than the 3-D system. To deduce the time scale we will first derive an equation for the evolution of mean kinetic energy.

2. Evolution equation for the mean kinetic energy

We assume that the flow is described by the Navier–Stokes equations under the Boussinesq approximation

(2.1)\begin{gather} \frac{{\mbox{D}} \boldsymbol{u}}{{\mbox{D}} t} ={-} \frac{1}{\rho} \boldsymbol{\nabla} p + g \alpha T {\boldsymbol{e}}_z + \nu \nabla^2 \boldsymbol{u} , \end{gather}

(2.2)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0 , \end{gather}

(2.3)\begin{gather}\frac{{\mbox{D}} T }{{\mbox{D}} t} = \kappa \nabla^2 T , \end{gather}

where $\rho, p, g, \nu$ and $\kappa$ are density, pressure, acceleration due to gravity, kinematic viscosity and diffusivity, ${\boldsymbol {e} }_z$ is the vertical unit vector, $T$ is temperature and $\alpha$ is the thermal expansion coefficient. We locate a lower boundary at $z = -L/2$, and an upper boundary at $z = L/2$, lateral boundaries at $x = -X_0$ and $x = X_0$. In three dimensions we also introduce lateral boundaries at $y = -X_0$ and $y = X_0$. We assume that $X_0 \gg L$, so that horizontal mean values will be well converged. We apply constant temperature boundary conditions at the lower and upper boundaries, with $T = \Delta T/2$ and $T = -\Delta T/2$ at the lower and upper boundaries, respectively. At the lateral boundaries we apply adiabatic boundary conditions. For the velocity field we apply stress-free or no-slip boundary conditions. We assume that the initial temperature is linear, $T = -\Delta T z/L$, and that the initial velocity field is very close to zero. The non-dimensional input parameters of the problem are the Rayleigh and Prandtl numbers, defined as

(2.4a,b)\begin{equation} Ra = \frac{g \alpha \Delta T L^3}{\nu \kappa} , \quad Pr = \frac{\nu}{\kappa} , \end{equation}

while the output parameters are the time-dependent Nusselt and Reynolds numbers defined as

(2.5a,b)\begin{equation} Nu =\left. -\frac{{\mbox{d}} \bar{T} }{{\mbox{d}} z}\right|_{z={-}L/2} / {( \Delta T /L )} , \quad Re = \frac{E^{1/2} L}{\nu}, \end{equation}

where $E$ is the domain mean value of the kinetic energy per unit mass, and the bar denotes a horizontal mean.

Using Cartesian tensor notation, the kinetic energy equation can be written as

(2.6)\begin{equation} \frac{1}{2}\,\frac{\partial u_i u_i}{\partial t} ={-} \frac{\partial }{\partial x_j} \left (u_j \left (\frac{1}{2} u_i u_i + \frac{p}{\rho} \right ) \right ) + g \alpha T w - 2\nu S_{ij} S_{ij} + 2\nu \frac{\partial }{\partial x_j} (u_i S_{ij} ) , \end{equation}

where $w$ is the vertical velocity and $S_{ij}$ is the strain rate tensor. In turbulence theory, dissipation is often expressed in terms of vorticity, ${\boldsymbol {\omega }}$, rather than strain. For an incompressible fluid such a formulation is perfectly consistent, which can be seen from the identity

(2.7)\begin{equation} 2\nu S_{ij} S_{ij} = \nu \omega_i \omega_ i + 2\nu \frac{\partial }{\partial x_j} \left ( u_i \frac{\partial u_j}{\partial x_i} \right ) . \end{equation}

Since the last term will integrate to zero over a volume with no-slip or stress-free boundary conditions, dissipation can be defined as $\nu \omega ^2$ instead of $2 \nu S_{ij} S_{ij}$. This definition is preferable in the context of 2-D Rayleigh–Bénard convection, for two reasons. First, conservation of enstrophy, $\omega ^2/2$, is central in the theory of 2-D turbulence. To be able to link dissipation to enstrophy has therefore certain theoretical advantages. Second, linking dissipation to vorticity will clarify the difference between stress-free and no-slip boundary conditions. With stress-free conditions, vorticity is zero at the boundaries in two dimensions which is not generally true for $S_{ij} S_{ij}$. A vorticity based definition will thus guarantee that boundary layer dissipation is small with stress-free conditions as opposed to the case with no-slip conditions.

To derive the expression for the evolution of $E$ we integrate the temperature equation (2.2) to obtain

(2.8)\begin{equation} \overline{wT} ={-} \int_{{-}L/2}^{z} \frac{\partial \bar{T}}{\partial t} \, {\mbox{d}}z + \kappa Nu \frac{\Delta T}{L} + \kappa \frac{\partial \bar{T}}{\partial z} . \end{equation}

Integrating (2.6) over the whole domain, using (2.7) and (2.8), assuming that $\bar {T}$ remains an odd function of $z$ during the evolution of the flow and integrating in time with given initial conditions, we find

(2.9)\begin{equation} E(t) = \alpha g \int_{{-}L/2}^{L/2} \frac{z}{L} \left ( \bar{T}(z, t) + \frac{z}{L} \Delta T \right ) {\mbox{d}} z + \int_0^t \left ( \frac{\kappa^2 \nu}{L^4} Ra (Nu(t)-1) - \epsilon(t) \right ) {\mbox{d}}t , \end{equation}

where $\epsilon$ is the mean dissipation, which in two dimensions can be expressed as

(2.10)\begin{equation} \epsilon = \frac{1}{2LX_0} \int_{{-}X_0}^{X_0} \int_{{-}L/2}^{L/2} \nu \omega^2 \, {\mbox{d}} z \,{\mbox{d}}\kern0.7pt x , \end{equation}

with a corresponding expression in three dimensions. The first term on the right-hand side of (2.9) arises from the first term on the right-hand side of (2.8) in the following way:

(2.11)\begin{align} - \int_0^t \frac{\alpha g}{L} \int_{{-}L/2}^{L/2} \int_{{-}L/2}^{z} \frac{\partial \bar{T}}{\partial t} \, {\mbox{d}}z^\prime\, {\mbox{d}}z\,{\mbox{d}} t &={-} \frac{\alpha g}{L} \int_{{-}L/2}^{L/2} \int_{{-}L/2}^{z} (\bar{T}(z^\prime, t)- \bar{T}(z^\prime, 0)) \, {\mbox{d}}z^\prime \,{\mbox{d}}z \nonumber\\ &={-} \frac{\alpha g}{L} \left [ z \int_{{-}L/2}^{z } (\bar{T}(z^\prime, t)- \bar{T}(z^\prime, 0)) \, {\mbox{d}}z^\prime \right ]_{z ={-}L/2}^{z = L/2} \nonumber\\ &\quad + \frac{\alpha g}{L} \int_{{-}L/2}^{L/2} z (\bar{T}(z, t)- \bar{T}(z, 0))\,{\mbox{d}} z \nonumber\\ &=\alpha g \int_{{-}L/2}^{L/2} \frac{z}{L} \left ( \bar{T}(z, t) + \frac{z}{L} \Delta T \right){\mbox{d}} z , \end{align}

where it has been assumed that $\bar {T}(z,t)$ is an odd function of $z$ and that $\bar {T}(z,0) = - \Delta T z/L$. Introducing the free-fall velocity, $u_f = \sqrt {gL\alpha \Delta T}$, and the non-dimensional variables

(2.12a–e)\begin{equation} \tilde{E} = \frac{E}{u_f^2} , \quad \tilde{\epsilon} = \frac{ L \epsilon}{u_f^3} , \quad \tilde{t} = \frac{u_f t} {L} , \quad \tilde{\bar{T}} = \frac{\bar{T}}{\Delta T} , \quad \tilde{z}= \frac{z}{L} , \end{equation}

(2.9) can be written in non-dimensional form as

(2.13)\begin{equation} \tilde{E}(\tilde{t}) = \int_{{-}1/2}^{1/2} \tilde{z} ( \tilde{\bar{T}}(\tilde{z}, \tilde{t})+ \tilde{z}) \, {\mbox{d}} \tilde{z} + \int_0^{\tilde{t}} [ Pr^{{-}1/2} Ra^{{-}1/2} (Nu(\tilde{t})-1) - \tilde{\epsilon}(\tilde{t}) ] \, {\mbox{d}} \tilde{t} . \end{equation}

From (2.13) it follows that, in the stationary state, we will have

(2.14)\begin{equation} \tilde{\epsilon} = Pr^{{-}1/2} Ra^{{-}1/2} (Nu -1) , \end{equation}

which previously has been shown in many studies.

3. Reynolds-number scaling and convergence time scale

The key property distinguishing 2-D from 3-D turbulence is the conservation of enstrophy by the nonlinear term. The equation for mean enstrophy, $\varOmega = \langle \omega ^2 \rangle /2$, can be written as

(3.1)\begin{equation} \frac{\partial \varOmega }{\partial t} = Q-\epsilon_{\omega} + \frac{1}{2L X_0} \int_{boundaries} \nu {\boldsymbol{n}} \boldsymbol{\cdot}\boldsymbol{\nabla} \left ( \frac{\omega^2}{2} \right ) {\mbox{d}} s , \end{equation}

where $Q = -g \alpha \langle \omega \partial _x T \rangle$ is mean enstrophy production by buoyancy, $\epsilon _{\omega } = \nu \langle \partial _i \omega \partial _i \omega \rangle$ is mean enstrophy dissipation and the last term (where $\boldsymbol n$ is the outwards pointing normal unit vector) is enstrophy production at the boundaries.

With stress-free boundary conditions, the last term is zero, because $\omega$ is zero at the boundaries. For the same reason, kinetic energy dissipation, defined as $\nu \omega ^2$, is zero at the boundaries and the contribution to total dissipation from the boundaries is negligible. As shown by Fjørtoft (1953) and Kraichnan (1967) there can be no downscale energy cascade in a system where enstrophy is conserved by the nonlinear term. Instead, energy will tend to cascade to larger scales. There is probably no extended inverse cascade range in the 2-D Rayleigh–Bénard system, since the energy injection scale is most likely proportional to $L$. In three dimensions, the kinetic energy spectrum peaks at a wavenumber that is inversely proportional to $L$ and falls off as $\sim k^{-5/3}$ at higher wavenumbers (e.g. Pandey et al. 2022). Most likely, the energy injection scale in the 2-D system is also proportional to $L$ but the energy spectrum falls off as $\sim k^{-3}$ at high wavenumbers, as predicted by Kraichnan (1967). Even though there is no extended inverse cascade range, energy and enstrophy will accumulate in the lowest-order modes. Convection rolls that extend over the whole domain can be seen as a manifestation of this. If the production term, $Q$, in (3.1) is dominated by the large-scale modes, mean enstrophy will, in the absence of a downscale energy cascade, scale as $\varOmega \sim E/L^2$. This is the only assumption that is needed in order to derive the Reynolds-number scaling. Kinetic energy dissipation then scales as

(3.2)\begin{equation} \epsilon = 2 \nu \varOmega \sim \nu \frac{E}{L^2} \quad \Rightarrow \quad \frac{\epsilon}{E^{3/2}/L} \sim Re^{{-}1} . \end{equation}

In the 2-D system, the scaling of dissipation is thus strikingly different from the scaling (1.3) in the 3-D system. Although the argument for (3.2) seems very strong in the case with stress-free boundary conditions it can be questioned that (3.2) also holds in the case with no-slip conditions. The production of enstrophy at the boundaries by the last term in (3.1) could give rise to a considerably larger dissipation in the boundary layers than in the central region. However, DNS data by Zhang, Zhou & Sun (2017, figure 7) show that the ratio between the boundary layer dissipation and the central region dissipation is of the order of unity and is independent of $Ra$ at $Ra \in [10^6, 10^{10}]$. It seems unlikely that this ratio should increase dramatically for higher $Ra$. We therefore assume that (3.2) also holds in the case with no-slip conditions. The data of Zhang et al. (2017, figure 5a,b) show that dissipation is more or less constant in the central region, with a sharp increase at the edge of the boundary layers. Exact solutions with stress-free boundary conditions (Chini & Cox 2009) also show that enstrophy is constant in the central region but with a sharp decrease at the edge of the boundary layers. In both cases, it can be assumed that the thermal boundary layer width, $\delta _T$, is determined only by $\kappa$ and $\varOmega$. Dimensional considerations then give $\delta _T \sim \kappa ^{1/2} \varOmega ^{-1/4} \sim \kappa ^{1/2} \nu ^{1/4} \epsilon ^{-1/4}$, which is the Batchelor scale (Batchelor 1959) calculated using the mean dissipation of the system. As pointed out by Lindborg (2023) this is exactly the scale of $\delta _T$ that corresponds to classical Nusselt-number scaling.

Using (2.14), with $Nu-1 \approx Nu$, (3.2) and $Nu \sim Ra^{1/3}$, we obtain

(3.3)\begin{gather} \tilde{E} \sim Pr^{{-}1} Nu \sim Pr^{{-}1} Ra^{1/3} , \end{gather}

(3.4)\begin{gather}Re \sim Pr^{{-}1} Ra^{1/2} Nu^{1/2} \sim Pr^{{-}1} Ra^{2/3} , \end{gather}

which are two equivalent expressions of the same identity. The Reynolds-number scaling (3.4) is identical to the expression derived for exact solutions with stress-free boundary conditions by Wen et al. (2020). Substituting (3.3) into the left-hand side of (2.13) we can determine a minimum time it will take for the system to settle. An upper limit of the first term on the right-hand side of (2.13) can be estimated by putting $\tilde {\bar {T}} = 0$ in the integral, which in this case will be equal to $1/12$. This estimate is not crucially dependent on the assumption that the initial temperature profile is linear. The closer the initial temperature profile is to the final stationary profile, the smaller is this term. The first term on the right-hand side of (2.13) is thus negligible in comparison with the left-hand side which is of the order of $Pr^{-1} Nu \gg 1$. Assuming that $Nu(\tilde {t}) < C Nu_{st}$ where $C$ is a constant and $Nu_{st}$ is the Nusselt number in the stationary state, we obtain

(3.5)\begin{equation} \tilde{\tau} \gtrsim Pr^{{-}1/2}Ra^{1/2} . \end{equation}

It is noteworthy that (3.5) is derived using only (2.13), (3.2) and the assumption $Nu(\tilde {t}) < C Nu_{st}$, but no assumption regarding the Nusselt number scaling. In dimensional form (3.5) is expressed as $\tau \gtrsim L^2/\nu$, which is a general expression for the convergence time scale of a 2-D system undergoing an inverse energy cascade.

For the 3-D system it is not as straightforward to estimate $\tilde {\tau }$. Most likely, $\tilde {\tau }$ increases with $Ra$ also in the 3-D system (private communication with J. Schumacher), although not as fast as in the 2-D system. In three dimensions the left-hand side of (2.9) is of the order of $Ra^{-1/9} Pr^{-1/3}$. Formally, it is therefore a subleading term in the limit of high $Ra$ since the first term on the right-hand side is independent of $Ra$ and therefore, formally, is of the order unity, although it is smaller than $1/12$. The Prandtl-number dependence complicates the matter and we therefore only make an estimate for $Pr = 1$. In this case stationarity cannot be reached until the dissipation term is of the order of unity. Assuming that mean dissipation, $\tilde {\epsilon }$, during the evolution of the flow is limited by the value it takes in the final state we obtain $\tilde {\tau } \gtrsim Ra^{1/6}$ for the 3-D system.

The total computational cost in a simulation is proportional to the number of grid points multiplied by the number of time steps. For $Pr \sim 1$, the temperature and velocity fields should be resolved at the Kolmogorov scale, which will require that the number of grid points in each direction is proportional to $Ra^{1/3}$. Assuming that $\tilde {\tau } \sim Ra^{1/2}$ in two dimensions and $\tilde {\tau } \sim Ra^{1/6}$ in three dimensions the total computational cost will scale as $Ra^{7/6} \times N_{\Delta t}$ in both cases, where $N_{\Delta t}$ is the number of time steps per non-dimensional unit time, that surely is as large in two dimensions as in three dimensions. As a matter of fact, it is likely that $N_{\Delta T}$ will be larger in two dimensions than in three dimensions. A Courant condition based on the magnitude of the velocity will require a smaller time step in two than in three dimensions. The total cost will thus scale at least as fast with $Ra$ in two dimensions as in three dimensions. At $Ra = 10^{14}$, it may actually be more costly to run a fully resolved simulation to a stationary state in two dimensions than in three dimensions.

4. Comparison with DNS data

We first consider the $Nu$ and $Re$ scaling for the case with stress-free boundary conditions and then move to the case with no-slip conditions. In each case, we consider (i) $\gamma$ in $Nu \sim Ra^{\gamma }$, (ii) whether $Nu$ is independent of $Pr$, (iii) $\beta$ in $Re \sim Pr^{\sigma } Ra^{\beta }$ and (iv) $\sigma$ in the same expression. Finally, we consider evidence of the scaling of the convergence time scale.

4.1. Stress-free boundary conditions

Wang et al. (2020a) made an extensive DNS study of 2-D Rayleigh–Bénard convection with stress-free boundary conditions and studied the flow evolution for different roll states, quantified by the roll aspect ratio, $\varGamma _r \in [1.6, 8]$, with $Ra$ and $Pr$ systematically varied in the ranges $Ra \in [10^{7}, 10^{9}]$ and $Pr \in [1, 100]$.

1. (i) The value of $\gamma$ was slightly increasing with $\varGamma _r$, within the range $\gamma \in [0.302, 0.321]$.

2. (ii) The value of $Nu$ was virtually independent of $Pr$.

3. (iii) The value of $\beta$ was slightly increasing with $\varGamma _r$, within the range $\beta \in [0.657, 0.675]$.

4. (iv) The value of $\sigma$ was slightly increasing with $\varGamma _r$, within the range $\sigma \in [-1.078, -1.043]$.

In conclusion, the data of Wang et al. (2020a) strongly suggest that the 2-D system approaches classical Nusselt-number scaling and Reynolds-number scaling (3.4) at $Ra \sim 10^{9}$, as also pointed out by Wen et al. (2020).

4.2. No-slip boundary conditions

1. (i) Johnston & Doering (2009) report $\gamma = 0.285$ at $Pr = 1$ and $Ra \in [10^{7}, 10^{10}]$. Zhang et al. (2017) report $\gamma = 0.3$ for $Pr = 0.7$ and $5.3$ at $Ra \in [10^{6}, 10^{10}]$. Wang et al. (2020b) plot $Nu/Ra^{1/3}$ at $Pr = 10$ in a lin–log plot and find that $\gamma$ is slightly increasing from $\gamma = 0.262$ in $Ra \in [10^{7}, 10^{8}]$ to $\gamma = 0.289$ at $Ra \in [10^9, 10^{10}]$. van der Poel et al. (2013) plot $Nu/Ra^{1/3}$ at $Pr = 4.38$ in a lin–log plot showing a slightly convex curve at $Ra \in [10^{7}, 10^{10} ]$. If extrapolated to higher $Ra$, the curve would approach a straight line at $Ra \sim 10^{14}$. Zhu et al. (2018) claim that they have determined $Nu$ up to $Ra= 10^{14}$ and Zhu et al. (2019) that they have even reached $Ra = 4.64 \times 10^{14}$. They find that $\gamma = 0.357$ at $Ra > 10^{13}$, which they interpret as evidence of a transition to the ‘ultimate regime’ of heat transfer predicted by Kraichnan (1962) and Grossmann & Lohse (2011). However, all the data of Zhu et al. (2018) and Zhu et al. (2019) for $Ra > 10^{10}$ were evaluated in states that were very far from stationarity (private communication with D. Lohse and X. Zhu). Lohse & Zhu claim (private communication) that the Nusselt number is converged, although the mean kinetic energy is very far from having reached stationarity. As evidence they point out that in a simulations at $Ra = 10^{11}$ that was ran to stationarity after publication of Zhu et al. (2018), $Nu$ increased only by three per cent during the last phase of the simulation. However, in the author's opinion, there is no other way to test whether $Nu$ is converged to the accuracy that is needed in order to distinguish $\gamma = 1/3$ from $\gamma = 0.357$ at $Ra > 10^{13}$ than running the simulations at $Ra > 10^{13}$ for a considerably longer time.

2. (ii) van der Poel et al. (2013) present a plot showing that $Nu$ is virtually independent of $Pr$ at $Ra = 10^{8}$ and $Pr \in [0.3, 100]$. The two curves of $Nu$ by Zhang et al. (2017) for $Pr = 0.7$ and $5.3$ are indistinguishable.

3. (iii) Zhang et al. (2017) report $\beta = 0.6$ for $Pr = 0.7$ and $5.3$. Wang et al. (2020b) observed a slight increase of $\beta$ at $Ra \in [10^{7}, 10^{10} ]$ and $Pr = 10$ with $\beta = 0.565$ for $Ra \in [10^{7}, 10^8]$ and $\beta = 0.595$ for $Ra \in [10^{9}, 10^{10}]$.

4. (iv) van der Poel et al. (2013) present a log–log plot of $Re/Pr^{3/4}$ at $Pr \in [0.1, 60]$. From the slope of the curve it can be estimated that $\sigma \approx -0.9$, which is not too far from $-1$. From the Nusselt-number plot given by Zhang et al. (2017) at $Pr = 0.7$ and $Pr = 5.3$ it can also be estimated that $\sigma$ is not very far from $-1$.

In conclusion, we find it likely but not certain that $Nu$ will approach $Ra^{1/3}$, independent of $Pr$, in the limit of high $Ra$, also in the case with no-slip boundary conditions. Despite the difference between $\beta = 0.6$ and $\beta = 2/3$ it also seems likely that $Re$ will approach $Pr^{-1} Ra^{2/3}$. It can be noted from (3.4) that $\beta$ is not expected to reach $2/3$ until $\gamma$ reaches $1/3$. In the limit of high $Ra$ we should have $\beta = (1+\gamma )/2$. Thus, the observed deviation of $\beta$ from $2/3$ can partly be explained by the observed deviation of $\gamma$ from $1/3$. The observation by Wang et al. (2020b) of a slight increase of $\beta$ with $Ra$ also supports the hypothesis that $\beta$ will approach $2/3$.

4.3. Convergence time scale

The prediction (3.5) of the slow convergence is supported by data communicated to the author by X. Zhu from simulations at $Ra \in [10^{10}, 10^{11} ]$ and $Pr =1$. The simulation is the same as reported by Zhu et al. (2018) that were continued after publication to investigate the convergence of $\tilde {E}$. For the convergence of $\tilde {E}$ we refer to figure 1, which has been prepared using data communicated to the author by X. Zhu. For $Ra = 10^{10}$ and $Ra = 10^{11}$, $\tilde {E}$ reaches approximate stationarity at $\tilde {\tau } \approx 1000$ and $\tilde {\tau } \approx 3000$, with stationary values $\tilde {E} \approx 0.25$ and $\tilde {E} \approx 0.48 \approx 0.5$, in each case respectively, which roughly means that $\tilde {\tau }$ tripled and $\tilde {E}$ doubled when $Ra$ was increased by a factor of ten. These numbers suggest that $\tilde {\tau } \sim Ra^{0.47}$ and $\tilde {E} \sim Ra^{0.30}$. It is clear that the data are generally consistent with (3.5) and (3.3) even if it is impossible to make any detailed comparison. It should be noted that the original simulation at $Ra = 10^{11}$ by Zhu et al. (2018) was far from stationarity when it was ended at $\tilde {t} = 1000$, and the higher $Ra$ simulations were even farther away when they were ended. Assuming that $\tilde {\tau }$ continues to triple and $\tilde {E}$ continues to double when $Ra$ is increased by a factor of ten, the simulation at $Ra = 10^{14}$ would reach stationary first at $\tilde {\tau } \approx 80\,000$ with $\tilde {E} \approx 4$. Since this simulation was ended at $\tilde {t} = 250$ with $\tilde {E} \approx 0. 2$, the Nusselt number was, indeed, calculated in a state very far from stationarity. How close $Nu$ at $\tilde {t}= 250$ was to its stationary value is not known.

Figure 1. Non-dimensional kinetic energy, $\tilde {E}$, vs non-dimensional time, $\tilde {t}$, at $Ra \in [10^{10}, 10^{11}]$, from four simulations that were continued after publication of Zhu et al. (2018). The data from the simulations were communicated to the author by X. Zhu.

5. Conclusion

Assuming that classical Nusselt-number scaling holds and that dissipation, in the absence of a downscale energy cascade, scales as (3.4) in the 2-D Rayleigh–Bénard system, we have shown that the Reynolds-number scales as $Re \sim Pr^{-1} Ra^{2/3}$. This is a much stronger scaling with $Ra$ than the corresponding 3-D relation (1.4). If it is assumed that classical scaling holds for stress-free as well as no-slip boundary condition in two and three dimensions, the Reynolds-number scaling will, most likely, also be equal with stress-free and no-slip conditions, in two dimensions as well as three dimensions.

Using the scaling $\tilde {E} \sim Pr^{-1} Nu$ and the equations of motion we deduced the convergence time scale $\tilde {\tau } \gtrsim Pr^{-1/2}Ra^{1/2}$, without making any assumption regarding the Nusselt-number scaling. The slow convergence is confirmed by data communicated to the author by X. Zhu. From a computational point of view the slow convergence is, of course, disappointing. The general motivation for carrying out 2-D simulations is that they are supposed to reach the same $Ra$ as 3-D simulations, but at a lower computational cost. If this is not true, there is a risk that very few 2-D DNS will be carried out at $Ra > 10^{10}$ in the future, which would be a pity. It is quite possible that strategies can be developed by which the convergence can be substantially improved. Such strategies would include carefully designed initial conditions as well as running the simulations at low resolution until a stationary state is reached after which the resolution is increased. It is the hope of the author that we will see such a development in the near future, so that the predictions of the present paper can be tested in simulations up to at least $Ra = 10^{14}$.