Temporal dynamics of large-scale structures for turbulent Rayleigh–Bénard convection in a moderate aspect-ratio cylinder

Abstract We investigate the spatial organization and temporal dynamics of large-scale, coherent structures in turbulent Rayleigh–Bénard convection via direct numerical simulation of a 6.3 aspect-ratio cylinder with Rayleigh and Prandtl numbers of $9.6\times 10^7$ and $6.7$, respectively. Fourier modal decomposition is performed to investigate the structural organization of the coherent turbulent motions by analysing the length scales, time scales and the underlying dynamical processes that are ultimately responsible for the large-scale structure formation and evolution. We observe a high level of rotational symmetry in the large-scale structure in this study and that the structure is well described by the first four azimuthal Fourier modes. Two different large-scale organizations are observed during the duration of the simulation and these patterns are dominated spatially and energetically by azimuthal Fourier modes with frequencies of 2 and 3. Studies of the transition between these two large-scale patterns, radial and vertical variations in the azimuthal energy spectra, as well as the spatial and modal variations in the system's correlation time are conducted. Rotational dynamics are observed for individual Fourier modes and the global structure with strong similarities to the dynamics that have been reported for unit aspect-ratio domains in prior works. It is shown that the large-scale structures have very long correlation time scales, on the order of hundreds to thousands of free-fall time units, and that they are the primary source for a horizontal inhomogeneity within the system that can be observed during a finite, but a very long-time simulation or experiment.

convection is considered an ideal problem for investigating the complex phenomenon of turbulent thermal convection because the simple boundary conditions make it more manageable to study without sacrificing the core complexity of thermal convection. The canonical form of turbulent RBC is defined by a domain with a fixed height that extends infinitely in the horizontal directions, thus creating a flow field that is statistically homogenous in the horizontal direction. The infinite domain makes mathematical analysis more tractable, but is obviously not achievable in experiments or discrete computational analysis. Thus, the aspect ratio (Γ ), or the ratio of the horizontal and vertical length scales is a parameter of choice in experimental or numerical investigations. Smaller aspect-ratio domains are easier to simulate numerically or access during the measurements. However, many RBC applications are well represented by wide horizontal layers (Adrian, Ferreira & Boberg 1986), and so the large Γ cases are of significant interest to several different scientific communities.
Rayleigh-Bénard convection in large Γ domains have been primarily studied to investigate the structure and behaviour of pattern formation and a spiral defect chaos at the onset of convection (Meyer, Ahlers & Cannell 1987;Bodenschatz, Pesch & Ahlers 2000). Recently, attention has been turned to studying Rayleigh-Bénard convection in a fully turbulent regime with large (Fernandes & Adrian 2002;Du Puits, Resagk & Thess 2007;Hardenberg et al. 2008;Xia, Sun & Cheung 2008;Bailon-Cuba, Emran & Schumacher 2010;Sakievich, Peet & Adrian 2016), and very-large (Hartlep, Tilgner & Busse 2003;Pandey, Scheel & Schumacher 2018;Stevens et al. 2018;Krug, Lohse & Stevens 2019) aspect-ratio domains. Studies reaching aspect ratios as large as Γ = 128 were able to comment on the size of structures not influenced by the boundary conditions (Stevens et al. 2018), and the natural sizes of these structures were reported to be of six to seven times the height of the domain, consistent with the previous studies of Hartlep et al. (2003) with varying Pr numbers. Despite these emerging studies devoted to wider aspect-ratio cells, there are still a large number of gaps in the community's knowledge of the aspect-ratio affect on the flow field. This is because the main body of knowledge of RBC is still drawn from smaller, primarily unit aspect-ratio studies, which constitute the vast majority of efforts in the field.
One such gap is how Γ affects the structure of the flow field because the organization of structures within the flow field changes dramatically as the domain size is increased (Du Puits et al. 2007;Sakievich et al. 2016;Pandey et al. 2018;Stevens et al. 2018). These structures contain a large portion of the field's kinetic and thermal energy and play an important role in the transport of these quantities. Furthermore, understanding the flow field structure creates a fundamental framework for reasoning and comprehension of the physics.
A prime example is unit Γ turbulent RBC. Much of the work in turbulent RBC over the last several decades has been focused on the unit Γ case wherein a single large-scale circulation (LSC) dominates the flow structure. Identification and dissection of this structure has led to important theories, models and scaling correlations.
One of the first major studies to focus on the LSC was the pioneering work of Krishnamurti & Howard (1981) who observed a persistent, large-scale circulation in unit Γ RBC experiments when the Rayleigh number was greater than 1 × 10 6 . Repetition of the experiments showed that the direction of the LSC changed between realizations, but the LSC was a consistent structure in the flow field. The LSC has often been referred to as the 'wind of turbulence', or a 'mean wind', and this puzzling phenomenon has garnered much attention over the last several decades. The experimental study by Zocchi, Moses & Libchaber (1990) showed that the LSC is generated by the organization of small-scale plumes that gather and cross the layer depth on opposing sides of the convection cell.
Plumes that form in the central region of the cell are generally swept up by the momentum of the LSC leading to a behaviour that is similar to a canonical boundary layer. This similarity was harnessed by Grossmann & Lohse (2000, 2001 to derive a semi-empirical theory for predicting the scaling of heat transfer and Reynolds number as a function of Rayleigh and Prandtl numbers by assuming that the LSC generates a Prandtl-Blasius boundary layer. Another field of interest that has stemmed directly from the LSC is its dynamics. The experimental work by Cioni, Ciliberto & Sommeria (1997) was one of the earliest studies to identify a dynamic nature in the LSC. Through measuring the horizontal temperature distribution, which is strongly correlated with the LSC, Cioni et al. (1997) demonstrated that the LSC wanders dynamically within a cylindrical domain. Mechanisms contributing to LSC wandering have mainly been attributed to rotations and cessations. The distinction between these two mechanisms is that a rotation is the gradual change in azimuthal angle, and a cessation is when the LSC suddenly dies and then reappears with a totally new orientation (Brown, Nikolaenko & Ahlers 2005;Brown & Ahlers 2006). The rotation process sees a wide range of azimuthal drift speeds, and the slowest speeds have been attributed to a diffusive meandering that is driven by turbulent fluctuations within the flow field (Brown et al. 2005;Brown & Ahlers 2006). Mishra et al. (2011) studied LSC dynamics via DNS and found evidence of partial reversals and double cessations in addition to the standard rotations and cessations, which was also reported by Xi, Zhou & Xia (2006) via experiments. Mishra et al. (2011) also found that cessation is marked by a distinct rise in the amplitude of the second Fourier mode when the LSC dies down. Additional studies have also shown that (a) the LSC experiences a torsional mode that causes its flow near the top and bottom plates to rotate out of phase with one another (Funfschilling & Ahlers 2004;Resagk et al. 2006;Funfschilling, Brown & Ahlers 2008); and (b) a sloshing mode that causes the entire structure to shift back and forth with respect to the LSC symmetry plane (Brown & Ahlers 2009;Xi et al. 2009;Zhou et al. 2009).
Many of the studies cited above have been conducted in unit-Γ cylinders, but studies are also routinely performed in box domains. In boxes the LSC tends to lock into opposing corners, and it will periodically switch between the two pairs of corners in the domain (Bai, Ji & Brown 2016). Stochastic models have been shown to successfully predict the dynamics of the LSC in the Γ = 1 case for boxes and cylinders (Brown & Ahlers 2007, 2008aBai et al. 2016); these models are derived with a primary assumption that the 'wind of turbulence' or the single LSC is present.
While current understanding of the LSC has been fruitful, it is not the only large-scale structure that can reside within turbulent RBC. Du Puits et al. (2007) showed that as Γ increases the wind of turbulence concept breaks down. Numerical studies of RBC in moderate Γ cylinders have shown that the large-scale structure organizes into a series of three-dimensional roll cells (Bailon-Cuba et al. 2010;Sakievich et al. 2016). Qualitatively, these patterns show clear signs of organization that are similar to the structures seen at the onset of thermal convection, but the patterns vary with Γ and Rayleigh number. Furthermore, the dynamic events such as net rotation, cessations and sloshing remain to be analysed and quantified in the larger Γ cylinders. For example, Vogt et al. (2018) recently discovered an appearance of a completely new LSC dynamic mode which they called a 'jump rope vortex', different from both torsional and sloshing modes, in a Γ = 2 cylindrical container. The spatial organization, length scales and time scales of wider Γ cases must be quantified before the predictive powers of the current low-order models can be extended to more general cases of RBC in wider domains.
Recent works (Stevens et al. 2018;Krug et al. 2019) performed direct numerical simulation (DNS) studies at very large Γ in turbulent RBC, and these studies identified persistent structures whose horizontal length scales are several times larger than the height of the domain, which were recently termed as 'superstructures'. While these studies provide a valuable insight into statistical properties of the superstructures, the analysis in these papers was largely concerned with the temporally averaged flow fields, filtering out the temporal dynamics of the structures. Pandey et al. (2018) have devoted significant attention to the determination of the averaging time scales to properly capture the properties of the superstructures, since, evidently, an infinite time average would smooth out all the large-scale structures and result in a horizontally homogeneous mean field. They also examined the large-scale drift of the structure patterns by measuring the phase change of the temporally filtered angular wavenumber in a sliding time average, and identified a slow evolution of the superstructures on a time scale which is of the same order as the filtering time scale that they identified, i.e. tens to hundreds of free-fall time units, depending on the flow parameters, such as Pr and Ra numbers. A similar conclusion was reached in Sakievich et al. (2016) who showed that the time scales required to observe a change in the large-scale pattern in a Γ = 6.3 cylinder with Pr = 6.7 and Ra = 9.6 × 10 7 were of the order of 600 free-fall time units.
Different from the previous work on RBC in large-aspect-ratio domains, which was mainly concerned with the statistical properties of the structures, this work focuses on temporal dynamics and evolution of the different structure modes in a moderately high aspect-ratio RBC domain with Γ = 6.3, as calculated from DNS using a spectral element approach. Unlike the previous works, we focus on an individual structure rather than their statistical representation, so that the dynamical events that are ultimately responsible for the large-scale properties of the turbulent superstructures can be understood from the bottom-up perspective. The analysis of temporal dynamics of the large-scale structure is performed via studying both the temporal evolution and the statistics of the different azimuthal Fourier modes in a cylindrical RBC domain. The Fourier decomposition approach allows us to highlight the dynamical processes at different scales, and their interaction, that accompany the large-scale structure formation and evolution in a cylindrical RBC cell. In this sense, this is the first study which bridges the gap between our understanding of the temporal processes in a large-aspect-ratio domain with a well-studied subject of a temporal dynamics of RBC in a unit aspect-ratio cell, and highlights the similarities and differences between the two cases. The aspect ratio we consider, Γ = 6.3, is about the size of the superstructures naturally found in larger aspect-ratio domains (Hartlep et al. 2003;Stevens et al. 2018). While the effect of the boundary conditions is present in this study, some of the features of the large-scale mode organization and dynamics that we observe resonate remarkably well with both the time scales (Pandey et al. 2018) and the statistical properties (Krug et al. 2019) of the superstructures found in larger aspect-ratio domains.
The current study is performed at a single value for the Rayleigh and Prandtl numbers. It is known that these two parameters have a significant influence on the fluid and heat transfer dynamics in a Rayleigh-Bénard convection, including the effect on both length and time scales (Hartlep et al. 2003;Pandey et al. 2018), and the spatial structure of the flow (Malevsky 1995;Verzicco & Camussi 1999;Breuer et al. 2004). The evidence that exists of the influence of Ra and Pr on length and time scales is unfortunately not enough to form a comprehensive picture at this point. For example, previous studies have shown that the time scales generally increase with Prandtl number while they seem to be relatively unaffected by Rayleigh number at a fixed Pr = 0.7 (Pandey et al. 2018). It was shown that the length scales monotonically grow with Ra at a Pr = 0.7-1 (Pandey et al. 2018;Krug et al. 2019), while this tendency might be reversed at very high values of Pr (Malevsky 1995). The dependence of length scales on Prandtl number at a fixed Rayleigh number seem to exhibit a growth with Pr followed by a decay, with the maximum value around Pr ∼ 7-10 (Hartlep et al. 2003;Pandey et al. 2018). In terms of structure, the large-scale flow organization is dominated by rolls at lower Pr, and by cellular structures at higher Pr, with the transition between the two regimes happening around Pr = 7 (Malevsky 1995;Verzicco & Camussi 1999;Hartlep et al. 2003;Breuer et al. 2004). The current study can thus be viewed as an in-depth exploration of one particular regime, while larger parametric studies are required to achieve a broad understanding of the temporal dynamics of large-scale structures across the parameter range.
The paper is organizes as follows. In § 2 we present the problem formulation and the numerical method. In § 3 we comment on the mean field of RBC achievable with the finite-time DNS simulations. In § 4 the large-scale structure organization is described via a temporal filtering approach. In § 5 we analyse the temporal dynamics and the integral times scales of the azimuthal Fourier modes. In § 6 statistical properties of the Fourier modes and their spatial variability are discussed, while in § 7 concluding remarks are presented.

Problem formulation and numerical method
The purpose of this section is to provide a description of the problem formulation, notation, governing equations and numerical methodology used throughout this study. This work relies heavily on Fourier decomposition to analyse the structure of the flow field, and so a small primer on Fourier decomposition is included at the end of this section.

Equations, computational domain and scaling
The computation domain Ω in this study is a cylinder with height H and diameter D. We can express Ω in cylindrical coordinates that are normalized by H and symmetrized about the mid-plane (H/2 → z = 0) such that the normalized Ω is defined as where Γ is the aspect ratio (Γ = D/H) of the cylinder. Here Ω is also aligned with the gravitational vector (g) such that g/|g| = −ê z , whereê z is the unit normal in the z-direction. Velocity and temporal units are normalized by the 'free-fall' velocity (w f = √ βgΔTH) and time (t f = H/w f ), where β is the coefficient of thermal expansion, g is the gravitational constant and ΔT is the temperature difference between the top and bottom plates of the convection cell. Here ΔT is used to normalize the temperature field and the Boussinesq approximation is applied to the incompressible Navier-Stokes equations for computation and analysis. The reference temperature for the Boussinesq approximation is taken to be the average mid-plane temperature such that Utilizing these scales, the non-dimensional form of the Boussinesq equations for RBC can be expressed as where u, p and ϑ are the dimensionless velocity, pressure and temperature. The Rayleigh (Ra) and Prandtl (Pr) numbers in (2.4) and (2.5) are defined as where α and ν are the thermal diffusivity and kinematic viscosity, respectively. In this study Γ = 6.3, Ra = 9.6 × 10 7 , Pr = 6.7. When presenting the results, the velocity field is expressed in terms of cylindrical coordinates such that u = {u r , u θ , u z } and u i represents any of the components {u r , u θ , u z }.
Two primary time scales are used in this work. The free-fall time scale (t f ), and an eddy turnover (t ). Here t approximates the mean up and down times for large-scale motions that span the depth of the domain. It is defined as where V indicates a spatial average over the entire domain volume V and t is an average in time. Eddy-turnover time can be thought of as the average time it will take for a particle to cross the layer depth twice driven by a turbulent diffusion. For reference, t is approximately 31t f in this study.

Numerical method
The data in this study is obtained by DNS using the open source spectral element code Nek5000. Nek5000 is a thoroughly validated research code that has been used extensively in scientific literature (see Fischer, Lottes & Kerkemeier 2008). The spectral element method uses high-order polynomial approximation of the solution within each element, while local element solutions are assembled globally through gather-scatter operations (Deville, Fischer & Mund 2002). A P n − P n−2 (Fischer 1997) spectral element formulation is employed herein. Temporal integration is of second order, utilizing an implicit backward-differentiation formula for the viscous terms, and an explicit extrapolation for nonlinear and Boussinesque terms. Pressure decoupling is performed using the operator splitting method (Karniadakis, Israeli & Orszag 1991), and a resulting pressure Poisson equation is solved by GMRES (Saad & Schultz 1986) with a multigrid preconditioning. The computational grid is discretized with hexahedral elements and a marginal amount of biasing toward the upper and lower plates is applied to the element distribution. The spectral element method (SEM) used in this simulation also applies a Gauss-Lobatto-Legendre (GLL) quadrature within each element which clusters points toward the boundaries of each element and greatly improves resolution at the walls. z y x z y x FIGURE 1. Images of the computational grid using spectral elements (44.6 million grid points with a ninth-order polynomial basis in the spectral elements) (a) and an example of the cylindrical coordinates post-processing grid with a resolution of [60,512,32] points in the r, θ and z directions, respectively (b). The sampling resolution in (b) is reduced from the actual resolution used for analysis [160,2048,64] to improve image quality in the figure.
Ninth-order polynomials were used for the quadrature resulting in roughly 44.6 million grid points. A schematic of the computational domain and the spectral element grid employed is provided in figure 1(a). In the experimental work of Fernandes (2001) it is reported that the Kolmogorov length for RBC at this Γ , Ra and Pr is approximately 1.2 × 10 −2 H, and our simulation's grid had five points within this range at the wall. We also determined that this grid satisfies the spatial resolution criteria of Grötzbach (1983). The temporal resolution for each time step was approximately t × 10 −4 with a corresponding Courant-Friedrichs-Lewy range of ∼0.6-0.7. In our prior work we conducted a a-posteriori analysis to evaluate the resolution of our results utilizing the techniques outlined by Scheel, Emran & Schumacher (2013) and confirmed that our simulation met the requirements for dissipation continuity across elements, as well as the resolution with respect to the height dependent Kolmogorov and Batchelor scales (see Sakievich et al. 2016). Boundary conditions are specified as follows: no-slip conditions for velocity are set at all bottom, top and side walls. Temperature is set to a constant value T = T hot at the bottom wall, and T = T cold at the top wall, with the temperature difference ΔT = T hot − T cold used for non-dimensionalization in (2.2) and (2.6), while adiabatic boundary conditions are used for the temperature at the side walls. Additional details regarding resolution, convergence and comparison with experiments for the specific computations in this work can be found in the prior work of Sakievich et al. (2016). The total run time of the simulations in the current work is 3054 free-fall time units, or close to 100 eddy-turnover times, which makes it one of the longest DNS studies of turbulent Rayleigh-Bénard convection up to date (Sakievich et al. 2016;Pandey et al. 2018), and perhaps the longest for the considered Rayleigh number. Note that while the current study undoubtedly pushes the limit in terms of the simulation time, it still only covers less that one viscous time unit t v = √ Ra/Pr t f = 3800t f , and only 10 % of a thermal diffusion time unit t d = √ Ra Pr t f = 25 460t f for the current Pr number, once again highlighting the challenges of accessing the longest possible time scales in the RBC studies at high Ra number. While the total span of the simulations is more than enough to capture the superstructures (Pandey et al. 2018) and their dominant dynamics, it might not be enough to capture the effect of slow diffusive processes on their evolution.
The data for calculating statistical quantities was sampled every three free-fall times. This choice of a sampling rate is a compromise between fine resolution in time and a storage requirement for a very long temporal study, such as the one performed in the current work. The sampling rate is more than adequate for studying the long time scales and slow energetic processes associated with the evolution of the large-scale structures, which is the focus of the current work. For post-processing, each snapshot is projected onto cylindrical coordinates using spectral interpolation routines native to Nek5000, and the velocity components are transformed from Cartesian to cylindrical. This was previously done on a smaller scale (Sakievich et al. 2016), but in this work the transformation has been extended to the entire domain. Cylindrical coordinates are the logical choice for analysing the current dataset and facilitate operations along the domain's periodic, azimuthal direction. The DNS snapshots are resampled with [160,2048,64] points in r, θ and z, respectively, to generate the cylindrical grids used for analysis. Non-uniform, Gauss-Legendre (GL) quadrature is used to sample in the r and z directions, but the θ direction uses equispaced sampling points to facilitate Fourier transforms. Gauss-Legendre quadrature does not include the endpoints and is defined on the standard interval x ∈ (−1, 1). Gauss-Legendre quadrature is selected to facilitate high accuracy numerical integration and to remove unnecessary sampling at the walls of the cell where the system is well defined. Boundaries in the z direction are constrained with Dirichlet boundary conditions, so that sampling on them for Fourier transforms is trivial. Points along the central axis (r = 0) are at a spatial singularity in the cylindrical coordinate representation and provide no additional data when Fourier transforms in θ and integration over the r-z plane are applied. The points along r = Γ /2 have Neumann boundary conditions in the temperature field but virtually no information is lost since the gradient at the wall is zero (adiabatic), and the GL quadrature samples very close to the boundaries. An example image of a post-processing grid in cylindrical coordinates is shown in figure 1(b).
The post-processing grid has a little less than half the number of grid points when compared with the computational grid, but the change in coordinate system and quadratures leads to non-uniform sampling ratios with respect to the original domain Ω. In the vertical direction this causes the post-processing grid to have approximately twice as many points in the boundary layer region, and about half as many in the bulk region. The horizontal sampling is not as straightforward to compare because the coordinate systems are fundamentally different. However, it is definitely the case that the centre of the post-processing grid is sampled more finely than the computational grid. The spectral interpolation used during resampling ensured that no wavenumber information was lost on a post-processing grid according to its resolution, but neither is gained, compared to the original computational grid. Figure 1 provides an example of the two different grids (computational and post-processing), however, the number of depicted grid points is reduced in figure 1(b) with respect to the actual post-processing grid used for analysis to make the visualization comprehensible. The spectral element grid in figure 1(a) does not require a reduction in sampling for visualization purposes because the GLL points within each element can be given a different contrast.

Fourier decomposition
Fourier decomposition in the azimuthal direction renders insights into the structure of the flow field. Fourier modes are a natural choice because the azimuthal direction is geometrically periodic. Fourier decomposition provides additional benefits in this study that extend beyond the mathematical significance of the modes. For example, azimuthal motions for RBC in cylinders tend to evolve on extremely long time scales, and the azimuthal velocity signals are relatively weak (Brown et al. 2005;Mishra et al. 2011).
Performing an analytical decomposition such as Fourier analysis allows the azimuthal evolution of the flow to be studied in a well understood format with precise measurements.
For a general signal u (r, θ, z, t), azimuthally periodic with a period 2π, the Fourier series decomposition is given by where j = √ −1 and k is the integer Fourier mode number. Fourier coefficients are given by the inverse operationû (2.10) The Fourier series representation can be approximated by a finite number of modes, N, and a uniform convergence of a truncated series to the original signal u(x) is guaranteed given the signal is (a) smooth, and (b) periodic. A discrete Fourier transform in this case can be defined asû (2.11) Note that in this formulation, k can be understood as the integer azimuthal mode number, and, at each particular radius r, the angular wavenumberk could be defined as so that the wavelength, that gives a relation to a physical length of the structures, for each particular azimuthal mode k at each particular radius is given as where the dependence on r denotes that the wavelength for a certain k mode has been calculated at a given radius. Throughout this work, Fourier coefficients are indicated by theû accent, and the Fourier operator (computed with its discrete representation (2.11)) is indicated by F [u]. For each complex Fourier coefficientû k , its amplitude (2.14) and phase can be defined, where Re(û k ) and Im(û k ) correspond to the real and imaginary parts of u k , respectively. All averages are noted by the brackets and subscripts are listed by the order in which the averaging operations are applied. For instance, u z (r, θ, z, t) θ,t (r, z) is the vertical profile of the vertical velocity field after averaging in the azimuthal direction and in time.
Throughout the paper, we also look at the scaled volume integrated values of the Fourier coefficients defined as Note that the integral defined in (2.16) is also equivalent to (2.17)

The mean field
The primary interest of this study is to investigate the spatial and temporal properties of the large-scale structures in the flow field. These structures all have a finite life span and, therefore, reside in the fluctuating field with respect to Reynolds decomposition. However, fluctuations must be defined with respect to mean values. Therefore, it is essential to define the mean field and the averaging operators that create the mean field about which the fluctuations occur. Turbulence analysis typically employs the assumption of ergodicity to define the steady-state mean field as a time average. In terms of dynamical systems ergodicity means that every point in state space is sampled during the system's evolution, though in practice it is often sufficient to sample a large number of statistically independent instances. This has been a challenge in large Γ RBC studies (Bailon-Cuba et al. 2010;Emran & Schumacher 2015;Sakievich et al. 2016) because the large-scale patterns evolve over very long periods of time, rendering the available samples statistically correlated, at least with respect to large scales. Adrian, Sakievich & Peet (2017) define these large-scale organization as 'super-coherent states' because the strong spatio-temporal coherence persists over many eddy-turnover times. A super-coherent state can be thought of as a deep basin of attraction in state space where the realizations within the basin have a strong similarity (i.e. they are highly correlated) over a very long period of time. It can seem like the system is converging to a steady state within one of these basins when in reality the total state space can contain other deep basins, and transitions between these basins may be triggered by some perturbation events. Adrian et al. (2017) show that in the case of moderate and large Γ RBC in cylindrical domains there are multiple states that have the potential for super-coherency and that some of these states can be identified by the symmetries of the domain and the boundary conditions. For example, in our previous work (Sakievich et al. 2016) a large, long-lived updraft was observed in the central region of the cylinder, and this updraft biased the statistics in an otherwise stationary system. This updraft could just as likely have been a downdraft in another realization of the flow, just as an evenly balanced coin has equal probability of landing on heads or tails over a sufficient number of samples. Another example of a deep basin is highlighted by the azimuthal symmetry of a cylindrical RBC cell. Since there is nothing to constrain the orientation of the flow's structure the system has an equal probability of assuming any azimuthal orientation, and so ergodicity demands that over a long enough period every orientation must be sampled. Even if the large-scale structures are not observed to rotate about the cylindrical domain's central axis during a simulation or experiment within the state space that defines the Reynolds average, these different orientations must be sampled. The key mechanisms that have been identified for reorienting the large-scale structures are cessations and rotations, but up to this point they have only been observed in Γ ≤ 1 systems (Brown et al. 2005;Mishra et al. 2011).
The challenge for turbulent RBC is that obtaining a sufficient number of samples, or sufficient averaging time to sample many of these deep basins of attraction is a non-trivial task since these dynamic events tend to evolve over long periods of time and state perturbations occur infrequently.
Luckily the symmetries of the RBC system can be employed to account for orientation changes and directional bias of the large-scale structures. Conventionally employed azimuthal averaging can be used to account for all possible azimuthal orientations of an observed structure, while Adrian et al. (2017) define an additional transformation to account for the vertical antisymmetry imposed by the thermal boundary conditions and the direction of gravity. This transformation can provide a mean field that accounts for the equal probability of updraft and downdraft. An unbiased mean is defined as a mean computed based on the sum of the current state u (r, z, θ, t), and a new state S [u(r, z, θ, t)], where S[ ] is the symmetry transformation operator defined by Adrian et al. (2017), and is described in more details in appendix A. Due to a commutativity of the azimuthal and time averaging operators with the symmetry transformation, a new unbiased mean can be defined as where the superscript 'SM' stands for a symmetrized mean. Fluctuating quantities are then defined with respect to this new mean field, The mean operator SM θ,t can also be expressed in terms of Fourier coefficients, as given below since the contribution of any non-zero azimuthal mode will be zero due to azimuthal averaging. Conversely, this also means that the fluctuating field contains the higher-order modes k > 0, and the mode k = 0 fluctuations. In this work we have applied the transformation defined by (A 1) to the azimuthally and time-averaged mean fields to construct an approximation based on (3.1) for the Reynolds averaged mean fields from the available finite-time DNS sampled data originally containing 'super-coherent states' (Sakievich et al. 2016;Adrian et al. 2017). Figure 2 shows the symmetrized mean statistics for the temperature and the azimuthal velocity field, with radial and vertical mean velocities superimposed as vector plots in figure 2(b). The mean field in figure 2 displays several interesting characteristics. We note that the transformation (A 1), (3.1) is designed to remove the bias in preferential statistically significant thermal updraft and downdraft congregations, which primarily targets a symmetrization of the temperature field and the associated vertical velocity field. Judging by a relative quiescence in both temperature and vertical velocity non-uniformities in the central region of the cell in figure 2, one can see that the transformation did successfully achieve that goal. However, starting at the sidewalls (r = Γ /2), two counter-rotating roll cells can be observed with stagnation point at z = 0 (horizontal midplane) where the two roll cells meet. Additionally, a thermal boundary layer can be seen along the adiabatic sidewalls. These roll cells and the accompanying boundary layers are in-line with the vertical antisymmetry of the mean temperature field and are not expected to be removed FIGURE 2. Azimuthal and temporally averaged mean fields. The colour plot in (a) corresponds to ϑ θ,t and in (b) it corresponds to u θ θ,t , while the vector field in plot (b) is the two-dimensional vector of {u r , u z } θ,t where the peak vector magnitude is 0.06 w f . All length scales are normalized with H.
by a proposed transformation , and, thus, they are likely to be present in the true Reynolds averaged flow field due to the effect of side walls. The mean azimuthal velocity component shows that a preferential direction for rotation or drift is not consistently present across the entire time series, i.e. net rotation is not a product of the mean that has been constructed from this dataset. However, local patches of weak non-zero azimuthal velocities can be observed in the mean field presented in figure 2(b). As noted in the appendix A, the current transformation is not meant to mitigate a potential bias in the azimuthal velocities due to a structure drift. Accounting for the azimuthal symmetry with yet another transformation could further improve the mean statistics in the azimuthal velocity field shown in figure 2(b). On the other hand, the computed mean azimuthal motions are weak and will likely asymptotically approach zero 'in a natural way' as the number of samples is progressively increased. Finally, we would like to note that the observed corner structures in temperature and velocity fields are solely due to the effect of side walls, no-slip or stress free, and, thus, are expected to vanish in the situations with periodic boundary conditions and with infinite aspect ratio Γ → ∞ cells. A mean flow that is identically zero everywhere would thus be expected for the infinite aspect ratios. However, one has to be careful with periodic boundary conditions, since they still impose a non-physical domain truncation and, thus, unphysically affect the length scales of the structures that have to be 'squeezed' into a finite-size domain. It is hard to say what effect, if any, it will have on mean flow in the RBC problem, but unphysical modifications to mean flow due to periodic boundary conditions were previously reported in the simulations of channel flows due to a 'structure locking' phenomenon (

Global description of the large-scale structure
In this section the largest scales of the flow field are investigated using azimuthal Fourier decomposition. They are of interest because of the important role they play in transporting energy. It will be shown that these large-scale structures contain the majority of the energy in the flow field, persist for long periods of time, and are responsible for much of the flow inhomogeneity. Figure 3 shows the scaled volume integrated (see (2.16) for the definition) and time-averaged energy spectra of the three velocity components, {|û (r, k, z, t)| 2 } V/2π t , and temperature (defined analogously).
The spectra in figure 3 indicate that the k = 2 Fourier mode is the most dominant mode over the range of the simulation. The peak is very pronounced in the temperature and azimuthal velocity fields, but more subtle in the radial and vertical velocity components. Even though the spectra in figure 3 indicates that the dominant structure over the life span of the simulation was the k = 2 mode, it does not provide a clear indication of the spectra's evolution in time. In the authors' previous work (Sakievich et al. 2016) and the work of Bailon-Cuba et al. (2010) no significant evolution of the large-scale structures was observed on the time scale of the simulations which was about O(10 2 )t f at similar aspect ratios. However, in the work of Emran & Schumacher (2015) it was estimated that the large-scale structures would drift on time scales of O(10 3 )t f , albeit for a higher aspect ratio, lower Pr and a lower Ra flow field. In the present work the simulation time has been extended to the order where a global drift of the large-scale structures has been predicted by Emran & Schumacher (2015). Note that for the current parameter range of Ra = 9.6 × 10 7 , Pr = 6.7, the drift time scales might be even higher (Brown et al. 2005;Pandey et al. 2018).
Observations of the evolution in the flow field are provided in this section through temporal filtering of the dataset. Temporal filtering is defined as a box cut filter (1/T) where the time period T = 600t f was chosen for the filtering duration, and the starting times of filtering t i were, respectively, 0, 600t f , 1200t f , 1800t f and 2400t f . Figure 4 contains the scaled volume integrated energy spectrum of the temporally filtered temperature field for these time intervals, {|F [ ϑ (r, θ, z, t) t ]| 2 } V/2π . Temporal filtering removes the majority of the small-scale structures leaving the highly correlated large-scale structures clearly visible, and it is a good technique for observing the slowly evolving large-scale dynamics (Sakievich et al. 2016). The temperature field's spectrum is selected for comparison because it contains the most distinguished peak in figure 3. Visualizations of the temporally filtered temperature field as it evolves in time are provided in figure 5, and the instances in figure 5 correspond to the energy spectra plots in figure 4. Figure 4 shows that over the first 600t f the dominant mode is k = 3, which corresponds to the structure observed in Sakievich et al. (2016). However, in the next 600t f , the dominant mode transitions to k = 2. The first instance of the filtered field also shows a larger distribution of energy in the other low-order modes, including the zeroth and the first mode (see figure 4), but by k = 12 the energy content is about the same for all instances of the filtered field. The second instance of the filtered field shows higher energy content in modes k = 1 and 3 (which are likely the remnants of the previous state of the structure), but by the third instance the energy has concentrated predominantly in k = 2. One possible interpretation of this transition is that the observed k = 3 dominant structure is less stable than the structure corresponding to k = 2 because the turbulent thermal energy is distributed among a larger number of low-order modes.
Looking at the individual modes can help explain their contribution to the overall flow field. The first few low-order modes and their cumulative summations corresponding to temperature fields in figures 5(a) and 5(e) are provided in figures 6 and 7.
. Scaled volume integrated energy spectrum for the temporally filtered temperature field. Filtering is performed by applying a temporal average with a period of 600t f . The legend entries refer to the averaging period of each instance in multiples of t f .
The modes in figure 6 can be interpreted with the following roles: k = 0 establishes a central, warm column, k = 1 and 2 shift the central column and bias the structure away from the centre breaking its symmetry, and k = 3 finalizes the hub-and-spoke like pattern that was outlined in Sakievich et al. (2016). A qualitative comparison of figures 6(h) and 5(a) show that the total structure is well described by the first 4 (k = 0 : 3) Fourier modes and figures 6(i) and 6( j) show that three-dimensional representation of these modes takes the form of large-roll cells that span the entire layer depth.
However, examination of the modes displayed in figure 7 shows that as the simulation evolves the structure becomes almost fully described by k = 2. This convergence of energy and structure toward a single mode seems to indicate a stabilization for the system as a whole. It could be the case that k = 2 is the long-term structure of the flow field at this Γ , Ra and Pr, but because the k = 3 structure remained coherent for approximately 1/5th of the total simulation time, nothing definitive can be determined. It is still possible that the system could undergo yet another transition or modulate back to a k = 3 dominated structure. For future studies of RBC, it is also worth noting the length of time in which this transient evolved in moderate to large Γ where multi-roll cell structures persist, especially since the time scale of this transition is larger than the entire simulation time in many previous studies. We should also note that even though we are using the dominant modes to describe the overall structure observed in this study, it is entirely possible for the modes to decay and emerge independent of one another. For example, it is possible that the k = 0 and/or k = 1 modes will reappear while the k = 2 mode is dominant and create a different state which could include global updrafts and downdrafts observed with the current k = 3 structure, or cause a reversion back to the k = 3 state. In fact, a slight downward reversal of the k = 0 mode can already be noticed in figure 5(e), where a weak connection of cold plumes, as opposed to hot plumes, can be observed, and in figure 7(e), where slight negative temperatures in the central region resulting form k = 0 mode contribution can be seen.
To draw a similarity with a turbulent pipe flow, a concentration of energy in a few low-order azimuthal modes was observed experimentally by Bailey & Smits (2010), and computationally by Baltzer, Adrian & Wu (2013), with the dominant azimuthal mode being k = 3 in both studies.
and relating this wavelength Λ k to a typical size of the superstructures reported to be around 6-7H in terms of their spectral wavelength for convection in air (Hartlep et al. 2003;Pandey et al. 2018;Stevens et al. 2018), yields k = 2.8-3.3 for the Γ = 6.3 case. However, we should bear in mind that the current simulations are performed for convection in water, with Pr = 6.7, versus convection in air, where Pr ∼ 0.7. Pandey et al. (2018) showed that the structure sizes are supposed to increase for water versus air by approximately 1.5 times, while Busse (1994) measured the sizes to be around 9-10H for his experiments with Pr ∼ 7 fluids at Ra = 10 5 -10 6 . For Λ k = 10, (4.1) recovers k = 2 precisely, thus hinting that both k = 2 and 3 modes are consistent with the previously reported structure sizes. Moreover, for Γ = 1, this relationships gives Λ 1 ∼ 3 being the largest possible wavelength fitting into a Γ = 1 cell, which is still smaller than the natural size of the superstructures that would want to be formed, explaining why the k = 1 mode is dominant in the Γ = 1 case. Note that the wavelength Λ k defined in (4.1) is essentially equivalent to λ k (r) in (2.13) evaluated at r = Γ /2 for a given RBC cell.

Temporal evolution of the flow field
In the previous section it was shown that the large-scale flow transitioned from a structure dominated by a k = 3 Fourier mode to one dominated by a k = 2 Fourier mode. In this section the temporal evolution of the transition will be investigated in greater detail using visualization of the unfiltered mid-plane temperature field. Figure 8 presents snapshots of ϑ(r, θ, z = 0, t) at times that bracket the transition each separated by 90t f . In figure 8(a) the flow is dominated by the k = 3 mode, albeit with the warm spokes located at two o'clock being rather weaker than the other two spokes. In figure 8(b) the one o'clock spoke weakens more as the upward plume at the centre of the cylinder becomes stronger. The process continues in figure 8(c) wherein the warm plume at one o'clock has almost vanished, and the central updraft plume has moved toward the x-axis x-axis -5.0 × 10 -2 5.0 × 10 -2 -2.5 × 10 -2 2.5 × 10 -2 0 A more detailed investigation of the transition is performed by plotting the scaled volume integrated Fourier coefficients for a given mode defined by (2.16).
Volume integration removes the localized spatial variations of the mode and allows the temporal evolution to be investigated from a macro perspective. Even though the volume integrated Fourier coefficients only depend on wavenumber and time, they are still complex variables. The phase and amplitude of the volume integrated coefficient can simultaneously change with time.
Figures 9 and 10 display the temporal evolution ofû r andθ fluctuations for the first five (k = 0 : 4) volume integrated Fourier modes in terms of their phase and amplitude. Note that the k = 0 modes represent an azimuthal mean of the corresponding physical quantities, and their further integration with (2.17) results in volume averaging in physical space, which is identically zero forû z andû r due to mass conservation. The reader is cautioned that while interpreting figures 9 and 10 to remember that any phase jumps of 2π in the plots are continuous in phase space. In general, the modes in figures 9 and 10 share some common behaviour. The variables and plots are divided into two highly correlated groups,û z withθ in figure 9 andû r withû θ in figure 10. A strong degree of correlation betweenθ andû z dynamics testifies that temperature and vertical velocity are the signatures of the same large-scale structure, consistent with the findings of Krug et al. (2019). Theθ andû z coefficients for the low-order modes except for k = 2 show larger mean amplitudes and smaller fluctuations in phase during the initial part of the simulation. During the transition from a k = 3 to k = 2 in the interval 500 < t/t f < 1000, the amplitude of the dominant structures tends to decrease, and the phase fluctuations tend to increase.
While the signature of the transition from a k = 3 to k = 2 dominant structure around 500-1000t f can be seen in each of the low wavenumber modes, the effects are most clearly displayed in the k = 2 and k = 3 modes. Inspection of the amplitude and phase for k = 2 and k = 3 shows that the amplitude of k = 3 starts out at a relatively large value and, beginning at around 500t f , decays steadily to a negligible value over another 500t f , completing the process around 1000t f . Furthermore, while the amplitude is large, the phase remains relatively constant, but once the amplitude dies down the phase fluctuations increase. The opposite effects are seen in the k = 2 mode where the amplitude is initially low but then gradually increases over the same period that k = 3 decays.
A simple analysis of the time scales involved with the dynamical processes occurring in Rayleigh-Bénard convection can be presented by estimating a 'mode turnover time', which is the time it takes a particle to travel one full circuit in a particular mode k. The full circulation length for each mode (normalized with H) can be defined as where Λ k is the kth mode wavelength given by (4.1), which yields The 'mode turnover' time therefore can be defined through a classical eddy-turnover time as Note that in the eddy-turnover time definition, (2.8), the variance of the vertical velocity fluctuations is used as a velocity scale. In the present study the magnitude of the horizontal and vertical velocity fluctuations is of the same order, consistent with the observations in which 'accounts for the fact that an individual parcel is not perfectly circulating around in such a roll when the flow is turbulent'. We choose not to introduce any empirical factors, but admit that the proposed time scale is only an order of magnitude estimate, and it might take several of such time scales for a particular event to happen, as discussed below.
Since the mode turnover time reflects the time it takes an information (disturbance) to propagate across the entire mode, it should be representative of the time scales associated with the mode destabilization. The time scale analysis with (5.3) predicts, for an aspect ratio Γ = 1, where a single k = 1 mode dominates, a time scale of t (1) 1 = 2.57t .
For an aspect ratio Γ = 1, Mishra et al. (2011) identified a time scale of 10t as the time scale associated with the LSC reversal, while Brown et al. (2005) reported 30t as the average time between reorientations. However, these time scales refer to the time between the reorientation events, and not to the duration of the events themselves. If one closely examines the data presented in, for example, Brown et al. (2005), Mishra et al. (2011) and Zürner et al. (2019) one can see that the time scale associated with the process of transition of the global structure into a state with a new LSC reorientation, i.e. from the beginning of the mode destabilization to the time when it stabilizes again, is on the order of 2t in Brown et al. (2005) and Zürner et al. (2019), and on the order of 2.63t in Mishra et al. (2011), well in line with the predicted t (1) 1 = 2.57t from (5.3). Thus, while cessations in unit aspect-ratio cells are often perceived as instantaneous events, they are, in fact, events of a finite, albeit short, duration. While a phase reversal happens almost instantaneously during the LSC reorientation, the other accompanied events, such as a decrease of the k = 1 mode amplitude preceding the phase reversal (or partial reversal), and a subsequent increase back to its original value, happen more gradually (Brown et al. 2005;Brown & Ahlers 2007). For a current aspect ratio of Γ = 6.3, the time scales associated with the destabilization of the modes are about six times larger for a given k. The corresponding time scales for the first three modes k = 1, 2, 3 would be t (6.3) 1 = 10.90t , t (6.3) 2 = 6.00t and t (6.3) 3 = 4.30t , which, if scaled with the free-fall time units, are t (6.3) 3 ∼ 130t f . Note that the duration of the observed 3 to 2 mode transition correlates well with these time scales.
One can observe that as the mode number k → ∞, the mode turnover time converges to a constant value of t , irrespective of Γ . This perhaps can explain a self-similarity of higher-order modes, and an apparent lack of dependence on the problem geometry, as will be illustrated below. This also shows that the time scale of an eddy turnover is still a valid measure to describe the fast processes in the RBC cell.
The 3 to 2 mode transition in the global structure in the current Γ = 6.3 domain occurs on a time scale of approximately 500t f , or 16 eddy turnovers, but it involves not only modes k = 2 and k = 3, but the other low-order modes, for example, k = 1. A careful observation of the k = 1 mode behaviour in figure 9(c) shows a rather rapid decrease in its amplitude around the time t = 500t f followed by a relatively rapid increase back to its original value, with another event of a vanishing amplitude for k = 1 at t = 1000t f . A rapid decay in amplitude remarkably resembles the cessations observed in Mishra et al. (2011) where it was shown that the k = 1 mode amplitude rapidly vanishes during cessation. A close inspection of the phase diagram in figure 9(d) reveals that this rapid decrease is indeed accompanied by a phase shift close to 180 • , pointing to a reorientation in a mode 1. These mode 1 cessations happen to fall onto the interval corresponding to the mode k = 3 to 2 transition, which suggests that the two events might be associated with each other. Interestingly, the duration between the consecutive cessations correlates well with t (6.3) 1 predicted by (5.3), which might signify that they are associated with the same k = 1 destabilization event and might resemble double cessations identified in Γ = 1 cells (Xi et al. 2006;Mishra et al. 2011), albeit on longer time scales here than in unit aspect-ratio domains, commensurate with the 6.3 times difference in the aspect ratio. Note that mode 2 also exhibits an event similar to a cessation at a time of approximately 300t f where its amplitude essentially vanishes and its orientation rapidly changes. While such cessations in a dominant mode would lead to a complete reorientation of the structure, they do not result in observable changes when the modes are not dominant. This once again testifies of a relative complexity of a wide aspect-ratio system, where a global reorganization is a more complex process that involves the mode interaction.
FIGURE 11. Temporal evolution of the scaled volume integrated Fourier coefficients plotted on the complex plane for k = 2 (a,c) and k = 3 (b,d), where (a,b), • =û r , and (c,d), =û θ . Markers are colour-coded according to their time, from blue (start of the simulations) to red (finish).
In general, the k = 2 and k = 3 modes exhibit similar characteristics during the time periods when each respective mode is the most dominant mode in the flow field. However, there is one notable difference between the temporal evolution of the phases for these two modes. While the phase of k = 3 remains virtually constant during its period of dominance (albeit showing signs of low-amplitude, low-frequency fluctuations), the k = 2 mode's phase changes at a relatively constant rate. This indicates a steady rotation event for the k = 2 mode that becomes strikingly clear when the scaled volume integrated coefficients are plotted on the complex plane. Plotting on the complex plane (see figures 11 and 12) is another intuitive way to view the changes in amplitude and phase for a given wavenumber.
The trajectory of the k = 2 mode for the temperatureθ in figure 12(c) begins near the origin and as time progresses it tracks up along the imaginary axis and then begins to drift into quadrant two of the real-complex plane. The rate of rotation manifested by a steady increase in the phase angle Φ is measured to be 1.1 degrees per eddy turnover by As a general observation, when the mode is the dominant mode in the large-scale structure, its complex Fourier coefficient drifts away from the origin, and when it loses its dominance it becomes centred around the origin with a rapidly changing amplitude scattered between zero and some threshold value. The same can be said about the behaviour ofû r andû θ coefficients, even in the modes 2 and 3, shown in figure 11, which reflects the fact that these variables do not play an active role in the global structure dynamics. The key observation here is that active components of the dominant modes display persistence in terms of phase and amplitude of their Fourier coefficients, while the supporting modes display chaotic behaviour with zero mean in all their Fourier coefficients.
Even though there appears to be no net rotation when the k = 3 mode is dominant on a time scale while it persisted, a rotation event is occurring in the k = 1 mode. Figure 9(d) shows a steady change in phase for the first 1500t f of the simulation, after which it suddenly begins to see large fluctuations in phase like the k = 3 mode.
During the first 500t f , the k = 1 mode rotates approximately 60 • with a least squares fit providing a rotation speed of approximately 3 degrees per eddy turnover. The phase jumps up and down during the period 500-1000 t f between the two cessations (marking the period of mode 3 to 2 transition), when at the time of 1000t f , a sudden acceleration occurs, and a higher value of rotation velocity persists until 1500t f , at which time the phase randomizes. Interestingly, the moment of sudden acceleration in k = 1 rotation coincides with a second cessation of this mode, and a completion of k = 3 to k = 2 transition. Brown & Ahlers (2007) reported a similar phenomenon of a rapid acceleration of motion in the azimuthal direction accompanying a cessation predicted by their stochastic LSC model which was explained by the fact that when the amplitude becomes small, the azimuthal motion becomes fast, due to a reduction in the LSC angular momentum. The above observations indicate that there may be a connection between the dynamics of the k = 1 mode, i.e. rotation and cessation, and the transition of the flow's global structure. The k = 1 mode rotation and cessation supports our observation that the transition is marked by the movement of the central column, while the phase fluctuations of the k = 3 mode can explain the azimuthal shifts of the large-scale structure lobes along the edge of the domain (see figure 8).
The rotation rate of 3 degrees per eddy turnover corresponds to a time scale of rotation of approximately 40-60 eddy turnovers (per 1/2 revolution). A rotation rate of 1.1 degree per eddy turnover observed for mode 2 corresponds to the time scale of 160 eddy turnovers (per 1/2 revolution). Interestingly, very different time scales associated with the rotation were also observed in Brown et al. (2005). Fast rotation typically classified as a reorientation seemed to scale with the time of about 10 eddy turnovers for Γ = 1, which is approximately 3.6t (1) 1 , while much slower drift that was not associated with a reorientation event would occur on time scales of an order of magnitude larger. Similar time scale discrepancy in a duration between different reorientation events was observed with Γ = 1 in Sreenivasan, Bershadskii & Niemela (2002). The current Γ = 6.3 cell is shown to exhibit similar dynamics, where the rotation of mode 1 accompanied by a mode transition scales on a mode turnover time 3.6-5.5t (6.3) 1 , while a slow rotation in a persistent mode 2, not undergoing a transition, is approximately 26.6t (6.3) 2 . Brown & Ahlers (2006 attributed slow rotations, or drifts, to diffusive processes. Viscous time scale is given by (Pandey et al. 2018), which, for given flow parameters, is equal to 3800t f , or 120t , while a thermal diffusion time scale is t d = √ Ra Pr t f = 25 460t f = 804t for the current Pr = 6.7 case. It is seen that the slow mode rotation events ('azimuthal jitter') are expected to occur on very long time scales in the current case, and a slow rotation observed in mode 2 indeed falls in between these time scales.
So far virtually all of the discussion for the volume integrated Fourier coefficients has centred on theθ andû z fields. These two fields are highly correlated, and they tend to describe the events that are oriented in the inhomogeneous, vertical direction, and provide a profound amount of information regarding the time scales and dynamics that are associated with the large-scale structure's transition. Figure 10 shows that the amplitude ofû r changes rapidly in time, but the magnitude of these fluctuations stays within a specific band for each wavenumber. Meanwhile, the phase ofû r for each wavenumber tends to evolve at a slower pace with several low frequency components. These time scales are measured by performing a fast Fourier transform (FFT) analysis of theû r andû θ signals in figure 10 (using time as the abscissa). This analysis allows us to identify the temporal frequencies and their accompanying periods that are associated with the consistent phase oscillations in these fields. Since these temporal signals are not strictly periodic, a Blackman window function is multiplied with the signal prior to taking the FFT to improve peak detection by limiting spectral leakage (Blackman & Tukey 1958). The measured peaks are not consistent between wavenumbers, but the general trend is that theû r andû θ amplitudes fluctuate with time scales on the order of 1 eddy-turnover time, while the phases change at a much slower rate with time scales of O(10) eddy turnovers. This is consistent with the behaviour that is exhibited byû z andθ in the non-dominant modes. Non-dominant modes are characterized by the mode numbers k below 2, and k greater than 3. The k = 2 mode is non-dominant at a time below 500t f , and the k = 3 mode is non-dominant at a time above 1000t f . The primary difference between {û z ,θ} and {û r ,û θ } groups of variables is that the {û r ,û θ } coefficients of the non-dominant as well as the dominant modes show no noticeable change in the behaviour of their phase and amplitude throughout the course of the simulation, and, thus, seem to be unaffected by the large-scale mode transition, as was previously observed in figure 11.
The distinction between the coefficients of the horizontal velocity components and u z andθ becomes less clear at higher wavenumbers. The temporal evolution of two additional modes (k = 5 and k = 10) is plotted in figure 13 to show how the phase and amplitude of the volume integrated coefficients are affected by increasing Fourier wavenumber. Figure 13 shows that as the wavenumber increases, so do the frequencies associated with the temporal evolution of the coefficients, consistent with the t k prediction in (5.3). The rapid oscillations and lack of distinction between the behaviour of the individual variables indicates that the high-wavenumber fields are becoming increasingly random in nature, and less descriptive of the physical system as a whole. Also, it is worth noting that the amplitudes of the high-frequency modes are higher for the vertical and radial velocities as opposed to the temperature and the azimuthal velocity. Stronger high-frequency contributions to the energy content for the vertical velocity rather than the temperature were also previously observed by Krug et al. (2019).

Integral time scale
In this section we attempt to quantify the temporal behaviour of the modes by measuring the classical integral time scale, T . Normally T is interpreted to be the coherence time of the flow field. It is defined in terms of the autocorrelation function, where ψ is a vector containing variables of interest and τ is the temporal offset between the two instances of the flow field, or snapshots. Usually ψ is taken to be the turbulent velocity field (ψ = {u r , u θ , u z }), where the prime indicates the fluctuating portion of the velocity. An autocorrelation based on this particular vector will determine a correlation based on the turbulent kinetic energy. A summation over the indices i = 1 . . . d, where FIGURE 13. Temporal evolution of the scaled volume integrated Fourier coefficients ofû z ,θ (left pane) andû r ,û θ (right pane) for k = 5 (a-d), k = 10 (e-g) plotted by amplitude (left) and phase (right). The phase plots have been rescaled to cover a period of 4π to highlight the low frequency cycles that occur in the temporal evolution of these modes.
d is the dimension of the vector, is then implied in the definition of the autocorrelation function in (5.4). The interest of this work includes the global correlation times as well as the time scales of the individual Fourier modes, since the large-scale structures have been shown to contain a relatively small number of Fourier modes. In terms of the Fourier coefficients averaging (5.5) over θ gives us R ii (r, θ, z, τ ) where δ k,−k is the Kronecker delta function, in which statistical stationarity in θ requires k = −k . Since all the flow variables are real signals, the negative Fourier mode −k can be expressed as the complex conjugate of the positive Fourier mode k. Therefore, the Dirac delta function in (5.6) shows that all wavenumbers will contribute to the correlation when multiplied by their complex conjugates. This also ensures that the correlation will be comprised entirely of real numbers which is required since the dependent variables are all defined in real space. The discrete representation of (5.6) is where N θ is the number of samples for the Fourier transform in the θ direction and * indicates the complex conjugate. Defining an autocorrelation function per wavenumber and recognizing that R ii (r, −k, z, τ ) = R * ii (r, k, z, τ ), (5.9) which follows from the corresponding properties of the Fourier coefficients (Canuto et al. 1988), one can see that (5.7) can be rewritten as where Re{R ii (r, k, z, τ )} is the real part of the autocorrelation (5.8). Utilizing (5.7) and (5.8), the total integral time scale (T ) and the integral time scales for each wavenumber (T k ) are (r, z, 0) dτ, (5.11) (5.12) which is equivalent to The quantity under the integral in (5.11) can also be referred to as a normalized autocorrelation,R (5.14) These integral time scales depend on the field ψ that is chosen. ψ = {u r , u θ , u z }, We used ψ = {ϑ } and ψ = {u r , u θ , u z , ϑ } to perform proper orthogonal decomposition in the study of Bailon-Cuba et al. (2010) on data from turbulent RBC simulations at various Γ . We shall also, likewise, be referring to the norms of these fields as the turbulent kinetic energy, thermal energy and total energy, respectively. Thus, in the definition of the autocorrelation functions in (5.4)-(5.14), the values i = 1 : 3 are assumed for the turbulent kinetic energy, i = 4 for the thermal energy and i = 1 : 4 for the total energy. Figure 14(a-c) shows T for the entire field when ψ is defined as the turbulent kinetic energy, the turbulent thermal energy and the total turbulent energy, respectively. The r-z plots of T also give insight into the structure of the flow field by indicating which regions of the flow field have longer correlation times, and also by how much the correlation times vary. Figures 14(a) and 14(b) show very different behaviour between the correlation of the turbulent kinetic and thermal energy fields. The two fields have little overlap between regions with very long correlation times. The kinetic energy field has its longest correlation times in the boundary layer while the turbulent thermal energy field has its longest correlation times in the bulk region. Regions where the fluctuations change sign frequently tend to have lower correlation times while regions where fluctuations maintain the same sign for long periods of time have longer correlations.
Further understanding of these correlation times can be gained by reviewing the variance of the various components that comprise the correlation metrics. Variance of the velocity and temperature fluctuations is plotted in figure 15 to show where the strongest fluctuations most frequently appear. Variance is defined as azimuthally and temporally averaged turbulence fluctuations σ u i (r, z) = u i (r, θ, z, t) 2 θ,t . The peak variance of the velocity components are aligned with the shape of the large-scale structures. The horizontal components (u r and u θ ) have a strong variance in the boundaries where the large-scale roll cells will have the strongest horizontal velocity contributions. Likewise the vertical velocity (u z ) has a larger variance in the bulk region where the large-scale structures are principally updrafts or downdrafts. This is where the opposing plumes pass one another as they cross the layer depth (see figure 15d), and this appears to be the factor that is driving the correlation times of the turbulent kinetic energy field down in the bulk region, due to a prevalence of the opposite sign fluctuations in this region.
The variance of u z also shows a peak in the central region near r = 0 which can be attributed to the time period when k = 3 was the dominant structure, and the uncertainty as the lobes of the k = 2 structure dance back and forth across the centre of the convection cell.
The peak variance of the temperature near the boundaries in figures 15(a) and 16(b) highlight how well mixed the thermals are in the bulk region. Note that the variance of the total energy (figure 15f ) has its peak values in the thermal boundary layers, the same as the variance of temperature, showing a noticeable contribution of the thermal fluctuations into a variance of the total energy field. However, in the bulk region, the variance of the total turbulent energy is similar to the variance in kinetic energy. This also explains while integral time scales shown in figure 14 are similar for the total and kinetic turbulent energy fields in the bulk region, but different from the thermal energy field: the large integral time scales for the temperature in the bulk are generated by relatively small temperature fluctuations divided by a small variance (the denominator of (5.11)), which account for a negligible contribution to the time scales when added to the velocity field with larger values of the fluctuations and the larger variance in the bulk. However, an increase of integral times scales for the total turbulent energy on the side walls due to a thermal field contribution is pronounced. Figure 15 also shows the skewness of the temperature field defined as which is a measure of asymmetry of the probability distribution of a fluctuating temperature field about its mean. It can be seen that the skewness of temperature peaks just as the variance begins to increase, and the two quantities are anti-correlated as they approach the wall (see figure 16b,c). The skewness profile shows that even though the fluctuations about the mean are very small in the bulk region, the lower half of the domain is still biased toward warmer fluid and the upper half toward colder. When the integral time scales are averaged over the volume, the differences between the three metrics narrows (see table 1). Following (2.16), volume integrated integral time scales (scaled with 2π) are defined as The turbulent kinetic energy shows the correlation time averaged across all Fourier modes that is about 50 % longer than the turbulent thermal energy, but the total turbulent energy shows about the same correlation time as the turbulent thermal energy. This means that there are regions where correlation among the velocity components cancels out correlation from the thermal region. Table 1 and figure 17 further demonstrate the effect of the dominant Fourier modes where the scaled volume integrated integral time scales are shown per mode number, {T k } V/2π . Here it can be seen that correlation times for the low-order modes, and specifically the k = 2 mode increase in magnitude in all the presented metrics. It is also noteworthy to see that the correlation time of the k = 2 mode is dramatically larger than the other modes with the second longest mode (k = 0 for thermal energy and k = 3 for kinetic and total energy) being almost five times shorter. Table 1 and figure 17 clearly show that the low frequency modes are responsible for the large integral time scales, and that the majority of the correlation time can be attributed to the k = 2 mode. Further evidence of this can be seen by comparing the spatial distribution   ). Figure 17 also shows that T k decays to a value of approximately 3t f for all three energy vectors after the first 12 Fourier modes. This is the minimum limit that can be obtained with this dataset since the snapshots were sampled 3t f apart. Shorter T k 's are probable for the higher wavenumbers.

Effects of the spatial inhomogeneity
In the previous sections the effects of spatial inhomogeneity have been seen in the mean flow and integral time scales. In this section the effects of spatial inhomogeneity will be analysed more carefully by looking at the r-z variations in the normalized autocorrelation, the integral time scale, T , and the Fourier spectra.

Spatial variability of the integral time scales
In the previous section we presented the integral time scales for three different quantities. Moving forward we will restrict our analysis to the integral time scale of the total turbulent energy since it is an aggregate of the other two. The normalized autocorrelation for the total turbulent energy field,R ii (r, z, τ ), i = 1 : 4, is plotted versus snapshot spacing in figure 18 at a selection of points in the r-z plane. While the entire simulation time is over 3000t f , only about half of that period is used to calculate the autocorrelations, due to the large values of the maximum time separations.
By definition (see (5.11)), T (r, z) is equal to the area under each of the plots in figure 18(a). Figure 18(a) shows that while fluctuations at points C and D in the highly correlated viscous boundary layer monotonically decay, they remain correlated over the entire dataset. Linear extrapolation of the lines in figure 18 can be used to estimate the time it will take forR ii (r, z, τ ) at points C and D to reach a value of zero. This time turns out to be approximately 90-100 eddy turnovers (≈2700-3000t f ). Based on the measured rotation rate of 1.1 degrees per eddy turnover, a rotation during the decorrelation time is ≈90 • , corresponding to the k = 2 mode changing from positive to negative vertical velocity. The updraft and downdraft of the k = 2 mode will exactly cancel one another. This suggests that decorrelation of the most persistent structures in this flow is caused by the observed, global rotation. The other two probes are taken at the mid-plane. The probe at point B is at a local minimum in T and shows sufficient decay inR ii (r, z, τ ) to indicate that the values become uncorrelated during this computation. The other probe at point A is near a local maxima in T . It shows signs of a weak long-lived transient as the correlation decays to zero with a separation time of approximately 800t f , but then begins to grow again. Some possible sources for this transient include the low-amplitude, low-frequency transient in the magnitude ofθ for k = 0 (see figure 9a), which can potentially be associated with a time scale of updraft and downdraft reversals (Sakievich et al. 2016). Note also the higher-frequency oscillations in the correlation function at all the probes A, B, C and D. The time scale of these higher-frequency processes stays very close between all the probes and varies between 84-94 t f , which corresponds to roughly 3 eddy turnovers. This time correlates well with the mode turnover time for the high-order modes that converges to a scale of an eddy turnover. Since the correlation functions here include the representation of all the modes, it is probable that higher-order modes are responsible for these high-frequency oscillations. Similar oscillations in a correlation function have been observed in Xi et al. (2006) and Mishra et al. (2011) in unit aspect-ratio cells with a time scale close to an eddy-turnover time. This suggests a similar origin of these oscillations in different aspect-ratio cells, coming from the contribution of high-order modes which are independent of geometry. Additional insight into the spatial variance of T k (r, z) can be found by investigating the contribution from the individual Fourier modes. Recall that R ii (r, z, τ ) can be defined as a summation of R ii (r, k, z, τ ) over all k's. A T k (r, z) field (5.13) can be calculated for each Fourier mode giving an indication as to how the individual modes contribute in the total correlation. Plots of the total turbulent energy based T k (r, z) for a selection of Fourier modes is provided in figure 19. Note that the values of the integral time scales for some modes can be negative, which explains why the modal values at some locations exceed their cumulative value in figure 14(c). Only low wavenumber plots are included in figure 19 due to the short correlation times of modes k > 10 (see figure 17 and table 1).
One observation of the subplots in figure 19 is that the globally dominant, highly correlated modes (subplots (c) and (d), and also subplot (a)) show a high level of symmetry about the mid-plane with long correlation times covering a substantial portion of the r-z plane, but the other modes do not. Subplots (b), (e) and ( f ) also show much smaller peak values for T k . The long correlation times in the dominant modes support our previous observations that the k = 2 and k = 3 modes are the primary sources for the long correlation times, and the high level of symmetry indicates that the shapes of these modes do not see much spatial variation over time. It is interesting to see that the peak in the correlation time T k for the mode k = 0 (subplot (a)) is the same as for the mode k = 3, but the spatial locations of the highly correlated events are different. This supports an earlier observation that the modes k = 0 and k = 3 coexisted when the mode k = 3 was strong, where mode k = 0 was dominant in the central region and the mode k = 3 was responsible for the formation of the side plumes. The lack of spatial symmetry and smaller range of T k in the non-dominant modes indicate that these modes see a large amount of variation in their spatial structure and could contain rare energetic events in the flow field which have life spans much less than the length of the simulation, but much longer than the sampling rate of 3t f .

Statistics of the Fourier modes and their spatial variability
In this section the temporally averaged statistics is presented for the Fourier modes which allows one to judge about the spatial variability of modes and the contribution of different length scales into the overall flow structure. This is achieved by evaluating the time-averaged energy spectra of the modes at different r-z locations. Up to this point in the paper all data has been presented with respect to the azimuthal Fourier modes, which, as discussed in § 2.3, are simply the integer mode indicators and are not directly related to the structure sizes. Therefore, seven different locations is provided iexamining Fourier coefficients at different radii corresponds to different physical length scales and energy densities per unit length. A more consistent way to compare the flow structure at various locations in the flow field is to normalize the energy spectra and frequency with respect to a geometric length scale, which is the wavelength λ k (r) = 2πr/k defined in (2.13). This is done by premultiplying the energy spectra with the radial location and plotting against the inverse of the wavelength 1/λ k (r). A sampling of the spectra at seven different locations is provided in figure 20. These locations are at various points within the boundary layers (bottom plate and side walls), and bulk region of the flow field; regions where different physical phenomena dominate. Here z = −0.45 and r = 3.1 are within the viscous boundary layers for the bottom and side walls, respectively, while z = −0.4 is just outside the viscous boundary layer in the vertical direction.

Variations in radial location
Spectra of ϑ, u θ and u z evaluated at various radial locations with a fixed height z = −0.4 (figure 20a, c, e and g) show a good collapse across virtually all length scales. For length scales that are greater than Γ = 6.3 (k/2πr 2 · 10 −1 ), ϑ and u θ collapse less well, and this behaviour is also seen in the u r plot. From the plots, it can be concluded that the effect of the side wall boundary conditions is rather small and is mostly prominent in the u r variable. It is not unexpected, since u r analytically must tend to zero at the wall. In fact, the same can be said about the centreline location, where u r is, again, analytically zero, and a lack of collapse in u r is again observed at low wavenumbers. The same behaviour is manifested in u θ . Other than u r , the side walls influence the large length scales in all the variables but u z , which shows an almost perfect collapse across all radii.
Failure to collapse in the larger length scales can be attributed to the dominance of low-order Fourier modes that describe the flow field's large-scale structure. Since Fourier mode k = 2 contains a large amount of energy throughout the entire domain it will disrupt the collapse of the spectra by affecting different length scales at each radii. 10 -2 10 -2 10 -1 10 -1 10 0 10 0 10 1 10 2 10 1 10 2 10 3 10 4 10 -2 10 -2 10 -1 10 -1 10 0 10 0 10 1 10 2 10 1 10 2 10 3 10 4 10 -2 10 -2 10 -1 10 -1 10 0 10 0 10 1 10 2 10 1 10 2 10 3 10 4 10 -2 10 -2 10 -1 10 -1 10 0 10 0 10 1 k /(2πr) k /(2πr) 10 2 10 1 10 2 10 3 10 4 In fact, if the spectra were to collapse across all length scales for all variables then it would be horizontally homogenous as in the canonical form of RBC with infinite Γ . In a sense the side walls of the convection cell act as a high pass filter because they limit the size of the largest length scales that can be observed in the flow. The fact that the k = 2 mode dominates the energy spectra at multiple length scales indicates  that the underlying structure has a modal nature, and that it is the principle cause for radial inhomogeneity. This is most likely due to the confining, geometric effects of the cylinder. Table 2 illustrates this point by providing the time-averaged energy values for the first several Fourier coefficients at several different radial locations at z = −0.4. For a reference, the scaled volume integrated value of |ϑ k | 2 corresponding to the values in figure 3 is also provided. This data complements the data in figure 20 for z = −0.4. Table 2 shows that the k = 2 mode is indeed energetically dominant at two of the three tabulated radii, and is responsible for the peaks in energy observed in figure 20. The k = 0 mode is approximately 40 % more energetic at r = 1.0, z = −0.4. This is most likely due to the residual effects of the central updraft early in the time series combined with geometric effects. Careful observation of figure 20(a) shows that as the radius decreases, the strength of the energy peak associated with the k = 2 mode decays. Due to a singular nature of approaching r = 0, an azimuthal alignment of the most dominant energetic structures might break, contributing its energy to an azimuthally invariant k = 0 mode. However, it is worth noting that the k = 2 mode is still larger than the k = 1 and k = 3 modes at r = 1.0, z = −4.0, thus breaking a monotonic decay of the energy spectrum and displaying a level of dominance as the second most energetic mode.

Variations in vertical location
When spectra are sampled at various heights at a fixed radius, the behaviour is virtually opposite to the fixed height, varying radius situation (see figure 20b, d, f and h). In the previous case ϑ and u z 's spectra showed the best collapse, but when the vertical location is varied their collapse is considerably worse than u r and u θ . Additionally, u r and u θ show the best collapse at the lowest frequencies, and a poorer collapse at higher frequencies.
In fact, divergence at high frequencies is seen for all three velocity components and the temperature, and their energy content decreases as the vertical position approaches the mid-plane. This is because there are more small-scale fluctuations for all the variables in the boundary layer region.
The spectrum of ϑ shows a strong collapse at the frequency associated with the k = 2 Fourier mode, indicating the dominance of the k = 2 mode in the temperature field across the entire depth of the domain but a decay with increasing distance from the wall for all other frequencies, commensurate with the presence of stronger high-frequency fluctuations in the near-wall region. A persistence of a dominant low-order mode all the way down to the wall presents evidence of the influence of the large-structure organization on the boundary layer flow, observed in previous studies (Sakievich et al. 2016;Pandey et al. 2018;Krug et al. 2019).
The spectrum for u z in figure 20(h) shows some special characteristics that deserve a discussion of their own. Perhaps the most notable is that at the mid-plane the smallest energy among the heights occurs at high frequencies and the largest energy among the heights at lower frequencies. The point in the spectrum where the energy in the mid-plane is no longer smallest occurs at a non-dimensional frequency of approximately 8.5 corresponding to a physical length scale of 0.118H. For lower frequencies (larger length scales), there is a region where the spectrum collapse for the vertical positions that are outside the viscous boundary layer. This region of collapse starts to break apart at a non-dimensional frequency of 2 corresponding to a length scale of 0.5H, and the energy in frequencies lower than 2 increases with the distance from the wall. The absence of a clear peak in the vertical velocity spectrum for the k = 2 mode can perhaps be explained by more effective mechanisms of energy transfer from large to small scales in this field, which also amounts to a larger amplitude of higher-order modes in the vertical velocity (as opposed to, e.g. temperature) observed in figure 13. A similar phenomenon of a higher small-scale contribution to the energy associated with the vertical velocity fluctuations as opposed to the temperature is discussed in Krug et al. (2019). In general, a collapse of the temperature spectrum in the large scales across z planes, and the vertical velocity spectra in the intermediate scales, but not large or small, is commensurate with the findings of Krug et al. (2019).

Discussion and conclusions
It has been shown that Γ > 1 turbulent RBC has dynamics that occur on much longer time scales and affects more spatial Fourier modes than RBC in a Γ = 1 cell. A general explanation for the increase of time scales in wider aspect-ratio cells can be provided that is due to the increase of the length scales of the coherent motions that are able to settle in larger aspect-ratio domains. The long correlation time scales of coherent structures observed in the current study resonate well with the recent studies of Pandey et al. (2018) who examined the evolution of times scales of turbulent superstructures in the square domain with Γ = 25. While previous works on the large-scale patterns in Rayleigh-Bénard convection in wide aspect-ratio cells were primarily concerned with the statistical analysis of the flow field, and, thus, the properties of the 'average' flow structure (Hartlep et al. 2003;Pandey et al. 2018;Stevens et al. 2018;Krug et al. 2019), the current study focuses on the individual structure and a detailed analysis of its temporal dynamics. This is achieved via investigating each spatial Fourier mode independently. The individual modal analysis shows that the correlation times of the dominant Fourier modes are substantially longer than the other modes and that these individual modal correlation times can exceed the total correlation time of the system (see table 1). Furthermore, the spatial variation of integral time scales shows an alignment between the location of maximum temporal correlation for the entire system (figure 14) and the individual modes ( figure 19). This is a noteworthy observation since the integral time scales have been shown to vary by three orders of magnitude across the r-z plane ( figure 14). Examination of the temporal correlation in the regions where the integral scale is largest provides a further connection between the system's behaviour and the dominate modes. The rate at which temporal correlation decays was found to match the measured rate of rotation for the most energetic Fourier mode (k = 2) at 1.1 degrees per eddy turnover. From these observations we can conclude that the dominant Fourier modes leave a strong signature on the temporal scales of the system.
While mode 2 has clearly emerged as being a dominant mode after approximately 1000 t f of the simulations, it is not the only mode that influences the overall organization of the structure. In general, the dynamics of the large-scale structure in the current study exhibits a complex interplay between several low-order modes. For example, the transition of dominant mode from k = 3 to k = 2 occurs gradually over approximately 1000t f . This event can be observed in the time-averaged field (figure 5) and through careful observation of instantaneous fields (figure 8), but becomes blatantly obvious when analysing the dynamics of the individual Fourier modes (see figures 6,7,9). Interestingly, dynamics of the first, k = 1, Fourier mode also play a role in the overall organization of the structure and the mode interplay observed in the current study, at least for the first 1000t f , including the rotations and the cessations of the k = 1 mode, similar to k = 1 dynamics in the unit-Γ cells (Brown et al. 2005;Mishra et al. 2011). However, this mode never dominates the overall structure dynamics, as opposed to the Γ = 1 case.
The dominance of modes 2 and 3 in a Γ = 6.3 aspect-ratio cell can potentially be explained by considering the 'natural' sizes of the superstructures, reported to be of the order of 6-7 of the domain height as deduced from numerical studies in the domains with Γ = 10-60 and Pr ∼ 0.7 (Hartlep et al. 2003;Pandey et al. 2018;Stevens et al. 2018), while the sizes up to 10H can be expected for Pr = 6.7 (Busse 1994;Pandey et al. 2018). The sizes reported in these studies, however, correspond to the spectral wavelength, which equals Λ k = πΓ /k (normalized with H) according to (4.1), thus giving k ∼ 2-3 when the sizes of Λ k = 6-10 are substituted into this equation for Γ = 6.3. Following the same logic, in a unit aspect-ratio cell, the longest mode that can settle (k = 1) gives the wavelength of Λ 1 = πΓ ∼ 3, which is still smaller than the size of a natural superstructure, thus explaining why the mode k = 1 is clearly dominant in a unit aspect-ratio case, since higher-order modes would have even smaller wavelengths. In this sense, the existence of the k = 1 mode in the current Γ = 6.3 case is interesting, since it manifests an existence of the correlated structures of even larger length scales than an average size of the superstructures.
The concept of a 'mode turnover time' has been introduced in this study as an initial hypothesis for explaining the interactions between the mode number, Γ , and the time scales. Conceptually, a mode turnover is the time it would take a particle moving at the r.m.s. velocity to traverse a modal structure at a given Γ . The observed time scales of the mode transition, on the order of 20-30 eddy turnovers, were shown to correlate well with the introduced concept of a mode turnover time. Moreover, the concept of the 'mode turnover time' also predicts similar time scales associated with the duration of destabilization events associated with reorientations observed in Γ = 1 cells; see Brown et al. (2005) and Mishra et al. (2011). We note that the time scales are not expected to grow indefinitely as the aspect ratio increases, due to a saturation of the sizes of the superstructures. Indeed, time scales identified in larger aspect-ratio studies (Emran & Schumacher 2015;Pandey et al. 2018) are similar to the time scales observed here, and not 5 to 10 times larger which would be commensurate with the considered aspect-ratio sizes.
We also propose to separate the influence of 'fast' and 'slow' time scales on the processes observed in turbulent RBC at Ra = 9.6 × 10 7 and Pr = 6.7. 'Fast' time scale correlates with the mode turnover time, while 'slow' time scale is based off the diffusion, or viscous, time scale. Viscous time scale can be estimated as t v = √ Ra/Pr t f (Pandey et al. 2018), equal to t v ∼ 3800t f ∼ 120t in the current case, while the diffusion time scale t d = √ Ra Pr t f ∼ 25 460t f ∼ 804t is even larger. The current results indicate that the effects associated with the mode transition, as well as the fast rotation observed in the first mode both scale with the mode turnover time ('fast' time scale), while the slow azimuthal drifts, reported previously for unit aspect-ratio cells (Brown et al. 2005;Brown & Ahlers 2006) and also observed here in a mode 2 once it stabilized, occur on time scales that are at least an order of magnitude larger ('slow' time scale). We hypothesize that the difference in the time scales can be explained by the difference in the physical mechanisms that cause these events. The slow rotations, or drifts, are associated with the slow diffusive processes (Brown & Ahlers 2006, while the fast rotations, as well as the mode cessations and transitions, are related to the destabilization processes, which has been previously linked to the interaction between buoyancy and friction (Sreenivasan et al. 2002;Brown & Ahlers 2007). In a conceptual model of an LSC reversal by Sreenivasan et al. (2002), the reversal is explained as a loss of equilibrium in the dynamics associated with the ascending and descending plumes within the LSC circulation cycle. Their model fits very well with the current hypothesis that such a process would occur on a mode circulation time scale. Previous studies demonstrated a stochastic nature of the reorientation processes, caused by the perturbations. These perturbations, likely appearing locally, need to propagate through an entire mode to cause a mode destabilization, which would require a time scale of the order of the mode turnover time. We would like to stress once again that the mode turnover time scales, of the order of 2.75t for k = 1 in Γ = 1 cells, correspond to the duration of the transition events (such as an LSC reorientation), from destabilization to restabilization, following the origination of the perturbation. The origin of the perturbations themselves is stochastic, so the time scales between events are much larger, e.g. reported to be around 10 t to 30 t in Γ = 1 cells (Sreenivasan et al. 2002;Brown et al. 2005;Mishra et al. 2011). Brown & Ahlers (2007) elaborated on the physical processes accompanying the dynamics of the LSC reversal by showing that a destabilization can also lead to an onset of a fast azimuthal motion, ultimately responsible for the structure reorientation. The reason for the increased rotation rate during the mode destabilization is a reduction of an angular momentum of a weakened mode associated with a reduction in the mode's amplitude (Brown & Ahlers 2007). Both cessations and rotations were shown to be accompanied by this fast azimuthal drift (Brown et al. 2005;Brown & Ahlers 2007;Mishra et al. 2011), the difference being that the mode amplitude essentially vanishes during the cessation, but stays finite (although often reduced) during the rotation-led reorientations. Similar cessation-like events accompanied by fast rotations of the modes were observed twice in mode 1 and once in mode 2 in this study, all during the global structure transition process. Even though cessations are perceived as instantaneous events, the processes that lead to it and follow it until a mode stabilization is completed are finite scale, evidence of this can be seen in time series presented in Brown et al. (2005), Mishra et al. (2011) and Zürner et al. (2019), as well as in the current DNS results. It is conjectured here that both cessation-led and rotation-led reorientation events have a similar origin and operate on similar time scales, and are, in fact, just different manifestations of the same destabilization process. Moreover, it is likely that, during the process of destabilization, multiple rotations, cessations (double cessations), as well as the mode transitions in wider aspect-ratio systems, can occur simultaneously as a signature of the same destabilization process that leads to an eventual reorganization of the structure. Since destabilization events are caused by stochastic perturbations, the duration between them is random (Sreenivasan et al. 2002;Brown & Ahlers 2007). We were only able to observe one destabilization event over the 3000 t f duration of the simulation. Note, once again, that due to a high aspect ratio, the relevant times scales are several times larger in our system than in unit aspect-ratio cells, which significantly reduces the probability of observing the reorganization events in higher aspect-ratio systems, experimentally, or numerically.
The current study is also in a position to contribute to the discussion initiated in the previous works of Pandey et al. (2018), Stevens et al. (2018) and Krug et al. (2019) on whether the vertical velocity and temperature fields correspond to the signature of the same large-scale structure, or whether there are structures of different sizes that are formed by the temperature and the vertical velocity components. The results in the current study and in our prior work (Sakievich et al. 2016) support the conclusion of Krug et al. (2019) that ϑ and u z correspond to the same structure for several reasons: (i) temporal dynamics of the low-order Fourier modes of temperature and vertical velocity is highly correlated, which illustrates that they go through the same processes of transition and reorganization testifying their link to the same large-scale structure; (ii) the imprint of the low-order modes on the spatially and temporally averaged spectra of ϑ and u z is similar and seems to produce a clear peak at the k = 2 wavenumber, albeit this peak is stronger in temperature than in vertical velocity; (iii) visualizations in our previous work (Sakievich et al. 2016) show a high degree of coherence between the spatial location of the thermal updrafts and downdrafts, and the velocity roll cells identified by a visualization of the three-dimensional velocity field, illustrating that both temperature and velocity fields are effected by the same large-scale organization. Similar to the study of Krug et al. (2019), we also find that there are more small-scale fluctuations in the vertical velocity rather than in the temperature field, manifested by the fact that the amplitudes of the higher-order modes are larger for the vertical velocity than for the temperature manifested in figure 13. This abundance of small-scale energy contribution in the vertical velocity is responsible for the shift of the spectral peak towards higher frequencies compared to the temperature in the works of Pandey et al. (2018), Stevens et al. (2018) and Krug et al. (2019). Similarly, it results in weakening of the k = 2 peak in the volume integrated spectra of u z in the current work (see figure 3), and eliminating it from the local (in r, z) azimuthal spectral plots (see figure 20). Again, commensurate with the findings of Krug et al. (2019), we see a good collapse of the spectra of the temperature in the low-order modes across the vertical planes in the RBC cell (also showing a clear k = 2 peak across all vertical locations), while the vertical velocity collapses in the intermediate scales, but not in large and small scales. This striking similarity of the current data and the results of Krug et al. (2019) testifies to the similar principles of the spatial organization of structures between the current Γ = 6.3 case and the superstructures found in the larger domains.
While the length scales, time scales and the principles of spatial organization similar to the properties of the superstructures have been observed in the current work, cylindrical geometry and the side wall boundary conditions in the current Γ = 6.3 case do influence the organization of structures observed in this study. It is mostly manifested via the fact that the structures organize themselves in line with the azimuthal Fourier modes. It also influences the length scales in the core of the cylindrical cell, which do not correspond to the same size motions, but rather to the same wavenumber motions, which shortens the physical length scales as one tends toward the centre of the cell (see figure 20). Similar azimuthal mode organizations have been observed for large-scale-motions and very-large-scale-motions (VLSM), otherwise known as superstructures, in pipe flow. For example, studies of Bailey & Smits (2010) and Baltzer et al. (2013) showed that VLSM have large streamwise scales that concentrate around a single azimuthal mode, with the dominant azimuthal mode being k = 3 in both studies. It remains to be answered whether this concentration of energy towards the same azimuthal modes found in pipe flow and in the current RBC case has far-reaching consequences, or whether it is a pure coincidence, given that the k = 2, 3 modes are likely emerged in the current study due to a spatial fit of the natural RBC superstructures into the given cylinder size. Further studies of the mode dynamics in cylindrical and other shape domains in the turbulent Rayleigh-Bénard convection with even larger aspect ratios would be of interest in this respect, so that the principles via which superstructures are organized geometrically and are evolved dynamically can be investigated with a minimum influence of the confining geometry.
Finally, since only a single regime of Ra and Pr numbers has been investigated in the current study, it naturally leads to a question: what is expected to change, and what will likely stay the same, for different parameter values? It can be said almost for certain, that the time scales when defined with respect to a free-fall time will change when other Ra and Pr are considered (Pandey et al. 2018). However, the mode destabilization time scales as defined through the eddy-turnover time, see (5.3), and the role of the aspect ratio and the mode number in the global structure dynamics across different Ra and Pr is expected to be robustly represented by this scaling. The reason is that the dependence on Ra and Pr is already encoded in the r.m.s. value of the turbulent velocity fluctuations u 2 z V,t , while the role of the individual mode dynamics with respect to u 2 z V,t is solely reflected by the mode circulation length, which is only a function of k and Γ . Indeed, a robust scaling of ∼2 eddy turnovers for the duration of the LSC reorientation process in Γ = 1 cells was observed, at least, for the ranges of Pr numbers from 0.029 (Zürner et al. 2019) to 0.7 (Mishra et al. 2011) to 4.4 (Brown et al. 2005. The time scale of the azimuthal meandering (slow azimuthal rotation) is likely to change with Pr, both in free fall and eddy-turnover scaling, favouring longer rotation time scales at either very low or very high Pr, due to a large difference between viscous and diffusion time scales (Pandey et al. 2018;Zürner et al. 2019). A spatial organization of the modes in a container of a given size will also likely change with Ra and Pr. This is due to different length scales of the favoured structures observed at different Ra and Pr regimes (Hartlep et al. 2003;Pandey et al. 2018). Due to a difference in length scales, the mode numbers which will settle in a given container will be different, which will determine the overall appearance and dynamics of the global structure. Additionally, apart from the length scales, a striking difference of the spatial convection patterns at different Ra and Pr numbers might also play a role. For example, at a low Pr, the Rayleigh-Bénard convection was found to be dominated by rolls, or elongated 'rivers', while at high Pr, the convection pattern takes a form of cells connected by ridges (Malevsky 1995;Breuer et al. 2004;Pandey et al. 2018). Interestingly, the current Pr = 6.7 case is near a transitional regime, where the convection pattern changes from rolls to cells, and, according to Hartlep et al. (2003), both styles of convection may coexist. It would be interesting to explore if the 3 to 2 mode transition observed here and the disappearance of the central column (k = 0 mode) might be related to a potential switch in a convection pattern.
Extension of the presented study to other parameter regimes and other aspect ratios would be a natural important direction for future work. It must be realized, however, that the time scales of interest, considering the findings of the current study and similar works, are extremely long. For example, it would be of interest to increase the temporal duration of the current DNS by another decade, to observe even longer-term dynamics of the structure which we just started to uncover, such as: will the structure be destabilized again, and how a newly formed structure will look? Will it transition back to a 3 mode? Will a central column reappear or reverse its direction? Will a structure rotate around a full circle? Similar long-time studies must also be performed at different Γ , Ra and Pr regimes. In this context it would be important to evaluate whether other, lower-fidelity, but computationally more efficient approaches, such as large eddy simulations, would be able to predict the important dynamical events in Rayleigh-Bénard convection. While formulation of reliable subgrid closure models for turbulent heat transfer problems remains challenging, the current physical problem presents a clear motivation and the need for these efforts to be undertaken on a large-scale basis.
original field. That is, if the original field's statistics had a bias due to preferential large-scale updrafts, the simulations starting from a transformed field will give statistics commensurate with the prevalence of the downdrafts (over the same run time). The sum of the two statistics produced an essentially unbiased statistics in very good agreement with carefully conducted long-time experiments (Fernandes 2001;Fernandes & Adrian 2002). Since the transformation defined by (A 1) commutes with the azimuthal and time averaging operators, in this paper, instead of running new simulations starting with the transformed field, we are applying a symmetry transformation onto a computed θ,t averaged field directly, so that the sum of the two statistics presented in § 3 is unbiased to the extent possible given the finite time of the simulations. Note that the transformation described here is geared towards removal of a bias due to a formation of the preferential thermal plumes and the associated large-scale vertical motions. There are other potential symmetries in the RBC problem, for example, associated with a preferential direction of an azimuthal rotation of the structure in the RBC cell, which the current transformation would not be able to mitigate.