# Boundary layers in thermal convection are fluctuation-dominated

Roshan J. Samuel  Institute of Thermodynamics and Fluid Mechanics, Technische
Universität Ilmenau, D-98684 Ilmenau, Germany  Mathis Bode  Institute for
Advanced Simulation, Jülich Supercomputing Centre, D-52428 Jülich, Germany
Janet D. Scheel  Department of Physics, Occidental College, Los Angeles, CA
90041, USA  Katepalli R. Sreenivasan  Tandon School of Engineering, New York
University, New York, NY 11201, USA  Courant Institute of Mathematical
Sciences, New York University, New York, NY 10012, USA  Center for Space
Science, New York University, Abu Dhabi, 100565, UAE  Jörg Schumacher
Institute of Thermodynamics and Fluid Mechanics, Technische Universität
Ilmenau, D-98684 Ilmenau, Germany  Tandon School of Engineering, New York
University, New York, NY 11201, USA

###### Abstract

We study the dynamics of thermal and momentum boundary layers in three-
dimensional direct numerical simulations of Rayleigh-Bénard convection for the
Rayleigh number range $10^{5}\leq Ra\leq 10^{11}$ and $Pr=0.7$. Using a
Cartesian slab with horizontal periodic boundary conditions and an aspect
ratio of 4, we obtain statistical homogeneity in the horizontal $x$\- and
$y$-directions, thus approximating an infinitely extended system. We observe
that upon canonical use of long-time and area averages, a coherent mean flow
is practically absent. Instead, the velocity field close to the wall is a
collection of differently oriented local shear patches interspersed by shear-
free, incoherent flow regions. These shear patches occupy an area fraction of
approximately 40% for all $Ra$. Rather than resulting in a pronounced mean
with small fluctuations about it, the velocity field is dominated by strong
fluctuations of all three components around a non-existent or weak mean. This
feature is particularly pronounced for $Ra\geq 10^{9}$. We discuss the
consequences of these observations for convection layers with larger aspect
ratios, including boundary layer instabilities and the resulting turbulent
heat transport.

##  I Introduction

Rayleigh-Bénard convection (RBC) is one of the fundamental flow configurations
in fluid turbulence research. The fluid in this configuration is nominally
confined to an infinitely extended layer enclosed between two parallel,
horizontal and impermeable plates separated by a vertical distance $H$
(Rayleigh, 1916). When the fluid layer is heated sufficiently strongly from
below (and cooled from above), buoyancy forces initiate a turbulent fluid
motion that has a statistically preferred state with respect to the direction
of gravity, ${\bm{g}}={g\bm{e}}_{z}$. A central question concerns the amount
of heat and momentum carried through the layer, and their dependencies on the
imposed temperature difference between the top and bottom plates, $\Delta
T=T_{\rm bot}-T_{\rm top}$. The temperature difference is expressed by the
dimensionless Rayleigh number $Ra=g\alpha\Delta TH^{3}/(\nu\kappa)$, where
$\alpha$ is the isobaric expansion coefficient, $\nu$ the kinematic viscosity,
and $\kappa$ the temperature diffusivity (Ahlers et al., 2009a; Chillà and
Schumacher, 2012; Verma, 2018).

Since the RBC system is enclosed by walls at the top and bottom, the viscous
and thermal boundary layers formed on these walls pose a bottleneck for the
global transport of both heat and momentum. Their composition and dynamics at
very high Rayleigh numbers still need to be better understood as emphasized
recently (Iyer et al., 2020; Lindborg, 2023; Shishkina and Lohse, 2023;
Creyssels and Martinard, 2024). For Rayleigh numbers $Ra\gtrsim 10^{10}$, no
laboratory experiment to-date has resolved the dynamic interplay of the
boundary layers and their fluctuations with the basic structural elements,
namely thermal plumes and shear layers. Direct numerical simulations (DNS) are
thus the only way to compare their structure and statistical properties with
predictions from theories for canonical laminar and turbulent boundary layers
(Schlichting and Gersten, 2016). Furthermore, the closed-cell geometry of
high-Rayleigh-number studies (Castaing et al., 1989; Chavanne et al., 1997;
Niemela et al., 2000; Ahlers et al., 2009b; Urban et al., 2012) breaks the
horizontal translation symmetry of the statistics, except possibly when the
aspect ratio is very large (Pandey et al., 2018). Small aspect ratio enforces
a dominant large-scale circulation (LSC) in the cell (Kadanoff, 2001),
manifesting as a relatively coherent shear flow connecting the top and bottom
plates, fluctuating only moderately in its mean orientation (Sreenivasan et
al., 2002; Stevens et al., 2011; Shi et al., 2012; Scheel and Schumacher,
2017). Both aspects take us away from the original question on the heat and
momentum transfer in an infinitely extended plane layer.

In this work, we focus on a configuration that is closer to the original RBC
model of convection between a pair of infinitely extended planes, by using
periodic boundaries in horizontal directions. Simulations with similar
boundary conditions and Prandtl numbers have been done by Kerr (1996) for
$\Gamma=6$ and $Ra\leq 2\times 10^{7}$, by Hartlep et al. (2003) for
$\Gamma=10$ and $Ra\leq 1\times 10^{7}$, by van Reeuwijk et al. (2008) and van
Reeuwijk and Jonker (2008) for $\Gamma=4$ and $Ra\leq 1\times 10^{8}$, and by
De et al. (2018) for $\Gamma=6$ but $Ra\leq 2\times 10^{6}$. Our DNS spans
Rayleigh numbers of six orders of magnitude up to $Ra=10^{11}$ for long
periods of time (see table 1 for details). A canonical mean flow analysis
reveals practically no mean flow; instead, strong velocity fluctuations
dominate the flow at all $Ra$. Fits to the mean velocity profiles result in
very small free-stream velocities $U_{\infty}\sim 10^{-3}$ in terms of the
free fall velocity $U_{f}=\sqrt{g\alpha\Delta TH}$, and thicknesses
$\delta_{\infty}\sim 10^{-2}$ in terms of $H$, resulting in small shear
Reynolds numbers $Re_{\rm shear}\lesssim 1$–10 even for the largest $Ra$.
Furthermore, we analyse fluctuations of the velocity components, determine the
distances of maximum mean-square fluctuation from the wall, and discuss the
resulting Reynolds numbers. A first impression of the complex boundary layer
dynamics can be obtained by the streamline plots in figure 1 at
$z=\delta_{T}/2$ for two $Ra$. The figures indicate a prominent patchiness of
the boundary layer viewed as a whole from the top. We will quantify the
corresponding area fraction and relate it to the large-scale patterns, which
are the turbulent superstructures of convection (Pandey et al., 2018).

##  II Numerical simulations

Figure 1: Streamline snapshots of the horizontal velocity field
$(u_{x},u_{y})$ at $z=\delta_{T}/2$. (a) $Ra=10^{8}$. (b) $Ra=10^{11}$. The
entire cross section of $4H\times 4H$ is shown.

We solve the three-dimensional Boussinesq equations of RBC (Verma, 2018) by
the GPU-based spectral element code nekRS (Fischer et al., 2022). Length,
velocity, and temperature are expressed in units of $H$ and $U_{f}$, and the
outer temperature difference $\Delta T$. No-slip boundary conditions apply for
the velocity field at the plates at $z=0$ and $H$. Table 1 summarizes 7
simulations, all at a Prandtl number $Pr=\nu/\kappa=0.7$ and aspect ratio
$\Gamma=L/H=4$, where $L$ is the horizontal length. The number of collocation
points inside the thermal boundary layer (based on the temperature fluctuation
profiles) is always $N_{BL}\geq 11$. Furthermore, we verified that the Nusselt
numbers $Nu_{\rm vol}$ and $Nu_{\rm wall}$, which are given by combined
volume-time $\langle\cdot\rangle_{V,t}$ and area-time averages
$\langle\cdot\rangle_{A,t}$,

$Nu_{\mathrm{vol}}=1+\sqrt{RaPr}\langle
u_{z}T\rangle_{V,t}\quad\mbox{and}\quad
Nu_{\mathrm{wall}}=-\frac{\partial\langle T\rangle_{A,t}}{\partial
z}\big{|}_{z=0,H}\,,$ (1)

result in practically the same values (table 1).

Figure 2(a,b) plot the Nusselt and Reynolds numbers, $Nu=Nu_{\rm wall}$ and
$Re=U_{\rm rms}\sqrt{Ra/Pr}$ versus the Rayleigh number $Ra$, compensated by
the high-$Ra$ scaling result of Iyer et al. (2020); here $U_{\rm
rms}=\langle{\bm{u}}^{2}\rangle_{V,t}^{1/2}$. We also compare our results with
those of Scheel and Schumacher (2017) in a closed cylinder at $\Gamma=1$.
While the Nusselt numbers collapse fairly well (figure 2(a)) and are in
agreement with previous simulations for similar parameters and those for
$Ra\leq 10^{8}$ from van Reeuwijk et al. (2008), the Reynolds number shows a
strong geometry dependence (figure 2(b)): they show good agreement in scaling
exponents at high $Ra$ but not in the prefactor. Figure 2(c) shows the thermal
fluctuation boundary layer thickness-based Rayleigh number $Ra_{\delta,{\rm
rms}}=g\alpha(\Delta T/2)\delta^{3}_{T,{\rm rms}}/(\nu\kappa)$ versus $Ra$,
together with the power law fit. Statistics in the all runs are obtained for
more than 100 free-fall times $T_{f}=H/U_{f}$; see table 1. We also verified
that the vertical profiles of the temperature fluctuation, fitted by
polynomial orders $p=7$ and 9 at $Ra=10^{10}$ and $p=5,7$ and 9 at
$Ra=10^{11}$, collapse on each other. These observations confirm that the DNS
are sufficiently well-resolved.

Figure 2: (a) Compensated turbulent heat transfer $Nu/Ra^{1/3}$ versus $Ra$. (b) Compensated turbulent momentum transfer $Re/Ra^{0.458}$ versus $Ra$. (c) Thermal fluctuation boundary layer thickness-based Rayleigh number $Ra_{\delta,{\rm rms}}$ versus $Ra$, together with a power law fit. Panels (a,b) contain data from the present DNS and from Scheel and Schumacher (2017) and Iyer et al. (2020), compensated by the fits to the data from Iyer et al. (2020); see also table 1. $Ra$ | $N_{e}$ | $p$ | $N_{BL}$ | $\Delta t$ | $Nu_{\mathrm{vol}}$ | $Nu_{\mathrm{wall}}$ | $Re$ | $Ra_{\delta,{\rm rms}}$
---|---|---|---|---|---|---|---|---
$10^{5}$ | $100\times 100\times 64$ | 5 | 71 | 1000 | $4.27\pm 0.24$ | $4.27\pm 0.15$ | $93.0\pm 2.6$ | $102\pm 19$
$10^{6}$ | $100\times 100\times 64$ | 7 | 57 | 1000 | $8.15\pm 0.35$ | $8.15\pm 0.19$ | $296\pm 5$ | $115\pm 15$
$10^{7}$ | $100\times 100\times 64$ | 9 | 42 | 1000 | $15.6\pm 0.9$ | $15.6\pm 0.3$ | $892\pm 17$ | $144\pm 14$
$10^{8}$ | $150\times 150\times 96$ | 7 | 24 | 1000 | $30.4\pm 1.9$ | $30.4\pm 0.4$ | $2571\pm 47$ | $175\pm 15$
$10^{9}$ | $150\times 150\times 96$ | 9 | 16 | 400 | $60.9\pm 4.4$ | $60.8\pm 0.9$ | $7300\pm 150$ | $209\pm 15$
$10^{10}$ | $250\times 250\times 128$ | 9 | 13 | 200 | $123.5\pm 9.7$ | $123.1\pm 1.4$ | $20510\pm 470$ | $258\pm 36$
$10^{11}$ | $500\times 500\times 256$ | 5 | 11 | 100 | $255\pm 19$ | $256\pm 3$ | $57960\pm 870$ | $277\pm 22$
Table 1: Details of the simulation series. The Rayleigh number $Ra$, the
number of spectral elements $N_{e}$, the polynomial order $p$ on each element
with respect to the space direction, the number of collocation points within
the thermal boundary layer $N_{BL}$, the total averaging time in free-fall
units $\Delta t$, the volume averaged and the wall-averaged Nusselt numbers,
$Nu$ and $Nu_{\mathrm{wall}}$, the Reynolds number $Re$, and the Rayleigh
number based on the thermal fluctuation boundary layer thickness,
$Ra_{\delta,{\rm rms}}$. The total number of mesh cells is $N_{e}\times
p^{3}$. It increases to almost $2.2\times 10^{10}$ collocation mesh cells for
$Ra=10^{11}$. Mean values in the last four columns are accompanied by the
standard deviations.

Figure 3: (a) Planar-averaged profiles of the horizontal velocity component
$u_{x}$ for different averaging intervals for $Ra=10^{9}$. (b) Match of the
profiles from panel (a) to the Blasius profile $f^{\prime}(\eta)$ (dashed
line). Notice the very small values of $U_{\infty}$ and the non-monotonic
convergence of the profiles. (c) Variation of the free-stream velocity
$U_{\infty}$ and boundary layer height $\delta_{\infty}$ obtained from the
fits to the Blasius profile versus averaging time. (d) Time-averaged profiles
$\langle u_{x}\rangle_{A,t}$ for all 7 Rayleigh numbers. (e,f) Combined plane-
time average profile (thick red line) and instantaneous plane average profiles
(80 thin grey lines) of the horizontal velocity fluctuations, $U^{h}_{\rm
rms}(z)$ and $U_{\rm rms}(z)$ for $Ra=10^{9}$.

##  III Mean temperature and velocity profiles

The mean velocity profiles for the horizontal components are obtained by a
combined average over the area $A=L^{2}$ and $N_{t}=20$ statistically
independent realizations of the turbulent flow separated from each other by at
least 5$T_{f}$ as

$\langle u_{i}\rangle_{A,t}(z)=\frac{1}{A\Delta
t}\sum_{i=1}^{N_{t}}\int_{A}u(x,y,z,t_{i})dA\,,$ (2)

for $i=x,y$. Figure 3 displays the result of this analysis for $Ra=10^{9}$. In
panel (a), the mean profile of the $x$-velocity component is shown as a
function of the averaging time which was varied from $\Delta t=400\,T_{f}$ (=
20 snapshots) to $1600\,T_{f}$ (= 80 snapshots). The profile converges
steadily to zero, though not uniformly. A fit of the mean profiles to the two-
dimensional Blasius solution $f(\eta)$ has been applied subsequently; the
solution $f(\eta)=\psi(x,z)/\sqrt{xU_{\infty}\nu}$ and $\eta=z/\delta(x)$
(Schlichting and Gersten, 2016). Here, $\psi(x,z)$ is the stream function and
$u_{x}=\partial\psi/\partial z$ and $u_{z}=-\partial\psi/\partial x$. Recall
that at distance $\eta=5$ the Blasius profile reaches a streamwise velocity
magnitude of $f^{\prime}(\eta)=0.99U_{\infty}$.

In the absence of a definable leading edge distance $x$, we match $\langle
u_{x}\rangle_{A,t}(z)/U_{\infty}$ and $z/\delta_{\infty}$ to
$f^{\prime}(\eta)$ which is reported in panels (b,c) of figure 3. Therefore,
the numerical profiles are rescaled such that the first local maximum of
$\langle u_{x}\rangle_{A,t}$ corresponds to $U_{\infty}$. The fits are shown
again for different time intervals. It is seen in panel (c) that the velocity
scale $U_{\infty}$ drops to a very small magnitude $\sim 10^{-3}U_{f}$ while
$\delta_{\infty}\lesssim 0.1H$ in case of $Ra=10^{9}$. In table 2, the results
for all Rayleigh numbers and both horizontal components are listed. No clear
trend of the velocity $U^{x,y}_{\infty}$ with Rayleigh number is detectable,
the magnitude is between ${\cal O}(10^{-3})$ and ${\cal O}(10^{-2})$. The
boundary layer thickness parameters vary as well when $x$\- and $y$-directions
are compared at fixed $Ra$. They decrease with increasing Rayleigh number.
Furthermore, we calculate the corresponding shear Reynolds numbers $Re_{\rm
shear}=U_{\infty}\delta_{\infty}/\nu$, which are found to be very small for
all cases. We have also repeated this analysis for a combination of both
components (not shown), though not for the magnitude which would become
insensitive to different directions of local flow patterns. Panel (d) of
figure 3 shows the mean profiles for all 7 simulation runs.

Panels (e) and (f) show the vertical profiles of root-mean-square velocities
in red,

$U_{\rm rms}^{h}(z)=\sqrt{\langle
u_{x}^{2}+u_{y}^{2}\rangle_{A,t}}\quad\mbox{and}\quad U_{\rm
rms}(z)=\sqrt{\langle u_{x}^{2}+u_{y}^{2}+u_{z}^{2}\rangle_{A,t}}\,,$ (3)

together with profiles obtained in the individual snapshots in gray.

$Ra$ | $U^{x}_{\infty}$ | $U^{y}_{\infty}$ | $\delta^{x}_{\infty}$ | $\delta^{y}_{\infty}$ | $Re^{x}_{\rm shear}$ | $Re^{y}_{\rm shear}$
---|---|---|---|---|---|---
$10^{5}$ | $4.8\times 10^{-4}$ | $1.7\times 10^{-3}$ | $9.1\times 10^{-1}$ | $6.5\times 10^{-1}$ | $0.16$ | $0.43$
$10^{6}$ | $2.5\times 10^{-3}$ | $8.4\times 10^{-3}$ | $1.7\times 10^{-1}$ | $1.5\times 10^{-1}$ | $0.51$ | $1.5$
$10^{7}$ | $3.0\times 10^{-3}$ | $2.7\times 10^{-3}$ | $8.2\times 10^{-2}$ | $8.1\times 10^{-2}$ | $0.93$ | $0.82$
$10^{8}$ | $1.2\times 10^{-3}$ | $6.6\times 10^{-3}$ | $5.8\times 10^{-2}$ | $3.7\times 10^{-2}$ | $0.80$ | $2.9$
$10^{9}$ | $2.5\times 10^{-4}$ | $9.3\times 10^{-4}$ | $5.1\times 10^{-2}$ | $2.7\times 10^{-2}$ | $0.48$ | $0.94$
$10^{10}$ | $3.3\times 10^{-3}$ | $1.8\times 10^{-3}$ | $1.3\times 10^{-2}$ | $1.1\times 10^{-2}$ | $5.1$ | $2.4$
$10^{11}$ | $3.6\times 10^{-3}$ | $1.4\times 10^{-4}$ | $5.8\times 10^{-3}$ | $7.0\times 10^{-3}$ | $7.9$ | $0.38$
Table 2: Fit parameters $U^{x}_{\infty}$ and $U^{y}_{\infty}$ as well as the
corresponding thicknesses $\delta^{x}_{\infty}$ and $\delta^{y}_{\infty}$ of
the combined plane-time averaged mean profiles of the horizontal velocity
components $\langle u_{x}\rangle_{A,t}$ and $\langle u_{y}\rangle_{A,t}$ to
the Blasius profile $U_{\infty}f^{\prime}(\eta)$ with
$\eta=z/\delta_{\infty}$. The corresponding shear Reynolds numbers
$Re^{x,y}_{\mathrm{shear}}=U^{x,y}_{\infty}\delta^{x,y}_{\infty}/\nu$ are also
given.

Figure 4 summarizes mean profiles for all 7 simulation runs, the mean
temperature profile and the root-mean-square profiles of temperature,
horizontal velocity components, and all three velocity components. The
temperature fluctuation profile is obtained by $T_{\rm
rms}(z)=\langle\theta^{2}\rangle^{1/2}_{A,t}$ with
$\theta({\bm{x}},t)=T({\bm{x}},t)-\langle T\rangle_{A,t}(z)$. The
correspondingly related characteristic scales are indicated by horizontal
dashed lines and detailed in table 3. It is seen that the thermal boundary
layer thickness $\delta_{T}=1/(2Nu)$ is slightly larger than the distance from
the wall of the maximum of the temperature fluctuation profile,
$\delta_{T,{\rm rms}}$. Increasingly larger are distances from the wall to the
maxima of the velocity fluctuation profiles, obeying a ratio of
$\delta_{U,{\rm rms}}/\delta_{T,{\rm rms}}\approx 2$ for $Ra=10^{5}$ up to
approximately 14 for $Ra=10^{11}$; see panel (d).

Figure 4: Mean temperature profiles compared to each other for all 7 runs. We show the mean profile of temperature $\langle T\rangle_{A,t}(z)$ in panel (a), root-mean-square profile of temperature $T_{\rm rms}(z)$ in (b), root-mean-square profile with respect to the horizontal velocity components, $U^{h}_{\rm rms}(z)$ in (c), and with respect to all three components $U_{\rm rms}(z)$ in (d). The corresponding height scales of these profiles are indicated by dashed horizontal lines and listed in table 3. $Ra$ | $\delta_{T}$ | $\delta_{T,{\rm rms}}$ | $\delta^{h}_{U,{\rm rms}}$ | $\delta_{U,{\rm rms}}$
---|---|---|---|---
$10^{5}$ | ($1.17\pm 0.04)\times 10^{-1}$ | ($1.27\pm 0.08)\times 10^{-1}$ | ($1.23\pm 0.07)\times 10^{-1}$ | ($2.5\pm 1.4)\times 10^{-1}$
$10^{6}$ | ($6.13\pm 0.15)\times 10^{-2}$ | ($6.13\pm 0.27)\times 10^{-2}$ | ($8.80\pm 0.54)\times 10^{-2}$ | ($1.8\pm 1.2)\times 10^{-1}$
$10^{7}$ | ($3.21\pm 0.06)\times 10^{-2}$ | ($3.07\pm 0.10)\times 10^{-2}$ | ($6.99\pm 0.58)\times 10^{-2}$ | ($1.07\pm 0.51)\times 10^{-1}$
$10^{8}$ | ($1.64\pm 0.02)\times 10^{-2}$ | ($1.52\pm 0.04)\times 10^{-2}$ | ($5.08\pm 0.35)\times 10^{-2}$ | ($6.51\pm 0.92)\times 10^{-2}$
$10^{9}$ | ($8.22\pm 0.12)\times 10^{-3}$ | ($7.47\pm 0.18)\times 10^{-3}$ | ($3.52\pm 0.28)\times 10^{-2}$ | ($4.3\pm 2.0)\times 10^{-2}$
$10^{10}$ | ($4.06\pm 0.05)\times 10^{-3}$ | ($3.72\pm 0.17)\times 10^{-3}$ | ($2.32\pm 0.20)\times 10^{-2}$ | ($2.82\pm 0.37)\times 10^{-2}$
$10^{11}$ | ($1.96\pm 0.02)\times 10^{-3}$ | ($1.77\pm 0.05)\times 10^{-3}$ | ($2.12\pm 0.24)\times 10^{-2}$ | ($2.61\pm 0.37)\times 10^{-2}$
Table 3:  Differently defined boundary layer heights of temperature and
velocity according to the vertical mean profiles plotted in figure 4. The
standard deviation accompanies each mean value.

##  IV Scaling laws of fluctuations with Rayleigh number

Figure 5: Scaling of the root-mean-square velocities and temperature with
respect to Rayleigh number. All values are obtained by a combined average with
respect to volume and time. (a) Horizontal velocity components. (b) Vertical
velocity component. (c) All three velocity components. (d) Temperature. Dashed
lines denote power law fits to the data. Blue open circles stand for data
obtained in the bulk volume $V_{b}$. In panels (b–d), we add data from Iyer et
al. (2020) with separate fits. The corresponding temperature data are obtained
for the cylindrical cell interior, $r\leq 0.03H$ and $Ra\geq 10^{9}$.

Figure 5 summarizes the root-mean-square fluctuations of the three velocity
components and the temperature. They are obtained by a combined average with
respect to the full volume $V=L^{2}H$ and time, e.g. $u_{x,{\rm rms}}=\langle
u_{x}^{2}\rangle^{1/2}_{V,t}$. The quantity $U_{\rm rms}$ denotes again the
fluctuations with respect to all three velocity components. It is seen that
the dependence of the velocity fluctuations on the Rayleigh number is very
weak with $\beta\lesssim 0.042$. The temperature fluctuations drop with a
smaller power law exponent, $T_{\rm rms}\sim Ra^{-\beta}$, which is found to
be $\beta=0.118$ for the present data. This exponent is slightly smaller in
magnitude than those reported in experiments in cylindrical cells of aspect
ratio 1/2. For comparison, Castaing et al. (1989), Niemela et al. (2000) and
Wu and Libchaber (1992) report exponents of $\beta\approx 0.145$. We also
analysed the temperature fluctuations in the bulk of the layer, which takes a
volume average with respect to $V_{b}=L^{2}\times[0.4,0.6]$ and time. The
exponent changes to $\beta=0.140$ which is closer to the experiments. We have
verified that a variation of the thickness of the bulk volume $V_{b}$ does not
alter the results significantly.

In panels (b–d) of figure 5, we added data from the DNS of Iyer et al. (2020)
for comparison, which were obtained in a closed slender cylindrical cell of
aspect ratio $\Gamma=0.1$. It is seen that exponents of the power law fits are
close to those of the present simulation series. The prefactors differ as
expected, because the former DNS data were obtained for a geometrically
constrained convection flow. This finding of nearly the same scaling exponents
suggests a robust geometry-independent trend of all fluctuations with respect
to the Rayleigh number. Geometry-specific aspects mostly affect the prefactor.

##  V Decomposition into coherent and incoherent boundary layer regions

The orientation of the boundary layer flow varies strongly as shown in figures
6 (a) and (d), where we plot the pointwise orientation angle of the horizontal
velocity $\varphi=\arctan(u_{x}/u_{y})\in[-\pi,\pi]$ for a snapshot at
$Ra=10^{10}$ at $z=\delta_{T}$ and $z=H-\delta_{T}$, respectively. At both
these heights, we break down the horizontal plane into a cover of $10^{4}$
disjoint square boxes of area content $A_{i}=A/10^{4}$, where $A=L^{2}$. We
then calculate the mean horizontal velocity $\bar{\bm{u}}_{h}(A_{i})$ in each
of the $A_{i}$ and decompose the cross section into coherent and incoherent
boundary layer regions for $|\bar{\bm{u}}_{h}(A_{i})|>U^{h}_{\rm
rms}(\delta_{T})$ and $|\bar{\bm{u}}_{h}(A_{i})|\leq U^{h}_{\rm
rms}(\delta_{T})$, respectively. Panels (b, e) of figure 6 show that coherent
shear-dominated patches separated by incoherent flow regions (in gray). The
superposed streamlines indicate the different flow orientations of the shear-
dominated regions. Panels (c,f) show the corresponding snapshots of the
temperature field $T$ at $z=0.1$ above the bottom and $z=0.9$ below the top
which is a distance of 25 $\delta_{T}$ away from the walls.

It is clearly seen that the hotter regions at $z=0.1$ and the colder regions
at $z=0.9$ coincide with the incoherent flow regions at the edge of the
thermal boundary layer. This is the imprint of the turbulent superstructure of
convection in this setup at $\Gamma=4$. We have also determined that the area
fraction of the incoherent regions remains nearly constant at approximately
60% of $A$ for the whole Rayleigh number range, which further supports its
connection to the turbulent superstructure that consists of a few large-scale
circulation rolls which extend over the periodic horizontal boundaries. The
physical interpretation is that the incoherent regions correspond to the
upwelling (downwelling) motion of clustered plumes outside shear dominated
patches. It suggests a vertical coarsening from the plate into the bulk, which
remains to be analysed in the future.

Figure 6: Boundary layer flow structure for $Ra=10^{10}$ at the top and bottom
walls. (a) Snapshot of the pointwise orientation angle of the horizontal
velocity vector at $z=\delta_{T}$. (b) Decomposition of the same snapshot into
coherent and incoherent (in gray) flow regions at $z=\delta_{T}$. For the
coherent shear-dominated regions, we overlay horizontal streamlines and
indicate their local flow direction. (c) Corresponding temperature field at
$z=0.1\approx 25\delta_{T}$. (d) Same snapshot with the pointwise orientation
angle of the horizontal velocity vector at $z=H-\delta_{T}$. (e) Decomposition
similar to (b) for $z=H-\delta_{T}$. (f) Corresponding temperature field at
$z=0.9\approx 25\delta_{T}$.

##  VI Conclusions

Our DNS of the turbulent Rayleigh-Bénard convection encompasses a Cartesian
domain with $\Gamma=4$ and periodic boundary conditions in the horizontal
directions. These simulations up to $Ra=10^{11}$ are aimed at being close to
the original canonical case of a plane convection layer between a pair of
infinitely extended rigid plates. We demonstrated that a standard mean flow
analysis with a subsequent match to laminar boundary layer profiles gives no
systematic profiles; to the extent that we can define them, they are very
small in shear Reynolds numbers. In the long-time limit, which we could follow
for 1600 $T_{f}$ at $Ra=10^{9}$, the velocity mean profiles have to converge
to $\mbox{lim}_{t\to 0}\langle u_{i}\rangle_{A,t}(z)\to 0$ due to homogeneity
in $x$ and $y$. The simulations by Hartlep et al. (2003) (for $Ra\leq 1\times
10^{7}$) also found that the mean flow contained very little energy, but De et
al. (2018) found a long time periodicity to the mean flow for $Ra\leq 2\times
10^{6}$.

Rather than having a mean flow profile with small fluctuations, we are faced
with small mean velocity amplitudes with fluctuations that are up to 3 orders
of magnitude larger when the statistics are taken over finite time intervals
$\Delta t$, as seen from comparisons of table 2 with the data in figure 4. The
strong fluctuations cause the boundary layer thickness scales of temperature
and velocity to differ by an order of magnitude for the highest Rayleigh
numbers, as summarized in table 3. These results hold for all Rayleigh numbers
of the series, but become particularly pronounced for $Ra\geq 10^{9}$, a range
to which most previous larger-aspect-ratio DNS studies did not advance.

Our analysis also revealed that the velocity boundary layer in the present
configuration is a carpet of differently oriented shear dominated patches
interspersed by regions of incoherent mean flow, the latter of which occupies
about 60% of the plate area for all Rayleigh numbers. This heterogeneous
composition crystallizes particularly for $Ra\gtrsim 10^{9}$, underlying again
the importance of DNS with larger aspect ratios and high Rayleigh numbers. The
present results also raise many questions on the possible transition
mechanisms of the boundary layer to a turbulent regime and the possible
consequences for the global heat transfer. The heterogeneity of the velocity
boundary layer detected here suggests to us the prevalence of local, rather
than global, instability mechanisms, which would bring us back to the marginal
stability concept of the boundary layer, see e.g. Howard (1966). However,
figure 2(c) shows a power law fit of $Ra_{\delta,{\rm rms}}=ARa^{\gamma}$ with
a very small exponent $\gamma=0.077$ and $A\approx 41$. The resulting
$Ra_{\delta,{\rm rms}}$ are by at least a factor of 4 smaller than Howard’s
critical Rayleigh number of $Ra_{\delta}\sim{\cal O}(10^{3})$. These open
points will form the subject of further studies including higher Rayleigh
numbers in this configuration.

###### Acknowledgements.

The work of RJS is funded by the European Union (ERC, MesoComp, 101052786).
Views and opinions expressed are however those of the authors only and do not
necessarily reflect those of the European Union or the European Research
Council. The work of JDS was supported by a Mercator Fellowship of the
Deutsche Forschungsgemeinschaft within the Priority Programme DFG-SPP 1881 on
Turbulent Superstructures. The authors gratefully acknowledge the Gauss Centre
for Supercomputing e.V. (https://www.gauss-centre.eu) for funding this project
by providing computing time through the John von Neumann Institute for
Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing
Centre (JSC).

## References

  * Ahlers et al. [2009a] G. Ahlers, S. Grossmann, and D. Lohse.  Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection.  _Rev. Mod. Phys._ , 81(2):503–537, April 2009a.
  * Ahlers et al. [2009b] Guenter Ahlers, Eberhard Bodenschatz, Denis Funfschilling, and J. Hogg.  Turbulent Rayleigh-Bénard convection for a Prandtl number of 0.67.  _J. Fluid Mech._ , 641:157–167, November 2009b.
  * Castaing et al. [1989] B. Castaing, G. Gunaratne, F. Heslot, L. P. Kadanoff, A. Libchaber, S. Thomae, X.-Z. Wu, S. Zaleski, and G. Zanetti.  Scaling of hard thermal turbulence in Rayleigh-Bénard convection.  _J. Fluid Mech._ , 204:1–30, July 1989.
  * Chavanne et al. [1997] X. Chavanne, F. Chillà, Bernard Castaing, B. Hebral, B. Chabaud, and J. Chaussy.  Observation of the ultimate regime in Rayleigh-Bénard convection.  _Phys. Rev. Lett._ , 79(19):3648–3651, 1997.
  * Chillà and Schumacher [2012] F. Chillà and J. Schumacher.  New perspectives in turbulent Rayleigh-Bénard convection.  _Eur. Phys. J. E_ , 35(7):58, 2012.
  * Creyssels and Martinard [2024] M. Creyssels and D. Martinard.  Stability analysis of sheared thermal boundary layers and its implication for modelling turbulent Rayleigh-Bénard convection.  _Eur. J. Mech. B/Fluids_ , 105:97–103, 2024.
  * De et al. [2018] A. K. De, V. Eswaran, and P. K. Mishra.  Dynamics of plumes in turbulent Rayleigh–Bénard convection.  _Eur. J. Mech. B/Fluids_ , 72:164–178, 2018.  doi: 10.1016/j.euromechflu.2018.05.007.
  * Fischer et al. [2022] P. F. Fischer, S. Kerkemeier, M. Min, Y.-H. Lan, M. Phillips, T. Rathnayake, E. Merzari, A. Tomboulides, A. Karakus, N. Chalmers, and T. Warburton.  NekRS, a GPU-accelerated spectral element Navier–Stokes solver.  _Parallel Comput._ , 114:102982, 2022.
  * Hartlep et al. [2003] T. Hartlep, A. Tilgner, and F.H. Busse.  Large scale structures in Rayleigh-Bénard convection at high Rayleigh numbers.  _Phys. Rev. Lett_ , 91:064501, 2003.  doi: 10.1103/PhysRevLett.91.064501.
  * Howard [1966] L. N. Howard.  Convection at high Rayleigh number.  _Applied Mechanics, 11th Congress of Applied Mechanics, Munich_ , pages 1109–1115, 1966.
  * Iyer et al. [2020] K. P. Iyer, J. D. Scheel, J. Schumacher, and K. R. Sreenivasan.  Classical 1/3 scaling of convection holds up to Ra=$10^{15}$.  _Proc. Natl. Acad. Sci. U.S.A._ , 117:7594–7598, February 2020.
  * Kadanoff [2001] L. P. Kadanoff.  Turbulent heat flow: Structures and scaling.  _Phys. Today_ , 54(8):34–39, August 2001.
  * Kerr [1996] R. M Kerr.  Rayleigh number scaling in numerical convection.  _J. Fluid Mech._ , 310:139–179, January 1996.
  * Lindborg [2023] E. Lindborg.  Scaling in Rayleigh-Bénard convection.  _J. Fluid Mech._ , 956:A34, 2023.
  * Niemela et al. [2000] J. J. Niemela, L. Skrbek, K. R. Sreenivasan, and R. J. Donnelly.  Turbulent convection at very high Rayleigh numbers.  _Nature_ , 404:837–840, January 2000.
  * Pandey et al. [2018] A. Pandey, J. D. Scheel, and J. Schumacher.  Turbulent superstructures in Rayleigh-Bénard convection.  _Nat. Commun._ , 9:2118, May 2018.
  * Rayleigh [1916] O. M. Lord Rayleigh.  LIX. On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side.  _Philos. Mag._ , 32:529–546, 1916.
  * Scheel and Schumacher [2017] J. D. Scheel and J. Schumacher.  Predicting transition ranges to fully turbulent viscous boundary layers in low Prandtl number convection flows.  _Phys. Rev. Fluids_ , 2(12):123501, December 2017\.
  * Schlichting and Gersten [2016] H. Schlichting and K. Gersten.  _Boundary-layer theory_.  Springer Berlin, Heidelberg, 2016.
  * Shi et al. [2012] N. Shi, M. S. Emran, and J. Schumacher.  Boundary layer structure in turbulent Rayleigh–Bénard convection.  _J. Fluid Mech._ , 706:5–33, June 2012.
  * Shishkina and Lohse [2023] O. Shishkina and D. Lohse.  Ultimate turbulent thermal convection.  _Phys. Today_ , 76:26–32, 2023.
  * Sreenivasan et al. [2002] K. R. Sreenivasan, A. Bershadskii, and J. J. Niemela.  Mean wind and its reversal in thermal convection.  _Phys. Rev. E_ , 65:056306, 2002.
  * Stevens et al. [2011] R. J. A. M. Stevens, D. Lohse, and R. Verzicco.  Prandtl and Rayleigh number dependence of heat transport in high Rayleigh number thermal convection.  _J. Fluid Mech._ , 688:31–43, 2011.
  * Urban et al. [2012] P. Urban, P. Hanzelka, T Kralik, V. Musilová, A Srnka, and Ladislav Skrbek.  Effect of Boundary Layers Asymmetry on Heat Transfer Efficiency in Turbulent Rayleigh-Bernard Convection at Very High Rayleigh Numbers.  _Phys. Rev. Lett._ , 109(15):154301, October 2012\.
  * van Reeuwijk and Jonker [2008] M. van Reeuwijk and H. J. J. Jonker.  Wind and boundary layers in Rayleigh-Bénard convection. II. Boundary layer character and scaling.  _Phys. Rev. E_ , 77:036312, 2008.
  * van Reeuwijk et al. [2008] M. van Reeuwijk, H. J. J. Jonker, and K. Hanjalic.  Wind and boundary layers in Rayleigh-Bénard convection. I. Analysis and modeling.  _Phys. Rev. E_ , 77(3):036311, 2008.
  * Verma [2018] M. K. Verma.  _Physics of Buoyant Flows: From Instabilities to Turbulence_.  World Scientific, Singapore, 2018.
  * Wu and Libchaber [1992] X. Z. Wu and A. Libchaber.  Scaling relations in thermal turbulence: The aspect ratio dependence.  _Phys. Rev. A_ , 45:842–845, 1992.