Fixed-flux Rayleigh–Bénard convection in doubly periodic domains: generation of large-scale shear

Abstract

This work studies two-dimensional fixed-flux Rayleigh–Bénard convection with periodic boundary conditions in both horizontal and vertical directions and analyses its dynamics using numerical continuation, secondary instability analysis and direct numerical simulation. The fixed-flux constraint leads to time-independent elevator modes with a well-defined amplitude. Secondary instability of these modes leads to tilted elevator modes accompanied by horizontal shear flow. For $Pr=1$, where $Pr$ is the Prandtl number, a subsequent subcritical Hopf bifurcation leads to hysteresis behaviour between this state and a time-dependent direction-reversing state, followed by a global bifurcation leading to modulated travelling waves without flow reversal. Single-mode equations reproduce this moderate Rayleigh number behaviour well. At high Rayleigh numbers, chaotic behaviour dominated by modulated travelling waves appears. These transitions are characteristic of high wavenumber elevator modes since the vertical wavenumber of the secondary instability is linearly proportional to the horizontal wavenumber of the elevator mode. At a low $Pr$, relaxation oscillations between the conduction state and the elevator mode appear, followed by quasi-periodic and chaotic behaviour as the Rayleigh number increases. In the high $Pr$ regime, the large-scale shear weakens, and the flow shows bursting behaviour that can lead to significantly increased heat transport or even intermittent stable stratification.

1. Introduction

Fixed-flux temperature boundary conditions describing adiabatic boundaries appear in a wide range of geophysical and astrophysical applications, including convection in the Earth's mantle (McKenzie, Roberts & Weiss 1974; Chapman, Childress & Proctor 1980; Hewitt, McKenzie & Weiss 1980; Sakuraba & Roberts 2009; Long et al. 2020) and models of solar supergranulation (Vieweg, Scheel & Schumacher 2021; Vieweg et al. 2022; Käufer et al. 2023). Such boundary conditions also model Rayleigh–Bénard convection (RBC) between poorly conducting horizontal plates, suggesting an explanation for discrepancies between experimentally measured heat transport and numerical simulations (Verzicco 2004; Verzicco & Sreenivasan 2008). Fixed-flux temperature or mass boundary conditions are also used to model the free surface in ocean circulation models (Huck, de Verdière & Weaver 1999; Abernathey, Marshall & Ferreira 2011) and for understanding nutrient productivity in the oceans (Pasquero, Bracco & Provenzale 2005). A low enough constant injection rate of ${\rm CO}_2$ concentration in saline aquifers can also be modelled as a fixed-flux problem (Amooie, Soltanian & Moortgat 2018). Moreover, fixed-flux boundary conditions are also relevant for concentration transport within an enclosure with impermeable boundaries (Mamou, Vasseur & Bilgen 1998; Mamou & Vasseur 1999).

Rayleigh–Bénard convection provides a canonical set-up for understanding the effect of fixed-flux conditions; see e.g. the monograph by Goluskin (2016). With fixed flux at both the top and bottom boundaries, the critical horizontal wavenumber at the onset of convective instability vanishes (Sparrow, Goldstein & Jonsson 1964; Hurle, Jakeman & Pike 1967; Chapman & Proctor 1980), and the resulting near-onset evolution is described by the Cahn–Hilliard equation (Novick-Cohen 2008; Miranville 2019). In the weakly supercritical Rayleigh number regime, each convection cell is thus unstable to perturbations with ever longer wavelength in a process referred to as coarsening (Chapman & Proctor 1980; Chapman et al. 1980). This large scale manifests itself in the turbulent regime of three-dimensional (3-D) fixed-flux RBC, and organizes the resulting flow. For example, recent 3-D direct numerical simulations (DNS) with aspect ratio $\varGamma = 60$ show that convection cells aggregate gradually into a single large cell that eventually fills the whole domain, thereby providing insight into the aggregation of granules into a supergranule in the solar convection zone (Vieweg et al. 2021, 2022), a process also confirmed in experiments (Käufer et al. 2023). There is also evidence that fixed-flux boundary conditions influence the generation and reversals of large-scale shear. For example, experimental studies found that a configuration with constant flux at the bottom and constant temperature at the top exhibits less frequent reversals of the large-scale circulation than in a configuration with constant temperature on both surfaces (Huang et al. 2015).

Heat transport in fixed-flux RBC has also been analysed widely. In two-dimensional (2-D) domains of modest aspect ratio, fixed-flux and fixed-temperature RBC display essentially identical heat transport, overall flow dynamics and mean temperature profiles at Rayleigh number $Ra_T=10^{10}$ based on temperature difference (Johnston & Doering 2009), despite possible differences in the dominant scale. In 3-D cylindrical geometry, Stevens, Lohse & Verzicco (2011) investigated the difference in heat transport introduced by replacing the bottom plate by a fixed-flux condition, and showed that this difference decreases with increasing Rayleigh number. Different heat transport scaling predictions may be realized depending on the details of the thermal forcing (Lepot, Aumaître & Gallet 2018; Bouillaut et al. 2019), as shown in experiments using radiative heating in a thermally insulating container to control the heat flux. However, the current upper bound on the Nusselt number with fixed-flux boundary conditions displays the same scaling law with the Rayleigh number based on temperature difference as the fixed-temperature configuration with either no-slip (Otero et al. 2002; Wittenberg 2010) or stress-free velocity boundary conditions (Fantuzzi 2018).

Fixed-flux temperature boundary conditions are also studied in related set-ups. For example, the difference between Neumann and Dirichlet boundary conditions on the buoyancy field in moist convection decreases as the Rayleigh number increases (Weidauer & Schumacher 2012), while fixed-heat-flux and fixed-temperature boundary conditions are shown to be asymptotically equivalent in rapidly rotating convection (Calkins et al. 2015). For rotating convection with no-slip boundaries, the Nusselt number increases significantly when a fixed heat flux is imposed instead of a fixed temperature difference (Kolhey, Stellmach & Heyner 2022).

The top and bottom boundaries are often absent in geophysical applications, suggesting that periodic boundary conditions in the vertical are more appropriate. The resulting homogeneous RBC problem driven by a constant temperature gradient has been employed to understand bulk RBC (Borue & Orszag 1997; Lohse & Toschi 2003; Calzavarini et al. 2005, 2006; Ng et al. 2018; Pratt, Busse & Müller 2020; Barral & Dubrulle 2023). Similar homogeneous configurations are also commonly employed to analyse double-diffusive convection (Stellmach et al. 2011; Radko 2013; Garaud 2018) and different shear flows (Rogers & Moin 1987; Sekimoto, Dong & Jiménez 2016). Periodic boundary conditions in the vertical within these homogeneous set-ups have the benefit of eliminating inessential but computationally expensive thermal or viscous boundary layers. Within homogeneous RBC, the Nusselt number $Nu$ scales with the Rayleigh number $Ra$ according to the ultimate regime prediction (Lohse & Toschi 2003; Calzavarini et al. 2005), a prediction supported by experimental evidence from a cylindrical cell (Schmidt et al. 2012). Moreover, Calzavarini et al. (2005) showed that $Nu$ scales with the Prandtl number $Pr$ according to mixing length theory (Spiegel 1963), and attributed this fact to more frequent appearances of elevator modes at high $Pr$. An exponentially growing elevator mode is an exact nonlinear solution of the homogeneous RBC problem, whose growth in DNS is arrested only by secondary numerical noise with a resolution-dependent instability ultimately leading to statistically steady transport (Calzavarini et al. 2006). The appearance of elevator modes also leads to high intermittency in the heat transport (Borue & Orszag 1997; Calzavarini et al. 2005, 2006; Barral & Dubrulle 2023), thereby affecting the flow statistics adversely, by increasing the sensitivity to round-off noise and discretization error (Calzavarini et al. 2006). In DNS, these polluting elevator modes can be suppressed by removing explicitly the mean flow parallel to gravity at each time step (Pratt et al. 2020), or by introducing an artificial horizontal buoyancy field (Xie & Huang 2022) or large-scale friction (Barral & Dubrulle 2023), but the modes remain a major source of contention.

This work formulates a fixed-flux homogeneous RBC problem that is not only relevant to a wide range of geophysical and astrophysical applications but also avoids the exponentially growing elevator modes that plague homogeneous RBC driven by a constant temperature gradient. Our study is motivated by a recent formulation imposing fixed-flux salinity conditions on homogeneous inertia-free salt-finger convection (IFSC) (Xie, Julien & Knobloch 2020). This fixed-flux constraint saturates the elevator mode in IFSC, and it does so in the present case as well. In both cases, the resulting formulation leads to differential–integral equations with time-varying mean salinity or temperature gradients that adjust the system response. Moreover, fixed-flux conditions result in a more potent restoring mechanism towards the statistically steady state that can be used as a diagnostic to determine whether in situ salt-finger convection is flux-driven or gradient-driven (Xie et al. 2020).

This work thus focuses on the underlying dynamics of fixed-flux homogeneous RBC using numerical continuation, secondary instability analysis, and DNS. Secondary instabilities of the elevator mode lead to tilted elevator modes accompanied by horizontal jet formation. At $Pr=1$, this state is in turn unstable to a subcritical Hopf bifurcation displaying hysteresis behaviour between this state and the resulting direction-reversing state. A subsequent global bifurcation of Shilnikov type (Shilnikov & Shilnikov 2007) leads to modulated travelling waves without flow reversal. Single-mode equations that severely truncate the horizontal flow structure reproduce this moderate Rayleigh number behaviour rather well.

At high Rayleigh numbers, chaotic flow dominated by modulated travelling waves appears, and resembles no-slip instead of stress-free boundary conditions in RBC with fixed temperature. The vertical wavenumber of the secondary instability of steady elevator modes leading to these transitions is linearly proportional to the horizontal wavenumber of the elevator mode, leading to its suppression when the vertical extent of the domain precludes its presence. Thus the domain aspect ratio requires adjustment as the parameters are varied.

At low Prandtl numbers, relaxation oscillations between the conduction state and the elevator mode appear, followed by quasi-periodic and chaotic behaviour as the Rayleigh number increases. Since the secondary and Hopf bifurcation points shift closer to the primary instability as $Pr$ decreases, the single-mode description works well in this regime. In contrast, at high $Pr$ the large-scale shear weakens, and the flow exhibits bursting behaviour resembling that in homogeneous RBC driven by a constant temperature gradient (Borue & Orszag 1997; Calzavarini et al. 2005, 2006) and resulting in significantly increased heat transport or even intermittent stable stratification.

The remainder of this paper is organized as follows. Section 2 formulates the fixed-flux homogeneous RBC problem. Bifurcation analysis at moderate Rayleigh numbers performed via numerical continuation is described in § 3 and confirmed by DNS. Section 4 analyses the high Rayleigh number dynamics arising from the secondary instability of the elevator modes. The effects of changing the Prandtl number are discussed in § 5. The paper concludes with a summary and suggestions for future work in § 6.

2. Fixed-flux homogeneous RBC

We consider a layer of fluid of depth $h$ with a constant upward heat flux $-k q$ through it, where $k$ is thermal conductivity, and $q<0$ is the associated vertical temperature gradient. The equation of state $(\rho _*-\rho _{r*})/\rho _{r*}=-\alpha (T_*-T_{r*})$ is linear in the Boussinesq approximation, with constant expansion coefficient $\alpha$, constant reference density $\rho _{r*}$, and constant reference temperature $T_{r*}$. The subscript $*$ denotes a dimensional variable. In the following, we non-dimensionalize the temperature $T_*$ by the temperature gradient $|q|$ associated with imposed constant heat flux, $\underline {T}=T_*/|hq|$. Spatial variables are normalized by the depth $h$ of the layer, while time and velocity are normalized using the thermal diffusion time $h^2/\kappa _T$ and the corresponding speed $\kappa _T/h$, respectively. Here, $\kappa _T=k/(\rho _{r*} c_p)$ is the thermal diffusivity, with $c_p$ denoting specific heat capacity. In homogeneous double-diffusive convection, lengths are usually normalized by the expected finger width (Stellmach et al. 2011; Radko 2013), while here we normalize lengths by the layer depth $h$ for consistency with the usual procedure for RBC.

We introduce the velocity field $\boldsymbol {u}:=(u,v,w)$ in Cartesian coordinates $(x,y,z)$ with $z$ in the upward vertical direction. Under the Boussinesq approximation, the system is governed by

(2.1a)$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u}={-}\boldsymbol{\nabla} p+Pr\,Ra_{T,q}\underline{T}\boldsymbol{e}_z+ Pr\,\nabla^2 \boldsymbol{u}, \end{gather}$$

(2.1b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}$$

(2.1c)$$\begin{gather}\frac{\partial \underline{T}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \underline{T}=\nabla^2 \underline{T}, \end{gather}$$

where $\boldsymbol {e}_z$ is the unit vector in the vertical. The governing parameters are the flux Rayleigh number $Ra_{T,q}$ and the Prandtl number $Pr$:

(2.2a,b)\begin{equation} Ra_{T,q}:=\frac{\alpha g\,|q|\,h^4}{\nu \kappa_T},\quad Pr:=\frac{\nu}{\kappa_T}. \end{equation}

A similar flux Rayleigh number is also employed by Otero et al. (2002), Johnston & Doering (2007, 2009), Verzicco & Sreenivasan (2008) and Goluskin (2016).

We decompose the total temperature $\underline {T}(x,y,z,t)$ as

(2.3)\begin{equation} \underline{T}(x,y,z,t)=1+\bar{\mathcal{T}}_{z,q}z+T(x,y,z,t), \end{equation}

where $\bar {\mathcal {T}}_{z,q}$ is a spatially and temporally averaged temperature gradient. This decomposition allows us to impose vertically periodic boundary conditions on $T$. The velocity is taken as periodic in the vertical, and periodic conditions in the horizontal are imposed on all variables. We then define the volume-averaged heat flux

(2.4a)\begin{align} Q(t)&:=\langle (\boldsymbol{u}\underline{T}-\boldsymbol{\nabla}\underline{T})\boldsymbol{\cdot} \boldsymbol{e}_z\rangle_{h,v} \end{align}

(2.4b)\begin{align} &\,=\langle wT\rangle_{h,v}-\bar{\mathcal{T}}_{z,q}, \end{align}

where $\langle \cdot \rangle _{h,v}$ is the horizontal and vertical average. The equality in (2.4b) is obtained on assuming a vanishing homogeneous mode, $\langle w\rangle _{h,v}(t)=\langle T\rangle _{h,v}(t)=0$. We further assume that the instantaneous heat flux $Q(t)$ recovers the imposed value $Q_c$ exponentially rapidly at a rate $\beta$:

(2.5)\begin{equation} \frac{{\rm d}Q}{{\rm d}t}+\beta(Q-Q_c)=0. \end{equation}

Here, $Q_c=1$ because the temperature is normalized based on the imposed heat flux.

We can now write (2.1c) in terms of $T$ that is periodic in all spatial directions. This is obtained by substituting the decomposition (2.3) and (2.4) into (2.1c):

(2.6)\begin{equation} \frac{\partial T}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}T-Qw + w\langle wT\rangle_{h,v}=\nabla^2 T,\end{equation}

where $Q(t)$ is governed by (2.5). Note that (2.5) is taken to be independent of $(\boldsymbol {u},T,p)$. Thus setting $Q(t=0)=Q_c$ leads to $Q(t)=Q_c=1$. This corresponds to the $\beta \to \infty$ limit in which $Q(t)$ recovers the reference value $Q_c$ instantaneously. We show, moreover, that for $\beta =10^4$ and a random $Q(t=0)$, the results display the same behaviour as those for ${\beta =\infty}$ (see Appendix). Setting $Q(t)=Q_c=1$ and eliminating the hydrostatic pressure, we obtain the governing equations in the form

(2.7a)$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u}={-}\boldsymbol{\nabla} p+Pr\,Ra_{T,q}\,T\boldsymbol{e}_z+Pr\, \nabla^2 \boldsymbol{u}, \end{gather}$$

(2.7b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}$$

(2.7c)$$\begin{gather}\frac{\partial T}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}T-w + w \langle wT\rangle_{h,v}=\nabla^2 T. \end{gather}$$

The integral term $w\langle wT\rangle _{h,v}$ in (2.7c) is the new flux feedback term that does not appear in earlier formulations of the homogeneous RBC problem driven by a constant temperature gradient (Borue & Orszag 1997; Lohse & Toschi 2003; Calzavarini et al. 2005, 2006; Ng et al. 2018; Pratt et al. 2020; Xie & Huang 2022; Barral & Dubrulle 2023).

The response parameter is the instantaneous Nusselt number, which measures the ratio of the total convective transport to the conductive heat transport in the vertical:

(2.8)\begin{equation} nu(t):=\frac{ \langle (\boldsymbol{u}\underline{T}-\boldsymbol{\nabla}\underline{T})\boldsymbol{\cdot} \boldsymbol{e}_z\rangle_{h,v}}{\langle(-\boldsymbol{\nabla}\underline{T})\boldsymbol{\cdot} \boldsymbol{e}_z\rangle_{h,v}}=\frac{1}{1-\langle wT\rangle_{h,v}(t)}. \end{equation}

We also define a Nusselt number measuring the time-averaged heat transport as

(2.9)\begin{equation} Nu:=\frac{1}{1-\langle wT\rangle_{h,v,t}}, \end{equation}

where $\langle \cdot \rangle _{h,v,t}$ denotes spatio-temporal averaging. We can also obtain the mean temperature gradient by time-averaging (2.4):

(2.10)\begin{align} \bar{\mathcal{T}}_{z,q}&=\langle wT\rangle_{h,v,t}-\langle Q\rangle_t \end{align}

(2.11)\begin{align} &=\langle wT\rangle_{h,v,t}-1, \end{align}

where we assume $\langle Q\rangle _t=Q_c=1$, as appropriate for long-time averages. The Rayleigh number based on the mean temperature gradient is

(2.12)\begin{equation} Ra_T:=\frac{\alpha g\,|\bar{\mathcal{T}}_{z,q*}|\,h^4}{\nu \kappa_T}=\frac{\alpha g\,|q|\,h^4}{\nu \kappa_T}\,(-\bar{\mathcal{T}}_{z,q})=\frac{Ra_{T,q}}{Nu}, \end{equation}

where $\bar {\mathcal {T}}_{z,q*}=q\bar {\mathcal {T}}_{z,q}$ is the dimensional mean temperature gradient. A relation similar to (2.10) was also noted in RBC with a fixed imposed flux (Otero et al. 2002; Johnston & Doering 2007, 2009; Verzicco & Sreenivasan 2008; Goluskin 2016). As a result, a scaling law $Nu\sim Ra_{T,q}^{1/3}$ based on imposed flux corresponds to $Nu\sim Ra_T^{1/2}$ based on the mean temperature gradient, a relation similar to that between fixed-flux and fixed-temperature RBC (Otero et al. 2002).

3. Bifurcation analysis at moderate Rayleigh number

In this section, we analyse flow structures originating from the primary instability at moderate Rayleigh numbers and their subsequent destabilization by means of analytical calculation and numerical continuation as well as DNS. The nonlinear solutions and their stability determined from this analysis provide the pathway towards chaotic behaviour or even fully developed turbulent states, which generally visit neighbourhoods of (unstable) steady, periodic or travelling wave solutions, and these visits leave an imprint on the flow statistics; see e.g. Kawahara & Kida (2001), van Veen, Kida & Kawahara (2006) and the reviews by Kawahara, Uhlmann & Van Veen (2012) and Graham & Floryan (2021). We keep our analytical calculation general as appropriate for three dimensions, although our numerical results are confined to 2-D $(x,z)$ configurations.

3.1. Primary instability and the steady elevator mode

We start from the primary instability and the steady elevator mode that originates from this instability, which allows analytical progress. We linearize (2.7) around the conduction base state ($\boldsymbol {u}=\boldsymbol {0}$, $T=0$) by dropping the nonlinear terms. After eliminating the pressure in the vertical momentum equation by applying $-\boldsymbol {e}_z\boldsymbol {\cdot }\boldsymbol {\nabla }\times [\boldsymbol {\nabla }\times (\cdot )]$ to the momentum equation, we obtain

(3.1a)$$\begin{gather} \frac{\partial \nabla^2 w}{\partial t}=Pr\,\nabla^4 w+Pr\,Ra_{T,q}\,\nabla^2_\perp T, \end{gather}$$

(3.1b)$$\begin{gather}\frac{\partial T}{\partial t}=w+\nabla^2 T, \end{gather}$$

where $\nabla ^2_\perp :=\partial _x^2+\partial _y^2$. We use the normal mode assumption $\phi (x,y,z,t)=\hat {\phi }\exp [\text {i}(k_x x+k_y y+k_z z) +\lambda t]+\text {c.c.}$, where $\phi =w,T$, and $k_x$, $k_y$ and $k_z$ are the wavenumbers in the corresponding directions, and $\lambda$ is the (necessarily real) growth rate. Here, $\text {i}$ is the imaginary unit, and c.c. denotes the complex conjugate. This normal mode assumption yields

(3.2a)$$\begin{gather} \lambda \hat{w}=-Pr\,K^2\hat{w}+Pr\,Ra_{T,q}\,\frac{k_\perp^2}{K^2}\,\hat{T}, \end{gather}$$

(3.2b)$$\begin{gather}\lambda \hat{T}=\hat{w}-K^2\hat{T}, \end{gather}$$

where $K^2:=k_x^2+k_y^2+k_z^2$, and $k_\perp ^2:=k_x^2+k_y^2$. Solving this eigenvalue problem gives growth rate

(3.3)\begin{equation} \lambda={-}\frac{1}{2}\,(Pr+1)K^2\pm \frac{1}{2}\sqrt{(Pr+1)^2K^4+4\,Pr\left(Ra_{T,q}\,\frac{k_\perp^2}{K^2}-K^4\right)}. \end{equation}

When $k_z=0$, this growth rate is the same as that associated with the elevator mode in homogeneous RBC driven by a constant temperature gradient (Calzavarini et al. 2006, Eq. (9)).

The growth rate $\lambda$ vanishes at the onset of a steady bifurcation, leading to the neutral curve $Ra_{T,q}={K^6}/{k_\perp ^2}$. For a given $Ra_{T,q}$, the most unstable mode corresponds to $k_z=0$, i.e. to an elevator mode. The corresponding neutral curve then simplifies:

(3.4)\begin{equation} Ra_{T,q}=k_\perp^4, \end{equation}

again as for the case of constant temperature gradient forcing (Calzavarini et al. 2006). As a result, the critical horizontal wavenumber is $k_{\perp,c}=0$, as in RBC with fixed-flux boundary conditions; see e.g. Sparrow et al. (1964) and Chapman & Proctor (1980).

The resulting steady elevator mode ($k_z=0$) plays an important role in the subsequent behaviour of the system. The amplitude of the elevator mode is obtained by substituting

(3.5a)$$\begin{gather} w(x,y,z,t)=\hat{w}_e \exp[\text{i}(k_x x+k_y y)]+\text{c.c.}, \end{gather}$$

(3.5b)$$\begin{gather}T(x,y,z,t)=\hat{T}_e\exp[\text{i}(k_x x+k_y y)]+\text{c.c.} \end{gather}$$

into (2.7), which gives

(3.6a,b)\begin{equation} \hat{w}_e=\sqrt{\frac{Ra_{T,q}}{2k_\perp^2}-\frac{k_\perp^2}{2}}, \quad \hat{T}_e=\frac{k_\perp^2 \hat{w}_e}{Ra_{T,q}}. \end{equation}

Note that this is an exact solution of the nonlinear governing equation (2.7) and corresponds to the Nusselt number

(3.7)\begin{equation} Nu=\frac{Ra_{T,q}}{k_\perp^4}. \end{equation}

The steady elevator mode within the fixed-flux formulation thus has a unique time-independent amplitude (3.6a,b) for each Rayleigh number. In contrast, within homogeneous RBC driven by a constant temperature gradient, the steady elevator mode bifurcates from $Ra_{T,q}^{(p)}$ with arbitrary amplitude, and grows exponentially for $Ra_{T,q}>Ra_{T,q}^{(p)}$, leading to intermittent heat transport in DNS (Borue & Orszag 1997; Lohse & Toschi 2003; Calzavarini et al. 2005, 2006).

3.2. Numerical methods

For solution branches beyond the steady elevator mode, numerical computations are required. We compute each solution branch and associated bifurcation points by numerical continuation using pde2path (Uecker, Wetzel & Rademacher 2014; Uecker 2021a) with horizontal and vertical directions discretized by the Fourier collocation method (Weideman & Reddy 2000) following the implementation in Uecker (2021b). We use a streamfunction formulation of the full 2-D equations in (2.7) to reduce the number of variables, thereby facilitating computation. The horizontal and vertical directions use $N_x=N_z=32$ grid points, and doubling the number of grid points in each direction does not influence the results. The tolerance of the maximum absolute value of the residual at each vertical location ($L_\infty$ norm) is set to $10^{-6}$. We implement the phase condition associated with horizontal translation symmetry for elevator modes, and the phase conditions corresponding to both horizontal and vertical translations for all other 2-D solution branches (Rademacher & Uecker 2017). The stability of each branch is determined by computing a subset of the eigenvalues, and this subset is enlarged as necessary to ensure that instability and bifurcation points are identified correctly.

We also analyse time-dependent states through DNS in Dedalus (Burns et al. 2020) using a Fourier spectral method in both horizontal and vertical directions. We set the spatial homogeneous mode associated with $k_x=k_z=0$ to zero, which can be implemented by adding the constraint $\langle T\rangle _{h,v}(t)=0$. The quantity $\langle w\rangle _{h,v}(t)$ is conserved over time, and all of the results here are for $\langle w\rangle _{h,v}(t)=0$. A non-zero $\langle w\rangle _{h,v}(t)$ leads to vertically advected structures whose behaviour in the comoving frame is identical to that described below. We use $N_x=N_z=128$ grid points for moderate $Ra_{T,q}$, and set $Pr=1$.

3.3. Flow structures beyond elevator mode

Here, we choose the domain size $L_x=0.2{\rm \pi}$, unless otherwise mentioned, selected to accommodate the secondary instability of the elevator mode. With this domain size, the horizontal wavenumber of a domain-filling elevator mode corresponds to ${k_\perp =2{\rm \pi} /L_x=10}$, thus the critical Rayleigh number is $Ra_{T,q}=10^4$ according to (3.4). A larger domain will instead display a stable finite-amplitude elevator mode up to ${Ra_{T,q}=10^8}$ at $Pr=1$, as demonstrated below through secondary instability analysis (figure 12) and DNS (figure 13).

Figure 1 shows the resulting bifurcation diagram using thick (thin) lines for stable (unstable) states obtained from numerical continuation. The markers in figure 1 correspond to the final stable state as obtained from DNS. Figure 1(a) shows that the elevator mode (EM, black) bifurcates from the primary instability at $Ra_{T,q}^{(p)}=10^4$, consistent with (3.4), and that the Nusselt number of this mode displays a linear relation with $Ra_{T,q}$ according to (3.7). The secondary instability of the elevator mode at $Ra_{T,q}^{(s)}=19576.3$ leads to a branch of steady tilted elevator modes (TEMs, blue) accompanied by large-scale shear. Figure 2(a) shows the evolution of this shear from DNS at $Ra_{T,q}=3\times 10^4$, starting from an unstable elevator mode at this Rayleigh number. The figure shows that the large-scale shear $\langle u\rangle _h(z,t)$ becomes non-zero at $t\approx 2$ and then saturates in a horizontal flow with an approximately sinusoidal profile in the vertical. Figure 2(b) shows that the associated temperature deviation $T(x,z,t)$ at $t=10$ is tilted in the direction corresponding to the generated large-scale shear. This secondary bifurcation resembles the behaviour observed in RBC between fixed temperature boundaries, whereby a secondary bifurcation of steady convection rolls leads to tilted rolls accompanied by large-scale shear (Howard & Krishnamurti 1986; Rucklidge & Matthews 1996), as also observed in both experiments (Krishnamurti & Howard 1981) and DNS (Matthews et al. 1996; Goluskin et al. 2014; Von Hardenberg et al. 2015; Wang et al. 2020). This TEM branch terminates in another unstable steady state (red) that bifurcates from the conduction state without large-scale shear generation and resembles the two-layer (S2) solutions identified in salt-finger convection (Liu, Julien & Knobloch 2022).

Figure 1. (a) Bifurcation diagram with elevator mode (EM, black line), tilted elevator mode (TEM, blue line), direction-reversing state (DRS, magenta square) and modulated travelling waves (MTW, green cross). The bifurcation points include the primary bifurcation $Ra_{T,q}^{(p)}$, the secondary bifurcation $Ra_{T,q}^{(s)}$, the Hopf bifurcation $Ra_{T,q}^{(h)}$, and a global bifurcation at $Ra_{T,q}^{(g)}$. (b) Hysteresis diagram near the Hopf bifurcation point $Ra_{T,q}^{(h)}$ with TEM initial condition and increasing $Ra_{T,q}$ (blue star) or DRS initial conditions and decreasing $Ra_{T,q}$ (magenta square).

Figure 2. (a) Large-scale shear $\langle u\rangle _h(z,t)$ and (b) temperature deviation $T(x,z,t)$ at $t=10$ for the steady TEM at $Ra_{T,q}=3\times 10^4$, $Pr=1$ and $L_x=0.2{\rm \pi}$.

The steady TEM loses stability at a Hopf bifurcation at $Ra_{T,q}^{(h)}=32085.1$ leading to oscillations about the TEM state with frequency $\omega _h=43.1$. This Hopf bifurcation is subcritical, however, implying that the tilted oscillations are unstable. Computations indicate that the system instead evolves into a symmetric direction-reversing state (DRS) with associated hysteresis near $Ra_{T,q}^{(h)}$ as shown in figure 1(b). Here, we use $\langle nu(t)\rangle _t$ to distinguish the TEM state from the DRS, reached from TEM initial conditions upon increasing $Ra_{T,q}$ (blue star) or from DRS initial conditions upon decreasing $Ra_{T,q}$ (magenta square). When $Ra_{T,q}< Ra_{T,q}^{(h)}$, the DNS with TEM initial conditions show excellent agreement with numerical continuation results (thick blue line), while DNS with TEM initial conditions at $Ra_{T,q}>Ra_{T,q}^{(h)}$ evolve into DRS. Figure 3 shows $\langle u \rangle _h(z,t)$ and $nu(t)$ at $Ra_{T,q}=32\,100$, a value slightly larger than $Ra_{T,q}^{(h)}$, with a TEM initial condition from a lower $Ra_{T,q}$. The large-scale flow oscillates with frequency $\omega \approx 43.0$ in $t\in [0,10]$ (not shown in figure 3) that is close to the Hopf frequency $\omega _h=43.1$ but does not reverse. At $t\approx 74$, the flow abruptly transitions to a DRS as shown in figure 3(a). Figure 3(b) shows the corresponding instantaneous Nusselt number $nu(t)$. Similar DRS are also observed in RBC (Sugiyama et al. 2010; Chandra & Verma 2013; Winchester, Dallas & Howell 2021), magnetoconvection (Matthews et al. 1993; Proctor et al. 1994) as well as in salt-finger convection (Liu et al. 2022).

Figure 3. (a) Large-scale shear $\langle u\rangle _{h}(z,t)$ and (b) instantaneous Nusselt number $nu(t)$ at $Ra_{T,q}=3.21\times 10^4$, $Pr=1$ and $L_x=0.2{\rm \pi}$ with a TEM initial condition (supplementary movie 1, available at https://doi.org/10.1017/jfm.2023.1057, shows the corresponding temperature deviation $T(x,z,t)$).

At higher Rayleigh numbers, this DRS spends more time displaying flow structures close to an elevator mode. Figure 4 shows three snapshots of the temperature deviation $T(x,z,t)$ at $Ra_{T,q}=4.6\times 10^4$. At $t=2.21$ and $t=2.37$, the temperature deviation tilts in opposite directions. At $t=2.29$, $T(x,z,t)$ displays flow structures close to an elevator mode, followed at $t=2.37$ by a restored and approximately reflected tilted state. At yet higher Rayleigh numbers, the DRS collides with the unstable steady elevator mode leading to a global bifurcation at $Ra_{T,q}^{(g)}\approx 46892.03$, as indicated in figure 1(a). At this global bifurcation, the DRS transitions to modulated travelling waves (MTW, green) that do not reverse direction. Figure 5(a) shows the corresponding large-scale shear $\langle u\rangle _h(z,t)$ at $Ra_{T,q}=6\times 10^4$. Figure 5(b) displays the corresponding temperature deviation $T(x,z,t)$ at $z=0.1$, including the MTW that sets in at $t\approx 0.8$. This global bifurcation is illustrated in the phase diagram shown in figure 6 near $Ra_{T,q}^{(g)}$, revealing an abrupt change in topology before and after this global bifurcation. Note that both states pass through $\langle T\rangle _h (z_p, t)=\langle u\rangle _h(z_p, t)=0$ corresponding to the elevator mode. A similar global bifurcation must take place on the subcritical DRS branch in order to generate the DRS from the oscillating TEM state, but is inaccessible to DNS. Such gluing bifurcations are also seen in RBC, where oscillatory tilted convection rolls originating from a Hopf bifurcation of steady tilted convection rolls may collide with steady convection rolls and glue together in a global bifurcation; see e.g. Rucklidge & Matthews (1996).

Figure 4. Temperature deviation $T(x,z,t)$ associated with the DRS at three different times when $Ra_{T,q}=4.6\times 10^4$, $Pr=1$ and $L_x=0.2{\rm \pi}$: (a) $t=2.21$, (b) $t=2.29$, and (c) $t=2.37$.

Figure 5. (a) Large-scale shear $\langle u\rangle _h(z,t)$ and (b) temperature deviation $T(x,z,t)$ at $z=0.1$ for MTWs at $Ra_{T,q}=6\times 10^4$, $Pr=1$ and $L_x=0.2{\rm \pi}$.

Figure 6. Phase diagram showing $\langle T\rangle _h(z_p,t)$ as a function of $\langle u\rangle _h(z_p,t)$ at (a) $Ra_{T,q}=46\,892.0$ and(b) $Ra_{T,q}=46\,892.1$, with $z_p:={\text {arg\,max}}_{z} \langle T\rangle _h(z,t=10)$. The global bifurcation takes place in between.

3.4. Single-mode equations

The previous results show that the flow in the horizontal direction is dominated by a domain-filling mode. Moreover, in fixed-flux RBC, any long box will eventually contain a single pair of rolls (Chapman & Proctor 1980), and domain-filling modes also organize the flow in the turbulent regime at high Rayleigh number (Vieweg et al. 2021, 2022; Käufer et al. 2023). This motivates us to derive single-mode equations that have been used successfully in a wide range of convection problems (Herring 1963; Toomre, Gough & Spiegel 1977; Gough & Toomre 1982; Paparella & Spiegel 1999), especially for well-organized columnar structures in the presence of strong restraining body forces, including rapid rotation and strong magnetic field (Julien & Knobloch 2007), or large-scale damping in salt-finger convection (Liu et al. 2022), or convection in a porous medium (Liu & Knobloch 2022).

Single-mode equations are obtained from a severely truncated Fourier expansion in the horizontal, which reduces the governing equations from three spatial dimensions to equations for the vertical solutions profile associated with a prescribed horizontal planform. Here, we derive the single-mode equations by decomposing variables into a mean mode in the horizontal and horizontal harmonics:

(3.8a)$$\begin{gather} T(x,y,z,t)=\bar{T}_0(z,t)+\hat{T}(z,t)\exp({\text{i}(k_x x+k_y y)})+\text{c.c.}, \end{gather}$$

(3.8b)$$\begin{gather}\boldsymbol{u}(x,y,z,t)=\bar{U}_0(z,t)\,\boldsymbol{e}_x+\hat{\boldsymbol{u}}(z,t) \exp({\text{i}(k_x x+k_y y)})+\text{c.c.}, \end{gather}$$

(3.8c)$$\begin{gather}p(x,y,z,t)=\bar{P}_0(z,t)+\hat{p}(z,t)\exp({\text{i}(k_x x+k_y y)})+\text{c.c.} \end{gather}$$

We truncate the resulting equations at these harmonics to obtain the single-mode equations:

(3.9a)$$\begin{gather} \partial_t \hat{u}+\bar{U}_0\text{i}k_x \hat{u}+\hat{w}\,\partial_z \bar{U}_0={-} \text{i}k_x \hat{p}+Pr\,\hat{\nabla}^2 \hat{u}, \end{gather}$$

(3.9b)$$\begin{gather}\partial_t \hat{v}+\bar{U}_0\text{i}k_x\hat{v}={-}\text{i}k_y \hat{p}+Pr\,\hat{\nabla}^2 \hat{v}, \end{gather}$$

(3.9c)$$\begin{gather}\partial_t \hat{w}+\bar{U}_0\text{i}k_x\hat{w}={-}\partial_z \hat{p}+Pr\,\hat{\nabla}^2\hat{w}+Pr\,Ra_{T,q}\,\hat{T}, \end{gather}$$

(3.9d)$$\begin{gather}\text{i}k_x \hat{u}+\text{i}k_y\hat{v}+\partial_z \hat{w}=0, \end{gather}$$

(3.9e)$$\begin{gather}\partial_t\hat{T}+\bar{U}_0\text{i}k_x\hat{T}+\hat{w}\,\partial_z \bar{T}_0- \hat{w}+\hat{w}\int_0^1\left(\hat{w}^* \hat{T}+\hat{w}\hat{T}^*\right){\rm d}z= \hat{\nabla}^2\hat{T}, \end{gather}$$

(3.9f)$$\begin{gather}\partial_t \bar{U}_0+\partial_z\left(\hat{w}^*\hat{u}+\hat{w}\hat{u}^*\right)= Pr\,\partial_z^2\bar{U}_0, \end{gather}$$

(3.9g)$$\begin{gather}\partial_t \bar{T}_0+\partial_z\left(\hat{w}^*\hat{T}+\hat{w}\hat{T}^*\right)= \partial_z^2 \bar{T}_0. \end{gather}$$

Here, the integral term in (3.9e) represents the fixed-flux constraint originating from $w\langle wT\rangle _{h,v}$ in (2.7c). The horizontal wavenumber is chosen as $k_x=10$ and $k_y=0$ corresponding to a domain-filling mode within a 2-D domain with $L_x=0.2{\rm \pi}$. Numerical continuation of the single-mode equations (3.9) is performed using pde2path (Uecker et al. 2014; Uecker 2021a) with $N_z=128$, while DNS of (3.9) are conducted using Dedalus (Burns et al. 2020) with $N_z=128$.

Figure 7(a) shows the bifurcation diagram, while figure 7(b) shows the hysteresis diagram obtained from the single-mode equations in (3.9). Here, we can see that the single-mode equations reproduce the bifurcation and hysteresis diagrams obtained from the full equations in two dimensions, shown in figure 1. The hysteresis behaviour in figure 7(b) is present in a similar Rayleigh number range, $\Delta Ra_{T,q}\approx 200$, as in figure 1(b). Nevertheless, the bifurcation points are shifted slightly in the single-mode equations compared with the full equations, as shown in table 1. The success of the single-mode equations in predicting the Hopf frequency is perhaps in the same spirit as the real zero imaginary frequency (RZIF) ansatz that has shown success in predicting the oscillation frequency in nonlinear thermosolutal convection and shear flow (Turton, Tuckerman & Barkley 2015; Bengana et al. 2019; Bengana & Tuckerman 2021). Within the RZIF framework, the eigenvalues are computed based on dynamics linearized around a mean flow that can deviate from the laminar base flow, much as here the single-mode equations employ the large-scale modes $\bar {T}_0$ and $\bar {U}_0$ with superposed harmonics; see (3.9f) and (3.9g).

Figure 7. Same as figure 1 but obtained from the single-mode equations in (3.9).

Table 1. Comparison of the bifurcation points between the full 2-D equations in (2.7) with $L_x=0.2{\rm \pi}$, and the single-mode equations in (3.9) with $k_x=10$, including the primary bifurcation $Ra_{T,q}^{(p)}$, the secondary bifurcation $Ra_{T,q}^{(s)}$, the Hopf bifurcation $Ra_{T,q}^{(h)}$ with Hopf frequency $\omega _h$, and the global bifurcation $Ra_{T,q}^{(g)}$, all at $Pr=1$.

We further leverage the computational efficiency of single-mode equations to analyse the frequency scaling near the global bifurcation. We perform a bisection over $Ra_{T,q}$ using DNS of the single-mode equations in (3.9) to identify the global bifurcation point with more significant digits, $Ra_{T,q}^{(g)}=46\,761.0819762429$, than possible from the full equations. Figure 8 shows that the period $T_p$ of the direction reversals near $Ra_{T,q}^{(g)}$ diverges as $T_p=-0.0448\ln (Ra_{T,q}^{(g)}-Ra_{T,q})+0.6118$; cf. Knobloch & Proctor (1981) and Knobloch (1986).

Figure 8. (a) The reversal period $T_p$ as a function of $Ra_{T,q}$ near the global bifurcation $Ra_{T,q}^{(g)}$ obtained from the single-mode equations in (3.9) at $Pr=1$ and $k_x=10$. The black dashed line is $T_p=-0.0448\ln (Ra_{T,q}^{(g)}-Ra_{T,q})+0.6118$ and fits the DNS data with a relative residue of $0.4\,\%$, where $Ra_{T,q}^{(g)}=46\,761.0819762429$. (b) The leading eigenvalues of the (unstable) elevator mode at $Ra_{T,q}^{(g)}$ within the single-mode equations (3.9).

Since the DRS collides with an unstable elevator mode at the global bifurcation $Ra_{T,q}^{(g)}$, as indicated in figure 7(a), it is instructive to compute the eigenvalues of the elevator mode at this point within the single-mode equations (3.9). Figure 8(b) shows that the unstable eigenvalue $\lambda _1=43.69$ is real, and that the least stable eigenvalues are complex, with $\lambda _{2,3}=-\rho \pm \text {i}\omega =-49.27 \pm \text {i}117.11$. Thus this global bifurcation is associated with a saddle-focus equilibrium with $\delta \equiv \rho /\lambda _1=1.13>1$, i.e. the tame version of the Shilnikov bifurcation (Shilnikov 1965; Shilnikov & Shilnikov 2007). Here, we report these eigenvalues from the single-mode equations (3.9) to facilitate direct comparison of the logarithmic scaling law in figure 8(a) with theory. For the 2-D full equations, the corresponding eigenvalues of the elevator mode at the corresponding $Ra_{T,q}^{(g)}$ are $\lambda _1=41.11$, $\lambda _{2,3}=-\rho \pm \text {i}\omega =-56.67\pm \text {i}114.67$, leading to a Shilnikov bifurcation with $\delta =1.38$.

The logarithmic scaling law and associated coefficient can be predicted by constructing a Poincaré map near the global bifurcation point $Ra_{T,q}^{(g)}$ and the saddle-focus equilibrium (here the steady elevator mode) by composing a local map near this saddle focus and a global map (Shilnikov & Shilnikov 2007), as done by Glendinning & Sparrow (1984). For the local map, we consider the flow linearized around this elevator mode with $\mu :=Ra_{T,q}^{(g)}-Ra_{T,q}\ll 1$:

(3.10a)$$\begin{gather} \dot{\zeta}=\lambda_1 \zeta +\text{h.o.t.}, \end{gather}$$

(3.10b)$$\begin{gather}\dot{\theta}=\omega+\text{h.o.t.}, \end{gather}$$

(3.10c)$$\begin{gather}\dot{r}={-}\rho r +\text{h.o.t.}, \end{gather}$$

where $\zeta$ is the coordinate corresponding to the unstable eigenvalue, $(r,\theta )$ are the polar coordinates associated with the least stable eigenvalues, and h.o.t. refers to higher-order terms. We consider the Poincaré section $\varSigma ^{in}:=\{\theta =0\}$ and $\varSigma ^{out}:=\{\zeta =H\}$, and construct the local map $\varPi _{loc}:\varSigma ^{in}\rightarrow \varSigma ^{out}$ according to the linearized dynamics in (3.10):

(3.11a)$$\begin{gather} H=\zeta(T_f)=\zeta_0\,{\rm e}^{\lambda_1 T_f}, \end{gather}$$

(3.11b)$$\begin{gather}r(T_f)=r_0\,{\rm e}^{-\rho T_f}, \end{gather}$$

(3.11c)$$\begin{gather}\theta(T_f)=\omega T_f. \end{gather}$$

From (3.11a), the time of flight $\varSigma ^{in}:=\{\theta =0\}\to \varSigma ^{out}:=\{\zeta =H\}$ is

(3.12)\begin{equation} T_f={-}\frac{1}{\lambda_1}\ln \left( \frac{\zeta_0}{H}\right). \end{equation}

Substituting (3.12) into (3.11), we obtain the local map:

(3.13) \begin{equation} \varPi_{loc}:(r,\theta,\zeta)\rightarrow \biggl( r_0\left(\frac{\zeta_0}{H} \right)^{\delta},\frac{\omega}{\lambda_1}\ln\left(\frac{H}{\zeta_0}\right),H \biggr). \end{equation}

The global map $\varPi _{global}:\varSigma ^{out}\rightarrow \varSigma ^{in}$ is obtained from a Taylor series around the homoclinic orbit assumed to be present at $\mu =0$:

(3.14) \begin{align} \varPi_{global}: (r,\theta,h)\rightarrow (\bar{r}+a\mu+br\cos \theta +cr\sin\theta,0,d\mu+er\cos\theta+fr\sin\theta)+\text{h.o.t.}, \end{align}

where $a$, $b$, $c$, $d$, $e$ and $f$ are constants. By composing the local and global maps ($\varPi :\varSigma ^{in}\rightarrow \varSigma ^{in}=\varPi _{global}\circ \varPi _{loc}$), we obtain the Poincaré map

(3.15)\begin{equation} \varPi: \begin{bmatrix} r\\ \zeta \end{bmatrix}\rightarrow \begin{bmatrix} \bar{r}\\ 0 \end{bmatrix}+\begin{bmatrix} a\\ d \end{bmatrix}\mu+\begin{bmatrix} c_1 r \zeta^{\delta}\cos(k_1 \ln \zeta +\phi_1)\\ c_2 r \zeta^{\delta} \cos(k_2 \ln \zeta +\phi_2) \end{bmatrix}+\text{h.o.t.}, \end{equation}

where $c_i$, $k_i$ and $\phi _i$ ($i=1,2$) are constants. We may now search for a fixed point of the Poincaré map $\varPi$ that corresponds to the periodic orbit near $Ra_{T,q}^{(g)}$ in the original system. This point is approximated by the fixed point of the one-dimensional map

(3.16)\begin{equation} \zeta-d \mu=\varPhi(\zeta)\equiv c_2 r \zeta^{\delta}\cos(k_2\ln \zeta +\phi_2)+\text{h.o.t.} \end{equation}

When $\delta >1$, there is a unique fixed point of (3.16), which scales as $\zeta \sim d\mu$ near the global bifurcation $\mu \rightarrow 0$. Thus based on (3.12) and the assumption that the global return is much faster than the local passage past the fixed point, the period of the reversing orbit just before the global bifurcation scales as

(3.17)\begin{equation} T_p={-}\frac{2}{\lambda_1} \ln \mu +\text{const.}={-}\frac{2}{\lambda_1}\ln \left(Ra_{T,q}^{(g)}-Ra_{T,q}\right)+\text{const}. \end{equation}

Here, the factor 2 arises because the orbit makes two passes near the fixed point in each reversal period. Using $\lambda _1$ from figure 8(b), this calculation predicts that $2/\lambda _1=0.04578$, a coefficient that is almost exactly that obtained from the fit to the simulation data in figure 8(a).

4. Dynamics at high Rayleigh numbers

In this section, we study the dynamics at higher Rayleigh numbers, where chaotic behaviour appears. We first analyse the secondary instability of the elevator mode, which continues to play an important role in the high Rayleigh number regime. We focus on the 2-D elevator mode with a horizontal wavenumber $k_x=k_e$ in the $x$ direction,

(4.1a,b)\begin{equation} \bar{W}_e(x)=\hat{w}_e \exp(\text{i}k_e x)+{\rm c.c}.,\quad \bar{T}_e(x)=\hat{T}_e\exp(\text{i}k_e x)+{\rm c.c}., \end{equation}

and solution amplitude given by (3.6a,b). The decomposition

(4.2a,b)\begin{equation} \boldsymbol{u}=\bar{W}_e\boldsymbol{e}_z+\boldsymbol{u}',\quad T=\bar{T}_e+T' \end{equation}

leads to the linearized equations

(4.3a)$$\begin{gather} \frac{\partial \boldsymbol{u}'}{\partial t}+u'\,\frac{{\rm d}\bar{W}_e}{{\rm d}\kern 0.06em x}\, \boldsymbol{e}_z+\bar{W}_e\,\partial_z \boldsymbol{u}'={-}\boldsymbol{\nabla}p'+Pr\,Ra_{T,q}\,T'\boldsymbol{e}_z+Pr\, \nabla^2 \boldsymbol{u}', \end{gather}$$

(4.3b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}'=0, \end{gather}$$

(4.3c)$$\begin{gather}\frac{\partial T'}{\partial t}+u'\,\frac{{\rm d}\bar{T}_e}{{\rm d}\kern 0.06em x}+ \bar{W}_e\,\partial_z T'-w'+w'\langle \bar{W}_e\bar{T}_e\rangle_{h,v}=\nabla^2T'. \end{gather}$$

The cubic flux-feedback nonlinearity $w\langle wT\rangle _{h,v}$ in (2.7c) generates three linearized terms:

(4.4)\begin{equation} w'\langle \bar{W}_e\bar{T}_e\rangle_{h,v}+\bar{W}_e\langle w'\bar{T}_e\rangle_{h,v}+\bar{W}_e\langle \bar{W}_eT'\rangle_{h,v}, \end{equation}

with the latter two terms in (4.4) vanishing for $k_z\neq 0$. As a result, only the term $w'\langle \bar {W}_e\bar {T}_e\rangle _{h,v}$ originating from flux feedback appears in the linearized equation in (4.3c). The normal mode assumption in general 3-D form

(4.5a)$$\begin{gather} \boldsymbol{u}'=\tilde{\boldsymbol{u}}(x) \exp[\text{i}(k_y y+k_z z)+\lambda t]+\text{c.c.}, \end{gather}$$

(4.5b)$$\begin{gather}T'=\tilde{T}(x)\exp[\text{i}(k_y y+k_z z)+\lambda t]+\text{c.c.} \end{gather}$$

contains the coefficients $\tilde {\boldsymbol {u}}(x)$ and $\tilde {T}(x)$ that depend on $x$ because the base flow (elevator mode) also depends on $x$. In terms of the horizontal vorticity $\tilde {\omega }_x:=\text {i}k_y \tilde {w}-\text {i}k_z\tilde {v}$, we have the linear eigenvalue problem

(4.6)\begin{equation} \lambda \begin{bmatrix} \tilde{u}\\ \tilde{\omega}_x\\ \tilde{T} \end{bmatrix} = \begin{bmatrix} \mathcal{A}_{11} & \mathcal{A}_{12} & \mathcal{A}_{13}\\ \mathcal{A}_{21} & \mathcal{A}_{22} & \mathcal{A}_{23} \\ \mathcal{A}_{31} & \mathcal{A}_{32} & \mathcal{A}_{33} \end{bmatrix} \begin{bmatrix} \tilde{u}\\ \tilde{\omega}_x\\ \tilde{T} \end{bmatrix}=:\mathcal{A}\begin{bmatrix} \tilde{u}\\ \tilde{\omega}_x\\ \tilde{T} \end{bmatrix}, \end{equation}

where

(4.7a)$$\begin{gather} \mathcal{A}_{11}=\tilde{\nabla}^{{-}2}\left(-\text{i}k_z \bar{W}_e\,\tilde{\nabla}^2+ \text{i}k_z\,\frac{{\rm d}^2\bar{W}_e}{{{\rm d}\kern 0.06em x}^2}+Pr\,\tilde{\nabla}^4\right), \end{gather}$$

(4.7b)$$\begin{gather}\mathcal{A}_{12}=0, \end{gather}$$

(4.7c)$$\begin{gather}\mathcal{A}_{13}=\tilde{\nabla}^{{-}2}\,Pr\,Ra_{T,q}\,(-\text{i}k_z\,\partial_x), \end{gather}$$

(4.7d)$$\begin{gather}\mathcal{A}_{21}={-}\text{i}k_y\,\frac{{\rm d}\bar{W}_e}{{\rm d}\kern 0.06em x}, \end{gather}$$

(4.7e)$$\begin{gather}\mathcal{A}_{22}={-}\text{i}k_z \bar{W}_e+Pr\,\tilde{\nabla}^2, \end{gather}$$

(4.7f)$$\begin{gather}\mathcal{A}_{23}=Pr\,Ra_{T,q}\,\text{i}k_y, \end{gather}$$

(4.7g)$$\begin{gather}\mathcal{A}_{31}={-}\frac{{\rm d}\bar{T}_e}{{\rm d}\kern 0.06em x}+[1-\langle \bar{W}_e\bar{T}_e \rangle_{h,v}]\,\frac{\text{i}k_z\,\partial_x }{k_y^2+k_z^2}, \end{gather}$$

(4.7h)$$\begin{gather}\mathcal{A}_{32}=[1-\langle \bar{W}_e\bar{T}_e\rangle_{h,v}]\, \frac{-\text{i}k_y}{k_y^2+k_z^2}, \end{gather}$$

(4.7i)$$\begin{gather}\mathcal{A}_{33}={-}\text{i}k_z \bar{W}_e+\tilde{\nabla}^2, \end{gather}$$

with $\tilde {\nabla }^2:=\partial _x^2-k_y^2-k_z^2$, and $\tilde {\nabla }^4:=\partial _x^4-2(k_y^2+k_z^2)\,\partial _x^2+(k_y^2+k_z^2)^2$.

We compare the above formulation with that without the flux feedback, corresponding to setting the integral flux-feedback terms $\langle \bar {W}_e\bar {T}_e\rangle _{h,v}$ in (4.7g)–(4.7h) to zero, leading to a modified eigenvalue problem with $\mathcal {A}$ in (4.6) replaced by

(4.8)\begin{equation} \underline{\mathcal{A}}:=\begin{bmatrix} \mathcal{A}_{11} & \mathcal{A}_{12} & \mathcal{A}_{13}\\ \mathcal{A}_{21} & \mathcal{A}_{22} & \mathcal{A}_{23}\\ \underline{\mathcal{A}}_{31} & \underline{\mathcal{A}}_{32} & \mathcal{A}_{33} \end{bmatrix}, \end{equation}

where

(4.9a,b)\begin{equation} \underline{\mathcal{A}}_{31}:={-}\frac{{\rm d}\bar{T}_e}{{\rm d}\kern 0.06em x}+\frac{\text{i}k_z\, \partial_x }{k_y^2+k_z^2},\quad \underline{\mathcal{A}}_{32}:=\frac{-\text{i}k_y}{k_y^2+k_z^2}. \end{equation}

The horizontal direction is discretized using a Fourier collocation method with the horizontal derivative computed using a Fourier differentiation matrix (Weideman & Reddy 2000). The numerical implementation is validated against Floquet-based linear stability analysis (Holyer 1984; Garaud, Gallet & Bischoff 2015; Radko 2016; Garaud, Kumar & Sridhar 2019). We choose the horizontal domain $L_x$ to contain one or more wavelengths of the elevator wavelength $2{\rm \pi} /k_e$. For all the results reported here, we take $k_y=0$ corresponding to a 2-D configuration.

Figure 9(a) shows the growth rate $\max [{\rm Re}(\lambda )]$, comparing the fixed-flux case computed from $\mathcal {A}$ in (4.6) with the case without flux feedback computed from $\underline {\mathcal {A}}$ in (4.8). The flux feedback leads to $\lambda =0$ at $k_z=0$, and the instability is limited to a small range of $k_z$, a feature observed widely in systems with a conservation law (Matthews & Cox 2000). For the case without flux feedback, its growth rate is larger and decays to zero only at much higher wavenumbers $k_z$ (not shown in figure 9a). We further identify the wavenumbers $k_{z,0}$ and $k_{z,max}$ indicated in figure 9(a) corresponding, respectively, to zero growth rate and maximum growth rate:

(4.10a)$$\begin{gather} \lambda(k_{z,0})=0\quad \text{with } k_{z,0}\neq 0, \end{gather}$$

(4.10b)$$\begin{gather}k_{z,max}:=\underset{k_z}{\text{arg\,max}} \,{\rm Re}(\lambda). \end{gather}$$

Figure 9(b) displays $k_{z,0}$ and $k_{z,max}$, both of which increase as $Ra_{T,q}$ increases, with $k_{z,max}< k_{z,0}$. Here, we also plot $k_{z}=2{\rm \pi}$, which is the smallest non-trivial vertical wavenumber that fits in the $L_z=1$ domain. The secondary instability wavelength is not required to lie within $L_z=1$, thus a comparison of $k_{z,0}$ with $k_z=2{\rm \pi}$ determines whether the secondary instability can occur within the domain, or whether it is suppressed by its finite size. Figure 9(b) shows that $k_{z,0}>2{\rm \pi}$ when $L_x=0.2{\rm \pi}$ and $Ra_{T,q}\geq 3\times 10^4$, a result consistent with the presence of a secondary instability in an $L_z=1$ domain identified in § 3.

Figure 9. (a) The growth rate $\max [{\rm Re}(\lambda )]$ at $Ra_{T,q}=10^8$, $Pr=1$, $L_x=0.2{\rm \pi}$ and $k_e=10$, with flux feedback computed from $\mathcal {A}$ in (4.6), and without flux feedback as computed from $\underline {\mathcal {A}}$ in (4.8). (b) The wavenumbers $k_{z,0}$ and $k_{z,max}$ as functions of $Ra_{T,q}$ at $Pr=1$, $L_x=0.2{\rm \pi}$ and $k_e=10$. The red dashed line corresponds to $k_z=2{\rm \pi}$.

The secondary instability of the elevator mode and the MTWs generated through the sequence of bifurcations examined in § 3 continue to play an important role at larger Rayleigh numbers. Figure 10 displays DNS results at $Ra_{T,q}=10^8$ obtained with $N_x=N_z=256$ grid points ($N_x=N_z=512$ grid points generate the same behaviour). The large-scale shear $\langle u\rangle _{h}(z,t)$ in figure 10(a) is now dominated by MTWs, but displays chaotic behaviour. In addition, it slowly migrates in the vertical direction, behaviour that is permitted by the periodic boundary conditions in the vertical.

Figure 10. (a) Large-scale shear $\langle u\rangle _h(z,t)$, (b) total temperature $1+ \bar {\mathcal {T}}_{z,q} z+\langle T\rangle _{h,t}(z)$ (solid black line) compared with the linear profile $1+ \bar {\mathcal {T}}_{z,q} z$ (dashed blue line), and (c) root mean square vertical velocity $\sqrt {\langle ww\rangle _{h,t}}$, the latter two averaged over $t\in [0.275,0.465]$, at $Ra_{T,q}=10^8$, $Pr=1$ and $L_x=0.2{\rm \pi}$.

The mean total temperature averaged over $t\in [0.275,0.465]$ in figure 10(b) exhibits a deviation from a linear profile similar to canonical RBC. Moreover, in the region where the mean temperature deviation $\langle T\rangle _{h,t}(z)$ is close to zero, the corresponding large-scale shear $\langle u\rangle _{h}(z,t)$ also vanishes (white regions in figure 10a) and the root mean square (r.m.s.) vertical velocity $\sqrt {\langle ww\rangle _{h,t}}$ displays local minima with a zero vertical derivative, as shown in figure 10(c). The local maxima of the r.m.s. vertical velocity $\sqrt {\langle ww\rangle _{h,t}}$ correspond instead to the mixed mean temperature regions in figure 10(b). This behaviour resembles that of RBC in a bounded domain with constant temperature boundaries and no-slip instead of stress-free velocity boundary conditions at the top and bottom (van der Poel et al. 2014). The figure shows that the fixed-flux constraint suppresses bursting behaviour compared with homogeneous RBC driven by a constant temperature gradient (Borue & Orszag 1997; Lohse & Toschi 2003; Calzavarini et al. 2005, 2006). Figure 11 shows the Nusselt number scaling within $Ra_{T,q}\in [10^8,10^{10}]$ computed with $N_x=N_z=512$ grid points and $L_x=0.2{\rm \pi}$, $Pr=1$. Here, the fitted scaling law $Nu=0.189\,Ra_{T,q}^{0.217}$ is associated with an exponent lower than $Nu\sim Ra_{T,q}^{1/3}$ corresponding to ultimate regime scaling $Nu\sim Ra_T^{1/2}$. However, the flow structures associated with figure 11 are dominated by large-scale shear with a $\langle u\rangle _h(z,t)$ profile similar to that in figure 10(a), potentially reducing heat transport.

Figure 11. Plot of $Nu$ as a function of $Ra_{T,q}$ from DNS with $L_x=0.2{\rm \pi}$ and $Pr=1$.

The impact of the domain size and of the horizontal wavenumber of the elevator mode within the domain is analysed further in figure 12(a), which shows $k_{z,0}$ and $k_{z,max}$ as functions of $L_x$ when $k_e=2{\rm \pi} /L_x$, $Ra_{T,q}=10^8$ and $Pr=1$. As the horizontal domain and wavelength increase, both $k_{z,0}$ and $k_{z,max}$ decrease, suggesting a small vertical wavenumber of the secondary instability. With $L_z$ fixed at $L_z=1$, this requires that we choose a small enough horizontal domain such that $k_{z,0}\geq 2{\rm \pi}$. A similar requirement applies for mean flow generation in canonical RBC (Rucklidge & Matthews 1996; Fitzgerald & Farrell 2014; Wang et al. 2020) and is the reason for choosing $L_x=0.2{\rm \pi}$ for most of the results in this paper. Figure 12(b) then fixes $L_x$ at $L_x=0.4{\rm \pi}$, but varies $k_e$, allowing multiple elevator modes within the domain. Evidently, both $k_{z,0}$ and $k_{z,max}$ increase as $k_e$ increases, and both are fitted well by a linear scaling law. This trend can be confirmed also in DNS. The associated large-scale shear $\langle u\rangle _h(z,t)$ in figure 13(a) shows that $k_z=4{\rm \pi}$ for an initial elevator mode wavenumber $k_e=20$, while $\langle u\rangle _h(z,t)$ in figure 13(b) displays a $k_z=6{\rm \pi}$ instability associated with $k_e=30$. This observation corresponds directly to the $k_{z,max}$ value at $k_e=20$ and $30$ predicted by the secondary instability analysis in figure 12(b).

Figure 12. The wavenumbers $k_{z,0}$ and $k_{z,max}$ as functions of (a) $L_{x}$ when $k_e=2{\rm \pi} /L_x$, and (b) $k_e$ with $L_x=0.4{\rm \pi}$. The linear fit is $k_{z,0}=0.98k_e+0.10$ (black dash-dotted line) and $k_{z,max}=0.57k_e+0.20$ (blue dash-dotted line). The dashed red lines indicate $k_z=2{\rm \pi}$. Instability requires that $k_{z,0}>2{\rm \pi}$. Other parameters are $Ra_{T,q}=10^8$ and $Pr=1$.

Figure 13. Large-scale shear $\langle u\rangle _h(z,t)$ with initial elevator mode wavenumber (a) $k_e=20$ and (b) $k_e=30$, both at $Ra_{T,q}=10^8$, $Pr=1$ and $L_x=0.4{\rm \pi}$.

In a domain of horizontal size $L_x=0.4{\rm \pi}$, the final state is a stable domain-filling elevator mode with $k_e=5$ corresponding to the white region in figure 13 at $Ra_{T,q}=10^8$. A similar transition to a larger horizontal scale is also observed in fixed-flux RBC within both the weakly nonlinear regime (Chapman & Proctor 1980) and the turbulent regime (Vieweg et al. 2021, 2022; Käufer et al. 2023). The stable elevator mode that results suggests that vertical jets (i.e. elevator modes) are favoured in wide domains, while horizontal jets are found in narrow domains. This behaviour resembles that in rapidly rotating convection where jets parallel to the short side of an anisotropic domain are found (Julien, Knobloch & Plumley 2018), or in 2-D turbulence driven by stochastic forcing (Bouchet & Simonnet 2009). At higher Rayleigh numbers, the flow in wider domains does become turbulent despite the absence of a shear-generating instability. For example, when $L_x=0.4{\rm \pi}$, this instability is absent since $k_{z,0}<2{\rm \pi}$. However, when the Rayleigh number is increased to $Ra_{T,q}=10^{10}$, the flow nonetheless becomes turbulent or at least chaotic. Figure 14 shows two snapshots of a solution with $N_x=N_z=512$ grid points at these parameter values, displaying intermittent layering and vortex dipole generation. The mechanism responsible for destabilizing the stationary elevator state at these large Rayleigh numbers remains to be studied.

Figure 14. Snapshots of the temperature deviation $T(x,z,t)$ at (a) $t=0.03$ and (b) $t=0.032$. Parameters are $Ra_{T,q}=10^{10}$, $L_x=0.4{\rm \pi}$ and $Pr=1$.

5. Prandtl number effect

In this section, we investigate the effect of the Prandtl number using the full 2-D equations. A low Prandtl number is of interest in astrophysical applications, where the heat transport is dominated by photon diffusion (Garaud 2018, 2021). Prandtl number $Pr=7$ corresponds to thermal diffusivity in water appropriate to oceanographic applications. When the temperature field is exchanged for a concentration such as salinity, the corresponding thermal diffusivity is replaced by molecular diffusivity, leading to Prandtl numbers as large as $Pr=700$.

Figure 15 shows a bifurcation diagram similar to figure 1(a) but with ${Pr=0.1}$ and $Pr=7$ in figures 15(a,b), respectively. The onset of the primary instability $Ra_{T,q}^{(p)}$ and the amplitude of the resulting elevator mode branch are not affected by $Pr$, consistent with the analytical results in § 3.1. At low Prandtl numbers, the secondary bifurcation point $Ra_{T,q}^{(s)}$ and the Hopf bifurcation point $Ra_{T,q}^{(h)}$ are both shifted to much lower Rayleigh numbers, while increasing $Pr$ shifts these bifurcation points to a higher Rayleigh number. These bifurcation points and the associated Hopf frequency $\omega _h$ are listed in table 2. The Hopf frequency in table 2 obtained from numerical continuation is also very close to the oscillation frequency obtained from DNS at a parameter close to $Ra_{T,q}^{(h)}$ (not shown). Table 2 compares the bifurcation points computed from the full 2-D equations in (2.7) with the corresponding results from the single-mode equations in (3.9). The latter are in closer agreement with the full 2-D equations at ${Pr=0.1}$ than at $Pr=7$ because a low Prandtl number shifts these bifurcation points closer to the onset of primary instability at $Ra_{T,q}^{(p)}$. A related phenomenon is found in RBC, where at low $Pr$, a steady convection roll becomes immediately unstable to a large-scale (zonal) mode (Winchester, Howell & Dallas 2022). The orange markers in figure 15 show the Nusselt number of the time-dependent states obtained from DNS, where the low Prandtl number case follows closely that of the TEM branch.

Figure 15. Bifurcation diagram with an elevator mode (EM, black line), tilted elevator mode (TEM, blue) and time-dependent states (orange circles) at (a) ${Pr=0.1}$ and (b) $Pr=7$. The bifurcation points include the primary bifurcation $Ra_{T,q}^{(p)}$, the secondary bifurcation $Ra_{T,q}^{(s)}$, and a Hopf bifurcation $Ra_{T,q}^{(h)}$. The horizontal domain size is $L_x=0.2{\rm \pi}$.

Table 2. Comparison of bifurcation points between the full 2-D equations in (2.7) with $L_x=0.2{\rm \pi}$, and the single-mode equations in (3.9) with $k_x=10$, including the primary bifurcation $Ra_{T,q}^{(p)}$, the secondary bifurcation $Ra_{T,q}^{(s)}$, the Hopf bifurcation $Ra_{T,q}^{(h)}$, and the Hopf frequency $\omega _h$ at ${Pr=0.1}$ and $Pr=7$.

Figure 16(a) displays the large-scale shear $\langle u\rangle _h(z,t)$ found at $Ra_{T,q}=10\,800$ and ${Pr=0.1}$. The instantaneous Nusselt number grows exponentially from $Nu\approx 1$ corresponding to the conduction state, and saturates at $Nu\approx 1.08$ associated with the elevator mode based on (3.7). The Nusselt number then decreases rapidly to $Nu\approx 1$ as a result of shear generation. Four different snapshots corresponding to the four different markers in figure 16(b) are shown in figure 17. Here, the elevator mode grows from $t=160$ to $t=166$ and then tilts at $t=166.5$, followed by rapid decay back towards the conduction state at $t=167$. This behaviour resembles a relaxation oscillation between the conduction base state and the steady elevator mode associated with different time scales between the growth and decay of the elevator mode. A similar oscillation cycle between a sheared state and convection rolls is also observed in RBC (Matthews et al. 1996; Goluskin et al. 2014) and magnetoconvection (Matthews, Proctor & Weiss 1995). More generally, relaxation oscillations associated with disparate time scales are also observed widely in binary fluid convection (Batiste et al. 2006), double-diffusive convection (Beaume, Bergeon & Knobloch 2018) and convection in a rapidly rotating sphere (Busse 2002). Burst behaviour similar to the relaxation oscillation cycle observed here can be described by a 2-D ordinary differential equation model, where a burst is triggered by a symmetry-breaking instability of a growing symmetric state when it reaches a certain amplitude (Batiste et al. 2001; Bergeon & Knobloch 2002). The prevalence of large-scale shear at low $Pr$ is known to lead to suppression of convection, an effect observed widely in the low-$Pr$ regime of RBC; see e.g. Schumacher & Sreenivasan (2020, figure 9). The suppression of convection by large-scale shear resembles the suppression of turbulence by self-generated mean flows (Shats et al. 2007).

Figure 16. (a) Large-scale shear $\langle u\rangle _{h}(z,t)$ and (b) instantaneous Nusselt number $nu(t)$ displaying relaxation oscillations at $Ra_{T,q}=10\,800$, ${Pr=0.1}$ and $L_x=0.2{\rm \pi}$.

Figure 17. Evolution of the temperature deviation $T(x,z,t)$ during a burst corresponding to the four times indicated in figure 16(b): (a) $t=160$ (black star), (b) $t=166$ (blue circle), (c) $t=166.5$ (magenta square), and (d) $t=167$ (red cross). Parameters are $Ra_{T,q}=10\,800$, ${Pr=0.1}$ and $L_x=0.2{\rm \pi}$ (see supplementary movie 2).

Figure 18(a) shows $nu(t)$ at a higher Rayleigh number, $Ra_{T,q}=11\,900$. This time series no longer displays visible periodic behaviour, but its power spectrum density in figure 18(b) shows multiple peaks corresponding to frequencies $m\,f_1+n\,f_2$ with $f_1=0.1483$, $f_2=0.2044$ and $m,n\in \mathbb {Z}^+$. We conclude that the state in figure 18(a) is quasi-periodic, and likely generated by a torus bifurcation (equivalently, a Hopf bifurcation of a periodic orbit). Quasi-periodic orbits are also observed prior to the onset of non-periodic motion in RBC (Ahlers & Behringer 1978; Yahata 1982). At yet higher Rayleigh numbers, chaotic behaviour appears as shown in figure 19 for $Ra_{T,q}=13\,000$. The broad power spectrum density (PSD) is indicative of a chaotic state. The instantaneous Nusselt number does not reach the amplitude corresponding to the saturated steady elevator mode with $Nu=1.3$, as shown in figure 19(a). This is because the growing elevator modes are disrupted before they are able to saturate.

Figure 18. (a) The instantaneous Nusselt number $nu(t)$ and (b) its power spectrum density (PSD) within $t\in [200,450]$, displaying spectral peaks at frequencies $mf_1+nf_2$, where $f_1=0.1483$, $f_2=0.2044$ and $m,n\in \mathbb {Z}^+$, indicating quasi-periodic dynamics at $Ra_{T,q}=11\,900$, ${Pr=0.1}$ and $L_x=0.2{\rm \pi}$.

Figure 19. (a) The instantaneous Nusselt number $nu(t)$ and (b) its PSD at parameters $Ra_{T,q}=13\,000$, ${Pr=0.1}$ and $L_x=0.2{\rm \pi}$.

We next move on to the case $Pr=7$. Figure 20 displays $\langle u \rangle _h(z,t)$ and $nu(t)$ at $Ra_{T,q}=2.4\times 10^5$ and $Pr=7$. Here, MTWs are present, similar to those at $Pr=1$ in figure 5(a). The instantaneous Nusselt number shows irregular behaviour. At the higher Rayleigh number $Ra_{T,q}=3\times 10^5$, the large-scale shear can reverse its direction, as shown in figure 21(a). At this Rayleigh number, the instantaneous Nusselt number displays intermittent behaviour that can be temporarily much larger than that of the elevator mode ($Nu\approx 30$), as shown in figure 21(b). Similar intermittency in the heat transport is also observed in homogeneous RBC driven by a constant temperature gradient (Borue & Orszag 1997; Calzavarini et al. 2005, 2006) owing to the intermittent appearance of an elevator mode. Here, a larger Prandtl number suppresses the large-scale shear that disrupts elevator modes, leading to the observed intermittent bursting behaviour resembling that in homogeneous RBC driven by a constant temperature gradient.

Figure 20. (a) The mean horizontal flow $\langle u \rangle _h(z,t)$ and (b) instantaneous Nusselt number $nu(t)$ at $Ra_{T,q}=2.4\times 10^5$, $Pr=7$ and $L_x=0.2{\rm \pi}$.

Figure 21. (a) The mean horizontal flow $\langle u \rangle _h(z,t)$ and (b) instantaneous Nusselt number $nu(t)$ at $Ra_{T,q}=3\times 10^5$, $Pr=7$ and $L_x=0.2{\rm \pi}$.

Figure 22 shows four snapshots of the temperature deviation $T(x,z,t)$ close to the burst event in figure 21(b) at the four times indicated in the figure. Here, the local minima in the Nusselt number correspond to tilted states that can be associated with direction reversals, as shown in figure 22(a) at $t=3.62$ and figure 22(b) at $t=3.67$. The local maximum of the Nusselt number at $t=3.68$ (figure 22c) corresponds to a burst state with hot fluid moving upwards and cold fluid moving downwards without impediment. This burst significantly reduces the absolute value of the temperature gradient, which leads to a large Nusselt number based on their reciprocal relation in (2.8). This behaviour is similar to homogeneous RBC in the high Prandtl number regime, where the temperature field often leads to vertical jet formation associated with a strong influence on the vertical temperature gradient (Calzavarini et al. 2005).

Figure 22. Evolution of the temperature deviation $T(x,z,t)$ during a burst corresponding to the four times indicated in figure 21(b): (a) $t=3.62$ (black star), (b) $t=3.67$ (blue circle), (c) $t=3.68$ (magenta square), and (d) $t=3.69$ (red cross). Parameters are $Ra_{T,q}=3\times 10^5$, $Pr=7$ and $L_x=0.2{\rm \pi}$.

We next move on to $Pr=700$, corresponding to the molecular diffusivity of salt in water. At this parameter, the elevator mode remains ‘stable’ within DNS with $L_x=0.2{\rm \pi}$ up to $Ra_{T,q}=10^8$. This can be understood from the secondary instability of elevator mode as shown in figure 23, where the onset vertical wavenumber $k_{z,0}$ at $Pr=700$ is lower than $k_z=2{\rm \pi}$, which is the minimum wavenumber for the presence of a secondary instability in a domain with $L_z=1$. In order to incorporate the secondary instability of the elevator mode, we change the horizontal domain size to $L_x=0.1{\rm \pi}$, associated with domain-filling wavenumber $k_x=2{\rm \pi} /L_x=20$, and increase the Rayleigh number to $Ra_{T,q}=8\times 10^8$.

Figure 23. The wavenumbers $k_{z,0}$ and $k_{z,max}$ as functions of $Pr$ at $Ra_{T,q}=10^8$, $L_x=0.2{\rm \pi}$ and $k_e=10$. The red dashed line corresponds to $k_z=2{\rm \pi}$.

Figure 24(a) displays the mean horizontal flow $\langle u \rangle _h(z,t)$ at $Ra_{T,q}=8\times 10^8$, $L_x=0.1{\rm \pi}$ and $Pr=700$. This large-scale shear begins to drift vertically at $t\approx 0.3$, although its magnitude is in fact much reduced from that at $Pr=1$ (figure 10a). Moreover, as shown in figure 24(b), the mean temperature profile averaged over $t\in [0.02,0.26]$ exhibits apparent layer formation and several narrow regions displaying a stably stratified temperature profile. Such stably stratified regions are also present in extended domain DNS of RBC with $Pr=7$, but disappear when the Prandtl number decreases (Schumacher & Sreenivasan 2020, figure 9g). Related locally stably stratified regions are also observed in potential vorticity staircases (Read et al. 2006, figure 5) and rotating double-diffusive convection (Moll & Garaud 2016, figures 6–8). Note that due to the vertical drift, the mean temperature profile is close to a linear profile if averaged over longer times, as shown in figure 24(c).

Figure 24. (a) The mean horizontal flow $\langle u\rangle _h(z,t)$, and (b,c) the mean temperature $1+ \bar {\mathcal {T}}_{z,q} z+\langle T\rangle _{h,t}(z)$ (solid black line) compared with the linear profile $1+ \bar {\mathcal {T}}_{z,q} z$ (dashed blue line) averaged, respectively, over (b) $t\in [0.02,0.26]$ and (c) $t\in [0.5,1]$. The parameters are $Ra_{T,q}=8\times 10^8$, $Pr=700$ and $L_x=0.1{\rm \pi}$.

Figure 25 shows four snapshots of the temperature deviation $T(x,z,t)$ at $Pr=700$ and $Ra_{T,q}=8\times 10^8$. The temperature deviation in figures 25(a) and 25(d) displays tilted states associated with different tilting directions, with Nusselt numbers that are smaller than that of the elevator mode. Figure 25(b) shows a burst associated with the temporary creation of a positive or stable instantaneous mean temperature gradient $\bar {T}_{z,q}(t):=\langle wT\rangle _{h,v}-1$. At a slightly later instant, shown in figure 25(c), the temperature deviation is slightly tilted but still close to an elevator mode with an instantaneous Nusselt number close to that of an elevator mode ($Nu=Ra_{T,q}/k_x^4=5000$). In homogeneous RBC driven by a constant temperature gradient, elevator modes appear more frequently at high Prandtl numbers (Calzavarini et al. 2005), which is also the case here in the fixed-flux set-up, potentially due to a weaker large-scale shear that plays such an important role in the secondary instability of elevator modes.

Figure 25. Temperature deviation $T(x,z,t)$ at (a) $t=0.456$ with $nu(t)=1165.5$, (b) $t=0.457$ with $\bar {T}_{z,q}(t)=1.75\times 10^{-5}$, (c) $t=0.459$ with $nu(t)=4784.7$, and (d) $t=0.461$ with $nu(t)=961.5$, at $Ra_{T,q}=8\times 10^8$, $Pr=700$ and $L_x=0.1{\rm \pi}$, showing the generation of momentary stable stratification in case (b) (see supplementary movie 3).

6. Conclusions and future work

This work formulated a fixed-flux problem for homogeneous Rayleigh–Bénard convection (RBC) and analysed its underlying dynamics in detail using numerical continuation, secondary instability analysis and direct numerical simulations. The fixed-flux formulation leads to steady elevator mode solutions with a well-defined amplitude, something that is not the case in homogeneous RBC driven by a constant temperature gradient. The secondary instability of this elevator mode leads to tilted elevator modes (TEMs) accompanied by the generation of a horizontal shear flow or jet, provided that the elevator mode is sufficiently slender. We have chosen to admit this secondary instability by increasing the wavenumber $k_e$ of the elevator modes while keeping the vertical extent of the domain fixed. This procedure is equivalent to keeping $k_e$ fixed and increasing the domain height, provided that the Rayleigh number is adjusted accordingly. Thus narrow domains favour generation of horizontal jets, while wider domains favour stable elevator modes or equivalently vertical jets.

At $Pr=1$, a further increase in the Rayleigh number destabilizes the TEM state via a subcritical Hopf bifurcation leading to an interval of coexistence between the steady TEM state and a time-dependent state that we have called a direction-reversing state. With increasing Rayleigh number, this direction-reversing state encounters a global bifurcation of Shilnikov type (Shilnikov & Shilnikov 2007), leading to a modulated travelling wave state without flow reversal. Single-mode equations that severely truncate the horizontal structure reproduce this moderate Rayleigh number behaviour well, and confirm the tame (non-chaotic) nature of the Shilnikov bifurcation for the parameter values used. At higher Rayleigh numbers, chaotic behaviour appears but is dominated by modulated travelling waves in narrow domains, while in wider domains, simulations with $Ra_{T,q}=10^{10}$ display intermittent layering and vortex dipole generation. The correspondence between mean temperature and velocity profiles resembles behaviour encountered with no-slip instead of stress-free boundary conditions in RBC with fixed temperature. In contrast, the low Prandtl number regime displays relaxation oscillations between the conduction state and the elevator mode, and exhibits quasi-periodic and then chaotic behaviour as the Rayleigh number increases. At high Prandtl numbers, the large-scale shear generated by the secondary instability is much weaker, and the flow exhibits bursting behaviour resembling that in homogeneous RBC driven by a constant temperature gradient (Borue & Orszag 1997; Calzavarini et al. 2005, 2006). These bursts are associated with a significant increase in heat transport or even intermittent stable temperature stratification. Secondary bifurcation points are shifted closer to the primary instability at lower Prandtl numbers, leading to greater fidelity of our single-mode description. The relaxation rate $\beta$ of the heat flux does not influence the late-time flow behaviour of the system.

This work opens up several directions for future work. In particular, it is crucial to analyse the corresponding dynamics in three dimensions, where large-scale shears may form in an arbitrary horizontal direction or result from the excitation of the vertical vorticity mode. Inevitably, high Rayleigh number convection generates multiple scales, and the resulting interaction between different scales in the turbulent regime represents a continuing challenge to theory. Of particular interest is the recent discovery that temperature boundary conditions play a significant role in the multiscale structure of convection even in the turbulent regime (Vieweg et al. 2021, 2022) through a mechanism that remains elusive. The details of such high Rayleigh number flows in our configuration are beyond the scope of the present investigation, but a systematic study of the Nusselt number scaling with $Ra_{T,q}$ and $Pr$ for this case is clearly desirable. Moreover, the fixed-flux formulation can be extended naturally to double-diffusive convection or rotating convection set-ups, and to reduced equations valid in the low Prandtl number limit (Spiegel 1962; Thual 1992; Lignières 1999; Garaud 2021) or in the high Prandtl number limit (Constantin & Doering 1999; Wang 2004). The usefulness of these approaches for studying finite Prandtl number dynamics in fixed-flux systems merits detailed investigation. Moreover, a systematic comparison between single-mode equations and the full 2-D equations at different Prandtl numbers will provide further quantification of the validity of the single-mode description, and potentially provide further reduction in the low Prandtl number limit.