1. 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 Reference Rayleigh1916). 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, 
${\boldsymbol {g}} = {g \boldsymbol {e}}_z$. A central question concerns the amount and nature 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_{bot}-T_{top}$. The temperature difference is expressed by the dimensionless Rayleigh number 
$Ra=g\alpha \Delta T H^3/(\nu \kappa )$, where 
$\alpha$ is the isobaric expansion coefficient, 
$\nu$ the kinematic viscosity and 
$\kappa$ the temperature diffusivity (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009b; Chillà & Schumacher Reference Chillà and Schumacher2012; Verma Reference Verma2018).
A second important question is the structure of the velocity and thermal boundary layers on the horizontal walls and their effects on heat transport. 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. Reference Iyer, Scheel, Schumacher and Sreenivasan2020; Lindborg Reference Lindborg2023; Shishkina & Lohse Reference Shishkina and Lohse2023; Creyssels & Martinard Reference Creyssels and Martinard2024). 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 & Gersten Reference Schlichting and Gersten2016). Furthermore, the closed-cell geometry of high-Rayleigh-number studies (Castaing et al. Reference Castaing, Gunaratne, Heslot, Kadanoff, Libchaber, Thomae, Wu, Zaleski and Zanetti1989; Chavanne et al. Reference Chavanne, Chillà, Castaing, Hebral, Chabaud and Chaussy1997; Niemela et al. Reference Niemela, Skrbek, Sreenivasan and Donnelly2000; Ahlers et al. Reference Ahlers, Bodenschatz, Funfschilling and Hogg2009a; Urban et al. Reference Urban, Hanzelka, Kralik, Musilová, Srnka and Skrbek2012) breaks the horizontal translation symmetry of the statistics, except possibly when the aspect ratio is very large (Pandey, Scheel & Schumacher Reference Pandey, Scheel and Schumacher2018). Small aspect ratio enforces a dominant large-scale circulation (LSC) in the cell (Kadanoff Reference Kadanoff2001), manifesting as a relatively coherent shear flow connecting the top and bottom plates, fluctuating only moderately in its mean orientation (Sreenivasan, Bershadskii & Niemela Reference Sreenivasan, Bershadskii and Niemela2002; Stevens, Lohse & Verzicco Reference Stevens, Lohse and Verzicco2011; Shi, Emran & Schumacher Reference Shi, Emran and Schumacher2012; Scheel & Schumacher Reference Scheel and Schumacher2017). 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 more on the second question in a configuration that is closest to the original RBC model of convection between a pair of infinitely extended planes, by using periodic boundaries in both horizontal directions. Simulations with similar boundary conditions and Prandtl numbers have been done by Kerr (Reference Kerr1996) for an aspect ratio 
$\varGamma=6$ and 
$Ra\leq 2\times 10^7$, by Hartlep, Tilgner & Busse (Reference Hartlep, Tilgner and Busse2003) for 
$\varGamma = 10$ and 
$Ra\leq 1\times 10^7$, by van Reeuwijk, Jonker & Hanjalić (Reference van Reeuwijk, Jonker and Hanjalić2008a) and van Reeuwijk, Jonker & Hanjalić (Reference van Reeuwijk, Jonker and Hanjalić2008b) for 
$\varGamma = 4$ and 
$Ra\leq 1\times 10^8$, by De, Eswaran & Mishra (Reference De, Eswaran and Mishra2018) for 
$\varGamma = 6$ but 
$Ra \leq 2\times 10^6$ and by Stevens et al. (Reference Stevens, Blass, Zhu, Verzicco and Lohse2018) for 
$\varGamma \leq 32$ at 
$Ra=10^9$. Our DNS span Rayleigh numbers of six orders of magnitude up to 
$Ra=10^{11}$ for long periods of time (see table 1 for details). The choice of an aspect ratio of 4 for the present study provides a ‘sweet spot’. On the one hand, the domain is large enough that it does not generate strong LSCs, see Niemela & Sreenivasan (Reference Niemela and Sreenivasan2006). On the other hand, it is small enough to allow us to advance to very high Rayleigh numbers, here of up to 
$Ra=10^{11}$, since the required numerical resources grow with 
$\varGamma ^2$. Furthermore, this is the aspect ratio beyond which Nusselt and Reynolds numbers become independent of 
$\varGamma$, as discussed in Stevens et al. (Reference Stevens, Blass, Zhu, Verzicco and Lohse2018). We supplement these results by additional DNS runs at 
$\varGamma =2$ and 8 for 
$Ra=10^9$.
Details of the simulation series. Listed here are the Rayleigh number 
$Ra$, the aspect ratio 
$\varGamma$, the number of spectral elements 
$N_e$, the polynomial order 
$p$ on each element with respect to each space direction, the number of collocation points within the thermal boundary layer 
$N_{BL}$, the total averaging time in free-fall units 
$\tau _{total}$, the volume-averaged and wall-averaged Nusselt numbers, 
$Nu_{vol} (=Nu)$ and 
$Nu_{wall}$, the Reynolds number 
$Re$ and the Rayleigh number based on the thermal fluctuation boundary layer thickness, 
$Ra_{\delta,{rms}}$. The three DNS runs with dagger symbols at 
$Ra=10^9$ have a different vertical spectral element grid stretching than the fourth DNS run at this Rayleigh number. 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. The Prandtl number of all runs is 
$Pr=0.7$.
A canonical mean-flow analysis reveals practically no global mean flow; instead, strong velocity fluctuations dominate the flow at all 
$Ra$. Fits to the mean vertical 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 T H}$, and thicknesses 
$\delta _{\infty }\sim 10^{-2}$ in terms of 
$H$, resulting in small shear Reynolds numbers 
$Re_{shear}\lesssim 1\unicode{x2013}10$ even for the largest 
$Ra$. We further analyse fluctuations of the velocity components, determine the distances of maximum mean-square fluctuations from the wall and discuss the resulting Reynolds numbers.
A quick impression of the complex boundary layer dynamics is obtained by the streamline and contour plots in figure 1 close to the bottom wall for two Rayleigh numbers. The figures indicate a prominent patchiness of the whole velocity boundary layer viewed from the top. The boundary layer is composed of coherent shear-dominated and incoherent shear-free regions. This feature becomes less prominent for the contours of the temperature field, which display an increasingly dense skeleton of thermal plume ridges over the whole plate. Here, we quantify the corresponding area fraction, condition the fluctuations on the coherent and incoherent regions and relate the incoherent regions to the large-scale patterns in the bulk, which are the turbulent superstructures of convection (Pandey et al. Reference Pandey, Scheel and Schumacher2018; Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018). Only when the velocity fluctuations are conditioned on the coherent shear-dominated regions are the mean profiles close to those observed in closed cylinders for 
$\varGamma \sim 1$ (Scheel & Schumacher Reference Scheel and Schumacher2017).
Turbulent fields at two Rayleigh numbers. Streamline snapshots of the horizontal velocity field 
$(u_x, u_y)$ and contour plots of temperature field 
$T$ at the same instant. Panels show (a,c) 
$Ra=10^8$ and (b,d) 
$Ra=10^{11}$. The entire cross-section of 
$4H\times 4H$ at 
$z=\delta _T/2$ close to the bottom wall is shown. Here, 
$H$ is the height of the convection layer and 
$\delta _T$ the thermal boundary layer thickness. The colour legend holds for both temperature field panels (c,d). In panels (a,b) examples for local coherent shear-dominated and incoherent shear-free regions are indicated.
One important point needs to be made here. While a thermal boundary layer of the standard type is indeed present, no momentum boundary layer can be easily identified, as explained subsequently. For the velocity field, it is thus more accurate to merely discuss the flow near the wall instead of the boundary layer flow, but for convenience of identification and following convention, we continue to use the term boundary layer. There is no such ambiguity for the thermal boundary layer.
2. Numerical simulations and resolution analysis
We solve the three-dimensional Boussinesq equations of RBC (Verma Reference Verma2018) by the GPU-based spectral element code nekRS (Fischer et al. Reference Fischer2022) which combines an element decomposition of the computational domain with a spectral expansion in Lagrange polynomials of each involved field along each spatial dimension on each element. The equations are given in dimensionless form by
$$\begin{gather} \frac{\partial {\boldsymbol{u}}}{\partial t} + ({\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla}) {\boldsymbol{u}} =-\boldsymbol{\nabla}p + T \hat{{\boldsymbol{z}}} + \sqrt{\frac{Pr}{Ra}} \, \nabla^2 \boldsymbol{u}, \end{gather}$$
$$\begin{gather}\frac{\partial T}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) T = \frac{1}{\sqrt{PrRa}} \, \nabla^2 T, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0. \end{gather}$$
Here, 
$\boldsymbol {u}$, 
$p$ and 
$T$ are the velocity, pressure and temperature fields, respectively. Length, velocity and temperature are expressed in units of 
$H$, 
$U_f$ and the outer temperature difference 
$\Delta T$, respectively. No-slip boundary conditions apply for the velocity field at the plates at 
$z=0$ and 
$H$. Table 1 summarizes 10 simulations, all at a Prandtl number 
$Pr=\nu /\kappa =0.7$ and aspect ratios 
$\varGamma =L/H=2$, 4 and 8, 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 15$. Furthermore, we verified that the Nusselt numbers 
$Nu_{vol}$ and 
$Nu_{wall}$, which are given by combined volume–time 
$\langle {\cdot }\rangle _{V,t}$ and area–time averages 
$\langle {\cdot }\rangle _{A,t}$,
\begin{equation} Nu_{vol}=1+\sqrt{Ra Pr}\langle u_z T\rangle_{V,t} \quad\mbox{and}\quad Nu_{wall}=\left.-\frac{\partial \langle T\rangle_{A,t}}{\partial z}\right|_{z=0,H}, \end{equation}result in practically the same values (table 1).
As a first global result, we plot in figure 2(a,b) the Nusselt and Reynolds numbers, 
$Nu=Nu_{vol}$ and 
$Re=U_{rms} \sqrt {Ra/Pr}$ vs the Rayleigh number 
$Ra$, compensated by the high-
$Ra$ scaling result of Iyer et al. (Reference Iyer, Scheel, Schumacher and Sreenivasan2020); here, 
$U_{rms}=\langle \boldsymbol {u}^2\rangle _{V,t}^{1/2}$ is the root-mean-square velocity. Also compared are our results with those of Stevens, Verzicco & Lohse (Reference Stevens, Verzicco and Lohse2010) and Scheel & Schumacher (Reference Scheel and Schumacher2017) in closed cylinders at 
$\varGamma =1/2$ and 
$\varGamma =1$, respectively. Considering the large differences in aspect ratios, the Nusselt numbers collapse fairly well and follow the same trend for 
$Ra\leq 10^{11}$ (figure 2a); they are also in agreement with previous simulations for 
$Ra\leq 10^8$ from van Reeuwijk et al. (Reference van Reeuwijk, Jonker and Hanjalić2008a). The Reynolds number shows a strong geometry dependence, as visible in figure 2(b), although they tend to the same exponent towards high 
$Ra$; this suggests an agreement in scaling exponents at high 
$Ra$, but the prefactor seems to have a complicated Rayleigh-number dependence. Statistics in all runs are obtained for equal and more than 100 free-fall times 
$T_f=H/U_f$; see table 1.
(a) Compensated turbulent heat transfer 
$Nu/Ra^{0.331}$ vs 
$Ra$. (b) Compensated turbulent momentum transfer 
$Re/Ra^{0.458}$ vs 
$Ra$. The panels contain data from the present DNS and from Scheel & Schumacher (Reference Scheel and Schumacher2017), Stevens et al. (Reference Stevens, Verzicco and Lohse2010) and Iyer et al. (Reference Iyer, Scheel, Schumacher and Sreenivasan2020). All data are compensated by the fits to the data from Iyer et al. (Reference Iyer, Scheel, Schumacher and Sreenivasan2020); see also table 1.
We verified that the resolution of the boundary layers is sufficient. Figure 3 shows the vertical profiles of the temperature fluctuations for three different spectral element grids at 
$Ra=10^{10}$ and two at 
$Ra=10^{11}$ with different polynomial orders 
$p\geq 5$. It is seen that the profiles collapse well on each other, thus also demonstrating the convergence of the spectral method. This conclusion is further reinforced for 
$Ra=10^9$ where we have 4 DNS runs at different aspect ratio and vertical grid stretching. First, the Nusselt and Reynolds numbers in table 1 differ only slightly (of the order of a tenth of a per cent) between 
$\varGamma =4$ and 8. Secondly, they are very close for two element grids with different vertical grid stretching at 
$\varGamma =4$. The fast spectral convergence of spectral element methods in comparison with lower-order finite-difference schemes has been reported recently by Zahtila et al. (Reference Zahtila, Lu, Chan and Ooi2023) in comprehensive studies for turbulent channel flows.
Resolution study for the runs at the two highest Rayleigh numbers. The near-wall profile of the temperature fluctuations is plotted. The legend indicates the Rayleigh number together with the vertical numbers of spectral elements and the polynomial order. The runs, as listed in the legend from top to bottom, have 11, 14, 16, 11 and 16 points inside the thermal boundary layer. The horizontal lines are the thermal fluctuation boundary layer thicknesses.
3. Mean profiles of temperature and velocity
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
\begin{equation} \langle u_i\rangle_{A,t}(z)=\frac{1}{A N_t}\sum_{k=1}^{N_t}\int_A u_i(x,y,z,t_k)\,{\rm d}A, \end{equation}
for 
$i=x,y$ and 
$N_t$ the number of snapshots. For all runs, 
$N_t=20$ to obtain equidistant and statistically independent realizations of the flow. Figure 4 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 even longer averaging times, which were varied from 
$\tau _{total}=400\, T_f$ (
$=20$ snapshots) to 
$1600\, T_f$ (
$=80$ snapshots). The profile converges steadily to zero, although not uniformly. There is essentially no mean flow. If we insist upon fitting the near-wall mean profiles to the two-dimensional Blasius solution, for example, the result is shown in panel (b). 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 4. Recall that the Blasius solution 
$f(\eta )=\psi (x,z)/\sqrt {x U_{\infty } \nu }$ and 
$\eta =z/\delta (x)$ (Schlichting & Gersten Reference Schlichting and Gersten2016), where 
$\psi (x,z)$ is the streamfunction and 
$u_x=\partial \psi /\partial z$ and 
$u_z=-\partial \psi /\partial x$. 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 for different time intervals in (b). Recall that, at distance 
$\eta =5$, the Blasius profile reaches a streamwise velocity magnitude of 
$f^{\prime }(\eta )=0.99 U_{\infty }$. In this case, panel (a) shows that the maximum velocity reached is of the order 1 % or less of the free-fall velocity, which is the only characteristic velocity of the flow. As shown in panel (c), no clear trend of the velocity 
$U^{x,y}_{\infty }$ with averaging time is detectable, and the magnitude is between 
${O}(10^{-3})$ and 
${O}(10^{-2})$. In table 2, the results for all Rayleigh numbers and both horizontal components are listed.
(a) Planar-averaged profiles of the horizontal velocity component 
$u_x$ for different averaging intervals for 
$Ra = 10^9$. Notice the very small values of 
$U_{\infty }$ and the non-monotonic convergence of the profiles. (b) Match of the near-wall profiles from panel (a) to the Blasius profile 
$f^{\prime }(\eta )$ (dashed line). (c) Variation of the free-stream velocity 
$U_{\infty }$ and boundary layer height 
$\delta _{\infty }$ obtained from the fits to the Blasius profile vs averaging time at 
$Ra=10^9$. (d) Time-averaged profiles 
$\langle u_x\rangle _{A,t}$ for all 7 Rayleigh numbers. (e) Time-averaged profiles 
$\langle u_y\rangle _{A,t}$ for all 7 Rayleigh numbers. (f) Time-averaged profiles 
$\langle u_x\rangle _{A,t}$ for Rayleigh number 
$Ra=10^9$ and aspect ratios 
$\varGamma =2$, 4 and 8. These are the three runs which are indicated by a dagger in table 1.
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}_{shear}=U^{x,y}_{\infty }\delta ^{x,y}_{\infty }/\nu$ are also given. All runs are at 
$\varGamma =4$.
The boundary layer thickness parameters vary 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_{shear}=U_{\infty }\delta _{\infty }/\nu$, which are found to be very small for all cases. Panels (d,e) of figure 4 show the mean profiles for all 7 simulation runs and for both horizontal velocity components. They underline the very small mean-flow amplitudes for all Rayleigh numbers in this series. Panel (f) of the figure compares three runs at 
$Ra=10^9$, which are indicated with a dagger symbol in table 1, at aspect ratios 
$\varGamma =2$, 4 and 8. Again, the mean-flow amplitudes 
$\langle u_x\rangle _{A,t}(z)$ are comparable and very small such that an aspect-ratio dependence for this result can be excluded when periodic boundary conditions in the horizontal directions are applied. The non-monotonic behaviour is simply a reflection of the long averaging times required for convergence in convection studies.
4. Fluctuation profiles of temperature and velocity
The mean vertical profiles of the root-mean-square velocities are given by
\begin{equation} U_{rms}^h(z)=\sqrt{\langle u_x^2+u_y^2\rangle_{A,t}}\quad\mbox{and}\quad U_{rms}(z)=\sqrt{\langle u_x^2+u_y^2+u_z^2\rangle_{A,t}}, \end{equation}
where we distinguish between horizontal and full profiles. The fluctuation profiles for the velocity fields are obtained from the full components since the means are so small. We have verified that the differences in the procedure produce only very small changes. Figure 5 summarizes mean profiles for 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 similarly obtained by 
$T_{rms}(z)=\langle \theta ^2\rangle ^{1/2}_{A,t}$ with 
$\theta (\boldsymbol {x},t)=T(\boldsymbol {x},t)-\langle T\rangle _{A,t}(z)$. The corresponding 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/(2 Nu)$ is slightly larger than the distance from the wall of the maximum of the temperature fluctuation profile, which we term the thermal fluctuation boundary layer thickness or, for short, thermal fluctuation thickness, 
$\delta _{T,{rms}}$. Increasingly larger are distances from the wall to the maxima of the velocity fluctuation profiles, obeying a ratio of 
$\delta _{U,{rms}}/\delta _{T,{rms}}\approx 2$ for 
$Ra=10^5$ up to approximately 13 for 
$Ra=10^{11}$; see panel (d). The corresponding thicknesses are termed velocity fluctuation thickness.
Mean profiles are compared with 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_{rms}(z)$ in (b), root-mean-square profile with respect to the two horizontal velocity components, 
$U^h_{rms}(z)$, in (c) and with respect to all three components, 
$U_{rms}(z)$, in (d). The corresponding height scales of these profiles are indicated by dashed horizontal lines in the figures and listed in table 3. Note that, for the two lowest 
$Ra$, they exceed the displayed plot range of the figure in panel (d).
Differently defined boundary layer thicknesses of temperature and velocity according to the vertical mean profiles plotted in figure 5. The standard deviation accompanies each mean value.
Furthermore, we repeated the fluctuation profile analysis for the dependence on different aspect ratios, one on each side of 4, at 
$Ra=10^9$. These runs are indicated again by a dagger symbol in table 3, where we collect the corresponding thickness scales. The corresponding profiles are shown in figure 6. While the temperature profiles collapse close to the walls, thus displaying no sensitivity with respect to the aspect ratio 
$\varGamma$ in this range, the velocity profiles are affected by 
$\varGamma$. However, the resulting velocity fluctuation thicknesses are found to agree well for 
$\varGamma =4$ and 8 (by 1.5 % or less), in terms of both horizontal and full velocity fluctuations. The finding supports our considered view that 
$\varGamma \geq 4$ is sufficient to obtain horizontal homogeneity for the statistics already introduced.
Mean temperature profiles compared with each other for different aspect ratios at 
$Ra=10^9$. Similar to figure 5, we show the mean profile of temperature 
$\langle T\rangle _{A,t}(z)$ in panel (a), root-mean-square profile of temperature 
$T_{rms}(z)$ in (b), root-mean-square profile with respect to the horizontal velocity components, 
$U^h_{rms}(z)$ in (c) and with respect to all three components 
$U_{rms}(z)$ in (d). The corresponding height scales of these profiles are indicated by dashed horizontal lines and are also listed in table 3.
5. Scaling of combined volume–time-averaged fluctuation with Rayleigh number
Figure 7 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^2H$ and time, e.g. 
$u_{x,{rms}}=\langle u_x^2\rangle ^{1/2}_{V,t}$. The quantity 
$U_{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_{rms}\sim Ra^{-\beta }$, which is found to be 
$\beta =0.119$ 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. (Reference Castaing, Gunaratne, Heslot, Kadanoff, Libchaber, Thomae, Wu, Zaleski and Zanetti1989), Niemela et al. (Reference Niemela, Skrbek, Sreenivasan and Donnelly2000) and Wu & Libchaber (Reference Wu and Libchaber1992) 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.141$ 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.
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. (Reference Iyer, Scheel, Schumacher and Sreenivasan2020) with separate fits. The corresponding temperature data are obtained there for the cylindrical cell interior, 
$r\leq 0.03 H$ and 
$Ra\geq 10^9$.
In panels (b–d) of figure 7, we added data from the DNS of Iyer et al. (Reference Iyer, Scheel, Schumacher and Sreenivasan2020) for comparison, which were obtained in a closed slender cylindrical cell of aspect ratio 
$\varGamma =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 geometrically constrained convection. This finding of nearly the same scaling exponents suggests a robust geometry-independent trend of all thermal fluctuations with respect to the Rayleigh number. Geometry-specific aspects mostly affect the prefactor.
6. Decomposition into coherent and incoherent boundary layer regions
The orientation of the boundary layer flow varies strongly, as shown in figures 8(a) and 8(c), where we plot the orientation angle of the horizontal velocity 
$\varphi =\arctan (u_y/u_x) \in [-{\rm \pi},{\rm \pi} ]$ for a snapshot at 
$Ra=10^{10}$ at 
$z=\delta _T$ and 
$z=H-\delta _T$, respectively. At both these heights, we cover the horizontal plane into 
$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 {\boldsymbol {u}}_h(A_i)$ in each of 
$A_i$ and decompose the cross-section into coherent and incoherent boundary layer regions for 
$|\bar {\boldsymbol {u}}_h(A_i)| > U^h_{rms}(\delta _T)$ and 
$|\bar {\boldsymbol {u}}_h(A_i)|\leq U^h_{rms}(\delta _T)$, respectively. Similar decompositions have been applied to analyse the spatio-temporal intermittency of the transition to turbulence in shear flow turbulence in extended domains, see e.g. Hof (Reference Hof2022). Panels (a,c) of figure 8 show that coherent shear-dominated patches are separated by incoherent flow regions (in grey). See also figure 1. The superposed streamlines indicate the different flow orientations of the shear-dominated regions.
Boundary layer flow structure for 
$Ra=10^{10}$ at the top and bottom walls. (a) Decomposition of a snapshot into coherent and incoherent (in grey) flow regions at 
$z=\delta _T$. For the coherent shear-dominated regions, we overlay horizontal streamlines and indicate their local flow direction. (b) Corresponding temperature field in plane at 
$z=0.1\approx 25 \delta _T$. Grey area corresponds to 
$T>T_0=0.5$. (c) Same snapshot with the decomposition similar to (b) for 
$z=H-\delta _T$ close to the top wall. (d) Corresponding temperature field at 
$z=0.9\approx 25 \delta _T$. Grey area corresponds now to 
$T< T_0$. Note that all plots are shown at the coarse resolution of 
$100\times 100$ square boxes (which were used to calculate the local mean-flow magnitude and orientation) and not at the original spectral resolution of the DNS run.
Panels (b,d) of the same figure show the corresponding snapshots of the temperature field 
$T$ at 
$z=0.1$ above the bottom and 
$z=0.9$ below the top, which are distances of 25 
$\delta _T$ away from the walls at 
$Ra=10^{10}$. It is clearly seen that the hotter regions at 
$z=0.1$ and the colder regions at 
$z=0.9$, both of which are displayed in grey, coincide fairly well with the incoherent flow regions at the edge of the thermal boundary layer. We can define overlap factors 
$0\leq \tilde {O}\leq 1$ by
\begin{equation} \tilde{O}_{bot}=\frac{A_{incoh} \bigcup A_{T>T_0}}{\max(A_{incoh},A_{T>T_0})} \quad\mbox{and}\quad \tilde{O}_{top}=\frac{A_{incoh} \bigcup A_{T< T_0}}{\max(A_{incoh},A_{T< T_0})}, \end{equation}
with 
$T_0=0.5$. Here, we find mean overlaps of 
$\langle \tilde {O}_{bot}\rangle =0.60$ and 
$\langle \tilde {O}_{top}\rangle =0.63$, where the average is taken over the snapshots. The physical interpretation is as follows: the incoherent regions correspond to dominant hotter upwelling (colder downwelling) motions. These regions occur outside shear-dominated patches where the thermal plumes merge successively with growing distances from the walls. As one approaches the mid-plane of the convection cell, they tend to form the turbulent superstructure pattern of convection. We have determined that the area fraction of the incoherent regions remains nearly constant at approximately 60 % of 
$A$ for the whole Rayleigh-number range. The insensitivity of the volume fractions with respect to the Rayleigh number suggests that this skeleton of upwelling (downwelling) incoherent regions could be a relic from the weakly nonlinear regime of convection at much lower Rayleigh numbers, which itself arises from the onset of convection by a linear primary instability, filling the whole domain with convection rolls.
We have varied the threshold for this analysis from 
$0.5 U^h_{rms}(\delta _T)$ to 
$2 U^h_{rms}(\delta _T)$. While the incoherent fractions do depend on the threshold when its variations are large, they are practically independent of the Rayleigh number even for such large variations stated above. This supports our choice of 
$U^h_{rms}(\delta _T)$ as a physically meaningful threshold.
We can now return to the fluctuation analysis which is conditioned on coherent and incoherent regions in the following. Figure 9 replots the root-mean-square profiles of full and horizontal velocity and temperature profiles for Rayleigh numbers 
$Ra=10^8$, 
$10^9$ and 
$10^{10}$. We have chosen these three Rayleigh numbers of our series to provide a one-to-one comparison with DNS data in a closed cylindrical cell at 
$\varGamma =1$ of Scheel & Schumacher (Reference Scheel and Schumacher2017). They are also shown in the figure. Vertical profiles, which have been taken over the full cross-section (denoted as case G4 in the following), are shown in the left column of figure 9. Profiles conditioned on shear-dominated regions are displayed in the middle column (case G4C), while those for the cylindrical cross-section of the closed container (case G1) are shown in the right column. From the bottom row of the figure, it is clear that the temperature profiles of G4, G4C and G1 for all three 
$Ra$ agree. This suggests that the temperature boundary layers are alike in all cases. This is different for the velocity field, for which the horizontal velocity fluctuations (displayed in top row) show a clear trend. The thickness scale decreases from G4 to G4C and even more from G4C to G1. The close agreement of G4C and G1 clearly supports the dominance of shearing motion in the boundary regions in closed cylindrical cells, imposed by the prominent LSC. It is in line with a reduced fluctuation thickness. For fluctuations with respect to the full velocity field, we do detect a decrease of the thickness from G4 to G4C, but not from G4C to G1 for the two lower 
$Ra$. We suspect that this might be caused by prominent coherent up- and downwelling motions at the sidewalls for the lower 
$Ra$ which effectively enhance the thickness (Schumacher & Scheel Reference Schumacher and Scheel2016).
Vertical profiles of the velocity and temperature fluctuations for Rayleigh numbers 
$Ra=10^8$, 
$10^9$ and 
$10^{10}$. The figure compares averages over the full cross-section at 
$\varGamma =4$ (a,d,g) with averages over the coherent shear-dominated regions (b,e,h) and with averages over the full circular cross-section in a closed cylindrical cell at 
$\varGamma =1$ (c,f,i) from Scheel & Schumacher (Reference Scheel and Schumacher2017). Panels (a–c) are for 
$U^h_{rms}$, panels (d–f) for 
$U_{rms}$ and panels (g–i) for 
$T_{rms}$. The dashed lines in each panel mark the location of the first local maximum away from the wall, and indicate the corresponding fluctuation thickness. The colour coding, which is indicated in panel (c), holds for all panels.
Unlike G1 and G4C, for which the velocity and temperature boundary layers have comparable thicknesses, the G4 case shows that the velocity boundary layer is much thicker than the temperature thickness, suggesting a different mechanism in G4. We recall that the notion of the velocity boundary layer is only nominal in the sense that they are based on fluctuation profiles and the mean velocity variation within that region is quite small (see figure 4).
7. Final discussion
Our DNS of the turbulent Rayleigh–Bénard convection encompass a Cartesian domain with 
$\varGamma =4$, with no-slip horizontal walls and periodic boundary conditions for the side faces. These simulations of up to 
$Ra=10^{11}$ are aimed at approaching the original canonical case of a plane convection layer between a pair of infinitely extended rigid plates. We demonstrated that a standard mean-flow profile (obtained by combining averages with respect to time and the entire horizontal cross-sectional plane) have very small magnitudes, and that efforts to match them to laminar boundary layer profiles produced no conclusive results. To the extent that we can define the boundary layers, they give very small shear Reynolds numbers (see below). In the long-time limit, which we have followed for 1600 
$T_f$ at 
$Ra=10^9$, the velocity mean profiles have to converge to 
$\mbox {lim}_{t\to \infty }\langle u_i\rangle _{A,t}(z)\to 0$ due to statistical homogeneity in 
$x$ and 
$y$. The simulations by Hartlep et al. (Reference Hartlep, Tilgner and Busse2003) (for 
$Ra\leq 1\times 10^7$) also showed that the mean flow contained very little kinetic energy, but De et al. (Reference De, Eswaran and Mishra2018) found a long-time periodicity in the mean flow for low Rayleigh numbers, 
$Ra\leq 2\times 10^6$.
Rather than having a mean-flow profile with small velocity fluctuations, we are faced with small mean velocity amplitudes in the presence of velocity fluctuations that are up to 2 orders of magnitude larger when the statistics are taken over finite time intervals 
$\tau _{\rm total}$, as seen from comparisons of table 2 with the data in figure 5. This central result also holds when the aspect ratio of the simulation is varied. It is our view that fluctuations will be relevant for all configurations which includes closed cells of 
$\varGamma \lesssim 1$, see e.g. figure 9. But their relevance is strongest in the statistically homogeneous plane layer with periodic boundary conditions in the horizontal direction – the configuration that comes closest to the original physical problem of turbulent convection (Spiegel Reference Spiegel1962), as relevant for most geo- and astrophysical applications.
We also showed that the corresponding shear Reynolds numbers, which are based on mean-flow quantities for a finite averaging time, remain very small because the characteristic velocities 
$U_{\infty }$ are small. The strong fluctuations cause the fluctuation thicknesses of temperature and velocity (defined as the near-wall maxima of the root-mean-square profiles) to differ by an order of magnitude for the highest Rayleigh numbers, as summarized in table 3. This difference increases with Rayleigh number (although the Prandtl number is held fixed at order unity); it becomes particularly pronounced for 
$Ra\geq 10^9$, a range beyond which previous larger-aspect-ratio DNS studies rarely advanced.
Furthermore, our analysis revealed that the velocity boundary region in the present configuration is a carpet of differently oriented time-dependent and shear-dominated (coherent) regions interspersed by regions of incoherent flow. The latter regions occupy approximately 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 incoherent regions in the present flow can be as large as the entire cross-section of a cylindrical cell at 
$\varGamma \sim 1$. The coherent regions are the near-wall footprint of the circulation rolls which form the large-scale turbulent superstructure pattern (Pandey et al. Reference Pandey, Scheel and Schumacher2018; Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018). They change their orientation continually and thus result in a net zero mean flow, as stated above.
Finally, we showed that the velocity fluctuation thicknesses decrease when they are conditioned on shear-dominated patches in the near-wall region. They are then closer to those scales which are obtained in turbulent convection in closed cylindrical cells of aspect ratio 
$\varGamma \lesssim 1$. The geometry of the closed cell enforces a LSC which is mostly shear dominated in the vicinity of the walls, as already shown in Schumacher et al. (Reference Schumacher, Bandaru, Pandey and Scheel2016), where the time dependence of the orientation has been eliminated. This causes smaller velocity thickness scales that are, however, still larger than the thermal boundary layer thickness at 
$Pr\sim 1$. Nevertheless, the fluctuation thickness is the consistently definable velocity boundary layer scale for the present turbulent convection flow.
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. Differently from wall-bounded shear flows, we detect velocity fluctuations everywhere, even though at different strengths; see again figure 9. Furthermore, we do not observe a Rayleigh-number dependence of the ratio of coherent (‘laminar’) to incoherent (‘turbulent’) regions. The time scales, at which these complex spatio-temporal patterns change, become increasingly shorter with increasing Rayleigh number. A variation of the threshold for the decomposition into coherent and incoherent boundary regions practically did not affect this Rayleigh-number independence.
The spatio-temporal heterogeneity of the velocity boundary layer, which we 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 (Reference Howard and Goertler1966) and for a detailed boundary layer model with plume formation, Theerthan & Arakeri (Reference Theerthan and Arakeri1998). However, figure 10 shows a power-law fit of 
$Ra_{\delta,{rms}}=A Ra^\gamma$ with a very small exponent 
$\gamma =0.077$ and 
$A\approx 41$. The resulting 
$Ra_{\delta,{rms}}$ are, by at least a factor of 4, smaller than Howard's critical Rayleigh number of 
$Ra_{\delta }\sim {O}(10^3)$. Additional data from experiments show that, even at 
$Ra\sim 10^{17}$, a thickness-based Rayleigh number 
$Ra_{\delta }$ barely reaches a value of 
$10^3$. This challenges the original marginal stability concept. This question and higher Rayleigh-number simulation in the present configuration form the subject of further study.
Boundary layer thickness-based Rayleigh number 
$Ra_{\delta }$ vs Rayleigh number 
$Ra$. We have collected for comparison data from other simulations including Scheel & Schumacher (Reference Scheel and Schumacher2017) and Iyer et al. (Reference Iyer, Scheel, Schumacher and Sreenivasan2020); also plotted are experimental data from Fleischer & Goldstein (Reference Fleischer and Goldstein2002) and Niemela et al. (Reference Niemela, Skrbek, Sreenivasan and Donnelly2000). Note that all listed data from previous sources use 
$\delta _T=1/(2 Nu)$ to calculate the boundary layer thickness-based Rayleigh number 
$Ra_{\delta }$. Only the present data use 
$\delta _{T, {rms}}$ for the determination of 
$Ra_{\delta }$. The corresponding Rayleigh numbers 
$Ra_{\delta }=Ra_{\delta,{rms}}$ are found in table 1.
Acknowledgements
The authors thank J. Arakeri, B. Hof, H. Nagib and P. Schmid for helpful discussions.
Funding
The work of R.J.S. 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 the project nonbou by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS Booster at Jülich Supercomputing Centre (JSC).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available on reasonable request.
Author contributions
All authors designed the research. R.J.S., M.B. and J.S. carried out the supercomputer simulations. R.J.S., J.D.S and J.S. analysed the simulation data and generated the figures. All authors wrote the manuscript.
Open access
