2.1. Governing equations

We consider a layer of fluid sandwiched between two extended parallel plates, separated by depth $2h$ , with a uniform upper and lower plate temperature $T_U$ and $T_L$ . The horizontal size of the domain of interest is set with equal length and span, $L = L_x = L_z$ . The fluid is unstably stratified such that $\Delta T = T_L - T_U \gt 0$ . The fluid has a density of $\rho$ , kinetic viscosity, $\nu$ , and thermal diffusivity, $\kappa$ . The non-dimensionalised governing equations, composed of the Navier–Stokes equations with Boussinesq approximation, the energy equation and the divergence free constrain, are given by

(2.1a) \begin{equation} \frac {\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u} = -\boldsymbol{\nabla }\!p + \frac {1}{\textit{Re}}{\nabla} ^2 \boldsymbol{u} + \frac {Ra}{8\textit{PrRe}^2}\theta \boldsymbol{{j}}, \end{equation}

(2.1b) \begin{equation} \frac {\partial \theta }{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla })\theta = \frac {1}{\textit{RePr}}{\nabla} ^2 \theta , \end{equation}

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

with the following boundary conditions at the wall:

(2.2) \begin{equation} \boldsymbol{u}|_{y=\pm 1} = 0, \quad \theta |_{y = -1} = 1, \quad \theta |_{y=1} = 0, \end{equation}

and periodic boundary conditions in the planar $x$ and $z$ directions. The unit vector $\boldsymbol{j}$ is in the $y$ direction. The temporal coordinate $t$ denotes the time non-dimensionalised by the advective time scale, $U_c/h$ . Here, $\boldsymbol{u}(\boldsymbol{x},t)$ is the velocity vector non-dimensionalised by the laminar centreline velocity, $U_c$ (i.e. $U_{\textit{lam}}(y)= U_c(1-y^2)$ ), $p$ the pressure scaled by the dynamic pressure scale, $\rho U_c^2$ , and $\theta (\equiv (T-T_U)/\Delta T)$ the non-dimensionalised temperature with $T$ being the absolute temperature. We note that the rescaled term $Ra/8$ in the momentum forcing terms is equivalent to the Rayleigh number scaled based on the half-depth $h$ , since $Ra$ is scaled based on depth, $2h$ , as in classical RBC. The Rayleigh number $Ra$ , Reynolds number $\textit{Re}$ and Prandtl number $\textit{Pr}$ are defined in § 1. We set $\textit{Pr} = 1$ in this study. For the Rayleigh–Bénard convection problem, we note that (2.1) becomes singular at $\textit{Re} = 0$ . In such cases, we solve the non-dimensionalised incompressible Navier–Stokes equations based on thermal length, velocity and temporal scales shown in Appendix A.

2.2. Numerical methods

The governing equations are solved numerically using an open-source spectral/hp-element package, Nektar++ (Cantwell et al. Reference Cantwell2015; Moxey et al. Reference Moxey2020). The computational mesh consists of two-dimensional (2-D) quadrilateral elements in the $y{-}z$ plane generated using Gmsh (Geuzaine & Remacle Reference Geuzaine and Remacle2009) and then imported into Nektar++ using Nekmesh (Green et al. Reference Green, Kirilov, Turner, Marcon, Eichstädt, Laughton, Cantwell, Sherwin, Peiró and Moxey2024). We discretise the spatial domain based on the quasi-three-dimensional (quasi-3-D) approach, employing spectral/hp elements in the $y{-}z$ plane and Fourier expansions in the streamwise direction $x$ (Karniadakis Reference Karniadakis1989). The discretised equations are solved using a velocity correction scheme, based on a second-order stiffly splitting scheme, where the nonlinear advection and forcing terms are treated explicitly, while the pressure and diffusion terms are treated implicitly (Karniadakis, Israeli & Orszag Reference Karniadakis, Israeli and Orszag1991; Guermond & Shen Reference Guermond and Shen2003). The $3/2$ and polynomial dealiasing rule for the Fourier expansions and spectral/hp elements are applied during the evaluation of the nonlinear advection terms. We refer to the solution obtained at the end of the velocity correction scheme as the homogeneous velocity $\boldsymbol{u}_h$ . To sustain the flow for $\textit{Re}\gt 0$ , we use a Green’s function approach (Chu et al. Reference Chu, Henderson and Karniadakis1992) to impose a constant flow rate,

(2.3) \begin{equation} U_b = Q(\boldsymbol{u})=\frac {1}{2 L_z h}\int _{ y,z} \boldsymbol{u}(0,y,z) \; \mathrm{d}y\,\mathrm{d}z, \end{equation}

where $U_b$ and $Q(\boldsymbol{\cdot })$ refer to the desired flow rate and flow rate operator. A correction velocity, $\boldsymbol{u}_{\textit{corr}}$ , is obtained by solving the linear Stokes equations with unit forcing once and is stored for reuse. At the end of every time step, the final velocity field, $\boldsymbol{u}$ , is then updated by adding this correction velocity to the homogeneous velocity obtained from the velocity correction scheme,

(2.4) \begin{equation} \boldsymbol{u} = \boldsymbol{u}_h + \gamma \boldsymbol{u}_{\textit{corr}}, \end{equation}

where $\gamma$ , defined as

(2.5) \begin{equation} \gamma = \frac { U_b - Q(\boldsymbol{u}_h)}{Q(\boldsymbol{u}_{\textit{corr}})}, \end{equation}

is adjusted to satisfy the desired flow rate, $U_b$ . The flow rate, $U_b$ , is related to the laminar centreline velocity $U_c = 3/2 U_b$ , which defines the Reynolds number, $\textit{Re} = U_c h / \nu$ . For more details on the numerical method, the reader is referred to Hossain, Cantwell & Sherwin (Reference Hossain, Cantwell and Sherwin2021).

2.3. Parametric study in the $Ra$ – $\textit{Re}$ space

We consider fifty-two numerical simulations at $\textit{Re} = 0, 0.1, 1, 10, 100, 500, 750, 1000,{} 1050$ , $2000$ and $Ra = 0, 2000, 3000, 5000, 8000, 10000$ at $\textit{Pr} = 1$ with a large aspect ratio of $\varGamma = 16\pi$ . The initial conditions of all cases are sampled from a statistically stationary solution. Laminar solutions are obtained for $Ra = 0$ , $\textit{Re} \leqslant 1000$ , with the given set of parameters and are omitted from the analysis. For all of the cases considered here, we maintain the same spatial resolution of $(\Delta x, \Delta y|_{y=\pm 1}, \Delta z) = (0.25\pi , 0.0549, 0.1\pi )$ with polynomial order $P=4$ , except for the end case of $\textit{Re} = 2000$ , where the resolution in the streamwise direction is increased from $\Delta x = 0.25\pi$ to $\Delta x = 0.125\pi$ . This corresponds to a resolution in wall units of $(\Delta x^+, \Delta y^+|_{y=\pm 1}, \Delta z^+) = (37.27, 1.30, 7.49)$ . Then, we considered a mesh independence study for the end case of $Ra = 10000$ , $\textit{Re} = 2000$ , where increasing the resolution in streamwise direction or the polynomial order by 1 led to a ${\lt} 1\,\%$ change in near-wall transport properties defined by the Nusselt number, $\textit{Nu} = -2\langle \mathrm{d}\theta /\mathrm{d}y_{y=-1}\rangle _{x,z} h/\Delta T$ , and shear, $\langle \tau _w \rangle _{x,z} = \langle \mathrm{d}{u}/\mathrm{d}y|_{y=-1} \rangle _{x,z}$ , on the lower wall, where $\langle \boldsymbol{\cdot }\rangle _{x,z} = 1 / (L_xL_z) \int \; (\boldsymbol{\cdot }) \; \mathrm{d}x\,\mathrm{d}z$ refers to the plane-averaged operator. Their temporal numerical resolutions and time-integration horizon, $\zeta$ , are described in Appendix B. The temporal resolution between the numerical simulations differs due to time scales arising from the different flow physics, as we shall see later.

2.4. Linear stability analysis

In § 4, we will perform numerical experiments where small disturbances are added along the unstable manifolds of longitudinal rolls. To determine the unstable manifolds, we consider a small disturbance about the longitudinal roll state,

(2.6a) \begin{equation} \boldsymbol{u}(\boldsymbol{x}, t) = \boldsymbol{u}_{\textit{LR}}(\boldsymbol{x}) + \boldsymbol{\hat {u}}(\boldsymbol{x}, t), \end{equation}

(2.6b) \begin{equation} \theta (\boldsymbol{x}, t) =\theta _{\textit{LR}}(\boldsymbol{x}) + \hat {\theta }(\boldsymbol{x}, t), \end{equation}

(2.6c) \begin{equation} p(\boldsymbol{x}, t) = p_{\textit{LR}}(\boldsymbol{x}) + \hat {p}(\boldsymbol{x}, t), \end{equation}

where $\boldsymbol{q} = [\boldsymbol{u}, \theta , p]^T$ , $\boldsymbol{q}_{\textit{LR}}=[\boldsymbol{u}_{\textit{LR}}, \theta _{\textit{LR}}, p_{\textit{LR}}]^T$ and $\boldsymbol{\hat {q}}=[\boldsymbol{\hat {u}}, \hat {\theta }, \hat {p}]^T$ refer to the solution vector, the longitudinal state and the disturbance, respectively. We substitute (2.6) into (2.1) and neglect the nonlinear terms, leading to the linearised equations,

(2.7a) \begin{equation} \frac {\partial \boldsymbol{\hat {q}}}{\partial t} = \mathcal{A}(\boldsymbol{q}_{\textit{LR}}; \textit{Re}, Ra, \textit{Pr}) \boldsymbol{\hat {q}}, \end{equation}

where

(2.7b) \begin{align} \mathcal{A}(\boldsymbol{q}_{\textit{LR}}; \textit{Re}, Ra, \textit{Pr}) = \left ( \begin{array}{@{}c c c} - (\boldsymbol{u}_{\textit{LR}} \boldsymbol{\cdot }\boldsymbol{\nabla }) - (\boldsymbol{\nabla }\boldsymbol{u}_{\textit{LR}} \; \boldsymbol{\cdot }) + 1/\textit{Re} {\nabla} ^2 & \frac {Ra}{8\textit{Re}^2 \textit{Pr}}\boldsymbol{\hat {j}} & -\boldsymbol{\nabla }\\[9pt] -(\boldsymbol{\nabla }\theta _{\textit{LR}} \; \boldsymbol{\cdot }) & - (\boldsymbol{u}_{\textit{LR}} \boldsymbol{\cdot }\boldsymbol{\nabla }) +{\nabla} ^2 & 0 \\[9pt] \boldsymbol{\nabla }\boldsymbol{\cdot }& 0 & 0 \end{array} \right )\!. \end{align}

Consider the longitudinal rolls invariant along the $x$ -direction; the following form of normal-mode solution can be proposed:

(2.8) \begin{equation} \boldsymbol{\hat {q}} (\boldsymbol{x},t) = \boldsymbol{\breve {q}}( y,z)e^{i(\alpha x+\beta z)+\lambda t} + \text{c.c}, \end{equation}

where $\lambda , \alpha$ and $\beta$ are the complex frequency, the streamwise wavenumber and the spanwise wavenumber (or the Floquet exponent), respectively, and c.c. stands for complex conjugate. Using the periodic nature of $\boldsymbol{\breve {q}}(y,z)$ in the $z$ -direction, (2.8) can also be written as

(2.9) \begin{equation} \boldsymbol{\hat {q}} (\boldsymbol{x},t) = \left [\sum _{n=-\infty }^{\infty }\boldsymbol{\breve {\breve {q}}}_n(y)e^{i \frac {2\pi }{ L_z}(n+\sigma ) z} \right ] e^{i \alpha x +\lambda t} + \text{c.c}, \end{equation}

where $\sigma (=\beta L_z/(2\pi ))$ is the Floquet detuning parameter with $0 \leqslant \sigma \leqslant 1/2$ . In this study, we will only consider unstable manifolds of longitudinal rolls in a fixed computational domain; therefore, the fundamental mode, $\sigma = 0$ , is of sole interest. Substituting (2.9) into (2.7) results to a discretised eigenvalue problem with the eigenvalue $\lambda$ . The wavenumber $\alpha$ is restricted to discrete values of $\alpha =2 \pi m/L_x$ , and $m$ is a positive integer for the given computational domain. The resulting eigenvalue problems are solved using an iterative Arnoldi algorithm based on the time stepper implemented in Nektar++ (Tuckerman & Barkley Reference Tuckerman and Barkley2000; Barkley, Blackburn & Sherwin Reference Barkley, Blackburn and Sherwin2008).

4.1. Thermally assisted sustaining process (TASP) in small domains, $\varGamma = 2\pi$

Given the complex spatio-temporal dynamics of the newly identified intermittent rolls and the co-existence of longitudinal rolls with turbulent bands (see Appendices C.3 and C.4, respectively), we consider simulations in small domains defined by $\varGamma = 2\pi$ , where longitudinal rolls and localised turbulence could be viewed as spatially isolated. This domain size is chosen such that the spanwise wavelength of the most unstable longitudinal rolls ( $\lambda _z\simeq \pi$ ) is comparable to the horizontal length of the domain. Simulations for a smaller domain ( $\varGamma =\pi$ ) do not show any qualitative difference from the simulation results reported here (Huang Reference Huang2025). We start from a numerical simulation at $Ra = 10\,000$ and $\textit{Re}=1050$ in $\varGamma = 2\pi$ , integrated in time for $t \in [0, 3000]$ . The initial condition has been sampled from a statistically stationary turbulent field at $\textit{Re} = 2000$ for $Ra=10\,000$ , and the $\textit{Re}$ is then slowly lowered to $\textit{Re} = 1050$ . The time history for $t \in [0, 3000]$ of the two near-wall transport properties, the Nusselt number and the shear rate on the lower wall is presented in figure 2, together with snapshots of the temperature, $\theta (\boldsymbol{x})$ , and the near-wall streamwise and spanwise perturbation velocities, $u'|_{y^+= 15}, v'|_{y^+=15}$ at selected times. In this small domain, the dynamics of the system exhibits temporal intermittency, where the solution trajectory appears to wander between longitudinal rolls and highly disorganised chaotic flow fields, characterised by low- and high-near-wall transport properties, respectively.

Figure 2.

Intermittent dynamics in a small domain at $Ra = 10000$ , $\textit{Re} = 1050$ , $t\in [0,3000]$ , $\varGamma = 2\pi$ . (a) The time history of the Nusselt number and shear. Temporal snapshots of volumetric temperature, planar near-wall streamwise and spanwise perturbations at (b) $t = 1291.5$ , (c) $t = 1480.5$ , (d) $t = 1564.5$ , (e) $t = 1711.5$ . Longitudinal rolls and transient turbulence are observed in panels (b,d) and (c,e), respectively.

Starting from a longitudinal roll state of spanwise wavenumber of $\alpha d = 4$ at $t = 1291.5$ in figure 2(b), the solution erupts into a highly disorganised turbulent state at $t = 1480.5$ , marked by a disordered temperature field in figure 2(c). During this breakdown, the near-wall snapshots of streamwise perturbation velocity, $u'|_{y^+ = 15}$ , and wall-normal perturbation velocity, $v'|_{y^+ = 15}$ , illustrated in the bottom panels of figure 2(c), reveal three pairs of high- and low-speed streaks, each with an average spanwise wavelength of $\varLambda _z^+ \approx 339/3 = 113$ (where $\varLambda _z^+ = u_\tau \varLambda _z / \nu$ refers to non-dimensionalised wavelength), close to the mean streak spacing ( $\varLambda ^+_z \sim 100$ ) commonly reported in shear flow turbulence (Kline et al. Reference Kline, Reynolds, Schraub and Runstadler1967; Smith & Metzler Reference Smith and Metzler1983; Kim, Moin & Moser Reference Kim, Moin and Moser1987; Hamilton et al. Reference Hamilton, Kim and Waleffe1995). These streaks appear to be meandering, negatively correlated with wall-normal perturbation velocities, reminiscent of a streak breakdown process (Hamilton et al. Reference Hamilton, Kim and Waleffe1995), or a bursting event (Kim, Kline & Reynolds Reference Kim, Kline and Reynolds1971), where high- and low-speed streaks are brought close to and away from the wall, respectively, enhancing near-wall transport quantities. Indeed, this is reflected in large increments in the Nusselt number and wall shear rate of roughly $40\,\%$ at $t = 1480.5$ in figure 2(a). Subsequently, the solution trajectory returns to a longitudinal roll state at $t=1564.5$ , before erupting into turbulence at $t = 1711.5$ (see figure 2 d,e, respectively). This suggests that the turbulence has a finite lifetime, occurring transiently before decaying towards the laminar state at $\textit{Re} = 1050$ (Hof et al. Reference Hof, Westerweel, Schneider and Eckhardt2006; Schneider, Eckhardt & Yorke Reference Schneider, Eckhardt and Yorke2007), which is linearly unstable, leading to the onset of longitudinal rolls where transient turbulence could be re-excited again.

Figure 3.

Relaminarisation in a small domain at $Ra = 0$ , $\textit{Re} = 1050$ , $t\in [0,3000]$ , $\varGamma = 2\pi$ . The time history of the (a) Nusselt number and shear. Temporal snapshots of volumetric temperature at (b) $t = 31.5$ , (c) $t = 63$ , (d) $t = 157.5$ , (e) $t = 672$ .

To test this hypothesis, we consider a numerical simulation at $Ra = 0$ , $\textit{Re} = 1050$ , in $\varGamma = 2\pi$ , where longitudinal rolls cannot appear. The initial condition is taken from a stationary turbulent solution at a higher $\textit{Re}$ of $\textit{Re} = 2000$ , which is then slowly lowered to $\textit{Re} = 1050$ and then integrated in time for $t \in [0, 700]$ . The time history of the Nusselt number, $\textit{Nu}$ , and the wall shear rate, $\langle \tau _w \rangle _{x,z}$ , is reported in figure 3, with the temperature snapshots, $\theta (\boldsymbol{x})$ , at selected time units. An initial turbulent flow field decays towards the laminar solution in $t \in [0, 700]$ . Comparing the results between $Ra = 0$ and $Ra = 10\,000$ , we propose that the longitudinal rolls at $Ra= 10\,000$ could mediate a transition mechanism between the unstable laminar state and transient turbulence, so that the entire non-trivial flow dynamics could be sustained indefinitely.

Next, we investigate the impact of longitudinal rolls on this proposed mechanism at different $Ra$ . We capture the flow quantities for $Ra=10\,000$ and $\textit{Re}=1050$ at $t=850.5$ (just before the onset of longitudinal rolls, see figure 2), and use as an initial condition to perform four numerical simulations by lowering $Ra$ instantly to $Ra=8000,5000,3000,2000$ , and time-integrated further to $t \in [850.5, 5000]$ . The time history of the wall shear rate, $\langle \tau _w\rangle _{x,z}$ , and the temperature snapshots, $\theta (\boldsymbol{x})$ , of these experiments of ‘ $Ra$ -quenching’ are presented in figure 4. The time history of shear is visibly intermittent for $Ra = 8000, 5000$ , depicted as the orange and green trajectories in figure 4(a), similar to $Ra = 10\,000$ . At $Ra= 8000, 5000$ , the longitudinal rolls emerge approximately at $t = 1312.5$ (see figure 4 c,e), before erupting into turbulence at $t = 1743$ and $t = 3570$ in figure 4(d, f), respectively. This is then accompanied by a large spike in the wall shear rate before dipping briefly, shown as the orange and green trajectories of figure 4(a). As $Ra$ is lowered further to $Ra = 3000, 2000$ , the transients begin to decay into a longitudinal state from $t = 850.5$ to $t = 1312.5$ , which remains stable until $t = 4200$ illustrated by figure 4(h, j), accompanied by the red and purple asymptotic trajectories in figure 4(a). This suggests that the longitudinal rolls are likely to be linearly unstable for $Ra = 8000, 5000$ , leading to turbulence, while remaining stable for $Ra = 3000,2000$ . In particular, the longitudinal roll state at $Ra = 5000$ remained saturated for a longer period $t \in [1500, 3400]$ (green curve in figure 4 a), indicating that the growth rate of the linear instability is smaller than for $Ra = 8000$ . We note that the longitudinal rolls in figure 4 have a spanwise wavenumber of $\alpha d=4$ , which corresponds to the wavenumber of the dominant primary instability (see Appendix E), indicating that it is the most preferred wavenumber within a small domain.

Figure 4.

$Ra$ -quenching experiments for $Ra = 8000, 5000, 3000, 2000$ , at $\textit{Re} = 1050$ , $\varGamma = 2\pi$ , $t \in [850.5, 5000]$ . The time history of (a) shear and (b) volumetric temperature snapshots of the flow condition at $t = 850.5$ . Volumetric temperature snapshots for $Ra = 8000$ at (c,d) $t = 1312.5, 1743$ , and $Ra = 5000$ at (e, f) $t = 1312.5, 3570$ , revealing a longitudinal roll and a turbulent state, respectively. Stable longitudinal rolls emerge for $Ra = 3000$ at (g,h) $t = 1312.5, 4200$ , and $Ra = 2000$ at ( j,k) $t = 1312.5, 4200$ .

Figure 5.

Growth rates of infinitesimal perturbations linearised about longitudinal rolls, $\boldsymbol{q}_{\textit{LR}}$ , of spanwise wavenumber of $\beta h = 2$ , against (a) streamwise wavenumber $\alpha h$ and (b) $Ra$ for $\alpha h =1$ . The hatches in panel (a) refer to wavenumbers smaller than those admissible in $\varGamma = 2\pi$ . The dash-dotted line in panel (b) is a standard quadratic regression yielding $Ra_{s}\approx 4720$ .

To understand the stability characteristics of the longitudinal rolls, we perform a linear stability analysis about the longitudinal roll state ( $\alpha d= 4$ ), at $Ra = 10\,000, 8000,{} 5000, 3000, 2000$ . The details of linear stability analysis are described in § 2.4, where $\lambda$ and $\boldsymbol{\breve {q}}(y,z)e^{i\alpha x}$ refer to the eigenvalue and eigenmode. The longitudinal roll (base) states, $\boldsymbol{q}_{\textit{LR}}$ , are obtained by time integrating an initial condition consisting of the laminar (conduction) state, superimposed by the primary eigenmode, $\beta h = 2$ , at $Ra = 10\,000, 8000, 5000, 3000, 2000$ , in a two-dimensional $y{-}z$ plane, suppressing any three-dimensional perturbations numerically. The growth rates as a function of discrete streamwise wavenumbers, $1 \leqslant \alpha h \leqslant 2.5$ , are presented in figure 5. We note that the admissible streamwise wavenumbers within $\varGamma = 2\pi$ are $\alpha h = m$ , where $m$ is a positive integer, $m = 1, 2, ..$ , and $\alpha h = 1.5, 2.5$ are included for completeness. The longitudinal rolls are linearly unstable for $Ra \geqslant 5000$ , while they remain stable for $Ra \leqslant 3000$ , which confirms our earlier hypothesis. In particular, the growth rates between $Ra = 5000$ and $Ra = 10\,000$ differ by an order of magnitude, which could explain the prolonged period of saturation in the green curve of figure 4(a). The dominant secondary instability of longitudinal rolls in $\varGamma = 2\pi$ has a streamwise wavenumber of $\alpha h = 1$ . Using a standard quadratic regression, the critical Rayleigh number for disturbances with $\alpha h = 1$ is approximately $Ra_{s} \approx 4720$ , presented in figure 5(b).

Figure 6.

Time integration along the dominant unstable manifold, $\alpha h=1$ , of the longitudinal rolls at $Ra = 5000, \textit{Re} = 1050,$ $\varGamma = 2\pi$ for $t\in [0,8000]$ . Time history of the (a) Nusselt number and wall shear rate, and (b) midplane temperature space–time plot. This system oscillates between the longitudinal rolls ( $\textit{LR1}{-}4$ ) and turbulence ( $T1{-}4$ ) over four intervals. Snapshots of temperature at (c) $t = 105$ , (d) $t = 1512$ , (e) $t = 2446.5$ , ( f) $t = 2992.5$ , (g) $t = 3727.5$ , (h) $t = 5103$ , (i) $t = 6300$ , (j) $t = 6993$ .

Following this, we examine the dominant unstable manifold ( $\alpha h = 1$ ) of the longitudinal rolls ( $\beta h = 2$ ) at $Ra=5000$ , by considering an initial condition,

(4.1) \begin{equation} \boldsymbol{q}_0(\boldsymbol{x},t=0) = \boldsymbol{q}_{\textit{LR}}(y,z) + \boldsymbol{\breve {q}}(y,z)e^{i\alpha x}+c.c. \end{equation}

Here, $\breve {\boldsymbol{q}} e^{i \alpha x}$ is an eigenmode for the streamwise wavenumber $\alpha$ and its amplitude was scaled such that its total energy is defined by

(4.2) \begin{equation} \delta = \frac {1}{{V}} \int _\varOmega \Bigg (\boldsymbol{\hat {u}}(\boldsymbol{x})^T\boldsymbol{\hat {u}}(\boldsymbol{x}) + \frac {Ra}{8\textit{Re}^2\textit{Pr}}\hat {\theta }(\boldsymbol{x})^2 \Bigg )\; \mathrm{d}\varOmega = 10^{-3}. \end{equation}

We have also considered that $\delta = 10^{-2}, 10^{-4}$ , but $\delta = 10^{-3}$ was found to be small enough to ensure linear growth, while large enough without requiring a large number of time steps.

The initial condition is time-integrated for $t \in [0, 8000]$ and its time history of the near-wall transport properties, the space–time plot of the mid-plane temperature, $\theta |_{(x,y)=(\pi ,0)}$ , are presented in figure 6, with snapshots of temperature, $\theta (\boldsymbol{x})$ , and near-wall streamwise and spanwise perturbation velocities, $u'|_{y^+= 15}$ , $v'|_{y^+=15}$ at some selected times. The intermittent trajectory is visually apparent, oscillating between longitudinal rolls and transient turbulence over four cycles for $t \in [0, 8000]$ : for example, the regions of low and high near-wall transport quantities in figure 6(a) correspond well to the organised and disorganised longitudinal rolls in figure 6(b), respectively. As the solution emerges from the unstable manifold of the longitudinal roll state, $(\textit{LR1})$ in figure 6(c), the trajectory erupts into turbulence at $t = 1512$ , marked by a disordered volumetric temperature field in snapshot $(T1)$ in figure 6(d). Turbulence occurs transiently and the solution decays towards a longitudinal roll-like state at $t = 2446.5$ , shown by snapshot $(\textit{LR2})$ in figure 6(e), forming a single cycle. The intermittent cycle repeats over three subsequent intervals, where the transient turbulent state and longitudinal rolls emerge at $t = 2992.5, 5103, 6993$ and $t = 3727.5, 6300$ , as shown in the snapshots $(T2,3,4)$ and $(\textit{LR3},4)$ in figures 6( f,h, j) and 6(g,i), respectively. Lastly, here, we show that the dominant unstable manifold of the longitudinal rolls is connected to a transient turbulent state. Interestingly, a longitudinal roll with $\beta h = 1$ sometimes emerges after turbulence decays, shown as the snapshot $(\textit{LR3})$ . This suggests that other unstable manifolds may well be linked to the transition to transient turbulence.

Figure 7.

State space projection based on the planar averaged centreline velocity and shear, coloured by the volume normalised perturbation kinetic energy at $Ra = 5000$ , $\textit{Re} = 1050$ , $\varGamma = 2\pi$ , (a) $t \in [0,800]$ , (b) $t \in [0, 68750]$ . The open black circles represent the unstable equilibria of longitudinal rolls and the laminar state. Note that the black-crosses, labelled by (T1–4) and (LR1–4), refer to temporal snapshots in figure 6, not equilibria solutions.

To better visualise the temporal dynamics in figure 6, we project the solution trajectory onto the space composed of two state observables: planar averaged centreline velocity, $\langle u|_{y=0} \rangle _{x,z}$ , the wall shear rate, $\langle \tau _w\rangle _{x,z}$ , coloured by the volume-averaged perturbation kinetic energy, $1/(2V)||\boldsymbol{u}'||^2$ , in figure 7, where $\|\boldsymbol{u}'\|^2=\int _\varOmega \boldsymbol{u}'^H\boldsymbol{u}' \,{\rm d}V$ with $\varOmega$ being the flow domain. These observables are chosen because they are found to distinguish well the regions of turbulent states, longitudinal roll states and the laminar state that reside around $(0.82, 3.2)$ , $(0.90, 2.32)$ and at $(1, 2)$ , respectively. Indeed, the temporal snapshots of $(T1{-}4)$ appear around the representative location of turbulent states, $(0.82, 3.2)$ , and the snapshots of $(\textit{LR1}{-}4)$ and an equilibrium related to the longitudinal roll state with the spanwise wavenumber $\beta h=2$ (denoted by $\boldsymbol{q}_{\textit{LR}, \beta h = 2}$ ) are seen around $(0.90, 2.32)$ .

The solution trajectory emerges from the unstable manifold of the longitudinal roll state, $\boldsymbol{q}_{\textit{LR}_{\beta h = 2}}$ , evolving towards turbulent states around $(0.85,3.2)$ , characterised by a high wall shear rate. At this $\textit{Re}$ , the turbulence is transient in small domains, occurring with a finite lifetime, eventually decaying towards the laminar state (Zammert & Eckhardt Reference Zammert and Eckhardt2015; Tuckerman et al. Reference Tuckerman, Chantry and Barkley2020; Paranjape et al. Reference Paranjape, Yalnız, Duguet, Budanur and Hof2023). As the solution trajectory approaches the laminar solution at $(1,2)$ , it abruptly reverses towards the longitudinal roll state near $(0.95, 2.15)$ , $(\textit{LR3})$ , due to the buoyancy-driven linear instability of the laminar base state in the RBP flow. Subsequently, the solution trajectory could depart along an unstable manifold of the longitudinal rolls again, forming a cyclic process.

To see if this cycle could be sustained indefinitely, we consider a longer time horizon, $t \in [0, 68750]$ , illustrated in figure 7(b). The solution trajectory wanders between the ‘cloud’ of transient turbulence at the top left corner (in red), and longitudinal roll and laminar states (in blue) in the bottom right, forming a basin of attraction between the unstable longitudinal rolls, transient turbulence and the unstable laminar base state. This basin of attraction is likely established above a critical $Ra$ as the longitudinal rolls become linearly unstable (i.e. $Ra \gtrsim Ra_{s} \approx 4720$ , see figure 5 b) and the instability of the longitudinal rolls provides an intermediate pathway towards transient turbulence, which could be regenerated again – a ‘self-sustaining’ dynamical process. We refer to this process as the thermally assisted sustaining process (TASP), inspired by the self-sustaining process (SSP) from turbulent shear flows (Hamilton et al. Reference Hamilton, Kim and Waleffe1995). A schematic of the TASP is illustrated in figure 8.

Figure 8.

Thermally assisted sustaining process (TASP).

4.2. Variation of $Ra$ and $\textit{Re}$ on the TASP in small domains, $\varGamma = 2\pi$

Figure 9.

Behaviour of the unstable and stable longitudinal rolls at $Ra = 8000, 4000$ for (a,e) $\textit{Re} = 600$ , (b, f) $\textit{Re} = 700$ , (c,g) $\textit{Re} = 1000$ and (d,h) $\textit{Re} = 1400$ within $\varGamma = 2\pi$ . Each parameter regime consist of three panels from the top to bottom, depicting the midplane temperature space–time plot, $\theta |_{(x,y) = ( \pi ,0)}$ , time history of the Nusselt number and shear, and state space projection based on the planar averaged centreline velocity and shear, coloured by the volume normalised perturbation kinetic energy.

In this section, we further explore the behaviour of the TASP as $\textit{Re}$ and $Ra$ are varied. We consider eight different cases at $Ra = 8000, 4000$ and $\textit{Re} = 600, 700, 1000, 1400$ . The results of these eight cases, where longitudinal rolls are unstable at $Ra = 8000$ or stable at $Ra = 4000$ , are shown in figure 9, depicting the space–time plots of the midplane temperature, $\theta |_{y=0}(x,t)$ , the time history of the Nusselt number, $\textit{Nu}$ , and the wall shear rate, $\langle \tau _w \rangle _{x,z}$ , and the state space portrait using the plane-averaged centreline velocity, $\langle u|_{y=0} \rangle _{x,z}$ , the wall shear rate. For all cases, except $Ra = 4000$ , $\textit{Re} = 1000$ and $\textit{Re} = 1400$ , their initial conditions are prepared from the laminar state, superimposed by a random noise based on a Gaussian distribution with zero mean and unit variance, scaled to a total energy of $\delta = 10^{-3}$ (see (4.2) for the definition). For the exceptional cases at $Ra = 4000$ , $\textit{Re} = 1000$ and $\textit{Re} = 1400$ , where subcritical turbulence and stable longitudinal rolls are expected, their initial conditions are obtained by gradually lowering $\textit{Re}$ from a statistically stationary turbulent solution at $Ra = 4000, \textit{Re} = 2000$ .

We first consider $\textit{Re}=1000$ (figure 9 c,g). At $Ra = 8000$ (figure 9 c), the trajectory visits transient turbulence near $t=7200$ , decaying towards the longitudinal roll state, $\boldsymbol{q}_{\textit{LR}}$ with roll wavenumber $\beta h =2$ , at $t = 7400$ , before regenerating again, consistent with the TASP in § 4.1. As $Ra$ is lowered to $4000$ (figure 9 g), the solution trajectory decays towards the longitudinal roll state, $\boldsymbol{q}_{\textit{LR}}$ . In this case, the longitudinal rolls are linearly stable, confirming our hypothesis earlier that the TASP is only established when longitudinal rolls become linearly unstable above a certain $Ra$ -threshold (i.e. $Ra \gtrsim Ra_{s} \approx 4720$ ).

At $\textit{Re} = 1400$ , $Ra = 4000$ , the solution trajectory remains within the turbulent ‘cloud’ near $(0.8, 3.8)$ for $t \in [0,10000]$ , as illustrated in figure 9(h). This suggests that turbulence in this case might be an attractor sustaining indefinitely, although we have not investigated whether the turbulent chaotic saddle at $\textit{Re} = 1000$ truly transitioned into an attractor at $\textit{Re} = 1400$ . As $Ra$ is increased to $8000$ , the solution trajectory originating from the laminar state, evolves towards the unstable longitudinal roll state, $\boldsymbol{q}_{\textit{LR}}$ at $t = 1550$ , transitioning into turbulence at $t= 1800$ . Therefore, in this case, the linearly unstable longitudinal rolls serve as an intermediate transitional pathway between the laminar base state and subcritical turbulence, whereas at $Ra = 4000$ , a bistability between stable longitudinal rolls (not shown) and turbulence is expected.

Next, we examine the behaviour of TASP as $\textit{Re}$ decreases towards the intermittent regime at $\textit{Re} = 600, 700$ , where a periodic orbit emerges between the longitudinal roll and the laminar state (figure 9 a,b,e, f). At $\textit{Re} = 600, Ra = 8000$ in figure 9(a), the solution trajectory initially evolves towards the longitudinal roll state, $\boldsymbol{q}_{\textit{LR}}$ , which is linearly unstable and breaks down towards the laminar state at $t = 2200$ . This breakdown is evidenced by the trajectory’s proximity to the laminar state in state space and the presence of a narrow green patch in the midplane temperature space–time plot. The longitudinal roll state is regenerated again, forming a periodic orbit with a period of $T_{\textit{period}} = 8098-6108 =1990$ , oscillating between the longitudinal roll and laminar state over five intervals within $t \in [0, 10\,000]$ . As $\textit{Re}$ increases slightly to $700$ , the periodic orbit persists over a shorter period of $T_{\textit{period}} = 3889 - 3386 = 503$ . A notable difference is observed in the regenerated longitudinal rolls, which is continuously translated by $L_z/2$ in the spanwise $z$ -direction. Additionally, as $\textit{Re}$ increases from $600$ to $700$ , the trajectory moves further away from the laminar state during breakdown, suggesting an increasing attraction towards the longitudinal roll state, $\boldsymbol{q}_{\textit{LR}}$ (compare $t = 2200$ in figure 9 a and $t = 2750$ in figure 9 b). When $Ra$ is lowered to $Ra = 4000$ , the periodic orbit disappears and the trajectory stabilises into the longitudinal roll state, $\boldsymbol{q}_{\textit{LR}}$ , at $\textit{Re} = 600, 700$ (figure 9 e, f).

Figure 10.

A state space sketch of figure 9 at $Ra = 8000$ , (a) $\textit{Re} = 600$ , (b) $\textit{Re} = 700$ , (c) $\textit{Re} = 1000$ , (d) $\textit{Re} = 1400$ and $Ra = 4000$ at (e) $\textit{Re} = 600,700$ , ( f) $\textit{Re} = 1000$ , (g) $\textit{Re} =1400$ . The longitudinal roll is linearly unstable (saddle) at $Ra = 8000$ and is stable at $Ra = 4000$ , whereas the laminar state is always linearly unstable (saddle). The blue and orange solid arrows refer to the unstable manifold of longitudinal rolls and the laminar state. The red solid lines denote the chaotic trajectories of turbulence, likely forming a chaotic saddle at $\textit{Re} = 1000$ and a chaotic attractor at $\textit{Re} = 1400$ . The black-dashed trajectories refer to possible solution trajectories, forming a periodic orbit (P.O) at $Ra = 8000$ , $\textit{Re} = 600, 700$ , and a basin of attraction (B.o.A) at $Ra = 8000, \textit{Re} = 1000$ . We note that invariant states could exist at $Ra = 4000, \textit{Re} = 600,700$ , labelled as a saddle here (Paranjape et al. Reference Paranjape, Yalnız, Duguet, Budanur and Hof2023).

To summarise the dynamical processes identified in figure 9, we present their state space sketch in figure 10. At $Ra = 8000$ , $\textit{Re} = 600$ and $\textit{Re}= 700$ , the longitudinal rolls become linearly unstable, breaking down to the laminar state before being regenerated, forming a periodic orbit enclosed by black dotted paths in figure 10(a,b). For $\textit{Re} = 700$ , the regenerated longitudinal roll is continuously translated by $L_z/2$ , suggesting a possible merger of two periodic orbits into one as sketched in figure 9(b). Future bifurcation studies are required to establish this, providing an avenue for future work. As $Ra$ is lowered to $Ra = 4000$ , the laminar state stabilises into the longitudinal rolls in figure 10(e). However, the flow in this regime may contain some invariant solutions (e.g. see Paranjape et al. Reference Paranjape, Yalnız, Duguet, Budanur and Hof2023), denoted as saddle points here. Integrating along the unstable manifold of longitudinal roll states at $Ra= 8000, \textit{Re} = 1000$ leads to transient turbulence, which eventually decays to the laminar state before regenerating the longitudinal rolls again, forming the TASP in figure 10(c). In contrast, at $Ra = 4000, \textit{Re} = 1000$ , the longitudinal rolls become linearly stable, eliminating the intermediate (orange) pathway towards turbulence. Therefore, the transient turbulence stabilises into longitudinal rolls, as shown by the black-dashed trajectory in figure 10( f). For $Ra = 8000, \textit{Re} = 1400$ , the linearly unstable longitudinal rolls provide an intermediate pathway towards turbulence from the laminar state, as sketched in figure 10(d), breaking the bistability between the laminar state and subcritical turbulence seen at $Ra = 4000$ in figure 10(g). This difference highlights the contribution of unstable longitudinal rolls towards the transition to turbulence within $\varGamma = 2\pi$ .

We have examined the dynamics of unstable longitudinal rolls, as the Reynolds number, $\textit{Re}$ , and the Rayleigh number, $Ra$ , are varied, identifying three key dynamical processes: (1) periodic orbits between longitudinal rolls and the laminar state (figure 10 a,b); (2) the TASP, where transient turbulence can be sustained (figure 10 c); and (3) an intermediate transitional pathway towards sustained turbulence (figure 10 d). To establish a connection between these processes and refine their transitional boundaries, we conduct a further parametric study over fifteen numerical simulations at $Ra = 4000, 6000, 10\,000$ and $\textit{Re} = 600, 800, 900, 1200, 1400$ for $\varGamma = 2\pi$ . Figure 11 presents the midplane temperature space–time plot along with the time history of the wall shear rate, $\langle \tau _w\rangle _{x,z}$ , and the Nusselt number, $\textit{Nu}$ . For all simulations, the initial conditions are prepared from the laminar state, superimposed with a random noise based on a Gaussian distribution with zero mean and unit variance, scaled to a total energy of $\delta = 10^{-3}$ (see definition in (4.2)). Due to the subcritical nature of turbulence and expected stable longitudinal rolls, exceptions are made for $Ra = 4000$ , $\textit{Re} = 900, 1200, 1400$ , where initial conditions are taken from gradually lowering $\textit{Re}$ from a statistically stationary turbulent state at $Ra = 4000, \textit{Re} = 2000$ . The thermally assisted sustaining process is highlighted in green for $Ra \in [6000, 10\,000]$ and $\textit{Re} \in [900, 1200]$ , where temporally intermittent shear and Nusselt number fluctuations are observed, accompanied by a mixture of organised and disorganised flow structures in the temperature space–time plots. Here, longitudinal rolls provide an intermediate pathway towards transient turbulence. For the same $\textit{Re}$ and as $Ra$ lowers to $4000$ , transient turbulence eventually decays into stable longitudinal rolls labelled as ‘transient turbulence’ in figure 11. Periodic and quasi-periodic orbits between longitudinal rolls and the laminar state appear for $Ra = 10000$ , $\textit{Re} = 600, 800$ and $Ra = 6000$ , $\textit{Re} = 800$ , establishing threshold values of $Ra$ and $\textit{Re}$ shaded in yellow. Below this threshold, solutions stabilise into longitudinal rolls shaded in red. Although longitudinal rolls are linearly stable at $Ra = 4000, \textit{Re} = 1400$ , turbulence is sustained at least for a sufficiently long time, shaded blue at $\textit{Re} = 1400$ and labelled as ‘sustained turbulence’ in figure 11. In this case, a bistable system forms between longitudinal rolls and turbulence at $Ra = 4000$ , while longitudinal rolls provide an intermediate pathway to turbulence for $Ra \geqslant 6000$ .

Finally, we aim to establish the relevance of the local dynamical processes, such as the TASP identified in the small domain simulations of this section, to the results obtained in large domain simulations (see § 3), as discussed in Appendix D.

Figure 11.

Temperature space–time plots and time history of $(\tau _w, Nu)$ , for $Ra \in [5000, 10\,000]$ , $\textit{Re} \in [600, 1400]$ within $\varGamma = 2\pi$ . Unstable longitudinal rolls lead to the onset of (1) periodic orbits (yellow), (2) the thermally assisted sustaining process (green) and (3) sustained turbulence (blue), occurring beyond an $Ra{-}\textit{Re}$ boundary, below which longitudinal rolls remain stable (red).