1. Introduction
Natural convection is an important flow configuration with highly relevant applications in thermal convection in the Earth's atmosphere (Hartmann, Moy & Fu Reference Hartmann, Moy and Fu2001) and for stars (Cattaneo, Emonet & Weiss Reference Cattaneo, Emonet and Weiss2003) among many others. Fundamental aspects of natural convection are traditionally studied in the Rayleigh–Bénard set-up, in which the flow evolves between parallel horizontal plates that are heated from below and cooled from above (Siggia Reference Siggia1994; Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; Lohse & Xia Reference Lohse and Xia2010; Chilla & Schumacher Reference Chilla and Schumacher2012). In this case, the strength of the thermal driving is expressed by the dimensionless Rayleigh number 
$Ra$, while, the dimensionless heat flux, which is the most important system response, is quantified by the Nusselt number 
$Nu$.
A characteristic and very persistent feature of Rayleigh–Bénard convection (RBC) is the formation of a large-scale circulation (LSC). The LSC typically consists of one or several roll structures that fill the entire gap between the plates. Such structures continue to exist even in the highly turbulent state, at large 
$Ra$. It has further been amply demonstrated (e.g. Fitzjarrald Reference Fitzjarrald1976; Hartlep, Tilgner & Busse Reference Hartlep, Tilgner and Busse2003; Parodi et al. Reference Parodi, von Hardenberg, Passoni, Provenzale and Spiegel2004; von Hardenberg et al. Reference von Hardenberg, Parodi, Passoni, Provenzale and Spiegel2008; Pandey, Scheel & Schumacher Reference Pandey, Scheel and Schumacher2018; Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018; Krug, Lohse & Stevens Reference Krug, Lohse and Stevens2020) that the lateral extent of the LSC can reach multiple times the cell height 
$H$ in large aspect ratio domains. In laterally unconfined geometries, the LSC has therefore also been referred to as a ‘superstructure’ in the more recent literature.
The existence of such large-scale organization in an otherwise chaotic flow is arguably an exciting aspect of RBC, and is reminiscent of large-scale structures in other canonical flows (Hutchins & Marusic Reference Hutchins and Marusic2007a; Lee & Moser Reference Lee and Moser2018). Their slow spatio-temporal evolution renders these structures prime candidates for flow modelling approaches. Such efforts are indeed already under way, e.g. by Fonda et al. (Reference Fonda, Pandey, Schumacher and Sreenivasan2019), who used deep learning techniques to analyse superstructures up to a moderate 
$Ra$ of 
$10^7$. However, these authors also acknowledged that the relative contribution of the LSC to the total heat transport decreases with increasing 
$Ra$. Similarly, Krug et al. (Reference Krug, Lohse and Stevens2020) found that at 
$Ra = 10^8$, superstructures account for a maximum of 30 % of the total heat transport. It is therefore clear that the ‘background’ turbulence (i.e. small to intermediate scales) plays an essential role that needs to be accounted for in low-order models. This immediately raises the question on how structures of different sizes ‘interact’ in RBC or more precisely: How does a certain large-scale state affect local properties of the small-scale field? This is exactly what is going to be addressed in this study. One should note, however, that the interactions investigated here are not the same as the interscale transport of energy studied in e.g. Gayen, Hughes & Griffiths (Reference Gayen, Hughes and Griffiths2013), Togni, Cimarelli & De Angelis (Reference Togni, Cimarelli and De Angelis2015) or Green et al. (Reference Green, Vlaykov, Mellado and Wilczek2020).
Related findings on spatial variations in the small-scale statistics have recently been reported by He & Xia (Reference He and Xia2019). These authors found that in the centre regions of the LSC (i.e. in the wind-shear regions) the decay of temperature fluctuations away from the wall follows a power law dependence. Whereas in the plume emitting parts of the boundary layer (BL), temperature fluctuations drop off logarithmically with wall distance 
$z$. Relevant in the broader context are the studies by Shishkina et al. (Reference Shishkina, Horn, Wagner and Ching2015), Wang, He & Tong (Reference Wang, He and Tong2016) and Wang et al. (Reference Wang, Xu, He, Yik, Wang, Schumacher and Tong2018), who derived a BL equation for the mean and variance of the temperature profile, where they included the effect of the fluctuations into the otherwise laminar two-dimensional (2-D) BL equation. Further, Wang, He & Tong (Reference Wang, He and Tong2019) disentangled temperature fluctuations throughout the RBC cell into homogeneous background (Gaussian) and non-homogeneous plume contributions (exponential).
The approach we choose here is to identify and to quantify modulation effects between large and small scales in terms of amplitude and wavenumber. Naturally, this first and foremost necessitates a clear definition of what those ‘large’ (i.e. related to the LSC) and ‘small’ (the rest) scales are. Such a distinction can readily be made for a superstructure configuration (aspect ratio 
$\varGamma = 32$) at 
$Ra = 10^8$, as figure 1 demonstrates. The premultiplied temperature spectrogram 
$k\varPhi _{\theta \theta }$ shown there clearly separates into a large-scale peak at radial wavenumber 
$k \approx 1$ (normalized by the cell height 
$H$ as all length scales here) and one at significantly smaller scales (
$k \sim 10^2$). There is only little energy in the spectral gap between the two peaks, such that a scale decomposition can be achieved via spectral filtering. Here, we chose a cutoff wavenumber 
$k_{cut} = 2.5$ (indicated by the red line in the figure). Note, however, that the results are rather insensitive to the exact choice of 
$k_{cut}$ as long as the value lies within the spectral gap, since in this case small- and large-scale signals are dominated by energetic scales at wavenumber ranges far from 
$k_{cut}$. The normalized wavelength 
$\hat {l} = 2 {\rm \pi}/\hat {k}$ (where 
$\hat {k}$ is a characteristic wavenumber of the LSC based on a coherence metric, see Krug et al. Reference Krug, Lohse and Stevens2020) of the large-scale structures increases with increasing 
$Ra$ (Pandey et al. Reference Pandey, Scheel and Schumacher2018; Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018; Krug et al. Reference Krug, Lohse and Stevens2020), while the structure size corresponding to the small-scale peak scales with the thermal BL height 
$\delta _\theta$ (Krug et al. Reference Krug, Lohse and Stevens2020). As demonstrated in Blass et al. (Reference Blass, Verzicco, Lohse, Stevens and Krug2021), this leads to an increasing scale separation between the peaks at higher 
$Ra$, rendering their interaction even more clear and relevant. It should be noted in this context that others have used spatial filtering (Green et al. Reference Green, Vlaykov, Mellado and Wilczek2020) or time averaging (e.g. Pandey et al. Reference Pandey, Scheel and Schumacher2018; Fonda et al. Reference Fonda, Pandey, Schumacher and Sreenivasan2019) to extract the large-scale features in similar datasets.
Figure 1. Premultiplied temperature spectrogram 
$k\varPhi _{\theta \theta }$ at 
$Ra = 10^8$ and 
$\varGamma = 32$, where 
$k$ is the wavenumber and 
$\lambda =(2{\rm \pi} )/k$ the wavelength. The white cross located at 
$z = 0.85 \delta _\theta$ and 
$k = 34$ (
$\lambda /\delta _\theta \approx 11.4$) indicates the location of the small-scale or ‘inner’ peak. The black dashed-dotted line marks the superstructure scale 
$\hat {l}$ determined in Krug et al. (Reference Krug, Lohse and Stevens2020) based on the coherence metric and the red dashed line is located at the cutoff wavenumber 
$k_{{cut}} = 2.5$ used here.
For the reasons outlined above, large aspect ratio cells are convenient for the present purpose and are therefore considered in the following. However, we note that there is little reason to believe that the effects studied are restricted to such domains. In fact, we expect scale interactions similar to those uncovered here to also play a role in smaller cells. Yet, spectral overlap of the peaks and other complications, e.g. due to corner vortices, render their analysis much more cumbersome in these cases.
Finally, we point out that a considerable number of studies on scale interactions in other flows exist, tracing back to Brown & Thomas (Reference Brown and Thomas1977) and Bandyopadhyay & Hussain (Reference Bandyopadhyay and Hussain1984). Especially in the field of turbulent BLs, where a scale separation reminiscent of the situation displayed in figure 1 occurs (Hutchins & Marusic Reference Hutchins and Marusic2007b), this topic has attracted considerable attention over the last decade (Chung & McKeon Reference Chung and McKeon2009; Mathis, Hutchins & Marusic Reference Mathis, Hutchins and Marusic2009; Marusic, Mathis & Hutchins Reference Marusic, Mathis and Hutchins2010; Schlatter & Örlü Reference Schlatter and Örlü2010; Bernardini & Pirozzoli Reference Bernardini and Pirozzoli2011; Ganapathisubramani et al. Reference Ganapathisubramani, Hutchins, Monty, Chung and Marusic2012; Baars, Hutchins & Marusic Reference Baars, Hutchins and Marusic2017; Dogan et al. Reference Dogan, Örlü, Gatti, Vinuesa and Schlatter2019). We will make use of the mature framework and tools developed in this context to quantify amplitude modulation and wavenumber modulation effects in RBC.
The paper is structured as follows; in § 2 we present details of the numerical dataset along with relevant parameters. We then quantify modulation effects using two independent analysis tools. First, using a wavelet-based method (§§ 3.1 and 3.2) and then in § 3.3 on the basis of conditioned statistics. We conclude by discussing some physical implications of the results in § 4.
2. Governing equations and parameters
The simulations used here have previously been reported in Stevens et al. (Reference Stevens, Blass, Zhu, Verzicco and Lohse2018) and we include some relevant details here for completeness. The incompressible Navier–Stokes (NS) equations within the Boussinesq approximation are solved using a second-order finite difference scheme. The dimensionless equations are
\begin{gather} \frac{\partial \tilde{\boldsymbol{u}}}{\partial t} + \tilde{\boldsymbol{u}}\boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\boldsymbol{u}} = -\boldsymbol{\nabla} \tilde{p} + \sqrt{\frac{Pr}{Ra}}\nabla^2 \tilde{\boldsymbol{u}} + \tilde{\theta} \boldsymbol{e}_z, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{u}} = 0, \end{gather}
\begin{gather}\frac{\partial \tilde{\theta}}{\partial t} + \tilde{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\theta} = \frac{1}{\sqrt[]{Ra Pr}}\nabla^2 \tilde{\theta}, \end{gather}
where the characteristic length scale is the height of the RBC cell 
$H$, the free-fall time scale 
$T= \sqrt {H/(g\alpha \varDelta )}$ and the velocity scale is 
$\sqrt {g\alpha \Delta H}$, which is known as the free-fall velocity. Therein, 
$g$ is the gravitational acceleration (acting in the opposite direction to the unit vector 
$\boldsymbol {e}_z$), 
$\alpha$ is the thermal expansion coefficient and 
$\varDelta$ is the temperature difference between the top and the bottom plates. Two dimensionless parameters appear after normalizing the NS equations, namely the Prandtl number 
$Pr=\nu /\kappa$ and the Rayleigh number 
$Ra=\alpha g \Delta H^3/(\nu \kappa )$, where 
$\kappa$ is the thermal diffusivity and 
$\nu$ the momentum diffusivity (kinematic viscosity). Finally, the geometric control parameter is the aspect ratio 
$\varGamma =L/H$ of the square computational domain, with 
$L$ denoting the length of the cell in the periodic 
$x$ and 
$y$ directions. The simulations are performed with the second-order code AFiD (van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015) and parameters of the direct numerical simulation (DNS) are summarized in table 1. Importantly, the BLs are well resolved and comply with the criterion derived in Shishkina et al. (Reference Shishkina, Horn, Wagner and Ching2015). The present dataset was originally used to confirm grid convergence in Stevens et al. (Reference Stevens, Blass, Zhu, Verzicco and Lohse2018) and is therefore even better resolved than the data reported there.
Table 1. Parameters of the dataset. Input parameters are the 
$Ra$ number, the numerical resolution in the horizontal (
$N_x\times N_y$) and wall-normal (
$N_z$) directions, the aspect ratio 
$\varGamma$ and the runtime 
$t$ of the DNS normalized by the free-fall time scale 
$T$. Relevant results are the 
$Nu$ number, the superstructure size 
$\hat {l}$ which is calculated from the coherence spectrum in Krug et al. (Reference Krug, Lohse and Stevens2020), the thermal BL thickness 
$\delta _\theta$ and, finally, the number of grid points that are in the thermal BL 
$N_\theta$.
Instantaneous normalized values of pressure, temperature and the velocity vector are denoted by 
$\tilde {p}$, 
$\tilde {\theta }$ and 
$\tilde {\boldsymbol {u}}$, respectively and we refer to the components of 
$\tilde {\boldsymbol {u}}$ as 
$\tilde {v}$ (horizontal) and 
$w$ (vertical). The operator 
$\langle \cdot \rangle$ is used to indicate averaging over the statistically homogeneous horizontal plane and time. We use the notation 
$\tilde {\xi }(x,y,z,t) = \varXi (z) +\xi (x,y,z,t)$ to decompose an arbitrary instantaneous quantity 
$\tilde {\xi }(x,y,z,t)$ into its mean 
$\varXi (z) = \langle \tilde {\xi } \rangle$ and a fluctuating part 
$\xi (x,y,z,t)$.
3. Results
The essence of this study is best illustrated by the snapshot of the temperature field shown in figure 2(a). There, the colour contour represents 
$\theta _S(z=\delta _\theta )$, i.e. the small-scale temperature fluctuations at BL height. Note that we use subscripts 
$S,L$ to represent the small- (
$k>k_{cut}$) and large-scale (
$k<k_{cut}$) signals, respectively, which are obtained from spectral filtering with cutoff wavenumber 
$k_{cut}$. Even though the large scales – and hence the effect of a direct superposition of the LCS onto fluctuations at BL height observed before (Pandey et al. Reference Pandey, Scheel and Schumacher2018; Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018; Krug et al. Reference Krug, Lohse and Stevens2020) – have been removed, there is still a distinct large-scale ‘organization’ visible for 
$\theta _S$. This pattern stems from spatial variations of properties such as the local amplitude and the local wavenumber of 
$\theta _S$. Figure 2(a) further contains zero crossings of 
$\theta _L$ (dashed lines), which make it obvious that these spatial variations in the small-scale signal are highly correlated to the large-scale dynamics. In the next section, we will introduce a procedure to capture the large-scale ‘footprint’ in small-scale turbulence by extracting local amplitude and wavenumber. On this basis, we can then quantify the correlation of this footprint with the large scales.
Figure 2. (a) Snapshot of the small-scale temperature fluctuations 
$\theta _S$ at a distance 
$z = \delta _\theta$ above the hot plate. Insets show plume emitting (
$\theta _L>0$, blue box) and plume impacting (
$\theta _L<0$, red box) regions. (b,c) Show the large-scale amplitude 
$\sigma _{\theta L}$ and wavenumber 
$\kappa _{\theta L}$, respectively, for the snapshot shown in (a). The boxes indicate the positions of the insets in (a). Dashed black lines in all panels correspond to zero crossings in the large-scale signal 
$\theta _L$ at mid-height 
$z = 0.5$. The scales of the colour axes are 
${\pm }0.36$, 
${\pm }0.2$ and 
${\pm }0.4$ in (a), (b) and (c), respectively.
3.1. Reduced-order representation of small-scale fluctuations
To quantify how the large-scale Rayleigh–Bénard (RB) superstructures interplay with the small-scale turbulence (
$\xi _S$), we have to resort to a simplified representation of the small-scale fluctuations (
$\xi$ is a placeholder for the quantity concerned, e.g. 
$\theta$, 
$w$,…). These fluctuations are broadband (e.g. Krug et al. Reference Krug, Lohse and Stevens2020) and require reducing the wealth of information in 
$\xi _S(x,y)$ to a comprehensible set of representative quantities. Here, this is achieved using a wavelet procedure, which allows us to extract the large-scale spatial variations of the small-scale statistics. A similar procedure was applied to time series of velocity fluctuations by Baars et al. (Reference Baars, Talluru, Hutchins and Marusic2015). However, their methodology is limited to 1-D data and therefore needs to be extended to cope with 2-D fields to be applicable here. Our procedure yields two new variables: (i) the 2-D space-varying magnitude 
$\sigma _\xi$ and (ii) the 2-D space-varying characteristic wavenumber 
$k_\xi$ of the energy in 
$\xi _S(x,y)$. Through correlating these two representative variables of the small-scale turbulence with the large-scale superstructures in 
$\xi _L$, we study the interaction between the two energetic regions that were identified in spectral space (see the discussion alongside figure 1).
Our procedure to obtain the reduced-order representation of the small-scale fluctuations is illustrated in figure 3. A wall-parallel plane of 
$\xi$ forms the starting point in this diagram (top left of figure 3b). For ease of visualization, a 1-D subset taken along the 
$x$-direction of 
$\xi$ is shown in figure 3(a) for one particular 
$y = y_p$ location. Inspection yields a large-scale variation of the fluctuating magnitude 
$\sigma _{\xi L}$ (indicated by the dashed envelope) and a spatial variation in the characteristic wavenumber 
$k_{\xi L}$; this reflects the spatial variation in small-scale properties already observed in figure 2(a).
Figure 3. Diagram of the wavelet-based procedure, used in extracting local fluctuation properties of a turbulent quantity in a plane of RB turbulence, in particular the characteristic energy amplitude of the fluctuations, 
$\sigma _{\xi L}(x,y)$, and the characteristic wavenumber of the fluctuations, 
$k_{\xi L}(x,y)$.
The first step in providing a quantitative measure for these variations involves a 2-D wavelet transform. We restrict the discussion of this transform to a minimum here and refer to appendix A for details and to e.g. Farge (Reference Farge1992) and Addison (Reference Addison2002) for further mathematical intricacies. The wavelet transformation is defined by a convolution of a mother wavelet 
$\varPsi$ with 
$\xi (x,y)$, following
\begin{equation} A_{\xi}(x,y;l_s) = \int\int \xi(x^*,y^*)\bar{\varPsi}\left(\frac{x-x^*}{l_s},\frac{y-y^*}{l_s}\right)\textrm{d} x^*\,\textrm{d} y^*, \end{equation}
where the bar indicates the complex conjugate and 
$A$ is the so-called wavelet coefficient at a length scale 
$l_s$ of the self-similar wavelet. The wavelet 
$\varPsi$ itself is a localized basis function, and is here taken as the isotropic Morlet wavelet (a 2-D harmonic wave with a Gaussian envelope). This mother wavelet has no angular sensitivity, and is a suitable choice for RB turbulence, as there is no angular dependency in the statistics of this flow. The coefficient 
$A_{\xi }(x,y;l_s)$ is then a measure of the local amplitude (at point 
$(x,y)$) at the scale 
$l_s$.
A difficulty arises when relating wavelet decompositions to Fourier transforms because the wavelet scale 
$l_s$ does not have a precise counterpart in wavenumber space. This is because, being localized in space, a wavelet 
$\varPsi$ comprises energy at a range of wavenumbers instead of just a single one. The way to relate the wavelet scale to a wavenumber is then to determine the centre wavenumber of the wavelet's Fourier transform, denoted the ‘pseudo-wavenumber’ 
$k_s$. We will adopt this concept here in order to better connect to other work on RBC, which predominantly employs Fourier transforms. With this, coefficients 
$A_{\xi }(x,y;l_s)$ can be expressed as 
$A_{\xi }(x,y;k_s(l_s))$ and reflect the convolution magnitude as a function of a characteristic wavenumber. Finally, with the proper normalization, the square of the wavelet coefficients measures the energy level of a given field in terms of space and wavenumber, yielding the energy density
\begin{equation} E_\xi(x,y;k_s) = \frac{A_{\xi}(x,y;k_s)^2}{k_s}. \end{equation} The stack of 
$N$ planes in the top right of figure 3(b) represents quantity 
$E_\xi (x,y;k_s)$ computed at 
$N$ different scales. For a given sample location, denoted as 
$(x_i,y_i)$, the local wavelet spectrum 
$E_\xi (x_i,y_i;k_s)$ is shown as a function of 
$k_s$. A total of 
$N = 33$ logarithmically spaced scales were used in the small-scale range 
$k_{cut} < k < k_N$, with 
$k_N=4.5\times 10^3$ denoting the DNS Nyquist wavenumber. The local wavelet spectrum is a measure of the energy spectrum at one spatial location. It is influenced by a domain surrounding 
$(x_i,y_i)$, which is referred to as the cone-of-influence of the wavelet, i.e. the region spanned by the wavelet at a certain scale: a small-scale wavelet is more localized than a larger-scale wavelet. Since the boundaries of the computational domain are periodic, localized wavelet spectra exist for all points in the 
$(x,y)$ domain. When all localized spectra are averaged, one global wavelet spectrum is found for the full 
$(x,y)$ field which then matches the Fourier spectrum at one 
$z$ height in figure 1. So, the key feature of 
$E_\xi (x,y;k_s)$ is that it preserves spatial information, which is absent in Fourier spectra.
So far, the decomposition of 
$\xi (x,y)$ into energy densities 
$E_\xi (x,y;k_s)$ at 
$N$ different scales has increased the dimensions of the system. For a reduced-order representation of the small-scale turbulence, we now define the following quantities: first, we introduce a spatial field that captures the variation in the energy magnitude of 
$\xi (x,y)$. Integrating 
$E_\xi (x,y;k_s)$ over the small-scale wavenumber range 
$k_{cut} < k < k_N$ and taking the square root yields the space-varying standard deviation of 
$\xi _S(x,y)$ according to
\begin{equation} \varSigma_\xi + \sigma_\xi(x,y) = \left(\int_{k_{cut}}^{k_N} E_\xi(x,y;k_s)\,\textrm{d} k_s\right)^{0.5}. \end{equation}
Here, 
$\varSigma _\xi$ denotes the mean standard deviation of the small-scale fluctuations (i.e. 
$\varSigma _\xi = \langle \xi _S^2 \rangle ^{{1}/{2}}$), whereas 
$\sigma _\xi (x,y)$ is the zero-mean spatial fluctuation of the standard deviation. Only the large-scale variation of 
$\sigma _\xi$ is retained by low-wavenumber pass filtering with 
$k < k_{cut}$. The resulting spatial field 
$\sigma _{\xi L}(x,y)$ is interpreted as an envelope to the small-scale fluctuations and is plotted in figure 3(a). A second representative property of 
$\xi (x,y)$ is its local wavenumber, which we take as the first spectral moment of the local wavelet spectrum,
\begin{equation} K_\xi + k_\xi(x,y) = \frac{\displaystyle\int_{k_{cut}}^{k_N} kE_\xi(x,y;k_s)\,\textrm{d} k_s}{\displaystyle\int_{k_{cut}}^{k_N} E_\xi(x,y;k_s)\,\textrm{d} k_s} = \frac{1}{\left(\varSigma_\xi + \sigma_\xi(x,y)\right)^2} \int_{k_{cut}}^{k_N} kE_\xi(x,y;k_s)\,\textrm{d} k_s. \end{equation}
Similarly to the standard deviation, the mean wavenumber of the small-scale energy in 
$\xi (x,y)$ is 
$K_\xi$ and its zero-mean fluctuations are denoted as 
$k_\xi (x,y)$. Low-wavenumber pass filtering the latter yields 
$k_{\xi L}(x,y)$. In summary, the raw small-scale field 
$\xi _S(x,y)$ is represented by 
$\sigma _{\xi L}(x,y)$, which reflects the large-scale spatial deviation of the energy (or intensity) in the small-scale fluctuations from its mean energy. Further, 
$k_{\xi L}(x,y)$ describes the deviation of characteristic local wavenumber from the mean wavenumber. Thus, positive (negative) 
$k_{\xi L}$ means that the small-scale fluctuations are of a higher (lower) wavenumber locally, compared to an average over the entire domain. Results of 
$\sigma _{\theta L}$ and 
$k_{\theta L}$ computed for the snapshot in figure 2(a) are included in figure 2(b,c), respectively. To what degree these quantities are correlated with the large-scale signal can be gauged by comparing to the zero crossings of 
$\theta _L$ (shown as black lines). We will quantify this aspect next.
3.2. Correlations of large and small scales
To study the influence of the large scales on the small scales in RBC, we define an amplitude modulation coefficient similar to that used in work on turbulent boundary layers (e.g. Bandyopadhyay & Hussain Reference Bandyopadhyay and Hussain1984; Mathis et al. Reference Mathis, Hutchins and Marusic2009). This coefficient is the normalized correlation between the large-scale envelope of the small-scale turbulence 
$\sigma _{\xi ,L}(x,y,z)$ and the large-scale signal of the temperature at mid-height 
$\theta _L(x,y,H/2)$,
\begin{equation} R_{\xi,\sigma}(z) = \frac{\langle \sigma_{\xi L} (x,y,z) \theta_L(x,y,H/2) \rangle}{\sqrt{\langle \sigma_{\xi L}(x,y,z)^2\rangle} \sqrt{\langle \theta_L(x,y,H/2)^2\rangle}}. \end{equation}
It is important to emphasize that, in this definition, we always use the large-scale temperature signal at mid-height 
$\theta _L(x,y,H/2)$ as reference. Clearly, other choices are possible here. However, since the LSC is known to be highly coherent along the vertical direction and because 
$\theta _L$ correlates almost perfectly with the large-scale velocity signal (Krug et al. Reference Krug, Lohse and Stevens2020), this will have no significant impact on our results and conclusions. By definition, 
$R_{\xi ,\sigma }$ varies between 
$-1$ (perfect anti-correlation) and 
$1$ (perfect correlation and 
$R_{\xi ,\sigma } = 0$ implies uncorrelated signals.
Analogously, the wavenumber modulation coefficient 
$R_{k,\xi }(z)$ is defined by correlating the large-scale part of the wavenumber fluctuations 
$k_{\xi L}(x,y,z)$ to 
$\theta _L(x,y,H/2)$,
\begin{equation} R_{\xi,k}(z) = \frac{\langle k_{\xi L} (x,y,z) \theta_L(x,y,H/2) \rangle}{\sqrt{\langle k_{\xi L}(x,y,z)^2 \rangle} \sqrt{\langle \theta_L(x,y,H/2)^2\rangle}}. \end{equation}
We note that (3.5) and (3.6) define a normalized correlation. The physical relevance of such a correlation, however, also depends on how much (absolute) energy is contained in the fields of 
$\sigma _{\xi ,L}$ and 
$k_{\xi L}$. For the snapshot in figure 2, fluctuations in both 
$\sigma _{\theta ,L}$ and 
$k_{\theta L}$ are of the same order (tens of per cent) as their respective means and hence significant. We will revisit this point in § 3.3 in more detail.
Results for both modulation coefficients in the temperature field are presented as a function of the height 
$z$ in figure 4(a). The most striking observation from this plot is that the small-scale amplitude of 
$\theta$ is highly correlated to the large scales up to a distance of approximately 
$2\delta _\theta$ above the plate. The peak of this correlation at which remarkably 
$R_{\theta ,\sigma }<-0.9$ is found at 
$z\approx \delta _\theta$. This coincides with the location of the strongest small-scale fluctuations in the temperature field as is obvious from the right-hand plot of figure 4(a) as well as from the spectral peak of the energetic small scales in figure 1. Throughout this paper, all statistics are plotted with reference to the hot plate. Hence, the negative sign of 
$R_{\sigma ,\theta }$ in the near-wall region implies that in plume impacting regions (where 
$\theta _L<0$) the local small-scale energy is higher (
$\sigma _{\theta L}>0$). This is also confirmed visually from the insets in figure 2(a).
Figure 4. Wall-normal profiles of the amplitude modulation coefficient (solid dark blue), the wavenumber modulation coefficient (dash-dotted orange) and the small-scale energy 
$\overline {\sigma _\psi }^2$ (solid black) for (a) the temperature 
$\theta$, (b) the wall-normal velocity 
$w$ and (c) the small-scale turbulent heat transport term 
$(\theta _S w_S)$. The blue dotted line in (c) is the correlation coefficient for 
$(\theta _S w_S)_L$. (d) Illustrates the construction of the 
$(\theta _S w_S)$ signal. It contains a large-scale part (solid red), a non-zero mean (dashed black line) and the remaining small-scale signal, of which the respective variations with 
$z$ are shown in the right plot of panel (c).
Interestingly, 
$R_{\theta ,\sigma }$ crosses zero around 
$z=2\delta _\theta$ and attains somewhat lower, yet distinctly non-zero, positive values around 
$R_{\theta ,\sigma }\approx 0.4$ in the bulk. The subsequent drop to zero at mid-height is a symmetry requirement. The cross-over height of 
$z=2\delta _\theta$ is consistent with recent findings by He & Xia (Reference He and Xia2019), who found similar levels of temperature fluctuations in plume emitting and sheared regions at this height. Further away from the wall, He & Xia (Reference He and Xia2019) found that plume emitting regions contain significantly higher fluctuations of 
$\theta$, which is also consistent with the present observations.
Compared to 
$R_{\theta ,\sigma }$, the wavenumber modulation coefficient for 
$\theta$ (dashed red line) remains relatively low at all 
$z$. Also in this case, the strongest correlation encountered within the thermal BL with 
$R_{\theta ,k}\approx 0.3$. Unlike 
$R_{\theta ,\sigma }$, 
$R_{\theta ,k}$ does not change sign and remains positive throughout, which implies that structures in the temperature field are slightly more ‘compact’ in plume emitting regions.
In principle, correlation coefficients can be evaluated for all fluctuating quantities. Here, we will, in addition to 
$\theta$, only consider the wall-normal velocity 
$w$. This is motivated by the fact that 
$w$ also features in the turbulent heat transport 
$\theta w$. From figure 4(b), we find that the amplitude modulation coefficient for 
$w$ is negative at all heights at intermediate correlation levels of approximately 
$0.4$. On the contrary, 
$R_{w,k}$ is positive at all 
$z$ locations and reaches values up to 0.5. Both, 
$R_{w,\sigma }$ and 
$R_{w,k}$ peak in the range 1–10 
$\delta _\theta$, over which also the fluctuation level in 
$w_S$ builds up to its bulk value (see right-hand plot in figure 4b). The picture that emerges here is that 
$w_S$ fluctuations are of lower magnitude but higher wavenumber in plume emitting regions.
Besides the individual fields of 
$w_S$ and 
$\theta _S$, it is meaningful to study the modulation effects for the turbulent heat transport 
$\theta w$. Initially, we apply the same analysis leading to (3.5) and (3.6) also to 
$\theta w$. The corresponding results, 
$R_{\theta w,\sigma }$ (solid blue line) and 
$R_{\theta w,k}$ (dashed red), characterize zero-mean fluctuations of 
$(\theta w)_S$ and are shown in figure 4c. Both display substantial correlation with the large scales up to BL height, at which both the mean and the fluctuations in the normalized small-scale heat transport already reach the level of 
$Nu$, as can be seen from the right plot in 4(c). Consistent with observations for 
$\theta$ and 
$w$ individually, we find 
$R_{\theta w,\sigma }<0$ and 
$R_{\theta w,k}>0$ for 
$z\lessapprox \delta _\theta$, while at larger wall distances 
$R_{\theta w,\sigma }$ is slightly positive with no significant wavenumber modulation.
There is, however, another point to be made here. This relates to the fact that even though individually 
$\theta _S$ and 
$w_S$ are zero mean and small scale only, their product 
$\theta _Sw_S$ contains not only a significant mean (that contributes the bulk of the total heat transport above 
$z \approx \delta _\theta$, see figure 4(c), right plot), but also large-scale fluctuations. This is illustrated for a 1-D sample in figure 4(d). The large-scale variations 
$(\theta _S w_S)_L$ are not equally strong as those at smaller scales (note the factor 
$5$ for the rms of 
$(\theta _S w_S)_L$ in figure 4c), yet their standard deviation still amounts to approximately 10 % of 
$Nu$ around BL height, which is significant. At this wall distance, 
$(\theta _S w_S)_L$ displays a pronounced negative correlation with 
$\theta _L$ (dotted blue line in figure 4c), which extends even somewhat deeper into the bulk compared to 
$R_{\theta w,\sigma }$. This implies that the local heat transport carried by small-scale fluctuations is locally reduced (increased) in plume emitting (impacting) regions around BL height.
Correlation coefficients are very useful in quantifying the overall nature of scale interactions. However, they do not reveal all details about the spatial organization in the flow. We will therefore complement the results discussed so far by additionally considering conditioned averages in the following.
3.3. Conditioned averages
3.3.1. Averaging procedure
The results obtained so far indicate that the properties of the small-scale turbulence, especially in the vicinity of the BLs, vary significantly depending on the location within the LSC. In this section, we aim to investigate those spatial variations in more detail. To that end, we employ conditioned averages with respect to the zero crossing of the large-scale signal 
$\theta _L$. The definition of the new conditioned coordinate 
$d$ is illustrated in figure 5(a). For every point 
$P$ in a plane 
$d^*$ measures the smallest distance to a contour of 
$\theta _L =0$, i.e. 
$d^*$ is the smallest radius of a circle around 
$P$ that just touches a point on the zero contour. To distinguish large-scale positive from large-scale negative regions, 
$d$ then takes the sign of 
$\theta _L$, such that
\begin{equation} d=\textrm{sign}(\theta_L)d^*. \end{equation}
The 
$d$-field corresponding to the 
$\theta _L$ snapshot in figure 5(a) is shown in panel (b) of the same figure. The blow-up in figure 5(b) presents the large-scale horizontal velocity vector at BL height (with the same cutoff wavenumber as 
$\theta _L$). This plot demonstrates that the LSC flow at the hot plate is generally directed from small 
$d$ to large 
$d$. The largest large-scale velocity magnitudes are encountered around 
$d\approx 0$. Since the large-scale circulation is approximately aligned with 
$\boldsymbol {\nabla } d$, it is useful to decompose the horizontal velocity vector 
$\boldsymbol {v}$ into a component 
$v_p$ along 
$\boldsymbol {\nabla } d$ and one (
$v_n$) normal to it as illustrated in the blow-up of figure 5(a).
Figure 5. Illustration of the definition of the conditioned average. (a) Large-scale part of the temperature fluctuations at mid-height 
$\theta _L(z = 0.5)$. Black lines in all panels are iso-contours of 
$\theta _L =0$. Every point 
$P$ in the plane is assigned a distance 
$d$, where the magnitude of 
$d$ is determined by the distance to the nearest point on the iso-contour (see blow-up in (a)) and the sign is determined by that of 
$\theta _L$ at point 
$P$ . The resulting 
$d$-field corresponding to the snapshot of 
$\theta _L$ in (a) is shown in (b). The in-plane velocity vector at each point (see blow-up in (b) is decomposed into a component along (
$v_p$) and normal (
$v_n$) to the gradient 
$\boldsymbol {\nabla } d$ as shown in the blow-up in (a).
To retain the large-scale organization in the flow, we then compute conditional averages by taking the mean overall points with constant 
$d$. Mathematically, this results in a triple decomposition for any instantaneous quantity 
$\tilde {\xi }(x,y,z,t)$ according to
\begin{equation} \tilde{\xi}(x,y,z,t) = \varXi(z)+ \bar{\xi}(z,d) + \xi'(x,y,z,t), \end{equation}
where the overline denotes the conditional average, and 
$\xi '(x,y,z,t)$ is the fluctuation about the conditional and the temporal average.
Conceptually, the conditioned mean 
$\bar {\xi }(z,d)$ corresponds to the large-scale fluctuation field and represent the LSC, while the fluctuations 
$\xi '(x,y,z,t)$ around the conditioned mean are related to small-scale turbulence. It is hence of interest to analyse the spatial distribution of statistics related to 
$\xi '$ in the present context. Before doing so, we briefly present some results for the conditioned mean in order to orient the reader about basic properties of the LSC in the present configuration. For that purpose, the contour of the conditionally averaged temperature 
$\bar {\theta }$ is shown in figure 6(a) overlaid with vectors representing 
$(\bar {v}_p, \bar {w})$. It is clear from this figure that the ranges 
$d<-1$ and 
$d>1$ correspond to the plume impacting and emitting zones of the LSC, respectively, and the strongest horizontal flow is encountered for 
$d \approx 0$. Further, the histogram in figure 6(b) shows that the distribution of 
$d$ is already well converged from a single snapshot only with only few counts beyond half the typical superstructure size 
$\hat {l} \approx 3$ for this case. Relevant transport quantities that can be derived from the fields shown in figure 6(a) include the local wall shear stress 
$\bar {\tau }_w = |(\partial _z \overline {v_p}|_{z=0})|$, which is seen to peak at slightly negative 
$d$ (see figure 6c). Moreover, we observe a significant spatial variation of 
$\pm 0.3Nu$ for the local heat transport 
$\overline {Nu} = \partial _z(\varTheta +\bar {\theta })|_{z=0}$. As figure 6(c) shows, heat transport is enhanced in plume impacting regions and decreased in plume emitting regions. These observations are similar to those reported for RBC at 
$\varGamma =1$, 
$Pr=0.79$ and matching 
$Ra$, see Wagner, Shishkina & Wagner (Reference Wagner, Shishkina and Wagner2012). For a more detailed study on large-scale properties including their 
$Ra$ dependence, we refer to Blass et al. (Reference Blass, Verzicco, Lohse, Stevens and Krug2021).
Figure 6. (a) Superstructure convection roll as obtained from the conditional averaging procedure. The colour contours represent 
$\bar {\theta }$ and span the range 
${\pm }0.1$. Vectors indicate the field 
$(\bar {v}_p, \bar {w})$ with 15 data points skipped between adjacent vectors. (b) Number of points in a single plane as a function of the conditioning parameter 
$d$ for a total of 
$13$ snapshots. (c) Variation of the conditioned Nusselt number 
$\overline {Nu}$ and of the conditioned shear stress 
$\bar {\tau }_w$ as functions of 
$d$.
3.3.2. Conditioned variance
Several statistical properties are relevant when examining the fields of conditioned fluctuations 
$\xi '$ in terms of their magnitude. The most basic one is the variance 
$\langle \xi '^2\rangle (z)$, which indicates how energetic the small scales are at a given height (note that 
$\langle \cdot \rangle$ continues to imply horizontal averaging, in this case over 
$d$). The small-scale temperature variance 
$\langle \overline {\theta '^2}\rangle (z)$ is shown in figure 7(a) and indeed resembles 
$\langle \theta _S^2\rangle$ in figure 4(a) closely. The main interest here is in lateral (along 
$d$) variations of the variance, which we highlight by considering the difference 
$\Delta \overline {\xi '^{2} }(z,d) = \overline {\xi '^{2}}(z,d) - \langle \overline {\xi '^{ 2}}\rangle (z)$. This quantity is akin to the large-scale envelope 
$\sigma _{\theta L}$. We plot the normalized standard deviation of 
$\Delta \overline {\theta '^{2} }$ in figure 7(c) as a measure of how much the temperature variance varies with 
$d$ at a given height. Finally, we define the coefficient
\begin{equation} \rho_\xi(z,d) \equiv \frac{\overline{\xi'^{2}}(z,d) - \langle \overline{\xi'^{ 2}}\rangle(z)} {\langle (\overline{\xi'^{2}}(z,d) - \langle \overline{\xi'^{2}}\rangle(z))^2\rangle^{1/2}} = \frac{\Delta \overline{\xi'^{2}}(z,d) }{\textrm{{std}}\left( \Delta \overline{\xi'^{2} }\right)(z)}, \end{equation}
to illustrate how 
$\Delta \overline {\xi '^{2} }$ varies across the LSC. Results for 
$\rho _\theta$ are displayed in figure 7(b). Entirely consistent with figure 4(a), also this plot indicates enhanced (reduced) small-scale temperature fluctuations in the impacting (emitting) region below 
$z \approx 2\delta _\theta$. At this height, both 
$\langle \overline {\theta '^2} \rangle$ and 
$\Delta \overline {\theta '^{2} }$ are of significant magnitude. Also, the change in sign for 
$R_{\theta \sigma }$ towards positive values observed for 
$z \gtrsim 2\delta _\theta$ is echoed in the distribution of 
$\rho _\theta$. Generally, 
$\rho _\theta >0$ at 
$d>0$ (and vice versa) in the bulk. However, the zero crossing clearly does not coincide with 
$d =0$ in this case, which is testament to the additional information provided by considering the conditioned averages.
Figure 7. Statistics of the conditioned fluctuations of 
$\theta$ (a–c) and 
$w$ (d–f). Additionally, statistics of 
$(\theta w)'$ are shown in panels (g–i). Note that the standard deviation operator used in (c,f,i) implies horizontal averaging only, i.e. 
${\textrm {std}}(\cdot ) = \langle (\cdot )^2 \rangle ^{1/2}$ . The horizontal dashed line in (b,e,h) shows the local thermal BL height 
$\bar {\delta }_\theta (d) = 1/(2 \overline {Nu})$. The colour bar applies to all contour plots.
Results analogous to those presented for the temperature field in figure 4(a–c) are included for 
$w$ in panels (d–f) of the same figure. Also in this case, there is a striking agreement between the distribution of 
$\rho _w$, which indicates stronger fluctuations on the plume impacting side, and 
$R_{w \sigma }$ in figure 4(b).
Finally, we turn again to the turbulent heat transport. Here, we define the small-scale contribution according to
\begin{equation} (\theta w)'(z,d) = \overline{\theta w}(z,d) - \bar{\theta}\bar{w}(z,d), \end{equation}
that is, we subtract the part of the heat transport carried by the LSC from the overall turbulent transport. As figure 7(g) shows, 
$\langle (\theta w)' \rangle$ carries a bulk of the global transport at wall distances beyond the thermal BL height (note also the agreement with 
$\langle \theta _S w_S \rangle$ in figure 4c). The corresponding coefficient
\begin{equation} \rho_{\theta w}(z,d) \equiv \frac{(\theta w)'(z,d) - \langle (\theta w)' \rangle (z)} {\langle ((\theta w)'(z,d) - \langle (\theta w)' \rangle_d(z))^2\rangle^{1/2}} = \frac{\varDelta (\theta w)'(z,d)}{{\textrm{std}}\left( (\theta w)' \right)(z)}, \end{equation}
is analogous to (3.9) but provided explicitly for clarity. Importantly, it is clear from (3.10) and (3.11) that 
$(\theta w)'$ and hence also 
$\rho _{\theta w}$ describe how the net heat transport by the small scales varies spatially. The results presented here therefore approximately correspond to 
$(\theta _S w_S)_L$ in figure 4(c). The contour plot of 
$\rho _{\theta w}$ (figure 7h) then also confirms the observation that up to several BL heights, small-scale heat transport is stronger in plume impacting regions. In the bulk, the correlation coefficient for 
$(\theta _S w_S)_L$ is essentially zero (see dotted blue line in figure 4c). However, this does not imply that the small-scale heat transport is uniform there: as figure 7(h) reveals, small-scale transport is enhanced for both small and large 
$d$ and lower in the centre of the LSC around 
$d=0$. Even though a pronounced large-scale organization evidently exists, this pattern does not lead to a non-zero correlation coefficient due to its reflection symmetry around 
$d=0$. The horizontal standard deviation 
$\varDelta (\theta w)'$ amounts to typically 10 % to 20 % of 
$\langle (\theta w )\rangle$ (see figure 7i). These figures are comparable in terms of 
$Nu$ (cf. figure 7g), such that spatial variations of the small-scale heat transport are indeed significant on a global level.
3.3.3. Conditioned wavenumber
As a final point in our analysis, we reconsider spatial variations of the local wavenumber defined in (3.4). We normalize the conditioned wavenumber fluctuations 
$\bar {k}_\xi (z,d)$ by their mean 
$K_\xi (z)$ at the respective height to obtain
\begin{equation} \kappa_\xi(z,d) = \frac{\bar{k}_\xi(z,d)}{K_\xi (z)}. \end{equation}
The results for 
$\kappa _\theta$, 
$\kappa _w$ and 
$\kappa _{\theta w}$ are presented in figure 8(a,c,e), respectively, along with the corresponding means 
$K_\theta$, 
$K_w$ and 
$K_{\theta w}$ in figure 8(b,d,f) for reference.
Figure 8. Normalized conditioned wavenumber fluctuations of different fields; (a) temperature, (c) wall-normal velocity and (e) heat transport. The respective means as a function of 
$z$ used for normalization are shown in (b,d,f). The colour bar applies to all contour plots.
Largely, the plots of 
$\kappa _\theta$, 
$\kappa _w$ recount the observations made in figure 4(a,b), indicating higher wavenumbers on the plume emitting side for both quantities. This trend appears to be slightly more pronounced for 
$w$ compared to 
$\theta$.
The strongest wavenumber modulation is observed for 
$\kappa _{\theta w}$ in figure 8(e). The trend there agrees with 
$\kappa _\theta$ and 
$\kappa _w$ (and also with 
$R_{\theta w,k}$ in figure 4c) in indicating higher wavenumbers at 
$d>0$ approximately up to BL height 
$\delta _\theta$. Similar to 
$\rho _{\theta w}$, in the bulk a pronounced modulation pattern is also visible for 
$\kappa _{\theta w}$, which remains hidden in 
$R_{\theta w,k}$ due to its symmetry. In the case of 
$\kappa _{\theta w}$, this pattern is reversed for 
$\delta _\theta \lessapprox z \lessapprox 10\delta _\theta$, for which positive 
$\kappa _{\theta w}$ is found in the centre instead of at large 
$|d|$.
4. Concluding remarks
Our analysis on the basis of correlations as well as through conditioned averaging has revealed significant large-scale organization in small-scale turbulence fields in RBC. These modulation effects were seen to be strongest within – but by no means limited to – a layer close to the wall. Changes in the small-scale organization patterns occurred around a wall height of 
$z \approx 2\delta _\theta$. At this wall distance, the mean temperature profile approaches the bulk value and the location also coincides with the boundary of the ‘mixing zone’ described by Wang et al. (Reference Wang, He and Tong2019). They find that this mixing zone is characterized by the emergence of intermittent, strong, hot plumes, that lead to long exponential tails in the probability distribution function (PDF) of the temperature. In their case, the mixing region differs from the ‘boundary region’, where the temperature PDF is described by homogeneous, purely Gaussian, background fluctuations. This distinction does not appear to be relevant in the present context. Further, the cross-over height of 
$z \approx 2\delta _\theta$ links to a sharp decrease in vertical coherence with reference at BL height of the temperature field at large wavenumbers observed by Krug et al. (Reference Krug, Lohse and Stevens2020). This is indicative of the fact that, indeed, small structures emanating from the BL do not extend beyond this height, which may explain why the small-scale organization changes there.
In RBC, the small scales are typically most energetic close to the wall, making our results most relevant in this region; the strongest correlations are also encountered here with coefficients up to 0.9 for the temperature amplitude. Our findings regarding the near-wall zone are summarized graphically in figure 9. Up to a wall distance of 
$z \approx 2\delta _\theta$, we find that the amplitude of small-scale structures (indicted as plumes in figure 9) in the plume impacting region is higher, while the wavenumber (and hence the structure size) is higher there. Most importantly, we find that the turbulent heat transport carried by the small scales is higher (by about 10 % of 
$Nu$) in the plume impacting region and decreased by the same amount on the emitting side. Significant large-scale organization of the small scales was also identified in the turbulent bulk of the flow. However, the patterns are more complicated in this case and vary between different quantities, which precludes a straightforward (graphical) summary of the findings. Among the most notable results concerning the bulk flow is the fact that small-scale transport is enhanced at both sides of the LSC and is minimal in the centre around 
$d = 0$.
Figure 9. Schematic of the large-scale footprint in small-scale near-wall RBC turbulence. Small-scale structures are shown as plumes with their height indicating the local magnitude and variations in spacing the local size (i.e. the wavenumber) of the fluctuations. The qualitative mean temperature profile serves to put the different regions into perspective.
Thus far, we have not touched on the cause of the spatial variations in the small-scale statistics. The vertical localization of the peak in 
$\langle \theta _S^2\rangle$ correlates with the local maximum in the production term of temperature variance, which is proportional to 
$\theta w \partial _z \varTheta$ (see e.g. Deardorff & Willis Reference Deardorff and Willis1967; Togni et al. Reference Togni, Cimarelli and De Angelis2015). It appears therefore likely that also horizontal differences in 
$\langle \theta _S^2\rangle$ in the BL region might arise from locally varying production rates. Indeed, 
$\overline {Nu}$, and hence also the temperature gradient within the BL, is highest for 
$d<0$, which coincides with the strongest amplitude in 
$\langle \theta _S^2\rangle$. Statistics of 
$\theta$ and 
$w$ are interconnected since temperature fluctuations drive production of 
$\langle w^2\rangle$. Wall-normal velocity fluctuations then feed back to temperature fluctuations since the term 
$\theta w$ also features in the production term for 
$\langle \theta _S^2\rangle$. This picture is consistent at least with our results for 
$\rho _\theta$, 
$\rho _w$ as well as for 
$\rho _{\theta w}$ in the near-wall region. To what extent this theory actually applies will have to be evaluated in future studies. One way to address this would be to study the 
$Ra$ dependence of the effects described here. Recent studies in 2-D (Zhu et al. Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018) and 3-D RBC (Blass et al. Reference Blass, Verzicco, Lohse, Stevens and Krug2021) indicate that the 
$\overline {Nu}$ distribution changes with increasing 
$Ra$. This should lead to changes in the small-scale patterns if the above is true. A complete answer is likely even more complicated since the different scalings identified in He & Xia (Reference He and Xia2019) are linked to different leading-order balances (inertia–viscous vs. inertia–buoyancy) (Adrian Reference Adrian1996). Naturally, such a fundamental difference is also expected to have an impact on the quantities studied here.
A further relevant question is how the strength of the scale interaction varies as a function of 
$Ra$. Experience from flat plate turbulent boundary layers shows that in this case the modulation of small scales by the large superstructure scales intensifies with increasing shear Reynolds number 
$Re_\tau$ (e.g. Mathis et al. Reference Mathis, Hutchins and Marusic2009). It remains to be seen if a similar trend exists in RBC. This will also shed light on whether a ‘universal small-scale signal’ exists in pure convection. That is, if the properly normalized statistics of the small-scale signal at a given 
$z$ are invariant for different locations within the LSC and over at least a limited range of 
$Ra$ once they are stripped of the large-scale footprint by ‘demodulation’. Identifying such a universal signal will be crucial in representing small-scale turbulence in RBC effectively.
Finally, it would be interesting to uncover the 
$Pr$-dependence of the effects discussed here. Changes in 
$Pr$ affect the small-scale signal as well as the LSC, rendering it difficult to make any predictions on their interaction based on the present results.
Acknowledgements
We are grateful to R. Stevens for making the data available to us and thank R. Verzicco and D. Lohse for valuable feedback and support.
Funding
This project is supported by the Priority Programme SPP 1881 Turbulent Superstructures of the Deutsche Forschungsgemeinschaft and by the University of Twente Max-Planck Center for Complex Fluid Dynamics. Part of the work was carried out on the national e-infrastructure of SURFsara, a subsidiary of SURF cooperation, the collaborative ICT organization for Dutch education and research.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Details of the wavelet transform used in § 3.1
In decomposing the 2-D spatial field 
$\xi (x,y)$ in space-scale space, the 2-D continuous wavelet transform algorithm, , of the Matlab
Wavelet Toolbox was utilized. This algorithm efficiently performs the convolution (3.1) in the wavenumber domain, and is repeated for a range of spatial wavelet scales 
$l_s$. Following Plancherel's theorem the equivalent of the operation defined in (3.1) in wavenumber space is given by (for an isotropic wavelet, thus without angular dependency)
\begin{equation} A_{\xi}(x,y;l_s) = \int\int \hat{\xi}(k_x,k_y)l_s\hat{\bar{\varPsi}}\left(l_sk_x,l_sk_y\right)\textrm{d}k_x\,\textrm{d}k_y, \end{equation}
where the hat signifies the Fourier transform, and 
$k_x$ and 
$k_y$ are the horizontal wavenumbers, which are related to the radial wavenumber by 
$k \equiv \sqrt {k_x^2 + k_y^2}$. An isotropic and admissible Morlet wavelet was chosen, defined in wavenumber space as
\begin{equation} \hat{\varPsi}\left(k_x,k_y\right) = \exp({-0.5(\vert k_x+\textrm{i} k_y\vert - k_0)^2}), \end{equation}
with 
$k_0 = 6$. Alternative choices of isotropic wavelets (such as the Mexican hat wavelet) with different resolution characteristics in wavelet scale vs. physical space were found to not affect the conclusions made throughout the manuscript. This is related to the fact that the wavelet transform results are used in terms of the (large-scale) variation in the small-scale energy, being 
$E_\xi (x,y;k_s > k_{cut})$. Here, the scale separation between the large and small scales is large (recall figure 1 and its discussion) and is thus unaffected by a wavelet transform resulting in more or less resolution in wavelet scale or physical space, respectively (or vice versa).
Open access
