Lifetimes of metastable windy states in two-dimensional Rayleigh–Bénard convection with stress-free boundaries

Abstract

Two-dimensional horizontally periodic Rayleigh–Bénard convection between stress-free boundaries displays two distinct types of states, depending on the initial conditions. Roll states are composed of pairs of counter-rotating convection rolls. Windy states are dominated by strong horizontal wind (also called zonal flow) that is vertically sheared, precludes convection rolls and suppresses heat transport. Windy states occur only when the Rayleigh number $Ra$ is sufficiently above the onset of convection. At intermediate $Ra$ values, windy states can be induced by suitable initial conditions, but they undergo a transition to roll states after finite lifetimes. At larger $Ra$ values, where windy states have been observed for the full duration of simulations, it is unknown whether they represent chaotic attractors or only metastable states that would eventually undergo a transition to roll states. We study this question using direct numerical simulations of a fluid with a Prandtl number of 10 in a layer whose horizontal period is eight times its height. At each of seven $Ra$ values between $9\times 10^6$ and $2.25\times 10^7$ we have carried out 200 or more simulations, all from initial conditions leading to windy convection with finite lifetimes. The lifetime statistics at each $Ra$ indicate a memoryless process with survival probability decreasing exponentially in time. The mean lifetimes grow with $Ra$ approximately as $Ra^4$. This analysis provides no $Ra$ value at which windy convection becomes stable; it might remain metastable at larger $Ra$ with extremely long lifetimes.

1. Introduction

In the Rayleigh–Bénard convection (RBC) model, buoyancy-driven flow in a fluid layer is sustained by a destabilizing temperature drop from the bottom boundary to the top one. This system has been studied in laboratory experiments for over a century, and in recent decades it has also been the subject of many direct numerical simulations (DNS) in two and three dimensions with various boundary conditions for the velocity and temperature fields (Ahlers, Grossmann & Lohse 2009; Lohse & Xia 2010; Chillà & Schumacher 2012; Xia 2013; Shishkina 2021). In the two-dimensional (2-D) case, RBC most often forms convection rolls of alternating rotation direction, with a hot plume rising or a cold plume falling between adjacent rolls. A temperature snapshot for one such pair of rolls is shown in figure 1(e). Roll states are seen in simulations with all combinations of no-slip or stress-free velocity boundary conditions, and fixed-temperature or fixed-flux thermal boundary conditions, although the boundary conditions affect what width-to-height ratios the rolls can have (Wang et al. 2020a,b). Roll states are not the only type of RBC found in 2-D, however, at least with certain boundary conditions.

Figure 1. (a) Time series (in free-fall time units) of the Reynolds number ratio $Re_z/Re_x$ for one simulation with $(\varGamma,Pr,Ra)=(8,10,10^7)$. A transition from the windy state to the roll state occurs around time $\tau \approx 286$. (b–e) Temperature fields at the four instants $t=263.1,286.2,315,376.5$ indicated in (a). The supplementary material available at https://doi.org/10.1017/jfm.2023.875 includes a movie of the temperature field.

For 2-D RBC that is horizontally periodic and subject to stress-free velocity conditions at the top and bottom boundaries, some simulations have displayed a flow state dominated by a horizontal mean wind. The wind's strong vertical shear suppresses heat transport and precludes convection rolls. Figure 1(b) shows an example of such windy convection, sheared so that cold plumes move rightward along the top, while hot plumes move leftward. (The sign of the shear is an arbitrary breaking of symmetry.) The general phenomenon of windy convection, also called zonal flow or shearing convection, has been seen in 2-D simulations of various convection models at least as early as Thompson (1970) – see references in Goluskin et al. (2014, p. 363). This windy convection has features in common with strong zonal flows that arise in geophysical and astrophysical systems (Heimpel, Aurnou & Wicht 2005; Richards et al. 2006; Miyagoshi, Kageyama & Sato 2010; von Hardenberg et al. 2015; Kaspi et al. 2018) and in tokamak plasmas (Diamond et al. 2005). Although such applications have additional important physics, the 2-D RBC model may provide insight as an especially simple system in which convection drives strong zonal flow. Systematic parameter studies of windy convection in 2-D RBC were carried out by Goluskin et al. (2014) and Wang et al. (2020a).

The parameters in the equations modelling RBC can be reduced to two dimensionless numbers: the Rayleigh number $Ra$ that is proportional to the imposed temperature drop across the layer, and the Prandtl number $Pr=\nu /\kappa$, where $\kappa$ and $\nu$ are the fluid's thermal diffusivity and kinematic viscosity, respectively. The aspect ratio $\varGamma$ of the 2-D domain is the ratio of the horizontal period to the layer height. At $Ra$ just above the finite value where convection sets in, only roll states exist. For the small horizontal period $\varGamma =2$, roll states become unstable as $Ra$ is raised. There is then a narrow $Ra$ range where the flow seems to switch indefinitely between roll states and windy states with either wind direction (Winchester, Dallas & Howell 2021), and at larger $Ra$ only windy states can be found (Goluskin et al. 2014). However, the spontaneous transitions from rolls to wind were found to be a small-domain effect by Wang et al. (2020a), who simulated flows with horizontal periods as large as $\varGamma =128$. When $\varGamma \ge 4$, roll states appear stable for all combinations of $Pr\in [1,100]$ and $Ra\in [10^6,10^9]$ at which simulations were performed, never spontaneously undergoing a transition to windy convection. Windy states were also found in these larger domains at sufficiently large $Ra$; some initial conditions lead to roll states and others to windy states. When $Ra$ is barely large enough to find windy states, these states are transient and eventually undergo a transition to roll states.

In this work we study the spontaneous transition from windy states to roll states. Both states could be called 2-D turbulence, with the windy state being only metastable, whereas roll states are apparently stable. We fix $(\varGamma,Pr)=(8,10)$ and seven different $Ra$ values. At each $Ra$ we produce an ensemble of at least 200 simulations from slightly different initial conditions. Every simulation begins with windy convection but eventually undergoes a transition to a roll state with a single pair of rolls. Roll states with multiple pairs of rolls do not seem to arise from windy states, although they can develop from different initial conditions (Wang et al. 2020a).

Many other fluid systems also display metastable turbulence. Particularly well studied is the spatially localized turbulence in parallel shear flows, which decays to the laminar state at transitional values of the Reynolds number $Re$. Laboratory experiments and DNS have shown that localized ‘puffs’ in pipe flow and ‘spots’ in planar Couette flow and channel flow have survival probabilities that decrease exponentially in time (Bottin & Chaté 1998; Faisst & Eckhardt 2004; Shimizu, Kanazawa & Kawahara 2019), similar to what we report below for windy convection. The mean lifetime of a puff or spot, averaged over a large number of instances at each $Re$, appears to increase double-exponentially with $Re$ (Hof et al. 2008; Avila, Willis & Hof 2010; Shi, Avila & Hof 2013; Gomé, Tuckerman & Barkley 2020). This trend alone does not suggest that shear turbulence becomes truly stable at large $Re$, but a puff or spot has another possible fate: it can split in half, leading to two full-size puffs or spots. While the mean decay time increases superexponentially as $Re$ is raised, the mean splitting time decreases superexponentially. The $Re$ value at which splitting time crosses below decay time has been identified as the onset of sustained turbulence in pipe, Couette and channel flows (Avila et al. 2011; Shi et al. 2013; Gomé et al. 2020; Avila, Barkley & Hof 2023). In relation to the RBC model studied here, decay of a puff or spot is akin to decay of a windy state, but splitting has no analogue. Lacking a mechanism like splitting, if the mean lifetime of wind in RBC remains finite as $Ra$ is raised, then windy convection would not become truly stable, although it could be metastable with extremely long lifetimes.

Decay of metastable turbulence has also been studied in various systems beyond parallel shear flows. Rempel, Lesur & Proctor (2010) simulated turbulent-to-laminar decay in magnetized Keplerian shear flow, finding that mean lifetimes increase exponentially with the magnetic Reynolds number. Linkmann & Morozov (2015) simulated turbulent-to-simple-flow decay in isotropic turbulence forced by negative damping at large scales. They find that mean lifetimes increase superexponentially with that system's Reynolds number, as with puff or spot decay in shear flows, but there is no analogue of the splitting mechanism.

Several systems display transitions from one turbulent state to another, as in our present study, rather than decay to a simple state. Gayout, Bourgoin & Plihon (2021) report transitions between two turbulent states in wind-tunnel experiments where fluid interacts with a pendulum, and de Wit, van Kan & Alexakis (2022) report transitions into and out of a vortex condensate state in simulations of body-forced turbulence that is triply periodic with a small period in one direction. In both studies the direction of the transition depends on the control parameter. Gayout et al. (2021) suggest that lifetimes for each transition direction depend double-exponentially on the control parameter. On the other hand, de Wit et al. (2022) suggest that lifetimes diverge at a finite critical value of the control parameter, with the critical value and the rate of divergence differing between the two transition directions.

Metastable turbulence can manifest as switching back and forth between turbulent states, as opposed to permanent disappearance of one state. The present RBC configuration can switch between windy states and roll states in small domains (Winchester et al. 2021), as mentioned above, but lifetime statistics have not been studied. In RBC with sidewalls there is no windy state, but switching occurs between different large-scale circulation patterns in a 2-D or quasi-2-D square (Sugiyama et al. 2010; Chen et al. 2019) and in a 3-D cylinder (Brown & Ahlers 2006). Chen et al. (2019) suggest that mean switching times increase as a power of $Ra$.

The rest of this paper is organized as follows. Section 2 describes the governing equations and our method for simulating an ensemble of flows with transient windy convection. Results are presented in § 3, followed by conclusions in § 4. The Appendix describes additional computations that verify the robustness of our results.

2. Simulation methods

Rayleigh–Bénard convection can be modelled by the Boussinesq equations (Chandrasekhar 1981), in which the fluid's velocity is divergence-free, and buoyancy force in the vertical $z$ direction is created by the fluid's linear thermal expansion coefficient $\alpha$. In terms of a dimensionless velocity field $\boldsymbol {u}(x,z,t)$, temperature field $T(x,z,t)$ and pressure field $p(x,z,t)$, the equations are

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

(2.2)\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u} ={-}\boldsymbol{\nabla} p+ (Pr/Ra)^{1/2}\,\nabla^2\boldsymbol{u} + T{\boldsymbol{e}_z}, \end{gather}

(2.3)\begin{gather}\frac{\partial T}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} T = (PrRa)^{{-}1/2}\nabla^2 T. \end{gather}

The Rayleigh number is $Ra=g\alpha h^3\varDelta /(\kappa \nu )$, where $g$ is gravitational acceleration in the $-z$ direction, $h$ is the layer height and $\varDelta$ is the imposed temperature difference between the top and bottom boundaries. Here $\boldsymbol {e}_z$ is the unit vector in the $z$ direction. We have scaled length by $h$, so the dimensionless 2-D domain is $(x,z)\in [0,\varGamma ]\times [0,1]$, with the horizontal $x$ direction being periodic. The dimensionless time $t$ has been scaled by the free-fall time $h/U_f$, where $U_f=(g\alpha h \varDelta )^{1/2}$. For the boundary conditions, stress-free conditions on the velocity vector $\boldsymbol {u}=(u,w)$ are imposed at both boundaries by $w=0$ and $\partial u/\partial z=0$, and the dimensionless temperatures imposed are $T|_{z=0}=1$ and $T|_{z=1}=0$. We simulated (2.1)–(2.3) with these boundary conditions using the second-order staggered finite difference code AFiD, which has been extensively verified elsewhere; see Verzicco & Orlandi (1996) and van der Poel et al. (2015) for details.

We fix $Pr=10$ because this is the value for which Wang et al. (2020a) carried out a parameter study of $\varGamma$ and $Ra$. Based on this study we fix $\varGamma =8$ to safely avoid the spontaneous roll-to-wind transitions that occur only in small domains. Simulations are carried out at seven different $Ra$ values between $9\times 10^{6}$ and $2.25\times 10^7$. The minimum $Ra$ is large enough to induce windy convection, at least initially. The maximum $Ra$ leads to windy convection that can have very long lifetimes but eventually undergoes a transition to rolls in all simulations.

We use the following procedure to create an ensemble of initial conditions that all lead at first to windy convection. At each $Ra$ we start a simulation with temperature that is a random perturbation of the conductive profile $T=1-z$ and with horizontal velocity that is sheared as $u=2(z-\tfrac{1}{2})$, where we anticipate that developed flows will have order-unity velocities in our free-fall units. After windy convection develops but before it undergoes a transition to a roll state, we arbitrarily choose one snapshot of the flow. Results are not sensitive to the choice of snapshot (cf. the Appendix). For the snapshot chosen at each $Ra$, every initial condition in the ensemble is generated by perturbing the temperature at all interior grid points with pseudorandom numbers drawn uniformly from the interval $[-A,A]$. Our main results all use perturbation amplitude $A=10^{-4}$. The Appendix reports additional tests confirming that our results are not sensitive to increasing the grid resolution or decreasing the perturbation amplitude.

3. Results

Transitions from the windy state to the roll state occurred in all of our simulations. Each transition was detected as in Wang et al. (2020a) by using the vertical-to-horizontal ratio of Reynolds numbers, $Re_z/Re_x=\sqrt {\langle w^2\rangle _V/\langle u^2\rangle _V}$, where $\langle \boldsymbol {\cdot }\rangle _V$ denotes a volume average. Figure 1 shows one such time series from a simulation where the transition occurred relatively quickly (panel a), along with temperature fields before, during and after the transition (panels b–e). To precisely define the time $\tau$ of a transition, at each $Ra$ we time-averaged $Re_z/Re_x$ in the windy state and in the roll state to find the mean value of each state, then we used the average of these two values as the transition threshold. The first time $Re_z/Re_x$ crosses this threshold defines the lifetime $\tau$ of the simulation. Results are insensitive to the exact threshold; if we instead use the roll-state mean value of $Re_z/Re_x$ as the threshold, the ensemble-averaged lifetime increases by less than 2 % at the smallest $Ra$ (when lifetimes are shortest), and this percentage of increase is smaller at larger $Ra$.

To examine lifetime statistics of windy convection within each ensemble, for every time $\tau$ at which a simulation undergoes a transition, we calculate the fraction of simulations that survive longer than $\tau$. Figure 2(a) shows these fractions plotted vs $\tau$, with a different data series for each $Ra$. Each plotted series is close to a straight line, and the vertical axis scale is logarithmic, so this indicates that the fractions decrease exponentially in time. Such exponential decrease suggests that the wind-to-roll transitions behave as a memoryless random process, as in most other studies of metastable turbulence recalled in the introduction. For a memoryless random process with mean lifetime $\tau _m$, the probability of a chosen ensemble member surviving past time $t$ is exactly $S(t)=e^{-t/\tau _m}$. The straight lines in figure 2(a) show this $S(t)$ for the $\tau _m$ values estimated at the various $Ra$, meaning the lines have slopes of $-1/\tau _m$. Each $\tau _m$ is estimated simply as the average over all lifetimes $\tau$ in an ensemble, which is possible because we have run every simulation until it undergoes a transition. Estimating $\tau _m$ instead from the slope of the data in figure 2(a) gives similar values, as reported in the Appendix.

Figure 2. (a) Symbols show, for each lifetime $\tau$ of a windy state measured in free-fall times, the fraction of simulations at the same $Ra$ with longer lifetimes. (Some $\tau$ are beyond the plotted timespan.) For the mean lifetime $\tau _m$ at each $Ra$, a solid line shows the survival probability $S(t)=e^{-t/\tau _m}$. (b) Mean lifetimes $\tau _m$ of windy convection (${\square}$) for each $Ra$. Error bars show 95 % confidence intervals (see text). The best-fit power-law scaling (solid line) is $\tau _m\approx c\,Ra^{4.05}$ with $c=1.26\times 10^{-25}$. Also shown is the exponential relation obtained by linearly fitting $\log \tau$ to $Ra$ (blue dashed line), which does not fit the data. The inset shows the same plot with the vertical axis compensated by $Ra^{4.05}$; the range of this axis is $10^{-23}$ to $2\times 10^{-23}$.

Figure 2(b) shows how the mean lifetimes $\tau _m$ vary with $Ra$. Both axes are logarithmic, so a linear trend corresponds to a power-law scaling of $\tau _m$ with $Ra$. The best-fit line gives the scaling $\tau _m\approx c\,Ra^{4.05}$. Error bars plotted for each $\tau _m$ estimate are $\pm 1.96\tau _m/\sqrt {N}$, which is the 95 % confidence interval for an $N$-member ensemble from an exponential random process (cf. Avila et al. 2010). The best-fit scaling exponent of $4.05$ is indistinguishable from 4, according to the 95 % confidence interval $\pm 0.14$ that the MATLAB function confint estimates from our $\tau _m$ values. Figure 2(b) shows the line of this best-fit scaling law, along with the curve of an exponential fit. The exponential curve clearly does not fit the data well, and any fits with double-exponential or divergent dependence on $Ra$ would be even worse.

4. Discussion and conclusions

This first investigation of the lifetime statistics of windy convection in 2-D stress-free RBC has been specific to flow with $Pr=10$ and a horizontal period eight times the layer height. Over the studied range of $9\times 10^6\le Ra\le 2.25\times 10^7$, the wind-to-roll transition behaves as a memoryless random process, and the ratio of mean wind lifetime to free-fall time scales approximately like $Ra^4$. The corresponding ratio of mean lifetime to diffusive time scales like $Ra^{3.5}$. While $Ra$ varies only by a factor of 2.5 over the range of our simulations, the mean lifetimes vary by a factor of almost 40. This range is large enough to see clearly that the mean lifetimes are better described by power-law scaling with $Ra$ than by exponential or superexponential dependence on $Ra$. Whether this scaling differs at other Prandtl numbers and horizontal periods remains to be studied; both parameters are known to affect the $Ra$ values at which windy states occur in RBC (Wang et al. 2020a) as well as in penetrative convection (Fuentes & Cumming 2021). Power-law dependence on $Ra$ has been suggested also for mean switching times between different large-scale circulation patterns in RBC (Wang et al. 2018; Chen et al. 2019).

If power-law scaling of mean lifetimes continues as $Ra\to \infty$, this would mean that windy convection is never truly stable but is metastable with such long lifetimes at large $Ra$ as to be effectively stable. However, the history of transitional shear flow studies suggests caution when extrapolating our findings to asymptotically long lifetimes; for turbulent puffs and slugs, the double-exponential dependence of lifetimes on $Re$ is technically hard to determine, and earlier studies with less data and smaller domains suggested other dependence on $Re$ – see discussion in Avila et al. (2010). Extending our DNS approach to larger $Ra$ would be very expensive; nearly 0.47 million CPU hours were needed to observe transitions in all 200 of our highest-$Ra$ simulations. Somewhat larger $Ra$ could be reached if simulations are truncated at a maximum timespan rather than waiting for every one to undergo a transition; mean lifetime estimates can account for such truncation, as described in Avila et al. (2010). At even larger $Ra$ where mean lifetimes are extremely long, one might employ a rare-event algorithm where DNS is performed selectively for cases leading to a transition while pruning others, as has been done for shear flows (e.g. Gomé, Tuckerman & Barkley 2022). Another possibility is a laboratory experiment, which is well suited to long time spans, but 2-D convection with stress-free top and bottom would be hard to approximate. For this reason and others, it is of great interest to find a laboratory model that displays windy convection.

Regarding a theoretical explanation of how mean lifetimes of windy convection depend on $Ra$, further data may allow an argument that invokes extreme value statistics. For pipe flow, Goldenfeld, Guttenberg & Gioia (2010) suggested that puff decay is triggered when the maximum-over-space turbulent intensity drops below some threshold. They note that if this maximum intensity follows a Gumbel distribution, and if the threshold decreases linearly in Re as $Re$ is raised, then the mean time before falling below the threshold would increase double-exponentially with $Re$. For RBC, the wind-to-roll transition mechanism seems the opposite: rolls might be triggered when the maximum-over-space deviation from the mean wind rises above some threshold. If this threshold increases linearly in $Ra$, then arguments analogous to Goldenfeld et al. (2010) imply that the mean lifetime would increase exponentially with $Ra$, not double-exponentially. However, if the threshold increases only logarithmically with $Ra$, then the mean lifetime would increase as a power of $Ra$, as we observe in our data. Further DNS is needed to determine whether such extreme value arguments apply to the present model, as has been largely confirmed for pipe flow by Nemoto & Alexakis (2021). In particular, are wind-to-roll transitions triggered when maximum turbulent intensity rises above a threshold? If so, which extreme value distribution governs this maximum, and how does the threshold vary with $Ra$?

Another direction for simplified modelling of windy convection is as a stochastic predator–prey system. A model of this type has been proposed for pipe flow (Shih, Hsieh & Goldenfeld 2016), where it captures double-exponentially varying mean lifetimes of not only puff decay but also splitting. This model consists of stochastic predator–prey dynamics in a spatially extended domain, with the predator representing azimuthal zonal flow in the pipe and the prey representing turbulent fluctuations. A similar predator–prey analogy has been described for windy convection (Goluskin et al. 2014), again with the zonal flow as the predator, but there is no clear analogue to the spatial localization and splitting of puffs in pipe flow. As with the extreme value arguments, closer analysis of various quantities in windy convection may be needed to formulate an insightful stochastic model.