Flow organisation in laterally unconfined Rayleigh–Bénard turbulence

Abstract

We investigate the large-scale circulation (LSC) of turbulent Rayleigh–Bénard convection in a large box of aspect ratio $\varGamma =32$ for Rayleigh numbers up to $Ra=10^9$ and at a fixed Prandtl number $Pr=1$. A conditional averaging technique allows us to extract statistics of the LSC even though the number and the orientation of the structures vary throughout the domain. We find that various properties of the LSC obtained here, such as the wall-shear stress distribution, the boundary layer thicknesses and the wind Reynolds number, do not differ significantly from results in confined domains ($\varGamma \approx 1$). This is remarkable given that the size of the structures (as measured by the width of a single convection roll) more than doubles at the highest $Ra$ as the confinement is removed. An extrapolation towards the critical shear Reynolds number of $Re_s^{{crit}} \approx 420$, at which the boundary layer (BL) typically becomes turbulent, predicts that the transition to the ultimate regime is expected at $Ra_{{crit}} \approx {O}(10^{15})$ in unconfined geometries. This result is in line with the Göttingen experimental observations (He et al., Phys. Rev. Lett., vol. 108, 2012, 024502; New J. Phys., vol. 17, 2015, 063028). Furthermore, we confirm that the local heat transport close to the wall is highest in the plume impacting region, where the thermal BL is thinnest, and lowest in the plume emitting region, where the thermal BL is thickest. This trend, however, weakens with increasing $Ra$.

1. Introduction

Rayleigh–Bénard (RB) convection (Ahlers, Grossmann & Lohse 2009; Lohse & Xia 2010; Chilla & Schumacher 2012; Xia 2013) is the flow in a box heated from below and cooled from above. Such buoyancy driven flow is the paradigmatic example for natural convection which often occurs in nature, e.g. in the atmosphere. For that case, a large-scale horizontal flow organisation is observed in satellite pictures of weather patterns. Other examples include the thermohaline circulation in the oceans (Rahmstorf 2000), the large-scale flow patterns that are formed in the outer core of the Earth (Glatzmaier et al. 1999), where reversals of the large-scale convection roll are of prime importance, convection in gaseous giant planets (Busse 1994) and in the outer layer of the Sun (Miesch 2000). Thus, the problem is of interest in a wide range of scientific disciplines, including geophysics, oceanography, climatology and astrophysics.

For a given aspect ratio and given geometry, the dynamics in RB convection is determined by the Rayleigh number $Ra=\beta g\varDelta H^3 /(\kappa \nu )$ and the Prandtl number $Pr=\nu /\kappa$. Here, $\beta$ is the thermal expansion coefficient, $g$ the gravitational acceleration, $\varDelta$ the temperature difference between the horizontal plates, which are separated by a distance $H$, and $\nu$ and $\kappa$ are the kinematic viscosity and thermal diffusivity, respectively. The dimensionless heat transfer, i.e. the Nusselt number $Nu$, along with the Reynolds number $Re$ are the most important response parameters of the system.

For sufficiently high $Ra$, the flow becomes turbulent, which means that there are vigorous temperature and velocity fluctuations. Nevertheless, a large-scale circulation (LSC) develops in the domain such that, in addition to the thermal boundary layer (BL), a thin kinetic BL is formed to accommodate the no-slip boundary condition near both the bottom and top plates. Properties of the LSC and the nature of the BLs are highly relevant to the theoretical description of the problem. In particular, the unifying theory of thermal convection (Grossmann & Lohse 2000, 2001, 2011; Stevens et al. 2013) states that the transition from the classical to the ultimate regime takes place when the kinetic BLs become turbulent. This transition is shear based and driven by the large-scale wind, underlying the importance of the LSC to the overall flow behaviour.

So far, the LSC and BL properties have mainly been studied in cells featuring a small aspect ratio $\Gamma$, typically $\varGamma =1/2$ or $\varGamma =1$. Various studies have shown that the BLs indeed follow the laminar Prandtl–Blasius (PB) type predictions in the classical regime (Ahlers et al. 2009; Zhou & Xia 2010; Zhou et al. 2010; Shi, Emran & Schumacher 2012; Stevens et al. 2012; Shishkina et al. 2015; Schumacher et al. 2016; Shishkina et al. 2017a). Previous studies by, for example, Wagner, Shishkina & Wagner (2012) and Scheel & Schumacher (2016), have used results from direct numerical simulations (DNS) in aspect ratio $\varGamma =1$ cells to study the properties of the BLs in detail. Wagner et al. (2012) showed that an extrapolation of their data gives that for $Pr=0.786$ the critical shear Reynolds number of $420$ is reached at $Ra\approx 1.2\times 10^{14}$, while Scheel & Schumacher (2016) predict a value of $Ra\approx 3\times 10^{13}$.

Despite the wealth of studies in low aspect ratio domains, many natural instances of thermal convection take place in very large aspect ratio systems, as mentioned above. Previous research has demonstrated that several flow properties are significantly different in such unconfined geometries. Hartlep, Tilgner & Busse (2003) and von Hardenberg et al. (2008) performed DNS at $Ra={O}(10^7)$ and $\varGamma =20$. They observed large-scale structures by investigating the advective heat transport and found the most energetic wavelength of the LSC at $4H\text {--}7H$. Recently, DNS by Stevens et al. (2018) for $\varGamma =128$ and $Ra={O}(10^7\text {--}10^9)$ also reported ‘superstructures’ with wavelengths of $6\text {--}7$ times the distance between the plates. Similar findings were made by Pandey, Scheel & Schumacher (2018) over a wide range of Prandtl numbers $0.005\leq Pr \leq 70$ and $Ra$ up to $10^7$. It was shown that the signatures of the LSC can be observed close to the wall, which Parodi et al. (2004) described as clustering of thermal plumes originating in the BL and assembling the LSC. Krug, Lohse & Stevens (2020) showed that the presence of the LSC leads to a pronounced peak in the coherence spectrum of temperature and wall-normal velocity. Based on DNS at $\varGamma =32$ and $Ra={O}(10^5\text {--}10^9)$, they determined that the wavelength of this peak shifts from $\hat {l}/H \approx 4$ to $\hat {l}/H \approx 7$ as $Ra$ is increased.

Stevens et al. (2018) have shown that, in periodic domains, the heat transport is maximum for $\varGamma =1$ and reduces with increasing aspect ratio up to $\varGamma \approx 4$ when the large-scale value is obtained. They also found that fluctuation-based Reynolds numbers depend on the aspect ratio of the cell. However, other than the structure size, it is mostly unclear how the large-scale flow organisation and BL properties are affected by different geometries. Not only is the size of the LSC more than 2 times larger without confinement (note that $\skew2\hat{l}$ measures the size of two counter-rotating rolls combined), but also other effects, such as corner vortices, are absent in periodic domains. Therefore one would expect differences in wind properties and BL dynamics. It is the goal of this paper to investigate these differences. Doing so comes with significant practical difficulty due to the random orientation of a multitude of structures that are present in a large box. To overcome this, we adopt the conditional averaging technique that was devised in Berghout, Baars & Krug (2020) to reliably extract LSC features even under these circumstances. Details on this procedure are provided in section § 3 after a short description of the dataset in § 2. Finally, in §§ 4 and 5 we present results on how superstructures affect the flow properties in comparison to the flow formed in a cylindrical $\varGamma =1$ domain (Wagner et al. 2012) and summarise our findings in § 6.

2. Numerical method

The data used in this manuscript have previously been presented by Stevens et al. (2018) and Krug et al. (2020). 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 1996; van der Poel et al. 2015). 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 (1996), Verzicco & Camussi (1997, 2003), Stevens, Verzicco & Lohse (2010) and Kooij et al. (2018).

Table 1. Data from Stevens et al. (2018) and Krug et al. (2020) 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. 2010; Stevens et al. 2010). 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. (2010) 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 2008; Wagner et al. 2012) to intermediate (van Reeuwijk, Jonker & Hanjalić 2008) aspect ratio domains. To overcome this issue, we use a conditional averaging technique developed in Berghout et al. (2020), 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. (2020) 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. 2020). 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. 2020) (b) Coherence spectra of temperature and wall-normal velocity at mid-height, figure adopted from Krug et al. (2020). 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. (2020)), 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. (2020) 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. (2020) 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. (2012) 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. (2008) 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. (2008), 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. (2018) (blue squares). We also show the Pandey et al. (2018) 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. (2018) 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. (2018) 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. (2018). 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. (2018) 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. 2012), 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. (2012) 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. (2012) 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. (2008) and Wagner et al. (2012). 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. 2018) 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. (2012).

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. (2012) and Scheel & Schumacher (2016), 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. (2012) (open symbols, W.S.W.) and Scheel & Schumacher (2016) (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. (2012).

When interpreting results for the BL thicknesses, it should be kept in mind that different definitions exist in the literature (du Puits et al. 2007; Zhou & Xia 2010; Zhou et al. 2010; Schmidt et al. 2012; du Puits, Resagk & Thess 2013; Zhou & Xia 2013; Scheel & Schumacher 2014; Shishkina et al. 2015, 2017b; Ching et al. 2019). 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. (2012) and Scheel & Schumacher (2016). 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. (2012) obtained in a cylindrical $\varGamma =1$ sample with $Pr=0.786$ and with the data of Scheel & Schumacher (2016) for $Pr=0.7$. The $Re$ values obtained by Wagner et al. (2012) and Scheel & Schumacher (2016) 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 2000, 2001) are in good agreement with the data. The unifying theory, using the updated constants found by Stevens et al. (2013), 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. (2012) and Scheel & Schumacher (2016), both for a cylindrical $\varGamma =1$ domain. We also show the predictions from the unifying theory (Grossmann & Lohse 2000, 2001) using the updated prefactors (Stevens et al. 2013). 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. (2012), 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 2000, 2011) at a critical shear Reynolds number of $Re_s^{crit} \approx 420$ (Landau & Lifshitz 1987). 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. (2012), who find $Ra_{crit}\approx 1.2\times 10^{14}$, and Scheel & Schumacher (2016), who determined $Ra_{crit}=(3\pm 2) \times 10^{13}$, both for a cylindrical $\varGamma =1$ cell, and the results from Sun et al. (2008) 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. 2012, 2015). As is explained by Ahlers, Bodenschatz & He (2017) also measurements of the shear Reynolds number in low $Pr$ number simulations by Schumacher et al. (2016) support the observation of the ultimate regime in the Göttingen experiments.

6. Conclusions

We have used a conditional averaging technique to investigate the properties of the LSC and the boundary layers in $\varGamma =32$ RB convection for unit Prandtl number and Rayleigh numbers up to $Ra=10^9$. The resulting quasi-two-dimensional representation of the LSC allowed us to analyse the wind properties as well as wall shear and local heat transfer. We found the distribution of the wall-shear stress $\bar {\tau }_w$ to be asymmetric. The maximum of $\bar {\tau }_w$ is located closer to the plume impacting side and its value increases as $\bar {\tau }_w^{max} / (\rho V_{ff} ^2) \sim Ra^{0.24}$ with increasing $Ra$. The local heat transfer at the wall, represented by the conditioned Nusselt number $\overline {Nu}$, has its highest values in the plume impacting zone at all $Ra$ considered here. Going from the plume impacting towards the plume emitting region, $\overline {Nu}$ is seen to decrease consistently as is expected from the fluid near the hot wall heating up. However, as $Ra$ is increased, the differences in $\overline {Nu}$ even out more and more. For the plume emitting side in particular, we were able to connect this trend to increased advective transport in the wall-normal direction at higher $Ra$. When extrapolating the trends for $\overline {Nu}$ to $Ra$ higher than those available here, our results appear consistent with Zhu et al. (2018). These authors observed a reversal of the $\overline {Nu}$-distribution in two-dimensional RB turbulence above $Ra \gtrapprox 10^{11}$ with higher values of the heat transport in the emitting region.

Further, we examined the thermal and the viscous BLs. At low $Ra$, both increase along $d$ in an approximately linear fashion, whereas flat plate boundary layer theory would suggest a growth proportional $\sqrt {d}$ (Landau & Lifshitz 1987). As $Ra$ increases, and especially for $d>0$, the growth becomes successively weaker and stops entirely beyond $Ra \gtrapprox 10^8$. Again, this is likely a consequence of the increased convective mixing in this region. For increasing $Ra$, both $\bar {\lambda }_{\theta }$ and $\bar {\lambda }_u$ become thinner, with $\bar {\lambda }_{\theta }$ showing an effective scaling of $\langle \bar {\lambda }_{\theta } \rangle _J /H \sim Ra^{-0.3}$. At $Ra\gtrapprox 4\times 10^6$ we observed a cross-over point where the thermal BL becomes smaller than the viscous BL. It should be noted that the cross-over appears specific to the definition of $\bar {\lambda }$ since a similar behaviour was not observed when an alternative definition ($\bar {\ell }$, based on the location of the level) was employed. Nevertheless, the scaling behaviour of $\bar {\lambda }$ and $\bar {\ell }$ was seen to be very similar. When calculating instantaneous BL ratios, a convergence to $\bar {\varLambda }^{MP} \rightarrow 1$ for high enough $Ra$ can be observed as predicted by the PB theory for laminar BLs. As pointed out in Shishkina, Wagner & Horn (2014), the PB limit only strictly applies to wall parallel flow and the ratio is expected to be higher if the flow approaches the plate at an angle. This incidence angle is higher at smaller $\varGamma$ which can explain why at comparable $Ra$ the BL ratios reported in Wagner et al. (2012) are slightly higher than what is found here.

We expected to find significant differences in the LSC statistics obtained in a confined $\varGamma =1$ system and a large $\varGamma =32$ system. However, surprisingly, we find that the thermal BL thickness $\langle \bar {\lambda }_{\theta } \rangle _J$ obtained for both cases agrees very well. It turns out that the viscous BL thickness $\langle \bar {\lambda }_u \rangle _J$ is significantly larger (${\approx }55\text {--}65\,\%$) for the periodic $\varGamma =32$ case than in a $\varGamma =1$ cylinder. However, the wall shear and its scaling with $Ra$ are similar in both cases. Here we find that in a periodic $\varGamma =32$ domain, the shear Reynolds number scales as $Re_s \sim Ra^{0.243}$. This is a bit lower than the corresponding result for $\varGamma =1$, although one needs to keep in mind the slight difference in $Pr$ ($Pr = 0.786$ at $\varGamma = 1$ vs. $Pr = 1$ for $\varGamma =32$) is responsible for part of the observed difference. An extrapolation towards the critical shear Reynolds number of $Re_s^{{crit}} \approx 420$ when the laminar-type BL becomes turbulent predicts that the transition to the ultimate regime is expected at $Ra_{{crit}} \approx {O}(10^{15})$. This is slightly higher than the corresponding result for a $\varGamma =1$ cylinder, i.e. $Ra_{{crit}} \approx {O}(10^{14})$, by Wagner et al. (2012). However, it should be noted that considering inherent uncertainties and differences in $Pr$, the results for $\varGamma =32$ agree with the observed transition to the ultimate regime in the Göttingen experiments (He et al. 2012, 2015), as well as with previous measurements of the shear Reynolds number (Wagner et al. 2012; Scheel & Schumacher 2016). So surprisingly, we find that in essentially unconfined very large aspect ratio systems, in which the resulting structure size is significantly larger, the differences in terms of $Re_{{wind}}$ or $Re_s$ with respect to the $\varGamma =1$ cylindrical case are marginal.