Multifluids for Representing Subgrid‐Scale Convection

Traditional parameterizations of convection are a large source of error in weather and climate prediction models, and the assumptions behind them become worse as resolution increases. Multifluid modeling is a promising new method of representing subgrid‐scale and near‐grid‐scale convection allowing for net mass transport by convection and nonequilibrium dynamics. The air is partitioned into two or more fluids, which may represent, for example, updrafts and the nonupdraft environment. Each fluid has its own velocity, temperature, and constituents with separate equations of motion. This paper presents two‐fluid Boussinesq equations for representing subgrid‐scale dry convection with sinking and w = 0 air in Fluid 0 and rising air in Fluid 1. Two vertical slice test cases are developed to tune parameters and to evaluate the two‐fluid equations: a buoyant rising bubble and radiative convective equilibrium. These are first simulated at high resolution with a single‐fluid model and conditionally averaged based on the sign of the vertical velocity. The test cases are next simulated with the two‐fluid model in one column. A model for entrainment and detrainment based on divergence leads to excellent representation of the convective area fraction. Previous multifluid modeling of convection has used the same pressure for both fluids. This is shown to be a bad approximation, and a model for the pressure difference between the fluids based on divergence is presented.


Introduction
Slow progress has been made improving the parameterization of subgrid-scale convection since Arakawa and Schubert (1974), Tiedtke (1989), Gregory and Rowntree (1990), and Kain and Fritsch (1990). Convection is still one of the weakest aspects of large-scale models of the atmosphere (Holloway et al., 2014;Sherwood et al., 2013;Solomon et al., 2007). Improvements have been made by relaxing the quasi-equilibrium assumption (Gerard & Geleyn, 2005;Pan & Randall, 1998;Park, 2014), including stochasticity (Plant & Craig, 2008) and allowing net mass transport by convection (Kuell & Bott, 2008;Malardel & Bechtold, 2019). These developments have enabled limited simulation in the gray zone, where convection is partially resolved. However, there are still systematic biases implying that improvements should still be possible. Lappen and Randall (2001) introduced the "assumed-distribution higher-order closure" (ADHOC) parameterization of the boundary layer and convection, which solves prognostic equations for higher-order moments of the subgrid-scale flow, such as the triple correlation of the vertical velocity anomalies, and diagnoses updraft and downdraft means from the higher-order moments. Thus, all necessary statistics can be calculated from the grid box-mean values that are predicted by the dynamical core. Net mass transport by convection is included because it is part of the grid box mean. ADHOC was extended by Golaz et al. 10.1029/2019MS001966 (2002 to use bi-Gaussian probability distribution functions rather than top-hat distributions to represent subgrid-scale variability. This formed the basis of the "Cloud Layers Unified By Binormals" (CLUBB) parameterization (Larson et al., 2012), which was first used at gray zone resolutions. Thuburn et al. (2018) proposed multifluid modeling of convection in order to treat nonequilibrium and net mass transport consistently. The multifluids approach has a potential advantage over ADHOC and CLUBB in solving prognostic equations for first-order moments directly (mean properties of, e.g., updraft and environment) rather than solving prognostic equations for higher-order moments of the whole fluid in order to diagnose means of updrafts and the environment. Solving directly for first-order moments in each fluid could be less sensitive to closure assumptions. The multifluids approach has an advantage over prognostic plume models such as Gerard and Geleyn (2005) as multifluids enables the mass transported by convection to feed back onto the continuity equation. Yano (2014) derived similar equations to Thuburn et al. (2018) and showed how this is a generalization of the mass-flux formulation. The Tan et al. (2018) "extended eddy diffusivity mass flux" or extended EDMF scheme uses similar equations, but Tan et al. (2018) only solve the equations in a single column without mass transported by convection feeding back onto the continuity equation.
The disadvantage of the multifluids approach is that it requires a dedicated model rather than being a stand-alone parameterization for a dynamical core. Thuburn et al. (2019) and Weller and McIntyre (2019) presented solutions of two-fluid equations using the same pressure for both fluids. Thuburn et al. (2019) showed that these equations are in fact unstable and can be stabilized by diffusion of vertical velocity whereas Weller and McIntyre (2019) stabilized the equations with diffusion and drag between fluids. Neither of these stabilization techniques are acceptable; sufficient diffusion to stabilize the multifluid equations may smooth out real features of the flow, and the diffusion between fluids and drag used by Weller and McIntyre (2019) meant that both fluids tended to move together rather than through each other. In this paper, it is shown that significant and sustained pressure differences between the fluids exist (section 4) and that the use of separate pressures stabilizes the equations. This paper presents the multifluid Boussinesq equations for representing subgrid-scale convection including a parameterization of the pressure difference between fluids (section 2.2) and a new model of entrainment and detrainment (section 2.3). The numerical method for solving the equations is described in section 3. Results of two vertical slice test cases are presented in section 4. These are first simulated at high resolution in a single fluid model, and then statistics of the flow are reproduced solving the two-fluid equations in one column.

Derivation
Many of the challenges of representing convection with the multifluid Navier-Stokes equations (e.g., those described in Weller & McIntyre, 2019) carry over to the Boussinesq equations because convection is buoyancy dominated and close to incompressible. The resolved test cases presented in section 4 give very similar results for Navier-Stokes and Boussinesq equations. Therefore, the Boussinesq equations without background stratification are chosen for simplicity: where variables are defined in Table 1. Note that the Boussinesq pressure potential, P, is actually pressure divided by a constant reference density and is referred to as just the pressure when considering the Boussinesq equations. The pressure variations are determined up to a constant in order to ensure that the velocity is divergence free.
Velocity of the single-fluid equations and velocity averaged over all fluids of the multifluid equations Rate of volume fraction transfer from fluid i to j P (m 2 s −2 ) Boussinesq pressure potential of the single-fluid equations = p ′ / r and pressure potential averaged over all fluids of the multifluid equations Buoyancy of single-fluid equations defined as -g ′ ∕ r , where ′ is departures in density from a horizontally uniform reference, r and buoyancy averaged over all fluids of the multifluid equations b i (m s −2 ) B u o y a n c yo ff l u i di The single-fluid Boussinesq equations will be conditionally averaged following Thuburn et al. (2018) to give the multifluid Boussinesq equations. Lagrangian labels, I i , pick out different fluid types or regions. I i = 1 in fluid i and 0 elsewhere. I i satisfies Lagrangian conservation without transfers: It is more straightforward to introduce transfers after volume averaging because the transfer terms for (4) involve delta functions (this will be addressed in future publications). By applying a volume average (denoted by the overbar), the fluid i volume fractions, velocities, and buoyancies are defined by which implies ∑ i i = 1.
As the system is governed by the Boussinesq equations, it is assumed that variations in density only affect the definition of the buoyancy terms (Table 1) so u i and b i are defined without density weighting.
To derive a transport equation for i , (4) is combined with I i times (3) and volume averaged. It is necessary to assume that partial derivatives commute with the averaging (as described by Thuburn et al., 2018).The equation for i is then 10.1029/2019MS001966 where i S ij is the transfer rate from fluid i to j, which must satisfy These transfer terms account for the fact that the labels I i are not truly Lagrangian conserved (i.e., fluid can transform from fluid i to fluid j).
Next the conditional averaging of terms in the momentum (1) and buoyancy (2) equations is considered. The filter scale is assumed to be equal to the grid scale, but the derivation does not rely on this. Conditionally averaging the advection terms leads to subfilter-scale fluxes, F SF , since the velocity and buoyancy are not uniform within each fluid at the filter scale: The subfilter-scale fluxes are represented by constant viscosity diffusion terms for each fluid. Molecular diffusion is assumed to be of much smaller magnitude than this constant-viscosity turbulent diffusive mixing, so molecular diffusion is neglected hereafter. The diffusion terms are written so that, in the absence of transfers between fluids and if both fluids are initialized to be identical, the fluids remain identical and i acts as a passive tracer: The Boussinesq pressure potentials within each fluid are defined as The pressure gradient can then be decomposed as The terms in square brackets are the difference between the pressure force pushing outward from the fluid calculated using the filtered variables and the fully resolved version of the same force, both of which are related to drag. They will be represented by a drag, D ij . However, this drag is likely to be smaller than the drag that is used, for example, in multiphase modeling (e.g., Guelfi et al., 2007) as it is representing the difference between two representations of drag, rather than the whole drag. Sensitivity to D ij will be investigated in section 4.
Combining all terms gives the conditionally averaged momentum and buoyancy equations: These equations include mass transfers between fluids with u T i being the velocity of the fluid transferred from i to j and b T i being the buoyancy of the fluid transferred from i to j. Equations 17 and 18 are solved in advective form in which the mass transfer terms are more complicated: where D i ∕Dt = ∕ t + u i · ∇ is the total derivative with respect to fluid i. The instabilities reported by Weller and McIntyre (2019) showed growing divergence in a fluid with vanishing volume fraction. For single-fluid equations, divergence would lead to pressure anomalies, which would act as a negative feedback on the divergence. In the multifluid case with just one pressure, fluids with a sufficiently small volume fraction cannot force pressure anomalies big enough to stabilize the divergence. Use of a separate pressure for each fluid is therefore proposed to ensure stability.

Pressure of Each Fluid
It is assumed that the Boussinesq pressure potential in each fluid is given by a perturbation, p i , from the total, P: For stability and to respect ∑ i ∇ · i u i = 0, we propose four constraints on the pressure perturbations, p i : 1.
∑ i i p i = 0. 2. p i = 0 when u i = u so that pressure perturbations do not force individual fluid anomalies away from the mean. 3. p i ↛ 0 when i → 0 so that pressure perturbations control divergence in vanishing fluids. 4. p i ↛ ∞ as j → 0 for any combination of i, j.
In order to derive a suitable model for p i , we will first consider the Baer and Nunziato (1986) model of multiphase flow. Baer and Nunziato (1986) proposed a multiphase explosion model that has been widely used in mathematics and engineering for simulating, for example, fluidized beds of granular material (e.g., Embid et al., 1992) and two compressible fluids (e.g., Saurel & Abgrall, 1999). This model has a separate pressure for each phase, P i , and provides a transport equation for the phase fraction (equivalent of i ), which closes the system enabling diagnosis of P i . This model is not used here but is examined in order to derive a parameterization of pressure differences between fluids suitable for atmospheric convection. Saurel and Abgrall (1999) wrote the Baer and Nunziato (1986) model using transport equations for both the volume fraction, i , and for the mass fraction, i i , without transfers between fluids, which in our notation are

The Baer and Nunziato (1986) Model
where V is the velocity of the interface between fluids; is the compaction viscosity, which controls how quickly the pressure of each fluid relaxes to the mean pressure; and index j refers to the other fluid in a two-fluid system.
Different authors make different assumptions about the interface velocity. Assuming uniform density, the combination of Equations 23 and 24 (to eliminate i / t) shows that the pressure difference between fluids is related to divergence.

Proposed Model for Pressure Differences between Fluids
Thuburn et al. (2019) stabilized the multifluid equations by adding a diffusion term w 2 / z 2 to the RHS of the w equation. If a fluid perturbation pressure is set to p i = − w/ z, then this is equivalent to the diffusion term in one dimension. However, adding an artificial diffusion term will harm three-dimensional simulations of resolved convection in the free troposphere, whereas adding a fluid perturbation pressure will not. Inspired by both Baer and Nunziato (1986) and Thuburn et al. (2019) and obeying the constraints given above, the pressure in each fluid is chosen to be where is a volume viscosity with units m 2 s −1 . The "stabilization" term reduces to that of Thuburn et al. (2019) in one dimension and the "correction" term ensures that ∑ i i p i = 0. It should be noted that this parameterization destroys exact energy conservation. In future work solving the compressible equations, the energy lost due to this stabilization will be reinserted as turbulent kinetic energy or heat.
Scale analysis in section 2.6 will inform the choice of , and sensitivity to will be evaluated in section 4.

Entrainment and Detrainment
It is useful to compare the transfer terms, S ij , with models of entrainment, , and detrainment, used in convection parameterizations. If Fluid 0 is the environment and Fluid 1 is updrafts, then the relationship is Results using two formulations of entrainment are presented in section 4. The first is similar to the dynamic entrainment described by Houghton and Cramer (1951), Asai and Kasahara (1967), and de Rooy et al. (2013): This formulation prevents one fluid from accumulating when a fluid decelerates (at the top of a rising plume or at the bottom of descending air). A common form of lateral entrainment (e.g., de Rooy et al., 2013) is also tested in section 4: where r c is the cloud or plume radius.

Drag in the Momentum Equation
Pressure differences between the fluids can lead to form drag, which is parameterized following Simpson and Wiggert (1969), Romps and Charn (2015), and Weller and McIntyre (2019) as where C D is a drag coefficient and r c is a cloud radius. Sensitivity to C D /r c is explored in section 4.

The Buoyancy and the Momentum of the Mass That Is Transferred
Each fluid may not be well mixed, so the fluid transferred may not have the mean properties of the fluid it is leaving (as was assumed by Weller & McIntyre, 2019). A realistic estimate of the properties of the fluid transferred can only be based on knowledge of the subfilter-scale variability (i.e., the variability within each fluid). Modeling of subfilter-scale variability is beyond the scope of this paper. However, since conditional averaging in this study is based on the vertical velocity, it is possible to say something about the vertical velocity transfer. The subset of Fluid 0 that is about to rise (i.e., within the next time step) should be transferred from Fluids 0 to 1, and vice versa for the fluid about to fall in Fluid 1. The fluid transferred will therefore

10.1029/2019MS001966
have a vertical velocity of w = 0. For three-dimensional multifluid simulations one could assume that the horizontal velocity transferred is equal to the mean horizontal velocity of the fluid transferred from However, only one-dimensional multifluid simulations are presented in this paper.
The assumption that w T i = 0 does not imply anything about b T i . In the absence of knowledge of the subfilter-scale variability, the buoyancy transferred follows Thuburn et al. (2019) Results are presented using b = 1 2 and b = 1. b = 1 is described as an upstream approximation by Yano (2014) and is widely used for entrainment (e.g., Asai & Kasahara, 1967;Arakawa & Schubert, 1974).
Note that b ≠ 1 implies that the buoyancy of the fluid transferred depends on the properties of the receiving fluid, which is not logical and can lead to unbounded values of b i in the fluid that is losing mass (as demonstrated in section 4).

The model b T
i = 0 is also tested. The rationale is that if fluid zero contains negatively buoyant air and fluid one contains positively buoyant air, then the air that is transferred has zero buoyancy.

Scale Analysis
Scale analysis can guide our choice of for the parameterization of the pressure difference between. The scaling presented can be derived using fewer assumptions than those below, but these assumptions simplify the analysis: 1. The flow is buoyancy dominated. 2. The flow is close to inviscid. 3. The flow is slowly varying. 4. Fluid i properties are anomalies from a neutrally stable resting mean state implying that the average over all fluids of b i , u i , and P i are zero.

Only the vertical direction is considered where
8. There are two fluids. 9. Without loss of generality, the transfers S i = − w i z > 0 and S ji = 0 are used. 10. The drag term, D ij , is zero.
These assumptions can be used to simplify the multifluid w equation, (19): This is nondimensionalized using a length scale L, a buoyancy scale B, and a time scale T to get the nondimensional variables:w The nondimensional version of (32) is then 10.1029/2019MS001966 The flow is buoyancy dominated, so the buoyancy term should be O(1), which implies the scaling B = L∕T 2 ⇒ T = (L∕B) 1 2 . Therefore, the nondimensional momentum equation becomes The stabilization term must be large enough to smooth out oscillations inw due to advection and buoyancy but not too large to remove all variability inw. Therefore, the stabilization term should be O(1), which implies that This will be used to guide the choice of in section 4.

Numerical Method
The spatial discretization follows exactly Weller and McIntyre (2019), solving advective form equations using monotonic, finite volume advection and C-grid, Lorenz staggering. The time stepping is Crank-Nicolson with no off-centering and with deferred correction of explicitly solved variables. The method for calculating a separate pressure for each fluid is new because Weller and McIntyre (2019) solved the Euler equations with a single pressure. Two outer iterations per time step are used to update explicitly solved variables. For the first outer iteration, predicted values for time n + 1 are set to those from time level n and are given the superscript . For the second iteration, values at are set to those from the end of the first iteration. The implicit numerical method for applying the transfers between fluids is specific for two fluids because it involves the inversion of 2 × 2 matrices.
The prognostic variables are b i in cell centers, i in cell centers, and the velocity flux at cell faces, u i = u i · S, which is the velocity at cell faces dotted with the area vector for each face (normal to the face). The pressures, P and p i , are diagnostic and are located at cell centers.
The first calculation of each outer iteration is for the transfer terms, S ij using (27) or (28). Next, (9) is solved for each i , operator split: first applying advection, then correcting ∑ i i = 1, and finally applying mass transfers: where Δt is the time step and the superscripts ′ and ′′ denote intermediate values. n i is interpolated from cell centers onto cell faces using monotonic van Leer advection (operator ( n i ) vL ) as in Weller and McIntyre (2019) and values exclusively at time level n for best accuracy and guaranteed monotonicity. For calculating the divergence, the normal component of u i is needed at cell faces, which are prognostic variables of the C grid. The S ij is limited to ensure 0 < i < 1 ∀i.
Next, (20) is solved for the buoyancy in each fluid, first transporting buoyancy and then applying the transfer terms. The advection term, u i · ∇b i , is split into ∇ · (u i b i ) − b i ∇ ·u i for better conservation as calculated by the finite volume method. The transfer terms can be very large due to the presence of i S i in the transfer term for b i , so they are treated operator split and implicitly:

Journal of Advances in Modeling Earth Systems
The use of ( 0 / 1 ) ′′ from before the mass transfers is necessary to ensure bounded and conservative transfers, as described by McIntyre et al. (2020). Simultaneous solutions of (41) and (42) give Next, a Poisson equation is solved to calculate the total pressure, P, to ensure that the total velocity is divergence free. This is the standard approach for ensuring divergence constraints (e.g., Smolarkiewicz et al., 2014). The pressure equation also provides updates for the velocity flux, u i . First, an intermediate velocity flux is calculated with partial updates and updates from a previous value of p i but without updates from ∇P n + 1 : where operator () f means linear interpolation from cell centers onto cell faces and ∇ f P is a compact discretization of (∇ P) f · S: that is, the normal component of the pressure gradient calculated from just the values of pressure either side of the face. Once P n + 1 is calculated, the velocity flux (without momentum transfers) is updated from the back-substitution: The Poisson equation for P n + 1 is found by multiplying (47) by i , summing over all fluids and taking the divergence. Knowing that ∇ · ∑ i i u i = 0, which is solved to find P n + 1 .
Next, the momentum transfers due to mass transfer and due to drag are calculated implicitly, similar to the implicit transfers of buoyancy, calculating u n+1 i from u ′′ i . These do not influence the total P or total divergence because momentum is transferred conservatively. However, they do influence p i , which is why they are calculated before p i . The next step is to solve Helmholtz equations implicitly for each p n+1 The stages calculating P, p i , and u i are repeated twice per time step iteration. The Poisson and Helmholtz equations are solved to a tolerance of 10 −6 using the OpenFOAM conjugate gradient solver with incomplete Cholesky preconditioning.

Discrete Conditional Averaging
The multifluid equations will be evaluated in section 4 by comparison with high-resolution single fluid solutions. This has to be done carefully as even these reference solutions exist on a discrete grid. A single-fluid solution is conditionally averaged by If an entire grid box is assigned to be either in Fluid 0 or 1 based on w at the cell center, then can only take a small number of discrete values depending on the resolution of the conditional averaging. For example, if conditional averaging goes from n columns to one column with vertical resolution staying the same, then can only take values k n where k is an integer between zero and n. This means that is either uniform (for small n) or has jumps in the vertical (for larger n). An alternative is therefore proposed in order to approximate based on w and ∇w at cell centers. This is based on the intersection between a grid box and the surface of w = 0. Given a grid box with center x c , the surface of points, x w of w = 0, satisfies Rather than the expense of calculating exact grid box volume fractions either side of this surface, the signed distance of each grid box vertex, x v , to the surface is calculated: The volume fraction either side of the w = 0 surface is approximated by the average distance to vertices in one side of the surface divided by the average distance to vertices on both sides: This gives smooth fields that will be shown in section 4.

Evaluation Test Cases
Two test cases are developed to tune parameters and to evaluate the use of two-fluid equations to represent dry, two-dimensional, subgrid-scale convection. The first is transient; the standard rising bubble of Bryan and Fritsch (2002). The second is steady state with heat transfer at the ground and uniform radiative cooling (radiative convective equilibrium [RCE]). For both test cases, the reference solution is a single-fluid solution of the two-dimensional Boussinesq equations at high horizontal and vertical resolution. These will be compared with the two-fluid model run using one column. The aim is to reproduce the vertical heat transport, the mean velocity of the rising air, the buoyancy of the rising air, and the pressure differences between fluids in one column of the two-fluid model.

Rising Bubble
The two-dimensional test case of Bryan and Fritsch (2002) was designed for fully compressible inviscid models but similar solutions are obtained solving inviscid Boussinesq equations ( = 0, = 0). The Boussinesq version has no background stratification in a domain of height 10 km and width 20 km initially at rest. The buoyant perturbation is given by for ) 2 , x c = 10 km, z c = 2 km, and x r = z r = 2 km. The 100 m grid spacing is used in the x and z directions for the two-dimensional reference simulation, and Δz = 100 m is used in the vertical for the one-column simulations. All simulations use a time step of 2 s.

Single-Fluid Solutions
The solution of the Boussinesq equations in two dimensions for the rising bubble at high resolution is shown in Figure 1, at the initial time, 500 and 1,000 s, with buoyancy colored, pressure contoured, and the w = 0 contour dotted. The solution is very similar to the fully compressible solution (Bryan & Fritsch, 2002). The pressure gradients accelerate the flow beneath the bubble, decelerate the flow above the bubble, and pull in air horizontally. The w = 0 contour at the initial time shows the initial tendency of w, which demonstrates that there is no straightforward relationship between buoyancy and w; w is the solution of an elliptic problem involving the momentum equation and the divergence-free constraint. After 500 s the center of the bubble has risen further than the edges, and after 1,000 s there is a complex region of entrainment inside the bubble consisting of both buoyant and neutral rising air and buoyant and neutral sinking air. It is unlikely that the simple parameterizations proposed here will be able to reproduce statistics of this complex behavior. Thuburn et al. (2019) and Weller and McIntyre (2019) used the same pressure for both fluids with differences parameterized as drag. The high pressure above and low pressure below the bubble in Figure 1 are a manifestation of form drag and will act to slow the ascent. However, the pressure patterns in Figure 1 also show the limitations of this approach; the low pressure region at the bottom of the bubble will also act to accelerate the rising flow into the bubble, which is not accounted for by drag, and there are also widespread pressure anomalies in the sinking fluid, which will accelerate the descent whereas drag will always act to slow the fluid.
Large-scale models of the atmosphere often have coarse horizontal resolution and relatively fine vertical resolution. The high-resolution solution is horizontally averaged, conditioned based on the sign of w, to give the vertical profiles that coarse resolution models or parameterizations should reproduce. Sinking air and air with w = 0 are assigned to Fluid 0 and rising air to Fluid 1, and an average is taken over the x direction at time t = 1,000 s to give the profiles shown in the top row of Figure 2. 1 is the volume fraction of rising fluid at every height and b 0 and b 1 , w 0 and w 1 , and P 0 and P 1 are the average buoyancy, vertical velocity, and pressure of the rising and falling fluids. 1 varies little with height. The complex shape of the bubble is evident in b 0 , which has two peaks; the lower peak associated with re-entrained air. The fastest ascent is at the bottom of the bubble, colocated with a large drop in pressure.
In order to see what happens in a single fluid atmospheric model at modest horizontal resolution without a parameterization of convection, the bubble is simulated at coarser horizontal resolutions, using 9, 5, 3, and 1 columns instead of the 200 columns of the fine resolution to cover the 20 km wide domain and a vertical space of Δz = 100 m. The conditional and horizontal averages of these simulations are shown in Figure 2.
These simulations are initialized by conservative horizontal averaging from the full-resolution solution. As the horizontal resolution is reduced, the vertical transport reduces until there is no movement at all for the one-column simulation: There is no possibility of fluid rising because in one column with a single-fluid there can be no compensating subsidence. The lack of vertical transport at coarse resolution motivates convection parameterization.

Results of the Two-Fluid, One-Column Simulations
Two-fluid, one-column simulations are initialized by horizontally and conditionally averaging the high-resolution single-fluid initial conditions. Sinking air is put into fluid zero and rising air into fluid one. The velocity field after one time step is used for initialization as the air is stationary at t = 0. Results after 500 s are shown in Figures 3 and 4 for various assumptions about transfers between fluids and various parameter values.
The solutions in the top row of Figure 3 have no drag between fluids, almost no transfers, no diffusion, and the pressures of both fluids are assumed equal meaning that the multifluid equations should be unstable (Thuburn et al., 2019). Small transfers are applied to prevent either i from becoming negative. This is sufficient to stabilize the solution but does not prevent oscillations from developing in all fields. The buoyant fluid has risen little from the initial conditions. 1 oscillates between zero and one meaning that only one fluid is present at some locations, so there is no possibility for the fluids to move past each other. The two fluids have mixed due to the two-way stabilizing transfers, which have reduced the buoyancy in the rising fluid and increased the buoyancy in the sinking fluid.
Including entrainment and detrainment as S ij = − ∇ ·u i stabilizes the solution and leads to a realistic rising volume fraction (second row of Figure 3). pressure differences between the fluids. However, there is too much mixing between the fluids and a discontinuity arises at the bubble leading edge; the buoyancy force is producing vertical motion, but there is no pressure to make it smooth.
Using entrainment and detrainment set by divergence (S ij = − ∇ ·u i ) and pressure difference between the fluids controlled by divergence ( P i = P − ∇ · u i + ∑ ∇ · u ), the simulations are more realistic (bottom three rows of Figure 3). The first value, = 2.3 × 10 4 m 2 s −1 , is set from Equation 35 using length scale equal to the initial bubble radius and a buoyancy scale equal to the initial maximum buoyancy. This makes the velocity profile too smooth and the updraft too weak, so the bubble does not rise enough. It is reasonable to consider smaller values of because the buoyancy does not retain its initial maximum for very long. Using = 10 4 m 2 s −1 gives more accurate solutions whereas using = 4 × 10 3 m 2 s −1 leads to the leading edge of the bubble rising too quickly and does not fully remove the discontinuity in the updraft velocity at the leading edge of the bubble. The pressure differences between the fluids have some similarity with the resolved simulation, but none of the values of reproduce the pressure dip at the trailing edge of the bubble. = 10 4 m 2 s −1 gives the closest approximation to the resolved simulation.
The results in the top row of Figure 4 use divergence transfer and pressure differences between fluids but, different from previous simulations, also use drag with a coefficient satisfying C D /r c = 1/2,000 m −1 . The drag slows the updrafts without increasing the mixing between the fluids, which makes the solutions less like the resolved simulations. Although drag is known to be an important term in the momentum equations for buoyant thermals and plumes (e.g., Romps & Charn, 2015), the drag is related to pressure differences caused by the relative motion and to mixing with downdraft air. The mixing with downdraft air is already accounted for in the entrainment, and the pressure differences are already accounted for in the pressure parameterization. Therefore, the drag parameterization is not a useful contribution to the multifluid Boussinesq equations.
If pressure differences between fluids are included with = 10 4 m 2 s −1 without entrainment (S ij ≈ 0) (second row of Figure 4), the fluids move past each other, but there is a buildup of Fluid 1 at the leading edge of the bubble and very little Fluid 1 below the bubble. With one column and two fluids, all air must rise in one fluid and sink in the other. Continuity does not allow for any other solution. Pressure gradients can accelerate or decelerate this flow, but as the model is set up, with no transfers between fluids, no fluid can change direction. It is therefore essential to have transfers between fluids.
The other entrainment option is the commonly used fractional entrainment rate, = 0.2/r c . This is used in combination with pressure differences between fluids in the third row of Figure 4. The solutions are smooth, but the rising volume fraction is not right with a buildup of Fluid 1 at the bubble leading edge, the updraft is too weak, and too much is entrained into Fluid 1 above the bubble. Using this entrainment, the updraft will entrain more and more air until 1 = 1. The S ij = − ∇ ·u i model gives results closer to the reference solution.
Two alternative models for estimating the buoyancy of the fluid transferred between updraft and downdraft have been tested. The results in the fourth row of Figure 4 assume that the fluid transferred has zero buoyancy in order to keep the positively buoyant air in fluid one and the negatively buoyant air in Fluid 0. Thus, when zero buoyancy air leaves fluid one (the updraft), the mean buoyancy of fluid one increases. This leads to higher buoyancy in fluid one than when using b T i = b i and gives results close to the reference solution. The results in the bottom row of Figure 4 following Thuburn et al. (2019) so that the fluid transferred is influenced by the fluid it is moving toward as well as where it has come from. These solutions are also similar to the resolved solutions.
From the results at 500 s, we can conclude the following: • Divergence transfer (entrainment and detrainment) leads to realistic updraft area fractions. • Pressure differences between fluids are usefully modeled by divergence, with a coefficient close to the derived value from scale analysis. • Drag terms are not useful. • Entrainment calculated as the reciprocal of plume width is not useful.
In order to examine further the effects of different assumptions about the buoyancy of the fluid transferred, a selection of results at 1,000 s is presented in Figure 5. None of the simulations reproduce the complex buoyancy pattern at the trailing edge of the bubble associated with a deep low pressure. The rising air fraction ( 1 ) at and below the bubble is too low in all simulations. The top row, which uses b T i = b i (transferring the mean buoyancy), shows too much mixing between the fluids-too much of the most buoyant air has been transferred out of Fluid 1, but the leading edge of the bubble has risen about the right amount. In the other two simulations, the air that is transferred out of Fluid 1 has lower buoyancy than the mean in Fluid 1 (b T 10 < b 1 ), so the remaining rising fluid stays very buoyant and rises too far. The assumption that the buoyancy transferred is zero (middle row of Figure 5) means that the buoyant fluid increases its buoyancy above the initial conditions, which would not be possible in the resolved single-fluid solution.
It might be possible to improve aspects of this simulation by optimizing and b (the coefficient from (31) for calculating b T i ), but it is probably not helpful to over fit the two-fluid model to this dry, two-dimensional test case. There are clearly discrepancies between the single-fluid resolved rising bubble and the two-fluid one-column rising bubble but rather than trying to fit this test case closely, cases with multiple plumes in each column should be considered. Using the lessons learned from modeling the rising bubble, a statistically stationary case with many plumes is simulated next.

Radiative-Convective Equilibrium
A two-dimensional, dry RCE test case is devised to mimic some properties of atmospheric convection, in particular the difficulty resolving the flow at coarse resolution. The domain is 160 km wide and 10 km tall and is resolved by 640 cells in the horizontal and 40 in the vertical (Δx = Δz = 250m). The top and bottom boundaries are zero velocity, and the lateral boundaries are periodic. A heat flux of h = 10 −3 m 2 s −3 is imposed at the bottom boundary leading to a boundary condition of b/ z = −h/ , where = 100 m 2 s −1 is the buoyancy diffusivity. The top boundary has b/ z = 0. Diffusion is applied to the momentum equation with a coefficient = 70.7 m 2 s −1 so as to give a Prandtl number of 0.707, the value for dry air. Uniform cooling of −1 × 10 −7 m s −3 is applied to the domain to maintain equilibrium.
The fluid is initially stationary with b = 0. In order to quickly and reproducibly initialize instability, b = 10 −4 m 2 s −1 is imposed in a square of side length 2,000 m at the center in the x direction and at the ground. The simulation is run using a time step of 5 s. A quasi-steady, chaotic state is reached after about 3 × 10 4 s; the simulation is run for 2 × 10 5 s, and conditional averages are calculated over the final 10 5 s.
This test case is designed to have strong updrafts in narrow plumes and weak descent elsewhere to mimic atmospheric convection but without the complication of moisture, phase changes, and three spatial dimensions.

Single-Fluid Solutions
A snapshot of the buoyancy, vertical velocity, and pressure of the RCE test case at 2 × 10 5 s are shown in Figure 6 showing intense, narrow updrafts. The buoyant air at the ground rises in narrow plumes, Figure 6. Buoyancy, vertical velocity, and pressure at t = 2 × 10 5 s of the resolved radiative-convective equilibrium test case. Vertical velocity contours are black every 0.5 m/s with negative contours dashed and the zero contour dotted. The domain is 10 km high and 160 km wide and the z direction is stretched a factor of 2 in the plots. accelerating due to buoyancy and pressure gradients then decelerating due to pressure gradients before reaching the top. The plumes spread out, and the sinking air is a mixture of warm and cold with similar pressure gradients in the rising and falling air, accelerating then decelerating the flow so as to maintain continuity.
Horizontal averages conditioned on w and time averaged every 1,000 s between 10 5 s and 2 × 10 5 s are shown as dashed lines in Figure 7. The rising and falling volume fractions, 1 and 0 , remain close to 1 2 throughout the depth. The rising plumes look narrow in Figure 6, but there are wide regions of slowly rising air around the plume cores. The entrained rising cool air and warm air sinking that has hit the model top mean that the rising and falling air have similar average buoyancy away from the ground. The ascending air gradually accelerates toward the middle of the domain and then decelerates toward the model top. At the ground the acceleration is forced by buoyancy, and away from the ground the acceleration is controlled by pressure gradients. The descending air mean velocity is close to equal and opposite to the ascending air, as expected since 1 = 0.5.

Results of the Two-fluid, One-Column Simulations
The two-fluid model is used to simulate the RCE test case in one column of 40 cells with zero buoyancy gradient at the top and a buoyancy gradient of 10 −5 s −2 at the ground for both fluids and uniform radiative cooling. The two fluids are initialized using the time mean conditional horizontal averages from the single-fluid resolved simulation. The entrainment model S ij = − ∇ ·u i is used for all simulations to account for transfers throughout the domain to represent cloud base mass flux (at the bottom boundary), cloud top detrainment, and lateral entrainment. is set from Equation 35 using the depth of the high-resolution solution super adiabatic layer as the length scale (L = 800 m) and the buoyancy scale B = L b/ z = 8 × 10 −3 m s −2 giving = 2,000 m 2 s −1 . A simulation using S ij = − ∇ ·u i , no drag, and b T i = b i is shown in the top row of Figure 7. The entrainment model S ij = − ∇ ·u i keeps i at about 0.5 throughout the depth, as in the resolved simulation, with entrainment in the lower half of the domain (S 01 > 0) and detrainment in the upper half (S 10 > 0), where the air decelerates toward the top. The one-column model reproduces the superadiabatic layer near the ground and the near uniform buoyancy away from the ground, but the buoyancy difference between the fluids is too large. The updraft and downdraft velocities are captured well. The pressure differences between fluids are too large, which is consistent with the over estimated buoyancy difference between fluids. Figure 7. This model of the buoyancy of the transferred fluid is not helpful; it leads to strong cooling in fluid zero near the ground due to the large transfers. It is more accurate to assume that fluids take their mean properties when transferred,

Results of a simulation using
The two-fluid, one-column model of RCE gave results closer to the reference solution than for the rising bubble. This is likely to be because the rising and falling fluids in the RCE test case are averaged over a large area with plumes at different stages of development rather than trying to simulate the complex entrainment pattern behind a single bubble. The RCE test case is horizontally homogeneous and in equilibrium making it more straightforward.

Summary, Conclusions, and Further Work
This paper demonstrates the plausible utility of using multifluid equations to represent subgrid-scale convection. Solutions of the two-fluid Boussinesq equations using coarse horizontal resolution (one column) reproduce aspects of highly resolved single-fluid simulations of dry, two-dimensional convection. Two test cases are used to evaluate the multifluid model. The first is a Boussinesq version of the well-known two-dimensional dry rising bubble of Bryan and Fritsch (2002), and the second is a dry, two-dimensional version of RCE with heating at the ground and uniform radiative cooling. High-resolution single-fluid simulations are conditionally averaged based on the sign of vertical velocity (w) in order to define rising and falling fluids, which each have their own buoyancy, velocity, and pressure. The two-fluid, one-column Boussinesq model reproduces mean vertical heat and mass transport of these test cases.
The RCE test case complies with the assumptions of a traditional mass-flux convection scheme-an ensemble of plumes per column, horizontal homogeneity, and equilibrium. The rising bubble test does not satisfy these constraints. The multifluid model gives good solutions for the RCE test case and good solutions for the bubble for the first part of the simulation but does not represent the complex re-entrainment as the bubble continues to rise.
The multifluid Boussinesq equations can have divergence in each fluid while the divergence averaged over all fluids is constrained to be zero by the single-fluid Boussinesq equations. A model is proposed and evaluated for the pressure difference between the rising and falling fluids based on the divergence of each fluid. The use of a different pressure for each fluid is not only realistic but also leads to stable solutions of the multifluid equations.
Conditionally, averaging the high-resolution solutions based on w leads to rising and falling fluids with area fractions, i , which are approximately uniform with height. Models for entrainment and detrainment into and out of these regions need to be consistent with the definition of the convecting region. A model for entrainment and detrainment based on divergence is proposed that leads to excellent agreement with from the high-resolution solutions. A more conventional fractional entrainment rate based on the reciprocal of the plume width is also tested, which leads to unrealistic rising area fractions.

10.1029/2019MS001966
The properties of the fluid that are transferred from one fluid to another are unknown in the multifluid equations, and these require closure. An obvious first guess is that fluid takes its mean buoyancy and velocity with it when it is transferred. Improvements on this assumption will depend on knowledge of subfilter-scale variability. However, since conditional averaging is based on the sign of w, it is consistent to assume that the fluid transferred has zero vertical velocity (w T i = 0). Three models for the buoyancy of the fluid transferred are tested. The most successful is to assume that the fluid takes its mean buoyancy with it (b T i = b i ). Assuming that the buoyancy of the transferred fluid is an average of the buoyancy of both fluids can lead to unrealistically large or small, unbounded buoyancies. Assuming that one fluid has negative buoyancy and the other positive buoyancy and fluid of zero buoyancy is transferred (b T i = 0) also leads to unbounded buoyancy.
One-column multifluid modeling with moisture has already been demonstrated by Thuburn et al. (2019) and two-dimensional simulations of the multifluid Euler equations by Weller and McIntyre (2019). Next steps are to simulate convection in the gray zone (where it is partially resolved) and to represent subfilter-scale variability within each fluid, further extending the EDMF scheme (Tan et al., 2018). Cloud base mass flux could then be estimated based on higher-order moments of the vertical velocity and buoyancy within the environment fluid and vertical transport of turbulence by plumes can inform estimates of entrainment. This will lead to a unified parameterization of turbulence and convection and will take multifluid modeling of convection beyond existing paradigms.