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.


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 et al. 2009;Lohse & Xia 2010;Chillà & Schumacher 2012;Xia 2013;Shishkina 2021).In the two-dimensional (2D) 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) below.Roll states are seen in simulations with all combinations of noslip or stress-free velocity boundary conditions, and fixed-temperature or fixedflux thermal boundary conditions, although the boundary conditions affect what width-to-height ratios the rolls can have (Wang et al. 2020b,a).Roll states are not the only type of RBC found in 2D, however, at least with certain boundary conditions.
For 2D 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 2D simulations of various convection models at least as early as Thompson (1970) -see references in Goluskin et al. (2014, page 363).This windy convection has features in common with strong zonal flows that arise in geophysical and astrophysical systems (Richards et al. 2006;Miyagoshi et al. 2010;von Hardenberg et al. 2015;Heimpel et al. 2005;Kaspi 2018) and in tokamak plasmas (Diamond et al. 2005).Although such applications have additional important physics, the 2D RBC model may provide insight as an especially simple system in which convection drives strong zonal flow.Systematic parameter studies of windy convection in 2D 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 P r = ν/κ, where κ and ν are the fluid's thermal diffusivity and kinematic viscosity, respectively.The aspect ratio Γ of the 2D 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 Γ = 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 et al. 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 Γ = 128.When Γ 4, roll states appear stable for all combinations of P r ∈ [1, 100] and Ra ∈ [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 2D turbulence, with the windy state being only metastable, whereas roll states are apparently stable.We fix (Γ, P r) = (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 et al. 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 et al. 2010;Shi et al. 2013;Gomé et al. 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 et al. 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 et al. (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 et al. (2021) report transitions between two turbulent states in wind tunnel experiments where fluid interacts with a pendulum, and de Wit et al. (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 side walls there is no windy state, but switching occurs between different large-scale circulation patterns in a 2D or quasi-2D square (Sugiyama et al. 2010;Chen et al. 2019) and in a 3D 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 section 3, followed by conclusions in section 4. The appendix describes additional computations that verify the robustness of our results.

Simulation methods
Rayleigh-Bénard convection can be modeled 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 α.In terms of a dimensionless velocity field u(x, z, t), temperature field T (x, z, t) and pressure field p(x, z, t), the equations are The Rayleigh number is Ra = gαh 3 ∆/κν, where g is gravitational acceleration in the −z direction, h is the layer height, and ∆ is the imposed temperature difference between the top and bottom boundaries.Here e z is the unit vector in the z direction.We have scaled length by h, so the dimensionless 2D domain is (x, z) ∈ [0, Γ ] × [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αh∆) 1/2 .For the boundary conditions, stress-free conditions on the velocity vector u = (u, w) are imposed at both boundaries by w = 0 and ∂u/∂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 P r = 10 because this is the value for which Wang et al. (2020a) carried out a parameter study of Γ and Ra.Based on this study we fix Γ = 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 × 10 6 and 2.25 × 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−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.

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 = w 2 V / u 2 V , where • 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 τ 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 τ 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 τ at which a simulation undergoes a transition, we calculate the fraction of simulations that survive longer than τ .Figure 2(a) shows these fractions plotted versus τ , 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 τ m , the probability of a chosen ensemble member surviving past time t is exactly S(t) = e −t/τm .The straight lines in figure 2(a) show this S(t) for the τ m values estimated at the various Ra, meaning the lines have slopes of −1/τ m .Each τ m is estimated simply as the average over all lifetimes τ in an ensemble, which is possible because we have run every simulation until it undergoes a transition.Estimating τ m instead from the slope of the data in figure 2(a) gives similar values, as reported in the appendix.
Figure 2(b) shows how the mean lifetimes τ m vary with Ra.Both axes are logarithmic, so a linear trend corresponds to a power law scaling of τ m with Ra.The best-fit line gives the scaling τ m ≈ c Ra 4.05 .Error bars plotted for each τ m estimate are ±1.96τm / √ N , which is the 95% confidence interval for an N -member ensemble from an exponential random process (cf.Avila et al. 2010).The bestfit scaling exponent of 4.05 is indistinguishable from 4, according to the 95% confidence interval ±0.14 that the MATLAB function confint estimates from our τ 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.

Discussion and conclusions
This first investigation of the lifetime statistics of windy convection in 2D stressfree RBC has been specific to flow with P r = 10 and a horizontal period 8 times the layer height.Over the studied range of 9 × 10 6 Ra 2.25 × 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 → ∞, 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 doubleexponential 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é et al. 2022).Another possibility is a laboratory experiment, which is well suited to long time spans, but 2D 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 et al. (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 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 modeling of windy convection is as a stochastic predator-prey system.A model of this type has been proposed for pipe flow (Shih et al. 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.
The same simulations used to confirm that τ m is robust to increased resolution and to decreased perturbation amplitude also confirm that τ m is robust to the flow snapshot from which initial conditions are generated by random perturbations.This is because the ensembles with medium and high resolutions with A = 10 −4 used different snapshots, and so did the ensemble with medium resolution and A = 10 −5 .As reported above, all three snapshots led to similar τ m estimates.

Figure 1 :
Figure 1: (a) Time series (in free-fall time units) of the Reynolds number ratio Rez/Rex for one simulation with (Γ, P r, Ra) = (8, 10, 10 7 ).A transition from the windy state to the roll state occurs around time τ ≈ 286.(b-e) Temperature fields at the four instants t = 263.1,286.2, 315, 376.5 indicated in panel a.The supplementary material includes a movie of the temperature field.

Figure 2 :
Figure 2: (a) Symbols show, for each lifetime τ of a windy state measured in free-fall times, the fraction of simulations at the same Ra with longer lifetimes.(Some τ are beyond the plotted timespan.)For the mean lifetime τm at each Ra, a solid line shows the survival probability S(t) = e −t/τm .(b) Mean lifetimes τm of windy convection ( ) for each Ra.Error bars show 95% confidence intervals (see text).The best-fit power law scaling ( ) is τm ≈ c Ra 4.05 with c = 1.26 × 10 −25 .Also shown is the exponential relation obtained by linearly fitting log τ to Ra ( ), 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 × 10 −23 .