1. Introduction
Thermal convection is a pervasive phenomenon that governs heat and mass transfer across a variety of processes. Examples include the ocean circulation (Rahmstorf Reference Rahmstorf2006), the convection flows in the inner core of the Earth and its mantle (Glatzmaier & Roberts Reference Glatzmaier and Roberts1995; Bercovici Reference Bercovici2003; Coltice, Gérault & Ulvrová Reference Coltice, Gérault and Ulvrová2017) and the cooling process of liquid metal in fusion reactor (Tagawa & Ozoe Reference Tagawa and Ozoe1997; Wang & Zhou Reference Wang and Zhou2019). To understand the generic mechanisms of these buoyancy-driven flows, an idealised system is that of Rayleigh–Bénard convection (RBC), i.e. a layer of fluid heated from below and cooled from above (for reviews, see, e.g. Ahlers, Grossmann & Lohse (Reference Ahlers, Grossmann and Lohse2009), Lohse & Xia (Reference Lohse and Xia2010); Chillà & Schumacher (Reference Chillà and Schumacher2012), Xia (Reference Xia2013), Xia et al. (Reference Xia, Huang, Xie and Zhang2023)). This system is controlled by three dimensionless parameters: the Rayleigh number
$Ra=\alpha g\Delta \textit{TH}^3/\nu \kappa$
, measuring the strength of buoyancy forcing; the Prandtl number
$\textit{Pr}=\nu /\kappa$
, characterising the relative diffusion of momentum and heat; and the aspect ratio
$\varGamma =l/H$
, quantifying the geometrical influence. Here,
$\Delta T$
is the temperature difference between two horizontal plates with a vertical distance
$H$
,
$g$
is the acceleration due to gravity and
$\alpha$
,
$\kappa$
and
$\nu$
represent the coefficients of thermal expansion, thermal diffusivity and kinematic viscosity of the working fluid, respectively. The lateral dimension
$l$
of the system could be either the diameter
$D$
for cylindrical cells or the combination of length
$L$
and width
$W$
for rectangular cells.
One noteworthy flow structure in turbulent RBC is the large-scale circulation (LSC), which originates from self-organised thermal plumes and extends across the system (Krishnamurti & Howard Reference Krishnamurti and Howard1981; Xi, Lam & Xia Reference Xi, Lam and Xia2004). Its crucial role in heat and mass transfer has prompted numerous studies to explore its properties. A distinctive aspect of the LSC is its low-frequency oscillation, initially observed in the Chicago Experiment (Heslot, Castaing & Libchaber Reference Heslot, Castaing and Libchaber1987; Castaing et al. Reference Castaing, Gunaratne, Heslot, Kadanoff, Libchaber, Thomae, Wu, Zaleski and Zanetti1989). Specifically, based on local measurements, Heslot et al. (Reference Heslot, Castaing and Libchaber1987) reported the re-emergence of coherent oscillation in the temperature signals when RBC gets into the turbulent regime, indicating the presence of a large-scale coherent flow within the system. Later, Sano, Wu & Libchaber (Reference Sano, Wu and Libchaber1989) established a direct link between the temperature oscillation frequency and the turnover time of the mean flow, suggesting that the oscillation is associated with the LSC. This relationship was supported by subsequent studies of turbulent RBC using various working fluids and convection cells (Takeshita et al. Reference Takeshita, Segawa, Glazier and Sano1996; Ashkenazi & Steinberg Reference Ashkenazi and Steinberg1999; Niemela et al. Reference Niemela, Skrbek, Sreenivasan and Donnelly2001; Qiu & Tong Reference Qiu and Tong2001, Reference Qiu and Tong2002; Qiu et al. Reference Qiu, Shang, Tong and Xia2004; Sun & Xia Reference Sun and Xia2005; Brown & Ahlers Reference Brown and Ahlers2006b ; Niemela & Sreenivasan Reference Niemela and Sreenivasan2006; Zhou, Sun & Xia Reference Zhou, Sun and Xia2007).
Further studies employing multi-point velocity measurements and shadowgraph imaging revealed that the LSC exhibits a three-dimensional pattern: thermal plumes near the top and bottom plates oscillate at a frequency close to the mean-circulation frequency, but out of phase in the horizontal direction (Qiu et al. Reference Qiu, Shang, Tong and Xia2004; Funfschilling & Ahlers Reference Funfschilling and Ahlers2004). This oscillatory behaviour is known as the torsional mode, due to its twisting nature (Xi, Zhou & Xia Reference Xi, Zhou and Xia2006; Funfschilling, Brown & Ahlers Reference Funfschilling, Brown and Ahlers2008). Later, the torsional mode is found to coexist with the sloshing mode, which is characterised by a periodic off-centre motion in the bulk fluid (Xi et al. Reference Xi, Zhou, Zhou, Chan and Xia2009; Zhou et al. Reference Zhou, Xi, Zhou, Sun and Xia2009). These two oscillation modes are dominant in cylindrical convection cells with
$\varGamma \simeq 1$
. Recently, a novel three-dimensional oscillation mode, characterised by the orbiting motion of the vortex centre of the LSC, was observed in convection cells with
$\varGamma \gt \sqrt {2}$
(Vogt et al. Reference Vogt, Horn, Grannan and Aurnou2018). This oscillation mode, referred to as the jump rope mode, has been confirmed experimentally and investigated thoroughly by numerical simulations (Li et al. Reference Li, Chen, Xu and Xi2022; Horn, Schmid & Aurnou Reference Horn, Schmid and Aurnou2022; Cheng et al. Reference Cheng, Mohammad, Wang, Keogh, Forer and Kelley2022; Teimurazov et al. Reference Teimurazov, Singh, Su, Eckert, Shishkina and Vogt2023; Wondrak et al. Reference Wondrak, Sieger, Mitra, Schindler, Stefani, Vogt and Eckert2023).
Decades of research have led to several physical mechanisms being proposed to explain the oscillation of the LSC in turbulent RBC (see table 1). Since the LSC is formed by self-organised thermal plumes, its oscillation behaviour might be attributed to the plume emission dynamics. According to the model proposed by Howard (Reference Howard1966), a plume emits once a growing diffusive thermal boundary layer reaches a critical size, resulting in a plume emission frequency that scales with the
$Ra$
as
$Ra^{2/3}$
. Although this prediction is in good agreement with the behaviour of laminar thermal plumes (Qin et al. Reference Qin, Hou, He, Wei and Huang2024) and the finding in RBC without an LSC (Niemela et al. Reference Niemela, Skrbek, Sreenivasan and Donnelly2001), it deviates from the
$Ra$
-dependent scaling (close to 0.5) of the observed LSC oscillation frequency, indicating that plume emission alone cannot account for the oscillation of the LSC. In this regard, Villermaux (Reference Villermaux1995) proposed a model in which the LSC dynamics is coupled with boundary layers. In this scenario, the boundary layer instability is triggered by the incoming plumes carried by the LSC from the opposite boundary. As a result, there is a time delay between the emission of thermal plumes due to the cell-crossing time, leading to temporally alternating and periodic oscillation.
Physical mechanisms proposed to understand the oscillation of the LSC.
Although the Villermaux model was supported by some experimental observations (Qiu & Tong Reference Qiu and Tong2001, Reference Qiu and Tong2002; Qiu et al. Reference Qiu, Shang, Tong and Xia2004; Sun, Xia & Tong Reference Sun, Xia and Tong2005; Song, Villermaux & Tong Reference Song, Villermaux and Tong2011), the experiment conducted by Xi et al. (Reference Xi, Zhou, Zhou, Chan and Xia2009) in a cylindrical convection cell demonstrated that thermal plumes are emitted neither periodically nor alternately, but rather randomly and continuously, from the top and bottom boundary layers. This experimental finding challenged the key assumption of the Villermaux model. Inspired by the finding by Xi et al. (Reference Xi, Zhou, Zhou, Chan and Xia2009), Brown & Ahlers (Reference Brown and Ahlers2009) extended their previous stochastic model of the LSC (Brown & Ahlers Reference Brown and Ahlers2007, Reference Brown and Ahlers2008) to explain the observed torsional and sloshing modes. Based on the extended model, they further suggested that the torsional and sloshing modes can be interpreted as a pair of travelling waves consisting of counter-moving plumes, where the oscillations are advected along with the LSC. Yet, the travelling-wave pattern proposed in this advected-oscillation picture has not been revealed explicitly through experiments.
In addition to the plume-based models, there are other alternative viewpoints. For instance, by considering a dynamical model in the inviscid and unforced limit, Resagk et al. (Reference Resagk, du Puits, Thess, Dolzhansky, Grossmann, Fontenele Araujo and Lohse2006) found that both boundary layers and plumes are not necessary factors for the LSC oscillation. They argued that the oscillation properties are intrinsic to the bulk flow dynamics and can be regarded as ‘inertial oscillations’. However, their model only applies when the convection cell’s ellipticity exceeds a threshold, and does not account for the typical oscillation phenomena observed in other geometries.
Among the models discussed above, the Brown–Ahlers model and its extension (Brown & Ahlers Reference Brown and Ahlers2007, Reference Brown and Ahlers2008, Reference Brown and Ahlers2009) have achieved a great success in explaining most of the three-dimensional dynamics of the LSC (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009; Song et al. Reference Song, Brown, Hawkins and Tong2014; Ji & Brown Reference Ji and Brown2020). On the other hand, the recently discovered jump rope mode in larger
$\varGamma$
systems (Vogt et al. Reference Vogt, Horn, Grannan and Aurnou2018) suggests that the torsional and sloshing modes may not be requisite for observing the oscillation dynamics. Indeed, significant periodic oscillations in both temperature and velocity fields can be observed in a vertical thin disk, where the torsional and sloshing modes are impossible to excite in such a strongly confined geometry (Song et al. Reference Song, Villermaux and Tong2011). The persistence of oscillations, regardless of specific geometric configurations of the convection cells, coupled with the intricate interplay between the velocity and temperature fields, underscores the need for further investigation of these phenomena in different geometries.
In the present study, we investigate the LSC oscillation in a quasi-two-dimensional (quasi-2-D) turbulent RBC cell. Here, ‘quasi-2-D’ refers to the fact that the lateral confinement of the container suppresses the development of large-scale motions in the spanwise direction. This confinement reduces, but does not eliminate, the complexity associated with the three-dimensional dynamics of the LSC, and thereby enables a more focused examination of its oscillation mechanism. To reveal the flow structures responsible for the oscillation dynamics, we move beyond traditional single-point or two-point measurements that provide only limited access to the global organisation of the flow. Instead, we employ data-driven (DD) modal analysis to identify the dominant oscillatory structures. In particular, we employ spectral proper orthogonal decomposition (SPOD), which is an effective tool for identifying frequency-resolved coherent structures in complex datasets (Lumley Reference Lumley2007; Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Bres2018; Abreu et al. Reference Abreu, Cavalieri, Schlatter, Vinuesa and Henningson2020). Through SPOD analysis, we show that the structures associated with the dominant oscillation frequencies display an approximate azimuthal periodicity when viewed in polar coordinates. These coherent structures can be interpreted as travelling waves, in line with the advected-oscillation picture (Brown & Ahlers Reference Brown and Ahlers2009).
Furthermore, motivated by the links between SPOD, stochastic forcing models and resolvent analysis (Semeraro et al. Reference Semeraro, Jaunet, Jordan, Cavalieri and Lesshafft2016; Towne, Bres & Lele Reference Towne, Bres and Lele2017; Towne et al. Reference Towne, Schmidt and Colonius2018), we use a DD resolvent analysis (McKeon & Sharma Reference McKeon and Sharma2010; McKeon Reference McKeon2017; Herrmann et al. Reference Herrmann, Baddoo, Semaan, Brunton and McKeon2021) as a complementary input–output framework inferred from the data. Note that this approach is not used to extract spectral information that are unavailable from SPOD; rather, it is employed to examine whether the dominant oscillatory structures can be interpreted as the most amplified responses of an effective linear dynamics around the mean circulation. The applicability of this linear-amplification perspective for cylindrical RBC is also discussed.
The paper is structured as follows. Section 2 introduces the experimental set-up and the analytical methods. Section 3 is divided into four parts: § 3.1 examines the oscillation properties in the velocity field; § 3.2 assesses the spatial coherence of the flow through cross-correlation and cross-spectrum analysis; § 3.3 employs the SPOD method to identify the oscillatory structures; and § 3.4 applies the resolvent analysis to motivate the linear-amplification interpretation. Section 4 examines this interpretation further by using datasets from complementary numerical simulations of quasi-2-D RBC and velocity measurements in cylindrical RBC. Section 5 concludes with a summary of our findings and a discussion of future directions.
2. Experiments and methods
2.1. Experiments
The experiment was conducted in a rectangular convection cell, similar to those used in previous studies (Xia, Sun & Zhou Reference Xia, Sun and Zhou2003; Sun, Cheung & Xia Reference Sun, Cheung and Xia2008; Li et al. Reference Li, He, Tian, Hao and Huang2021). Briefly, the convection cell was constructed with a Plexiglas tank forming the sidewalls and two copper plates serving as the top and bottom boundaries. The dimensions of the cell are
$25.2$
cm in length (
$L$
),
$25.0$
cm in height (
$H$
) and
$7.6$
cm in width (
$W$
), resulting in aspect ratios of
$\varGamma _\perp = L/H \approx 1$
for the length to height and
$\varGamma _\parallel = W/H \approx 0.3$
for the width to height. The temperatures of both copper plates were monitored by thermistors embedded 2 mm away from the fluid–solid interface. Temperature regulation was achieved by using a refrigerated circulator for the top plate and a resistive wire heater for the bottom plate. By adjusting the temperature difference between the top and bottom plates,
$Ra$
was varied from
$2.5\times 10^9$
to
$1.2\times 10^{10}$
in the present study. Deionised and degassed water was used as the working fluid, with the mean temperature of the bulk flow kept at 40
$^\circ \text{C}$
approximately, corresponding to a Prandtl number of
$\textit{Pr} \approx 4.3$
. Before performing the velocity measurements, the convection cell was carefully levelled to an accuracy better than
$0.1^\circ$
and enclosed within a thermostat box to prevent the flow dynamics from being influenced by room temperature variations.
The velocity fields were measured in the vertical mid-plane of the convection cell, using a planar particle image velocimetry (PIV) system. Polyamide particles with a density of 1.03 g cm−3 and a diameter of 50 μm were seeded into the fluid and then illuminated by a green laser beam. The particle images were recorded by a front-view CMOS camera at a sampling rate of approximately 5 Hz over nearly two hours to ensure sufficient data for statistical analysis. The particle images were processed with PIVlab software (Stamhuis & Thielicke Reference Stamhuis and Thielicke2014; Thielicke & Sonntag Reference Thielicke and Sonntag2021), using sub-windows of 32 pixels
$\times$
32 pixels with a
$50\,\%$
overlap, resulting in
$93 \times 93$
vectors for each velocity field. The obtained velocity fields were used for analysis by following methods.
2.2. Methods
2.2.1. Cross-correlation and cross-spectrum
The time-correlation function
$\rho _{\textit{ab}}(\tau )$
is one of the commonly used methods to characterise the oscillation behaviours;
$\rho _{\textit{ab}}(\tau )$
is defined as
\begin{equation} \rho _{\textit{ab}}(\tau ) = \langle [a(t)-\langle a\rangle _t][(b(t+\tau )-\langle b\rangle _t]\rangle _t/\sigma _a\sigma _b. \end{equation}
Here,
$\langle \boldsymbol{\cdot }\rangle _t$
denotes temporal averaging, while
$\sigma _a$
and
$\sigma _b$
are standard deviations of time series
$a$
and
$b$
, respectively. When
$a=b$
,
$\rho _{\textit{ab}}$
is the auto-correlation function and we denote it as
$\rho _a$
. Under the Fourier transformation, the cross-spectrum can be obtained as
\begin{equation} \mathcal{S}_{\textit{ab}}(f) = \int \rho _{\textit{ab}}(\tau )e^{-i2\pi f\tau }d\tau = \mathcal{C}_{\textit{ab}}(f) + i\mathcal{Q}_{\textit{ab}}(f), \end{equation}
where
$\mathcal{C}_{\textit{ab}}(f)$
represents the cospectrum and
$\mathcal{Q}_{\textit{ab}}(f)$
denotes the quadrature spectrum (Stull Reference Stull2012). Then, the amplitude spectrum
$|\mathcal{S}_{\textit{ab}}(f)|$
, the coherence spectrum
$\gamma _{\textit{ab}}(f)$
, and the phase spectrum
$\phi _{\textit{ab}}(f)$
can be defined as
\begin{eqnarray} |\mathcal{S}_{\textit{ab}}(f)| = \sqrt {\mathcal{C}^2_{\textit{ab}}(f) + \mathcal{Q}^2_{\textit{ab}}(f)}, \end{eqnarray}
\begin{eqnarray} \gamma _{\textit{ab}}(f) = \sqrt {\frac {|\mathcal{S}_{\textit{ab}}(f)|^2}{\mathcal{S}_a(f)\mathcal{S}_b(f)}}, \end{eqnarray}
\begin{eqnarray} \phi _{\textit{ab}}(f) = \arctan (\mathcal{Q}_{\textit{ab}}(f)/\mathcal{C}_{\textit{ab}}(f)), \end{eqnarray}
where
$\mathcal{S}_{\{a,b\}}(f)$
represents the energy spectra of
$a,b$
. Physically, the phase spectrum
$\phi _{\textit{ab}}(f)$
measures the synchronisation between two time series, while the coherence spectrum
$\gamma _{\textit{ab}}(f)$
can be seen as the correlation coefficient in the frequency domain. If two time series tend to co-occur at a given frequency, for example, when the minima of
$a(t)$
coincide with the minima of
$b(t)$
, then
$\phi _{\textit{ab}}(f)$
would be close to zero. Otherwise, if the minima of
$a(t)$
frequently coincide with the maxima of
$b(t)$
, then
$\phi _{\textit{ab}}(f)$
would be close to
$\pm \pi$
. Besides, low values of the coherence spectrum
$\gamma _{\textit{ab}}(f)$
indicate that the underlying physical processes behind
$a(t)$
and
$b(t)$
may differ.
2.2.2. Spectral proper orthogonal decomposition
Spectral proper orthogonal decomposition, essentially a POD method in the frequency domain, aligns with our objective of identifying the flow modes associated with the dominant oscillation frequencies. In the following, we will briefly outline its practical applications without delving into mathematical intricacies and its connection with other analyses, which can be readily found in previous studies (Lumley Reference Lumley2007; Towne et al. Reference Towne, Schmidt and Colonius2018; Abreu et al. Reference Abreu, Cavalieri, Schlatter, Vinuesa and Henningson2020). Given a matrix
$\boldsymbol{X}=[\boldsymbol{x}_1,{\cdots} , \boldsymbol{x}_i,{\cdots} ,\boldsymbol{x}_{\!N}]$
with each column
$\boldsymbol{x}_i$
representing a PIV snapshot, we can divide it into
$N_b$
blocks
$\boldsymbol{X} = [\boldsymbol{Y}_1,{\cdots} ,\boldsymbol{Y}_{\!j},{\cdots} ,\boldsymbol{Y}_{N_b}]$
, where each block contains
$N_{\!f}$
snapshots, i.e.
$\boldsymbol{Y}_{\!j} = [\boldsymbol{x}_{1}^{j},{\cdots} ,\boldsymbol{x}_{q}^{j},{\cdots} ,\boldsymbol{x}_{N_{\!f}}^{j}]$
. Here,
$\boldsymbol{x}_{q}^{j}=\boldsymbol{x}_{q+(j-1)p}$
with
$p$
denoting the separation between adjacent blocks. After performing the discrete Fourier transform for each block, we obtain
$\boldsymbol{\hat {Y}}_{\!j} = [\boldsymbol{\hat {x}}_{1}^{j},{\cdots} ,\boldsymbol{\hat {x}}_{k}^{j},{\cdots} ,\boldsymbol{\hat {x}}_{N_{\!f}}^{j}]$
with each entry
$\boldsymbol{\hat {x}}_{k}^{j}$
as follows:
\begin{gather} \boldsymbol{\hat {x}}_{k}^{j} = \sum _{q=1}^{N_{\!f}}\boldsymbol{x}_{q+(j-1)p}e^{-i2\pi (q-1) f_k\Delta t}, \end{gather}
\begin{gather} f_k = \begin{cases} \frac {k-1}{N_{\!f}\Delta t} \quad \text{for} \quad k\leqslant N_{\!f}/2,\\[11pt] \frac {k-1-N_{\!f}}{N_{\!f}\Delta t} \quad \text{for} \quad k\geqslant N_{\!f}/2, \end{cases} \end{gather}
where
$\Delta t = 1/f$
is the sampling interval. Gathering the Fourier component at a given frequency
$f_k$
from each block, we obtain
$\boldsymbol{\hat {X}}(f_k) = [\boldsymbol{\hat {x}}_{k}^{(1)},{\cdots} ,\boldsymbol{\hat {x}}_{k}^{j},{\cdots} , \boldsymbol{\hat {x}}_{k}^{(N_b)}],{}j=1,2,{\cdots} ,N_b$
. Similar to the snapshot POD method (Sirovich Reference Sirovich1987; Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993; Holmes Reference Holmes2012), we perform the singular value decomposition for
$\boldsymbol{\hat {X}}(f_k)$
as follows:
\begin{align} \boldsymbol{\hat {X}}(f_k) = \boldsymbol{\varPhi }(f_k)\boldsymbol{\varSigma }(f_k)\boldsymbol{\varPsi }(f_k)^*, \end{align}
where
$\boldsymbol{\varSigma }(f_k) = diag[\sqrt {\lambda _1},\sqrt {\lambda _2},{\cdots} ,\sqrt {\lambda _{N_b}}]$
and
$\boldsymbol{\varPhi }(f_k)=[\boldsymbol{\phi }_1(f_k),\boldsymbol{\phi }_2(f_k),{\cdots} , \boldsymbol{\phi }_{N_b}{}(f_k)]$
. Each column vector
$\boldsymbol{\phi }_i(f_k)$
represents the SPOD mode ranked by its energy
$\lambda _i(f_k), i = 1, 2, {\cdots} , N_b$
.
For the experimental PIV data analysed in the main text, the state vector contains only the two in-plane velocity components sampled on a uniform grid, so the SPOD inner product reduces to the standard Euclidean one. For the direct numerical simulation (DNS) data discussed in § 4.1, the state vector
$(u,w,\theta )$
contains both velocity and temperature in the vertical mid-plane of the convection cell, and the SPOD inner product uses the local cell-volume weighting to maintain a physically meaningful and grid-independent modal energy; the relative temperature weighting is taken as unity, which is not a critical issue for the present purpose (Podvin & Sergent Reference Podvin and Sergent2012).
2.2.3. Resolvent analysis
For turbulent flows under statistically steady conditions, the governing equations can be linearised through a standard Reynolds decomposition, yielding
\begin{equation} \frac {\partial \boldsymbol{u}^\prime }{\partial t} + \mathcal{L}\boldsymbol{u}^\prime = \boldsymbol{q}. \end{equation}
Here,
$\mathcal{L}$
is the linearised flow operator,
$\boldsymbol{q}$
contains the remaining nonlinear terms and
$\boldsymbol{u}^\prime = \boldsymbol{u}-\boldsymbol{U}$
is the fluctuation velocity. Because the exact 3-D mean field is inaccessible to the present planar PIV measurements, we do not attempt to reconstruct the full physics-based resolvent operator. Instead, we adopt the DD procedure of Herrmann et al. (Reference Herrmann, Baddoo, Semaan, Brunton and McKeon2021) to infer an effective linear operator from the PIV snapshots. In the present paper, this operator is interpreted pragmatically as a reduced input–output model for the dominant measured dynamics, rather than as a complete surrogate for the linearised Boussinesq operator. Briefly, the method uses dynamic mode decomposition to approximate the linear evolution operator as
\begin{equation} \boldsymbol{u}^\prime _{k+1} = \exp (\boldsymbol{A}\Delta t)\boldsymbol{u}^\prime _{k}, \end{equation}
where
$\boldsymbol{u}^\prime _{k+1}$
is the fluctuation quantity at time
$t_k=k\Delta t$
, and
$\Delta t$
is the sampling time. By substituting
$\boldsymbol{A}$
into (2.9) and taking the temporal Fourier transform for both sides, (2.9) can be written in the frequency domain as
\begin{equation} \tilde {\boldsymbol{u}^\prime } = \boldsymbol{H}(\omega )\tilde {\boldsymbol{q}}, \end{equation}
where
$\boldsymbol{H}(\omega )=(i\omega \boldsymbol{I} + \boldsymbol{A})^{-1}$
is the resolvent operator. It serves as a transfer function that amplifies (or attenuates) the harmonic forcing input
$\tilde {\boldsymbol{q}}$
and maps it to the response
$\tilde {\boldsymbol{u}^\prime }$
. The basic object of resolvent analysis is to find the dominant direction along which the largest input–output gains can be obtained, i.e.
\begin{equation} \sigma _1^2(\omega ) = \textit{max}\frac {||\tilde {\boldsymbol{u}^\prime }||^2}{||\tilde {\boldsymbol{q}}||^2} = ||\boldsymbol{H}(\omega )||^2. \end{equation}
For the experimental data considered here, the norm
$||\boldsymbol{\cdot }||$
is the Euclidean norm associated with the two in-plane velocity components on the uniform PIV grid. Practically, it is accomplished by performing the singular value decomposition of the resolvent operator
\begin{equation} \boldsymbol{H}(\omega ) = \boldsymbol{U}(\omega )\boldsymbol{S}(\omega )\boldsymbol{V}^*(\omega ), \end{equation}
where
$\boldsymbol{S}$
is a diagonal matrix containing the gains
$\sigma _1\geqslant \sigma _2\geqslant {\cdots} \geqslant \sigma _n\geqslant 0$
, while
$\boldsymbol{U}(\omega )$
and
$\boldsymbol{V}(\omega )$
are unitary matrices whose columns are the resolvent forcing modes and resolvent response modes, respectively.
3. Results
3.1. Global angular momentum
As mentioned above, previous studies of the LSC oscillation in turbulent RBC frequently utilised local measurements of temperature or velocity, which are susceptible to the influence of local flow structures and might produce inconsistent results across different regions. To eliminate the limitations of localised measurements, here, we analyse the global angular momentum
$L(t)$
defined as
\begin{equation} L(t)=\frac {1}{\varOmega }\int _{\mathcal{D}} \boldsymbol{e}_y \boldsymbol{\cdot }(\boldsymbol{x}\times \boldsymbol{U}(\boldsymbol{x},t))\, \text{d}\varOmega , \end{equation}
which can serve as a more holistic indicator of the large-scale flow dynamics. Figure 1(a) displays the evolution of the global angular momentum
$L(t)$
for different
$Ra$
numbers. These time series indicate that the LSC maintained a consistent direction and no flow reversal occurred during the PIV experiments. The auto-correlation function
$\rho _L(\tau )$
of
$L(t)$
, as shown in figure 1(b), reveals well-defined oscillation signals that persist over multiple periods, suggesting that the oscillation in the flow is highly coherent. The oscillation period
$\tau _o$
can be determined by the first peak
$\tau _{peak}$
of
$\rho _L(\tau )$
, as indicated by the coloured arrows in figure 1(b). The values of
$\tau _{peak}$
plotted in figure 1(c) closely match the turnover time (
$\pi H/U_{\textit{max}}$
) of the LSC, and both time scales follow a scaling relation
$\tau _o \sim Ra^{-0.45}$
, consistent with the findings in previous studies (Qiu & Tong Reference Qiu and Tong2001; Xia et al. Reference Xia, Sun and Zhou2003; Sun & Xia Reference Sun and Xia2005; Niemela & Sreenivasan Reference Niemela and Sreenivasan2006; Vogt et al. Reference Vogt, Horn, Grannan and Aurnou2018; Li et al. Reference Li, Chen, Xu and Xi2022; Cheng et al. Reference Cheng, Mohammad, Wang, Keogh, Forer and Kelley2022). The correspondence between the oscillation period and the circulation time of the LSC supports the physical picture that the periodic oscillations are associated with the advection of the LSC (Sano et al. Reference Sano, Wu and Libchaber1989; Cioni, Ciliberto & Sommeria Reference Cioni, Ciliberto and Sommeria1997; Qiu & Tong Reference Qiu and Tong2001; Brown & Ahlers Reference Brown and Ahlers2009).
(a) Time series of the normalised global angular momentum
$L/\sigma _L$
. From bottom to top,
$Ra\in \{2.5\times 10^9, 3.9\times 10^9, 6.2\times 10^9, 9.8\times 10^9, 1.2\times 10^{10}\}$
. The time series is vertically shifted for clarity. (b) Normalised auto-correlation function
$\rho _L(\tau )$
of
$L (t)$
, with arrows indicating the peak time
$\tau _{peak}$
for each
$Ra$
. (c) Comparison between the oscillation period
$\tau _o$
estimated from the peak time
$\tau _{peak}$
, i.e.
$\tau _o = \tau _{peak}$
, and the turnover time of the LSC
$\pi H/U_{\textit{max}}$
. Here,
$U_{\textit{max}}$
represents the maximum mean velocity of the LSC. It is seen that
$\tau _{peak}$
and
$\pi H/U_{\textit{max}}$
are closely aligned and both follow a power-law relation:
$\tau _o\sim Ra^{-0.45}$
. (d) Normalised power spectral densities
$\mathcal{S}(f)/\mathcal{S}(f_o)$
as a function of
$f/f_o$
. The two dashed lines denote two peak frequencies at
$f_1 = f_o = 1/\tau _o$
and
$f_2 = 2f_o$
.
By taking the oscillation frequency as
$f_o = 1/\tau _o = 1/\tau _{peak}$
, we can further obtain the normalised power spectral density
$\mathcal{S}(f)/\mathcal{S}(f_o)$
as shown in figure 1(d). In the low-frequency range, there are two peak frequencies at
$f_1\approx f_o$
and
$f_2\approx 2f_o$
. The primary peak
$f_1\approx f_o$
is pretty strong, while its harmonic frequency
$f_2\approx 2f_o$
acts as a weak bump due to turbulent noise. Similar harmonic frequencies have also been observed in other studies of temperature signals (see, e.g. Takeshita et al. (Reference Takeshita, Segawa, Glazier and Sano1996), Ciliberto, Cioni & Laroche (Reference Ciliberto, Cioni and Laroche1996), Qiu & Tong (Reference Qiu and Tong2001), Song et al. (Reference Song, Villermaux and Tong2011)), whose origin was at first attributed to the alternating emission of thermal plumes between the upper and lower boundary layers, but latter recognised to result from the horizontal sloshing motion of the LSC in a cylindrical convection cell (Xi et al. Reference Xi, Zhou, Zhou, Chan and Xia2009; Zhou et al. Reference Zhou, Xi, Zhou, Sun and Xia2009; Brown & Ahlers Reference Brown and Ahlers2009). However, similar to the study by Song et al. (Reference Song, Villermaux and Tong2011), the existence of horizontal sloshing modes and their contribution to the harmonic frequency within the present quasi-2-D convection cell are unlikely. This raises the possibility that an alternative physical mechanism may be at play in generating this harmonic frequency. Since
$\mathcal{S}(f)/\mathcal{S}(f_o)$
shown in figure 1(d) share a remarkable consistency for different
$Ra$
numbers, which suggests that the underlying flow dynamics is governed by a universal mechanism, our subsequent analysis will focus on the case at
$Ra=1.2\times 10^{10}$
, without loss of generality.
3.2. Cross-correlation and cross-spectrum
In this subsection, we employ two different probing configurations to examine the flow oscillation characteristics. Figure 2(a) illustrates the conventional configuration, where we select symmetric points along the horizontal and vertical mid-lines that traverse the flow domain. These points are strategically chosen because, as depicted in the contour plot, they are located on the circulation path of the LSC. Consequently, the cross-correlation analysis at these points can effectively reflect the coherence of the LSC. Figure 2(b) introduces an alternative approach to examine the temporal correlation in the flow. This approach focuses on the flow regions with significant turbulent fluctuations, which are roughly demarcated into the lower-left and upper-right areas by the black lines with appropriate thresholds of velocity fluctuations. According to previous findings in similar convection cells (Xia et al. Reference Xia, Sun and Zhou2003; Huang et al. Reference Huang, Kaczorowski, Ni and Xia2013), these regions correspond to areas where rising and falling thermal plumes are abundant. Therefore, analysing the data in these regions can help to understand the role of thermal plumes in the LSC oscillation.
Two different probing configurations for examining the oscillation characteristics in the flow. (a) The conventional configuration is based on local point probes located symmetrically at the horizontal and vertical mid-lines, as illustrated by the red pentagrams and the blue triangles. (b) An alternative approach is based on the flow regions with energetic turbulent fluctuations, which are demarcated by the black lines with a threshold of
$1.45$
. The contours in panels (a) and (b) are the magnitude of the mean velocity
$\sqrt {\boldsymbol{U}^2}$
and the root-mean-square (r.m.s.) of the velocity fluctuation
$\sqrt {\langle \boldsymbol{u}^{\prime 2}\rangle _t}$
with
$\boldsymbol{u}^\prime = \boldsymbol{u} -\boldsymbol{U}$
, respectively. Both
$\sqrt {\boldsymbol{U}^2}$
and
$\sqrt {\langle \boldsymbol{u}^{\prime 2}\rangle _t}$
are normalised by their respective spatial r.m.s. values.
Contour plots of the cross-correlation between the symmetric points illustrated in figure 2: (a,b)
$\rho _{u_i}(\tau /\tau _o;-x,x)$
for probes on the horizontal mid-line; (c,d)
$\rho _{u_i}(\tau /\tau _o;-z,z)$
for probes on the vertical mid-line. Panels (a,c) correspond to the horizontal velocity component
$u_i = u$
, while panels (b,d) correspond to the vertical velocity component
$u_i = w$
. The black solid line in each panel represents the spatially averaged cross-correlation value, with black dashed lines marking the negative peaks around
$\tau /\tau _o = \pm 0.5$
. The blue vertical dashed lines in panel (a) indicate the positive peaks around
$\tau /\tau _o = 1/3$
.
Figure 3 examines the oscillation behaviour based on the conventional probing configuration. Specifically, figures 3(a) and 3(b) display the contour plots of the cross-correlation
$\rho _{u_i}(\tau ;-x,x)$
for probes on the horizontal mid-line, while figures 3(c) and 3(d) show
$\rho _{u_i}(\tau ;-z,z)$
for probes on the vertical mid-line. Here,
$u_i$
represents the horizontal velocity component
$u$
in (a, c) and the vertical velocity component
$w$
in (b, d). Although the detailed correlation patterns vary strongly with probe position, all four panels exhibit clear negative troughs around
$\tau /\tau _o=\pm 1/2$
. This is even more evident after spatial averaging, as shown by the black solid lines. Additionally, figure 3(a) shows positive peaks at
$\tau /\tau _o\approx \pm 1/3$
, suggesting a robust contribution from a component with
$f\tau _o\approx 3$
. In the travelling-wave picture revealed in § 3.3, the sign of the cross-correlation is determined by the phase relation between the two probe signals at a given lag: positive values correspond to an approximately in-phase alignment, whereas negative values correspond to an approximately anti-phase alignment. The negative extrema near
$\tau /\tau _o=\pm 1/2$
therefore indicate that a half-period shift tends to align the signals in anti-phase, consistent with a coherent propagating pattern.
Cross-correlation (a) and cross-spectra (b–d) between the angular momenta calculated from the lower-left region
$L_l(t) = \langle L\rangle _{\boldsymbol{x}\in \mathcal{A}_l}$
and the upper-top region
$L_u(t) = \langle L\rangle _{\boldsymbol{x}\in \mathcal{A}_u}$
, as illustrated in figure 2(b). Specifically, panels (b–d) are the amplitude spectrum
$|\mathcal{S}_{ul}(f)|$
, the coherence spectrum
$\gamma _{ul}(f)$
, and the absolute value of the phase spectrum
$|\phi _{ul}(f)|$
. The three vertical dashed lines correspond to normalised frequencies
$f\tau _o$
of
$1,2$
and
$3$
, respectively.
Figure 4 shows the results under the second probing configuration. Rather than using local velocity components, this configuration analyses the spatially averaged angular momenta of two distinct regions illustrated in figure 2(b). The cross-correlation
$\rho _{ul}(\tau )$
exhibits positive peaks around
$\tau /\tau _o\approx 1/3$
and
$1/2$
, together with a negative trough around
$\tau /\tau _o\approx 1$
. Again, these alternating signs can be naturally interpreted as a consequence of the phase shift accumulated by a propagating structure between the two regions, not as evidence for different physical mechanisms at different lags. Fourier transforming
$\rho _{ul}(\tau )$
yields the cross-spectrum
$\mathcal{S}_{ul}(f)$
, from which three dominant frequencies,
$f\tau _o\approx 1$
,
$2$
and
$3$
, can be identified in the amplitude spectrum
$|\mathcal{S}_{ul}(f)|$
shown in figure 4(b). The coherence spectrum
$\gamma _{ul}(f)$
in figure 4(c) confirms that these frequencies correspond to strongly correlated large-scale motions. The phase spectrum in figure 4(d) further shows that the diagonal regions are approximately out of phase at
$f\tau _o\approx 1$
and closer to in phase at
$f\tau _o\approx 2$
and
$3$
.
Taken together, the correlation analysis above demonstrates that the oscillation dynamics is coherent over the cell, which contains a fundamental low-frequency component and higher harmonics. Note that the appearance of harmonics is not necessary to be interpreted as several independent intrinsic time scales. Rather, it may imply that the flow oscillation is not purely sinusoidal and therefore projects energy onto integer multiples of the fundamental frequency, as expected for a travelling-wave-like dynamics. This motivates us to carry out a field-based modal analysis to extract the corresponding spatial structures.
3.3. Travelling-wave-like flow structures identified by SPOD
To reveal the underlying flow processes responsible for the observed oscillations, we employed the SPOD method outlined in § 2.2.2 to analyse the velocity field in this subsection. This method is particularly well suited for extracting intrinsic patterns from a large dataset of flow fields. The SPOD modes were estimated by using blocks that contain
$N_{\!f} = 1024$
snapshots, leading to
$N_b=51$
realisations. A standard Hann window was used to reduce the spectral leakage.
(a) Energy spectra
$\lambda _n(f)$
of the leading four SPOD modes. The black dashed lines denote four dominant oscillation frequencies
$f\tau _o=1,2,3$
and
$4$
, whose corresponding first SPOD modes are respectively shown in (b–e). The vortex structures are characterised by streamlines, with the contour indicating the vorticity distribution (scaled by the r.m.s. value of the entire field).
Figure 5(a) presents the SPOD eigenspectra
$\lambda _n(f)$
as a function of the dimensionless frequency
$f\tau _o$
. To facilitate interpretation, all eigenvalues are normalised by the total flow energy integrated over all frequencies, allowing each eigenvalue to represent the contribution of its respective mode. For clarity, only the leading four eigenvalues are displayed. The significant gap between the first and second eigenvalues highlights the low-rank character of the dynamics. Notably, clear spectral peaks are concentrated in the leading SPOD mode
$\lambda _1(f)$
. This does not imply that the higher modes should peak successively at
$f\tau _o=2,3,4$
, as SPOD modes are ranked independently at each frequency; when one coherent structure dominates a given frequency, most of the organised energy remains captured by the first mode, and the subleading modes mainly contain weaker or less coherent fluctuations. The first SPOD modes at the dominant frequencies are shown in figure 5(b–e). These modes manifest as counter-rotating vortex pairs, with the number of pairs increasing together with the oscillation frequency. Specifically, the flow domain shows one, two, three and four pairs of vortices for frequencies
$f\tau _o=1, 2, 3$
and
$4$
, respectively. A closer examination reveals that the vortices are approximately azimuthally distributed around the periphery of the convection cell. Here, ‘azimuthal’ means that the flow field is viewed in a polar coordinate system centred at the cell centre.
This observation offers a significant insight: even though the cell’s geometry is not azimuthally periodic, the flow modes can exhibit a form of azimuthal periodicity to a certain extent. Thus, we can assign an azimuthal wavenumber
$k_\varphi$
to each vortex pair, corresponding to the number of pairs:
$1, 2, 3$
and
$4$
. The azimuthally distributed patterns suggest that there is a linear relation between the azimuthal wavenumber and the dominant oscillation frequencies, i.e.
$k_\varphi \sim f\tau _o$
(see Appendix A for more discussion). Given that
$\tau _o \approx \pi H/U_{\textit{max}}$
(see figure 1
c), we can obtain a linear dispersion relation between the angular frequency
$\omega =2\pi f$
and the azimuthal wavenumber
$k_\varphi$
, i.e.
$\omega \sim {2U_{\textit{max}}}/{H}k_\varphi$
. This dispersion relation implies that the observed periodic oscillations can be regarded as a linear wave propagating around the convection cell in polar coordinates
$(r, \varphi )$
\begin{equation} \boldsymbol{u}^\prime (r,\varphi ,t) \approx \boldsymbol{u}_0^\prime (r)e^{i(k_\varphi \varphi -\omega t)}. \end{equation}
Equation (3.2) closely mirrors the advected-oscillation model proposed by Brown & Ahlers (Reference Brown and Ahlers2009). In their model, the LSC oscillations observed in cylindrical convection cells can be described by a pair of oppositely travelling waves, which obey the empirical equation
\begin{equation} \alpha _h = Ae^{i(nk_0z - \omega _n t)},\quad \alpha _c = Ae^{i(nk_0z - \omega _n t) + i\psi }. \end{equation}
Here,
$\alpha _h$
and
$\alpha _c$
represent the sloshing displacement angles of the hot and cold plumes, respectively;
$\psi = (n+1)\pi$
represents the relative phase difference between
$\alpha _h$
and
$\alpha _c$
; and
$k_0 = \pi /(2L)$
. Their analysis reveals that the integer value of
$n$
must be even order for the model to describe the experimentally observed torsion and sloshing modes. Besides, due to the damping effect, the lowest-order even mode tends to dominate, i.e.
$n=2$
.
Comparing (3.2) and (3.3) leads to several insights. Firstly, their consistent forms suggest that advection-induced waves are crucial for the observed periodic oscillations, regardless of the cell’s geometry. This implies that the concept of travelling waves is universally applicable to both quasi-2-D and cylindrical convection cells. Secondly, a key distinction between the Brown–Ahlers model and our case is the direction of the wavenumber. In Brown–Ahlers model, it is constrained to the
$z$
-axis, whereas in our scenario, it is oriented along the azimuthal direction, which aligns more naturally with the flow patterns as manifested by the first SPOD modes. As we will discuss in § 4, our scenario can better account for the circulation properties of thermal plumes, making it suitable for different geometric configurations.
3.4. Linear-amplification mechanism revealed by resolvent analysis
Upon identifying the travelling-wave-like oscillation structures, a natural question arises: Can these dominant structures be interpreted within an input–output framework? This motivates us to perform a DD resolvent analysis in this subsection. We stress that the purpose of this analysis is not to provide spectral content independent of SPOD, but to examine whether the same measured structures align with the most amplified responses of an effective linear dynamics.
Figure 6(a) displays the gain spectra
$\sigma _n^2(f)$
as a function of the dimensionless frequency
$f\tau _o$
. Again, we focus on the four most dominant resolvent modes for clarity. Similar to the SPOD eigenspectra in figure 5(a), there is a significant gap between the first and second gains, and the leading gain
$\sigma _1^2(f)$
exhibits four distinct peaks that closely align with the peaks of the principal SPOD mode
$\lambda _1(f)$
. The corresponding first resolvent response modes are shown in figure 6(b–e). Their qualitative similarity to the SPOD modes indicates that the measured travelling-wave-like structures are also the dominant responses of the inferred input–output operator. As noted in previous studies (Towne et al. Reference Towne, Colonius, Jordan, Cavalieri and Bres2015, Reference Towne, Bres and Lele2016, Reference Towne, Schmidt and Colonius2018; Semeraro et al. Reference Semeraro, Jaunet, Jordan, Cavalieri and Lesshafft2016), resolvent and SPOD modes become identical under white-noise forcing. Turbulent RBC is not expected to satisfy that ideal limit exactly, so we do not interpret the similarity as proof that the two analyses contain independent information. Instead, we take it as evidence that the dominant oscillatory structures identified statistically by SPOD are consistent with an amplification picture based on an effective linear response. More broadly, the agreement may also partly reflect the shared connections of DD-resolvent and SPOD with Koopman-based descriptions (Rowley et al. Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009; Mezić Reference Mezić2013; Budišić et al. Reference Budišić, Mohr and Mezić2012). This suggests that some degree of similarity between the two approaches may also be expected on theoretical grounds.
(a) Energy spectra
$\sigma ^2_n(f)$
of the leading four resolvent response modes. The black dashed lines denote four dominant oscillation frequencies
$f\tau _o=1,2,3$
and
$4$
, whose corresponding first resolvent response modes are respectively shown in (b–e). The large-scale vortex structures are characterised by streamlines, with the contour indicating the vorticity distribution (scaled by the r.m.s. value of the entire field).
Input–output description of the LSC oscillation mechanism.
To put this mechanism more clearly, let us write down the linearised Navier–stokes equation more explicitly and perform temporally Fourier transformation for both sides. The resulting governing equation is given by
\begin{align} \begin{bmatrix} \hat {\boldsymbol{u}} \\[3pt] \hat {\theta } \\[3pt] \hat {p} \end{bmatrix} = \underbrace { \left ( -i\omega \begin{bmatrix} \boldsymbol{I} & & \\ & 1 & \\ & & 0 \end{bmatrix} + \begin{bmatrix} \mathcal{A} & 1 & \boldsymbol{\nabla }\\[3pt] \mathcal{C} & \mathcal{B} & 0 \\[3pt] \boldsymbol{\nabla }& 0 & 0 \end{bmatrix} \right )^{-1} }_{\mathcal{H}} \begin{bmatrix} \hat {\boldsymbol{f}} \\[3pt] \hat {g} \\[3pt] 0 \end{bmatrix}\!, \end{align}
where
$\mathcal{A},\mathcal{B}$
and
$\mathcal{C}$
are linear operators, such that
\begin{align} \mathcal{A} = ()\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{U} &+ \boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla }() - 1/\sqrt {Ra/Pr}\boldsymbol{\nabla} ^2(), \nonumber \\ \mathcal{B} = \boldsymbol{U}\boldsymbol{\cdot }&\boldsymbol{\nabla }() - 1/\sqrt {\textit{RaPr}}\boldsymbol{\nabla} ^2(), \nonumber \\ &\mathcal{C} = ()\boldsymbol{\cdot }\boldsymbol{\nabla }T, \end{align}
while
$\boldsymbol{f} = -(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}+\overline {(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}}$
and
$g = -(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\theta +\overline {(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\theta }$
represent the nonlinear interactions within the velocity fluctuations and between the velocity and temperature fluctuations, respectively. If we regard these nonlinear terms as nonlinear forcing, then (3.4) can be interpreted as an input–output process, as depicted in figure 7.
Specifically, the nonlinear terms
$\boldsymbol{f}$
and
$g$
can be interpreted as effective forcing that feeds fluctuation energy back into the linear dynamics. The resolvent operator
$\mathcal{H}(\omega )$
then determines which forcing directions are preferentially amplified. Although the full operator cannot be recovered from the present data, a simplified interpretation can still be obtained from the observed bulk flow. In the central region occupied by the persistent LSC, the mean azimuthal advection is the dominant contribution, whereas the azimuthal variations of the mean quantities are comparatively weak. We therefore use
$\boldsymbol{\nabla }\boldsymbol{U}\approx 0$
and
$\boldsymbol{\nabla }T\approx 0$
only as bulk approximations, not as exact statements valid throughout the whole cell; likewise, the diffusive terms are subdominant at the high
$Ra$
considered here. This approximation is consistent with the well-mixed character of the turbulent bulk, where the strongest mean temperature gradients remain concentrated in the thermal boundary layers near the plates (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009). Under these approximations, the azimuthal Fourier transform of (3.4) yields
\begin{equation} \mathcal{H}(\omega ) \approx \begin{bmatrix} i(2k_\varphi U_\varphi /H-\omega )\boldsymbol{I} & 1 \\[5pt] 0 & i(2k_\varphi U_\varphi /H -\omega ) \end{bmatrix}^{-1}\!. \end{equation}
Equation (3.6) shows that when the frequency
$\omega$
is close to the advection frequency set by the azimuthal wavenumber
$k_\varphi$
, namely
$\omega \approx {2U_{\varphi }}/{H}k_\varphi$
, the resolvent norm becomes large. In the simplified operator, this amplification arises because the input–output map approaches a resonance-like singular response of an advective system; it should not be confused with modal instability. In practice, viscous effects neglected here, finite bandwidth and nonlinear forcing regularise the response, so the measured gain remains finite and appears as broadened spectral peaks rather than mathematical singularities. Thus, the predominantly imaginary spectrum reflects the advection-dominated dynamics, while the large resolvent gain indicates selective amplification of fluctuations whose frequency matches the mean-circulation advection rate (McKeon & Sharma Reference McKeon and Sharma2010; McKeon Reference McKeon2017).
This result is consistent with the linear dispersion relation
$\omega \sim {2U_{\textit{max}}}/{H}k_\varphi$
inferred earlier from the SPOD analysis (see § 3.3). In other words, the observed frequency–wavenumber relation is the natural signature of fluctuations being preferentially amplified when they are convected at the mean azimuthal speed of the LSC. Since the admissible azimuthal wavenumbers are discrete, the preferred frequencies are also discrete. We therefore interpret the oscillation as an advection-selected, linearly amplified response of a persistent mean circulation, rather than as a self-excited instability of the mean flow.
4. Discussion
Through SPOD analysis, we have identified travelling-wave-like structures underlying the oscillation dynamics. The complementary resolvent analysis suggests that these structures can be interpreted within a linear-amplification framework associated with the persistent mean circulation. Within this interpretive framework, we consider two corresponding corollaries in this section.
-
(i) If velocity and temperature fluctuations are governed by the same simplified amplification relation in (3.6), one expects their dominant oscillations to occur concurrently.
-
(ii) If this input–output interpretation captures the leading dynamics of the LSC rather than a geometry-specific detail, it may be also applicable to explain the oscillation behaviours observed in RBC with other geometries, such as cylindrical ones.
4.1. Co-occurrence of the temperature–velocity oscillations
To check the first corollary, we performed a complementary numerical simulation in a quasi-2-D domain with
$\varGamma _\perp =1$
and
$\varGamma _\parallel =1/3$
. The domain has no-slip boundaries and adiabatic sidewalls. The simulation was conducted at
$Ra=1\times 10^9$
and
$\textit{Pr}=4.3$
by employing a parallel DNS method under a second-order finite-difference scheme on staggered grids (Bao, Luo & Ye Reference Bao, Luo and Ye2018; Gao et al. Reference Gao, Luo and Bao2018, Reference Gao, Tao, Huang, Bao and Xie2024). To analyse the temperature and velocity fields jointly, we use the state vector
$(u,w,\theta )$
in the vertical mid-plane slice for the SPOD analysis. Following Podvin & Sergent (Reference Podvin and Sergent2012), the rescaling factor for the temperature relative to the velocity is set to
$1$
. Accordingly, the SPOD inner product for two DNS snapshots
$\boldsymbol{q}_1=(u_1,w_1,\theta _1)$
and
$\boldsymbol{q}_2=(u_2,w_2,\theta _2)$
is defined as
$\langle \boldsymbol{q}_1,\boldsymbol{q}_2 \rangle = \sum _{\!j} (u_{1,j}u_{2,j}+w_{1,j}w_{2,j}+\theta _{1,j}\theta _{2,j})\Delta V_{\!j}$
, where
$\Delta V_{\!j}$
is the local cell volume of the
$j$
th grid point in the slice; the associated norm is
$||\boldsymbol{q}||=\langle \boldsymbol{q},\boldsymbol{q}\rangle ^{1/2}$
. Thus, the essential weighting here is the cell-volume factor required for grid independence, whereas the unit temperature rescaling is kept for convenience and is not a critical parameter for the present analysis (Podvin & Sergent Reference Podvin and Sergent2012). A total of 5000 slices (corresponding to several turnover times of the LSC) were used for the analysis. These data were processed in blocks of
$N_{\!f} = 1024$
snapshots with
$50\,\%$
overlap, resulting in
$N_b = 8$
realisations.
The SPOD analysis results of the numerical simulation in a quasi-2-D RBC at
$Ra=1\times 10^9$
and
$\textit{Pr}=4.3$
; (a) SPOD energy spectra
$\lambda _n(f)$
. The black dashed lines denote three dominant oscillation frequencies
$f\tau _o=1,2$
and
$f\tau _o \approx 3$
, whose corresponding first SPOD modes are respectively shown in (b–d). The large-scale vortex structures are characterised by streamlines, with the contour indicating the temperature distribution.
Figure 8 presents the SPOD eigenspectra and the corresponding first SPOD modes. It is seen that the most distinct low-rank dynamics occurs at dominant frequencies around
$f\tau _o \approx 1, 2$
. The corresponding flow structures manifest as one-pair and two-pair counter-rotating vortices, each with azimuthal wavenumbers of 1 and 2, respectively. At higher frequencies, such as
$f\tau _o \gt 3$
, this dynamics is less distinct, and the flow structures shown in figure 8(d) exhibit less organisation. This reduced clarity may be attributed to the limited numerical dataset available for the SPOD analysis, as only
$N_b=8$
realisations might not provide sufficient data for the statistical convergence of high-frequency, low-energy modes. Despite this limitation, the azimuthal distribution of these flow structures still generally follows the travelling-wave pattern, in line with our experimental observations in the previous section. Besides, a strong correlation is observed when comparing the streamline patterns with the temperature distributions. The paths of the rising hot plume and the sinking cold plume closely align with the flow directions shown by the streamlines. This close correlation indicates that the temperature and velocity oscillations are intrinsically linked in the system, maintaining the coherence of flow structures.
To validate this synchronous behaviour, we have also performed SPOD analysis for the temperature and velocity fields separately. For the sake of conciseness, these additional results are displayed in the Appendix B. Although their eigenspectra show some slight variations, the peak frequencies and oscillation patterns are consistent with those observed in figure 8(b–d). Based on the resemblance between the velocity and temperature SPOD modes, it is evident that the temperature and velocity oscillations are interconnected and co-occur in the flow. This may explain why the previous studies found that the temperature oscillation frequencies and the velocity oscillation frequencies coincide with each other, e.g. Qiu & Tong (Reference Qiu and Tong2001, Reference Qiu and Tong2002); Zhou et al. (Reference Zhou, Sun and Xia2007), Brown & Ahlers (Reference Brown and Ahlers2009) and Xi et al. (Reference Xi, Zhou, Zhou, Chan and Xia2009).
4.2. Applicability to different geometric configurations
To explore whether the linear-amplification mechanism is applicable to other geometric configurations, we now examine the velocity fields measured in a cylindrical RBC with
$Ra=5.2\times 10^9$
and
$\textit{Pr} = 5.7$
. Details of the experiment have been introduced in Li et al. (Reference Li, Chen, Xu and Xi2022). Here, we only highlight the essential feature. The convection cell’s dimensions are 19 cm in diameter and 18 cm in height, thus its aspect ratio is around one. In order to focus on the oscillation dynamics, the convection cell was tilted at an angle of
$2^\circ$
to prevent the diffusive meandering of the LSC. The velocity fields were measured by a standard PIV technique and a total of 10 820 velocity maps were collected at a frequency of 1
$Hz$
. During the SPOD analysis, the dataset was segmented into
$N_b = 169$
blocks, each containing
$N_{\!f} = 128$
snapshots with
$50\,\%$
overlap.
The SPOD results are shown in figure 9(a–c). It is seen that the flow structures at the dominant frequencies
$f\tau _o \approx 1$
and
$2$
exhibit travelling-wave-like patterns with azimuthal wavenumbers
$k_\varphi = 1$
and
$2$
, respectively. For the SPOD mode at the primary frequency
$f\tau _o\approx 1$
(figure 9
b), the sign difference of the vorticity distribution suggests that the flow structures in the lower and upper parts of the convection cell move out of phase. This oscillation mode results in a torsional motion when superimposed on the mean velocity field. Similarly, the superposition of the mean flow and the SPOD mode at the harmonic frequency
$f\tau _o\approx 2$
, characterised by two-pair vortices, as manifested in figure 9(c), leads to a sloshing motion of the bulk flow. These flow patterns are consistent with the oscillation behaviours of the LSC in cylindrical RBC. Note that similar observations were reported in a recent study by Horn et al. (Reference Horn, Schmid and Aurnou2022). By employing dynamic mode decomposition on numerical data of a liquid metal, they showed that the torsional and sloshing oscillations originate from the combination of a tilted elliptical mean flow and a higher-order mode consisting of four rolls, aligning with our present findings revealed by SPOD analysis.
The SPOD and resolvent analysis results of the velocity fields measured in a cylindrical RBC with
$Ra=5.2\times 10^9$
and
$\textit{Pr} = 5.7$
. Panels (a,d) present the SPOD energy spectra
$\lambda _n(f)$
and the resolvent gain spectra
$\sigma ^2_n(f)$
, respectively. Panels (b,c) and panels (e, f) depict the corresponding first SPOD and resolvent response modes for frequencies
$f\tau _o \approx 1$
and
$f\tau _o\approx 2$
, where large-scale vortex structures are characterised by streamlines, with the contour indicating the vorticity distribution (scaled by the r.m.s. of the entire field).
To further assess the linear-amplification interpretation in a different geometry, we conducted a resolvent analysis for the cylindrical-cell data; the results are shown in figure 9(d–f). The gain spectra again exhibit dominant frequency peaks near those identified by SPOD, and the leading resolvent response modes remain qualitatively similar to their SPOD counterparts. The two analyses are not expected to coincide exactly, because their equivalence requires an idealised white-noise forcing assumption that is unlikely to hold in turbulent RBC (Towne et al. Reference Towne, Bres and Lele2017, Reference Towne, Schmidt and Colonius2018; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Bres2018). The discrepancies in both peak frequency and spatial detail are therefore informative: they indicate that the forcing statistics matter more strongly in the cylindrical cell, presumably because the additional 3-D dynamics makes the effective forcing more structured. Even so, the persistence of the same dominant travelling-wave-like responses supports the view that linear amplification remains a useful organising principle for the oscillations in cylindrical RBC.
While the above analysis suggests that the linear-amplification interpretation remains useful across different geometric configurations, the apparent differences between the results presented in figures 9 and 5 underscore the significant role of geometric factors in shaping the flow dynamics. Firstly, compared with the results obtained in the quasi-2-D RBC, the separation between the eigenvalues associated with the first and second SPOD modes is more pronounced in the cylindrical RBC. This suggests a more dominant contribution from the primary mode in the cylindrical geometry. Furthermore, the eigenvalue decay rate is significantly steeper in the cylindrical convection cell, with only the primary harmonic frequency standing out clearly. This is in stark contrast to the quasi-2-D configuration, where multiple harmonic frequencies are distinctly observable. The rapid decay of eigenvalues in the cylindrical case implies a more rapid energy transfer to higher modes, which are then quickly damped out. These differences hint at the fact that while the linear-amplification mechanism might be universally applicable, the geometric configuration of the system can significantly influence the flow properties and, as a result, the behaviour of oscillations. Specifically, the relatively loose lateral confinement allows more degrees of freedom in cylindrical cells, consequently, the growth of perturbations are more substantial compared with those in quasi-2-D cases. These stronger perturbations could have dual effects: on the one hand, they may amplify low-frequency oscillations due to the increased energy available for these modes; on the other hand, they might obscure high-frequency oscillations by integrating them into the background turbulence. This integration makes the high-frequency oscillations less distinct and more challenging to isolate from the overall flow dynamics.
As flow oscillations are phenomena ubiquitously observed in turbulent RBC with various geometries, future studies will need to delve deeper into the geometric influence to check the applicability and robustness of the linear-amplification mechanism. This will involve a more detailed examination of how different geometries alter the flow properties and the resulting oscillation behaviours, which is beyond the scope of present study.
5. Concluding remarks
To summarise, we have investigated the oscillation dynamics in turbulent RBC by combining experimental velocity measurements and DNSs. The velocity field in the laterally confined rectangular cell exhibits higher-order harmonic frequencies in addition to the circulation frequency of the mean flow. The SPOD identifies the coherent structures associated with these frequencies as travelling-wave-like patterns with approximate azimuthal periodicity in polar coordinates. A complementary DD resolvent analysis shows that the same dominant structures can be interpreted as selectively amplified responses of an effective linear dynamics about the mean circulation. This interpretation further suggests that the velocity and temperature oscillations should occur concurrently, which is supported by the SPOD analysis of the numerical data in a quasi-2-D RBC. The analysis of the data measured in cylindrical RBC also indicates that the linear-amplification perspective remains applicable in a different geometry, although the agreement between SPOD and resolvent analysis becomes less exact because the effective forcing is more strongly structured.
To the best of our knowledge, the spatially coherent structures underlying specific oscillation frequencies in turbulent RBC, especially in laterally confined geometries, have not previously been characterised as explicitly as done here. The present study identifies these structures and places them in an input–output framework, in which nonlinear interactions act as effective forcing and the mean circulation selects the most amplified responses. We emphasise, however, that the present DD resolvent analysis should be understood as an effective reduced description of the measured large-scale dynamics. It supports the interpretation that the dominant oscillations are governed to leading order by linear amplification about a persistent mean state, rather than proving that all aspects of the turbulent flow are captured by a linear operator. In that sense, the dynamics resembles an intra-well oscillator whose preferred frequencies are set by the local topology of a stable large-scale state, consistent with the strong suppression of flow reversals once a robust LSC is established at sufficiently high
$Ra$
(Brown & Ahlers Reference Brown and Ahlers2006a
; Ni, Huang & Xia Reference Ni, Huang and Xia2015).
However, the absence of high-frequency harmonics in cylindrical RBC indicates that the geometry of the convection cell significantly influences the oscillation characteristics by altering the properties of the mean flow. This observation has prompted us to investigate the impact of geometry on the linear-amplification mechanism. Experimental efforts are on the way and the findings will be reported in the future. Additionally, since the periodic oscillations could be attributed to travelling waves, it is intriguing to examine whether there exist triadic interactions among them, which is another topic for future exploration. A more systematic examination of the conditions under which DD-resolvent and SPOD agree, in view of their shared connections with Koopman-based descriptions, is also left for future work.
Acknowledgements
We thank H.-D. Xi and K.-Q. Xia for their helpful discussions and Y.-Z. Li for sharing us the PIV data measured in the cylindrical convection cell. The authors also gratefully acknowledge the support of computational resources by the Centre for Computational Science and Engineering of Southern University of Science and Technology.
Funding
This work was supported by the National Natural Science Foundation of China (Grant Nos. 11988102 and 12272162) and the Department of Science and Technology of Guangdong Province (Grant Nos. 2019QN01H691 and 2023B1212060001).
Declaration of interests
The authors report no conflict of interest.
Appendix A. The linear dispersion relation
The linear dispersion relation can be derived from the azimuthal Fourier transform of the velocity fluctuation. The derivation is grounded on two key hypotheses: Taylor’s frozen-flow hypothesis (Taylor Reference Taylor1938) and the Kraichnan–Tennekes random sweeping hypothesis (Kraichnan Reference Kraichnan1964; Tennekes Reference Tennekes1975). The mathematical representation is given by the following equation:
\begin{equation} \frac {\partial \hat {u}_\varphi (k,r,t)}{\partial t} + i(U_\varphi +v_\varphi )k\hat {u}_\varphi (k,r,t) = 0. \end{equation}
In this equation,
$k$
denotes the azimuthal wavenumber. The solution to this equation leads to the determination of the space–time energy spectrum, which is expressed as
\begin{equation} E(k,\omega ) = \frac {E(k)}{\sqrt {2\pi k^2v_\varphi ^2}}\exp \Bigg [\frac {-(\omega -U_\varphi k)^2}{2k^2v_\varphi ^2}\Bigg ]. \end{equation}
Here,
$E(k,\omega )$
is a product of a wavenumber spectrum
$E(k)$
and a Gaussian frequency distribution. The Doppler shift is attributed to the constant mean velocity
$U_\varphi$
, while the random advection contributes to the Doppler broadening of the spectrum.
(a) Interpolation of the velocity field in a polar coordinate system, with its origin located at the centre of the convection cell. (b) Frequency–wavenumber diagram with a red dashed line denoting
$f\tau _0=k_\varphi$
.
To obtain the space–time energy spectra, we first interpolate the velocity fields in polar coordinates, as depicted in figure 10(a). Subsequently, we perform a Fourier transformation on both the azimuthal direction and the temporal direction of the interpolated velocity fields. This process is followed by a radially weighted averaging to derive the space–time energy spectrum
$E(k,\omega )$
, which is illustrated in figure 10(b). It is found that the so-obtained energy spectrum aligns with the convection lines
$f\tau _0=k_\varphi$
. This observation confirms that the linear dispersion relation arises from a linear advection mechanism driven by the mean flow.
Appendix B. The SPOD analysis of the temperature and velocity fields obtained from numerical simulation of quasi-2-D RBC
The SPOD analysis results of the velocity and temperature fields obtained from numerical simulation of quasi-2-D RBC at
$Ra=1\times 10^9$
and
$\textit{Pr} = 4.3$
. The SPOD energy spectra
$\lambda _n(f)$
and the first SPOD mode for three oscillation frequencies
$f\tau _0 = 1, 2$
and
$f\tau _o \approx 3$
(dashed lines) of the velocity field (a–d) and temperature field (e–h). The vorticity and temperature are coded in colour and scaled by their r.m.s., respectively.
This appendix is included to support the specific corollary discussed in § 4.1, namely the co-occurrence of temperature and velocity oscillations in the DNS data. It is not intended as a second methodological section or as a separate resolvent analysis. Briefly, we performed SPOD analysis on the velocity and temperature fields obtained from numerical simulation of quasi-2-D RBC, and the results are presented in figure 11. The dominant oscillation frequencies and their corresponding flow patterns are consistent with those presented in figure 8. This consistency confirms that the velocity and temperature fields share similar dynamical properties and their oscillations occur concurrently.
