2. Numerical methodology

We perform DNS in a 2-D fluid layer with horizontal dimension $L$ , with an imposed temperature difference $\Delta T$ between the bottom and top plates separated by vertical dimension $H$ . For this work, the aspect ratio is $\Gamma = L/H = 1$ . The following Oberbeck–Boussinesq equations dictating the flow are solved using the Nek5000 solver (which has been used extensively for the simulation of turbulent convection – for details, see Scheel, Emran & Schumacher Reference Scheel, Emran and Schumacher2013):

(2.1) \begin{align}\boldsymbol{\nabla \cdot u} &= 0, \end{align}

(2.2) \begin{align}\frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u \cdot \nabla u} &= - \frac { \boldsymbol{\nabla} p}{\rho _0} + \alpha g (T-T_0) \hat {z} + \nu\nabla^2 \boldsymbol{u}, \end{align}

(2.3) \begin{align}\frac {\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} T &= \kappa\nabla^2 T. \end{align}

Here, $\boldsymbol{u} = (u_x, u_z)$ , $p$ and $T$ are the velocity, pressure and temperature, respectively; $\rho _0$ is the reference density, and $T_0$ is the reference temperature. We non-dimensionalise (2.1)–(2.3) using $H$ , $\Delta T$ , $u_f$ and $t_f$ as the scales for length, temperature, velocity and time, respectively, where $u_f = \sqrt {\alpha g\, \Delta \textit{TH}}$ is the free-fall velocity, and $t_f = H/u_f$ is the free-fall time. The result contains the Prandtl number $Pr$ and the Rayleigh number $Ra = \alpha g\, \Delta \textit{TH}^3/(\nu \kappa )$ , where $\alpha$ is the isobaric thermal expansion coefficient, and $g$ is the acceleration due to gravity. We explore two fluids with $Pr = 0.1$ and $1$ for $Ra$ between $10^6$ and $10^{12}$ . The square domain is decomposed into $N_e$ elements, and each element is further resolved using Lagrangian interpolation polynomials of order $N$ in both horizontal and vertical directions. Thus the entire flow is resolved using $N_e N^2$ mesh cells. The no-slip condition for the velocity field is imposed on all boundaries. Isothermal and adiabatic conditions for the temperature field are imposed on the horizontal plates and sidewalls, respectively. The flow consists of a single large-scale structure covering the entire domain, except for $Pr = 1,\ Ra \leqslant 10^7$ , where two vertically stacked structures occur. The flow morphology is similar to that observed by Labarre, Fauve & Chibbaro (Reference Labarre, Fauve and Chibbaro2023) and Samuel & Verma (Reference Samuel and Verma2024), so we do not show it here.

To resolve the thermal and viscous boundary layers, finer mesh is used near all boundaries. The spatial resolution in the flow is further ensured by computing the Kolmogorov length scale $\eta$ , which is nominally the finest scale in the velocity field, and stipulating that the local vertical grid spacing is $\varDelta_{z}(z)/\eta (z) \lt 1.5$ for all simulations. The kinetic energy and the thermal dissipation rates are defined, respectively, as

(2.4) \begin{eqnarray} \varepsilon _u (\boldsymbol{x}) &= \frac {\nu }{2} \sum _{l,m} \left ( \frac {\partial u_l}{\partial x_m} + \frac {\partial u_m}{\partial x_l} \right )^2 , \end{eqnarray}

(2.5) \begin{eqnarray} \varepsilon _T (\boldsymbol{x}) &=\kappa \left ( \frac {\partial T}{\partial x_l} \right )^2 , \end{eqnarray}

and represent the rates of loss of kinetic energy and thermal energy per unit mass. Here, $l, m = (x,z)$ , and the local Kolmogorov scale is given by $\eta = (\nu ^3/\varepsilon _u)^{1/4}$ . As the intermittent variation of $\varepsilon _u (\boldsymbol{x})$ in the flow leads to variations in the Kolmogorov scale as well, we estimate an average Kolmogorov scale in each horizontal plane using the horizontally and temporally averaged dissipation as

(2.6) \begin{equation} \eta (z) = \frac {\nu ^{3/4}}{\langle \varepsilon _u \rangle _{x,t}^{1/4}(z)}. \end{equation}

The finest length scale in the temperature field, the Batchelor scale $\eta /\sqrt {Pr}$ , is of the same order as $\eta$ , or coarser, in the present work. Thus it is always adequately resolved.

3. Other associated definitions

The Nusselt number in a horizontal plane is computed as (Chillà & Schumacher Reference Chillà and Schumacher2012)

(3.1) \begin{equation} Nu(z) = \frac { \langle u_z T \rangle _{x,t} - \kappa\partial \langle T \rangle _{x,t}/\partial z} { \kappa\kern-0.5pt\Delta T /H}, \end{equation}

where $\langle u_z T \rangle _{x,t}$ is the convective component of the heat flux, and $- \kappa\partial \langle T \rangle _{x,t}/\partial z$ is the diffusive component. As the vertical velocity vanishes at the horizontal plates, this relation yields the definition using the wall temperature gradient as

(3.2) \begin{equation} Nu_W = - \frac {H}{\Delta T} \left . \frac {\partial \langle T \rangle _{x,t}} {\partial z} \right |_{z=0,H} . \end{equation}

Averaging (3.1) along the vertical direction yields the following relation for the global heat transport across the convective layer:

(3.3) \begin{equation} Nu = 1 + \frac {H}{\kappa\kern-0.5pt\Delta T}\, \langle u_z T \rangle _{A,t} . \end{equation}

Here, $\langle \cdot \rangle _{A,t}$ stands for the average over the entire flow domain and simulation time covering the steady state.

In another family of relations, the exact relations in RBC (Howard Reference Howard1972; Shraiman & Siggia Reference Shraiman and Siggia1990) connect the Nusselt number $Nu$ with the globally averaged kinetic energy dissipation rate $\varepsilon _u$ and the thermal dissipation rate $\varepsilon _T$ , as

(3.4) \begin{align}\langle \varepsilon _u \rangle _{A,t} &= \frac {\nu ^3}{H^4}\, \frac {(Nu-1)Ra}{Pr^2} , \end{align}

(3.5) \begin{align}\langle \varepsilon _T \rangle _{A,t} &= \kappa\, \frac {(\Delta T)^2}{H^2}\, Nu. \end{align}

Thus the Nusselt number can also be estimated from mean dissipation rates:

(3.6) \begin{align}Nu_{\varepsilon _u} &= 1 + \frac {H^4}{\nu ^3} \frac {Pr^2}{Ra} \langle \varepsilon _u \rangle _{A,t} , \end{align}

(3.7) \begin{align}Nu_{\varepsilon _T} &= \frac {H^2}{(\Delta T)^2} \frac {\langle \varepsilon _T \rangle _{A,t}}{\kappa }. \end{align}

The agreement among simulated values from different definitions serves as a check on the accuracy and adequacy of spatial and temporal resolutions (Stevens, Verzicco & Lohse Reference Stevens, Verzicco and Lohse2010; Zhang et al. Reference Zhang, Zhou and Sun2017; Pandey et al. Reference Pandey, Krasnov, Schumacher, Samtaney and Sreenivasan2022a ). In table 1, we list the global heat flux estimated from all the methods, and note that the agreement among various estimates of the Nusselt number is indeed very good. Also listed are other crucial parameters of the flow. We discuss below the approach to the final state of the Nusselt number using these definitions and relations, with comments on the scaling of global averages. In (3.1)–(3.3) and (3.6)–(3.7), the Nusselt numbers are globally averaged quantities, though we do not show the averaging symbol explicitly, following standard usage of the past.

Table 1.

Important parameters of DNS in a 2-D box with $\Gamma = 1$ . We list the Prandtl number, the Rayleigh number, the total number of mesh cells in the entire flow domain $N_e N^2$ , the Nusselt numbers using (3.3), (3.6), (3.7) and (3.2), respectively, the Reynolds number, and the simulation time after the flow attains a steady state, $t_{sim}$ . Error bars in Nusselt and Reynolds numbers are the corresponding standard deviations.

In the following, we study the evolutions of the instantaneous domain-averaged heat fluxes, which are defined as

(3.8) \begin{align}\overline {Nu} &= 1 + \frac {H}{\kappa\kern-0.5pt\Delta T}\, \langle u_z T \rangle _{A} , \\[-10pt]\nonumber\end{align}

(3.9) \begin{align} \overline {Nu_{\varepsilon _u}} &= 1 + \frac {H^4}{\nu ^3}\, \frac {Pr^2}{Ra}\, \langle \varepsilon _u \rangle _{A} , \end{align}

(3.10) \begin{align} \overline {Nu_{\varepsilon _T}} = \frac {H^2}{(\Delta T)^2}\, \frac {\langle \varepsilon _T \rangle _{A}}{\kappa }. \end{align}

4. Transient characteristics

A common method for initiating high- $Ra$ simulations is to start them from the flow at a lower $Ra$ value. Simulations can also be performed ab initio from the conduction state with random perturbations. In both cases, global heat flux and kinetic energy evolve with time, and one needs to wait some time before a statistically steady state is attained. Figure 1 shows the temporal evolution of integral quantities for the two $Ra$ indicated and $Pr = 0.1$ , with simulations initiated from the conduction state. On the one hand, we observe that $\overline {Nu}$ and $\overline {Nu_{\varepsilon _T}}$ in figure 1(b,d) start oscillating about some mean value shortly after the simulation begins. On the other hand, $\overline {Nu_{\varepsilon _u}}$ in figure 1(c) initially increases and starts to fluctuate about a mean value only after some time has elapsed. The instantaneous domain-averaged non-dimensional kinetic energy $E = \langle (u_x^2 + u_z^2)/2 \rangle _A/u_f^2$ in figure 1(a) offers the best means to determine the time to the steady state. There is no ambiguity about this approach for the lower $Ra$ ; however, aside from the fact that it takes longer at the higher $Ra$ for the energy to achieve its steady state, the latter is somewhat nominal because fluctuations about the average are significant. The situation becomes more so at even higher Rayleigh numbers. Fluctuations in the domain-averaged energy $E$ follow from dynamical considerations; we will show this in § 6.

Figure 1.

Evolution of the integral quantities in the transient state for $Pr = 0.1$ , and $Ra = 3 \times 10^8$ (black curves) and $Ra = 3 \times 10^9$ (orange curves). (a) The domain-averaged kinetic energy $E$ increases slowly and takes a few thousand free-fall times to reach the steady state. The transient time is longer for higher $Ra$ . Dashed vertical lines show a quantitative measure of the transient time $t_{\it trns}$ , obtained from (4.3), to be discussed later. (b,c,d) Plots show that the Nusselt number fluctuates rapidly about its mean nearly from the start, but fluctuations have different characters depending on the definition of the Nusselt number. In (b), the fluctuations in $\overline {Nu}$ are strong and of high frequency with no well-defined transient state, and there is an overlap for the two $Ra$ values. In (c), one can approximately identify a transient state in $\overline {Nu_{\varepsilon _u}}$ , which exhibits very strong fluctuations containing both high and low frequencies. Plot (d) shows that $\overline {Nu_{\varepsilon _T}}$ fluctuates similarly to $\overline {Nu}$ in (b), but there is no overlap for the two Rayleigh numbers. The occasional appearance of negative values of $\overline {Nu}$ suggests the likelihood that a small parcel of the coldest fluid from the top wall registers directly at the bottom wall, and vice versa.

We define the average kinetic energy in the steady state, $E_{av}$ , as

(4.1) \begin{equation} E_{av} = \tfrac {1}{2}u_{\it RMS}^2/u_f^2 , \end{equation}

where $u_{\it RMS} = \sqrt { \langle u_x^2 + u_z^2 \rangle _{A,t}}$ is the root mean square (RMS) velocity of the flow in the steady state, and plot it in figure 2 as a function of $Ra$ . Before discussing differences between the two Prandtl numbers, we note the major difference between 2-D and 3-D fluctuations. In the 3-D case, for low and moderate Prandtl numbers, $E_{av}$ is a slowly decreasing function of $Ra$ ; for example, $E_{av}$ from a horizontally periodic cuboid of $\Gamma = 4$ for $Pr = 0.7$ (from Samuel et al. Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024), shown as green squares in figure 2, follows a $Ra^{-0.07}$ scaling. The corresponding behaviour for 2-D cases is that the fluctuations increase with $Ra$ . The behaviour is different for different $Pr$ when the Rayleigh numbers are low, but follows an approximately $Ra^{1/3}$ scaling at high Rayleigh numbers, as indicated by the blue and red dashed lines in figure 2. (Our objective is not a detailed study of differences between 2-D and 3-D convection. Some such differences have been pointed out e.g. by van der Poel et al. (Reference van der Poel, Stevens and Lohse2013) and Pandey et al. (Reference Pandey, Kumar, Chatterjee and Verma2016)). The challenges posed by this increasing trend will be discussed below.

Figure 2.

Global (area- as well as time-averaged) kinetic energy in the steady state $E_{av}$ as a function of $Ra$ . An increasing trend with $Ra$ is observed in 2-D RBC. In the turbulent regimes ( $Ra \geqslant 10^8$ for $Pr = 0.1$ , and $Ra \gt 10^{9}$ for $Pr = 1$ ), the data approximately follow $E_{av} \sim Ra^{1/3}$ , shown as dashed lines. In contrast, $E_{av}$ in 3-D RBC for $Pr = 0.7$ (taken from Samuel et al. Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024) shows a weakly decreasing trend. The transitions observed in this figure are related to transitions in the flow structure; see § 5.3. The statistical error bars in almost all the figures here are comparable to the thickness of the symbols. The exception is figure 6, for which the errors bars are shown explicitly.

It is useful to study the variation of $t_{\it trns}$ (which is the time required for $E$ to attain some stipulated fraction of $E_{av}$ ) with respect to $Ra$ and $Pr$ by some simple scheme. To this end, we plot the evolution of $E$ for $Ra = 3 \times 10^9$ and $10^{10}$ for $Pr = 0.1$ in figure 3. The growth of $E$ for both $Ra$ values is comparable for initial times but the curve for higher $Ra$ continues to grow for a longer period before a nominally steady state is reached, but with conspicuous fluctuations.

Figure 3.

Evolution of the domain-averaged kinetic energy $E(t)$ in the transient state for $Pr = 0.1$ and $Ra = 3 \times 10^9$ and $10^{10}$ can be described well by (4.2), and the dashed curves show $E_{\it fit}(t)$ for the growth rate. Vertical lines show that the transient time $t_{\it trns}$ is much longer for $Ra = 10^{10}$ than for $Ra = 3 \times 10^9$ .

Figure 3 suggests that $E(t)$ can be fitted with an equation of the form

(4.2) \begin{equation} [E_{av} - E(t)]/E_{av} = c \exp (-kt), \end{equation}

where $E_{av}$ is the ‘steady-state’ or ‘asymptotic’ mean energy of the flow; the growth rate $k$ depends on the Rayleigh and Prandtl numbers. The factor $c$ captures the finite energy at the initial instant in figure 3: though $E(t)$ grows rapidly only when the convective motion is established, there is a finite $E(t)$ from which it starts.

We fit the suitable segment of the growth curve $E(t)$ with (4.2). The segment at late times is not suitable for fitting the formula because the kinetic energy does not attain a constant value but fluctuates strongly for long periods of time. The scaling factor $c$ is calculated from the data as $c = (E_{av} - E(0))/E_{av}$ , where $E(0)$ is the energy at $t = 0$ . The dashed curves in figure 3 are fits to the data. Having determined the growth rate, we define the transient time $t_{\it trns}$ as the time when the fitted curve $E_{\it fit}(t)$ in figure 3 reaches 96 % of its ‘constant’ value, $E_{av}$ . Thus the transient time is estimated as

(4.3) \begin{equation} t_{\it trns} = \frac {1}{k} \log \frac {c}{0.04} . \end{equation}

The transient time thus estimated is plotted as a function of $Ra$ in figure 4(a). Note that $t_{\it trns}$ increases as a power law; for $Pr = 0.1$ , the best fit yields $t_{\it trns} \sim Ra^{0.59 \,\pm\, 0.03}$ , while we find $t_{\it trns} \sim Ra^{0.71\, \pm\, 0.08}$ for $Pr = 1$ .

We have experimented with different definitions of $t_{\it trns}$ (e.g. by requiring it to reach 90 % of $E_{av}$ ), and find the same trends, though the numbers are different. The fact that the $Ra$ exponent is non-trivial suggests an important role for the boundary layers and their relation to the large structure; those details (as well as the effect of the aspect ratio) are in need of further exploration.

Figure 4.

(a) The transient time $t_{\it trns}$ as a function of $Ra$ shows a power law. For a given $Ra$ , $t_{\it trns}$ is longer for lower $Pr$ . (b) The time $t_{\it trns}$ as a function of ${Re}$ is essentially the same for both values of $Pr$ , and exhibits a nearly linear trend. Moreover, $t_{\it trns}$ is nearly the same for both $Pr$ values when the Reynolds numbers are the same. Data correspond only to the turbulent regimes. Sparser data sets for $Pr = 0.021$ are consistent with the Prandtl numbers of these plots.

Figure 5.

Inverse growth rate as a function of (a) $Ra$ and (b) ${Re}$ . The trends with ${Re}$ are approximately the same for both Prandtl numbers, and are slightly below linear.

Figure 4(a) further shows that for a given $Ra$ , $t_{\it trns}$ is longer for $Pr = 0.1$ than for $Pr = 1$ . This is expected as $u_{\it RMS}$ – consequently the Reynolds number – is higher for smaller $Pr$ (Schumacher, Götzfried & Scheel Reference Schumacher, Götzfried and Scheel2015; Pandey & Verma Reference Pandey and Verma2016; Pandey et al. Reference Pandey, Krasnov, Schumacher, Samtaney and Sreenivasan2022a ). In figure 4(b), we plot $t_{\it trns}$ as a function of the Reynolds number ${Re}$ , which is defined as

(4.4) \begin{equation} {Re} = u_{\it RMS}H/\nu . \end{equation}

We find that the two $t_{\it trns}$ fall approximately on the same line for both Prandtl numbers, the best fit for which is nearly linear. Thus a 2-D RBC at high $Ra$ has to be simulated for very long times to achieve a statistically steady state. If a steady state is not achieved and $E$ is still growing, then one will find $Nu_{\varepsilon _u} \lt Nu$ (see § 6), even when spatial and temporal resolutions are adequate.

Shown in figure 5 is the variation of the inverse growth rate $k^{-1}$ with respect to $Ra$ (figure 5 a) and ${Re}$ (figure 5 b). It appears that $k^{-1}$ is approximately a linear function of ${Re}$ , and depends only weakly on the Prandtl number.

A longer transient at higher $Ra$ compels researchers to use some workaround to achieve the statistically steady state in a shorter time. For example, one might start with a much coarser resolution and ramp it up to the required level only after the kinetic energy has reached an approximate steady state. It is practically impossible with available computing power to conduct very high $Ra$ simulation with fully-resolved fields in the entire transient state. One can use a coarser mesh during the transient when all degrees of freedom have not yet been excited. However, it is important to ensure that simulations on the coarser mesh lead to the same steady state as that attained in a well-resolved simulation. A concern otherwise would be that the coarse-grid simulation results in kinetic energy different from the properly resolved case, ultimately affecting the scaling of ${Re}$ , particularly because of the presence of the low frequency modes (see § 6). Though we demonstrate in Appendix A that the mean kinetic energy in the steady state and the transient time are essentially the same in coarser simulations, it is not necessarily likely to be a general result. The data in figures 1 and 3 correspond to simulations at the coarser resolution.

5. Scaling exponents for integral transport in the steady state

5.1. Heat transport

We plot $Nu$ , computed from (3.2), as a function of $Ra$ in figure 6(a). The data for high Rayleigh numbers seem to closely follow similar power laws for both $Pr$ values. For $Pr = 1$ , the best fit for $Ra \geqslant 6 \times 10^7$ yields the scaling $Nu = 0.12\, Ra^{0.29\,\pm\,0.003}$ . The exponent is close to $2/7$ proposed for the so-called ‘hard turbulence’ in confined 3-D RBC (Castaing et al. Reference Castaing, Gunaratne, Kadanoff, Libchaber and Heslot1989), also explored in various later studies (Siggia Reference Siggia1994; Chillà & Schumacher Reference Chillà and Schumacher2012; Johnston & Doering Reference Johnston and Doering2009). We plot the normalised Nusselt number $Nu\, Ra^{-2/7}$ versus $Ra$ in figure 6(b), which reveals that the local exponent varies considerably and is only approximately $2/7$ . Furthermore, we note that the Nusselt numbers fluctuate significantly, being significantly larger when $Nu$ is computed using (3.3). To some extent, this result indicates that the local scaling exponents obtained by fitting heat transport data over short ranges of Rayleigh numbers could lead to misleading conclusions. Figure 6(a) also shows that the Nusselt numbers for $Ra \leqslant 10^7$ for $Pr = 1$ do not follow the higher- $Ra$ trend. This is a consequence of differing flow structures observed: the flow for $Ra \leqslant 10^7$ consists of two vertically stacked rolls, whereas for $Ra \gt 10^7$ , a single-roll structure is observed. These findings suggest that the flow is less efficient in transporting heat when the double-roll state occurs, which is in line with previous observations in both two and three dimensions (Xi & Xia Reference Xi and Xia2008; van der Poel, Stevens & Lohse Reference van der Poel, Stevens and Lohse2011; Weiss & Ahlers Reference Weiss and Ahlers2011).

Figure 6.

(a) The Nusselt number as a function of $Ra$ for $Pr = 0.1$ (red circles) and $Pr = 1$ (blue triangles). The scaling for high Rayleigh numbers differs only slightly from the low- $Ra$ behaviour. (b) The normalised Nusselt number $Nu\, Ra^{-2/7}$ shows that the $2/7$ scaling is only approximate for moderate Rayleigh numbers and moderate aspect ratios. Wall heat flux from (3.2) is shown here with the error bars representing the standard deviation.

For $Ra \leqslant 3 \times 10^7$ , the Nusselt numbers for $Pr = 0.1$ are higher than those for $Pr = 1$ . Unlike for $Pr = 1$ , the flow for $Pr = 0.1$ shows no double-roll state at lower Rayleigh numbers, which is why a larger difference appears in $Nu$ for the two Prandtl numbers. Figure 6(b) clearly shows that heat transport is smaller for $Pr = 0.1$ than for $Pr = 1$ in the turbulent regime, i.e. for $Ra \gt 3 \times 10^7$ . This feature is similar to that observed in 3-D RBC, where lower- $Pr$ fluids are less efficient at transporting heat when $Pr \lt 1$ (Verzicco & Camussi Reference Verzicco and Camussi1999; van der Poel et al. Reference van der Poel, Stevens and Lohse2013; Pandey & Sreenivasan Reference Pandey and Sreenivasan2021; Pandey et al. Reference Pandey, Krasnov, Sreenivasan and Schumacher2022b ). The Grossmann–Lohse model (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001) also suggests a similar trend.

Figure 7.

The Reynolds number based on $u_{\it RMS}$ scales nearly as $Ra^{2/3}$ in the turbulent regime (indicated by the blue solid line), which is distinctively different from scaling $Ra^{1/2}$ reported in 3-D RBC. The green dashed line indicates the ${Re} \sim Ra^{0.46}$ scaling observed for 3-D RBC in a $\Gamma = 4$ box by Samuel et al. (Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024).

5.2. Momentum transport

A variety of velocity scales can be defined in turbulent RBC, and can be used to define the Reynolds number ${Re}$ . In cylindrical or cubic domains with $\Gamma \approx 1$ , the most dominant eddy in the flow is in the form of a large-scale circulation, and its velocity is observed to scale with the free-fall velocity (Lam et al. Reference Lam, Shang, Zhou and Xia2002; Xia, Sun & Zhou Reference Xia, Sun and Zhou2003). In DNS, the Reynolds number is often obtained using the RMS velocity and the depth of the convective layer $H$ (see (4.4)) (Scheel & Schumacher Reference Scheel and Schumacher2016; Pandey et al. Reference Pandey, Krasnov, Sreenivasan and Schumacher2022b ). We plot ${Re}\,Pr$ against $Ra$ in figure 7. Also indicated by the dashed line is the power law ${Re}\,Pr = 0.33\, Ra^{0.46}$ , as found by Samuel et al. (Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024) for $Pr = 0.7$ . There is some overlap in the magnitude of the Reynolds numbers between 2-D and 3-D RBC, but the scaling exponents differ, especially for high Rayleigh numbers. In the high- $Ra$ regime, we find that the Reynolds number scales as ${Re} = 0.13\, Ra^{0.65\,\pm\,0.01}$ for $Pr = 0.1$ , and ${Re} = 0.02\, Ra^{0.65\,\pm\,0.01}$ for $Pr = 1$ . The data for lower Rayleigh numbers exhibit approximately the ${Re} \sim Ra^{0.55}$ scaling for both $Pr$ values.

As is well known, the Reynolds number decreases with increasing $Pr$ in 3-D RBC (van der Poel et al. Reference van der Poel, Stevens and Lohse2013; Pandey & Sreenivasan Reference Pandey and Sreenivasan2021); similarly, in 2-D RBC, the Reynolds number scales approximately as $Ra^{2/3}\,Pr^{-1}$ , consistent with the result found semi-analytically by Chini & Cox (Reference Chini and Cox2009) and Wen et al. (Reference Wen, Goluskin, LeDuc, Chini and Doering2020). This ${Re}$ scaling further suggests that

(5.1) \begin{equation} \frac {u_{\it RMS}}{u_f} = {Re}\, \sqrt {\frac {Pr}{Ra}} \sim Ra^{1/6}\, Pr^{-1/2} , \end{equation}

and consequently

(5.2) \begin{equation} E_{av} = \tfrac{1}{2} \frac {u^2_{\it RMS}}{u^2_f} \sim Ra^{1/3}\, Pr^{-1}. \end{equation}

The high- $Ra$ data in figure 2 are indeed consistent with this scaling.

5.3. The RMS temperature fluctuation

We compute the RMS temperature fluctuation as

(5.3) \begin{equation} T_{\it RMS} = \sqrt { \langle T^2 \rangle _{A,t} - \langle T \rangle _{A,t}^2}, \end{equation}

and plot $T_{\it RMS}/\Delta T$ as a function of $Ra$ in figure 8(a). It is clear that $T_{\it RMS}/\Delta T$ decreases with increasing $Ra$ , but various regions can be identified in figure 8(a). For $Pr = 0.1$ , $T_{\it RMS}$ for $Ra \leqslant 2 \times 10^7$ exhibits $Ra^{-0.07}$ scaling, but it starts, somewhat abruptly, to follow a steeper $Ra^{-0.12}$ scaling for higher $Ra$ . For $Pr = 1$ , too, we find that $T_{\it RMS}$ for $Ra \lt 10^9$ shows a scaling $Ra^{-0.11}$ , while a scaling $Ra^{-0.14}$ ensues for $Ra \geqslant 10^9$ . Again, the transition between these two regimes is nearly abrupt. The RMS fluctuations for $Ra \leqslant 10^7$ are larger and clearly depart from the trend for higher $Ra$ ; this occurs because of the double-roll state observed for weak thermal forcing at $Pr = 1$ . It is interesting that the scaling $Ra^{-0.14}$ for high $Ra$ agrees well with scalings of fluctuations at the centre of cylindrical RBC cells for $Pr \approx 0.7$ (Castaing et al. Reference Castaing, Gunaratne, Kadanoff, Libchaber and Heslot1989; Niemela et al. Reference Niemela, Skrbek, Sreenivasan and Donnelly2000), as well as with that in the bulk region of a horizontally periodic box for $Pr = 0.7$ , $\Gamma = 4$ (Samuel et al. Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024); see also Pandey & Verma (Reference Pandey and Verma2016). We also note that the magnitudes of $T_{\it RMS}$ for moderate Rayleigh numbers ( $10^7 \lt Ra \lt 10^9$ ) are quite similar for the two Prandtl numbers.

Figure 8.

(a) The RMS temperature fluctuation $T_{\it RMS}$ decreases with $Ra$ for both $Pr$ values, and has nearly the same magnitude for moderate Rayleigh numbers. Different exponents as well as prefactors are found for moderate and high Rayleigh numbers. (b) Fluctuation $T_{\it RMS}$ as a function of ${Re}$ shows that it scales approximately as ${Re}^{-0.2}$ in the turbulent regime (for ${Re} \gt 10^4$ ) for both Prandtl numbers.

The nearly abrupt transitions in $T_{\it RMS}$ are related to changes in flow morphology and statistics, and in the corresponding heat transport scaling. Data for $Pr = 0.1$ show that the scaling exponent $\gamma$ in the $Nu$ – $Ra$ relation changes from ${\approx}0.35$ to ${\approx}0.25$ for $Ra \geqslant 3 \times 10^7$ , with the latter value persisting for nearly two decades of $Ra$ (i.e. up to $Ra \approx 10^9$ ), beyond which it again increases to 0.30. Similarly, for $Pr = 1$ , $\gamma$ changes from ${\approx}0.31$ to ${\approx}2/7$ after the transition at $Ra \approx 10^9$ . Labarre et al. (Reference Labarre, Fauve and Chibbaro2023) recently noted that the ratio between the RMS fluctuations and the mean heat flux increases abruptly in two dimensions once $Ra/Pr$ increases beyond ${\approx}10^9$ . A somewhat similar transition was reported also by Gao et al. (Reference Gao, Tao, Huang, Bao and Xie2024), who found the transition $Ra$ to scale with the Prandtl number as $Pr^{1.41}$ . Our own data show that the relative fluctuations in the heat flux near the wall and in the bulk are enhanced significantly for $Ra \geqslant 3 \times 10^7$ at $Pr = 0.1$ , and for $Ra \geqslant 10^9$ at $Pr = 1$ (see table 1).

The decrease in $T_{\it RMS}/\Delta T$ with $Ra$ is related to the thermal boundary layer, which becomes thinner with $Ra$ (Scheel, Kim & White Reference Scheel, Kim and White2012; Pandey Reference Pandey2021). As $T_{\it RMS}$ represents an average measure of the temperature anomaly in the flow, the dominant contribution to $T_{\it RMS}$ arises from regions occupied by thermal plumes. This is because the temperature within the plumes varies slowly and differs strongly from the ambient temperature, which is approximately the mean temperature $\Delta T/2$ in the flow. As the fraction of the volume occupied by the plumes decreases with $Ra$ , so does their contribution to $T_{\it RMS}/\Delta T$ . A similar magnitude of RMS fluctuations for the two Prandtl numbers (and moderate Rayleigh numbers) is due to the similar nature in the two cases of heat transport (see figure 6), which determines the thickness of the thermal boundary layer.

In figure 8(b), we show $T_{\it RMS}/\Delta T$ as a function of ${Re}$ . Although the data for the two Prandtl numbers are distinct, the scaling regimes reveal themselves clearly. We observe that for large ${Re}$ , the scaling exponent of $T_{\it RMS}/\Delta T$ with respect to ${Re}$ is essentially the same for both $Pr$ values. For ${Re} \gt 10^4$ , temperature RMS exhibits the same scaling, $T_{\it RMS}/\Delta T \sim {Re}^{-0.2}$ , for both $Pr$ values. However, the exponents in moderate Reynolds numbers, to the extent that they can be defined at all, are different, with $T_{\it RMS}/\Delta T$ showing ${Re}^{-0.13}$ and ${Re}^{-0.21}$ for $Pr = 0.1$ and $Pr = 1$ , respectively. They are unlikely to be of fundamental significance.

6. Fluctuation of global quantities

We now discuss the fluctuation of integral quantities in the nominally steady state. Taking the dot product of (2.2) with $\boldsymbol{u}$ , and averaging over the entire domain, we obtain

(6.1) \begin{equation} \frac {\partial }{\partial t} \big\langle u_i^2 / 2 \big\rangle _A = \alpha g\, \langle u_z T \rangle _A - \langle \varepsilon _u \rangle _A . \end{equation}

Recalling that $E = \langle (u_x^2 + u_z^2)/2 \rangle _A/u_f^2$ is the domain-averaged kinetic energy, (6.1) takes the non-dimensional form

(6.2) \begin{equation} \sqrt {Ra\,Pr}\, \frac {\partial E}{\partial t} = \overline {Nu} - \overline {Nu_{\varepsilon _u}} , \end{equation}

where $\overline {Nu}$ and $\overline {Nu_{\varepsilon _u}}$ are the instantaneous domain-averaged heat fluxes defined in (3.8) and (3.9), respectively. Equation (6.2) states that $\overline {Nu}$ and $\overline {Nu_{\varepsilon _u}}$ are not equal to each other whenever $E$ varies with time.

Figure 9.

Temporal evolution of the integral quantities in statistically steady state for $Pr = 0.1,\ Ra = 10^{10}$ : (a) domain-averaged scaled kinetic energy $\sqrt {Ra\,Pr} \, E$ is dominated by a slow evolution; (b) $\overline {Nu}$ fluctuates rapidly about its mean; (c) $\overline {Nu_{\varepsilon _u}}$ , in addition to having rapidly fluctuating components, evolves slowly and is related to $E$ (see (6.2)); (d) $\overline {Nu_{\varepsilon _T}}$ fluctuates rapidly, but a weak slowly-varying trend is present. The horizontal dashed line in all the plots indicates the time-averaged quantity. Here, the origin is taken to be $3000 \, t_f$ of figure 3.

In figure 9, we show the temporal evolutions of $\sqrt {Ra\,Pr}\, E,\ \overline {Nu},\ \overline {Nu_{\varepsilon _u}},\ \overline {Nu_{\varepsilon _T}}$ in the nominally steady state for $Pr = 0.1,\ Ra = 10^{10}$ . Each quantity evolves differently from the others. Figure 9(a) shows that domain-averaged kinetic energy $E$ contains sizeable changes with respect to its mean value, indicated by a dashed horizontal line, occurring in 40–60 units of free-fall time. In figure 9(c), $\overline {Nu_{\varepsilon _u}}$ shows a slow variation superimposed by strong rapid fluctuations. In contrast, $\overline {Nu}$ in figure 9(b) fluctuates rapidly around its mean value. In figure 9(d), $\overline {Nu_{\varepsilon _T}}$ also fluctuates rapidly about its mean, with weak (but non-monotonic) trends. Data for $Ra \gt 10^8$ at $Pr = 0.1$ exhibit similar characteristics.

Figure 10.

Temporal evolution of the integral quantities in steady state for $Pr = 1,\ Ra = 10^{12}$ . The descriptions are the same as in figure 9.

Coming to the energy balance equation, if we average (6.2) over a finite interval of time, say between an initial (non-dimensional) time $t_{\it init}$ and a final time $t_{\it fin}$ , then we obtain

(6.3) \begin{equation} \sqrt {Ra\,Pr}\, \frac {\Delta E}{\Delta t} = \langle Nu \rangle _{\Delta t} - \langle Nu_{\varepsilon _u} \rangle _{\Delta t}, \end{equation}

where $\Delta E = E(t_{\it fin}) - E(t_{\it init})$ , $\Delta t = t_{\it fin} - t_{\it init}$ , and $\langle \cdot \rangle _{\Delta t}$ denotes an averaging over the time interval $\Delta t$ . As there exist long periods of growth or decay of $E$ (see figure 9 a), we can apply (6.3) to those intervals. For example, focusing on a segment where $E$ grows in figure 9 (the region highlighted by red shading), we find the left-hand side of (6.3) to be $\sqrt {Ra\,Pr}\, \Delta E/\Delta t = 23.8$ . During this period, the average values of the heat fluxes are found to be $\langle Nu \rangle _{\Delta t} = 80.3$ and $\langle Nu_{\varepsilon _u} \rangle _{\Delta t} = 56.4$ , yielding 23.9 for the right-hand side. Thus (6.3) is satisfied perfectly. Similarly, in the blue-shaded region in figure 9 where $E$ decays, we obtain $\sqrt {Ra\,Pr}\, \Delta E/\Delta t = -33,\ \langle Nu \rangle _{\Delta t} = 90.4$ and $\langle Nu_{\varepsilon _u} \rangle _{\Delta t} = 123.5$ ; thus the terms of (6.3) balance perfectly again. We note that the power spectra of the quantities shown in figure 9 show slow oscillations, or low-frequency modes, in $E$ and $Nu_{\varepsilon _u}$ , $Nu$ and $Nu_{\varepsilon _T}$ , although these modes diminish in power compared to the significantly strengthening high frequencies.

On the one hand, due to the rapid fluctuation of $\overline {Nu}$ about its mean, its short-term average does not differ much from the long-term average. For example, $\langle Nu \rangle _{\Delta t} = 80.3$ and $\langle Nu \rangle _{\Delta t} = 90.4$ in the same two intervals are not far from the average 84.7. On the other hand, the presence of a strong low-frequency component in $\overline {Nu_{\varepsilon _u}}$ causes short-term averages to differ significantly from the long-term value. For example, $\langle Nu_{\varepsilon _u} \rangle _{\Delta t} = 56.4$ and $\langle Nu_{\varepsilon _u} \rangle _{\Delta t} = 123.5$ in the growing and decaying periods of $E$ , respectively, differ by up to 50–60 % from $Nu_{\varepsilon _u} = 79.2$ . This applies to all the high- ${Re}$ data in 2-D RBC that we have explored. For example, figure 10 demonstrates the same picture for $Pr = 1,\ Ra = 10^{12}$ , where in the red- and blue-shaded regions, $\langle Nu \rangle _{\Delta t}$ is, respectively, larger and smaller than $\langle Nu_{\varepsilon _u} \rangle _{\Delta t}$ , and (6.3) applies perfectly.

As $Ra$ approaches very high values, the overwhelming resources required to simulate convective flows in two dimensions limit the total simulation time available to gather statistics. As a result, $Nu$ and $Nu_{\varepsilon _u}$ may not converge perfectly even if the sufficiency of spatial and temporal resolutions is ensured. For the simulation at $Pr = 0.1,\ Ra = 10^{10}$ (shown in figure 9), $Nu$ and $Nu_{\varepsilon _u}$ differ by more than 6 % (see table 1). Similarly, for $Pr = 1,\ Ra = 10^{12}$ (shown in figure 10), the difference is also approximately 5 %. For lower $Ra$ , on the other hand, there is better convergence, to within 1–2 %. The convergence of $Nu_{\varepsilon _T}$ and $Nu$ is much better because both $\overline {Nu_{\varepsilon _T}}$ and $\overline {Nu}$ oscillate with comparable rapidity about their long-term averages.

7. Summary and conclusions

Our goal here has been to study the nature of the transient evolution of the DNS of 2-D thermal convection, using the no-slip boundary condition on all the walls, along with isothermal bottom and top walls, and adiabatic sidewalls. We illustrate related features using a square box for Prandtl numbers 0.1 and 1, in the Rayleigh number range between $10^6$ and $10^{12}$ . We particularly study the temporal evolution of integral transport quantities – such as the Nusselt number (defined in three different ways) and the turbulent energy – and discuss their scaling. The nominally steady state is reached exponentially with substantial dependence on Rayleigh and Prandtl numbers. Although there is some degree of common behaviour of transients for all the conditions explored here, there is no strict universality to the details of the exponential approach. We find, perhaps not surprisingly, that the velocity field evolves more slowly than the thermal field. That the velocity field evolves more slowly is not surprising because this evolution involves a few intermediate steps from the onset of heating changes. In the relation ${Re} \sim Ra^{\zeta }$ , if the exponent $\zeta \gt 0.5$ as in 2-D RBC, then the ratio $u_{\it RMS}/u_f$ continues to increase with $Ra$ . It is evident from (5.2) (see also (4.2) and (4.3)) that $t_{\it trns}$ increases with increasing fluctuating energy.

We also call attention to large oscillations of the velocity field in what may be regarded effectively as the steady state. One main conclusion is that these low-frequency oscillations are related to differences between the Nusselt number defined by the correlation of $u_z$ and $T$ , and the Nusselt number based on the energy dissipation (see (6.2)). The time to saturation of the turbulent energy is presumably dependent on the formation of the large structure (Smith & Yakhot Reference Smith and Yakhot1993), which itself would depend on the aspect ratio. The relation between the formation of the large structure and the time to saturation remains unclear at present, but it appears that achieving the so-called ultimate state of convection for smooth boundaries is as elusive in two dimensions as in three. The naive expectation that the reduced dimensionality in two dimensions allows one to settle the question of the ultimate state by attaining high enough $Ra$ is thwarted to some degree by the appearance of strong, low-frequency fluctuations that demand excessively long simulation times for definitive scaling relations.