2. Numerical method

The data used in this manuscript have previously been presented by Stevens et al. (Reference Stevens, Blass, Zhu, Verzicco and Lohse2018) and Krug et al. (Reference Krug, Lohse and Stevens2020). A summary of the most relevant quantities for this study can be found in table 1; note that, there and elsewhere, we use the free-fall velocity $V_{ff} = \sqrt {g\beta H\varDelta }$ as a reference scale. In the following, we briefly report details on the numerical method for completeness. We carried out periodic RB simulations by numerically solving the three-dimensional incompressible Navier–Stokes equations within the Boussinesq approximation. They read:

(2.1)\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} =-\boldsymbol{\nabla} P + \nu \nabla^2\boldsymbol{u}+\beta g \theta \hat{z}, \end{gather}

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

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

Here, $\boldsymbol {u}$ is the velocity vector, $\theta$ the temperature and the kinematic pressure is denoted by $P$. The coordinate system is oriented such that the unit vector $\hat {z}$ points up in the wall-normal direction, while the horizontal directions are denoted by $x$ and $y$. We solve (2.1)–(2.3) using AFiD, the second-order finite difference code developed by Verzicco and coworkers (Verzicco & Orlandi Reference Verzicco and Orlandi1996; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015). We use periodic boundary conditions and a uniform mesh in the horizontal direction and a clipped Chebyshev-type clustering towards the plates in the wall-normal direction. For validations of the code against other experimental and simulation data in the context of RB we refer to Verzicco & Orlandi (Reference Verzicco and Orlandi1996), Verzicco & Camussi (Reference Verzicco and Camussi1997, Reference Verzicco and Camussi2003), Stevens, Verzicco & Lohse (Reference Stevens, Verzicco and Lohse2010) and Kooij et al. (Reference Kooij, Botchev, Frederix, Geurts, Horn, Lohse, van der Poel, Shishkina, Stevens and Verzicco2018).

Table 1. Data from Stevens et al. (Reference Stevens, Blass, Zhu, Verzicco and Lohse2018) and Krug et al. (Reference Krug, Lohse and Stevens2020) for the global Nusselt number, the grid resolution $(N_x,N_y,N_z)$ in the streamwise, spanwise and wall-normal directions, the location of the coherence spectrum peak $\hat {l}$, the root mean square velocity $v_{RMS}=\sqrt {\langle v_x^2 + v_y^2 + w^2 \rangle _V}$ non-dimensionalised with the free-fall velocity $V_{ff}=\sqrt {\beta g H \varDelta }$, the estimated thermal BL thickness $\lambda _{\theta }^*/H=1/(2Nu)$ and the amount of non-dimensional time units used for our statistical analysis $t V_{ff}/H$.

The aspect ratio of our domain is $\varGamma =L/H=32$, where $L$ is the length of the two horizontal directions of the periodic domain. The used numerical resolution ensures that all important flow scales are properly resolved (Shishkina et al. Reference Shishkina, Stevens, Grossmann and Lohse2010; Stevens et al. Reference Stevens, Verzicco and Lohse2010). We note that the grid resolution at $Ra=10^9$ still has $11$ grid points in the thermal and kinetic boundary layer, while the criteria by Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010) state that 8 grid points are sufficient in this case. In the Appendix we give further details on the simulations and show that the average of the horizontal velocity components is almost zero.

In this manuscript, we define the decomposition of instantaneous quantities into their mean and fluctuations such that $\psi (x,y,z,t)= \varPsi (z) + \psi ' (x,y,z,t)$, where $\varPsi = \langle \psi (x,y,z,t) \rangle _{x,y,t}$ is the temporal and horizontal average over the whole domain and $\psi '$ denotes the fluctuations with respect to this mean.

3. Conditional averaging

Extracting features of the LSC in large aspect ratio cells poses a significant challenge. The reason is that there are multiple large-scale structures of varying sizes, orientation and inter-connectivity at any given time. It is therefore not possible to extract properties of the LSC by using methods that rely on tracking a single or a fixed small number of convection cells, such as those applied successfully in analysing the flow in small (Sun, Cheung & Xia Reference Sun, Cheung and Xia2008; Wagner et al. Reference Wagner, Shishkina and Wagner2012) to intermediate (van Reeuwijk, Jonker & Hanjalić Reference van Reeuwijk, Jonker and Hanjalić2008) aspect ratio domains. To overcome this issue, we use a conditional averaging technique developed in Berghout et al. (Reference Berghout, Baars and Krug2020), where this framework was employed to study the modulation of small-scale turbulence by the large flow scales. This approach is based on the observation of Krug et al. (Reference Krug, Lohse and Stevens2020) that the premultiplied temperature power spectra $k\varPhi _{\theta \theta }$ (shown in figure 1a,c,e) is dominated by two very distinct contributions. One is due to the ‘superstructures’ whose size (relative to $H$) increases with increasing $Ra$ and typically corresponds to wavenumbers $kH \approx 1\text {--}1.5$. The other contribution relates to a ‘near-wall peak’ with significantly smaller structures whose size scales with the thickness of the BL (Krug et al. Reference Krug, Lohse and Stevens2020). This implies that this peak shifts to larger $k$ (scaled with $H$) as the BLs get thinner at higher $Ra$. Hence, there is a clear spectral gap between superstructures and small-scale turbulence, which widens with increasing $Ra$, as can readily be seen from figure 1(a,c,e). This figure also demonstrates that a spectral cutoff $k_{{cut}} = 2/H$ is a good choice to separate superstructure contributions from the other scales over the full $Ra$ range $10^5\leq Ra\leq 10^9$ considered here.

Figure 1. (a,c,e) Premultiplied temperature power spectra $k \varPhi _{\theta \theta }$ for $Ra=10^5;10^7;10^9$. The blue line indicates the cutoff wavenumber $k_{{cut}} = 2/H$ used for the low-pass filtering. The dashed black lines indicate alternative cutoffs ($k_{{cut}} = 1.8/H$ and $k_{{cut}} = 2.5/H$) considered in panel (d). The white plusses are located at $k=0.57/\lambda _{\theta }^*$ and $z=0.85 \lambda _{\theta }^*$ (with $\lambda _{\theta }^* =H/(2Nu)$) in all cases, which corresponds to the location of the inner peak (Krug et al. Reference Krug, Lohse and Stevens2020) (b) Coherence spectra of temperature and wall-normal velocity at mid-height, figure adopted from Krug et al. (Reference Krug, Lohse and Stevens2020). The black line illustrates the choice of $k_{{cut}} = 2/H$ and the legend of figure 4(a) applies for the $Ra$ trend. (d) Snapshot of temperature fluctuations for $Ra=10^7$ at mid-height. The black lines show contours of $\theta _L ' =0$ evaluated for different choices of $k_{{cut}}$.

The choice for $k_{{cut}} = 2/H$ is further supported by considering the spectral coherence

(3.1)\begin{equation} \gamma^{2}_{\theta w}(k)=\frac{|\varPhi_{\theta w}(k)|^2}{\varPhi_{\theta \theta}(k)\varPhi_{w w}(k)}, \end{equation}

where $\varPhi _{w w}$ and $\varPhi _{\theta w}$ are the velocity power spectrum and the co-spectrum of $\theta$ and $w$, respectively. The coherence $\gamma ^2$ may be interpreted as a measure of the correlation per scale. The results at $z = 0.5H$ in figure 1(b) indicate that there is an almost perfect correlation between $\theta '$ and $w'$ at the superstructure scale. At larger wavenumbers, this correlation is seen to drop significantly. For the lower $Ra$ values, the coherence rises again at the very small scales. However, almost no energy resides at the scales corresponding to the high-wavenumber peak in $\gamma ^{2}_{\theta w}$ (see figure 1(a,c,e) and for a more detailed discussion Krug et al. (Reference Krug, Lohse and Stevens2020)), such that the coherence there is of little practical relevance. The threshold $k_{{cut}} = 2/H$ effectively delimits the large-scale peak in $\gamma ^2_{\theta w}$ towards larger $k$ for all $Ra$ considered, such that this value indeed appears to be a solid choice to distinguish the large-scale convection rolls from the remaining turbulence. To confirm this, we overlay a snapshot of $\theta '$ with zero crossings of the low-pass filtered signal (with cutoff wavenumber $k_{{cut}}$) $\theta '_L$ in figure 1(d). These contours reliably trace the visible structures in the temperature field. Furthermore, it becomes clear that slightly different choices for $k_{{cut}}$ do not influence the contours significantly. This is consistent with the fact that only limited energy resides at the scales around $k\approx 2/H$, such that $\theta '_L$ only changes minimally when $k_{{cut}}$ is varied within that range. In the following, we adopt $k_{{cut}} = 2/H$ to obtain $\theta '_L$ except when we study the effect of the choice for $k_{{cut}}$.

We use $\theta '_L$ evaluated at mid-height to map the horizontal field onto a new horizontal coordinate $d$. To obtain this coordinate, first the distance $d^*$ to the nearest zero crossing in $\theta '_L$ is determined for each point in the plane. This can be achieved efficiently using a nearest-neighbour search. Then the sign of $d$ is determined by the sign of $\theta '_L$, such that $d$ is given by

(3.2)\begin{equation} d=\textrm{sgn}(\theta '_L){d^*}. \end{equation}

All results presented here are with reference to the lower hot plate. Hence $d<0$ and $d>0$ correspond to plume impacting and plume emitting regions, respectively. The averaging procedure is illustrated in figure 2(a,b). Another important aspect is a suitable decomposition of the horizontal velocity component $v$. Figure 2(c) shows how we decompose $v$ into one component ($v_p$) parallel the local gradient $\boldsymbol {\nabla } d$, and another component ($v_n$) normal to it. This ensures that $v_p$ is oriented normal to the zero crossings in $\theta '_L$ for small $|d|$, where the wind is strongest. However, at larger $|d|$, the orientation may vary from a simple interface normal, which accounts for curvature in the contours. It should be noted that the $d$-field is determined at mid-height and consequently applied to determine the conditional average at all $z$-positions. This is justified since Krug et al. (Reference Krug, Lohse and Stevens2020) showed that there is a strong spatial coherence of the large scales in the vertical direction. Therefore, the resulting zero contours would almost be congruent if $\theta '_L$ was evaluated at other heights. The time-averaged conditional average is obtained by averaging over points of constant $d$, while we make use of the symmetry around the mid-plane to increase the statistical convergence. Mathematically, the conditioned averaging results in a triple decomposition according to $\psi (x,y,z,t) = \varPsi (z) + \bar {\psi }(z,d)+ \tilde {\psi }(x,y,z,t)$, where the overline indicates conditional and temporal averaging. As bin size of the $d$-array we have used the horizontal grid spacing $dx=\varGamma /N_x$.

Figure 2. Illustration of the conditional averaging method based on simulation data for $Ra=10^7$. $(a)$ Temperature fluctuation field at mid-height and corresponding distance field (right). The black lines correspond to the zero crossings $\theta _L' =0$ relative to which the distance $d^*$ is defined (see blow-up in panel b). Note that by definition isolines $\theta _L '=0$ correspond to contours of $d=0$ in the distance field. $(b)$ Illustration of the distance definition; for every point $d^*$ is equal to the radius of the smallest circle around that point which touches a $\theta _L' =0$ contour. $(c)$ Illustration of the decomposition of the horizontal velocity $v$, here at boundary layer height, into the parallel $v_p$ and the normal $v_n$ component to the gradient vector $d$. The colour scheme in (b,c) indicates the $d$-field as in (a).

Applying the outlined method to our RB dataset results in a representative large-scale structure like the one depicted in figure 3 for $Ra= 10^7$. In general, we find $\bar {\theta }<0$ with predominantly downward flow for $d<0$, while lateral flow towards increasing $d$ dominates in the vicinity of $d = 0$. In the plume emitting region $d>0$ the conditioned temperature $\bar {\theta }$ is positive and the flow upward. In interpreting the results it is important to keep in mind that the averaging is ‘sharpest’ close to the conditioning location ($d = 0$) and ‘smears out’ towards larger $|d|$ as the size of individual structures varies. We normalise $d$ with $\hat {l}$ to enable a comparison of results across $Ra$. Based on the location of the peak in $\gamma ^2$, Krug et al. (Reference Krug, Lohse and Stevens2020) found that the superstructure size is $\hat {l}=5.9 H$ at $Ra= 10^7$. As indicated, the conditionally averaged flow field in figure 3 corresponds to approximately half this size.

Figure 3. Contour plot of the conditionally averaged temperature $\bar {\theta }/\varDelta$ for $Ra=10^7$. The arrows show $\bar {w}/V_{ff}$ and $\bar {v}_p/V_{ff}$ and are plotted every 24 and every 6 data points along $d$ and $z$, respectively. The white line is the streamline which passes through $z^*/H$ at $d=0$.

We present the probability density function (PDF) of the distance parameter $d$ in figure 4(a). The data collapse to a reasonable degree, indicating that there are no significant differences in how the LSC structures vary in time and space across the considered range of $Ra$. Visible deviations are at least in part related also to uncertainties in determining $\hat {l}$ via fitting the peak of the $\gamma ^2$-curve.

Figure 4. (a) PDF of the normalised distance parameter $d/ \skew2\hat {l}$ using all available snapshots. (b) Sample velocity profile, locally determined in $(x,y)$, to illustrate the slope method ($\lambda$) and the level method ($\ell$) used to determine the instantaneous BL thicknesses.

The LSC is carried by $v_p$, which is also supported by the fact that the velocity component normal to the gradient $\boldsymbol {\nabla } d$ averages to zero, i.e. $\bar {v}_n \approx 0$, for all $d$. The determination of the viscous BL thickness is therefore based on $v_p$ only. We use the ‘slope method’ to determine the viscous ($\lambda _u$) and thermal ($\lambda _{\theta }$) BL thickness. Both are determined locally in space and time and are based on instantaneous wall-normal profiles of $\theta$ and $v_p$, respectively. As sketched in figure 4(b), $\lambda$ is given by the location at which linear extrapolation using the wall gradient reaches the level of the respective quantity. Here the ‘level’ (e.g. $u_l(x,y) =\max _{z \in I}(v_p(x,y,z))$ for velocity) is defined as the local maximum within a search interval $I$ above the plate. In agreement with Wagner et al. (Reference Wagner, Shishkina and Wagner2012) we find that the results for both thermal and viscous BL do not significantly depend on the search region when it is larger than ${I=}4 \lambda _{\theta }^*$. Therefore, we have adopted this search region in all our analyses.

In figure 5(a) we present the conditionally averaged temperature $\bar {\theta }$ as a function of $z/H$ at three different locations of $d/\hat {l}$. Consistent with the conditioning on zero crossings in $\theta '_L =0$, we find that $\bar {\theta }\approx 0$ for all $z$ at $d = 0$. In the plume impacting ($d/\hat {l} = -0.25$) and emitting ($d/\hat {l} = 0.25$) regions, $\bar {\theta }$ is respectively negative and positive throughout. On both sides, $\bar {\theta }$ attains nearly constant values in the bulk, the magnitude of which is decreasing significantly with increasing $Ra$.

Figure 5. $(a)$ Conditioned temperature $\bar {\theta }/\varDelta$ at $d=0$ and in the plume impacting $(d/\hat {l}=-0.25)$ and in the plume emitting region $(d/\hat {l}=0.25)$ for various $Ra$, see legend in $(c)$. $(b)$ Wind velocity $\bar {v}_p/V_{ff}$ at $d=0$ versus $z/H$ at the same $Ra$. The inset shows the sensitivity of the results to different choices of $k_{{cut}}$ in the range $1.8\leq k_{{cut}} H \leq 2.5$ (same range used in figure 1) for $Ra=10^7$. $(c)$ Mean wind velocity at $d=0$ normalised by its maximum value for various $Ra$ (see legend). The dashed black lines in $(c)$ represent experimental data from Sun et al. (Reference Sun, Cheung and Xia2008) at $\varGamma =1$ for $Ra=1.25 \times 10^9$ (short dash) and $Ra=1.07 \times 10^{10}$ (long dash) and the dotted black line represents the Prandtl–Blasius profile.

Profiles for the mean wind velocity $\bar {v}_p(z)$ at $d = 0$ are shown in figure 5(b,c). These figures show that the viscous BL becomes thinner with increasing $Ra$, while the decay from the velocity maximum to $0$ at $z/H = 0.5$ is almost linear for all cases. We note that of all presented results the wind profile is most sensitive to the choice of the threshold $k_{{cut}}$. The reason is that the obtained wind profile depends on both the contour location and orientation. To provide a sense for the variations associated with the choice of $k_{{cut}}$, we compare the present result at $Ra = 10^7$ to what is obtained using alternative choices ($k_{{cut}} = 1.8/H$ and $k_{{cut}} = 2.5/H$) in the inset of figure 5(b). This plot shows that results within the BL are virtually insensitive to the choice of $k_{{cut}}$ while the differences in the bulk consistently remain below 5 %. In figure 5(c) we re-plot the data from figure 5(b) normalised with the BL thickness $\bar {\lambda }_u(d=0)$ and the velocity maximum $\bar {v}_p^{{max}}$. The figure shows that the velocity profiles for the different $Ra$ collapse reasonably well for $z\lessapprox \bar {\lambda }_u$. A comparison to the experimental data by Sun et al. (Reference Sun, Cheung and Xia2008), which were recorded in the centre of a slender box with $\varGamma =1$ and $Pr = 4.3$, reveals that, although the overall shape of the profiles is similar, there are considerable differences in the near-wall region. With their precise origin unknown, these discrepancies could be related to the differences in $Pr$ and $\varGamma$.

Another interesting question that we can address based on our results concerns the evolution time scale $\mathcal {T}$ of the LSC. We estimate $\mathcal {T}$ as the time it takes a fluid parcel to complete a full cycle in the convection roll obtained from the conditional average. To do this we compute the streamline that passes through the location $z^*/H$ of the velocity maximum $\bar {v}_p(z^*/H)=\bar {v}_p^{max}$ at $d=0$ as shown in figure 3. The integrated travel time $\mathcal {T}$ along this averaged streamline as a function of $Ra$ is presented in figure 6(a). We find $\mathcal {T}/T_{ff} \gg 1$, i.e. the typical time scale of the LSC dynamics is much longer than the free-fall time $T_{ff} = \sqrt {H / (\beta g \varDelta )}$. Up to $Ra =10^7$ the time scale $\mathcal {T}$ grows approximately according to $\mathcal {T}/T_{ff} = (7.7 \pm 1.5) \times Ra^{0.139 \pm 0.014}$, but the trend flattens out at $Ra$ beyond that value. For the determination of all uncertainties in this manuscript we have used a $95\,\%$ confidence interval.

Figure 6. $(a)$ Time scale $\mathcal {T}$ versus $Ra$ using different methods. The datasets are: the time needed to circulate the flow along a streamline, which passes through $z^*/H$ at $d=0$ (red circles), see figure 3; the time scale calculated with the EAM method of Pandey et al. (Reference Pandey, Scheel and Schumacher2018) (blue squares). We also show the Pandey et al. (Reference Pandey, Scheel and Schumacher2018) data, which were calculated for the smaller $Pr=0.7$ (black diamonds). The dashed line shows $\mathcal {T}/T_{ff} = (7.7 \pm 1.5) \times Ra^{0.139 \pm 0.014}$. $(b)$ Average velocity $v_{{wind}}$ determined along the streamline chosen in $(a)$, normalised with $v_{RMS}$. $(c)$ Comparison between the length of the streamline and the circumference ${\rm \pi} (0.25 \hat {l}+0.5 H)$ of the ellipse (EAM method), both used to calculate the respective time scales in $(a)$.

To compare our results to other estimates in the literature, we also adopt the method used by Pandey et al. (Reference Pandey, Scheel and Schumacher2018) to estimate $\mathcal {T}$. These authors assumed the LSC to be an ellipse, used $v_{RMS}$ as the effective velocity scale and introduced a empirical prefactor of $3$ (which is equivalent to assuming a velocity scale $v_{RMS}/3$). The results for the ‘elliptical approximation method (EAM)’, using $v_{RMS}/3$ as the velocity scale, are compared to the corresponding results by Pandey et al. (Reference Pandey, Scheel and Schumacher2018) in figure 6(a). Results are consistent between the two methods in terms of the order of magnitude. However, the actual values, especially at lower $Ra$, differ significantly, and also the trends do not fully agree. The streamline approach allows us to determine the average convection velocity along the streamline $v_{{wind}}\equiv \mathcal {L}/\mathcal {T}$, where $\mathcal {L}$ is the length of the streamline. Figure 6(b) show that $v_{wind}$ is indeed proportional to $v_{RMS}$ with $v_{wind} \approx 0.45 \, v_{RMS}$ in the considered $Ra$ number regime. In figure 6(c), we present $\mathcal {L}$ along with the ellipsoidal estimate used in Pandey et al. (Reference Pandey, Scheel and Schumacher2018). From this, it appears that an ellipse does not very well represent the streamline geometry. Further, it becomes clear that it is the difference in the length-scale estimate that leads to the different scaling behaviours for $\mathcal {T}$ in figure 6(a).

It should additionally be noted that the present approach provides information on the typical turnover time scale of the superstructure in an averaged sense. This is different from Schneide et al. (Reference Schneide, Pandey, Padberg-Gehle and Schumacher2018) who studied turnover times for individual fluid particles. Particles may linger for long times in either the core of the structures or within the boundary layers, leading to a very wide distribution of time scales in the latter case.

4. Wall-shear stress and heat transport

The shear stress $\bar {\tau }_w$ at the plate surface is defined through

(4.1)\begin{equation} \bar{\tau}_w/\rho = - \nu \langle \partial_z \bar{v}_p \rangle_t. \end{equation}

Here, $\partial _z$ is the spatial derivative in the wall-normal direction. In figure 7(a) we show that the normalised shear stress $\bar {\tau }_w / \bar {\tau }_w^{max}$ as a function of the normalised distance $d/\hat {l}$ is nearly independent of $Ra$. Similar to findings in smaller cells (Wagner et al. Reference Wagner, Shishkina and Wagner2012), the curves are asymmetric with the maximum ($d/\hat {l} \approx -0.05$) shifted towards the plume impacting region. The value of $\bar {\tau }_w / \bar {\tau }_w^{max}$ drops to approximately 0.25 in both the plume impacting ($d/\hat {l} = -0.25$) and the plume emitting region ($d/\hat {l} = 0.25$).

Figure 7. $(a)$ Normalised shear stress $\bar {\tau }_w$ as a function of $d/\hat {l}$ and $(b)$ mean shear stress $\langle \bar {\tau }_w \rangle _J$ and maximum shear stress $\bar {\tau }_w^{max}$ versus $Ra$. The filled symbols show data of the present study ($\varGamma =32$ periodic domain), where the data for $Ra\geq 4\times 10^6$ can be fitted as $\bar {\tau }_w^{max} /(\rho V^2_{ff})=(0.12 \pm 0.04) \times Ra^{0.242 \pm 0.020}$ and $\langle \bar {\tau }_w \rangle _J /(\rho V^2_{ff})=(0.10 \pm 0.04) \times Ra^{0.236 \pm 0.016}$. The open symbols represent the data of Wagner et al. (Reference Wagner, Shishkina and Wagner2012) for $\varGamma =1$ with a cylindrical domain. The blue symbols show the maximum shear stress and the red symbols the mean shear stress over the interval $J=\{d/\hat {l} | d/\hat {l} \in [-0.2 : 0.15] \}$.

We use the good collapse of the $\bar {\tau }_w / \bar {\tau }_w^{max}$ profiles across the full range of $Ra$ considered to separate regions with significant shear from those with little to no lateral mean flow. We define the ‘wind’ region based on the approximate criterion $\bar {\tau }_w / \bar {\tau }_w^{max} \gtrapprox 0.5$, which leads to the interval $J=\{d/\hat {l}| d/\hat {l} \in [-0.2 : 0.15] \}$ that is indicated by the blue shading in figure 7(a). We use the average over this interval to evaluate the wind properties and indicate this by $\langle \rangle _J$. In figure 7(b) the data for mean $\langle \bar {\tau }_w\rangle _J$ and for maximum $\bar {\tau }_w^{max}$ wall-shear stress are compiled for the full range of $Ra$ considered. Both quantities are seen to increase significantly as $Ra$ increases. Around $Ra=1\text {--}4\times 10^6$ we can see a transition point at which the slope steepens. For lower $Ra$ the scaling of $\langle \bar {\tau }_w \rangle _{J}$ is much flatter. A fit to the data for $Ra\geq 4\times 10^6$ gives

(4.2)\begin{equation} \frac{\bar{\tau}_w /\rho}{V_{ff}^2} \sim Ra^{0.24}, \end{equation}

for both $\langle \bar {\tau }_w\rangle _J$ and $\bar {\tau }_w^{max}$. Overall, we find that the shear stress at the wall due to the turbulent thermal superstructures (in the periodic $\varGamma =32$ domain with $Pr=1$) compares well with the shear stress in a cylindrical $\varGamma =1$ domain by Wagner et al. (Reference Wagner, Shishkina and Wagner2012) with $Pr=0.786$. Most importantly, the scaling with $Ra$ is the same for both cases. The actual shear stress seems to be somewhat higher in the cylindrical aspect ratio $\varGamma =1$ domain than in the periodic domain in which the flow is unconfined. In part this difference may be related to the difference in $Pr$. Besides that, as we will show in the next section, the shear Reynolds number is slightly lower for the periodic domain than in the confined domain.

Next, we consider the local heat flux at the plate surface, given by

(4.3)\begin{equation} \overline{Nu}(d) \equiv - \frac{H}{\varDelta}\partial_z \bar{\theta}(d), \end{equation}

which is plotted in figure 8(a) for the full range of $Ra$. In all cases $\overline {Nu}/Nu$ is higher than one on the plume impacting side ($d<0$). This is consistent with the impacting cold plume increasing the temperature gradient in the BL locally. The fluid subsequently heats up while it is advected along the plate towards increasing $d$ by the LSC. As a consequence, the wall gradient is reduced and $\overline {Nu}$ decreases approximately linearly with increasing $d/\hat {l}$, which is consistent with observations by van Reeuwijk et al. (Reference van Reeuwijk, Jonker and Hanjalić2008) and Wagner et al. (Reference Wagner, Shishkina and Wagner2012). This leads to the ratio $\overline {Nu}/Nu$ dropping below 1 for $d>0$. For increasing $Ra$, the local heat flux becomes progressively more uniform across the full range of $d$. To quantify this, we plot the mean local heat fluxes in the plume impacting and emitting regions, respectively, in figure 8(b). The former is decreasing while the latter is increasing with increasing $Ra$, bringing the two sides closer. Again, and in both cases, a change of slope is visible in the range of $Ra=1\text {--}4\times 10^6$. In this context it is interesting to note that in a recent study on two-dimensional RB convection at $\varGamma =2$ (Zhu et al. Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018) it was found that at significantly higher $Ra \gtrapprox 10^{11}$ the heat transport in the plume emitting range dominated, reversing the current situation. If we boldly extrapolate the trend for $Ra\geq 4\times 10^6$ in our data, we can estimate that a similar reversal may occur at $Ra \approx {O}(10^{12}\text {--}10^{13})$, see figure 8(b).

Figure 8. $(a)$ Local heat flux $\overline {Nu}$ at the wall normalised by the global heat flux $Nu$ as function of the normalised spatial variable $d/ \hat {l}$. $(b)$ Values in the impacting ($-0.3 \leq d/ \hat {l} \leq -0.2$) and emitting ($0.2 \leq d/ \hat {l} \leq 0.3$) range as a function of $Ra$.

A possible mechanism that might explain this behaviour is increased turbulent (or convective) mixing, which can counteract the diffusive growth of the temperature BLs. To check this hypothesis, we plot the heat transport term $\widetilde {(\theta w)} \equiv \overline {w \theta } - \bar {w} \bar {\theta }$ in figure 9. The normalisation in the figure is with respect to the total heat flux $Nu$, the plotted quantity reflects the fraction of $Nu$ carried by $\widetilde {(\theta w)}$. It is obvious from these results that the convective transport contributes significantly, even within the BL height $\langle \bar {\lambda }_{\theta } \rangle _J$. Moreover, this relative contribution is independent of $Ra$ (except for the lowest value considered) in the plume impacting region (see figure 9a). However, figure 9(b) shows that already around $d = 0$ the convective transport in the BL increases with increasing $Ra$. This trend is much more pronounced in the plume emitting region $d>0.2$, see figure 9(c). Hence, convective transport in the BL plays an increasingly larger role for $d\geq 0$ with increasing $Ra$. Its effect is to smooth out the near-wall region, thereby increasing the temperature gradient at the wall. It is conceivable that the increased convective transport in the near-wall region (provided the trend persists) eventually leads to a reversal of the $\overline {Nu} (d)$ trend observed at moderate $Ra$ in figure 8(a).

Figure 9. Large-scale turbulent heat transport term $\widetilde {(\theta w)} / ( Nu \kappa \varDelta /H)$ evaluated in $(a)$ the plume impacting region, $(b)$ for small $|d|$ around zero and $(c)$ in the plume emitting region.

5. Thermal and viscous boundary layers

Next, we study how the BL thicknesses $\lambda _{\theta }$ and $\lambda _u$ vary along the LSC. In figure 10(a) we present $\bar {\lambda }_{\theta }$, normalised by $\lambda _{\theta }^*$. As expected from figure 8, $\bar {\lambda }_{\theta }$ is generally smaller in the plume impacting region and then increases along the LSC. However, unlike $\overline {Nu}$, $\bar {\lambda }_{\theta }$ is not determined by the gradient alone but also depends on the temperature level (see figure 4(b) for the definition of the level) such that differences arise. Specifically, $\bar {\lambda }_{\theta }/\lambda _{\theta }^*$ is rather insensitive for $Ra \geq 4\times 10^6$ in the plume impacting region ($d/\hat {l}<-0.1$). Furthermore, for $Ra \geq 10^7$, the growth of the thermal BL with $d/\hat {l}$ comes to an almost complete stop beyond $d= 0$, which is entirely consistent with the conclusions drawn in the discussion on $\widetilde {(\theta w)}$ above. Finally, we note that $\bar {\lambda }_{\theta }$ is generally larger than the estimate $\lambda _{\theta }^*$, which agrees with previous observations by Wagner et al. (Reference Wagner, Shishkina and Wagner2012).

Figure 10. $(a)$ Thermal BL thickness $\bar {\lambda }_{\theta }$ normalised by the estimated thermal BL thickness $\lambda _{\theta }^*$ and $(b)$ viscous BL thickness $\bar {\lambda }_u$ normalised by the mean viscous BL thickness in the interval $d/ \hat {l} \in J$ versus normalised distance $d/ \hat {l}$. The colour indicates the Rayleigh number, see legend.

There is no obvious choice for the normalisation of the viscous BL thickness and we therefore present $\bar {\lambda }_u$ normalised with its mean value $\langle \bar {\lambda }_u \rangle _{J}$ in figure 10. Overall, these curves for $\bar {\lambda }_u$ exhibit a similar trend as we observed previously for $\bar {\lambda }_{\theta }$. The values of $\bar {\lambda }_u$ are smaller in the plume impacting region ($d<0$) and the variation with $Ra$ is limited. Also for $\bar {\lambda } _u/ \langle \bar {\lambda }_u \rangle _{J}$ the growth with increasing $d$ is less pronounced the higher $Ra$ and the curves almost collapse for $d>0$ at $Ra\geq 10^7$.

Figure 11(a) shows $\langle \bar {\lambda }_{\theta } \rangle _{J}$ and $\langle \bar {\lambda }_u \rangle _{J}$ as a function of $Ra$. For the thermal BL thickness, the scaling appears to be rather constant over the full range and from fitting $4 \times 10^6 \leq Ra \leq 10^9$ we obtain

(5.1)\begin{equation} {\langle \bar{\lambda}_{\theta} \rangle_{J}/H = (4.7 \pm 0.9) \times Ra^{-0.298 \pm 0.012}}. \end{equation}

The reduction of the viscous BL thickness $\langle \bar {\lambda }_u \rangle _{J}$ with $Ra$ is significantly slower than for the thermal BL thickness $\langle \bar {\lambda }_{\theta } \rangle _{J}$. For low $Ra$, $\langle \bar {\lambda }_u \rangle _{J} < \langle \bar {\lambda }_{\theta } \rangle _{J}$. However, due to the different scaling of the two BL thicknesses, $\langle \bar {\lambda }_u \rangle _{J} > \langle \bar {\lambda }_{\theta } \rangle _{J}$ for $Ra \approx 4\times 10^6$. Comparing the periodic $\varGamma =32$ data with the confined $\varGamma =1$ case reported in Wagner et al. (Reference Wagner, Shishkina and Wagner2012) and Scheel & Schumacher (Reference Scheel and Schumacher2016), we note that the results for $\langle \bar {\lambda }_{\theta } \rangle$ agree closely between the two geometries. The scaling trends for $\langle \bar {\lambda }_u \rangle$ also appear to be alike in both cases. However, the viscous BL is significantly thinner in the smaller box, with a slight difference between the two datasets of $\varGamma =1$, which may be due to the difference in $Pr$. This situation is similar, and obviously also related to, the findings we reported for the comparison of the wall-shear stress in figure 7(b).

Figure 11. $(a)$ Mean BL thicknesses versus $Ra$ for the present data at $\varGamma = 32$ (filled symbols) and those of Wagner et al. (Reference Wagner, Shishkina and Wagner2012) (open symbols, W.S.W.) and Scheel & Schumacher (Reference Scheel and Schumacher2016) (x, S.S.), both with $\varGamma = 1$. The dashed line shows $\langle \bar {\lambda }_{\theta } \rangle _{J}/H = (4.7 \pm 0.9) \times Ra^{-0.298 \pm 0.012}$. $(b)$ Most probable BL ratio $\bar {\varLambda }^{MP}$ versus normalised distance $d/ \hat {l}$ for various $Ra$.

We further computed the most probable instantaneous BL ratio $\bar {\varLambda }^{MP} (d) \equiv \textrm {mode} (\varLambda (d(x,y),t |d = \textrm {const.}))$, where the ‘mode’-operator returns the most common value of the instantaneous BL ratio $\varLambda =\lambda _{\theta } / \lambda _u$, and present results in figure 11(b). Since the statistics of $\varLambda$ were found to be quite susceptible to outliers, we decided to report the most probable value $\bar {\varLambda }^{MP}$ as this provides a more robust measure than the mean. The Prandtl–Blasius BL theory for the flow over a flat plate suggests that $\varLambda = 1$ for $Pr= 1$. The figure shows that $\bar {\varLambda }^{MP}$ is almost constant as function of $d/\hat {l}$. However, unexpectedly, $\bar {\varLambda }^{MP}$ turns out to depend on $Ra$. For $Ra =10^5$, $\bar {\varLambda }^{MP}\approx 2$, which is larger than the theoretical prediction, but similar to the ratio of the means reported in figure 11(a). $\bar {\varLambda }^{MP}$ decreases with $Ra$ and approaches the predicted value of $1$ for $Ra=10^9$. We note that, although this $Ra$ dependence is not expected, it was also observed by e.g. Wagner et al. (Reference Wagner, Shishkina and Wagner2012).

When interpreting results for the BL thicknesses, it should be kept in mind that different definitions exist in the literature (du Puits et al. Reference du Puits, Resagk, Tilgner, Busse and Thess2007; Zhou & Xia Reference Zhou and Xia2010; Zhou et al. Reference Zhou, Stevens, Sugiyama, Grossmann, Lohse and Xia2010; Schmidt et al. Reference Schmidt, Calzavarini, Lohse, Toschi and Verzicco2012; du Puits, Resagk & Thess Reference du Puits, Resagk and Thess2013; Zhou & Xia Reference Zhou and Xia2013; Scheel & Schumacher Reference Scheel and Schumacher2014; Shishkina et al. Reference Shishkina, Horn, Wagner and Ching2015, Reference Shishkina, Horn, Emran and Ching2017b; Ching et al. Reference Ching, Leung, Zwirner and Shishkina2019). We note that values may depend on the boundary layer definition that is employed. To get at least a sense for which of the observations transfer to other possible BL definitions, we compare the results for $\lambda$ (the slope method) to those obtained by the location of the temperature and velocity levels $(\bar {\ell })$ (level method, see figure 4b) in figure 12. We note that the scalings versus $Ra$ are very similar, albeit not exactly the same, for both definitions of the BL thickness. However, the offset between $\bar {\lambda }$ and $\bar {\ell }$ is not the same for velocity and temperature. As a consequence, there is no cross-over between $\bar {\ell _{\theta }}$ and $\overline {\ell _u}$ within the range of $Ra$ considered.

Figure 12. Comparison of mean BL thicknesses versus $Ra$ using the slope method and the location of the respective temperature and velocity levels (level method). The dashed lines show $\langle \bar {\lambda }_{\theta } \rangle _{J}/H = (4.7 \pm 0.9) \times Ra^{-0.298 \pm 0.012}$ and $\langle \bar {\ell }_{\theta } \rangle _{J}/H = (8.6 \pm 3.1) \times Ra^{-0.283 \pm 0.023}$.

In figure 13 we compare the wind Reynolds number, which we determined as follows,

(5.2)\begin{equation} Re_{{wind}} = \langle {\bar{u}_l} \rangle_{J} \, H / \nu , \end{equation}

with the results of Wagner et al. (Reference Wagner, Shishkina and Wagner2012) and Scheel & Schumacher (Reference Scheel and Schumacher2016). The figure shows that our $Re_{{wind}}$ obtained in a periodic $\varGamma =32$ domain with $Pr=1$ agree surprisingly well with the results from Wagner et al. (Reference Wagner, Shishkina and Wagner2012) obtained in a cylindrical $\varGamma =1$ sample with $Pr=0.786$ and with the data of Scheel & Schumacher (Reference Scheel and Schumacher2016) for $Pr=0.7$. The $Re$ values obtained by Wagner et al. (Reference Wagner, Shishkina and Wagner2012) and Scheel & Schumacher (Reference Scheel and Schumacher2016) are slightly higher than our values. We note that the lower $Pr$ results in slightly higher $Re_{wind}$. This means that the main finding in this context is that $Re_{wind}$ in the turbulent superstructures is almost the same, perhaps slightly lower, than in a confined $\varGamma =1$ sample. We note that the predictions for the wind Reynolds number obtained from the unifying theory for thermal convection (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001) are in good agreement with the data. The unifying theory, using the updated constants found by Stevens et al. (Reference Stevens, van der Poel, Grossmann and Lohse2013), namely predicts that for $Pr=1$ the wind Reynolds number scales as $Re_{GL}=0.40\times Ra^{0.44}$, while the data for $Ra \geq 4\times 10^6$ are well approximated by $Re_{wind} = (0.22 \pm 0.05) \times Ra^{0.468 \pm 0.012}$.

Figure 13. $(a)$ Wind Reynolds number $Re_{wind}$ versus $Ra$ obtained in a periodic $\varGamma =32$ domain compared to the corresponding values obtained by Wagner et al. (Reference Wagner, Shishkina and Wagner2012) and Scheel & Schumacher (Reference Scheel and Schumacher2016), both for a cylindrical $\varGamma =1$ domain. We also show the predictions from the unifying theory (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001) using the updated prefactors (Stevens et al. Reference Stevens, van der Poel, Grossmann and Lohse2013). The blue dashed line shows $Re_{wind} = (0.22 \pm 0.05) \times Ra^{0.468 \pm 0.012}$. $(b)$ Plot of $Re_s$ versus $Ra$ with estimations for $Ra_{crit}$. In (a,b), we have fitted our own data points only from $Ra=4\times 10^6$ onwards to achieve consistent comparisons with the data by Wagner et al. (Reference Wagner, Shishkina and Wagner2012), where only data from $Ra=3\times 10^6$ on are available. The blue dashed line shows $Re_s = (0.09 \pm 0.04) \times Ra^{0.243 \pm 0.025}$.

To estimate when the BLs become turbulent we calculate the shear Reynolds number

(5.3)\begin{equation} Re_s= \frac{ \left [ {\bar{u}_l} \times \lambda_u^{MP} \right ]^{max} }{2 \nu}. \end{equation}

We expect the BL to become turbulent and the ultimate regime to set in (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2011) at a critical shear Reynolds number of $Re_s^{crit} \approx 420$ (Landau & Lifshitz Reference Landau and Lifshitz1987). A fit to our data gives

(5.4)\begin{equation} {Re_s=(0.09 \pm 0.04)\times Ra^{0.243 \pm 0.025}}, \end{equation}

from which we can extrapolate that $Re_s^{crit}=420$ is reached at $Ra_{crit} \approx 1.3\times 10^{15}$. Of course, this estimate comes with a significant error bar as our data for $\varGamma =32$ are still far away from the expected critical $Ra$ number. Nevertheless, it agrees well with the result from Wagner et al. (Reference Wagner, Shishkina and Wagner2012), who find $Ra_{crit}\approx 1.2\times 10^{14}$, and Scheel & Schumacher (Reference Scheel and Schumacher2016), who determined $Ra_{crit}=(3\pm 2) \times 10^{13}$, both for a cylindrical $\varGamma =1$ cell, and the results from Sun et al. (Reference Sun, Cheung and Xia2008) who find from experiments that $Ra_{crit}\approx 2\times 10^{13}$. We emphasise that all these estimates are consistent with the observation of the onset of the ultimate regime at $Ra_{*}\approx 2\times 10^{13}$ in the Göttingen experiments (He et al. Reference He, Funfschilling, Nobach, Bodenschatz and Ahlers2012, Reference He, van Gils, Bodenschatz and Ahlers2015). As is explained by Ahlers, Bodenschatz & He (Reference Ahlers, Bodenschatz and He2017) also measurements of the shear Reynolds number in low $Pr$ number simulations by Schumacher et al. (Reference Schumacher, Bandaru, Pandey and Scheel2016) support the observation of the ultimate regime in the Göttingen experiments.