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 Reference Ahlers, Grossmann and Lohse2009; Lohse & Xia Reference Lohse and Xia2010; Chillà & Schumacher Reference Chillà and Schumacher2012; Xia Reference Xia2013; Shishkina Reference Shishkina2021). 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. Reference Wang, Verzicco, Lohse and Shishkina2020a,Reference Wang, Chong, Stevens, Verzicco and Lohseb). Roll states are not the only type of RBC found in 2-D, however, at least with certain boundary conditions.
(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 (Reference Thompson1970) – see references in Goluskin et al. (Reference Goluskin, Johnston, Flierl and Spiegel2014, p. 363). This windy convection has features in common with strong zonal flows that arise in geophysical and astrophysical systems (Heimpel, Aurnou & Wicht Reference Heimpel, Aurnou and Wicht2005; Richards et al. Reference Richards, Maximenko, Bryan and Sasaki2006; Miyagoshi, Kageyama & Sato Reference Miyagoshi, Kageyama and Sato2010; von Hardenberg et al. Reference von Hardenberg, Goluskin, Provenzale and Spiegel2015; Kaspi et al. Reference Kaspi2018) and in tokamak plasmas (Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005). 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. (Reference Goluskin, Johnston, Flierl and Spiegel2014) and Wang et al. (Reference Wang, Chong, Stevens, Verzicco and Lohse2020a).
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 Reference Winchester, Dallas and Howell2021), and at larger 
$Ra$ only windy states can be found (Goluskin et al. Reference Goluskin, Johnston, Flierl and Spiegel2014). However, the spontaneous transitions from rolls to wind were found to be a small-domain effect by Wang et al. (Reference Wang, Chong, Stevens, Verzicco and Lohse2020a), 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. Reference Wang, Chong, Stevens, Verzicco and Lohse2020a).
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é Reference Bottin and Chaté1998; Faisst & Eckhardt Reference Faisst and Eckhardt2004; Shimizu, Kanazawa & Kawahara Reference Shimizu, Kanazawa and Kawahara2019), 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. Reference Hof, De Lozar, Kuik and Westerweel2008; Avila, Willis & Hof Reference Avila, Willis and Hof2010; Shi, Avila & Hof Reference Shi, Avila and Hof2013; Gomé, Tuckerman & Barkley Reference Gomé, Tuckerman and Barkley2020). 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. Reference Avila, Moxey, de Lozar, Avila, Barkley and Hof2011; Shi et al. Reference Shi, Avila and Hof2013; Gomé et al. Reference Gomé, Tuckerman and Barkley2020; Avila, Barkley & Hof Reference Avila, Barkley and Hof2023). 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 (Reference Rempel, Lesur and Proctor2010) simulated turbulent-to-laminar decay in magnetized Keplerian shear flow, finding that mean lifetimes increase exponentially with the magnetic Reynolds number. Linkmann & Morozov (Reference Linkmann and Morozov2015) 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 (Reference Gayout, Bourgoin and Plihon2021) report transitions between two turbulent states in wind-tunnel experiments where fluid interacts with a pendulum, and de Wit, van Kan & Alexakis (Reference de Wit, van Kan and Alexakis2022) 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. (Reference Gayout, Bourgoin and Plihon2021) suggest that lifetimes for each transition direction depend double-exponentially on the control parameter. On the other hand, de Wit et al. (Reference de Wit, van Kan and Alexakis2022) 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. Reference Winchester, Dallas and Howell2021), 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. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010; Chen et al. Reference Chen, Huang, Xia and Xi2019) and in a 3-D cylinder (Brown & Ahlers Reference Brown and Ahlers2006). Chen et al. (Reference Chen, Huang, Xia and Xi2019) 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 Reference Chandrasekhar1981), 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
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}
\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}
\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 (Reference Verzicco and Orlandi1996) and van der Poel et al. (Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015) for details.
We fix 
$Pr=10$ because this is the value for which Wang et al. (Reference Wang, Chong, Stevens, Verzicco and Lohse2020a) 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. (Reference Wang, Chong, Stevens, Verzicco and Lohse2020a) 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.
(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. Reference Avila, Willis and Hof2010). 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. Reference Wang, Chong, Stevens, Verzicco and Lohse2020a) as well as in penetrative convection (Fuentes & Cumming Reference Fuentes and Cumming2021). Power-law dependence on 
$Ra$ has been suggested also for mean switching times between different large-scale circulation patterns in RBC (Wang et al. Reference Wang, Lai, Song and Tong2018; Chen et al. Reference Chen, Huang, Xia and Xi2019).
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. (Reference Avila, Willis and Hof2010). 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. (Reference Avila, Willis and Hof2010). 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 Reference Gomé, Tuckerman and Barkley2022). 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 (Reference Goldenfeld, Guttenberg and Gioia2010) 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. (Reference Goldenfeld, Guttenberg and Gioia2010) 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 (Reference Nemoto and Alexakis2021). 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 Reference Shih, Hsieh and Goldenfeld2016), 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. Reference Goluskin, Johnston, Flierl and Spiegel2014), 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.
Supplementary material and movie
Supplementary material and movie are available at https://doi.org/10.1017/jfm.2023.875.
Acknowledgements
We thank R. Verzicco, N. Goldenfeld and H.-Y. Shih for valuable conversations, and we thank the anonymous referees for helpful comments. Prior to our investigation, D.G. discussed preliminary plans with Nigel and Hong-Yan, and Hans Johnston carried out preliminary simulations (unpublished) with a horizontal period twice the layer height. The present simulations were a separate effort, carried out on the national e-infrastructure of SURFsara, a subsidiary of SURF cooperation, the collaborative ICT organization for Dutch education and research.
Funding
D.G. was supported by the Canadian NSERC Discovery Grants Program via award numbers RGPIN-2018-04263, RGPAS-2018-522657 and DGECR-2018-0037. D.L. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 10109442).
Declaration of interests
The authors report no conflict of interest.
Appendix
Table 1 summarizes the main ensembles of simulations on which we have reported above. The full lifetime data sets for these ensembles are given in the supplementary material. In addition to the 
$\tau _m$ values estimated by averaging all lifetimes in each ensemble, the table reports alternative estimates where 
$-1/\tau _m$ is the slope of lines fit to the data in figure 2(a). The latter 
$\tau _m$ values have a best-fit scaling of 
$Ra^{4.11}$ instead of 
$Ra^{4.05}$.
For each of the main simulation ensembles, columns from left to right give the Rayleigh number, the horizontal and vertical grid resolution, the number of simulations (
$N$), the mean lifetime 
$\tau _m$ estimated by the mean of 
$\tau$ in the ensemble, 
$\tau _m$ estimated by fitting a line to data in figure 2(a), and the maximum lifetime in each ensemble. Times are in free-fall units. In all cases, 
$Pr=10$, the horizontal period is eight times the layer height, and random perturbations of the initial temperature have amplitude 
$A=10^{-4}$.
In the 
$Ra=10^7$ case, where our main ensemble size of 
$N=600$ is the largest, we have verified that our grid resolution is adequate, and that results are insensitive to how initial conditions are generated. As reported in table 1, this main ensemble was simulated with a resolution of 
$1536\times 192$ using initial conditions generated from a single snapshot by random temperature perturbations of amplitude 
$A=10^{-4}$. We performed 
$N=500$ additional simulations on the finer grid 
$2048\times 256$. The average over these lifetimes gives 
$\tau _m\approx 2806$, which agrees with the main ensemble's value 
$\tau _m\approx 2886$ to within the 
$\pm 8\,\%$ margin of the 95 % confidence interval. This indicates that the resolution 
$1536\times 192$ is adequate. We performed 
$N=800$ additional simulations on the coarser grid 
$1024\times 128$. The average over these lifetimes gives 
$\tau _m\approx 2332$. This is 19 % shorter than 
$\tau _m$ from the main ensemble, which indicates that 
$1024\times 128$ is under-resolved and produces bias towards shorter lifetimes. Finally, at the same resolution as the main ensemble, we performed 
$N=300$ simulations whose initial temperature perturbations had the smaller amplitude of 
$A=10^{-5}$. The resulting lifetimes are exponentially distributed, as in all other cases, and their average gives the very similar estimate of 
$\tau _m\approx 2864$.
The same simulations used to confirm that 
$\tau _m$ is robust to increased resolution and to decreased perturbation amplitude also confirm that 
$\tau _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 
$\tau _m$ estimates.
Open access
