Unified Entrainment and Detrainment Closures for Extended Eddy‐Diffusivity Mass‐Flux Schemes

Abstract We demonstrate that an extended eddy‐diffusivity mass‐flux (EDMF) scheme can be used as a unified parameterization of subgrid‐scale turbulence and convection across a range of dynamical regimes, from dry convective boundary layers, through shallow convection, to deep convection. Central to achieving this unified representation of subgrid‐scale motions are entrainment and detrainment closures. We model entrainment and detrainment rates as a combination of turbulent and dynamical processes. Turbulent entrainment/detrainment is represented as downgradient diffusion between plumes and their environment. Dynamical entrainment/detrainment is proportional to a ratio of a relative buoyancy of a plume and a vertical velocity scale, that is modulated by heuristic nondimensional functions which represent their relative magnitudes and the enhanced detrainment due to evaporation from clouds in drier environment. We first evaluate the closures off‐line against entrainment and detrainment rates diagnosed from large eddy simulations (LESs) in which tracers are used to identify plumes, their turbulent environment, and mass and tracer exchanges between them. The LES are of canonical test cases of a dry convective boundary layer, shallow convection, and deep convection, thus spanning a broad rangeof regimes. We then compare the LES with the full EDMF scheme, including the new closures, in a single‐column model (SCM). The results show good agreement between the SCM and LES in quantities that are key for climate models, including thermodynamic profiles, cloud liquid water profiles, and profiles of higher moments of turbulent statistics. The SCM also captures well the diurnal cycle of convection and the onset of precipitation.


Introduction
Turbulence and convection play an important role in the climate system. They transport energy, moisture, and momentum vertically, thereby controlling the formation of clouds and, especially in the tropics, the thermal stratification of the atmosphere. They occur on a wide range of scales, from motions on scales of meters to tens of meters in stable boundary layers and near the trade inversion, to motions on scales of kilometers in deep convection. General circulation models (GCMs), with horizontal resolutions approaching tens of kilometers, are unable to resolve this spectrum of motions. Turbulence and convection will remain unresolvable in GCMs for the foreseeable future , although some deep convective motions, on scales of kilometers to tens of kilometers, are beginning to be resolved in short-term global simulations (Kajikawa et al., 2016;Stevens et al., 2019).
Unable to resolve turbulence and convection explicitly, GCMs rely on parameterization schemes to represent subgrid-scale (SGS) motions. Typically, GCMs have several distinct parameterization schemes for representing, for example, boundary layer turbulence, stratocumulus clouds, shallow convection, and deep convection. The different parameterization schemes interact via trigger functions with discontinuous behavior in parameter space, even though in reality the flow regimes they represent lie on a continuous spectrum (Xie et al., 2019). This fragmentary representation of SGS motion by multiple schemes leads to a proliferation of adjustable parameters, including parametric triggering functions that switch between schemes. Moreover, most existing parameterizations rely on statistical equilibrium assumptions between the SGS motions and the resolved scales. These assumptions become invalid as model resolution increases and, for example, some aspects of deep convection begin to be explicitly resolved (Dirmeyer et al., 2012;Gao et al., 2017). It is widely recognized that these issues make model calibration challenging and compromise our ability to make reliable climate predictions (Hourdin et al., 2017;Schmidt et al., 2017;Schneider et al., 2017).
Many known biases in climate models and uncertainties in climate predictions are attributed to difficulties in representing SGS turbulence and convection. For example, biases in the diurnal cycle and the continental near-surface temperature, especially in polar regions, have been traced to inadequacies in turbulence parameterizations for stable boundary layers (Holtslag et al., 2013). Across climate models, biases in how tropical cloud cover covaries with temperature and other environmental factors on seasonal and interannual time scales are correlated with the equilibrium climate sensitivity, thus revealing the important role the representation of tropical low clouds plays in uncertainties in climate predictions (Bony & Dufresne, 2005;Caldwell et al., 2018;Ceppi et al., 2017;Cesana et al., 2018;Dong et al., 2019;Lin et al., 2014;Nam et al., 2012;Schneider et al., 2019;Teixeira et al., 2011). Differences in moisture export from the mixed layer to the free troposphere by cumulus convection lead to differences in the width and strength of the ascending branch of the Hadley circulation (Neggers et al., 2007). And biases in the structure of the South Pacific Convergence Zone have been traced to biases in the intensity of deep convective updrafts (Hirota et al., 2014). It is evident from these few examples that progress in the representation of SGS turbulence and convection is crucial for progress in climate modeling and prediction. At the same time, it is desirable to unify the representation of SGS motions in one continuous parameterization scheme, to reduce the number of adjustable parameters and obtain a scheme that more faithfully represents the underlying continuum of physical processes.
Different approaches for a systematic coarse graining of the equations of motion, leading to a unified parameterization, have been proposed (de Rooy & Siebesma, 2010;Han & Bretherton, 2019;Lappen & Randall, 2001a;Park, 2014aPark, , 2014bRio et al., 2019;Suselj et al., 2019b;Tan et al., 2018;Thuburn et al., 2018;Yano, 2014). They typically entail a conditional averaging (or filtering) of the governing equations over severalsubdomains (Weller & McIntyre, 2019), or an assumed probability density function (PDF) ansatz for dynamical variables and generation of moment equations from the ansatz Lappen & Randall, 2001a;Larson & Golaz, 2005;Larson et al., 2012). For example, conditional averaging can lead to a partitioning of a GCM grid box into subdomains representing coherent ascending and descending plumes, or drafts, and a more isotropically turbulent environment. Unclosed terms arise that, for example, to representinteractions among subdomains through entrainment and detrainment. Such unclosed terms need to be specified through closure assumptions (de Rooy et al., 2013). Or, if moment equations are generated through an assumed PDF ansatz for dynamical and thermodynamic variables, unclosed interactions among moments and dissipation terms need to be specified through closure assumptions Lappen & Randall, 2001b). Our goal in this paper is to develop a unified set of closures that work across the range of turbulent and convective motions, within one specific type of parameterization scheme known as an eddy-diffusivity mass-flux (EDMF) scheme (Siebesma & Teixeira, 2000;Siebesma et al., 2007;Wu et al., 2020).
We build on the extended EDMF scheme of Tan et al. (2018), which extends the original EDMF scheme of Siebesma and Teixeira (2000) by retaining explicit time dependence (SGS memory) and treating subdomain second-moment equations consistently, so that, for example, energy exchange between plumes and their environment obeys conservation requirements. The explicit SGS memory avoids any statistical equilibrium assumption. This is a necessary ingredient for the scheme to become scale aware and be able to operate in the convective gray zone, where deep convective motions begin to become resolved.
In this and the companion paper Lopez-Gomez et al. (2020) we present a set of unified closures that allow the extended EDMF parameterization to simulate stable boundary layers, dry convective boundary layers, stratocumulus-topped boundary layers, shallow convection, and deep convection, all within a scheme with unified closures and a single set of parameters. This paper focuses on unified entrainment and detrainment closures that are essential for convective regime, and Lopez-Gomez et al. (2020) present a closure for turbulent mixing. To demonstrate the viability of our approach, we compare the resulting parameterization scheme against large-eddy simulations (LESs) of several canonical test cases for different dynamical regimes. This paper is organized as follows. In section 2, we present the general structure of the extended EDMF scheme, including the subdomain decomposition and the prognostic equations for subdomain moments. Section 3 introduces the entrainment and detrainment closures that are key for the scheme to work across different dynamical regimes. Section 4 describes the numerical implementation of this scheme in a single-column model (SCM). In section 5, we describe the LES used in this study and how we compare terms in the EDMF scheme against statistics derived from the LES. Section 6 compares results from the EDMF scheme against LES of canonical test cases of dry convective boundary layers, shallow, and deep convection. Section 7 summarizes and discusses the main findings.

Equations of Motion
The extended EDMF scheme is derived from the compressible equations of motion of the host model. As thermodynamic variables, we choose the liquid-ice potential temperature θ l and the total water specific humidity q t , but these choices can easily be modified and harmonized with the thermodynamic variables of the host model in which the scheme is implemented. The unfiltered governing equations are as follows: In the momentum equation, to improve numerical stability, we have removed a reference pressure profile p h (z) in hydrostatic balance with a density ρ h (z): where g is the gravitational acceleration. Therefore, the perturbation pressure appear in the momentum equations in place of the full pressure p and gravitational acceleration g. Otherwise, the notation is standard: ρ is density, q t is the total water specific humidity, T v is the virtual temperature, R d is the gas constant for dry air, and is the liquid-ice potential temperature, with liquid and ice specific humidities q l and q i and reference surface pressure p s ¼ 10 5 Pa. In a common approximation that can easily be relaxed, we take the isobaric specific heat capacity of moist air c p to be constant and, consistent with Kirchhoff's law, the latent heat of vaporization L v to be a linear function of temperature (Romps, 2008). The temperature T is obtained from the thermodynamic variables θ l , ρ, and q t by a saturation adjustment procedure, and the virtual temperature T v is computed from the temperature T and the specific humidities (Pressel et al., 2015). The horizontal velocity vector is u h , and w is the vertical velocity component; ∇ h is the horizontal nabla operator. The symbol S stands for sources and sinks. For the velocities, the sources S uh and S w include the molecular viscous stress and Coriolis forces, and for thermodynamic variables, the sources S θl and S q t represent sources from molecular diffusivity, microphysics, and radiation.
When implemented in a GCM, the host model solves for the grid-averaged form of Equations 1-(6). In the averaged equations, SGS fluxes arise from the application of Reynolds averaging to quadratic and higher-order terms. As is common, we make the boundary layer approximation and focus on the vertical SGS fluxes, neglecting horizontal SGS fluxes. The role of the parameterization in the host model is to predict these vertical SGS fluxes, in addition to cloud properties that are used by radiation and microphysics schemes. In the next section, a decomposition of grid boxes into subdomains expresses the vertical SGS fluxes as a sum of turbulent fluxes in the environment (ED) and convective mass fluxes in plumes (MF).
To compute the MF component of the fluxes, the EDMF scheme solves for first moments of the host model's prognostic variables (w, θ l , q t ) in each of its subdomains, as well as for the area fraction of the subdomains.
To compute the ED component, the EDMF scheme solves additionally for the turbulence kinetic energy in the environment. Finally, to compute cloud properties by sampling from implied SGS distributions of thermodynamic variables, the EDMF scheme also solves for variances and covariance of θ l and q t in the environment. A summary of the prognostic and diagnostic variables in the scheme is given in Table 1.

Domain Decomposition and Subdomain Moments
The extended EDMF scheme is derived from the equations of motion by decomposing the host model grid box into subdomains and averaging the equations over each subdomain volume. We denote by ⟨ϕ⟩ the average of a scalar ϕ over the host model grid box, with ϕ * ¼ ϕ − ⟨ϕ⟩ denoting fluctuations about the grid mean.
Note. In the right two columns, "upd," "env," and "gm" stand for updrafts, environment, and grid mean, respectively, and these indicate whether a variable is prognostic or diagnostic in that model subdomain.
Similarly, ϕ i is the average of ϕ over the ith subdomain, and ϕ ′ i ¼ ϕ − ϕ i is the fluctuation about the mean of subdomain i. The difference between the subdomain mean and grid mean then becomes ϕ i * ¼ ϕ i − ⟨ϕ⟩. Common terminology assigns an area fraction a i ¼ A i /A T to each subdomain, where A i is the horizontal area of the ith subdomain and A T is the horizontal area of the grid box. This a i is more precisely a volume fraction, since A i is the vertically averaged horizontal area of the ith subdomain within the grid box. We retain here the terminology using subdomain area fractions, which reflect the subdomain volume fractions, consistent with previous works (Siebesma et al., 2007). (11) Equations (8) and (9) are self-evident; the derivation of (10) and (11) from (8) and (9) is given in Appendix A. Equation 10 with ϕ ¼ w is the vertical SGS flux of a scalar ψ, which is one of the key predictands of any parameterization scheme: The divergence of this flux appears as a source in the equations for the resolved scales of the host model. The decomposition in (9)-(11) only applies in general if ð · Þ is a Favre average-an average weighted by the density that appears in the continuity equation. However, in the EDMF scheme we describe in what follows, we make the approximation of ignoring density variations across subdomains (except in buoyancy terms), so that Favre and volume averages coincide within a grid box.
The central assumption in EDMF schemes is that within-subdomain covariances such as ϕ ′ i ψ ′ i and higher moments are neglected in all subdomains except one distinguished subdomain, the environment, denoted by index i ¼ 0. In the environment, covariances ϕ ′ 0 ψ ′ 0 are retained, and third moments such as w ′ 0 ϕ ′ 0 ψ ′ 0 , which appear in second-moment equations, are modeled with closures. The intuition underlying this assumption is that the flow domain is subdivided into an isotropically turbulent environment (i ¼ 0) and into coherent structures, identified with plumes (i ≥ 1). The environment can have substantial within-subdomain covariances, whereas the plumes are taken to have comparatively little variance within them. Variance within plumes can be represented by having an ensemble of plumes with different mean values (Neggers, 2012;Neggers et al., 2002;Sušelj et al., 2012). For the case of only two subdomains, an updraft (i ¼ 1) and its environment (i ¼ 0), the second-moment Equation 10 then simplifies to where the approximation in the second line reflects the EDMF assumption of neglecting within-plume covariances. The first line states that the covariance on the grid scale can be decomposed into the sum of the covariances within subdomains and the covariance among subdomain means, as in the analysis of variance (ANOVA) from statistics (Mardia et al., 1979). The second line reflects the EDMF approximation to only retain the covariances in the environment. The first term on the right-hand side is closed by a downgradient eddy diffusion (ED) closure and the second term is represented by a mass flux (MF) closure, whence EDMF derives its name (Siebesma & Teixeira, 2000). Whenever ϕ and ψ are both thermodynamic prognostic variables, the within-environment covariance ϕ ′ 0 ψ ′ 0 is solved prognostically. Under the EDMF assumption, the third-moment Equation (11) for two subdomains, written for a single scalar, simplifies to That is, third moments (i.e., skewness) on the grid scale are represented through covariances within the environment and through variations among means across subdomains with differing area fractions.

EDMF Assumptions
The extended EDMF scheme is obtained by applying this decomposition of grid-scale variations to the equations of motion (1)-(6), making the following additional assumptions: 1. We make the boundary layer approximation for subgrid scales, meaning that we assume vertical derivatives to be much larger than horizontal derivatives. This in particular means that the diffusive closure for fluxes in the environment only involves vertical gradients, where K ϕ,i is the eddy diffusivity (to be specified) for scalar ϕ in subdomain i. Consistent with the EDMF assumptions, we assume K ϕ,i ¼ 0 for i ≠ 0. 2. We use the same, grid-mean density ⟨ρ⟩ in all subdomains except in the buoyancy term. This amounts to making an anelastic approximation on the subgrid scale, to suppress additional acoustic modes that would otherwise arise through the domain decomposition. For notational simplicity, we use ρ rather than ⟨ρ⟩ for the grid-mean density in what follows, and ρ i for the subdomain density that appears only in the buoyancy term: The grid-mean density ρ appears in the denominator, playing the role of the reference density in the anelastic approximation. The area fraction-weighted sum of the subdomain buoyancies is the grid-mean buoyancy, ensuring consistency of this decomposition: 3. We take the subdomain horizontal velocities to be equal to their grid-mean values, This simplification is commonly made in parameterizations for climate models (Larson et al., 2019). It eliminates mass-flux contributions to the SGS vertical flux of horizontal momentum.

EDMF Equations
The full derivation of the subdomain-mean and covariance equations from (1)-(6) is given in Appendix B. The derivation largely follows Tan et al. (2018), except for a distinction between dynamical and turbulent entrainment and detrainment following de Rooy and Siebesma (2010). The resulting extended EDMF equation for the subdomain area fraction is the equation for the subdomain-mean vertical momentum is 10.1029/2020MS002162

Journal of Advances in Modeling Earth Systems
and the equation for the subdomain-mean of a thermodynamic scalar ϕ is The dynamical entrainment rate from subdomain j into subdomain i is E ij , and the detrainment rate from subdomain i into subdomain j is Δ ij . In addition to dynamical entrainment, there is turbulent entrainment from subdomain j into subdomain i, with rate Ê ij . Turbulent entrainment differentially entrains tracers but not mass (see Appendix B).
The pressure and buoyancy terms in the vertical momentum Equation (19) are written as the sum of their grid-mean value and perturbations from their grid-mean value. These perturbations vanish when summed over all subdomains because ∑ i a i ϕ i * ¼ 0; hence, the grid-mean values of the pressure and buoyancy terms are recovered upon summing over subdomains. Following Pauluis (2008), the pressure gradient term in (19) is written with 1/ρ inside the gradient to ensure energy conservation in our SGS anelastic approximation; see Appendix C for details. The subdomain density ρ i that is essential for the subdomain buoyancy is computed from the subdomain virtual temperature T v; i using the ideal gas law with the grid-mean pressure ⟨p⟩: In analogy with the anelastic approximation of Pauluis (2008), this formulation of the ideal gas law ensures that ∑ i a i ρ i T v; i ¼ ρ⟨T v ⟩, while accounting for subdomain virtual temperature effects that play a key role in the buoyancy of updrafts in shallow convection.
The scalar Equation (20) is applied to any thermodynamic variable, with its corresponding subdomain-averaged source S ϕ; i on the right-hand side. The terms on the left-hand side represent the explicit time tendencies and fluxes of the subdomain-means, which can be viewed as forming part of the dynamical core of the host model. The terms on the right-hand side are sources and sinks that require closure. The covariance equation for thermodynamic scalars (i.e., when ϕ,ψ ∈ [θ l , q t ]) in the environment becomes 10.1029/2020MS002162

Journal of Advances in Modeling Earth Systems
Consistently with the EDMF assumption, we have assumed here that ϕ ′ i ψ ′ i ¼ 0 for i > 0. Covariance equations of this form are used for the thermodynamic variances θ ′2 l; 0 and q ′2 t; 0 and for the covariance θ ′ l; 0 q ′ t; 0 , which are needed in microphysics parameterizations. Note that some of the entrainment and detrainment terms are cross-subdomain counterparts of the vertical gradient terms. For example, the "dynamical entrainment," "turbulent entrainment," and "turbulent entrainment production" are the cross-subdomain counterparts of the "vertical transport," "turbulent transport," and "turbulent production," respectively. The "dynamical entrainment flux" lacks any vertical counterpart. This term arises as a flux across a variable boundary in the conditional averaging process.
The subdomain turbulence kinetic energy (TKE) is defined as i Þ, and the TKE equation for the environment is written as

Effect on Grid Mean and Constraints on Entrainment/Detrainment
The conservation of mass and scalars in the host model grid box requires that by summing the EDMF equations over all subdomains, the equations for the grid-mean variables are recovered. The horizontal flux divergence terms that are included in the EDMF equations, ∇ h · ðρa i ⟨u h ⟩ϕ i Þ , represent the fluxes across the boundaries of the host model grid (see Appendix B) and, when summed over all subdomains, recover their grid-mean counterpart. Additionally, mass conservation requires that between two subdomains i and j, the entrainment and detrainment rates satisfy For entrainment and detrainment of subdomain-mean properties, scalar conservation further requires that so that when summing over two interacting subdomains, the entrainment and detrainment terms cancel out. Similarly, scalar conservation requires symmetry, Ê ij ¼ Ê ji , for turbulent entrainment.
Taking these requirements into account, a summation of Equation (20) over all subdomains yields the grid-mean scalar equation 10.1029/2020MS002162

Journal of Advances in Modeling Earth Systems
This is the form of the equation solved by the dynamical core of the host model. Using the covariance decomposition (10), the SGS flux in (25) is written as the sum of the eddy diffusivity and mass flux components: This illustrates the coupling between the dynamical core equations and the EDMF scheme. Similarly, the grid covariance equation follows by using the subdomain continuity equation (18), scalar-mean equation (20), and the scalar covariance equation (22) in the covariance decomposition (10), which yields Here, vertical SGS fluxes are decomposed according to Equation (26), and the turbulent transport term is decomposed according to Equation (11). In general, Equation (27) does not need to be solved by the host model. However, the consistency of the summation over subdomains to produce it ensures that the second moments are conserved within the EDMF scheme.
The subdomain equations in the EDMF scheme require closures for dynamical entrainment and detrainment, turbulent entrainment, perturbation pressure, eddy diffusivity, for the various sources, and for covariance dissipation. The following section focuses on closures for dynamical and turbulent entrainment and detrainment. The perturbation pressure closure is given by the sum of a virtual mass effect, momentum convergence, and pressure drag, see equation (11)

Closures
Entrainment and detrainment closures are a topic of extensive research (de Rooy et al., 2013). Following de Rooy and Siebesma (2010), we distinguish dynamical and turbulent entrainment and detrainment components. Turbulent entrainment is typically represented by a diffusive horizontal flux, while diverse closures for dynamical entrainment and detrainment are in use. It is common to write the dynamical entrainment and detrainment rates as a product of the vertical mass flux ρa i w i and fractional entrainment/detrainment rates ϵ ij and δ ij and Closures are then derived for the fractional rates ϵ ij and δ ij per unit length (they have units of 1/length).
Various functional forms for the fractional rates ϵ ij and δ ij have been proposed in the literature. For example, • Based on experiments on dry thermals, Morton et al. (1956) suggested ϵ ij to be inversely proportional to the updraft radius. This relation has been used in several closures (Bretherton et al., 2004;Kain & Fritsch, 1990).  Tian and Kuang (2016). In the steady equations, this entrainment functional form also ensures that the mass flux and the vertical velocity simultaneously go to zero at the top of updrafts; see Appendix E and Romps (2016). Alternative derivations of this functional form are based on a balance of sources and sinks of total kinetic energy in updrafts (Savre & Herzog, 2019), or on the dynamics of dry thermals (McKim et al., 2020).
Similar closures are often used for both entrainment ϵ ij and detrainment δ ij . Enhanced detrainment can occur in cloudy conditions: When the evaporation of cloud condensate after mixing with drier environmental air produces a buoyancy sink for an updraft, negatively buoyant air can detrain rapidly from the updraft (Kain & Fritsch, 1990;Raymond & Blyth, 1986). Various approaches for representing this enhanced detrainment owing to "buoyancy sorting" have been used, ranging from adding a constant background detrainment rate (Siebesma & Cuijpers, 1995;Tan et al., 2018), over explicitly modeling buoyancies of mixtures of cloudy and environmental air (Bretherton et al., 2004;Kain & Fritsch, 1990), to enhancing detrainment by functions of updraft-environment relative humidity differences (Bechtold et al., 2008(Bechtold et al., , 2014Böing et al., 2012;Savre & Herzog, 2019).
Here we combine insights from several of these studies into a new closure for entrainment and detrainment.

Dynamical Entrainment and Detrainment
We propose closures for dynamical entrainment and detrainment that are in principle applicable to many interacting subdomains (e.g., multiple updrafts, or updrafts and downdrafts). Our point of departure are dry entrainment and detrainment rates which are symmetric for upward and downward motions. To those we then add the contribution of evaporation, which is asymmetric between upward and downward motions. We first write our closures for the rates E ij and Δ ij , which facilitate ensuring mass and scalar conservation. In the end, we give the corresponding formulations in terms of the fractional rates ϵ ij and δ ij .

General Form of Entrainment and Detrainment Rates
The rates E ij and Δ ij have units of density divided by time and hence depend on a flow-dependent time scale, as well as on functions of nondimensional groups in the problem. Following Gregory (2001) This vertical velocity scale is taken to be representative of the vertical velocity difference across the updraft boundary, which we approximate as the difference between the subdomain means in convective conditions. In cases of strong environmental turbulence and weak updraft velocities, the environmental turbulent velocity scale ē 1=2 0 is a better representation of this velocity difference. This is the case in conditions of weak surface heating, such as those encountered in stratocumulus-topped boundary layers (Lopez-Gomez et al., 2020). Thus, the velocity scale w is taken as the maximum of the previously described scales. Considerations of symmetry and mass and tracer conservation lead to the inverse time scale Here, λ ij ¼ λ ji , c λ is a nondimensional fitting parameter, and s min is the smooth minimum function defined in Lopez-Gomez et al. (2020). The smooth minimum function ensures that the strongest characteristic velocity defines the entrainment rate. The inverse time scale λ ij depends on the buoyancy difference b i − b j between subdomains i and j, as is physical. Similarly, λ ij depends only on the mean vertical velocity difference w i − w j , as is required by Galilean invariance. In terms of this inverse time scale, the entrainment and detrainment rates are then written as and Mass and tracer conservation demand that E ij ¼ Δ ji (see Equation 24). This is satisfied by this formulation: The inverse time scale λ ij is symmetric under reversal of the i and j indices by construction. Conservation

Journal of Advances in Modeling Earth Systems
constraints are satisfied by the choice of the, as yet unspecified, nondimensional functions D ij and M ji in the entrainment rate (31) and, with inverted indices, D ji and M ij in the detrainment rate (32). The coefficients c ϵ and c δ are nondimensional fitting parameters. The functions D ij and M ij in principle can depend on all nondimensional groups of the problem. Once sufficient data are available, be they from high-resolution simulations or observations, they can be learned from data.
To demonstrate the viability of the EDMF closure, we use physically motivated and relatively simple functions for D ij and M ij .

Function D ij
We use the function D ij to estimate the relative magnitudes of entrainment and detrainment for a subdomain i in dry convection, in which case the subdomain buoyancy is linearly mixed. We consider the buoyancy b mix of a mixture, composed of a fraction χ i of air from subdomain i, and a fraction χ j of air from subdomain j (with χ i + χ j ¼ 1). We define an inverse time scale based on the mixture buoyancy as is the area-weighted mean buoyancy of subdomains i and j, such that a i + a j ¼ 1 implies b ij ¼ ⟨b⟩. (Note that we are assuming dry conditions here, so buoyancy averages linearly.) Here μ ij ¼ −μ ji , and its sign reflects the correlation between the sign of the velocity difference w i − w j and the sign of the mixture buoyancy b mix relative to the mean buoyancy b ij . The mixture buoyancy is defined as so that the buoyancy difference in (33) becomes which follows by using Thus, we assumed that the more rapidly rising subdomain entrains air if the mixture buoyancy is positive relative to the mean of the two interacting subdomains, and vice versa. This means that we expect entrainment from subdomain j into i if μ ij > 0, and we expect detrainment otherwise. This could be modeled by choosing D ij ¼ maxðμ ij ; 0Þ. However, we find that using a smooth sigmoidal function, between 0 and 1, improves our results, so we define Here, μ 0 is an inverse time scale, a fitting parameter that controls the smoothness of the sigmoidal function. We estimate μ 0 ¼ 4 × 10 −4 s −1 from examining various LES test cases. The fact that this is a dimensional coefficient is a shortcoming of the current model; we aim to replace this by a function of grid-mean quantities in future work. The fraction of air in the mixture, χ i , is typically taken from an assumed probability distribution (Bretherton et al., 2004;Kain, 2004). Here we choose a constant χ i for updrafts interacting with their environment, based on a heuristic assumption of an elliptical updraft in a surrounding mixing shell. If the mixing eddies at the updraft edge have similar radial extent in the updraft and in the shell, it implies that χ i is proportional to the ratio between the updraft area and the combined updraft and shell area; that is, χ i ¼ 0.25. For interactions between two updrafts (or downdrafts), the corresponding choice would be

Function M ij
In moist conditions, the function M ji represents the enhancement of detrainment from the rising subdomain i (and entrainment into the sinking subdomain j) by evaporation of liquid water when i is cloudy (saturated).
In dry conditions, we expect M ji ¼ M ij ¼ 0. Similar to Savre and Herzog (2019), the evaporative potential of the drier subdomain j is approximated here by an ad hoc function of the difference between the relative humidities RH i and RH j of the subdomains, conditioned on the saturation of subdomain i: Here, β is a nondimensional parameter that controls the magnitude of the evaporative potential for a given relative humidity difference. With this closure, a saturated updraft i detrains when the environment j ¼ 0 is subsaturated, and the detrainment rate increases with increasing subsaturation of the environment.

Fractional Entrainment and Detrainment Rates
Given the relationships (28) and (29) between the entrainment rates E ij and D ij and their fractional counterparts ϵ ij and δ ij , the fractional rates are and The relationship E ij ¼ Δ ji required for scalar and mass conservation in terms of the fractional rates implies The difference between the fractional rates, which is the source of ρa i , is The function D ij − D ji appearing here is a sigmoidal function between −1 and 1.
For the situation where entrainment is only considered between an updraft i and the environment j ¼ 0, and if the environmental mean vertical velocity w 0 and turbulent kinetic energy ē 0 are neglected, this closure reduces to a closure of the form b i =w i 2. It is heuristically modulated by the nondimensional functions D ij and M ij , which approximate the relative magnitudes of entrainment and detrainment while accounting for enhanced detrainment owing to evaporation of condensate.

Turbulent Entrainment
We assume that turbulent entrainment takes place only between the plumes (updrafts and downdrafts) and their environment, where second moments are not neglected. Therefore, we assume it depends on the turbulent velocity scale of the environment, ffiffiffiffi ffi ē 0 p , and the radial scale of a plume R i . The turbulent entrainment rate is related to the flux across the subdomain boundary via where A sg and V i are the updraft's interface area and volume (see the derivation of (B10)). We assume here that the updraft is cylindrical with a circular cross section, so that the ratio between its interface area and its volume is A sg /V i ¼ 2/R i . Following de Rooy and Siebesma (2010), Asai and Kasahara (1967), and Kuo

Journal of Advances in Modeling Earth Systems
(1962) the outward pointing turbulent flux across the boundary of the ith updraft, d ϕ ′ u ′ r; n , is modeled by downgradient eddy diffusion HereK i0 is the entrainment eddy diffusivity between the environment and the ith subdomain. The cross-subdomain gradient is discretized using the difference in the mean values of the two interacting subdomains and the radial scale of the updraft R i . The latter is written in terms of updraft height H i and an aspect ratio γ as R i ¼ γH i . The updraft height H i is taken to be the maximal height at which a i > 0 in the previous time step, but at least 100 m to avoid division by zero in the initial stages of the simulation. For the entrainment eddy diffusivity, we assume the form where R i is used as a mixing length and c t is a nondimensional fitting parameter.
Combining Equations 42- (44), we obtain the turbulent entrainment rate where c γ ¼ c t /γ is a fitting parameter that combines c t and γ ( Table 2). The middle term in (45) shows that Ê ij ∝ 1=R i , in agreement with laboratory experiments of dry plumes (Morton et al., 1956;Turner, 1963). It is also useful to define a fractional counterpart for turbulent entrainment,

Numerical Implementation
The model equations and closures are implemented in the SCM used in Tan et al. (2018), where a detailed description of the implementation of the initial and boundary conditions is given. The model solves for first moments of the prognostic variables fa i ; w i ; θ l; i ; and q t; i g in updrafts using (18), (19), and (20), respectively, and for the grid mean variables {⟨θ l ⟩,⟨q t ⟩} using equations of the form of (25), in which prescribed large-scale tendencies are applied as sources.
We consider a single updraft and its turbulent environment. The mean environmental properties are computed diagnostically as the residual of updraft and grid-mean quantities using (8) and (9). Prognostic equations for the second moments (θ ′2 l; 0 , q ′2 t; 0 , θ ′ l; 0 q ′ t; 0 ē 0 ) in the environment are solved using (22) and (23). The grid-scale second moments are diagnosed from (10), using the EDMF assumption of neglecting second moments in the updraft. Grid-scale third moments are diagnosed using (11), neglecting third moments in all individual subdomains. Thus, from a probability density function perspective, we are using a closure model that assumes a Gaussian environment and a delta distribution updraft (Lappen & Randall, 2001a).
The parameters we use in the entrainment and detrainment closures are shown in Table 2. The parameters in this study and in Lopez-Gomez et al. (2020) were chosen sequentially: We first calibrated a subset of parameters associated with turbulent mixing based on stable boundary layer simulations (Lopez-Gomez et al., 2020). We then searched for a combination of parameters related to dry convection (c ϵ , c t , c λ ) so that the EDMF scheme captures the DCBL and the subcloud layer in moist convective cases. Finally, we optimized the moisture-related parameters (β,c δ ) based on the EDMF scheme's ability to capture cloud layer properties and the cloud top height. The initial conditions, surface fluxes, and large-scale forcing are case specific. They are taken from the papers describing the cases, are linearly interpolated to the model resolution, and are implemented identically in the SCM and LES.
The SCM implementation of the EDMF scheme makes several assumptions because the SCM does not solve for the density, pressure, or vertical velocity of the grid-mean. In the SCM, it is assumed that ⟨w⟩ ¼ 0 and ρ ¼ ρ h in the EDMF equations, and consequently that ρ ¼ ρ h in the denominators of the buoyancy definitions (15) and (16). Furthermore, the grid-mean anelastic approximation requires the use of the reference pressure (p h ) in the ideal gas law (21) for consistency (Pauluis, 2008). The SCM is therefore fully anelastic, in contrast to the SGS anelastic approximation described in Appendix C. Since ⟨w⟩ ¼ 0, the balance in the ⟨w⟩ equation is reduced to thus removing from the subdomain equations the dependence on the grid-mean pressure.
All SCM simulations use a uniform vertical resolution of 50 m, with results from a resolution sensitivity test at 100 and 150 m shown for the first three moments in the grid. Other implementation details, such as how cloud properties are computed via numerical quadrature over implied SGS distributions, are described in Lopez-Gomez et al. (2020).

LES and Diagnosis of EDMF Subdomains
To assess the performance of the extended EDMF scheme, we compared it with LES in four convective test cases. We use PyCLES (Pressel et al., 2015), an anelastic LES code with weighted essentially nonoscillatory (WENO) numerics. We use an implicit LES strategy, which uses the dissipation inherent to WENO schemes as the only subgrid-scale dissipation. Such an implicit LES has been shown to outperform explicit SGS closures in simulations of low clouds Schneider et al., 2019). We use passive tracers that decay in time to diagnose updrafts and their exchanges with the environment in the LES (see Appendix D).
Four standard convective test cases are considered here: dry convective boundary layer, maritime shallow convection, continental shallow convection, and continental deep convection.
1. The Dry Convective Boundary Layer test case (DCBL, blue lines in all figures) is based on Soares et al. (2004). In this case, convection develops through 8 hr from an initially neutral profile below 1,350 m (which is stable above it) with prescribed sensible and latent heat fluxes and negligible large-scale winds. We use an isotropic 25 m resolution in a 6.4 km × 6.4 km × 3.75 km domain. 2. The marine shallow convection test case is based on the Barbados Oceanographic and Meteorological Experiment (BOMEX, orange lines) described in Holland and Rasmusson (1973). In this case, large-scale subsidence drying and warming and fixed surface fluxes are prescribed, and subtropical shallow cumulus convection evolves over 6 hr, with a quasi-steady state maintained in the last 3 hr (Siebesma et al., 2003). We use an isotropic 40 m resolution in a 6.4 km × 6.4 km × 3 km domain.

The continental shallow convection test case is based on the Atmospheric Radiation Measurement
Program at the United States's Southern Great Plains (ARM-SGP, green lines) described in Brown et al. (2002). This case exhibits a diurnal cycle of the surface fluxes, with cumulus convection first developing and then decaying between 5:30 and 20:00 local time. We use 100 m × 100 m × 40 m resolution in a 25 km × 25 km × 4 km domain. The large surface fluxes of latent and sensible heat erode the initial inversion as convection penetrates into the free atmosphere (Brown et al., 2002). 4. The continental deep convection test case is based on the Large-scale Biosphere-Atmosphere experiment with data from the Tropical Rainfall Measurement Mission (TRMM-LBA, red lines) observed on 23 February 1999 in Brazil (Grabowski et al., 2006). In this case, prescribed time-varying surface fluxes and radiative cooling profiles force a diurnal cycle, during which shallow convection transitions into deep convection in the 6 hr between 7:30 and 13:30 local time. We use 200 m × 200 m × 50 m resolution in a 51.2 km × 51.2 km × 24 km domain. No subsidence drying or warming are prescribed in this case. In our simulations of the TRMM-LBA case, microphysical rain processes are modeled by a simple warm-rain cutoff scheme that removes liquid water once it is 2% supersaturated. This simple scheme is 10.1029/2020MS002162

Journal of Advances in Modeling Earth Systems
implemented in the LES for a direct comparison with the same simple microphysics scheme in SCM. In future work, we will implement a more realistic microphysics scheme.
The different cases span a wide range of conditions that allow us to examine the different components of the unified entrainment and detrainment formulation presented in section 3. The DCBL case allows us to examine the dry formulations for dynamic and turbulent entrainment irrespective of the moisture related detrainment. The differences in environmental humidity between the shallow and deep convection cases allows us to test the moisture-dependent detrainment closure. For instance, we found the bulk detrainment value used in previous parameterization evaluated with BOMEX (Siebesma & Cuijpers, 1995;Tan et al., 2018) to be excessive for TRMM-LBA.
The diagnosis of the direct estimates of entrainment and detrainment and comparison with the closures (39) and (40) relies on decaying tracers with a surface source, which uniquely identify each LES grid box as either updraft or environment. Here we use the tracer scheme described in Couvreux et al. (2010), which labels a grid cell as updraft if its vertical velocity, tracer concentration, and liquid water specific humidity (above cloud base) exceed given thresholds. The net of entrainment minus detrainment [right-hand side of (18)] is diagnosed using the area and vertical velocity of updrafts identified with the help of the tracer scheme.
Fractional entrainment is diagnosed based on an advective form of the scalar equation, see Equation D1. Further information on the diagnosis is found in Appendix D.

Results
A comparison of the closures for the fractional turbulent and dynamic entrainment and detrainment rates with direct estimates of these terms from LES is shown in Figure 1. In this comparison, the profiles of the

Journal of Advances in Modeling Earth Systems
EDMF closures are based on diagnosing all EDMF components (area fractions, first and second moments) from LES and using those in the EDMF closures described in section 3. The profiles of the closures for entrainment and detrainment are similar to the direct estimates from LES. The role of the environmental moisture deficit in enhancing detrainment in the cloud layer is consistent with the directly diagnosed detrainment in ARM-SGP, in which convection penetrates into a dry layer with RH ≈ 50%.
When implemented in the SCM, these closures perform in a similar manner (Figure 2). Dynamic entrainment prevails in the subcloud layer, while dynamic detrainment prevails in the cloud layer, owing to the large environmental moisture deficit. The value of ϵ − δ predicted by the closures in the EDMF scheme is in agreement with direct estimates of this value from LES (solid gray lines). Turbulent entrainment is about half the dynamic entrainment in the boundary layer and vanishes above it. A discrepancy between the SCM and LES is found between the entrainment and detrainment profiles for the DCBL case. The LES updrafts detrain from mid levels and upward, whereas the SCM updrafts detrain mostly at their tops. This could indicate of a downside of the current closure that uses the subdomain mean buoyancy and does not detrain from buoyant updrafts. A more sophisticated scheme, in which entrainment dependents on second moments, could improve the performance at the cost of computing second moment in all subdomains.
We now turn to compare the performance of the EDMF scheme with LES. First, second, and third moments of θ l and q t are compared in Figures 3, 4, and 5. These show overall good matches between the SCM and LES, with a few notable mismatches. For example, in first moments in the sub-cloud layer in the ARM-SGP case, at cloud top in the BOMEX case, and at the top of the DCBL; in second moments (⟨θ * 2 l ⟩) throughout the DCBL; and in the third moments at the overshoots. Moreover, mismatches in sign are seen for ⟨θ * 3 l ⟩ in SCM simulations of TRMM-LBA at mid levels, and for ⟨q * 3 t ⟩ at the top of the DCBL. The sensitivity test at The grid-mean SGS fluxes, whose divergence is a source in the host model equations, are shown in Figure 6. We find good agreement in the fluxes except for ⟨w * θ * l ⟩ in TRMM-LBA at midlevels, where the SCM shows a strongly positive flux while the LES has a negligible flux there. The ED and MF components of the SCM   Figure 3 but for the third moments ⟨θ * l θ * l θ * l ⟩ and ⟨q * t q * t q * t ⟩. The DCBL spike in the ⟨q * t q * t q * t ⟩ profile (blue) has an amplitude of −1.5 (g 3 /kg 3 ). The comparison of updraft and cloud properties in Figure 7 shows good agreement with LES above cloud base. Below cloud base and in the DCBL, large disagreements in the mass flux and updraft fractions are found. However, in the boundary layer, the diagnosis of updrafts in the LES can be misleading because lateral turbulent mixing makes the distinction between updrafts and their environment ambiguous. We did not attempt to implement a more sophisticated scheme, such as (Efstathiou et al., 2020) in this work. However, the key predictions of the EDMF scheme (the SGS vertical fluxes and the mean profiles on the host model grid) are in good agreement with the LES (Figure 6). This implies that the net of ED and MF effects in the SCM reproduces the well-mixed boundary layer, even though the decomposition into updrafts and environment may not be exact.
Diurnal cycles of shallow and deep convection are shown in Figure 8. The onset of convection in the SCM is found to be about half an hour delayed compared with the LES, while cloud top height is in good agreement between the models. In the decay stage in the ARM-SGP case, the cloud in the SCM shuts off abruptly, unlike the gradual decline in the LES. This may result from the EDMF assumption that neglects variance in the (single) updraft, which cannot cross cloud base when its buoyancy right below cloud base is too low. Good agreement is found in the liquid water path (LWP) between the SCM and the LES in both cases. In the TRMM-LBA case, this agreement includes the effect of precipitation on the column integrated q t . The precipitation sink is used to compute rain rates in the cutoff microphysics scheme as the vertically integrated amount of q t removed at a model time step per unit area. The EDMF rain rates peak at nearly twice their LES counterparts in the TRMM-LBA case (Figure 9). This overestimation is consistent with the overestimation of w upd (Figure 7). Tuning the maximum supersaturation in the cutoff microphysics could improve both the vertical velocity and the rain rates, although this was not explored here. The coarse-graining of the convective plumes into a single updraft in the EDMF scheme may indicate that a different supersaturation should be applied in the SCM compared with the LES. Mean profiles of cloud properties over the last 2 hr (Hours 9-11 in ARM-SGP). Top to bottom rows correspond to DCBL, BOMEX, ARM-SGP. and TRMM-LBA, with SCM following the color-coding in Figure 3 and corresponding LES in gray. Left to right columns correspond to updraft massflux, updraft fraction (dashed) and cloud fraction (solid), updraft vertical velocity, and liquid water specific humidity, respectively.

Discussions and Conclusions
We have presented entrainment and detrainment closures that allow the extended EDMF scheme to simulate boundary layer turbulence, shallow convection, and deep convection, all within a unified physical framework. The results demonstrate the potential of the extended EDMF scheme to serve as a unified parameterization for all SGS turbulent and convective motions in climate models (other SGS motions such as gravity waves require additional parameterizations). The choice of parameters used to produce these results is uniform across all cases, as well as across all cases shown in Lopez-Gomez et al. (2020). We view these results as a proof of concept, which we will improve further using automated model calibration techniques and a larger LES data set in the future.
The dynamic entrainment/detrainment closures are based on a combination of a b/w 2 scaling and physically motivated nondimensional functions, which can in principle be learnt from data. At the moment, these functions are based on arguments from buoyancy sorting and relative humidity differences between clouds and their environment. The addition of turbulent entrainment, which only affects scalars, allows us to regulate the mass flux by reducing the vertical velocity without increasing the area fraction below cloud base, where detrainment is negligible.
The extended EDMF scheme produces good agreement with LES in key properties needed for climate modeling. The successful simulation of high-order moments and vertical fluxes justifies the EDMF assumption of a negligible contribution from updraft covariance to the grid scale covariance. It would be straightforward to include multiple updrafts (Neggers, 2012;Neggers et al., 2002;Sušelj et al., 2012), which can further improve the results. Using multiple updrafts would also open up the opportunity to include stochastic components either in the updrafts' boundary conditions or in the entrainment and detrainment closures (Romps, 2016;Suselj et al., 2013Suselj et al., , 2014Suselj et al., , 2019a, with the nonlinearity of the model ensuring that the stochastic effect will not average out in the grid mean. Nonetheless, the use of multiple updrafts results in a higher computational overhead of the parameterization in climate simulations. This added cost may be ameliorated harnessing the power of parallel architectures. There is a growing interest in using artificial neural networks as SGS models for turbulence and convection (e.g., O'Gorman & Dwyer, 2018;Rasp et al., 2018). It is worth noting that the extended EDMF scheme with multiple updraft and downdraft has a network structure: The subdomains play the role of network nodes, which interact through sigmoidal activation functions (entrainment/detrainment). Each node has memory (explicitly time-dependent terms), somewhat akin to long short-term memory (LSTM) networks (Hochreiter & Schmidhuber, 1997). Unlike artificial neural networks whose architecture is not tailor-made for the physical problem at hand, the architecture of the extended EDMF scheme ensures physical realizability and conservation of energy. Like for neural networks, the activation functions and other parameters in the extended EDMF scheme can be learnt from data. Our results, which required adjustment of only a handful of parameters, show that only a small fraction of the data typically required to train neural networks is needed to calibrate the extended EDMF scheme.
The explicitly time-dependent nature of the extended EDMF scheme makes it well suited to operate across a wide range of GCM resolutions  and under time-varying large-scale conditions that may include diurnal cycles and variability on even shorter time scales (Tan et al., 2018).
Without loss of generality, the subdomain boundary ∂Ω i can be expressed as the union the part of the subdomain Ω i boundary that coincides with the grid box Ω T boundary. The domain and subdomain boundaries are related through ∑ i ∂Ω g i ¼ Ω T . The subgrid boundary ∂Ω sg i is a free-moving surface with velocity u b , while boundary ∂Ω g i is fixed. Using the Reynolds transport theorem for the transient term, the Gauss-Ostrogradsky theorem for the divergence, and rearranging the surface integrals yields where n is the outward pointing unit vector normal to the surface over which the integration is performed. The first term on the right-hand side is the flux out of subdomain Ω i into other subdomains within the same grid box, and the second term on the left-hand side is the flux out of subdomain Ω i to a neighboring grid box. The total grid-scale divergence equals the sum of fluxes from all subdomains across the grid box, where the commutativity of the gradient and the volume average is exact for uniform grids and results in a small error otherwise (Fureby & Tabor, 1997). Using the domain decomposition in (9), the leftmost term in (B3) can be written in terms of the sum of the subdomain-mean values, where V i is the volume of subdomain Ω i , and (B4) holds generally. Note that the divergence in (B4) acts on the grid scale. The diagnosis of the contribution of each subdomain to the grid-mean divergence requires an assumption regarding the fraction of ∂Ω T covered by each ∂Ω g i . Here, we assume that A g i ¼ a i A g T , where A g i and A g T are the areas of surfaces ∂Ω g i and ∂Ω T , respectively. We further assume that for each Ω i the average over ∂Ω g i equals the subdomain mean. From this it follows that Z Note that (B5) cannot be obtained from the divergence theorem, since ∂Ω g i is not a closed surface. Using (B5) and dividing by the grid box volume V T , we can rewrite (B2) as where u r ¼ u − u b . Since the vertical extent of the volumes is fixed at the model vertical resolution, V i / V T ¼ ⟨A i ⟩/A T ¼ a i , with a i as the area fraction.
The net entrainment flux can be written in terms of a contribution from net mass entrainment and a contribution due to the subfilter-scale flux of ϕ: Here,ð · Þ represents the average over interface ∂Ω sg i , u r,n ¼ u r · n, and A sg is the total area of surface ∂Ω sg i . The two terms on the right-hand side of (B7) are denoted as net dynamical and turbulent entrainment 10.1029/2020MS002162

Journal of Advances in Modeling Earth Systems
fluxes, respectively. The net dynamical entrainment flux is taken to be the sum of two terms. For mass, it is written as and for a scalar as where the entrainment E ij and the detrainment Δ ij are positive semidefinite. We make the upwind approximation that the exchanged air mass carries with it the property of the subdomain from which it emanates, as is common in parameterizations (de Rooy et al., 2013).
The turbulent entrainment flux does not involve mass exchange between subdomains, and it is modeled as shown in section 3.2: Here, Ê ij is the turbulent entrainment rate from the jth subdomain into the ith subdomain. Using (B9) and (B10), decomposing the divergence term into vertical and horizontal components, and applying the eddy diffusivity assumption for the vertical turbulent flux, (B6) is written in the form (20). By setting ϕ ¼ 1 in (20), the mass continuity (i.e., area fraction) Equation (18) follows.
The second-moment equations can be derived by first writing (B6) for the product of two scalars ϕψ. Using (B7), and decomposing the divergence term into vertical and horizontal components, we obtain The subdomain covariance equation can then be obtained from (20), (18), and (B11) as which leads to Here, terms of the form ðψ i −ψÞ d u ′ r; n ϕ ′ are written as ψ i * d u ′ r; n ϕ ′ to ensure conservation of second moments on the host model grid. The last term in (B13) follows from (B12), given that The dissipation of covariance is represented by D ϕ ′ ψ ′ ; i . The vertical subgrid covariance flux is written as downgradient and proportional to the eddy diffusivity K ϕψ,i :

Journal of Advances in Modeling Earth Systems
Substituting (B9), (B10), and (B15) in (B13), we obtain (22). The extended EDMF scheme only makes use of covariance equations for thermodynamic variables θ l and q t and for the turbulence kinetic energy. Subgrid-scale covariances between thermodynamic variable and momentum are modeled diffusively following (14). ϵ i0 þε i0 ¼