## 1. Introduction

Rayleigh–Bénard (RB) convection (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; Lohse & Xia Reference Lohse and Xia2010; Chilla & Schumacher Reference Chilla and Schumacher2012; Xia Reference Xia2013) 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 Reference Rahmstorf2000), the large-scale flow patterns that are formed in the outer core of the Earth (Glatzmaier *et al.* Reference Glatzmaier, Coe, Hongre and Roberts1999), where reversals of the large-scale convection roll are of prime importance, convection in gaseous giant planets (Busse Reference Busse1994) and in the outer layer of the Sun (Miesch Reference Miesch2000). 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 Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001, Reference Grossmann and Lohse2011; Stevens *et al.* Reference Stevens, van der Poel, Grossmann and Lohse2013) 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.* Reference Ahlers, Grossmann and Lohse2009; Zhou & Xia Reference Zhou and Xia2010; Zhou *et al.* Reference Zhou, Stevens, Sugiyama, Grossmann, Lohse and Xia2010; Shi, Emran & Schumacher Reference Shi, Emran and Schumacher2012; Stevens *et al.* Reference Stevens, Zhou, Grossmann, Verzicco, Xia and Lohse2012; Shishkina *et al.* Reference Shishkina, Horn, Wagner and Ching2015; Schumacher *et al.* Reference Schumacher, Bandaru, Pandey and Scheel2016; Shishkina *et al.* Reference Shishkina, Emran, Grossmann and Lohse2017*a*). Previous studies by, for example, Wagner, Shishkina & Wagner (Reference Wagner, Shishkina and Wagner2012) and Scheel & Schumacher (Reference Scheel and Schumacher2016), 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.* (Reference Wagner, Shishkina and Wagner2012) 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 (Reference Scheel and Schumacher2016) 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 (Reference Hartlep, Tilgner and Busse2003) and von Hardenberg *et al.* (Reference von Hardenberg, Parodi, Passoni, Provenzale and Spiegel2008) 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.* (Reference Stevens, Blass, Zhu, Verzicco and Lohse2018) 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 (Reference Pandey, Scheel and Schumacher2018) 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.* (Reference Parodi, von Hardenberg, Passoni, Provenzale and Spiegel2004) described as clustering of thermal plumes originating in the BL and assembling the LSC. Krug, Lohse & Stevens (Reference Krug, Lohse and Stevens2020) 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.* (Reference Stevens, Blass, Zhu, Verzicco and Lohse2018) 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 (Reference Berghout, Baars and Krug2020) 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.* Reference Wagner, Shishkina and Wagner2012) and summarise our findings in § 6.

## 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:

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).

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 1*a*,*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.

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

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

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$.

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.

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.

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$.

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.

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

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$).

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

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

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*).

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 9*a*). 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*).

## 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).

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

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*).

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 Ching2017*b*; 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 4*b*) 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.

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

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}$.

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

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

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.

## 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.* (Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018). 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 Reference Landau and Lifshitz1987). 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 (Reference Shishkina, Wagner and Horn2014), 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.* (Reference Wagner, Shishkina and Wagner2012) 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.* (Reference Wagner, Shishkina and Wagner2012). 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.* Reference He, Funfschilling, Nobach, Bodenschatz and Ahlers2012, Reference He, van Gils, Bodenschatz and Ahlers2015), as well as with previous measurements of the shear Reynolds number (Wagner *et al.* Reference Wagner, Shishkina and Wagner2012; Scheel & Schumacher Reference Scheel and Schumacher2016). 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.

## Acknowledgements

We greatly appreciate valuable discussions with O. Shishkina. This work is supported by NWO, the University of Twente Max-Planck Center for Complex Fluid Dynamics, the German Science Foundation (DFG) via program SSP 1881, and the ERC (the European Research Council) Starting Grant No. 804283 UltimateRB. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre (www.lrz.de).

## Declaration of interests

The authors report no conflict of interest.

## Appendix. Contribution of the horizontal mean velocity

In figure 14 we show that at mid-height the contribution of the horizontal mean velocity $v_{h_{mean}}=\langle (v_x+v_y)/2 \rangle _A$ on the horizontal root mean square velocity $v_{h_{RMS}}=\sqrt {\langle v_x^2 + v_y^2 \rangle _A}$ is approximately one percentage point. This is in agreement with Hartlep *et al.* (Reference Hartlep, Tilgner and Busse2003), who find that the energy contained in the mean flow is less than $0.8\,\%$ of the total kinetic energy.