1. Introduction
An outstanding question in thermal convection is the dynamical interaction of the flow close to the boundaries with the rest of the flow, leading to the emission of plumes. For large aspect ratios, thermal boundary layers can be identified even at very high Rayleigh numbers, whereas a unidirectional momentum boundary layer does not seem to exist (Samuel et al. Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024; Shevkar et al. Reference Shevkar, Samuel, Zinchenko, Bode, Schumacher and Sreenivasan2025; Samuel & Schumacher Reference Samuel and Schumacher2025). For small aspect ratios, the emergence of a strong mean wind (Sreenivasan, Bershadskii & Niemela Reference Sreenivasan, Bershadskii and Niemela2002; Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; Shi, Emran & Schumacher Reference Shi, Emran and Schumacher2012; Lohse & Shishkina Reference Lohse and Shishkina2024) maintains a semblance to the momentum boundary layer – although, because of the small ‘fetch’ and the up- and down-welling regions near the bounding vertical walls, the flow bears little resemblance to the canonical form. Moreover, the direction of the mean wind drifts in cylindrical cells (Sreenivasan et al. Reference Sreenivasan, Bershadskii and Niemela2002; Xi, Zhou & Xia Reference Xi, Zhou and Xia2006; Mishra et al. Reference Mishra, De, Verma and Eswaran2011; Schumacher et al. Reference Schumacher, Bandaru, Pandey and Scheel2016) as well as horizontally periodic (van Reeuwijk, Jonker & Hanjalić Reference van Reeuwijk, Jonker and Hanjalić2008) domains.
However, if the flow has a longer fetch in a horizontal direction, the near-wall flow could be closer to canonical wall-bounded flows, at least in the mid-region. If a flow configuration with a well-defined mean flow can be identified, it would be worthwhile understanding its dynamics. To this end, Rayleigh–Bénard convection (RBC) in a rectangular cuboid of length
$L_x=2.4H$
, width
$L_y=0.8H$
and height
$H$
is studied here, in the configuration shown in figure 1. As identified schematically in figure 2, the flow is characterised by ejection regions in which the thermal plumes emanate, moving essentially upwards into the bulk of the flow, impact regions that receive vertically moving thermal plumes traversing from the opposite plate and shear regions in which the plumes after impact travel mostly parallel to the plate.
It is already known that the boundary layer properties in the ejection, shear and impact regions are different (He & Xia Reference He and Xia2019). Two-dimensional (2-D) models of RBC also reveal that the properties of thermal boundary layer (TBL) in the ejection region differ substantially from those in the shear and impact regions (van der Poel et al. Reference van der Poel, Ostilla-Mónico, Verzicco, Grossmann and Lohse2015), with temperature profiles in the ejection region exhibiting a logarithmic behaviour, which is absent in the shear and impact regions (Pandey Reference Pandey2021; He, Bao & Chen Reference He, Bao and Chen2024). The ejection and impact regions at the horizontal plates affect flow properties within the boundary layers as well as those deeper into the bulk (He, Tong & Xia Reference He, Tong and Xia2007; Chong & Xia Reference Chong and Xia2016; Huang & Xia Reference Huang and Xia2016; Xu, Zhang & Xia Reference Xu, Zhang and Xia2024).
Instantaneous flow realisation for
$\textit{Ra} = 10^{10}$
in the closed rectangular cuboid with the red and blue colours representing, respectively, hot upwelling and cold downwelling structures. Fine thermal structures are revealed by temperature isosurfaces: hot and cold structures correspond to
$T = 0.65$
and
$T = 0.35$
, respectively. Vertical velocity contours, shown in the vicinity of the sidewall at
$x \approx L_x$
and on the mid-vertical plane at
$x = L_x/2$
, as well as temperature contours on the sidewall at
$y \approx L_y$
are consistent with the double-roll structure spanning the width of the domain.
In this paper, we perform direct numerical simulations (DNS) of 3-D convection for
$\textit{Pr} = 1$
and
$\textit{Ra}$
up to
$10^{10}$
in a rectangular box of aspect ratios
$\varGamma _x = L_x/H = 2.4$
and
$\varGamma _y = L_y/H = 0.8$
. As the first task, we explore the large-scale structure, and find it to be two horizontally stacked cells over the entire domain, with fixed orientations for the entire duration of simulations (see figure 1). This flow morphology results in nearly well-defined and spatially fixed ejection, shear and impact regions on the bottom and top plates of the cell. Further, deeper in the bulk, the region near the mid length experiences continuous movement of thermal plumes and is thus an active region, whereas regions approximately a quarter length from the walls are relatively quiet; see figure 2. As a second task, we explore the characteristics of the flow in the bulk and the boundary layer corresponding to the various conditions stated above, especially the scaling characteristics of velocity and thermal fluctuations. One of our findings in this respect is that, despite these differences, the global transports of heat and momentum are practically identical, within error bars, to the one in horizontally unconstrained layer configurations (Samuel et al. Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024).
An illustration depicting various regions of the flow in the present configuration, using a representative mean temperature field in the
$y = L_y/2$
plane for
$\textit{Ra} = 10^9$
. Ejection, impact and shear regions are observed at the horizontal plates, while active and quiet regions exist deeper in the bulk of the convective layer.
The outline of the manuscript is as follows. Section 2 is an account of the numerical method – the account is brief because it is by now standard – and a summary of simulation parameters. Section 3 considers global transport in regions of high and low plume activity, while § 4 considers boundary layers. The paper concludes with summary remarks in § 5. A major finding of the present work is that the scaling properties depend strongly on whether or not one considers regions of high or low plume activity, with the temperature fluctuations and the normalised dissipation rates decaying more slowly with the Rayleigh number in regions of higher activity.
2. Numerical method
The following non-dimensionalised governing equations of RBC are integrated in time:
\begin{eqnarray} \frac {\partial {\boldsymbol u}}{\partial t} + {\boldsymbol u} \boldsymbol{\cdot }\boldsymbol{\nabla }{\boldsymbol u} & = & - \boldsymbol{\nabla }\!p + T \hat {z} + \sqrt {\frac {\textit{Pr}}{Ra}} {\nabla} ^2 {\boldsymbol u}, \end{eqnarray}
\begin{eqnarray} \frac {\partial T}{\partial t} + {\boldsymbol u} \boldsymbol{\cdot }\boldsymbol{\nabla }T & = & \frac {1}{\sqrt {\textit{RaPr}}} {\nabla} ^2 T, \end{eqnarray}
\begin{eqnarray} \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol u} & = & 0. \end{eqnarray}
Here,
${\boldsymbol u} = (u_x, u_y, u_z)$
,
$T$
and
$p$
are the velocity, temperature and pressure fields, respectively. These equations are non-dimensionalised using
$H$
,
$\Delta T$
, the free-fall velocity
$u_{\!f} = \sqrt {\alpha g \Delta T H}$
and the free-fall time
$t_{\!f} = H/u_{\!f}$
as the length, temperature, velocity and time scales, respectively. The Rayleigh and Prandtl numbers are defined as
$\textit{Ra} = \alpha g \Delta T H^3/(\nu \kappa )$
and
$\textit{Pr} = \nu /\kappa$
;
$\Delta T$
is the temperature difference between the bottom and top plates;
$\alpha , \nu , \kappa$
are respectively the isobaric thermal expansion coefficient, kinematic viscosity and thermal diffusivity;
$g$
is the acceleration due to gravity.
We conduct DNS of RBC for
$\textit{Pr} = 1$
and
$\textit{Ra}$
up to
$10^{10}$
using an open-source solver Nek5000 (Fischer Reference Fischer1997; Scheel, Emran & Schumacher Reference Scheel, Emran and Schumacher2013), which utilises the spectral element scheme. The rectangular box is decomposed into a finite number of elements
$N_e$
and each element is refined further by using a spectral interpolation with Legendre polynomials of order
$N$
in all directions. Thus, the flow domain is discretised using
$N_e N^3$
mesh cells. The velocity field satisfies the no-slip condition on all of the walls. For the temperature field, an adiabatic condition is imposed on the sidewalls and an isothermal condition on the horizontal plates. Mesh clustering near all the walls is used to capture rapid spatial variation of the temperature and velocity fields. The flow is allowed to evolve under the Oberbeck–Boussinesq approximation.
Crucial parameters of DNS are listed in table 1. We ensure that the ratio of local vertical grid spacing and the Kolmogorov scale
$\eta (z)$
remains lower than 1.8 (Pandey et al. Reference Pandey, Krasnov, Schumacher, Samtaney and Sreenivasan2022a
), where
$\eta =(\nu ^3 / \varepsilon _u)^{1/4}$
is nominally the finest scale in the velocity field and
$\varepsilon _u$
is the kinetic energy dissipation rate computed as
\begin{equation} \varepsilon _u (\boldsymbol x) = \frac {\nu }{2} \sum _{l,m=1}^3 \left ( \frac {\partial u_l}{\partial x_m} + \frac {\partial u_m}{\partial x_l} \right )^2 \!, \end{equation}
with
$l, m \in \lbrace x,y,z \rbrace$
. The thermal dissipation rate is computed as
\begin{equation} \varepsilon _T (\boldsymbol x) = \kappa \sum _{j=1}^3 \left ( \frac {\partial T}{\partial x_j} \right )^2 \! . \end{equation}
Parameters of DNS conducted in a domain of dimensions
$L_x:L_y:H = 2.4:0.8:1$
for
$\textit{Pr} = 1$
. Total number of mesh cells
$N_e \times N^3$
in the flow domain; Nusselt numbers
$\textit{Nu}$
,
$\textit{Nu}_{\varepsilon _u}$
,
$\textit{Nu}_{\varepsilon _T}$
and
$\textit{Nu}_W$
computed using (2.6), (2.7), (2.8) and (2.9), respectively; Reynolds number
$\textit{Re}$
computed using (3.2); and the total integration time
$t_{\textit{sim}}$
in the statistically steady state. Nusselt numbers computed using various methods agree to within 1 % for all the simulations, indicating a good spatial and temporal convergence (Pandey & Sreenivasan Reference Pandey and Sreenivasan2025). Error bars in global quantities represent the standard deviation.
Global dissipation rates, i.e. the volume- and time-averaged values, are related to global heat transport in RBC (Howard Reference Howard1972). The Nusselt number
$\textit{Nu}$
, which is the ratio of the total heat transport to that only by conduction, is computed as
\begin{equation} \textit{Nu} = 1 + \sqrt {\textit{RaPr}} \, \langle u_z T \rangle _{V,t}, \end{equation}
where
$\langle \boldsymbol{\cdot }\rangle _{V,t}$
is the combined volume–time average. Nusselt numbers computed from the dissipation rates are defined as
\begin{eqnarray} \textit{Nu}_{\varepsilon _u} & = & \sqrt {\textit{RaPr}} \, \langle \varepsilon _u \rangle _{V,t} + 1 , \end{eqnarray}
\begin{eqnarray} \textit{Nu}_{\varepsilon _T} & = & \sqrt {\textit{RaPr}} \, \langle \varepsilon _T \rangle _{V,t} . \end{eqnarray}
The Nusselt number at the walls is computed from the vertical derivative of the temperature
\begin{equation} \textit{Nu}_W = - \left . \frac { \partial \langle T \rangle _{A,t}} {\partial z} \right | _{z=0,H} . \end{equation}
The heat transports computed using different methods agree to within 1 % (see table 1). This concurrence strongly indicates that the resolution is adequate for faithful DNS (Stevens, Verzicco & Lohse Reference Stevens, Verzicco and Lohse2010; Shishkina et al. Reference Shishkina, Stevens, Grossmann and Lohse2010). Note that the computational constraints limit the integration time that could be realised for simulations at higher Rayleigh numbers. However, to ensure the adequacy of the duration of simulations, we have checked the convergence of global quantities by computing the difference between the two halves of the dataset. We find that the differences in results remain of the order of 2 % (in some cases a tenth as small).
3. Global transport, flow structure and fluctuations in high and low plume activity regions of the bulk flow
3.1. Global transports of heat and momentum
The variation of
$\textit{Nu}$
with
$\textit{Ra}$
in the rectangular box is shown in figure 3(a). The best fit to the data in this range of
$\textit{Ra}$
follows
$\textit{Nu} = 0.15 \textit{Ra}^{0.29}$
. We also plot the Nusselt numbers from DNS studies in horizontally periodic square cuboids of
$\varGamma = 32$
(Blass et al. Reference Blass, Verzicco, Lohse, Stevens and Krug2021) and
$\varGamma = 4$
(van Reeuwijk et al. Reference van Reeuwijk, Jonker and Hanjalić2008), both for
$\textit{Pr} = 1$
, as well as data for
$\textit{Pr} = 0.7$
in cylindrical cells of
$\varGamma = 1$
(Scheel, Kim & White Reference Scheel, Kim and White2012; Wagner, Shishkina & Wagner Reference Wagner, Shishkina and Wagner2012) and in a horizontally periodic cuboid of
$\varGamma = 4$
(Samuel et al. Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024). Here,
$\varGamma = L/H$
is the aspect ratio. This power-law exponent is close to
$2/7$
obtained in the mixing-zone model (Castaing et al. Reference Castaing, Gunaratne, Kadanoff, Libchaber and Heslot1989) and to that by the turbulent boundary layer model (Shraiman & Siggia Reference Shraiman and Siggia1990; Siggia Reference Siggia1994). We conclude that the heat transport is largely insensitive to the aspect ratio (for
$\varGamma \geqslant 0.5$
) and the shape of the convection cell, and also to whether the domain is closed or horizontally periodic. The mixing model of Castaing et al. (Reference Castaing, Gunaratne, Kadanoff, Libchaber and Heslot1989) is linked to boundary layer only indirectly by being the source of plumes, whereas the Shraiman–Siggia model (Shraiman & Siggia Reference Shraiman and Siggia1990) relies strongly on our empirical understanding of the turbulent boundary layer. The two models do not share much physics at first look, yet they produce the same result on the
$\textit{Nu}-Ra$
relation. A plausible inference to draw from this observation is that the dynamical aspects of boundary layers and the mixing regions are intimately intertwined for these Rayleigh numbers.
(a) Global heat transport,
$\textit{Nu}$
, increases with the thermal driving as
$\textit{Ra}^{0.29}$
in the rectangular box. Present results agree excellently with those for
$\textit{Pr} = 1$
from van Reeuwijk et al. (Reference van Reeuwijk, Jonker and Hanjalić2008) and Blass et al. (Reference Blass, Verzicco, Lohse, Stevens and Krug2021) obtained respectively in
$\varGamma = 32$
and
$\varGamma = 4$
horizontally periodic square cuboids as well as those for
$\textit{Pr} = 0.7$
from Samuel et al. (Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024). (b) Reynolds number based on
$u_{\textit{rms}}$
increases almost as
$\sqrt {Ra}$
. Dashed line stands for the scaling of the Reynolds number based on
$U_{ {max}}$
. Inset shows the normalised Reynolds number
$\textit{Re} \textit{Ra}^{-0.49} \textit{Pr}^{2/3}$
as a function of
$\textit{Ra}$
. The scaling exponents agree well with those observed in literature but the prefactor depends on the shape and size of the convective layer.
We now consider the Reynolds number
$\textit{Re}$
, which is commonly estimated using the root mean square (r.m.s.) velocity defined as
\begin{equation} u_{\textit{rms}} = \sqrt {\langle u_x^2 + u_y^2 + u_z^2 \rangle _{V,t}} . \end{equation}
The Reynolds number based on
$u_{\textit{rms}}$
and
$H$
\begin{equation} \textit{Re} = u_{\textit{rms}} \sqrt {\textit{Ra}/\textit{Pr}} \end{equation}
is plotted in figure 3(b) as a function of
$\textit{Ra}$
. The best fit to the data yields
$0.19 \textit{Ra}^{0.49}$
. We also compute the Reynolds number based on
$U_{\textit{max}}$
, the maximum horizontal velocity observed close the horizontal plates (Niemela et al. (Reference Niemela, Skrbek, Sreenivasan and Donnelly2001), Lam et al. (Reference Lam, Shang, Zhou and Xia2002); see § 4.1 for a concrete definition), and find good agreement with (3.2). The best fit to the data shows that the
$U_{\textit{max}}$
-based
$\textit{Re}$
follows
$0.17 \textit{Ra}^{0.50}$
. A discernible difference between the two Reynolds numbers occurs only at high Rayleigh numbers, for which
$U_{\textit{max}}$
is slightly higher than
$u_{\textit{rms}}$
.
For comparison, data from van Reeuwijk et al. (Reference van Reeuwijk, Jonker and Hanjalić2008) and Blass et al. (Reference Blass, Verzicco, Lohse, Stevens and Krug2021) are included in figure 3(b), for which the best fits are
$0.29 \textit{Ra}^{0.48}$
and
$0.18 \textit{Ra}^{0.49}$
. Reynolds numbers from Samuel et al. (Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024) are systematically larger – presumably because their configuration is horizontally periodic (with the aspect ratio of 4). We also show the normalised Reynolds number
$\textit{Re} \textit{Ra}^{-0.49} \textit{Pr}^{2/3}$
in the inset of figure 3(b). Although the Reynolds numbers in the closed cuboid are smaller than those in the horizontally periodic cuboids, the data from Blass et al. (Reference Blass, Verzicco, Lohse, Stevens and Krug2021) and Samuel et al. (Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024) are close. The scaling exponent of
$\textit{Re}$
in different convection domains is approximately 0.50, but the prefactors are different: the magnitude of the Reynolds number, not its scaling, depends on the shape and size of the container. The
$\textit{Re}-Ra$
scaling in the present work also agrees with those found for moderate-
$\textit{Pr}$
fluids (Scheel et al. Reference Scheel, Kim and White2012; Wagner et al. Reference Wagner, Shishkina and Wagner2012; Pandey & Verma Reference Pandey and Verma2016; Pandey et al. Reference Pandey, Krasnov, Sreenivasan and Schumacher2022b
).
3.2. Flow organisation
We compute the global kinetic energy corresponding to each velocity component,
$E_x = \langle u_x^2/2 \rangle _{V,t}, E_y = \langle u_y^2/2 \rangle _{V,t}, E_z = \langle u_z^2/2 \rangle _{V,t}$
, and plot in figure 4(a) their ratios with respect to the total energy
$E = E_x + E_y + E_z$
. We find that
$E_x$
is systematically larger than
$E_y$
, suggesting that the flow is stronger in the
$xz$
-plane; this results from our convection domain being longer in the
$x$
-direction than in
$y$
. The ratio
$E_x/E_y$
plotted in figure 4(b) shows that
$E_x \gg E_y$
at the lowest
$\textit{Ra}$
(it is approximately 30 for
$\textit{Ra} = 10^5$
, not shown here), suggesting that the flow is nearly two-dimensional at low Rayleigh numbers, as in extended and weakly forced convection (Busse Reference Busse1978; Pandey, Scheel & Schumacher Reference Pandey, Scheel and Schumacher2018). The flow becomes approximately three-dimensional in the
$\textit{Ra}$
-range between
$10^6$
and
$3 \times 10^6$
, as can be evidenced by the evolution of
$E_x$
,
$E_y$
and
$E_z$
with
$\textit{Ra}$
in this region. The asymptotic ratio
$E_x \approx 3 E_y$
is in the same proportion to the ratio of the configurational geometry.
(a) Fractions of the kinetic energy corresponding to velocity components as functions of
$\textit{Ra}$
exhibit non-monotonic behaviour. The weakest component is
$E_y$
, whereas
$E_x \approx E_z$
, which indicates that the flow evolves primarily in the
$xz$
-plane. (b) Ratio
$E_x/E_y$
as a function of
$\textit{Ra}$
; after a rapid initial decline,
$E_x/E_y$
grows from
$2.5$
to 4 for
$\textit{Ra}$
between
$10^7$
and
$10^{10}$
, implying that the flow becomes increasingly horizontally anisotropic. The horizontal line indicates that
$E_x$
remains more or less close to
$3 E_y$
for a wide range of
$\textit{Ra}$
.
In the following, we restrict comments to the 3-D range of Rayleigh numbers, even though we present data for all of the Rayleigh numbers considered.
A snapshot of the flow shown in figure 1 reveals two counter-rotating large-scale structures stacked along the
$x$
-direction. In the region
$x \approx L_x/2$
, the flow at the bottom plate comprises hot rising plumes for some Rayleigh numbers, and cold falling plumes for others. Thus, in this study, the ejection region at the bottom plate is located at
$x \approx L_x/2$
and the impact region at
$x \approx 0$
or
$ L_x$
, depending upon the specific
$\textit{Ra}$
. Irrespective of
$\textit{Ra}$
, regions
$x \approx L_x/4$
and
$x \approx 3 L_x/4$
on both bottom and top plates correspond to large shear; here, plumes are primarily carried along with the large-scale structures parallel to the horizontal plates. See also figure 2. The stable two-roll structure establishes a mean flow in the rectangular cuboid having a dominant
$x$
-component of velocity.
3.3. Fluctuations in high and low activity bulk regions
Hot or cold plumes are continuously emitted in the region
$x \approx L_x/2$
(see figures 1 and 2) and their emission affects the flow through the bulk of the flow. This leads us to explore the position-dependent flow characteristics in the bulk. As the flow evolves primarily in the
$xz$
-plane, we choose three cubic regions of volume
$V_b = (0.4 H)^3$
nested in the bulk, with the regions centred at fixed
$y_b = L_y/2, z_b = H/2$
, but varying
$x_b = \lbrace L_x/4, L_x/2, 3L_x/4 \rbrace$
. We compute the instantaneous r.m.s. velocity and temperature as well as mean kinetic energy and thermal dissipation rates in the bulk regions and find that all the quantities in quiet regions (at
$x_b = L_x/4$
and
$x_b = 3L_x/4$
) fluctuate about nearly the same mean. However, those in active region at
$x_b = L_x/2$
fluctuate about a higher mean. The local flow states at
$x_b = L_x/4$
and
$x_b = 3L_x/4$
are similar and differ from that at
$x_b = L_x/2$
. This behaviour is true for all Rayleigh numbers, so we combine the data for
$x_b = L_x/4$
and
$x_b = 3L_x/4$
. For reference, we call attention to figure 2 for the locations of the quiet and active regions.
Scaling relations in quiet (red) and active (blue) bulk regions. (a) Reynolds numbers based on the
$x$
-component (filled symbols) and
$z$
-component (open symbols) of velocity. The r.m.s. temperature fluctuation (b), kinetic energy dissipation rate (c) and thermal dissipation rate (d) are larger, and decrease slowly with
$\textit{Ra}$
, in the active region than those in the quiet region. Green diamonds in panels (b), (c), (d) are the corresponding quantities in the entire domain.
In figure 5(a), we show the Reynolds numbers in the bulk regions based on the
$x$
- and
$z$
-components of velocity, computed as
\begin{equation} \textit{Re}_x^{\textit{bulk}} = 0.4 \, \sqrt { \langle u_x^2 \rangle _{V_b,t} \, \textit{Ra}/\textit{Pr}} \quad \mbox{and} \quad \textit{Re}_z^{\textit{bulk}} = 0.4 \,\sqrt { \langle u_z^2 \rangle _{V_b,t} \, \textit{Ra}/\textit{Pr}} . \end{equation}
(The factor 0.4 arises because the side length of
$0.4 H$
is used to define the Reynolds number.) The relative magnitudes of these Reynolds numbers depend on the region being considered. In quiet regions,
$\textit{Re}_x^{\textit{bulk}}$
remains larger than
$\textit{Re}_z^{\textit{bulk}}$
, whereas, in the active region,
$\textit{Re}_z^{\textit{bulk}}$
is significantly larger than
$\textit{Re}_x^{\textit{bulk}}$
. The best fits for
$\textit{Re}_x^{\textit{bulk}}$
are
$\textit{Ra}^{0.46}$
and
$\textit{Ra}^{0.49}$
in the quiet and active regions, respectively. Similarly,
$\textit{Re}_z^{\textit{bulk}}$
scales as
$\textit{Ra}^{0.47}$
and
$\textit{Ra}^{0.48}$
, respectively, in quiet and active regions. These findings suggest that the scaling of
$\textit{Re}^{\textit{bulk}}$
is nearly the same in the entire domain, and the increase of
$\textit{Ra}$
affects the velocity fluctuations in all regions of the flow similarly (Xu et al. Reference Xu, Zhang and Xia2024).
We now compute the scalings of r.m.s. temperature for the entire domain as
\begin{equation} \overline {\sigma }_T = \sqrt {\langle T^2 \rangle _{V,t} - \langle T \rangle _{V,t}^2} \end{equation}
and show its variation with
$\textit{Ra}$
in figure 5(b):
$\overline {\sigma }_T$
decreases as
$0.97\textit{Ra}^{-0.137}$
, consistent with Castaing et al. (Reference Castaing, Gunaratne, Kadanoff, Libchaber and Heslot1989), Wu & Libchaber (Reference Wu and Libchaber1992), Niemela et al. (Reference Niemela, Skrbek, Sreenivasan and Donnelly2000) and also with Pandey & Sreenivasan (Reference Pandey and Sreenivasan2025) in a 2-D square for
$\textit{Pr} = 1$
. This feature seems to be common for all known constrained cases. The reason for decreasing
$\overline {\sigma }_T$
with
$\textit{Ra}$
is that the plumes become finer and occupy an increasingly smaller volume, which implies that increasingly larger fractions of fluid volume possess temperatures close to the mean and contribute weakly to
$\overline {\sigma }_T$
.
The r.m.s. temperature in the bulk regions,
$\overline {\sigma }_T^{\textit{bulk}}$
, is computed by using
$\langle \rangle _{V_b}$
in place of
$\langle \rangle _{V}$
in (3.4). Figure 5(b) shows that it remains larger in the active region, for which the best fit to the data yields
$\textit{Ra}^{-0.12}$
, similar to that found by Samuel et al. (Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024). The best fit to the data in quiet region is
$\overline {\sigma }_T \sim \textit{Ra}^{-0.21}$
, decaying substantially faster than in the active region.
If the observed trends for temperature fluctuation in various regions are extrapolated, the r.m.s. will be dominated increasingly by the contribution from active regions. This implies that the thermal plumes generated near the horizontal plates are swept away by the dominant double-roll structure more effectively as
$\textit{Ra}$
increases. This is plausible because the strength of the dominant flow structure increases with
$\textit{Ra}$
(Lui & Xia Reference Lui and Xia1998; Niemela & Sreenivasan Reference Niemela and Sreenivasan2003; Chandra & Verma Reference Chandra and Verma2013; Pandey Reference Pandey2021) in both two and three dimensions. The faster decrease of
$\overline {\sigma }_T^{\textit{bulk}}$
in the quiet region is also compatible with this picture as the thermal structures arriving in this region become less probable.
The kinetic energy dissipation rates in the entire cell as well as in the bulk regions are shown in figure 5(c). The best fit for the former is
$\langle \varepsilon _u \rangle _{V,t} = 0.10 \textit{Ra}^{-0.19}$
, consistent with the Nusselt number scaling shown in figure 3(a) (using the exact relation (2.7)). The dissipation in the active region scales the same as that of the global value,
$\varepsilon _u^{\textit{bulk}} = 0.05 \textit{Ra}^{-0.18}$
, although it is smaller by a factor of approximately two. In the quiet region,
$\varepsilon _u^{\textit{bulk}}$
decays more rapidly with
$\textit{Ra}$
, the best fit being
$0.09 \textit{Ra}^{-0.23}$
. Taken together, figure 5(c) suggests that, with increasing
$\textit{Ra}$
,
$\varepsilon _u^{\textit{bulk}}$
in the active region slowly approaches the global dissipation. Li et al. (Reference Li, Huang, Ni and Xia2021) found, using Lagrangian measurements, that the normalised local dissipation rate in plume-abundant regions scales as
$\textit{Ra}^{1.34}$
. Our data yield
$\textit{Ra}^{1.32}$
in the active region and
$\textit{Ra}^{1.27}$
in the quiet region.
We show the global and bulk thermal dissipation rates in figure 5(d), and note that bulk value is smaller by nearly an order of magnitude in the active region than the global value: the global thermal dissipation is dominated by contributions from the TBLs (Emran & Schumacher Reference Emran and Schumacher2012; Bhattacharya, Samtaney & Verma Reference Bhattacharya, Samtaney and Verma2019). The global dissipation follows a
$0.15 \textit{Ra}^{-0.21}$
scaling, consistent with that of
$\textit{Nu}$
from (2.8). The dissipation in the active region, in addition to being weaker, vanishes more rapidly with
$\textit{Ra}$
; the best fit yields
$\varepsilon _T^{\textit{bulk}} = 0.04 \textit{Ra}^{-0.28}$
. In the quiet region,
$\varepsilon _T^{\textit{bulk}}$
varies even more rapidly with
$\textit{Ra}$
, according to
$0.34 \textit{Ra}^{-0.48}$
. These trends confirm our intuition that vanishingly lower thermal dissipation occurs in the quiet region.
Mean thickness of the viscous layer at the horizontal plates, estimated using the slope method (filled red circles) and the maximum method (open red circles), decreases with
$\textit{Ra}$
. Orange pentagons represent the
$\delta _u^{\textit{max}}$
in a
$\varGamma = 4$
square cuboid from Samuel et al. (Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024). The maximum method yields a thicker VBL as well as a gradual decrease with
$\textit{Ra}$
. Red solid and dashed curves are the best fits to
$\delta _u^{\textit{sl}}$
and
$\delta _u^{\textit{max}}$
, respectively.
4. Boundary layer properties
4.1. Boundary layer thicknesses
We compute vertical profiles of the wall-parallel velocity as
\begin{equation} U = \sqrt {\langle u_x^2 + u_y^2 \rangle _{A,t}}, \end{equation}
and note that
$U$
displays a maximum before decreasing further away into the bulk region. Here,
$\langle \rangle _{A,t}$
stands for the averaging over the entire horizontal cross-section and time. This maximum horizontal velocity, denoted as
$U_{\textit{max}}$
, is used to define the
$U_{\textit{max}}$
-based Reynolds number in § 3.1. The position of
$U_{\textit{max}}$
from the wall is one estimate for the thickness of the viscous boundary layer (VBL) (Lam et al. Reference Lam, Shang, Zhou and Xia2002; Samuel et al. Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024; Pandey & Sreenivasan Reference Pandey and Sreenivasan2024), and this method of estimating the VBL thickness
$\delta _u^{\textit{max}}$
is referred to as the ‘maximum method’. We plot the averaged
$\delta _u^{\textit{max}}$
over both the horizontal plates against
$\textit{Ra}$
in figure 6 and find that it follows a power law
$\delta _u^{\textit{max}} = 0.59 \textit{Ra}^{-0.15}$
over the entire range of
$\textit{Ra}$
. On the other hand, the slope thickness,
$\delta _u^{\textit{sl}}$
, computed as (Qiu & Xia Reference Qiu and Xia1998a
)
\begin{equation} \delta _u^{\textit{sl}} = \frac {U_{\textit{max}}}{(\mathrm{d}U/\mathrm{d}z)_{z=0}}, \end{equation}
and shown in figure 6 as the average from top and bottom plates, yields a power law
$\delta _u^{\textit{sl}} = 0.68 \textit{Ra}^{-0.23}$
. We thus have the approximate Reynolds number scaling as
$\sqrt {Ra}$
and the mean VBL thickness varying as
$\delta _u^{\textit{sl}} \sim \sqrt {\textit{Re}}$
. This scaling suggests that the VBL retains laminar scaling at least up to
$\textit{Ra} = 10^{10}$
(Wagner et al. Reference Wagner, Shishkina and Wagner2012; Blass et al. Reference Blass, Verzicco, Lohse, Stevens and Krug2021).
In figure 7(a), we plot the normalised velocity profiles
$U/U_{\textit{max}}$
as functions of the normalised wall-normal distance
$z/\delta _u^{\textit{max}}$
, and note that the profiles do not collapse within the VBL. Thus, a dependence on
$\textit{Ra}$
remains even after the normalisation, which implies that the wall-parallel velocity
$U$
does not scale self-similarly in the VBL region.
(a) Normalised wall-parallel velocity profile
$U/U_{\textit{max}}$
as a function of the normalised distance
$z/\delta _u^{\textit{max}}$
from the plate. The profiles do not exhibit self-similarity within the VBL. (b) Normalised temperature fluctuation profiles as a function of the normalised distance from the plate collapse well in the TBL region. The profiles are additionally averaged over the bottom and top halves of the cell due to the up–down symmetry of Boussinesq convection.
We now show that, in contrast to velocity profiles, the temperature profiles collapse well within the TBL. We define the thickness of the TBL,
$\delta _\sigma$
by the vertical profiles of the r.m.s. temperature fluctuation (Belmonte et al. Reference Belmonte, Tilgner and Libchaber1993, Reference Belmonte, Tilgner and Libchaber1994; Lui & Xia Reference Lui and Xia1998) as
\begin{equation} \sigma _T = \sqrt {\langle T^2 \rangle _{A,t} - \langle T \rangle _{A,t}^2} . \end{equation}
Figure 7(b) shows that the normalised temperature fluctuation
$\sigma _T/\sigma _T^{\textit{max}}$
as a function of the normalised distance
$z/\delta _\sigma$
collapses well within the TBL.
We plot the thermal fluctuation thickness
$\delta _\sigma$
, averaged over the bottom and top plates, as a function of
$\textit{Ra}$
in figure 8(a) and find that it follows
$\delta _\sigma \sim 3.74 \textit{Ra}^{-0.30}$
. The mean TBL thickness
$\delta _T$
computed using the slope method is nearly the same as the fluctuation thickness. The best fit is
$\delta _T = 3.4 \textit{Ra}^{-0.29}$
, shown as a dashed curve in figure 8(a). Recalling the relation
$\delta _T = 0.5/Nu$
, the variation of
$\delta _T$
with
$\textit{Ra}$
is consistent with that of the measured
$\textit{Nu}$
. The robustness of
$\delta _T$
and
$\delta _\sigma$
, on the one hand, and differences between
$\delta _u^{\textit{sl}}$
and
$\delta _u^{\textit{max}}$
, on the other, show that the TBL can be identified less ambiguously than the VBL, probably because the temperature scale is imposed on the system, while the velocity scale results from the flow dynamics.
The peak temperature fluctuation with
$\textit{Ra}$
is plotted in figure 8(b); the best fit yields a scaling
$\sigma _T^{\textit{max}} = 0.22 \textit{Ra}^{-0.033}$
. We also include
$\sigma _T^{\textit{max}}$
from the
$\varGamma = 4$
square cuboid (Samuel et al. Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024) in figure 8(b), which follows a scaling
$\sigma _T^{\textit{max}} = 0.22 \textit{Ra}^{-0.032}$
. Thus, the observed scaling of
$\sigma _T^{\textit{max}}$
is quite robust. The decreasing
$\sigma _T^{\textit{max}}$
implies that the plume temperature becomes increasingly similar to that of the ambient fluid as
$\textit{Ra}$
increases (Pandey Reference Pandey2021).
(a) Thermal fluctuation thickness
$\delta _\sigma$
(symbols) decreases with
$\textit{Ra}$
with the best fit (solid curve) suggesting a
$3.74 \textit{Ra}^{-0.30}$
scaling for the entire range. The value of
$\delta _\sigma$
remains slightly smaller than the slope thickness
$\delta _T$
that decreases as
$3.4 \textit{Ra}^{-0.29}$
(dashed line). The present data agree well with the thermal fluctuation thickness observed for
$\textit{Pr} = 0.7$
in a square cuboid of
$\varGamma = 4$
(Samuel et al. Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024), shown as orange pentagons. (b) Amplitude of the peak temperature fluctuation,
$\sigma _T^{\textit{max}}$
, decreases as
$\textit{Ra}^{-0.033}$
in the rectangular cuboid and agrees very well with that observed in the
$\varGamma = 4$
square cuboid (Samuel et al. Reference Samuel, Bode, Scheel, Sreenivasan and Schumacher2024). The decreasing trend of
$\sigma _T^{\textit{max}}$
is also similar to that observed in 2-D convection for
$\textit{Pr} = 0.021$
(Pandey Reference Pandey2021).
4.2. Boundary layer evolution on horizontal walls
The flow near the horizontal plates is inhomogeneous and so we study the profiles extracted at various uniformly spaced locations along the length of the cell at the mid-width
$y = L_y/2$
. These local profiles are computed as
\begin{equation} U^{\textit{loc}}(x_0) = \sqrt { \langle u_x^2 + u_y^2 \rangle _{A_0,t} } \quad \mbox{and} \quad T^{\textit{loc}}(x_0) = \langle T \rangle _{A_0,t} \, . \end{equation}
Here,
$A_0$
is a horizontal cross-section of area
$0.05H \times 0.05H$
centred at
$(x_0,y_0 = L_y/2)$
with
$x_0/L_x = \lbrace 2/32, 3/32,\ldots , 30/32 \rbrace$
. (We have verified by varying the area of the averaging window between
$A_0/4$
and
$4 A_0$
that the following results remain robust.) The local profiles were found to be symmetric with respect to the mid-vertical plane
$x = L_x/2$
, as expected, so we average and discuss them at
$x_0 = \lbrace 2 L_x/32, 3 L_x/32,\ldots , L_x/2 \rbrace$
. Further, as the ejection region at the bottom plate is located at
$x \approx L_x/2$
for some
$\textit{Ra}$
, but at
$x \approx 0$
or
$L_x$
for others, to make the discussion clear, we transform the flow such that
$x = L_x/2$
corresponds at all Rayleigh numbers to the ejection and impact regions at the bottom and top plates, respectively. The location
$x = L_x/4$
continues to be the shear region on both plates (see figure 2). The horizontal flow is strong in the shear region and the maximum horizontal velocity decreases as one moves away from there (Qiu & Xia Reference Qiu and Xia1998b
). The local profiles are asymmetric with respect to the mid-horizontal plane except those in the shear region.
Structure of the VBL at the bottom plate. (a) Local thickness estimated using the maximum method normalised with that in the shear region increases in the direction of the wind (indicated by orange arrows). (b) Thickness estimated using the slope method exhibits an increasing trend up to
$x_0 \approx 0.4 L_x$
, beyond which it decreases.
Structure of the TBL at the bottom plate. Both the local fluctuation thickness (a) and the slope thickness (b) normalised with that in the shear region increase in the direction of the wind (indicated by orange arrows). However, the increasing trend breaks down at larger Rayleigh numbers beyond
$x_0 \approx 0.4 L_x$
.
We show in figures 9 and 10 the local boundary layer thicknesses computed from these profiles, normalised by that in the shear region, for normalised distances
$x_0/L_x$
. To avoid clutter, we show data only at the bottom plate. Figure 9 shows the variation of the VBL thickness estimated using the maximum method in panel (a), and that using the slope method in panel (b). Both methods yield local VBL thicknesses that increase, up to
$x_0/L_x \approx 0.4$
in the wind direction, indicated using orange arrows. See also Qiu & Xia (Reference Qiu and Xia1998a
) and Wagner et al. (Reference Wagner, Shishkina and Wagner2012) for similar observations. Figure 9(b) shows that
$\delta _u^{\textit{sl}}(x_0)$
starts to decrease as the ejection region is approached, but not
$\delta _u^{\textit{max}}(x_0)$
. As the ejection region in the present configuration is far from the sidewalls, it is possible that the decline of the local VBL thickness is an artefact of the slope method. Clusters of thermal plumes with locally smaller spacing detach at the edge of the ejection region rather than from its centre, where a stagnation region forms due to colliding mean flows. This behaviour can be probed only by the slope method that builds on vertical temperature derivatives at the wall. With increasing
$\textit{Ra}$
, a larger number of finer-scale plumes aggregate, causing broader and less coherent clustering regions, thus decreasing thickness ratios with increasing
$\textit{Ra}$
.
Scaling of (a) local VBL thickness and (b) TBL thickness as functions of
$\textit{Ra}$
. Red circles are data for the impact region (
$x_0 = L_x/2$
at the top plate), blue triangles are for shear region (
$x_0 = L_x/4$
and
$3 L_x/4$
) and green diamonds are for ejection region (
$x_0 = L_x/2$
at the bottom plate). Best fits reveal that the exponents are the same for the impact and shear regions and smaller than that for the ejection region. Black dashed curve in both (a) and (b) stand for the variation of the mean boundary layer thickness.
As shown in figure 10, TBL thicknesses from both methods increase all the way to
$L_x/2$
up to
$\textit{Ra} = 3 \times 10^6$
but decrease with
$x_0$
in the region
$x_0 \geqslant 0.4 L_x$
for higher
$\textit{Ra}$
. It is also clear that the variation of the TBL thickness over the plate becomes milder as
$\textit{Ra}$
increases, consistent with earlier observations (Lui & Xia Reference Lui and Xia1998; Pandey Reference Pandey2021). Figure 9 further suggests that the VBL thickness varies over the plate more strongly than the TBL thickness.
In figure 11(a), we plot the VBL thicknesses
$\delta _u^{\textit{sl}}(x_0)$
in the impact, shear and ejection regions as functions of
$\textit{Ra}$
, and observe that, as expected,
$\delta _u^{\textit{sl}}(x_0)$
decreases with
$\textit{Ra}$
everywhere. The best fit results in a scaling
$\delta _u^{\textit{sl}}(x_0) = 0.6 \textit{Ra}^{-0.23}$
in the impact region (
$x \approx L_x/2$
at the top plate). The thicknesses in the shear region (blue triangles) are slightly larger but the scaling exponent is nearly the same. In this region, the VBL thickness scales as
$\delta _u^{\textit{sl}}(x_0) = 0.9 \textit{Ra}^{-0.23}$
. Note that the thickness in the shear region corresponds to the average thickness at
$x_0 = L_x/4$
and
$x_0 = 3L_x/4$
both at the bottom and top plates. For comparison, we also show the scaling of the mean VBL thickness
$\delta _u^{\textit{sl}} = 0.68 \textit{Ra}^{-0.23}$
as a black dashed curve in figure 11(a). The VBL thickness in the impact (shear) region remains consistently smaller (larger) than the mean thickness at the entire horizontal plate. The thickness in the ejection region (at
$x \approx L_x/2$
at the bottom plate) declines more rapidly with
$\textit{Ra}$
, with the best fit yielding
$1.9 \textit{Ra}^{-0.28}$
.
Scalings of the local TBL thicknesses are shown in figure 11(b). As for the VBL thicknesses, scalings in the impact and shear regions are similar and differ from that in ejection region. The best fits result in the following scalings:
$\delta _T(x_0) = 2.4 \textit{Ra}^{-0.28}$
in the impact region,
$\delta _T(x_0) = 3.2 \textit{Ra}^{-0.28}$
in the shear region and
$\delta _T(x_0) = 11.8 \textit{Ra}^{-0.35}$
in the ejection region. Local thermal fluctuation thicknesses,
$\delta _\sigma (x_0)$
, scale very similarly to
$\delta _T(x_0)$
in these regions and thus also exhibit a strong difference between the scaling exponents in the ejection and shear regions. The scaling of TBL thickness in the shear region is nearly the same as that for the mean thickness over the entire plate, which is shown as a black dashed curve in figure 11(b). A rapid decrease of TBL thickness in ejection region, compared with those in shear and impact regions, is similar to that observed in low-
$\textit{Pr}$
convection, both in two (Pandey Reference Pandey2021) and three dimensions (Kim et al. Reference Kim, Schindler, Vogt and Eckert2024). As the heat transport is inversely proportional to the TBL thickness, this difference also implies that the local heat transport in the ejection region grows more rapidly with
$\textit{Ra}$
than that in the shear region (Zhu et al. (Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018)) for
$\textit{Pr} = 1$
in 3-D RBC.
5. Conclusions
In convection studies, it is customary to separate the flow into a boundary layer region and the bulk. There is, in general, no ambiguity to the notion when the aspect ratio is very large, which was the initial intent. When applying the concepts to small aspect ratios, the division is not conceptually clear cut because the flow near the boundary is strongly influenced by the vertical walls of the container. We have studied convection where the vertical walls are far enough away from each other so as to reduce their influence on the boundary layer. In this configuration, it is possible to study regions dominated by a strong activity of thermal plumes, quiet regions, as well as regions of shear and impact. The configuration we have chosen is a rectangular cuboid of length
$2.4H$
and width
$0.8H$
, and DNS were performed for
$\textit{Pr} = 1$
and
$3 \times 10^5 \leqslant Ra \leqslant 10^{10}$
. For these conditions, the walls do not strongly influence the plume dynamics, so that the flow characteristics can be explored in plume-abundant and plume-depleted regions, as well as regions of impact and shear.
What we learn is that, while the TBL can be unambiguously identified and is self-similar with respect to
$\textit{Ra}$
, the velocity boundary layer is not; its characteristic thickness depends on the definition. We observed, further, that the horizontal velocity changes gradually or rapidly from the wall depending on whether one considers the ejection, impact or shear region.
We explored the properties of velocity and temperature fluctuations and dissipation rates in active and quiet regions, located deeper in the bulk of the convective layer. The velocity fluctuations exhibit nearly similar characteristics in all these regions, this property being true for the kinetic energy dissipation rate as well. On the other hand, temperature fluctuations, as well as the thermal dissipation rate, do depend on the region of the flow where they are measured. We find that the scalings of all thermal quantities in the active region are shallower (i.e. the corresponding scaling exponents have a smaller magnitude) than those in the quiet region.
We further found that both the viscous and thermal boundary layers become thicker in the direction of the mean wind, although their structure in the ejection region is not monotonic. The behaviour in the ejection region was observed to depend on the definition of the boundary layer thickness in the case of the viscous layer, whereas only on the Rayleigh number in case of the thermal layer. We also studied the scalings of the local boundary layer thicknesses with respect to
$\textit{Ra}$
and observed that both VBL and TBL become thinner more rapidly in the ejection region than those in the shear and impact regions. These results form the bulk of the paper.
Our study also led us to two general qualitative lessons. First, the boundary layers and the bulk are intimately intertwined for the parameter ranges explored here, even though both boundary layers get increasingly thinner with increasing
$\textit{Ra}$
. Thus, theories that consider the bulk or the boundary layers as more important than the other cannot be distinguished for their accuracy. Second, despite these differences of detail, the global heat transport laws are practically the same as those in other configurations of low to moderate aspect ratios. In other words, focusing merely on the Nusselt number characteristics does not help us understand the nature of the flow adequately.
Acknowledgements
This research was carried out on the High Performance Computing resources at New York University Abu Dhabi. The authors also gratefully acknowledge Shaheen II of KAUST, Saudi Arabia (under project numbers k1491 and k1624) for providing computational resources.
Funding
This material is based upon work supported by Tamkeen under the NYU Abu Dhabi Research Institute grant G1502, and by the KAUST Office of Sponsored Research under Award URF/1/4342-01. A.P. also acknowledges financial support from ANRF (formerly SERB), India under the grant SRG/2023/001746. NYU supports the research of K.R.S.
Declaration of interests
The authors report no conflicts of interest.
Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.
