2.1. Physical model and simulation database

We consider RB convection confined between two parallel plates, with heating applied at the bottom and cooling at the top. The system is governed by the non-dimensional incompressible Navier–Stokes (N–S) equations under the Boussinesq approximation, written as

(2.1) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0, \end{equation}

(2.2) \begin{equation} \frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}=-\boldsymbol{\nabla }\!p+\sqrt {\frac {\textit{Pr}}{Ra}}{\nabla} ^{2}\boldsymbol{u}+\theta \boldsymbol{e}_{\!j}, \end{equation}

(2.3) \begin{equation} \frac {\partial \theta }{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\theta =\frac {1}{\sqrt {\textit{PrRa}}}{\nabla} ^{2}\theta , \end{equation}

where $\boldsymbol{u}$ is the velocity vector, $p$ denotes the kinematic pressure and $\boldsymbol{e}_{\!j}$ represents the wall-normal unit vector. Non-dimensional temperature, pressure and time are denoted by $\theta$ , $p$ and $t$ , respectively. For non-dimensionalisation, the temperature scale is the temperature difference between both plates $\varDelta$ , the length scale is the height of the fluid layer $H$ , the velocity scale is the free-fall velocity $U_{f}=\sqrt {\alpha g H \Delta }$ and the time scale is the free-fall time unit $t_{f}=\sqrt {H/(\alpha g\Delta )}$ . No-slip and impermeability conditions ( $\boldsymbol{u}=0$ ) are enforced at the plates, which are kept at two distinct constant temperatures, while periodic boundary conditions are imposed in both horizontal directions.

The statistical quantities used throughout this work are obtained from fully developed direct numerical simulation (DNS) studies. The numerical simulations were conducted with the open source AFiD code (Verzicco & Orlandi Reference Verzicco and Orlandi1996; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015; Zhu et al. Reference Zhu2018), which has been extensively validated in RB convection studies (Verzicco & Orlandi Reference Verzicco and Orlandi1996; Stevens et al. Reference Stevens, Verzicco and Lohse2010, Reference Stevens, Lohse and Verzicco2011; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015; Zhu et al. Reference Zhu2018). The code employs a finite-difference discretisation that is second-order accurate and energy-conserving, with velocity components defined on a staggered grid. Temporal marching is realised through a third-order Runge–Kutta method for the explicit terms together with a Crank–Nicolson scheme for the implicit contribution. The mesh is uniform in the horizontal directions, whereas the wall-normal direction is refined near the plates using clipped Chebyshev-type clustering to resolve the boundary layers.

The DNS cases considered in this work are grouped into two sets. The first set consists of the large-aspect-ratio cases 1E7 and 1E8, both with $\varGamma =24$ , for which turbulent superstructures are fully developed. Unless otherwise noted, all statistical quantities used in the linear analysis, including the mean temperature profile, turbulent thermal diffusivity and the characteristic scales of the structures, are obtained from this set. The second set consists of the cases 1E7-8, 1E7-1, 1E7-05 and 1E7-025, which are used specifically to examine the influence of aspect ratio. The corresponding computational parameters are summarised separately in tables 1 and 2. The DNS results have been validated against Stevens et al. (Reference Stevens, Blass, Zhu, Verzicco and Lohse2018).

Table 1.

Computational parameters of the large-aspect-ratio DNS cases. Here, $\varGamma =24$ , $N_{y}\times N_{x}\times N_{z}$ denotes the three-dimensional grid resolution, $t_{avg}/t_{f}$ is the sampling time used for statistical averaging and $N_{BL}$ is the number of grid points within the thermal boundary layer.

Table 1. Long description

The table presents computational parameters for two large-aspect-ratio DNS cases, labeled 1E7 and 1E8. Each row lists values for variables such as Ra, Pr, grid resolution denoted as Ny x Nx x Nz, the sampling time ratio tavg/tf, Nusselt number Nu, the ratio λ0/H, and the number of grid points within the thermal boundary layer NBL. The first case, 1E7, has an Ra value of 10⁷, a Pr value of 1.0, a grid resolution of 192 x 6144 x 6144, a sampling time ratio of 200, a Nusselt number of 15.85, a λ0/H ratio of 0.032, and 17 grid points within the thermal boundary layer. The second case, 1E8, has an Ra value of 10⁸, a Pr value of 1.0, a grid resolution of 192 x 6144 x 6144, a sampling time ratio of 200, a Nusselt number of 30.94, a λ0/H ratio of 0.016, and 11 grid points within the thermal boundary layer.

Table 2.

Computational parameters of the DNS cases at different aspect ratios for $Ra=10^{7}$ . The definitions of all variables are the same as in table 1.

Table 2. Long description

The table presents computational parameters for direct numerical simulation (DNS) cases at varying aspect ratios for a Rayleigh number of 10⁷ and a Prandtl number of 1.0. It includes four rows and nine columns. The columns are labeled as Cases, Ra, Pr, N_y x N_x x N_z, Γ, t_avg/t_f, Nu, λ_θ/H, and N_BL. The first row lists the case 1E7-8 with Ra 10⁷, Pr 1.0, grid resolution 192 x 2048 x 2048, Γ 8, t_avg/t_f 1000, Nu 15.70, λ_θ/H 0.032, and N_BL 17. The second row lists the case 1E7-1 with Ra 10⁷, Pr 1.0, grid resolution 192 x 256 x 256, Γ 1, t_avg/t_f 10000, Nu 17.12, λ_θ/H 0.029, and N_BL 17. The third row lists the case 1E7-05 with Ra 10⁷, Pr 1.0, grid resolution 192 x 128 x 128, Γ 0.5, t_avg/t_f 10000, Nu 16.88, λ_θ/H 0.030, and N_BL 17. The fourth row lists the case 1E7-025 with Ra 10⁷, Pr 1.0, grid resolution 192 x 64 x 64, Γ 0.25, t_avg/t_f 10000, Nu 15.00, λ_θ/H 0.033, and N_BL 17. The table provides a detailed comparison of the computational parameters used in the DNS cases to examine the influence of aspect ratio.

2.2. Relevance of the non-normality

We consider the dynamics of small-amplitude perturbations (denoted with primes) to the mean temperature profile $\varTheta (y)$ , which is governed by the non-dimensional linearised N–S equations within the Boussinesq approximation, written as

(2.4) \begin{equation} \partial _{t}\boldsymbol{u}^{\prime }=-\boldsymbol{\nabla }\!p^{\prime }+\sqrt {\frac {\textit{Pr}}{Ra}}{\nabla} ^{2}\boldsymbol{u}^{\prime }+\theta ^{\prime }\boldsymbol{e}_{\!j}, \end{equation}

(2.5) \begin{equation} \partial _{t}\theta ^{\prime }=-\frac {{\rm d}\varTheta (y)}{{\rm d}y}v^{\prime }+\frac {1}{\sqrt {\textit{PrRa}}}\boldsymbol{\nabla }\boldsymbol{\cdot }\left [\left (1+\widetilde {\kappa _{\theta }}(y)\right )\boldsymbol{\nabla }\theta ^{\prime }\right ], \end{equation}

(2.6) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^{\prime }=0. \end{equation}

Here, $\boldsymbol{u}= [u,v,w ]^{T}$ is the non-dimensional velocity with the coordinates $\boldsymbol{x}= [x,y,z ]^{T}$ with $y$ normal to the wall, and $x$ , $z$ spanning horizontal directions. Additionally, $\boldsymbol{e}_{\!j}$ denotes the unit vector in the wall-normal direction. Perturbation components of flow variables are marked using prime notation as follows:

(2.7) \begin{equation} \boldsymbol{u}^{\prime }(\boldsymbol{x},t)=\boldsymbol{u}(\boldsymbol{x},t), \quad \theta ^{\prime }(\boldsymbol{x},t)=\theta (\boldsymbol{x},t)-\varTheta (y). \end{equation}

The turbulent thermal diffusivity, $\widetilde {\kappa _{\theta }}(y)$ , is defined by

(2.8) \begin{equation} \langle v^{\prime }\theta ^{\prime }\rangle =-\frac {\widetilde {\kappa _{\theta }}(y)}{\sqrt {\textit{PrRa}}}\frac {{\rm d}\varTheta (y)}{{\rm d}y}, \end{equation}

where $\left \langle \varphi \right \rangle$ denotes the statistical averaging of an arbitrary variable $\varphi$ in horizontal and time directions. The relation in (2.8) follows from a Boussinesq-type eddy-diffusivity hypothesis. Unlike classical stability analyses based on laminar conductive states, the present formulation considers perturbations around a turbulent mean profile and the turbulent diffusivity is introduced to model the corresponding averaged turbulent transport provided by the fluctuations. No analogous eddy-viscosity term is introduced in the momentum equation, because the mean velocity is zero under horizontal and temporal averaging and thus no mean shear is present.

Following the derivation of the Orr–Sommerfeld and Squire equations (Pope Reference Pope2001), the linearised system (2.4)–(2.6) can be written in matrix form as $ (\partial _{t}\boldsymbol{B}+\boldsymbol{L} )\boldsymbol{X}=0$ , with $\boldsymbol{X}= [v^{\prime },\eta ^{\prime },\theta ^{\prime } ]^{T}$ , where $\eta ^{\prime }$ is the wall-normal vorticity. Detailed operator definitions, and in particular $\boldsymbol{B}$ , are given in Appendix A. The linear operator takes the form

(2.9) \begin{equation} \boldsymbol{L}=\left [\begin{array}{ccc} L_{v} & 0 & \partial _{y}^{2}-{\nabla} ^{2}\\[5pt] 0 & L_{\eta } & 0\\[5pt] \frac {{\rm d}\varTheta }{{\rm d}y} & 0 & L_{\theta } \end{array}\right ]\!, \end{equation}

where the off-diagonal coupling elements, ${\rm d}\varTheta /{\rm d}y$ and $ (\partial _{y}^{2}-{\nabla} ^{2} )$ , represent the contributions of the mean temperature gradient and the buoyancy term in RB convection, respectively, and act to enhance the non-normality of the linearised N–S system.

The importance of linear mechanisms in sustaining shear turbulence was highlighted by Kim & Lim (Reference Kim and Lim2000), who showed in turbulent channel flow that removing the lift-up term, a dominant non-modal component of transient growth, led to the decay of turbulence. In turbulent RB convection, the governing equations contain an analogous coupling term associated with the mean temperature gradient ${\rm d}\varTheta /{\rm d}y$ . In the temperature equation, this term acts as a source-like contribution, while in the full linearised system, it constitutes an off-diagonal coupling between the velocity and temperature perturbations. To examine its role, we conduct a numerical experiment in which the off-diagonal elements linked to ${\rm d}\varTheta /{\rm d}y$ are removed directly from the full nonlinear N–S equations (see Appendix A for details). In this modified system, the temperature equation becomes

(2.10) \begin{equation} \partial _{t}\theta ^{\prime }=-\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{\nabla }\theta ^{\prime }+\frac {1}{\sqrt {\textit{PrRa}}}{\nabla} ^{2}\theta ^{\prime }+\frac {{\rm d}\varTheta }{{\rm d}y}v^{\prime }, \end{equation}

where the additional term $({\rm d}\varTheta /{\rm d}y)v^{\prime }$ cancels the original coupling through the mean temperature gradient. This modification reduces the non-normality of the linearised system and suppresses linear forcing induced by the turbulent temperature profile, while keeping all nonlinear terms unchanged. Therefore, the experiment isolates the role of this linear coupling term from the remaining nonlinear source-like contributions.

The initial fields for the numerical experiments are based on the fully developed turbulent fields obtained from DNS. Simulations are performed using the open-source code AFiD (Verzicco & Orlandi Reference Verzicco and Orlandi1996; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015; Zhu et al. Reference Zhu2018), with computational parameters summarised in table 1. The computational details can be found in Appendix A.

In the numerical experiment, we modified the N–S equations in the AFiD code as given in (2.10), where the mean temperature profile $\varTheta$ was updated at each time step. The temporal evolution of $\textit{Nu}$ following the equation modification is shown in figure 2. Upon removal of the off-diagonal term corresponding to ${\rm d}\varTheta /{\rm d}y$ , $\textit{Nu}$ drops sharply to approximately $1$ , indicating rapid flow laminarisation. This demonstrates that the off-diagonal term and its induced non-normality are essential for sustaining nonlinear turbulence. In RB convection, neither superstructures nor small-scale structures can be sustained solely by nonlinear effects without the off-diagonal term ${\rm d}\varTheta /{\rm d}y$ .

Figure 2.

Time evolution of $\textit{Nu}$ on the wall. The time is non-dimensionalised by the free-fall time unit $t_{f}=\sqrt {H/(\alpha g\Delta )}$ . $t=0$ corresponds to the instant when the temperature equation is modified. The black dashed line denotes $\textit{Nu}=1$ .

Figure 2. Long description

A line graph showing the time evolution of Nusselt number on the wall for two different Rayleigh numbers. The x-axis represents time non-dimensionalized by the free-fall time unit, ranging from 0 to 300. The y-axis represents the Nusselt number, ranging from 0 to 30. Two data lines are present: one for Ra equals 10 to the power of 7 in blue and another for Ra equals 10 to the power of 8 in magenta. Both lines start at a high Nusselt number and decrease over time, approaching a steady state. The black dashed line denotes a Nusselt number of 1. All reported values are computed numerically.

2.3. Linear energy amplification

Building on the demonstrated crucial role of the off-diagonal term ${\rm d}\varTheta /{\rm d}y$ in turbulence maintenance, its enhancement of non-normality in the N–S system warrants further investigation, particularly regarding possible linear energy amplification of some certain disturbances. In wall-bounded shear flows, linear energy amplification theory has proven highly effective in uncovering the genesis of both small-scale and large-scale structures (Butler & Farrell Reference Butler and Farrell1993; del Álamo & Jiménez Reference del Álamo and Jiménez2006; Pujals et al. Reference Pujals, García-Villalba, Cossu and Depardon2009; McKeon & Sharma Reference McKeon and Sharma2010). del Álamo & Jiménez (Reference del Álamo and Jiménez2006) established that optimal linear modes, exhibiting maximal transient growth within turbulent mean profiles and eddy viscosities, well describe both near-wall streaks and large-scale motions. These analyses further revealed counter-rotating roll cells that closely match the characteristic scales of structures widely documented in turbulent shear flows.

Motivated by these insights, we compute the growth rate $G$ in turbulent convection, defined as

(2.11) \begin{equation} G(k_{x},t)=\underset {\forall \boldsymbol{q}(t=0)\neq 0}{\max }\left [\frac {\left \Vert \boldsymbol{q}(k_{x},t)\right \Vert }{\left \Vert \boldsymbol{q}(k_{x},t=0)\right \Vert }\right ]\!, \end{equation}

where the norm $\left \Vert \boldsymbol{q}(k_{x},t)\right \Vert$ characterises the amplitude of perturbation $\boldsymbol{q}$ at the horizontal wavenumber $k_{x}$ .

Given the horizontal symmetry of RB convection, subsequent analyses consider exclusively two-dimensional configurations. The wavenumber $k_x$ here denotes the horizontal wavenumber. Due to the inherent horizontal symmetry of the system, such two-dimensional perturbations may be regarded as disturbances with an arbitrary horizontal orientation in the corresponding three-dimensional case. Consequently, the perturbation $\boldsymbol{q}= [\hat {u}(y),\hat {v}(y),\hat {\theta }(y) ]^{T}\times \exp [\mathrm{i}(k_{x}x-\sigma t) ]$ , where $\hat {\varphi }$ denotes the Fourier coefficient of an arbitrary variable $\varphi$ . The norm of $\boldsymbol{q}$ is defined as

(2.12) \begin{equation} \left \Vert \boldsymbol{q}\right \Vert =\frac {1}{H}\intop _{0}^{H}\left (\left |\hat {u}\right |^{2}+\left |\hat {v}\right |^{2}+\gamma \left |\hat {\theta }\right |^{2}\right )\mathrm{d}y, \end{equation}

where $\gamma \gt 0$ weights the amplitude of thermal fluctuation relative to the kinetic energy and the bottom plate is chosen as the baseline $y=0$ . Unlike wall-bounded shear flows that typically consider only velocity fluctuations (Butler & Farrell Reference Butler and Farrell1993; del Álamo & Jiménez Reference del Álamo and Jiménez2006), RB convection also requires examining temperature fluctuations. Here, we include the squared temperature fluctuation term in the growth rate, which plays a role analogous to the potential energy contribution in strongly stratified shear layers (Kaminski et al. Reference Kaminski, Caulfield and Taylor2014). Such a combined kinetic–thermal energy measure provides a consistent framework for evaluating the total perturbation energy gain (Kaminski et al. Reference Kaminski, Caulfield and Taylor2014; Bose, Kannan & Zhu Reference Bose, Kannan and Zhu2024; Secchi et al. Reference Secchi, Gatti, Piomelli and Frohnapfel2025). By varying $\gamma$ , we could study modes that are dominated by either the velocity or the temperature field.

The growth rate is computed via the linearised N–S equations (2.4)–(2.6) using the methodology established by Reddy & Henningson (Reference Reddy and Henningson1993). The computational framework leverages the open-source Dedalus package (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020), employing spectral discretisation in horizontal directions and Chebyshev polynomials in the wall-normal direction. Analogous to the transient growth work in channel flow (del Álamo & Jiménez Reference del Álamo and Jiménez2006), we incorporate the mean temperature gradient and turbulent thermal diffusivity into (2.5). The gradient ${\rm d}\varTheta /{\rm d}y$ and the effective diffusivity $\widetilde {\kappa _{\theta }}$ defining the turbulent background are derived from statistical averaging of DNS data for RB convection, with computational parameters summarised in table 1. Figure 3(a,b) displays the distributions of the mean temperature gradient ${\rm d}\varTheta /{\rm d}y$ and turbulent thermal diffusivity $\widetilde {\kappa _{\theta }}$ across wall-normal heights, correspondingly. Additionally, ${\rm d}\varTheta /{\rm d}y$ remains negative under all conditions, with steeper gradients concentrated near the walls. These near-wall gradients intensify with increasing $Ra$ . Towards the flow centreline away from walls, the gradient diminishes towards nearly zero, exhibiting variations exceeding three orders of magnitude across different heights. Meanwhile, $\widetilde {\kappa _{\theta }}$ increases monotonically with wall-normal distance, peaking at the flow centre, and grows with $Ra$ throughout the domain.

Figure 3.

Turbulent profiles and corresponding linear energy amplification. (a, b) Wall-normal distributions of the mean temperature gradient ${\rm d}\varTheta /{\rm d}y$ and turbulent thermal diffusivity $\widetilde {\kappa _{\theta }}$ . (c–f) Growth rate $G$ as a function of time $t$ and horizontal wavenumber $k_{x}$ for different cases: (c) $Ra=10^{7}$ , $\gamma =0.01$ ; (d) $Ra=10^{8}$ , $\gamma =0.01$ ; (e) $Ra=10^{7}$ , $\gamma =100.0$ ; and (f) $Ra=10^{8}$ , $\gamma =100.0$ .

Figure 3. Long description

The image contains six graphs. The first two graphs (a and b) show wall-normal distributions of the mean temperature gradient and turbulent thermal diffusivity for two different Rayleigh numbers. The next four graphs (c, d, e, and f) display the growth rate as a function of time and horizontal wavenumber for different cases. Graphs c and d represent cases with Rayleigh numbers of 10⁷ and 10⁸, respectively, while graphs e and f show cases with different parameters. The color scales indicate the magnitude of the growth rate, with warmer colors representing higher values. The graphs illustrate the coexistence of superstructures and small-scale structures in turbulent Rayleigh-Bénard convection.

Figure 3(c–f) displays the maximum energy amplification $G$ for two typical weighting coefficients ( $\gamma =0.01$ and $100.0$ ), plotted against time $t$ and horizontal wavenumber $k_{x}$ . Critically, across all positive $\gamma$ , the growth rate $G$ exhibits two primary distribution types. Here, $\gamma$ is a weighting parameter in the perturbation norm, and the values $\gamma =0.01$ and $\gamma =100$ are shown as representative examples of two broader regimes rather than isolated cases. The evolution of the growth rate over the range from $\gamma =0.001$ to $1000$ is provided in Appendix B. For small $\gamma$ values, velocity perturbations dominate the weighting in $\left \Vert \boldsymbol{q}\right \Vert$ , producing the characteristic $G$ distribution shown in figure 3(c,d). These plots reveal pronounced large-scale linear energy amplification near $k_{x}=O(1)$ persisting beyond $20$ free-fall time units. The distribution of $G$ remains consistent for $\gamma \leq 0.1$ , indicating a persistent physical mechanism. Exemplary views of corresponding optimal solutions in physical space are displayed in figure 4(a,b). The perturbations extend across the entire wall-normal domain, consisting of normally penetrating large-scale cold and hot regions along with counter-rotating roll cells. In particular, this coupled organisation of large-scale temperature regions and embedded counter-rotating circulations is highly consistent with the conditionally averaged superstructures reported by Blass et al. (Reference Blass, Verzicco, Lohse, Stevens and Krug2021). Hot regions align with upward flows in the large-scale circulation, while cold regions correspond to downward flows. Between these zones, the circulation continuously transports fluid from plume-impacting to plume-emitting sections, sustaining a closed convective loop. A quantitative comparison between the characteristic scales of these optimal patterns and those of the superstructures obtained from DNS is provided in table 3, which will be discussed further later.

Table 3.

Structural length scale comparisons between linear optimal modes and DNS statistics. We display the optimal modes corresponding to the maximum linear energy amplification at time $t_g$ . The wavenumbers and wavelengths of the structures in DNS are identified from the peak positions of the premultiplied coherence spectra, in agreement with the results of Blass et al. (Reference Blass, Verzicco, Lohse, Stevens and Krug2021).

Table 3. Long description

The table presents a comparison of structural length scales between linear optimal modes and direct numerical simulation (DNS) statistics. It includes data for two Rayleigh numbers (10⁷ and 10⁸) and two scales (large and small). The columns are labeled with Ra, Scale, t_g, k_z H (linear), k_z H (DNS), λ_x/H (linear), and λ_x/H (DNS). The table shows that for Ra 10⁷ and large scale, t_g is 4.1, k_z H (linear) is 1.1, k_z H (DNS) is 1.1, λ_x/H (linear) is 5.7, and λ_x/H (DNS) is 5.9. For Ra 10⁷ and small scale, t_g is 1.1, k_z H (linear) is 17, k_z H (DNS) is 18, λ_x/H (linear) is 0.37, and λ_x/H (DNS) is 0.35. For Ra 10⁸ and large scale, t_g is 5.6, k_z H (linear) is 1.1, k_z H (DNS) is 1.0, λ_x/H (linear) is 5.9, and λ_x/H (DNS) is 6.3. For Ra 10⁸ and small scale, t_g is 1.0, k_z H (linear) is 42, k_z H (DNS) is 35, λ_x/H (linear) is 0.15, and λ_x/H (DNS) is 0.18. The table highlights the differences in structural length scales for various conditions, providing insights into the energy amplification and flow characteristics.

When $\gamma$ increases, temperature perturbations gain greater weight, shifting the linear energy amplification to a distinct physical mechanism. Figure 3(e,f) reveals energy amplification concentrated at small scales ( $k_{x}=10{-} 100$ ) within $10{-}20$ free-fall time units. This distribution pattern is also widely observed across cases with $\gamma \geq 10$ . Typical corresponding optimal solutions in physical space are displayed in figure 4(c,d). Unlike the large-scale results, the perturbations corresponding to small-scale energy amplification are exclusively localised near the wall, with scales of the order of the thermal boundary layer thickness, resembling the small-scale structures in turbulent RB convection. Small-scale perturbations similarly feature paired hot/cold regions and roll cells. Near the lower wall, these cells emit hot flow upward, while sweeping cold fluid towards the wall; vice versa at the upper wall.

Figure 4.

Linear optimal modes. Developed optimal solutions for the maximum energy amplification at structure generation time $t_g$ . The results represent the state after evolving the initial optimal solutions for time $t_g$ . (a) Case at $Ra=10^{7}$ , $\gamma =0.01$ , $t_g=4.1t_{f}$ for large-scale part; (b) case at $Ra=10^{8}$ , $\gamma =0.01$ , $t_g=5.6t_{f}$ for large-scale part; (c) case at $Ra=10^{7}$ , $\gamma =100.0$ , $t_g=1.1t_{f}$ for small-scale part; (d) case at $Ra=10^{8}$ , $\gamma =100.0$ , $t_g=1.0t_{f}$ for small-scale part.

Figure 4. Long description

The image consists of four heat maps labeled (a), (b), (c), and (d), each representing different cases of linear optimal modes for maximum energy amplification in fluid dynamics. The heat maps display color gradients ranging from red to blue, indicating varying levels of energy amplification. The x and y axes represent spatial coordinates, with the color scale on the right indicating the magnitude of energy amplification. In (a) and (b), the large-scale part of the fluid dynamics is shown, while (c) and (d) focus on the small-scale part. The color intensity varies across the maps, with red areas indicating higher energy amplification and blue areas indicating lower amplification. The heat maps illustrate the distribution and interaction of large-scale circulations and small-scale turbulent fluctuations in fluid dynamics.

Here, significant linear energy amplification occurs under both mechanisms. We emphasise that these two patterns should not be interpreted as purely velocity structures at large scales and purely temperature structures at small scales; in both cases, the optimal perturbations remain coupled velocity–temperature structures and the presence of both components is essential for the observed transient growth. The reason why different values of $\gamma$ highlight these two distinct regimes requires further investigation. Derived from the linearised N–S equations (2.4)–(2.6), the budget equations for velocity and temperature perturbations can be expressed as

(2.13) \begin{equation} \partial _{t}\left (\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{u}^{\prime }\right )=\underbrace {2v^{\prime }\theta ^{\prime }}_{P_{u}}+\underbrace {\left (-\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{\nabla }\!p^{\prime }\right )}_{\textit{VP}}+\underbrace {\sqrt {\frac {\textit{Pr}}{Ra}}\boldsymbol{u}^{\prime }\boldsymbol{\cdot }{\nabla} ^{2}\boldsymbol{u}^{\prime }}_{\varepsilon _{u}}, \end{equation}

(2.14) \begin{equation} \partial _{t}\left (\theta ^{\prime }\theta ^{\prime }\right )=\underbrace {\left (-2\frac {{\rm d}\varTheta }{{\rm d}y}v^{\prime }\theta ^{\prime }\right )}_{P_{\theta }}+\underbrace {\frac {\theta ^{\prime }}{\sqrt {\textit{PrRa}}}\boldsymbol{\nabla }\boldsymbol{\cdot }\left [\left (1+\widetilde {\kappa _{\theta }}\right )\boldsymbol{\nabla }\theta ^{\prime }\right ]}_{\varepsilon _{\theta }}. \end{equation}

Here, $P_{u}$ and $P_{\theta }$ denote the production terms for the velocity and temperature perturbations associated with the perturbation $v^{\prime }\theta ^{\prime }$ , respectively. Additionally, $\textit{VP}$ represents the velocity–pressure gradient term, while $\varepsilon _{u}$ corresponds to the dissipation term of the velocity perturbation and $\varepsilon _{\theta }$ denotes the effective dissipation term of the temperature perturbation, accounting for nonlinear effects. For the same $v^{\prime }\theta ^{\prime }$ perturbation, the temperature mode’s production term is scaled by $-{\rm d}\varTheta /{\rm d}y$ relative to the velocity mode. As shown in figure 3(a), the $-{\rm d}\varTheta /{\rm d}y$ profile in turbulent environments exhibits substantial height-dependent variations. This gradient remains strong near the thermal boundary layer height, but rapidly diminishes to near-zero values in the elevated flow core region, exhibiting over $1000$ -fold variations between these regions. This characteristic feature of turbulent RB convection consequently generates significant height-dependent disparities in the relative magnitudes of velocity and temperature modes.

In the flow core region away from walls, the mean temperature gradient approaches zero. Consequently, the production term for velocity modes significantly exceeds that of temperature modes, establishing velocity mode dominance. This vertically extensive core region enables perturbations to fully develop into larger scales spanning the entire wall-normal domain. Such dynamics is more readily captured at smaller $\gamma$ values, where velocity modes carry greater weight. Conversely, near the wall, the substantial mean temperature gradient causes the temperature mode production term to dominate over its velocity counterpart. Here, temperature modes prevail, with perturbations confined to scales characteristic of the thermal boundary layer where strong temperature gradients exist. This regime is better resolved at larger $\gamma$ values, which amplify the weighting of temperature modes.

Overall, the sharply varying mean temperature gradient accounts for the observed scale separation across different $\gamma$ values. The large- and small-scale perturbations identified in linear energy amplification qualitatively resemble superstructures and small-scale structures in turbulent RB convection, though specific horizontal wavelengths require further validation against DNS data. Critically, figure 3(c–f) demonstrates that the optimal mode, which exhibits maximal growth rate $G$ at different $t$ , progressively increases in scale over time. However, in actual DNS environments, nonlinear effects ultimately limit structural growth. Therefore, following the work by Butler & Farrell (Reference Butler and Farrell1993), we introduce a structure generation time scale $t_g$ derived from DNS statistics and consider optimal modes maximising growth rate at this specific $t_g$ .

Here, the structure generation time $t_g$ is defined as $t_g=E/P$ , where $E$ denotes the fluctuation energy at the corresponding scale and $P$ represents the energy input per unit time. Leveraging the distinct scale separation between superstructures and small-scale structures, we implement a spectral cutoff at horizontal wavelength $\lambda _{x}=\pi H$ to isolate these components, following Blass et al. (Reference Blass, Verzicco, Lohse, Stevens and Krug2021)’s methodology, where $\lambda _{x}=2\pi /k_{x}$ . Subsequent notation uses subscript ‘ $L$ ’ for large-scale structures ( $\lambda _{x}\gt \pi H$ ) and ‘ $S$ ’ for small-scale structures ( $\lambda _{x}\lt \pi H$ ). We first examine linear energy amplification for large-scale components, which is dominated by velocity modes. Consequently, the fluctuation energy primarily corresponds to large-scale turbulent kinetic energy. For these scales, energy input consists solely of the production term, while inter-scale transport and dissipation act as energy output (Bose et al. Reference Bose, Kannan and Zhu2024). Thus, $E_L$ and $P_L$ can be written into

(2.15) \begin{equation} E_{L}=\int _{0}^{H}\left \langle \boldsymbol{u}_{L}^{\prime }\boldsymbol{\cdot }\boldsymbol{u}_{L}^{\prime }\right \rangle {\rm d}y, P_{L}=\int _{0}^{H}2\left \langle v_{L}^{\prime }\theta _{L}^{\prime }\right \rangle {\rm d}y. \end{equation}

Conversely, for small-scale components where temperature modes govern energy amplification, fluctuations predominantly represent temperature variances. Here, energy input includes both production and inter-scale transport, balanced exclusively by dissipation (Bose et al. Reference Bose, Kannan and Zhu2024; Secchi et al. Reference Secchi, Gatti, Piomelli and Frohnapfel2025). This yields the following equations for $E_S$ and $P_S$ :

(2.16) \begin{equation} E_{S}=\int _{0}^{H}\left \langle \theta _{S}^{\prime }\theta _{S}^{\prime }\right \rangle \,{\rm d}y,\quad P_{S}=\int _{0}^{H}\frac {2}{\sqrt {\textit{PrRa}}}\left \langle \boldsymbol{\nabla }\theta _{S}^{\prime }\boldsymbol{\cdot }\boldsymbol{\nabla }\theta _{S}^{\prime }\right \rangle\, {\rm d}y. \end{equation}

Based on these, we obtain $t_g$ for superstructures and small-scale structures, while the former is $4{-}6$ free-fall time units and the latter is approximately $1$ free-fall time unit. As demonstrated in table 3, the structural length scales of the linear optimal perturbations (corresponding to maximum energy amplification at $t_g$ ) show agreement with the DNS statistics across both the large-scale and small-scale branches, for both $Ra = 10^7$ and $10^8$ cases. This quantitative agreement is a key support for our conclusion that the dominant large- and small-scale structures in turbulent RB convection are captured well by the optimal linear perturbations. Furthermore, Appendix C presents an additional DNS comparison based on the normalised wall-normal distribution of the temperature-fluctuation variance for both branches. While the optimal perturbations represent the initial linear organisation and DNS reflects the subsequent nonlinear development, the main wall-normal trends remain similar.

Figure 5.

Turbulent profiles and linear energy amplification at different aspect ratios. $Ra=10^{7}$ . (a–c) Wall-normal distributions of the mean temperature gradient ${\rm d}\varTheta /{\rm d}y$ and turbulent thermal diffusivity $\widetilde {\kappa _{\theta }}$ . In panel (b), the vertical axis is shown on a logarithmic scale to facilitate comparison of the temperature gradient in the bulk region. (d–k) Growth rate $G$ as a function of time $t$ and horizontal wavenumber $k_{x}$ for different cases. The two columns from left to right correspond to $\gamma =0.01$ and $\gamma =100.0$ , while the rows from top to bottom show results for aspect ratios $\varGamma =8$ , $\varGamma =1$ , $\varGamma =0.5$ and $\varGamma =0.25$ . The black lines represent the minimum wavenumber attainable by the structures at the corresponding aspect ratio.

Figure 5. Long description

The image contains multiple graphs showing turbulent profiles and linear energy amplification at different aspect ratios. The first set of graphs (a-c) displays wall-normal distributions of the mean temperature gradient and turbulent thermal diffusivity. Panel (b) uses a logarithmic scale on the vertical axis to compare the temperature gradient in the bulk region. The second set of graphs (d-k) shows the growth rate as a function of time and horizontal wavenumber for different cases. The two columns correspond to different values, while the rows show results for aspect ratios 24, 8, 1, 0.5, and 0.25. The black lines represent the minimum wavenumber attainable by the structures at the corresponding aspect ratio. All reported values are computed numerically.

Furthermore, the influence of the aspect ratio $\varGamma$ on linear energy amplification warrants further discussion. As the aspect ratio decreases, the characteristic scales of the flow structures also change. For example, the integral length scale present in the flow domain decreases continuously once $\varGamma \lt 8$ (Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018). To examine this effect, we consider the case $Ra=10^{7}$ with aspect ratios $\varGamma =0.25$ , $0.5$ , $1$ , $8$ and $24$ . The mean temperature gradient ${\rm d}\varTheta /{\rm d}y$ and turbulent thermal diffusivity $\widetilde {\kappa _{\theta }}$ used for the linear energy amplification analysis are obtained from statistics of fully developed DNS results, with the corresponding computational parameters listed in table 2. The results are shown in figure 5(a–c). Across different $\varGamma$ , the mean temperature gradient profiles generally collapse, with only small variations near the wall. In the bulk region, however, the temperature gradient gradually departs from zero as the aspect ratio decreases. This trend becomes particularly pronounced for the cases with $\varGamma \leq 0.5$ , for which the temperature gradient at the flow centre is more than ten times larger than that of the $\varGamma =24$ case. In contrast, $\widetilde {\kappa _{\theta }}$ exhibits pronounced differences. As $\varGamma$ decreases, the turbulent thermal diffusivity in the flow centre drops rapidly. Compared with the baseline case of $\varGamma =24$ , the $\varGamma =8$ profile largely collapses onto it except in the central region. However, when the aspect ratio is reduced to $\varGamma \leq 1$ , the central peak vanishes and is replaced by two off-centre peaks, leading to increasingly larger deviations both near the boundary layers and in the central region.

The growth rates $G$ at different aspect ratios are shown in figure 5(d–k). At $\varGamma =8$ , the growth rate distribution remains similar to that at $\varGamma =24$ , with wavenumber truncation cutting off part of the large-scale amplification (dominant wavelength of order $\varGamma$ ), while leaving the small-scale amplification nearly unaffected. When the aspect ratio is reduced to $\varGamma =1$ , wavenumber truncation together with the mean profile variations confines large-scale amplification to fluctuations whose wavelength equals the aspect ratio, where its magnitude is already weak; for $\varGamma \lt 1$ , it essentially vanishes. Meanwhile, small-scale amplification is less sensitive to truncation at $\varGamma =1$ , but differs from the $\varGamma =24$ case due to mean profile variations. When $\varGamma \lt 1$ , truncation effects and profile modifications jointly exert a strong influence on the growth rate; in particular, at $\varGamma =0.25$ , the peak of small-scale amplification shifts to disturbances of order $\varGamma$ . Overall, for large aspect ratios, the variation in linear energy amplification is primarily governed by wavenumber truncation, whereas at small aspect ratios, changes in the mean profile also play a significant role.