A two‐fluid single‐column model of turbulent shallow convection. Part 1: Turbulence equations in the multifluid framework

The multifluid equations are derived from the compressible Euler equations (or any of the usual approximate equation sets used in meteorology) by conditional filtering. They have the potential to provide the basis for an improved representation of cumulus convection and its coupling to the boundary layer and larger scale flow in numerical models. The present article derives the prognostic equations for subfilter‐scale turbulent second moments in the multifluid framework, along with certain systematic simplifications of them, thus providing a multifluid analogue of the well‐known Mellor and Yamada hierarchy of turbulence closures. As well as enabling a more accurate calculation of subfilter‐scale fluxes and the effects of subfilter‐scale variability on cloud fraction, liquid water, and buoyancy, the second moment information can be used to obtain a more accurate parameterization of entrainment and detrainment. A subset of the turbulence equations derived here is employed in the two‐fluid single‐column model described in Part 2 and applied to the simulation of shallow cumulus cases in Part 3.

the model as well as its possible future extensions and generalizations.
Atmospheric models are often run at grid resolutions too coarse to resolve individual convective updrafts explicitly. In this case, some of the effects of convection can be taken into account-parameterized-by including additional equations in the model that describe the dynamics and thermodynamics of the updrafts and their effects on the resolved-scale flow. Such parameterizations generally involve a number of simplifying assumptions and approximations. Mass flux convection schemes (e.g., Arakawa and Schubert, 1974;Tiedtke, 1989;Gregory and Rowntree, 1990;Zhang and McFarlane, 1995;Kain, 2004;Gerard and Geleyn, 2005;Park, 2014) are a widely used example.
By applying conditional filtering (see Section 2) to the continuous governing equations, Thuburn et al. (2018) derived a mathematical framework with a full set of prognostic equations for the convective updrafts and for their environment. If desired, multiple updraft types and downdrafts can be included as additional fluid types with their own equation sets. Thuburn et al. (2018) referred to the resulting framework as the "multifluid" equations because of their resemblance to the equation sets used in engineering to model multiphase flows, though their interpretation is somewhat different. Because they pick out flow structures such as updrafts, the multifluid equations are conceptually related to mass flux schemes (see the discussion in Thuburn et al. (2018)) as well as the more recently developed eddy-diffusivity mass-flux (EDMF) schemes (e.g., Soares et al., 2004;Siebesma et al., 2007;Neggers, 2009;Angevine et al., 2010;Witek et al., 2011;Sušelj et al., 2012Sušelj et al., , 2013 when subfilter-scale fluxes are included. However, because the multifluid framework involves virtually no approximation, it avoids several assumptions and approximations that are inherent in traditional mass flux and EDMF schemes and may limit their applicability. Such approximations include neglecting transience in the updraft equations and hence any memory of the state of the updraft from one time step to the next, neglecting horizontal advection of updrafts, and imposing environmental subsidence in the same grid column, so requiring the parent dynamical core to compensate if the subsidence should be remote from the updraft. The present work is motivated by the possibility of developing a three-dimensional multifluid scheme that can represent cumulus convection in weather and climate models more accurately and across a wider range of model resolutions than traditional schemes. The multifluid framework includes various terms that need to be parameterized in any practical application, namely subfilter-scale fluxes, effects of subfilter-scale pressure fluctuations, and terms representing the relabelling of fluid as it is entrained and detrained. The subfilter-scale fluxes may be modelled in a simple way as a downgradient eddy diffusion with a specified or stability-dependent profile; however, such formulations can lead to behaviour that is either unphysical or difficult to handle numerically (e.g., Thuburn et al., 2019). It may be advantageous to model such terms using a higher-order closure, for example via a scheme that derives eddy diffusivity from predicted turbulent kinetic energy (TKE: e.g., Witek et al., 2011;Tan et al., 2018). Our initial experiments with a TKE-based scheme led to better behaviour than with a specified eddy diffusivity profile, and encouraged us to pursue this route. Predicted or diagnosed subfilter-scale turbulence information can potentially be used in several ways to improve parameterized processes, including, for example, countergradient contributions to subfilter-scale fluxes such as those associated with the fallback of overshooting thermals, the effects of subfilter-scale variances on cloud fraction, liquid water amount, and buoyancy, and the "sorting" effects of entrainment and detrainment. In order to underpin such developments, the purpose of this article is to derive and document the prognostic equations for second-order turbulence quantities in the multifluid framework and certain systematic simplifications of them, thus providing a multifluid analogue of the Mellor and Yamada (1982) hierarchy of turbulence closures.
Over a series of articles, Mellor (1973); Yamada (1974, 1982) derived their well-known hierarchy of turbulence closure models (see also subsequent improvements by Helfand and Labraga, 1988;Nakanishi, 2001;Janjić, 2001, Nakanishi and Niino, 2004, 2006, 2009. They ensemble-averaged and manipulated the continuous governing equations to obtain prognostic equations for all second-order moments of velocity and conservative scalars, such as potential temperature in the dry case or liquid water potential temperature and total specific water in the moist case. By assuming a Gaussian joint distribution for liquid water potential temperature and total specific water, they were able to estimate the effect of condensation on the buoyancy and its correlations (Sommeria and Deardorff, 1977;Mellor, 1977b;Chen, 1991). They made some reasonable assumptions about how the third-order terms, dissipation terms, and terms involving subfilter-scale pressure fluctuations that arise in the equations may be expressed approximately in terms of predicted first-and second-order quantities; in this way the equation set is closed except for a number of scalar coefficients that may be estimated from observations and experiments and the space and time distribution of a certain master turbulence length-scale. The authors then made a series of simplifying approximations, the justification of which depends on the isotropy of the flow, to obtain a hierarchy of turbulence models of varying complexity.
The multifluid approach and the high-order turbulence approach provide two different ways of accounting for unresolved variability in models, each likely to be most useful in different regimes; the multifluid decomposition, like the mass flux approach, is intended to capture the largest-scale coherent turbulent structures such as updrafts, while the high-order turbulence approach should be most accurate for smaller-scale quasi-isotropic turbulence. A multifluid analogue of the Mellor-Yamada hierarchy would provide a unified approach with the potential to work well across a wider range of regimes. In particular, it might facilitate the construction of a "scale-aware" model that is applicable at a range of different resolutions encompassing the grey zone; the "grey zone" refers to those resolutions at which convective updrafts are marginally resolved, which are extremely challenging for current models. The decomposition of the flow into different fluid types in the multifluid approach explicitly takes into account some of the contributions of high-order moments (Appendix A, see also Lappen and Randall, 2001). Thus, a relatively simple high-order turbulence scheme within a multifluid model may be able to achieve comparable accuracy to a more complex high-order closure in a single-fluid model. As such, there appears to be considerable scope to explore the cost versus accuracy trade-offs for different numbers of fluids and different high-order closure approximations.
In this article we generalize the work of Mellor and Yamada in three ways. First, Mellor and Yamada restricted attention to the Boussinesq equations, both for simplicity and because their focus was on the atmospheric boundary layer. Their turbulence closures are usually employed in Boussinesq form, even in large-scale models that solve compressible governing equations. Many modern atmospheric models are based on the fully compressible Euler equations in order to be valid for a wide range of spatial scales, from convection-resolving to global. To ensure that the turbulence model may be used in a fully consistent way with a fully compressible parent model, here we contruct fully compressible analogues of Mellor and Yamada's equations. The compressible case involves mass-weighted filtering or Favre filtering (Favre, 1969) of certain terms, but otherwise the resulting equations for first and second moments have a very similar mathematical structure to those in the Boussinesq case. Second-moment equations or kinetic energy equations for the compressible case have been presented previously (e.g. Adumitroaie et al., 1999;Aluie, 2013), though that work was largely aimed at modelling high Mach number turbulence; here our focus is on buoyancy-driven low Mach number meteorological flows and their consistent coupling with the larger scale.
Second, Mellor and Yamada focused on a limited class of filters, including ensemble averages and global horizontal averages, that satisfy the Reynolds rules where the overbar indicates a filtered field and ′ = − . Hence they derived equations for second-order quantities of the form ′ ′ . One of the motivations for the present development is its possible application to the convective grey zone, in which turbulence is highly inhomogeneous on the model grid scale. In this case it is more appropriate to consider a spatial filter of finite width, comparable to the grid scale, in which case the Reynolds rules do not apply. Nevertheless, second moments such as ′ ′ may be replaced by generalized second moments such as − (with analogues for higher moments); the resulting equations then retain the same algebraic structure as those considered by Mellor and Yamada (Leonard, 1975;Germano, 1992, see also Appendix A below), though the parameters in the turbulence model need to be modified when the dominant turbulent eddies become partly resolved (Ito et al., 2015).
Third, and most significantly, we extend Mellor and Yamada's hierarchy to the conditionally filtered case. All of the terms in the Mellor-Yamada formulation have analogues in the multifluid case. In addition, new terms arise representing the effects of relabelling of fluid from one type to another, that is, entrainment and detrainment. New terms also arise because of the requirement to keep a single filtered pressure field for all fluid types. The multifluid second-moment equations are derived in Appendix B and presented in Section 3.
In Section 4, closures are presented for the terms in the second-moment equations that need to be modelled. The closures are based closely on those of Mellor and Yamada (1982), with minor modifications to account for changes in density, the fractional volumes occupied by different fluid types, and the fact that the conditionally filtered velocity may be divergent even for incompressible flow. The fact that the flow regime of interest here is the same as that of interest to Mellor and Yamada (1974);Mellor and Yamada (1982) suggests that no major changes should be needed to the Mellor and Yamada modelling to account for compressibility. On the other hand, it is not at all clear that the length-scales and distributions of the fluid labels that identify updrafts, downdrafts, and environment (see Section 2) can safely be ignored in formulating these closures. For now, in order to make progress, we ignore these factors. However, justification for this assumption, or modifications to it, must come from experience in application of multifluid higher-order closures in practice, and/or by direct conditional filtering of high-resolution reference simulations. The closures presented in Section 4 are then systematically simplified in Section 5 to give a multifluid analogue of the Mellor-Yamada hierarchy.
In the absence of convective updrafts, for example in a stable boundary layer, a multifluid scheme might reasonably allow several fluid types to exist, but with identical vertical profiles of fluid properties for each of those types. The requirement to keep those profiles identical puts some useful constraints on the form of certain terms. First, it shows that the approximation to neglect transience and advection applied in Section 5 must be applied to the advective form, rather than the flux form, of the second-moment equations. Second, it shows that subfilter-scale fluxes of any quantity must be accompanied by appropriate relabelling terms. Section 6 derives the required form of those relabelling terms. Such terms can be important, for example, during a spin-up phase prior to the onset of convective overturning.

CONDITIONAL FILTERING
In this section we briefly recall the idea of conditional filtering and, for later reference, present the conditionally filtered governing equations for first-order quantities. For more details see Thuburn et al. (2018). To aid the reader, Tables E1-E5 in Appendix E summarize the key notation used in this article.
Here we consider the fully compressible Euler equations. (Adapting the procedure for any of the usual approximate equation sets used in meteorology is straightforward.) Here, is the total fluid density, u = (u, v, w) is the fluid velocity, p is pressure, and Φ is geopotential. For simplicity, the governing equations have been expressed in terms of "conservative" variables, the specific entropy and q the total specific water content, with sources S and Q, respectively. Coriolis terms have also been omitted, but it is straightforward to include them. The equation of state has been written in the generic form (Equation 6), which accounts implicitly for any latent heating effects provided moisture is assumed to be in local thermodynamic equilibrium (i.e., no supersaturation or supercooled water, no condensed water in subsaturated air). Conditional filtering combines a spatial filtering operation, indicated by an overbar, with a set of n quasi-Lagrangian indicator functions I i , i = 1, … , n. The filter is assumed to commute with space and time derivatives. At any point in the fluid, one of the I i is equal to 1 while the others are equal to 0. The I i are assumed to satisfy The I i may be used to pick out different regions of the fluid, such as updrafts or their surrounding environment. Because I i is discontinuous, the relabelling term r i will have the form of a Dirac delta function, though this will be smoothed by the application of the filter. 1 Define i to be the volume fraction of the ith fluid type on the filter scale: Then Also, define the average density of the ith fluid type on the filter scale i by Combining the mass continuity equation (Equation 2) with the indicator function equation (Equation 7) gives Application of the filter and use of Equation (10) then gives the conditionally filtered density equation: where the mass-weighted filter-scale velocity u i is defined by i i u i = I i u. Here, for later convenience in parameterizing these terms, we have written where  i is the filter-scale rate per unit volume at which mass is relabelled from type to type i. In a similar way, application of the conditional filter to the scalar, momentum, and thermodynamic state equations (Equations 3-6) gives Here, i and q i are the mass-weighted filter-scale specific entropy and total specific water in fluid i, defined by i i i = I i , i i q i = I i q, with their filtered sources S i and Q i defined analogously. The subfilter-scale flux of entropy F i SF is related to the subfilter-scale u-covariance C u i and is defined by with an analogous expression for the subfilter-scale flux of water, where ( ) i is defined by i i ( ) i = I i . Similarly, the subfilter-scale flux of momentum is The pressure gradient term is treated slightly differently in order that a single filter-scale pressure p appears in the conditionally filtered equations. All departures of p from p are accounted for by the  i terms. (See Thuburn et al. (2018); Thuburn and Vallis (2018), but note also that Weller et al. (2020); Shipley et al. (2022) make an argument for parameterizing pressure differences between fluids.) Thus, is the subfilter-scale contribution to the pressure gradient force on fluid i. Note that the definition implies ∑ i  i = 0, as expected for momentum conservation (Thuburn et al., 2018). Subfilter-scale contributions to the equation of state are written as P i In large-scale models, subfilter-scale contributions to the equation of state are usually neglected. In high-order turbulence modelling, on the other hand, it is usual to account for subfilter-scale variations in q and to determine the cloud fraction, amount of liquid water, and filter-scale buoyancy (e.g., Sommeria and Deardorff, 1977;Mellor, 1977b, see also Appendix D below). Finally, hatted variables such aŝi represent the average values of that quantity in fluid that is relabelled from type to type i. It is common in convection schemes and EDMF schemes to approximate entrained and detrained quantities bŷi ≈ for a generic variable . However, to maintain generality we do not make this approximation here. Indeed, it is widely accepted that there is preferential detrainment of the least buoyant air from an updraft (e.g., Raymond and Blyth, 1986), while mass-flux closures based on convective inhibition (e.g., Fletcher and Bretherton, 2010) implicitly assume preferential detrainment of the lowest-w air at the boundary layer top. There is growing evidence that other quantities are also partially "sorted" by entrainment and detrainment (e.g., Romps, 2010;Thuburn et al., 2019). In Part 3 (McIntyre et al., 2022), we find significant sensitivity to how the quantitieŝi are parameterized in a single-column model.
In the conditionally-filtered first-moment equations (Equations 12 and 14-16), the terms  i ,̂i , and  i generally need to be parameterized. The subfilter-scale fluxes may be parameterized, or they may be predicted or diagnosed using the equations given in Sections 3-5 below. The specific formulations for all these terms used in our single-column model are detailed in Part 2.

MULTI-FLUID EQUATIONS FOR SECOND-ORDER MOMENTS
In this section we present prognostic equations for second-order quantities such as subfilter-scale momentum fluxes, scalar fluxes, and scalar variances, for the fully compressible equations. The derivations neglecting molecular viscosity and diffusion are given in Appendix B, while Appendix C outlines how to include those processes. The notation is also defined in full there, 2 but the main new 2 In Mellor and Yamada (1974); Mellor and Yamada (1982) and much of the other literature on high-order closures, vector quantities such as u and tensor quantities such as uu are represented by their components denoted by subscripts i, , k and so forth. Here, since subscripts are already used to indicate fluid types, further indices to indicate vector and tensor components would be unwieldy and risk making the presentation unintelligible. We therefore stick to bold font notation for vectors and sans serif for tensors. However, for the reader interested in details, some care is needed to keep track of the implied tensor indices in quantities such as T uu i . notation introduced is as follows: T uab i , T uua i , and T uuu i are third-order turbulent fluxes (Equations B12 and A15); B b is a modified definition of the covariance C ab when the first field is , andB pb is a modified definition of C ab when the first field is p (see Equation B16 and accompanying text);  ab i is a relabelling term (Equation B13);  ab i is a molecular dissipation term (Equation C3).
Dissipation terms associated with molecular viscosity or diffusion are included here. However, molecular diffusive contributions to the fluxes of second-order moments are expected to be negligible and so are omitted (e.g., Mellor and Yamada, 1974). In the conditionally filtered framework, further molecular diffusion terms arise due to diffusion across the boundary of fluid i (Appendix C); however, these terms are also expected to be negligible and are omitted here.
The equations are presented here in advective form, since the advective form is the basis for some of the simplifying approximations made in Section 5. If needed, the corresponding flux forms are easily obtained with the aid of the conditionally filtered mass Equation (12).
We first present the second-moment equations, and then below discuss their physical meaning and similarities to and differences from the single-fluid case. The scalar covariance equation, taking and q as prototypical scalars, is Here D i ∕Dt = ∕ t + u i ⋅ ∇ is the "material" rate of change following the velocity u i . Replacing q by gives a prototypical scalar variance equation: There is an analogous equation for the variance of q.
The equation for the subfilter-scale flux of a scalar such as entropy is There is an analogous equation for the subfilter-scale flux of q.
The equation for the subfilter-scale momentum flux is Here [ ] indicates the transpose of the immediately preceding term. Contraction of the tensor equation (Equation 24) gives the equation for the subfilter-scale kinetic energy: Here contracted on the last two indices. Also, while  k i = 1 2 Tr( uu i ). It will be convenient for later use to write down the deviatoric (i.e., trace-free) part of Equation (24): Here I is the 3 × 3 identity matrix. Many of the terms in these second-order moment equations resemble those in the single-fluid second-order moment equations used by Mellor (1973); Mellor and Yamada (1974), with straightforward modifications to allow for density variations (Adumitroaie et al., 1999;Aluie, 2013). The terms involving D i ∕Dt represent material advection of the second-order moments. Terms involving T represent the divergence of third-order turbulent fluxes. There are flux-gradient terms involving F or F that represent a source of the second-order moment when the flux is downgradient and a sink otherwise. For certain second-order moments, the molecular dissipation, represented by the  terms, cannot be neglected. Terms such as representing correlations with entropy and moisture sources, are included straightforwardly.
For the variable density case, the equations appear most simple if a slightly different definition is used for correlations involving pressure (B) or density (B) compared with other correlation terms (C) (e.g., Adumitroaie et al., 1999;Aluie, 2013, see Appendix B for details). The terms involving B always appear multiplied by (1∕ i )∇p, so in fact they represent buoyancy-correlation terms. Thus, these terms also have analogues in the single-fluid case.
The terms involvingB represent gradients of pressure-correlation terms, or correlations of pressure with gradients. Again these terms have analogues in the single-fluid case; however, for a term such as i ∇B ), the factor of i means that this term cannot be written as the divergence of some quantity and hence cannot be merged with the divergence of the third-order flux, as is done in the single-fluid case.
Some terms arise in the multifluid case that have no analogues in the single-fluid case. Terms involving  i result from fluid relabelling. These terms may be simplified if the properties of relabelled fluid may be assumed equal to the mean properties of the source fluid, so that ( ) i = ( ) (e.g., Equation B13), but for completeness here we keep the more general form; see the discussion at the end of Section 2. Finally, terms denoted  arise to account for certain subfilter-scale fluctuations, because we retain a single resolved-scale pressure p for all fluid types.

CLOSURE ASSUMPTIONS
As noted in Section 1, our interest is in low Mach number buoyancy-driven meteorological flows, such as the boundary layer and dry and moist convection. It is therefore expected that the modelling of subfilter-scale terms proposed by Mellor (1973); Mellor and Yamada (1974); Mellor and Yamada (1982) should not require any major modification to account for compressibility (e.g., Bilger, 1975, and references therein). It is much less clear whether their approach might require retuning, or even a more fundamental reformulation, to take into account the distribution of fluid labels in the conditionally filtered case. For example, Ito et al. (2015) found that it was necessary to reduce the length-scales in the Mellor-Yamada model when the grid resolution and filter scale were fine enough to partially resolve the dominant eddies. This suggests that a similar reduction of length-scales might be needed in a multifluid Mellor-Yamada model when the variations in fluid labels begin to be resolved. For the rest of this article, we assume that the Mellor and Yamada framework can be used, with possible recalibration of length-scales but otherwise minimal modification. However, this assumption should be tested.
To begin with, we examine the term B u i . When multiplied by ∇p∕ i , it gives, essentially, the subfilter-scale buoyancy flux in fluid i. In the single-fluid dry Boussinesq system, the subfilter-scale buoyancy flux is usually taken to be proportional to the subfilter-scale potential temperature flux, which is itself predicted or diagnosed, and so does not need to be modelled. In the moist case, the effects of condensation on the buoyancy flux must be taken into account. Following Sommeria and Deardorff (1977), it was shown by Mellor (1977b) (see also Chen, 1991) that if the joint distribution of entropy (or liquid water potential temperature) and total specific water can be approximated as Gaussian (see Lopez-Gomez et al., 2020, for an alternative log-normal approximation) and the saturation specific humidity q s is linearized about ( , q s ( )), then the integrals required to compute the buoyancy flux can be evaluated. Essentially the same derivation carries through here for a more general equation of state, provided we assume − g is proportional to q l , where g is the density the air would have if all water remained in vapour form, and q l is the specific liquid water. The multifluid nature of the equations does not affect this calculation. However, when the filter is of finite spatial scale, the Gaussian joint distribution assumption should be interpreted as taking into account the spatially varying filter weight.
Under these assumptions and approximations (see Appendix D), it may be shown that with = q or = , where  is the cloud fraction and Equation (28) defines D q i and D i . The cloud fraction is given by where q = Q s ( ) defines the saturation curve at a given pressure, and 2 s is the standard deviation of supersaturation in fluid i given by Equation (D3). (Some subscripts i have been omitted for clarity.) Equation (28) also holds with = w, provided the joint distribution of w and supersaturation is assumed Gaussian (and similarly for = u and = v). Now B u i is not directly predicted or diagnosed, but scale analysis shows that where | ′ | is the scale of the subfilter-scale density fluctuations in fluid i. Thus, to an excellent approximation, we can replace B u i by C u i , and similarly for B uq i . Thus B u i may be expressed in terms of quantities that are predicted or diagnosed: The third-order turbulent fluxes are modelled as downgradient diffusion terms: ) , ) , Here, 1 , 2 , 3 are turbulence length-scales, all assumed proportional to a certain master length-scale L, which may be time-and space-dependent. (In general, 1 , 2 , 3 , and L, as well as the length-scales l 1 , l 2 , Λ 1 , and Λ 2 introduced below, will also depend on the fluid type i; an index i to indicate this dependence has been suppressed.) U i = √ 2k i is a turbulent velocity scale for fluid i. "perm" indicates the two tensors similar to ∇C uu i , but with tensor indices permuted to make the overall expression symmetric. [ ] indicates the transpose of the immediately preceding term.
In the momentum flux equation, the pressure-velocity correlation ∇B Here, l 1 and l 2 are turbulence length-scales, again assumed proportional to the master length-scale. C 1 is another tunable constant, in this case dimensionless. In the pressure-strain correlation, Mellor and Yamada (1974) also include a possible contribution proportional to the buoyancy flux times gravity, but then set the coefficient to zero. However, Nakanishi (2001) found a better fit to LES data by retaining that term, and also included terms proportional to the momentum flux times the resolved strain. Similarly, in the correlation of pressure and scalar gradient, Mellor and Yamada (1974) also include a possible contribution proportional to (in the present notation) the correlation of the scalar with the buoyancy, but then set the coefficient to zero. Again, Nakanishi (2001) found a better fit to LES data by retaining that term, and also included a term proportional to the scalar flux times the resolved strain. The form of the first term on the right of Equation (33) is dictated by the requirement that, in the single-fluid incompressible case, both sides of the equation should be trace-free. For single-fluid low Mach number compressible flow, there is no longer a strict requirement for the left-hand side of the equation to be trace-free, but this incompressible pressure-strain model should still be suitable. In the multifluid case, on the other hand, ∇ ⋅ u i can be significant even for an incompressible fluid. We have therefore modified the second term on the right-hand side to ensure that it remains trace-free even when ∇ ⋅ u i is nonzero.
Terms involving  do not arise in the single-fluid second-order moment equations. The term  i arises in the multifluid (first-moment) momentum equation (Equation 16). Following Romps and Charn (2015), Thuburn et al. (2019) model this term as a drag on updrafts. The form of the term  a i (see Equation B16) resembles a i times  i minus a term similar in form to  i but with an additional factor a under the overbar. It is far from obvious how such a term may be modelled. However, if the fluctuations in a are small compared with its mean, which may be the case for a = or a = q, then it may be a reasonable approximation to write On the other hand, for updrafts u i may differ significantly from u; for example, the term u i ⋅  i may represent a significant sink of resolved kinetic energy and source of subfilter-scale kinetic energy, which is unlikely to cancel to leave  k i ≈ 0 in the TKE equation (Equation 25). In this case it may be possible to approximate  k i ≈ u i ⋅  i , but this needs to be investigated theoretically or in high-resolution reference simulations.
The molecular dissipation terms are modelled by Here, Λ 1 and Λ 2 are further turbulence length-scales, once again assumed proportional to the master length-scale. The molecular dissipation terms for scalar fluxes are assumed negligible, since the absence of an isotropic first-order tensor means there can be no dissipation term of the expected form (Mellor, 1973). The correlations between fluctuations and sources such as C S i and C uS i are likely to be small if the source terms are slowly varying. We will neglect such terms in the following, as do Mellor and Yamada (1982). However, a term such as C qQ i might well become important in a precipitating updraft, so further consideration should be given to how it may be approximated.
Substituting these modelling assumptions into the second-moment equations gives the multifluid analogue of the Mellor and Yamada level-4 closure.
Scalar covariance: Scalar variance: There is an analogous equation for the variance of q. Scalar flux: There is an analogous equation for the subfilter-scale flux of q. Momentum flux: Cohen et al. (2020) derive an equation very similar to Equation (36), with an eddy diffusive closure for the subfilter-scale fluxes, in their extended EDMF scheme. As noted by Mellor and Yamada (1982), planetary rotation or the presence of a solid lower boundary could suggest additional contributions to the closures discussed above. However, since results based on the above closures appeared satisfactory, and in the interest of avoiding unnecessary complexity, they did not include such additional contributions. We will do the same here.

SIMPLIFIED MODELS
Using an ordering analysis based on the anisotropy of the momentum flux tensor, Mellor and Yamada (1974); Mellor and Yamada (1982) systematically simplified their level-4 model by neglecting transience, advection, and turbulent transport of various second moments. When simplified in this way, the corresponding budget equations reduce to balances between local source and sink terms. Without repeating their ordering analysis here, we simply write down analogous expressions for the multifluid case. It is useful, first, to consider the general structure of the second-moment equations with Mellor-Yamada-type closures. Using the relation between subfilter-scale fluxes and second moments in Equations (18) and (19) and ignoring entropy and water source terms, the full set of second-moment equations may be written concisely in the form DC Dt Here C = (C , C q , C qq , … , C u , … , C uu , … , C ww ) T is a vector of the complete set of second moments, T represents the corresponding set of turbulent fluxes, and we have dropped the subscript i for clarity. The matrix S depends on the kinetic energy through the U i terms and may depend on other second moments, for example through relabelling terms; thus this system is not quite linear. Nevertheless, when the left-hand side terms are neglected for any subset of second moments, as in the simplified models discussed below, the corresponding equations effectively become a system of diagnostic equations for those second moments that are no longer predicted; in this way we always have the same number of equations as unknowns.
In the single-fluid case, the right-hand side of any second-moment equation is independent of whether that equation is written in flux form or advective form. Thus, neglecting transience and advection in the advective form is exactly equivalent to neglecting transience and advection in the flux form. In the multifluid case, however, the relabelling terms differ between the advective form and flux form. Therefore, it matters whether we work with the advective form or flux form when neglecting transience and advection.
To demonstrate this point, consider a case in which convective updrafts are absent, for example, a stable boundary layer. It would be reasonable, in this case, to require that all of the fluid types have identical properties i , i , q i and so forth and that the model should maintain this state of affairs. Now consider a generic relabelling term in the advective form equation for C ab i (Equation B13): Provided we choose a i = a =â i =â i and similarly for b, andĈ ab i =Ĉ ab i = C ab i , then the net relabelling effect on C ab i ( and likewise C ab ) will vanish, whatever the values of  i . In the flux form equation for C ab i , on the other hand, the relabelling term acquires an additional contribution, This additional contribution is exactly what is required to account for changes in mass and fluxes of mass in the flux form transience and advection terms, However, if we were to neglect these flux form transience and advection terms, then the additional relabelling terms would not be balanced and would disrupt the state of identical fluid properties in an unphysical way. For this reason, we must work with the advective form when neglecting transience and advection terms.

Level-3 model
In the level-3 model, transience, advection, and turbulent transport are neglected in the deviatoric part of the momentum flux equation and in the scalar flux equations. The isotropic part of the C uu i

equation (Equation 39) becomes
while the deviatoric part reduces to Also, the C u i equation (Equation 38) reduces to There is an analogous equation for the subfilter-scale flux of q.
The scalar variance and covariance equations (Equations 36 and 37) remain the same as at level 4.

Level-2.5 model
In the level-2.5 model, transience, advection, and turbulent transport are also neglected in scalar variance and covariance equations. Equations (36) and (37) then reduce to with an analogous equation for the variance of q. An important limitation of the level-2.5 model has been recognized by several authors (e.g., Helfand and Labraga, 1988;Janjić, 2001;Nakanishi, 2001).
The only remaining prognostic second moment in the second-moment system (Equation 40) is the TKE; all other second moments are diagnosed from the subsystem of diagnostic equations. Under certain conditions, namely when turbulence is convectively driven and is growing towards its equilibrium intensity, this subsystem can become singular, leading to unphysical values of some quantities and even causing the model to "blow up". These authors have suggested various modifications to the level-2.5 model to remove the problematic behaviour. Our experience to date suggests that similar modifications are needed in the multifluid case.

Level-2 model
In the level-2 model, transience, advection, and turbulent transport are neglected in the TKE equation (Equation 25), giving

ENTRAINMENT AND DETRAINMENT LINKED TO SUBFILTER-SCALE FLUXES
The subfilter-scale flux of a scalar, such as entropy, is given by ). The appearance of the factor i in the flux, and inside the divergence operator in the first-moment entropy equation (Equation 14), is crucial for ensuring that the flux is conservative. However, it means that the divergence of the flux has a dependence on the gradient of i .
The problem is most obvious in the scenario discussed at the start of Section 5, in which the properties of two fluids have identical profiles. If the corresponding i profiles are not also identical then the flux divergence terms will give rise to different entropy tendencies in the the two fluids and their properties will not remain identical. Figure 1 shows the problem schematically and suggests that the solution is to allow an additional relabelling that transfers the scalar conservatively between fluid types.
Let us propose a relabelling term R i of the following form: where, for clarity, the i equation has been written in advective form and all terms have been omitted except the divergence of the subfilter-scale flux and the proposed new relabelling term. The relabelling term can be thought of as arising from an equal and opposite mass relabelling  i =  i but witĥi ≠̂i.
To ensure conservation, the relabelling terms must satisfy ∑ Also, we require that, when fluids i and have identical properties, For the case of two fluids, it is straightforward to verify that has the required properties. Clearly R 1 + R 2 = 0. Also, if the two fluids have identical properties then F i SF = i F 0 for some F 0 . Then (using 1 + 2 = 1 and 1 1 + 2 2 = ), while Hence as required.
For n fluids it may be verified that the generalization has the required properties. Analogous relabelling terms are needed in the component momentum equations linked to the subfilter-scale momentum fluxes. The discussion here shows that, in the presence of nonzero subfilter-scale fluxes, the so-called turbulent entrainment term (e.g., de Rooy et al., 2013;Cohen et al., 2020), that is, the scalar transfer not carried by the mean mass entrainment  i −  i , must acquire a contribution given by Equations (51) or (55) above.
These relabelling terms will also affect the second-moment equations. Moreover, we should check that when two fluids have identical properties they have identical second-moment tendencies. The issue is complicated because the modification needed to the second-moment equations can depend on whether or not we make some of the simplifications discussed in Section 5. Let us illustrate this idea for the entropy variance; the treatment of the other second moments is analogous.
Rederiving the entropy variance equation with the new term R i in the entropy equation, and retaining only the relevant terms, gives Here we have included a possible additional variance relabelling term R i . The additional relabelling terms should be conservative, so we require First consider level-2 and 2.5 simplifications in which the left-hand side of Equation (56) is neglected. For the case of two fluids, it is easily verified that the choice ensures that identical fluid properties implies identical variance tendencies. The generalization to n fluids is . (59) If we do not neglect the left-hand side of Equation (56), as in level-3 and 4 schemes, or the kinetic energy equation at level 2.5, then gradients in i affect the ∇ ⋅ T u term; the problem is analogous to that associated with the ∇ ⋅ F i SF term in Equation (48), and so is the solution:

CONCLUSION
Prognostic equations for turbulent second-moment quantities have been derived in the multifluid framework. The equations presented allow for compressibility of the fluid (although compressibility effects are not expected to be large in the flow regimes of interest) and for the use of a filter of finite width that does not necessarily satisfy the Reynolds rules. The resulting equations contain terms that are analogues of all the terms that arise in the single-fluid case. In addition, they contain terms accounting for certain subfilter-scale pressure fluctuations that arise because of the use of a single filter-scale pressure for all fluid types. They also contain terms accounting for the possible relabelling of fluid parcels to represent entrainment and detrainment.
Following the work of Mellor and Yamada (1982) and others, the equations have been closed by modelling the terms that are not predicted, and a hierarchy of simplifications constructed by neglecting transience and transport terms in certain second-moment equations. In the single-fluid case it does not matter whether we work with flux form or advective form equations for the second moments when neglecting transience and transport. However, for the multifluid case we must work with the advective form to avoid picking up a spurious contribution from relabelling terms.
We have shown that, for consistency in the multifluid case, subfilter-scale fluxes in the presence of gradients in volume fraction must be accompanied by certain relabelling terms, which represent the effect of transport by those subfilter-scale fluxes across the boundary between fluid types.
The work so far thus establishes a consistent mathematical framework. In order for this framework to be useful for modelling, some further steps are needed.
• A variety of methods have been proposed in the literature, both prognostic and diagnostic, for obtaining a turbulent length-scale for the single-fluid case. The work of Ito et al. (2015) indicates that the turbulent length-scale must be modified in the grey zone, where the filter scale becomes comparable to or smaller than the scale of the largest turbulent structures. It seems very likely that the turbulent length-scale must also be modified in the multifluid case if the scale on which the fluid labels change is comparable with the scale of the largest turbulent structures. A parameterization for a turbulent length-scale that should be valid across a range of boundary-layer and convective-flow regimes is proposed by Lopez-Gomez et al. (2020) for their extended EDMF model.
• Suitable forms of the relabelling terms for second-moment quantities need to be formulated. A reasonable starting point would be to take the covariance in the relabelled air to equal the covariance in the "upstream" fluid:Ĉ i ≈ C for any second-moment property of the relabelled fluid. However, it is known that such a simple approach is not optimal for first-moment quantities, and improvements should be possible for second moments too.
• Correlations with water and entropy source terms are likely to be small in many situations. However, this might not be the case in precipitating cumulus clouds, and some quantitative estimates are needed.
• Some estimates for values of the pressure correlation terms denoted  (Equation B16) are needed. In particular, it would be useful to know whether  k i ≈ u i ⋅  i would be an adequate approximation in the TKE budget. It would be convenient if it could be assumed that  a i ≈ 0 for a = and a = q.
Future work should address these questions and others through a combination of quantitative analysis of high-resolution reference simulations, theoretical insights, and evaluation of models based on the framework presented here.
Part 2 of this set of articles presents the formulation of a two-fluid single-column model in which the ideas discussed here are applied to predict or diagnose a set of turbulent second-order moments; these second-order moments couple to the first-order moments through entrainment and detrainment, subfilter-scale fluxes, and the effects of subfilter-scale condensation on buoyancy.
Part 3 presents results from the single-column model for two shallow cumulus cases.

APPENDIX A. MULTIFLUID DECOMPOSI-TION OF MOMENTS
In this Appendix we introduce some of the notation used in the main text. We give expressions for the multifluid analogues of first, second, and third-order moments in compressible flow and for a general filter, and show how the single-fluid moments are given by a sum of contributions from each fluid in the multifluid case.

A.1 First-order moments
Let an overbar indicate a spatial filter of finite width-see Thuburn et al. (2018) for discussion-and let f * be the Favre, that is, density-weighted, average of some quantity f : it immediately follows from the definition of f i that Thus, the first-order moment f * is made up of contributions from each of the fluid types.

A.2 Second-order moments
The Favre average of fg is given by Defining the subfilter-scale contribution to the second moment (i.e. the generalized central second moment, Germano, 1992) by gives a decomposition into filter-scale and subfilter-scale contributions: In the multifluid case, where (fg) i is defined by i i (fg) i = I i fg. Now define the subfilter-scale contribution to the second moment from fluid i: Then we have the multifluid decomposition If needed, the decomposition of the total covariance can be obtained by using Equations (A9) and (A3) to eliminate fg, f * , and g * from Equation (A6), giving in agreement with Cohen et al. (2020).

A.3 Third-order moments
The Favre average of fgh is given by Defining the subfilter-scale contribution to the third moment (i.e. the generalized central third moment, Germano, 1992) by gives the decomposition ) .
(A13) In the multifluid case, where (fgh) i is defined by i i (fgh) i = I i fgh. Now define the subfilter-scale contribution to the third moment from fluid i: This gives the multifluid decomposition ) . (A16)

APPENDIX B. DERIVATION OF PROGNOS -TIC EQUATIONS FOR SECOND-ORDER MOMENTS
This Appendix gives details of the derivation of the multifluid prognostic equations for second-order moments presented in Section 3. Consider two generic variables (scalars or components of vector) a and b satisfying Combine Equation (B1) with the continuity Equation (2) to obtain the corresponding advective form equations: Combine Equation (B2) using the product rule to obtain D Dt (ab) = bA + aB.
With the aid of Equation (2) again, Equation (B3) becomes Next define (ab) i by i i (ab) i = I i ab; it is the total contribution to ab from fluid i. A prognostic equation for (ab) i is obtained by applying the conditional filter to Equation (B4), giving where F The conditionally filtered forms of Equation (B1) are Converting Equation (B7) to advective form using the density equation for fluid i (Equation 12), combining using the product rule, and converting back to flux form gives a prognostic equation for a i b i , the resolved contribution to ab from fluid i: Finally, we require a prognostic equation for C ab i = (ab) i − a i b i . C ab i is the unresolved or subfilter-scale contribution to ab from fluid i. It is the conditional filtering analogue, for a general filter, of the second-order quantity In the notation of Appendix A, implying that T uab is the contribution from fluid i to the third-order subfilter-scale turbulent flux. Now use the density equation for fluid i (Equation 12) to convert Equation (B9) to advective form. Since the relabelling terms take a similar form in all of the second-moment equations, it is convenient to introduce a shorthand notation. If we defineĈ ab i =(ab) i −â ibi to be the covariance of the relabelled fluid, then the relabelling terms become Thus, our final prognostic equation for C ab i is Taking a = and b = q and including the molecular dissipation term (Appendix C) immediately gives the prognostic equation (Equation 21) for the subfilter-scale covariance of and q within fluid i. Taking a = b = gives the prognostic equation (Equation 22) for the subfilter-scale variance of within fluid i. The equation for the variance of any other scalar, such as specific humidity, is analogous.
In the derivation of analogous prognostic equations for subfilter-scale scalar fluxes or momentum fluxes, there appear to be several alternative ways to rearrange the terms involving pressure to obtain equations of the desired form. It is highly desirable that the subfilter-scale terms should be term-by-term Galilean-invariant (e.g., Germano, 1992). Also, we desire that a buoyancy flux term should appear in the kinetic energy equation, as it does in the Boussinesq case. We therefore treat terms involving pressure in an analogous way to the u ⋅ ∇p term in Aluie (2013). (For an alternative, see e.g., Adumitroaie et al., 1999, but note that they assume a filter that obeys the Reynolds rules and also their expressions do not simplify so neatly for a general filter.) In the single-fluid case, Aluie's treatment of the pressure terms would correspond to the following in a scalar flux equation: The last term on the right-hand side is essentially the subfilter-scale covariance of a with the buoyancy. We wish to generalize this treatment to the conditionally filtered case, but in doing so we wish to avoid terms involving gradients of I i or i , as such terms may become large yet not have a simple and useful interpretation. The fact that we have a single mean pressure field seen by all fluid types also introduces a little extra complexity; see below. We therefore rearrange the pressure terms, for example, as they appear in Equation (B23) below, as follows: Here,  a i is defined by  a i = a i  i − I i a∇p + i a∇p. It accounts for differences between averages over fluid i and averages over all fluid for certain terms involving pressure gradients. An overbar with superscipt i indicates a (non-mass-weighted) filtered value over fluid i, for example, i a i = I i a.
Conditionally filtering then gives The advective form of the momentum equation for fluid i, Equation (16) Taking a i times Equation (B20) plus u i times the advective form of Equation (B7) gives Adding u i a i times the density equation for fluid i (Equation 12) and rearranging gives the corresponding flux form equation: Taking Equation (B19) minus Equation (B22) and converting to advective form gives a prognostic equation for the subfilter-scale contribution to the flux of a in fluid i: By analogy with Equations (B11) and (B12), Then, using our proposed form of the pressure terms in Equation (B16), Equation (B23) becomes Taking a = gives the C u i equation (Equation 23). Finally, consider the subfilter-scale momentum fluxes. Here, care is needed to keep track of the (unwritten) indices on tensor quantities, particularly which indices are contracted when a divergence is taken.
Taking u times the momentum equation (Equation 5) plus Equation (5) times u (where "times" here means the exterior product) gives Converting to flux form with the aid of the mass continuity equation (Equation 2) gives and applying the conditional filter then gives where Taking u i times Equation (B20) plus Equation (B20) times u i gives Here and below, [ ] indicates the transpose of the immediately preceding term. Adding u i u i times the density equation for fluid i (Equation 12) and rearranging gives the corresponding flux form equation: Also, it is to be understood that the divergence operator is contracted with the F u i SF in both the u i F u i SF and F u i SF u i terms.
Taking Equation (B28) minus Equation (B31), converting to advective form, and including the molecular dissipation term  uu i (Appendix C) then gives a prognostic equation for the contribution to the subfilter-scale momentum flux from fluid i: Keeping careful track of indices, it can be seen that . Also, by analogy with Equation (B16), This expression is notable because it makes the buoyancy flux term (third on the right) manifest, and also because the terms on the right-hand side are individually Galilean-invariant. Thus, Equation (B32) becomes

APPENDIX C. EFFECT OF MOLECULAR DIF-FUSION
In the single-fluid case, molecular diffusion and viscosity give rise to dissipation terms for some of the second-order moments. We must check that no significant additional complications arise in the conditionally filtered case.
Here we illustrate what happens in the case of a scalar variance; other cases are analogous.
Consider an advected scalar a subject to molecular diffusion: Repeating the steps of the derivation of Appendix B while carrying along the molecular diffusion term leads to a prognostic equation for the contribution from fluid i to the unresolved variance: The molecular diffusivity gives rise to several terms. In the first line is an additional molecular-diffusive contribution to the flux of variance. The analogous terms in the single-fluid case are neglected by Mellor (1973); Mellor and Yamada (1974) and we will do the same here. The fourth line contains contributions corresponding to molecular diffusion of variance across the boundary of fluid i. We will assume that the boundary of fluid i is sufficiently simple for these contributions to remain negligible. (Analogous terms are already implicitly neglected in the first-moment equations.) Finally, the third line contains the molecular dissipation of variance: This term will not generally be negligible, since, over long times and integrated over space, molecular dissipation must balance the various source terms. The dissipation can become significant despite the smallness of , because ∇a can become large on the subfilter scale. See Aluie (2013) for a discussion of the compressible case. Because the disspation is significant, it must be modelled in the variance equations. Analogous molecular dissipation terms arise in the scalar flux equations. In that case it is common to assume that, at small scale, velocity gradients are uncorrelated with scalar gradients and therefore the molecular dissipation of scalar flux is negligible: Moreover, the absence of an isotropic first-order tensor means there can be no dissipation term of the expected form in terms of filter-scale quantities (Mellor, 1973). Finally, analogous molecular dissipation terms arise in the momentum flux equations. It is common to assume that small-scale turbulence is isotropic, so that only TKE is dissipated:

APPENDIX D. EFFECT OF CONDENSATION
To account for the effects of subfilter-scale variability in and q on condensation, only very minor modifications are needed to the method of Mellor (1977b) to allow for compressibility. We therefore omit most of the details of the derivation. The results are assumed to apply separately in each fluid; the subscript i is therefore suppressed in many of the expressions below, except where needed to indicate a filter-scale mean.
The joint distribution G( , q) of entropy and total water, as sampled by the conditional filter at some location in space, is assumed Gaussian, with C i and C Source term for indicator function  (12)