Two-fluid single-column modelling of Rayleigh-B\'{e}nard convection as a step towards multi-fluid modelling of atmospheric convection

Multi-fluid models have recently been proposed as an approach to improving the representation of convection in weather and climate models. This is an attractive framework as it is fundamentally dynamical, removing some of the assumptions of mass-flux convection schemes which are invalid at current model resolutions. However, it is still not understood how best to close the multi-fluid equations for atmospheric convection. In this paper we develop a simple two-fluid, single-column model with one rising and one falling fluid. No further modelling of sub-filter variability is included. We then apply this model to Rayleigh-B\'{e}nard convection, showing that, with minimal closures, the correct scaling of the heat flux (Nu) is predicted over six orders of magnitude of buoyancy forcing (Ra). This suggests that even a very simple two-fluid model can accurately capture the dominant coherent overturning structures of convection.


Introduction
Despite being an important part of the global circulation and local variability, atmospheric convection is one of the weakest aspects of Numerical Weather Prediction (NWP) and climate models (Stephens et al. 2010;Sherwood et al. 2014;Stein et al. 2015;Clark et al. 2016). These difficulties are at least in part due to the "grey zone" problem: the resolution of current models is such that a typical grid spacing is neither much smaller, nor much larger, than a typical convective scale (say O(1 km) for shallow convection), meaning that neither traditional parametrizations, nor so-called "explicit convection", adequately represent the flow (Wyngaard 2004;Holloway et al. 2014;B. Zhou et al. 2014;Clark et al. 2016). If a separation of scales exists between the clouds scale(s) and the grid spacing, simplifying assumptions may be made to aid in parametrization of sub-grid processes. In atmospheric convection, traditional closures have assumed some form of balance between the large-scale forcing and the convective response, and that the area of a "grid box" taken up by cloud is small (Arakawa and Schubert 1974;Plant 2010). Both of these assumptions break down in the grey zone of current model resolutions, leading to unrealistic behaviour of models with traditional parametrizations.
At the other end of the scale, for grid spacings much smaller than the convective length scale(s), techniques of large-eddy simulation (LES) become applicable (see Mason 1994, for a review). However, true LES requires very high resolution -typically O(10 m) for the dry convective boundary layer (Sullivan and Patton 2011) -which is beyond the computational capabilities of NWP and climate models for the foreseeable future. Many current operational forecasting models (e.g. the Met Office UKV configuration of the MetUM, DWD's ICON-D2 ) use "explicit" convection, where the convection scheme is turned off. Some form of 'turbulence' scheme is still required, often an LES-like eddy viscosity/diffusivity scheme. While these perform better at grey zone resolutions, there are still undesirable effects, in particular the prediction of incorrect length scales (cloud size and inter-cloud spacing) typically larger than observed scales, even in cases where the model should be able to resolve the smaller scale (Lean et al. 2008). There is thus both a need for parametrization well into the future, and a need for new parametrization approaches in the grey zone.
Multi-fluid modelling has recently been proposed as an approach to representing convection in the grey zone (Yano 2014;Thuburn et al. 2018;Tan et al. 2018); similar equation sets are used for the modelling of multi-phase flows in engineering (e.g. Städtke 2007). In the convection context, this takes inspiration from traditional mass-flux parametrizations in splitting the fluid into multiple components, which may represent updrafts, environment, downdrafts etc.. The split is applied directly to the governing equations, which are then spatially filtered, allowing a fully 3D and time-dependent framework to be derived (Thuburn et al. 2018;Shipley et al. 2021). Neither quasi-equilibrium nor small updraft fraction are assumed in the derivation. Each "fluid" evolves according to its own prognostic equations, interacting with other fluids via the pressure gradient, and terms involving the exchange of mass, momentum, energy, and tracers. These exchange terms are the analogue of entrainment, detrainment, and cloud-base mass-flux in traditional models, and must be parametrized. Convection is inherently a part of the dynamics in this framework: there is no separate convection scheme which is called by the dynamical core.
The skewness of (joint) probability distribution functions of variables in convective flows is well known to be important (Larson et al. 2002;P. Zhu and Zuidema 2009) and is often poorly treated in first or second-order turbulence closures; one approach to modelling this variability is assuming bi-Gaussian joint probability distributions in PDF-based convective closures (Larson et al. 2012;Fitch 2019). Each Gaussian can be thought of as a different component of the fluid. A potential advantage of the multi-fluid approach is that even the simplest possible multi-fluid model, a two fluid model, intrinsically captures information about odd-order moments. It is therefore possible that the multi-fluid method can provide a better low-order approximation for flows with bimodal distributions, or large skewness.
In order to build a multi-fluid model of atmospheric convection, the multi-fluid equation set must be closed. The form of these closures directly depends on the definition of the fluid partitions (de Rooy et al. 2012;Shipley et al. 2021). For example, the single-column 2-fluid model of Thuburn et al. (2019) contains entrainment and detrainment closures designed to capture coherent structures in the convective boundary layer, whereas the closures in Cohen et al. (2020) are designed to model a second fluid in the cloud layer only. Perturbation pressure closures for the latter approach were suggested in He et al. (2020). Entrainment and detrainment closures based on velocity divergence, and a bulk viscous parametrization for the perturbation pressure, were proposed and tested in Weller et al. (2020), but the test cases used for comparison were non-turbulent, unlike the real atmosphere. All of these multi-fluid models have been single-column, and used standard atmospheric test cases (e.g. dry rising bubble, dry convective boundary layer, oceanic and continental shallow cumulus, diurnal deep convection) for verification. While prior work shows the considerable promise of the multi-fluid method, little work has been done testing the response of a specific multi-fluid scheme to a variety of forcings, or suggested how the closure constants should scale with that forcing. Such investigation could lead to more consistent results compared to tuning a model to a handful of test cases.
To gain a better understanding of the multi-fluid equations, and how some of the new closure terms affect the solution, we present a single-column model of dry Rayleigh-Bénard convection (RBC) with one rising and one falling fluid. RBC is the simplest relevant convection problem: the equations and boundary conditions are as simple as possible while still allowing for a fully turbulent convective solution. RBC has been extensively studied, and a wealth of experimental, numerical, and theoretical results make it a well-constrained starting point (Chandrasekhar 1961;Ahlers et al. 2009;Chillà and Schumacher 2012). In particular, the scaling of bulk buoyancy and momentum transport with the applied buoyancy forcing is well understood over at least ten orders of magnitude.
It is important to understand the response of the model in a fully-parametrized equilibrium setting before moving to the grey zone. This will help pin down the physics of a multi-fluid model of convection, free of the complexities -especially microphysics and phase changesof the real atmosphere.
The paper begins with an overview of Rayleigh-Bénard convection in section 2, motivating its use as a reasonable testbed for developing insights into "real-world" convection. Results from 2D direct numerical simulations (DNS) of dry RBC are presented, and shown to agree with reference results. In section 3 a multi-fluid Boussinesq equation set is presented, along with a discussion of how and why this equation set differs from previous papers on multi-fluid convection parametrization. Closures for one rising and one falling fluid which attempt to capture the large-scale overturning circulation are presented in section 3.1, and a scaling argument is presented for the magnitude of the pressure differences between the fluids. The numerical method is then described in section 4. In section 5, results of the two-fluid single-column model (section 3.1) are compared with horizontally-averaged results from the DNS (section 2.1) over a range of buoyancy forcing spanning seven orders of magnitude (10 3 ≤ Ra ≤ 10 10 ), and the sensitivity of the model to its two dimensionless closure constants is investigated. The paper concludes with a summary of its results and their relevance to convection parametrizations, and a discussion of avenues for future research.

Rayleigh-Bénard convection (RBC)
The Rayleigh-Bénard problem is the simplest fluid dynamical model of convection. First studied experimentally by Bénard (1900), the problem was given a theoretical treatment by Rayleigh (1916) which has been the basis of over a century of investigation. Rayleigh (1916) studied the motion of a Boussinesq fluid confined between two perfectly conducting horizontal plates of infinite extent, each held at a constant uniform temperature. For mathematical tractability he considered stress-free velocity boundary conditions at the plates; the no-slip case was tackled by Jeffreys (1926) and Jeffreys (1928). RBC has long been of interest to the meteorological community, being the basis of the Lorenz (1963) seminal discovery of deterministic chaos, and a key component of our understanding of convective systems (Emanuel 1994, ch. 3). Moist extensions of the model have been considered to gain insight into moist convection, though far less work has been performed on moist versions of the problem than on the dry case (Bretherton 1987;Bretherton 1988;Pauluis and Schumacher 2010;Weidauer and Schumacher 2012;Vallis et al. 2019). In this section, the classical results relevant to this paper are collected. The canonical text covering stability and the onset of convection is Chandrasekhar (1961); recent reviews covering fully turbulent convection are Ahlers et al. (2009) and Chillà and Schumacher (2012).
The setup of the Rayleigh-Bénard problem is as follows. A Boussinesq fluid is confined between two smooth, flat, horizontal plates, a fixed distance H apart. Each of these is held at a fixed buoyancy, ± ∆B/2, with no-slip, no-normal flow velocity boundary conditions. For both analytical and numerical simplicity we choose the lateral boundaries to be periodic in all fields. The motion of the fluid is described by the following Boussinesq equations of motion: Here u denotes the velocity field of the fluid; b := g(ρ ref − ρ)/ρ ref its buoyancy 1 ; P := p/ρ ref its pressure potential; ν its kinematic viscosity; κ its buoyancy diffusivity; and k is a unit vector antiparallel to gravity, defining the vertical (z) direction. All variables are defined relative to a resting, uniformly constant-density, hydrostatically-balanced pressure reference state. A diffusive nondimensionalization of equations (1)-(3) (as Chandrasekhar 1961; Emanuel 1994, but choosing a diffusive rather than viscous time-scale) by the external parameters, x := x/H ,b := b/∆B ,t := tκ/H 2 ,û := uH/κ ,P := H 2 /κν , shows that two dimensionless parameters govern the flow (the boundary conditions are given for completeness): u(ẑ = 0, 1) = 0.
1 For our desired application to atmospheric convection, ( where θ is the potential temperature, though much previous work on Rayleigh-Bénard convection (RBC) is performed in terms of temperature, using the approximation ( The equation set retains the same form.
Nondimensionalized variables are denoted by a hat, and the dimensionless parameters are defined by: The Rayleigh number, Ra, is the ratio of buoyancy forcing (∆B) to viscous diffusion (κν/H 3 ); and the Prandtl number, Pr, is the ratio of the diffusion of momentum (ν) to the diffusion of buoyancy (κ). The former can thus be seen as measure of the applied forcing in RBC, whereas the latter is an intrinsic property of the fluid. This nondimensionalization shows that any two RBC systems with the same Ra and Pr support the same solutions, i.e. are self-similar. 2 It is worth noting that this nondimensionalization specifically singles out the diffusive regime as the regime of interest, relevant for considerations of stability. For consideration of the convective solutions, a nondimensionalization based on the buoyancy forcing is more useful. This "free fall" or "free convective" scaling gives velocity and time scales U B := √ ∆B H, T B := H/∆B, and is ubiquitous in the CBL literature (where U B is denoted w * , see, e.g., Garratt (1994). Such a scaling also gives an a priori estimate for the Reynolds number, Re ∝ Ra 1/2 Pr −1/2 . 3 This approximate Re(Ra) scaling is observed for the regimes applicable to this paper.
The equation set (4)-(8) has a unique stationary zero-flow solution, with a linear buoyancy gradient between the plates and a quadratic pressure profile: This solution is both linearly and nonlinearly unstable to perturbations if and only if the Rayleigh number exceeds a critical value, Ra c ; importantly, the stability does not depend on the Prandtl number (see, for instance, Chandrasekhar 1961;Joseph 1966;Lindsay and Straughan 1990). Below Ra c , solutions are purely diffusive; above Ra c , a circulation develops which increases the heat transport. This circulation can either be steady, periodic, quasiperiodic, or turbulent, depending on the governing parameters (Ra, Pr). The precise value of Ra c depends on the velocity boundary conditions at the top and bottom boundaries, but not on the dimensionality of the domain; for our chosen no-slip conditions, Ra c ≈ 1708, and the wavelength of the most unstable mode is λ c ≈ 2.02H (Chandrasekhar 1961 , table 3). The domain-and time-averaged dimensionless buoyancy flux is given by the Nusselt number: which is the ratio of the actual buoyancy flux to the buoyancy flux of the purely diffusive solution. Averaging the buoyancy equation (5) over a horizontal plane and over time (denoted . . . A,t ) shows that the Nusselt number is independent of height in a statistically stationary flow.
Exact results for the domain-and time-averaged kinetic and thermal dissipation rates, ε u and ε b , are given by (Siggia 1994;Chandrasekhar 1961, appendix 1): 2 A third parameter, the aspect ratio of the domain, Γ := L/H, enters via the lateral boundary conditions; however, the dependence upon the aspect ratio is generally weak so long as Γ > 1 -see Ahlers et al. 2009, section 3E; also Johnston and Doering 2009;Bailon-Cuba et al. 2010;Q. Zhou et al. 2012 -and the dependence is weaker for periodic boundaries than for rigid boundaries.
3 The a priori scaling for the Nusselt number that this predicts, Nu ∝ Ra 1/2 -the so-called "ultimate scaling" -is steeper than observed to date in experimental or numerical dry RBC, because the non-turbulent surface layers next to the boundaries prevent a thermal shortcut.
Here the "double dot product" denotes the complete contraction of two rank-two tensors, following the convention A : B := a,b A ab B ab . Thus the vertical buoyancy flux is the only quantity that characterizes the stationary-state global energetic response of the system to the applied forcing (Ra, Pr) 4 . The statistically steady-state Rayleigh-Bénard problem can then be framed as asking the question: if we apply a buoyancy forcing Ra to a Boussinesq fluid characterized by Pr, what is the resulting Nu? Scaling theories for Nu as a function of Ra and Pr are well-developed, and there is good agreement between the theory and numerical and experimental results until at least Ra = 10 11 for Pr = O(1) (Ahlers et al. 2009;Chillà and Schumacher 2012). It is therefore a strong test of any dynamical low-order model of RBC to reproduce these scalings.

2D direct numerical simulation of RBC
To provide a reference "truth" for later sections in the paper, results from two dimensional direct numerical simulations of Rayleigh-Bénard convection over a wide range of Ra are presented. These simulations also serve to illustrate the phenomenology of RBC, and to indirectly validate the numerical methods via comparison with reference results.
While the restriction to two dimensions may seem like too great a simplification, global and large-scale results of Rayleigh-Bénard convection in two and three dimensions are remarkably similar so long as the Prandtl number is not too small.
The classical results regarding the critical Rayleigh number, critical wavelength, and onset of convection are unaffected (see Chandrasekhar 1961, ch. 2; though not explicitly stated, the stability analysis does not depend on the dimensionality of the domain). After the onset of convection, for O(1) Pr and greater, the scalings of global parameters such as the Nusselt and Reynolds numbers, as well as the boundary layer depths, are virtually the same in 2D as in 3D (although the magnitudes differ slightly)see Schmalzl et al. (2004). Many theoretical analyses of the problem have either included two dimensions as a special case, or actually assumed only two dimensions, the successful Grossmann and Lohse (2000) scaling theory for the Nusselt and Rayleigh numbers being a prime example of the latter. Therefore we choose to perform 2D simulations, given the similarity between 2D and 3D results and the vastly reduced computational requirements for 2D calculations.
Our simulation suite runs from fully diffusive (Ra 10 2 ) to well into the turbulent regime (Ra 10 10 ). Rayleigh numbers have been chosen such that there is at least one simulation per factor of ten of Ra, with extra simulations run in the vicinity of Ra c . The Prandtl number is fixed to be Pr = 0.707, the value for dry air at STP. Reviews of RBC suggest that qualitative results remain similar so long as the asymptotic range of Pr is the same, i.e. Pr = O(1) rather than Pr → 0 or ∞ (Ahlers et al. 2009;Chillà and Schumacher 2012). In particular, the scaling exponent Nu ∝ Ra β is not strongly Prandtl-number dependent.

Choice of resolution
• By "resolution", ∆ r , we mean the smallest length scale at which structures of the flow are well captured by the model.
• By "filter scale", ∆ f , we mean the length scale(s) associated with any filter applied to the flow, whether to the solutions or to the governing equations.
• By "grid scale" (alternatively, "grid length" or "grid spacing"), ∆ g , we mean the actual distance between points (or cell centres) within a discretized model.
A direct numerical simulation of a fluid ("DNS") must "resolve" all dynamically relevant scales of the fluid flow in order to justify the assumption that no small-scale processes need to be parametrized. But there are various metrics by which we can test whether a flow is "resolved".
To fully resolve a turbulent flow, the grid spacing must resolve at least a factor of ten into the viscous subrange (Kerr 1985), which is very computationally expensive. However, to get the majority of the statistics right the requirements are less extreme: the Kolmogorov dissipation length, η := H(Pr 2 /ε u ) 1/4 , must be resolved (Grötzbach 1983). Within fully-developed turbulence in the bulk of the fluid the exact result for the global kinetic energy dissipation rate, (12), may be used to estimate the smallest dynamically relevant scale: Towards the boundaries, the kinetic and thermal boundary layers must be resolveddissipation is typically higher in these regions, reducing the smallest dynamically relevant length scale. Shishkina et al. (2010) estimated local dissipation lengths based on dissipation rates defined within the boundary layers, using these to estimate the minimum number of points N u , N b required within each boundary layer (thickness δ u , δ b ) in order to adequately resolve the flow. This estimate is for 10 6 < Ra < 10 10 , so for Ra ≤ 10 6 we use the values of N u , N b estimated for Ra = 10 6 . Note that this extra resolution is only required in the vertical direction.
At any point in the flow the smallest of {η, δ b /N b , δ u /N u } must be resolved. Collecting the results of Grötzbach (1983) and Shishkina et al. (2010), the grid spacing is required to satisfy ∆x η < 2η, ∆x b < δ b 0.35 Ra 0.15 , ∆x u < δ u 0.31 Ra 0.15 to be adequate to resolve each respective scale.
To make use of the resolution requirements, the boundary layer thicknesses must be estimated. Since the centre of the domain will be statistically well-mixed after the onset of convection, we must have δ b /H ∼ 1/2 Nu. For the parameter regimes of this study, Nu ∼ Ra 2/7 and so δ b /H ∼ Ra −2/7 (Castaing et al. 1989;Shraiman and Siggia 1990;Ahlers et al. 2009). Prandtl-Blasius boundary layer theory suggests that the kinetic boundary layer thickness should scale as δ u /H ∼ Ra −1/4 , and δ u < δ b is expected over the entire Rayleigh number range here considered (Ahlers et al. 2009, fig. 3). To estimate the prefactors, an over-resolved simulation with ∆x/H = ∆z/H = 0.01 was run at Ra = 10 5 , finding δ u ≈ 0.56 Ra −1/4 , δ b ≈ 2.8 Ra −2/7 ; these prefactors do indeed ensure that δ u < δ b for the Rayleigh number regime of the study.
For each Ra we construct an orthogonal, rectangular grid such that the grid spacing is always smaller than the smallest of these length scales. This grid consists of, in the z-direction: a uniform grid with spacing ∆z (0) = ∆x u for 0 ≤ |z − z boundary | ≤ δ u ; a uniform grid with spacing ∆z (1) : In the horizontal direction, grid spacing is uniformly equal to 2η throughout the domain. Details of the grid for each simulation are given in table 1.
In principle, we could directly check that the resolution is sufficient post-hoc by refining the grid and re-computing all of the statistics; if they do not change as the resolution increases, then the lower resolution "fully resolves" the flow. In practice, for this paper we note that the grid spacings of our simulations are comparable to those in similar DNS of 2D RBC (e.g. Johnston and Doering 2009). Details of the numerical method are given in section 4 as a special case of the multi-fluid solver.

Calculation of Nu, Re, δ b
The Nusselt number, Reynolds number, and boundary layer depths are calculated as follows:   (14)) and the kinetic boundary layer thickness, δ u /H ≈ 0.56 Ra −1/4 . The Ra = 10 9 and 10 10 simulations were spun up on a coarser grid (the Ra = 10 8 grid), then after reaching equilibrium the grid was refined. The simulation time on the finer grid is given, followed by, in parentheses, the total simulation time on both grids for that Rayleigh number.
Nu: The most direct way of calculating Nu is to integrate the (dimensionless) heat flux over the entire domain, then take a time average: Nu = wb − ∂b/∂z V,t . However, if the flow is statistically stationary, then the time-averaged horizontally averaged (dimensionless) heat flux is independent of height, so calculating the time-averaged vertical buoyancy gradient averaged over the top and bottom boundaries gives a second estimate, Nu w := − ∂b/∂z A,t;z=0,H . The equivalence of these two expressions for Nu provides an extra check for the statistical steadiness of the numerical solutions. Another check for statistical stationarity is provided via the kinetic and thermal dissipation rates (calculated using equations (12)- (13)). Thus for a statistically stationary state, convergence of Nu = Nu w = ε b = 1 + ε u /Ra is required.
Re: The calculation of a Reynolds number based on the definition Re := U L/ν requires the choice of a velocity scale and a length scale. For RBC, the only length scale we can reasonably choose for a bulk Reynolds number must be the domain height H, as this is the only external length scale in the problem. However, what is a reasonable representative velocity scale, U ? Several possible choices are suggested in Kerr (1996) and Ahlers et al. (2009); we shall consider velocity scales based on the turning points of the velocity variance profile: An a priori estimate of Re can be found by assuming free-convective scaling, U = U B := √ ∆B H, implying Re = ∆B H 3 /ν 2 = Ra 1/2 Pr −1/2 . δ b : If the flow is statistically stationary, the buoyancy will be well-mixed in the interior of the domain, the time-averaged buoyancy profile must be approximately constant outside of the boundary layers, and approximately linear within due to the fixed buoyancy boundary conditions. Thus one measure of the thermal boundary layer thickness is Following Kerr (1996), we also estimate the thermal boundary layer thickness from the locations of the maxima of the buoyancy variance profile: Both the upper and lower boundary layer thicknesses should be the same.
The above time averages are calculated over at least 5 eddy turnover times (T e ≈ 4 T B ). Timeaverages are also calculated over twice and three times this minimum averaging time, and all simulations show convergence between the averages taken over these three different times. The total simulation time for each Rayleigh number is given in table 1.

The relevance of RBC to atmospheric flows
While RBC is a valuable test problem in its own right, it is worth considering similarities with and differences from atmospheric flows, in particular the dry atmospheric convective boundary layer (CBL). Besides the complexities of moisture, the dry RBC problem differs from even dry atmospheric convection in a few important ways. Firstly, the Boussinesq approximation is of questionable validity even on the scale of the atmospheric boundary layer; in practice however, it has long been used in the LES community with excellent results (Sullivan and Patton 2011, e.g. ). Furthermore, the Boussinesq form has been used to facilitate analysis; experiments using a non-Boussinesq (fully compressible) version of the same code show little qualitative or quantitative differences from their Boussinesq counterparts.
Secondly, the lower boundary in the CBL is neither smooth, nor uniformly heated. Recent results show that neither nonuniform heating (Bakhuis et al. 2018) nor rough boundaries (X. Zhu et al. 2019;Toppaladoddi et al. 2021) drastically change the dynamics of RBC, though the latter does tend to increase the heat flux towards the so-called "ultimate regime", equivalent to the free-convective regime which dominates discussion of scaling in the atmospheric convective boundary layer.
Thirdly, the fixed buoyancy boundary conditions are quite different to CBL conditions, where the lower boundary is closer to (and is often modelled as) a fixed buoyancy flux, and there is no fixed upper boundary for the convection (instead there is a stable atmospheric layer). In practice, solutions of RBC with fixed flux vs. fixed value boundary conditions are similar, especially in 2D (Verzicco and Sreenivasan 2008;Johnston and Doering 2009) (as are LES simulations of the CBL). It is thus only the upper boundary that introduces a major difference between RBC and the CBL. Even in that case there has been recent progress on studying modified Rayleigh-Bénard convection with the compensating heat flux provided by radiation in a layer of finite thickness (Lepot et al. 2018;Doering 2019), which the first authors note "spontaneously achieves the "ultimate" regime of thermal convection".
We thus consider the classical Rayleigh-Bénard problem to be sufficiently close to atmospheric convection to provide a reasonable testbed for investigating the behaviour of a multifluid model of turbulent convection. There remains the question of the applicable parameter regime, discussed in the next section.
2.2.1 An analogy between constant-viscosity RBC and large-eddy simulation of higher-Ra RBC.
Atmospheric flows generally involve very high Reynolds number; for example, the CBL might have depth ≈ 1000 m, and (even without a mean wind) velocities in convective updraughts ≈ 1 m s −1 . With kinematic viscosity of air ≈ 10 −5 m 2 s −1 we have Re ≈ 10 8 . In the context of RBC, this would lead to Ra ≈ 10 16 , i.e. much larger than in our simulations. Given the above considerations of resolution, a DNS of this problem is computationally impossible in 3D with current computing power and would be a challenge even in 2D. The atmospheric science community address this problem using 'Large-Eddy Simulation' (LES) as reviewed by Mason (1994). LES is based upon spatially low-pass filtering the equations of motion with a filter with characteristic scale ∆ f chosen such that the unfiltered flow remains well within the self-similar Kolmogorov 'inertial sub-range' (ISR) of scales. In this case, the sub-filter contribution to turbulent fluxes is small and can be represented by a simple eddyviscosity. In practice, the eddy-viscosity proposed by Smagorinsky (1963) has been found to give good results provided the simulation actually is well within the ISR.
In fact, Mason (1994) points out that acceptable results are obtained from a simulation of the CBL in which a constant viscosity is used at each level based upon the horizontal average of the Smagorinsky value. One might go further and suggest that the height-dependence is required primarily close to the surface and boundary-layer top where eddy length-scales are restricted. In this case, a simple view of LES is as follows. With a well-developed ISR the flow is essentially independent of Re. For sufficiently large Re we can choose an artificial larger viscosity such that the range of scales in the flow is smaller as the Kolmogorov microscale (i.e. the eddy scale at which Re = 1) is larger. The turbulent kinetic energy dissipation rate (ε) remains the same as does the flow at larger scales. Essentially the same argument applies to use of wind-tunnels with scale models.
Smagorinsky provides us with a method to estimate this artificially large viscosity, but let us take a more basic view. In the ISR the energy spectrum E(k) = K 0 ε 2 3 k − 5 3 , with K 0 a constant and k the wave number. Suppose we choose a filter scale with wave number k f , then the 'turbulent' kinetic energy (TKE) in the subfilter flow has a velocity scale U f given by Prandtl argues, by analogy with the kinetic theory of gases, that the eddy diffusivity is a product of the turbulent velocity scale and the 'mean free path' or 'mixing length'. It is possible to show this more rigorously using the dynamical equation for stress. We assume that the mixing length scales with k −1 f and hence the eddy viscosity is given by This assumption is precisely the same as stating that the eddy Reynolds number, the filter scale is proportional to the Kolmogorov microscale of the filtered flow, η f ∝ k −1 f . Indeed, if we absorb the constants in eq. (18) into η f , the Kolmogorov microscale for the filter, then f , Re f ≡ U f η f /ν f = 1 and all of the Kolmogorov scales apply. Note that with this viscosity the kinematic deviatoric stress is given by τ = U f η f ∇u + ∇u T . The spatially filtered equation for the TKE, in steady state and ignoring the transport term (both assumptions being appropriate for the homogeneous isotropic turbulence the ISR is considered to represent) leads to a simple balance between shear production and dissipation: The TKE is given by 1 2 U 2 f . The scaling above gives ε = U 3 f /η f (so the timescale for dissipation is η f /(2U f )) then this balance becomes: from which U f = η f 1 2 ∇u + ∇u T : ∇u + ∇u T 1 2 , leading to the Smagorinsky (1963) formulation of viscosity.
To give a simple example, suppose we have a convective boundary layer with H the depth of the layer, say 1000 m and convective velocity scale U = 2 m s −1 , corresponding to ∆B ≈ 4 × 10 −3 m s −2 . Then Re = U H/ν = 1.33 × 10 8 and Ra = 1.8 × 10 16 . The 'outer' mixing length is often taken to be L = 0.15H. Then a crude estimate of ε is ε = U 3 /L = 8/150 m 2 s −3 = 5.3 × 10 −2 m 2 s −3 . The Kolmogorov microscale is thus η = (ν 3 /ε) 1 4 ≈ 0.5 mm. If we choose a filter scale such that η f = 1 m, then U f = 0.376 m s −1 , the eddy viscosity is 0.376 m 2 s −1 and the Reynolds number of the whole flow is reduced to Re ≈ 5300. This should still be turbulent and is likely to be within the Reynolds number independent regime. In fact, the convective boundary layer is very amenable to LES because the large coherent structures with scales of order H dominate, and Sullivan and Patton (2011) show that "the majority of the low-order moment statistics (means, variances, and fluxes) become grid independent when the ratio z i /(Cs∆ f ) > 310". Here z i is the inversion depth (i.e. H), Cs is the Smagorinsky constant (≈ 0.2), and ∆ f essentially the grid length (so the actual filter scale is a multiple of this). This implies ∆ f < z i /(310Cs) ≈ 16 m in this case. Hence our notional 1 m resolution should be very well-converged LES.
Thus, provided the solution remains in the Re-independent turbulent regime, the relatively low Ra runs (Ra ≥ 10 8 ) may be interpreted as reasonable approximations to LES of much higher Ra (and hence Re) flows encountered in the atmosphere. Indeed, a similar argument is made by Mellado et al. (2018) in a study of Stratocumulus convection, except that they seem to have reversed the semantics of the conclusion by describing their artificially large viscosity runs as DNS. This 'DNS in a Re-independent regime' is, we would argue, more correctly described as a form of LES as its basis is precisely the same.
A slight note of caution may arise from consideration of the boundary conditions, as the turbulence length scale collapses as one approaches the boundary and buoyancy effects on turbulence become more dominant. (The same concerns apply to LES). With a fixed heat-flux boundary condition, the concern is less as the surface exchange serves merely to transport the given surface flux into the fluid where large eddies can start to transport it. In practice, our results are similar for fixed temperature and fixed heat flux boundary conditions, suggesting that, so long as the thermal boundary-layer is adequately resolved the solutions remain applicable to higher Re.

Phenomenology of RBC
Direct numerical simulations of 2D, dry, Boussinesq Rayleigh-Bénard convection were performed for the range 10 2 ≤ Ra ≤ 10 10 for a fluid with Prandtl number 0.707 (the value for dry air at standard temperature and pressure). For each Ra, the fluid was initialized from the hydrostatically-balanced resting state (10), with small random perturbations to the buoyancy field |δb pert | ≤ 0.01∆B drawn from a uniform distribution. The aspect ratio of the domain was set equal to the critical wavelength: Γ = L x /L z = λ c /H ≈ 2.02 (for 10 2 ≤ Ra ≤ 10 8 , simulations were also run with Γ = 10, which gave the same results for Nu, Re etc.; therefore only the smaller aspect ratio results are reported). Each simulation was run until a statistically-steady equilibrium was reached, determined by the convergence of the time-mean values of Nu, Re, δ b , and the equivalence of the four methods of estimating Nu. Since Nu ≈ Nu w ≈ ε b ≈ 1 + ε u /Ra for all simulations (not shown), verifying statistical steadiness, only Nu is discussed hereafter. All three methods of estimating the Reynolds number also produce very similar results ( fig. 2b), and the free-convective scaling (with proportionality factor ≈ 0.4) gives good agreement with the observed scaling, especially for Ra 10 6 . Figure 1 shows single-time snapshots of the 2D buoyancy field in fully developed RBC at various Rayleigh numbers. The solutions show several characteristic regimes. For Ra < Ra c , diffusion damps out any motion and the solution is entirely diffusive (not shown). As Ra increases above Ra crit the solution exhibits first steady convection (a), then transitional turbulence (b), and finally fully-developed convective turbulence (c-d). This broad phenomenology is valid in both 2D and 3D, so for the remainder of the paper we restrict to 2D. Reproducing this phenomenology serves both to demonstrate the usefulness of RBC as a model of convection, and to validate the chosen numerical method.
The scalings of Nu, Re, and δ b are shown in Fig. 2, along with snapshots of buoyancy fields from representative simulations in each phenomenological regime in figure 1. A transition from diffusive to convective behaviour is observed in both the Nusselt (figure 2a) and Reynolds (figure 2b) numbers at Ra 1700, in agreement with the prediction Ra c = 1708. A transition to turbulence follows between 10 7 Ra 10 8 , as expected given that Re ≈ 2000 for Ra ≈ 2×10 7 . This can be seen in the qualitative nature of the flow: figure 1a is steady, representative of all flows with Ra c Ra 10 7 ; while above Ra 2 × 10 7 the flow is intermittent and exhibits patterns on multiple scales, characteristic of turbulence, as seen in Figs. 1b-d.
The Nusselt number obeys a power law close to Ra 2/7 , and the Reynolds number a power law close to Ra 1/2 , from shortly after the onset of convection to the highest Rayleigh number considered. These are the expected exponents within this parameter regime (Ahlers et al.  (15); the theoretical scaling, Re ∝ Ra 1/2 , is shown as a solid blue line. In c), the solid blue line shows the theoretical scaling, δ b ∝ Ra −2/7 . 2009; Chillà and Schumacher 2012). The three different possibilities for the velocity scale in the Reynolds number calculation give similar results. A reduction in the prefactor of the power law for Nu is observed between 10 7 < Ra < 10 8 , which coincides with the onset of turbulence. A similar transition is seen in the results of Johnston and Doering (2009) for finite-difference DNS of 2D dry RBC with Pr = 1. Above Ra 10 7 they observe a power law relationship between Ra and Nu of Nu = 0.138 Ra 0.285 , which our data are in excellent agreement with.
Any two-fluid parametrization of RBC should therefore aim to capture the described scaling behaviour of Nu and Re with Ra.
3 Multi-fluid equation set and closure choices As a first step towards building a multi-fluid parametrization of convective turbulence, we motivate and present a two-fluid single-column model of Rayleigh-Bénard convection. The full viscous multi-fluid Boussinesq equation set is (Shipley et al. 2021): Here an overbar denotes a spatial filter (Germano 1992); i ∈ {0, 1, . . . , n} indexes the fluid partitions; I i is an indicator function for fluid i; σ i := I i is the fraction of fluid i contained within a characteristic filter volume; u i := I i u σ i and b i := I i b σ i are the velocity and buoyancy fields of fluid i; p i := I i P σ i − P is the difference between the conditionally-filtered pressure in fluid i and the unconditionally filtered pressure P ; S ± i , uS ± i , bS ± i are respectively sources and sinks of fluid fraction, momentum, and buoyancy in fluid i arising from the relabelling of fluid. The unconditionally-filtered pressure, P , ensures the incompressibility of the mean flow, equation (24).
Equations (21)-(25) are derived by conditionally spatially filtering the Boussinesq equations (1)-(3) in the manner set out by Thuburn et al. (2018); however, here viscous terms and sources and sinks of fluid fraction are retained from the outset. The only terms neglected here are those arising from possible non-commutation of the spatial filter with the partial derivatives. For a full derivation and discussion of the terms requiring closure, see Shipley et al. (2021).

Closures
The terms in equations (21)-(23) which require closure can be split into: • p i , the difference between the conditionally-filtered pressure in fluid i and the unconditionally filtered pressure; • P ∇σ i − P ∇I i , −ν∇ · (u ⊗ ∇I i ) T − ν∇I i · (∇u) T , and −κ∇I i · ∇b − κ∇ · b∇I i , which arise from conditionally-filtering the pressure gradient, viscous diffusion, and buoyancy diffusion terms; • −∇ · I i u ⊗ u − σ i u i ⊗ u i , and −∇ · I i ub − σ i u i b i , which are often termed "subfilter fluxes" and are akin to the Reynolds stress and subfilter buoyancy flux, respectively, in normal higher-order modelling of turbulence; • S ± i , uS ± i , bS ± i , which arise from filtering the re-labelling of fluid parcels.
We present closures that attempt to model the dominant coherent overturning structures of RBC (seen in the DNS, figure 1). For this study, differences between conditionally-filtered and unconditionally-filtered pressures are parametrized as p i = j σ j γ i ∇ · u j − γ i ∇ · u i , where γ i is a volume (or "bulk") viscosity. This has successfully been used by Weller et al. (2020), where it was argued that such a form is plausible since in the underlying Boussinesq flow, the pressure is simply a Lagrange multiplier to enforce the divergence-free condition. It is also possible to derive this form by analogy with the "bulk viscous pressure" which arises in compressible fluid dynamics, as in Batchelor (1967) -see Shipley et al. (2021) for details.
Residual terms arising from conditionally-filtering the pressure gradient and diffusion terms are closed via a mean-field approximation: These choices retain the correct sum over all fluids for the entire pressure, viscous, and diffusive terms, respectively. They also cause the fluid fractions to behave passively in the case of two fluids, and in the absence of transfers: the Eulerian derivatives for u i and b i do not depend on σ i if the two fluids have the same u i and b i .
The resolved velocities and buoyancies of the multi-fluid split are assumed to dominate the single-fluid subfilter fluxes, such that I i u ⊗ u ≈ σ i u i ⊗ u i and I i ub ≈ σ i b i u i . This is the same as assuming that the multi-fluid split captures all of the subfilter variability in the momentum and buoyancy fluxes, i.e. neglecting the residual subfilter fluxes of momentum and buoyancy, ∇ · I i u ⊗ u − σ i u i ⊗ u i and ∇ · I i ub − σ i u i b i . While this will never be exactly true, it is instructive to see how well a multi-fluid model with no extra subfilter modelling can perform when simulating a fully turbulent flow. In the single column context this requires the vertical grid to adequately resolve the boundary layers, as in the DNS.
To proceed further, we must decide what the labels I i represent. The simplest choice is to restrict to two fluids; the symmetries of the Rayleigh-Bénard problem suggest choosing one falling and the other rising: let i = 0 denote fluid with w ≤ 0, and i = 1 denote fluid with w > 0 (as in Weller et al. 2020). Then fluid 1 represents "updrafts" while fluid 0 represents "downdrafts". This choice of definitions for the two fluids, coupled with the discrete symmetry of the unfiltered equations under the simultaneous transformations z → H 2 − z, b → −b, forces D σ i dV = 1 2 . This constraint can be used as a "sanity check" for both the initial conditions and the transfer terms S ± i . The discrete symmetry of the fluids under exchange also forces γ i = γ j if γ is not a function of z.
Specializing to two fluids allows the sources of fluid fraction i to be written as S + i = σ j S ji , where σ j S ji is the rate of transfer of fluid fraction from j to i. A similar relation follows for the sinks. We choose to model the exchanges of momentum and buoyancy from fluid i to j as a characteristic value, u T ij or b T ij , times the rate of transfer of fluid fraction from i to j, σ i S ij . This aligns with the modelling approach taken in other recent works on multi-fluid modelling (Thuburn et al. 2018;Thuburn et al. 2019;Weller and McIntyre 2019;Weller et al. 2020;McIntyre et al. 2020).
Partitioning the flow based on the sign of w forces w T ij = 0. For a single-column model, it remains only to specify the form of the fluid fraction transfer rate, S ij , and the transferred buoyancy, b T ij (for a 2D or 3D model, the horizontal components of the transferred velocity would also need to be specified). For the fluid fraction transfer rate we choose which in 1D is similar to dynamical entrainment, and follows the successful implementation of the same divergence-based transfer in Weller et al. (2020). This aims to capture the largescale overturning circulation, and is exactly correct for the first normal mode of RBC with stress-free boundaries (Shipley et al. 2021). McIntyre (2020, chapter 2) also shows that this choice of transfer rate removes the problematic Kelvin-Helmholtz-like instability for a two-fluid Boussinesq system (Thuburn et al. 2019).
The transferred buoyancy must depend on the distribution of buoyancy within each fluid, and on the detailed dynamics of the relabelling. In the absence of this information, we choose a simple model: with some dimensionless constant C ≥ 0. That is, the buoyancy of fluid parcels relabelled from i to j is modelled as the mean buoyancy within the fluid i plus or minus some constant times the magnitude of the buoyancy, to crudely approximate the subfilter buoyancy variability. The signs are chosen to model the fact that the fluid transferred from the falling (0) to the rising (1) fluid is expected to be more buoyant than the average falling fluid parcel for that height, while the reverse should be true for transfers from the rising (1) to the falling (0) fluid. This is a similar formulation to that used by Thuburn et al. (2019), though in theirs the transferred value depends on both the initial and destination fluids, rather than just the initial fluid. Making these closure assumptions reduces the equation set to: with i ∈ {0, 1}, and the specific parametrization choices: The equations are given in vector form because of the desire to eventually create a 3D greyzone convection parametrization; to that end the subsequent numerical method is also threedimensional. Note, however, that in the form (28)-(34), the horizontal components of the transferred velocity still require closure.

Boundary conditions
Conditionally filtering the boundary conditions for RBC gives u i (z = 0, H) = 0, b i (z = 0, H) = ± ∆B 2 . The Neumann boundary condition for the unconditionally filtered pressure (required for the numerical solution, which solves elliptic equations for the pressures) is hydrostatic, dP dz = b(z = 0, H). Boundary conditions on the perturbation pressures are chosen to be zerogradient, dp i dz = 0. Because the σ i equation is a transport equation with no diffusion, boundary values of σ i are not in the domain of dependence of its solution. The asymptotic boundary behaviour of σ i is thus entirely dependent on the asymptotic behaviour of the transfer terms as the boundaries are approached. Boundary values of σ i are however required for the momentum and buoyancy equations, which do contain second derivatives of σ i . These boundary values should be set by extrapolated values of σ i from the interior of the domain. However, for this study we choose zero-gradient conditions for σ i for better numerical behaviour. Heuristically this means that we are imposing no creation of fluid in either partition at the boundary.

Scaling of pressure differences between fluids
In single-column form, equations (28)-(34) contain two free parameters: γ and C. C is dimensionless and should be O(1), but γ has the dimensions of (bulk) viscosity and does not have an obvious magnitude. In this section we present a scaling argument for γ with the external dimensionless control parameters Ra, Pr, thus reducing the model to the choice of two dimensionless constants which should both be O(1).
In convection, a distinction is often made between filamentary plumes and a well-mixed environment; this distinction is clearly seen in the example RBC buoyancy fields of figure 1, and is the basis of the conceptual "updraft"-"environment" partition. We assume that such a plume has a length O(H), a width δ, and the along-plume flow scales with the large-scale circulation U ∼ U B = √ ∆B H. Orienting a local Cartesian co-ordinate system such thatx points parallel to the plume andẑ points normal to it, the scaled continuity equation gives: Splitting the buoyancy equation similarly into its plume-parallel and -normal parts gives: Note the buoyancy scaling cancels here. The simplest choice of the time scale is T b = δ 2 /κ, which makes the coefficient of the final term on the RHS one, consistent with filamentary plumes being diffusion-limited in well-developed turbulent flows. Scaling the plume-parallel momentum equation with time scale T m = T b / Pr, buoyancy with ∆B and pressure with P ∼ U 2 (Bernoulli scaling), leads to: where Re = U H/ν = Pr −1/2 Ra 1/2 . The pressure gradient and buoyancy terms are assumed to drive the flow, and so Re δ 2 /H 2 = O(1) and: Hence the across-plume pressure contrast -i.e. the difference in pressure between the plume and the bulk -may be scaled as P z = P δ/H = ∆Bδ.
These results are the standard Prandtl-Blasius results with δ the boundary-layer depth, consistent with the presumption that plumes in RBC are simply detached from the boundary layers. This is a standard assumption for the kinetic boundary layer depth in scaling analysis of RBC, for example in the successful theory of Grossmann and Lohse (2000) for the Nusselt and Reynolds number scalings. The Re ∝ Ra 1/2 result is also expected for RBC in the parameter regimes under study in this paper (Ahlers et al. 2009, table 2).
We wish to parametrize the difference between the conditionally-filtered pressure in partition i, and the unconditionally-filtered pressure, as a bulk viscous stress: (34). Assuming that the multi-fluid split is dominated by a plume vs. bulk contrast, then u i scales with the velocity of the plumes, √ ∆B H, and the divergence within each fluid should then scale as ∇ · u i = (U/H)∇ ·ũ i (so long as the filter width is O(H)). Collecting the nondimensionalized expressions for the pressure and the bulk viscous stress gives: introducing the O(1), dimensionless constantγ 0 . This scaling law for γ(Ra, Pr) reduces the model for the pressure perturbation to the specification of an O(1) constant,γ 0 . Althoughγ 0 must be determined empirically, this determination need only be performed at one Rayleigh number. Since Pr = 0.707 is constant throughout our experiments, we choose to subsume the factor of Pr −1/2 ≈ 1.19 into the definition ofγ 0 from now on.

Single-fluid solver
Single-fluid reference solutions (section 2.1) were computed using the single-fluid Boussinesq finite volume code boussinesqFoam (available at www.github.com/AtmosFOAM/AtmosFOAM). This solves the single-fluid Boussinesq equation set (5)-(6) using precisely the same numerical method as detailed below for the multi-fluid equation set, but with only one fluid. This singlefluid solver gives statistically identical results to the multi-fluid solver when the latter is run with no coupling terms between the fluids, such that the σ i are simply passive tracers.

Two-fluid solver
The two-fluid Boussinesq equation set (30)-(34) is solved in advective form using the finite volume solver multiFluidBoussinesqFoam; this is part of the AtmosFOAM library of CFD codes for atmospheric fluid dynamics, based on the OpenFOAM open-source CFD library. The code is available at www.github.com/AtmosFOAM/AtmosFOAM-multiFluid. The method is similar to that detailed in section 3 of Weller et al. (2020); an overview, and choices specific to this paper, are presented below.
The spatial discretization uses Arakawa C-grid staggering in the horizontal and Lorenz staggering in the vertical. Temporal discretization is Crank-Nicolson with off-centring coefficient α = 0.55. Prognostic variables are b i and σ i at cell centres, and the volume flux φ i := u i · S f at cell faces, where S f is the outward-pointing area vector of face f . Advection of b i and σ i is total variation-diminishing (with a van Leer limiter) to preserve boundedness, while advection of φ i is linear upwind. Thus the spatial discretization is (almost) second-order accurate.
The transfer terms S ij are handled explicitly, while the momentum and buoyancy transfers are implicit and operator-split, as in (Weller and McIntyre 2019;McIntyre et al. 2020;Weller et al. 2020).
Diagnostic variables are the pressures P and p i at cell centres. Solutions for both P and p i are implicit but not simultaneous: first a Poisson equation is solved for P , which maintains a divergence-free mean velocity field (i.e. it ensures eq. (24) is satisfied), followed by a Helmholtz equation for each p i . These solutions are then iterated to convergence. The generalized Geometric-Algebraic MultiGrid (GAMG) method is used for the implicit pressure solves, with an absolute tolerance of 10 −6 .
Two outer iterations (for the whole of the above method) and two inner iterations (for the implicit pressure solves) are performed per time-step.
Apart from the transfer terms, this method is suitable for an arbitrary number of fluids, in up to 3 spatial dimensions. However, the transfer terms, and their inclusion into the algorithm, are currently specific to two fluids.   For 10 2 ≤ Ra ≤ 10 10 , single-column two-fluid simulations were run with the same vertical resolution as the reference DNS (see table 1) for various values ofγ 0 and C. The qualitative nature of the solutions is described in section 5.1, followed by an analysis of sensitivity to the choice ofγ 0 and C in section 5.2. In section 5.3 the global buoyancy and momentum transport, Nu and Re, is examined as a function of the buoyancy forcing Ra.

Two-fluid single-column model results
For all simulations, the initial state was constructed from a resting hydrostatically-balanced solution with a linear buoyancy profile and uniform σ i = 0.5 in both fluids. Small non-zero velocities equal to ±10 −3 U B were added to ensure correct labeling, and random perturbations of magnitude |δb| ≤ 0.0008 ∆B drawn from a uniform distribution were added to the initial linear profile to seed instability 5 . Simulations were run until a steady state was reached (9 − 12T e ); the steady-state profiles of buoyancy, pressure, vertical velocity, and fluid fraction, were then compared with the corresponding statistically steady-state time-mean conditionally horizontally averaged DNS profiles. Resolutions, time-step size, and total simulation run time for each simulation are given in table 2. The single column model spins up to equilibrium in a remarkably similar manner to the horizontally-averaged DNS; this is demonstrated in figure 3, which shows the Nusselt number vs. time for both DNS and single-column simulations at Ra = 10 5 and 10 8 . For these simulations, γ 0 = 1.861, and C = 0.5, 0 for Ra = 10 5 , 10 8 , respectively (see section 5.3). At each Ra, convection initiates at a similar time (≈ 2T e ) in both the single-column and DNS flows, seen in the sharp increase in Nu above the purely diffusive value of 1. This initial convective surge causes a strong peak in the Nusselt number (slightly overestimated by the single-column model), before the system gradually settles down towards equilibrium with decaying Nusselt number under-and overshoots. The under-and overshoots appear stochastic for the DNS, whereas they are periodic for the single-column model; that the single-column model appears less chaotic than the DNS is unsurprising.
The same steady state was reached when initializing from other initial conditions (e.g. initializing from the DNS reference profiles), provided the identities of the fluids were initialized correctly and the initial column-integrated fraction of fluid in each fluid was equal to 0.5. This suggests that the steady state is robust. Similar qualitative spin-up behaviour is also observed with different values ofγ 0 and C. Thus, for the remainder of the paper we consider only the steady state, and not the spin-up.
We begin our study of the two-fluid single-column model steady-state by looking at the qualitative behaviour of the equilibrium profiles in different Rayleigh number regimes. We then investigate the sensitivity of those profiles to the two closure constants, C andγ 0 . Finally we examine the scaling of the global parameters Nu and Re with Ra produced by the model.

Phenomenology
For each of the characteristic Rayleigh numbers Ra = 10 5 , 10 8 , 10 10 (as in figure 1), we present and discuss an example two-fluid single column simulation. Rather than use the fixed value used above for discussion of the spin-up, the values ofγ 0 and C in the example simulations were chosen to have the best qualitative fit to the conditionally horizontally averaged DNS for all profiles. The discussion for each of these examples qualitatively applies to all simulations within the characteristic Rayleigh number regime.
Laminar (Ra = 10 5 ) At Ra = 10 5 , the DNS exhibits laminar convective rolls (see Fig. 1a). This solution is qualitatively characteristic of the flow for all laminar Ra, Ra c < Ra 10 7 . Steady state results of a two-fluid single-column model governed by equations (28)-(34) withγ 0 ≈ 0.75, C = 0.5 are shown in Fig. 4. The mean buoyancy (a) and pressure (b) profiles match closely between the DNS and the single column model; in particular the model correctly predicts a well-mixed buoyancy in the fluid interior, with a sharp buoyancy gradient close to the top and bottom boundaries. The shape of the pressure profile is also correct, though the maxima are slightly too high close to the boundaries.
Good agreement is also seen between the DNS and two-fluid single column model for the individual fluid buoyancy profiles: the overall shape is correct, though the profiles are too far apart in the middle of the domain, leading to surplus buoyancy transport for a given velocity profile. Experiments varying C (see section 5.2.2) demonstrated C > 0 was required to reproduce a buoyancy overshoot at the top (bottom) of the rising (falling) fluid. By overshoot, we mean the part of the buoyancy profile at the interface between the bulk and the buoyancy boundary layer where db i dz changes sign. These overshoots can be seen in the 2D buoyancy field of the DNS flow of figure 1a and are a general feature of O(1) Prandtl number laminar RBC. (For Pr > 1, the overshoots become so strong that they begin to be seen even in the mean buoyancy profile; such profiles can be seen in e.g. Fig. 4b of Schmalzl et al. 2004.) The value C = 0.6 gives the best shape for b i (z) for Ra = 10 5 , but C ≈ 0.5 works for all laminar Ra.
The individual fluid velocity profiles are roughly the correct shape; the slight asymmetry in the location of the maxima in each fluid in the DNS is due to the gradient of the volume fraction profile in the DNS (i.e. forcing the correct gradient of σ i reproduces the asymmetry in the vertical velocity profiles).
The pressure profiles with each fluid are captured by the scheme, suggesting that to leading order p i ∝ −γ∇ · u i is an appropriate model of the pressure differences. The model is particularly good close to the boundaries, but the fluids are better mixed in the interior of the domain in the DNS, causing the pressure differences there to be smaller than predicted by the single column model. This could possibly be remedied by using a z-dependent γ parametrization, which would fit well with the discussion of LES in section 2.2.1.
The two-fluid model keeps area fractions, σ i (z), close to 0.5. This is expected as the divergence-based transfer is known to keep σ i (z) roughly constant . In contrast, the area fractions diagnosed from the DNS diverge from 0.5 either side of the centre (where symmetry demands equal fractions), reaching a maximum close to the boundaries approaching 0.3 and 0.7. Transition to turbulence (Ra = 10 8 ) Between 10 7 < Ra 5 × 10 8 , the DNS solutions transition from laminar flow to fully developed turbulence. The buoyancy field of figure 1b is characteristic of this transitional regime. Besides the solutions becoming intermittent and transient rather than (quasi-)periodic, the plume separation from the boundary layer fundamentally changes: above Ra ≈ 10 7 , regions of recirculation develop at the base of the plumes.
Results of a two-fluid single-column model withγ 0 ≈ 0.47, C = 0 are compared with those from the horizontally-averaged DNS in Fig. 5. As with the Ra = 10 5 results, the values ofγ 0 and C were chosen to give the best qualitative agreement for all profiles. Better prediction of the pressure differences between the fluids near the boundaries is achieved by increasingγ 0 by a factor of ≈ 2; however this degrades the agreement of the mean pressure profile with the DNS profile. This again suggests that γ should be a function of z, either directly or through dependence on other properties of the flow, for instance the TKE.
Comparisons with the DNS reference profiles are mostly the same as for the laminar case, except that the additional mixing caused by the recirculation regions at the base of the plumes modifies the profiles in the near-boundary regions. This has the most obvious effect on the buoyancy profiles within each fluid, which no longer overshoot, and on the volume fraction profile, which is no longer monotonic. The lack of overshoots is reproduced by transferring the mean buoyancy, C = 0, a suitable model for well-mixed turbulent flow. The detailed differences to the profiles caused by these recirculation regions are however not reproduced by this simple parametrization: better representation of the mass exchanges S ij is required. The recirculation is counter to the large-scale circulation, and hence is not captured either by our arguments for the scaling of γ, or by the divergence-based mass transfer.
Fully developed turbulence (Ra = 10 10 ) Above Ra 5 × 10 8 , the DNS flow is fully turbulent, exhibiting structures on many scales from the domain depth down to the exceptionally thin boundary layers, shown in Figs. 1cd for Ra = 10 10 . The recirculations at plume base first exhibited in the transitional regime divide into multiple small plumes which organize into a larger-scale circulation. The bulk of the domain is statistically well-mixed.
Results from a two-fluid single-column model withγ 0 ≈ 0.44, C = 0 are shown in Fig. 6. Qualitative agreement with the buoyancy and vertical velocity profiles is still good, but the mean pressure profile predicted by the model now has too little curvature in the centre of the domain, and does not get the gradient correct close to the boundaries. Again, the complex mixing of the turbulent flow has strong effects on the volume fraction profile, causing the volume fraction of rising (falling) fluid to be less than 0.5 close to the lower (upper) boundary.
These larger discrepancies between the DNS and the two-fluid model model are possibly because the w = 0 interface is now very complex. Figure 7 shows the w = 0 interface superimposed on the DNS buoyancy fields at Ra = 10 8 and Ra = 10 10 . Although the dominant rising/falling two-fluid split is still into columns of falling and rising air with an approximately vertical interface even in the higher Ra case, the simple split is increasingly complicated by the complex vortical motions in the bulk of the fluid, and especially close to the base of the plumes. The intricate dynamics of these interfaces are not accounted for by our single-column model.
While there are quantitative discrepancies, for all three Rayleigh numbers the overall the agreement between horizontally-averaged DNS and the two-fluid single-column model is good. Approximately the correct profiles are captured even in the highly turbulent regime of Ra = 10 10 . The model performs remarkably well given it has no representation of sub-filter variability beyond the two-fluid split, showing that the model captures the essential coherent overturning structures of Rayleigh-Bénard convection in all three characteristic regimes.

Sensitivity toγ 0 and C
In this section, the sensitivity of the model to the dimensionless closure parametersγ 0 and C is investigated. The effects of changingγ 0 and C are similar at all Rayleigh numbers, so for brevity only Ra = 10 5 is presented. Figure 8 shows the effect on the two-fluid single-column steady-state of varyingγ 0 from 10 −1 γ 0 10 1 , along with examples in the asymptotically-large and -smallγ 0 regimes. The experiments were performed with C = 0.5 at fixed Ra = 10 5 , but the results are similar for all Ra.

Sensitivity toγ 0
The best qualitative match between the single-column and DNS profiles is found when γ 0 ≈ 0.75, as discussed earlier, while the correct heat flux is predicted atγ 0 ≈ 1.861. These values are both O(1), as expected. Agreement with the reference profiles degrades sharply aŝ γ 0 moves away from this range.
Increasingγ 0 increases the buoyancy difference between the fluids, and damps the vertical velocities -which makes sense since in 1D this parametrization of p i is similar to diffusion of the vertical velocity within a fluid, even though the sum correction means no extra viscous term is added to the mean momentum budget. This effect is already clear atγ 0 = 2, where the vertical velocities are only ≈ 2/3 of those in the DNS, and the pressure profile is much shallower, though still with the correct number of turning points. Byγ 0 = 10, the pressure profile loses the minimum in the centre of the domain, and the vertical velocities are almost zero. At asymptotically largeγ 0 , the system becomes subcritical and the solution is purely diffusive. Decreasingγ 0 rapidly increases the pressure gradient, and deepens the minimum of the mean pressure in the centre of the domain. This drastically increases the vertical velocitiesbyγ 0 = 10 −1 , the maximum vertical velocities are over three times those of the DNS, and over twice those of the simulations withγ 0 = 0.75 discussed in detail earlier. Decreasingγ 0 further only slightly changes these results, as seen for the asymptotically-small case ofγ 0 = ×10 −5 . Figure 9 shows the steady-state effect of varying C from 0 (mean buoyancy is transferred: b T ij = b i ) to 1 (zero buoyancy is transferred over most of the domain: b T ij = 0 wherever b i = |b i |). Transfers with C > 1 amount to transferring buoyancies with magnitude greater than ∆B close to the boundaries, which causes the solution to become unstable at C ≈ 1.3. The main effect of increasing C is to generate the aforementioned overshoots in the withinfluid buoyancy profiles; this also steepens the pressure gradient, deepens the central pressure, and increases the magnitude of the vertical velocities in each fluid. These effects are small compared to the order-of-magnitude effects associated with varyingγ 0 : for example, the maximum velocity increases monotonically from 0.3 to 0.45 as C increases from 0 to 1. These effects are qualitatively similar at all Ra, but for Ra 10 7 , the individual fluid buoyancy profiles no longer exhibit overshoots, so C = 0 provides a better fit with the DNS buoyancy profiles.

Scaling of Nusselt number with Rayleigh number
To investigate the performance of the two-fluid single column model more systematically, the scaling of the Nusselt number for single-column models across the Rayleigh number range 10 2 ≤ Ra ≤ 10 10 is compared with the DNS results. The scaling γ/ν ∝ Ra 1/4 (section 3.2) is evaluated, along with two choices of the transferred buoyancy, C = 0 and C = 0.5. For each transferred buoyancy, the dimensionless proportionality factorγ 0 was fixed by finding the value which gave the correct Nusselt number at Ra = 10 5 . Fixing this constant at different Rayleigh numbers changes the prefactor of the Nu(Ra) scaling, but does not change the scaling itself. Figure 10a shows Nu against Ra for the different values of C and scalings for γ. The DNS results are shown for comparison, along with results from the single column model run with both tunable parameters set to zero, C =γ 0 = 0. All models withγ 0 > 0 perform significantly better than the model withγ 0 = 0, which becomes supercritical for Ra < 10 3 and follows a Nu(Ra) scaling with exponent everywhere > 0.33.
Models with γ/ν ∝ Ra 1/4 show exceptional agreement with the DNS heat fluxes for Ra ≥ 10 4 , giving Nu ∼ Ra 2/7 with both C = 0 (green curve) and C = 0.5 (purple curve). This shows that the Nusselt number scaling exponent depends on γ/ν but not on C; this makes sense since C is a crude parametrization for how the flow produces a given heat flux, and should not affect the scaling of the heat flux itself. Below Ra = 10 4 , the models with different values of C produce slightly different behaviour: the C = 0 solutions become supercritical below Ra = 10 3 , inconsistent with the known Ra c ≈ 1708. While the C = 0 simulations are still subcritical at Ra = 10 3 , the heat flux at Ra = 2 × 10 3 is roughly 30% too high. These discrepancies suggest that the scaling used for γ/ν is not quite correct in the low Ra regime; unsurprising since the scaling argument assumed Re 1. For the intended application to highly turbulent atmospheric convection, however, this does not present a severe problem.
The single-column model does not naturally capture the drop in the prefactor of the Nusselt number scaling which occurs as the flow transitions to turbulence around Ra ≈ 10 7 . The drop in the Nusselt number scaling prefactor may not be a robust feature of the convective flow, so it is far more important to get the scaling exponent correct. Such drops in the scaling prefactor are found in other RBC experiments (see Johnston and Doering 2009;Roche et al. 2004, for a 2D numerical and a 3D experimental example, respectively), but appear to be dependent directly on the nature of the flow, rather than global in nature like the scaling exponent. However, this drop can be accurately reproduced by using C = 0.5 for Ra ≤ 10 7 and C = 0 for Ra > 10 7 , retaining the value ofγ 0 ≈ 1.861. With this parametrization, the Nusselt number is correctly predicted to within 5% across six orders of magnitude of buoyancy forcing, 10 4 ≤ Ra ≤ 10 10 , and approximately the correct transitional behaviour is found for Ra < 10 4 . This could be diagnostically incorporated into the parametrization by, for instance, reducing C to 0 whenever the vertical velocity maximum gives a turbulent Re 2 × 10 3 .
The Reynolds number in the single-column simulations was estimated from the maximum magnitude of the vertical velocity; this should scale with the large-scale circulation, so makes sense for a bulk Reynolds number. The scaling behaviour of the Reynolds number is also wellcaptured (figure 10b), in particular giving the same scaling exponent as the DNS. Notably, the change in C required to capture the correct behaviour of Nu does not cause a corresponding kink in the Reynolds number scaling. This suggests that C really is just a crude measure of the flow state. Future work would hope to capture these flow states dynamically through representing the sub-filter scale variability of the variables within each fluid.
6 Summary, conclusions, and future work In this paper we have shown that the simple two-fluid single column model (28)-(34) can qualitatively reproduce horizontal-mean DNS buoyancy, vertical velocity, and pressure profiles in all three characteristic regimes of Rayleigh-Bénard convection. A scaling argument for the pressure differences between the fluids allows the model to predict the correct power-law scaling (a) Nusselt number vs. Rayleigh number for two fluid single column models with various values ofγ 0 and C.
The dashed blue curve shows the reference DNS results (as 2a), while the dashed brown curve shows the results of running the single column model with C = 0 (mass exchanges transfer the mean buoyancy) and γ = 0 (no pressure differences between the fluids). The green, purple, and red curves show the results for γ ∼ Ra 1/4 , with different values of γ 0 ; all give scalings of Nu ≈ Ra 2/7 . Single-column Nusselt numbers are calculated from the buoyancy gradient at the boundaries, and checked against the column-integrated buoyancy flux.
(b) Reynolds number vs. Rayleigh number for 10 2 ≤ Ra ≤ 10 10 . The dashed blue curve shows the reference DNS results (as 2b), while the solid orange curve shows results from a two fluid single column model obeying equations (28)-(34), with γ/ν = 1.861 Ra 1/4 and C = 0.5 for Ra ≤ 10 7 , C = 0 for Ra > 10 7 ; these constants give the best fit for Nu as a function of Ra (figure 10a). Both curves exhibit scalings of Re ≈ Ra 1/2 for Ra 10 4 . Single-column Reynolds numbers are calculated using the maxima of the individual fluid vertical velocity profiles for the velocity scale. Figure 10 of Nu ∼ Ra 2/7 , and after measuring a dimensionless constant at one Rayleigh number the magnitude of Nu can be predicted to within 5% over 6 orders of magnitude of Ra. The model also captures approximately the correct spin-up behaviour, and approximately the correct critical Rayleigh number. The closure set is minimal, requiring only two constants to be set; and not finely-tuned, as both closure constants may be varied significantly from their central values without destroying the solution.
Although we use a similar equation set and identical fluid definitions to Weller et al. (2020), this is the first such study to model a fully turbulent regime with these fluid definitions. It is also the first multi-fluid convection study to considerably vary the applied forcing, testing the robustness of the parametrization.
This demonstrates the essential validity of the multi-fluid concept: the model directly captures the dominant overturning circulation of convection, present even in the fully turbulent regime, by allowing for a circulation even in a single column. It is important to note that this performance is achieved without even a minimal treatment of fluxes due to variability within each fluid (i.e. conventional 'turbulent' or 'subfilter' fluxes) apart from the fixed viscosity and Prandtl number of the fluid.
With the current model the mean buoyancy profile (and therefore the Nusselt number), the vertical velocity maxima in each fluid (and therefore the implied Reynolds number), and the pressure profile, cannot all simultaneously have the correct magnitude. It is unclear whether this is due to neglected subfilter variability (in the form of exchanged buoyancy or neglected subfilter stresses, for example), or due to inadequate representation of the fluid fraction transfers. A more accurate and flexible representation of these transfers is essential to progressing beyond single-column modelling.
Future work will test the two-fluid model of this paper in the grey zone of RBC, investigating how the closures scale with resolution, and noting what flow features are missed by the simple closures in a higher dimensional setting. Improvements could arise from a partition which better selects the coherent structures, and from representation of within-fluid variability by consideration of higher moments of the flow. In particular, DNS data may be used to diagnose S ij , b T ij , u T ij for various filter scales and fluid definitions. Possible closures could be informed by direct analysis of the interactions between coherent structures, boundary layers, and homogeneous, isotropic bulk (Togni et al. 2015;Berghout et al. 2021).
All of the above will develop fundamental understanding of the multi-fluid equations for convection. A thorough understanding of the dry convective grey zone, and of possible multifluid approaches to its parametrization, will help sharpen the questions for the much thornier problem of moist convection.