A Generalized Mixing Length Closure for Eddy‐Diffusivity Mass‐Flux Schemes of Turbulence and Convection

Because of their limited spatial resolution, numerical weather prediction and climate models have to rely on parameterizations to represent atmospheric turbulence and convection. Historically, largely independent approaches have been used to represent boundary layer turbulence and convection, neglecting important interactions at the subgrid scale. Here we build on an eddy‐diffusivity mass‐flux (EDMF) scheme that represents all subgrid‐scale mixing in a unified manner, partitioning subgrid‐scale fluctuations into contributions from local diffusive mixing and coherent advective structures and allowing them to interact within a single framework. The EDMF scheme requires closures for the interaction between the turbulent environment and the plumes and for local mixing. A second‐order equation for turbulence kinetic energy (TKE) provides one ingredient for the diffusive local mixing closure, leaving a mixing length to be parameterized. Here, we propose a new mixing length formulation, based on constraints derived from the TKE balance. It expresses local mixing in terms of the same physical processes in all regimes of boundary layer flow. The formulation is tested at a range of resolutions and across a wide range of boundary layer regimes, including a stably stratified boundary layer, a stratocumulus‐topped marine boundary layer, and dry convection. Comparison with large eddy simulations (LES) shows that the EDMF scheme with this diffusive mixing parameterization accurately captures the structure of the boundary layer and clouds in all cases considered.


Introduction
Turbulence is ubiquitous in the planetary boundary layer. Small-scale chaotic air motions enhance mixing of energy and moisture in the lower troposphere. Under statically unstable conditions, convective updrafts and downdrafts further increase the vertical transport of energy and moisture between the surface and the air aloft. Together, turbulence and convection shape the vertical distribution of temperature and water vapor that sustains clouds. However, these processes act on scales far too small to be resolved in global climate models (GCMs), with resolutions constrained by current computational power . Although the unabated increase in processing power will make globally resolving deep convective processes routine in the coming years (Kajikawa et al., 2016), resolving turbulent mixing and shallow convection will remain an intractable problem for decades. Instead, parameterizations have to be used to approximate the average effect of these subgrid-scale processes on the grid scale.
Conventional parameterizations consider atmospheric turbulence and convection as independent processes, neglecting interactions that alter their combined effect on the large scale. These parameterizations are often regime dependent, leading to models that artificially split the spectrum of atmospheric conditions into a discrete number of cases. Examples of such case-dependent approaches include parameterizations of cumulus (Arakawa, 2004) and stratocumulus clouds (Lilly, 1968;Schubert, 1976). However accurate, the use of disparate schemes for different conditions complicates a seamless representation of subgrid-scale processes in the lower troposphere.
Several approaches to obtain a unified model of turbulence and convection have been proposed (Lappen & Randall, 2001;Park, 2014;Thuburn et al., 2018). Here we focus on the extended formulation of an eddy-diffusivity mass-flux (EDMF) scheme developed in Tan et al. (2018), which in turn built on work by Siebesma and Teixeira (2000), Soares et al. (2004), Siebesma et al. (2007), and Angevine et al. (2010), among others. In the EDMF framework, the flow within each grid cell is decomposed into several distinct subdomains, representing coherent convective structures and their relatively isotropic turbulent environment. Convective transport is captured by mass flux terms that depend on differences between subdomain mean properties; more isotropic turbulent transport, associated with small-scale fluctuations within each subdomain, is captured by eddy diffusion closures.
The extended EDMF framework uses additional prognostic equations for subdomain variables, such as the environmental turbulence kinetic energy (TKE), and it requires closures for local turbulent fluxes and for the mass exchange between subdomains (Tan et al., 2018). Even though the EDMF framework arises from the need for a unified model of turbulence and convection, the parameterizations used for entrainment and turbulent mixing are usually defined differently for each regime (Suselj et al., 2013;Witek et al., 2011). The development of regime-independent parameterizations for the required closures is the last step in the construction of a unified model of atmospheric turbulence and convection.
Here, a regime-independent closure for turbulent mixing within the EDMF framework is proposed. Section 2 reviews the decomposition of subgrid-scale fluxes in the extended EDMF scheme. Section 3 introduces the formulation of the closure. Section 4 illustrates the performance of the EDMF scheme with the turbulent mixing closure in boundary layer regimes where vertical transport is strongly dependent on the turbulence closure used: the stable boundary layer (SBL), the stratocumulus-topped boundary layer (STBL), and dry convection. The performance of the extended EDMF scheme with this closure in moist convective cases is demonstrated in a companion paper (Cohen et al., 2020). Finally, section 5 summarizes the results and conclusions.

EDMF Framework
In the EDMF framework, each grid cell volume is decomposed into n updrafts or downdrafts (labeled by index i ¼ 1; …; n) and an environment (labeled by index i ¼ 0) in which they are embedded. Following this decomposition, the grid-mean value of variable ψ may be written as Here, angle brackets ⟨·⟩ denote the grid mean, ψ i denotes the Favre average of ψ over subdomain i, and a i is the mean horizontal cross-sectional area covered by subdomain i within the grid cell. This partition is motivated by the anisotropy of turbulent convective flows, in which isotropic turbulent eddies coexist with coherent columnar structures that induce a strong vertical transport (Bjerknes, 1938). The subdomain decomposition is simplified for the horizontal velocity vector u h , which is taken to have the same mean value for all subdomains, u h; i ¼ ⟨u h ⟩. Applying the subdomain decomposition to higher-order moments introduces additional terms associated with the difference between grid and subdomain means. For the vertical subgrid-scale flux of ψ, this leads to Here, w is the vertical velocity, ψ * ¼ ψ − ⟨ψ⟩, ψ ′ i ¼ ψ − ψ i , and ψ The subdomain-mean terms can be explicitly solved for by introducing n prognostic subdomain equations for each variable and an additional equation for each plume area fraction a i , which may be diagnostic or prognostic. Cohen et al. (2020) derive the subdomain equations used in the EDMF framework, starting from the Navier-Stokes equations. The use of prognostic subdomain equations means that convective fluxes such as w * i ψ * i in (2) are explicitly solved for, while turbulent fluxes like w ′ i ψ ′ i must be modeled. Turbulent fluxes within each subdomain are modeled as downgradient and proportional to an eddy diffusivity K ψ, i , where ψ is the property being transported. For the vertical turbulent flux in (2), this gives The eddy diffusivity K ψ, i is proportional to a characteristic velocity scale and the length scale of the eddies driving the transport, both of which must be parameterized.
Proposed closures for the eddy diffusivity vary from simple diagnostic expressions to second-order models that introduce prognostic equations for both scales (Umlauf & Burchard, 2003). The 1.5-order TKE model is a particularly popular choice due to its balance between accuracy and computational efficiency (Mellor & Yamada, 1982). The 1.5-order model, also referred to as the Level 2.5 model in the Mellor-Yamada hierarchy, makes use of a prognostic equation for TKE and a diagnostic expression for the mixing length. In the EDMF framework, the grid-mean TKE ⟨e⟩ can be decomposed following Expression (2) for second-order moments as where ē i is the TKE of subdomain i and the second term represents the corresponding convective kinetic energy. This expression can be simplified by assuming that for the updrafts and downdrafts (i > 0), the contribution to the grid-mean TKE from small-scale turbulence is negligible compared to the convective term, an assumption commonly made in EDMF schemes: Thus, grid-mean TKE is given by the sum of environmental TKE and convective TKE. The TKE decomposition (5) can also be obtained by assuming a small updraft and downdraft area fraction and similar turbulence intensity in all subdomains (Siebesma et al., 2007). However, the equations derived for the subdomain second-order moments with these two approaches differ in the source terms that appear due to entrainment processes between subdomains. The former approximation is favored here to allow for the use of this framework in high-resolution models, where the assumption of slender updrafts may become inadequate (Randall, 2013).
Given an updraft area fraction a i , which may be diagnostic or prognostic (Tan et al., 2018), the grid-mean TKE is determined by the environmental TKE ē 0 and the subdomain-mean vertical velocities w i . The subdomain-mean vertical velocity equation for subdomain i is where ∇ h is the horizontal gradient operator, Ψ ¼ p=ρ is the pressure potential, and the turbulent transport terms on the right-hand side are negligible for all subdomains except the environment ( i ¼ 0 ). Subgrid density changes are only considered in the buoyancy term, such that ρ ¼ ⟨ρ⟩ in the previous equation, in order to avoid creation of spurious acoustic modes through the subdomain decomposition (Cohen et al., 2020). The buoyancy b i and the pressure potential anomaly Ψ † i are defined with respect to a reference hydrostatic pressure profile p h (z) and density ρ h (z), related by ∂ z p h ¼ −ρ h g: Here, p i is the subdomain-mean pressure. Density appears inside the pressure gradients in (6) and (7) to ensure thermodynamic consistency of the subgrid-scale anelastic approximation (Cohen et al., 2020).
Interactions between subdomains are captured by entrainment and detrainment fluxes. In the vertical velocity Equation 6, Δ ij is the dynamical detrainment of air mass from subdomain i into subdomain j, and E ij and E ij are the dynamical and turbulent entrainment from subdomain j into subdomain i, respectively. It is assumed that entrainment events occur over timescales much shorter than the eddy turnover rate K ψ; i =ē i , so that entrained air carries the properties of the subdomain it detrains from. In addition, for now we assume that entrainment occurs only between convective plumes and the environment, not among plumes.
The prognostic equation for environmental TKE can be written in non-conservative form as (Cohen et al., 2020) Here, ⟨u⟩ and ⟨v⟩ are the components of ⟨u h ⟩, P is the velocity-pressure gradient correlation, and D is the turbulent dissipation. All sources and sinks of ē 0 account for unresolved processes on the grid scale, so they must be parameterized. Subdomain covariances in (8) are modeled diffusively, with the environmental eddy diffusivity K ψ defined as where l is the mixing length and c ψ is a fitting parameter. The subscript 0 in the eddy diffusivity is dropped to simplify notation. The coefficient c ψ is taken to be equal to c h for the diffusion of all fields except for momentum, for which c ψ = c m . The eddy viscosity K m is related to K h through the turbulent Prandtl number Pr t , such that K m ¼ Pr t K h .
Under the assumption that subgrid-scale pressure work on the grid mean is negligible, P is taken as opposite to the pressure work on the plumes (Tan et al., 2018), The last term in (10) appears as a sink term in the convective TKE balance, which is derived in Appendix B. Hence, P acts as a return-to-isotropy term on the full grid, transferring momentum from the strongly anisotropic coherent structures into the relatively isotropic eddies in the environment. The pressure work on the plumes is formulated in terms of contributions from a virtual mass term (Gregory, 2001), an advective term, and a drag term (Romps & Charn, 2015), yielding the following expression for the velocity-pressure gradient correlation: where α a and α d are constant parameters, H i is the plume height, and α b is a function of the aspect ratio of the plume. Finally, assuming statistical equilibrium at scales l (Vassilicos, 2015), turbulent dissipation can be estimated from the spectral transport relation that follows from Kolmogorov's theory of inertial turbulence, giving Taylor's dissipation surrogate Here, c d is an empirical coefficient and l is the dissipation length, taken to be equal to the mixing length in our model. Expressions (3) and (5)-(12) provide closure to a 1.5-order model of turbulence within the EDMF framework, given diagnostic expressions for the mixing length l and for entrainment and detrainment.

Mixing Length Formulation
We seek to obtain a regime-independent eddy diffusivity closure that provides an accurate representation of turbulent subgrid-scale fluxes, over a wide range of host model resolutions. Thus, the eddy diffusivity should reduce to a large eddy simulation (LES)-type closure at high resolution, while being able to account for the processes that modify turbulent fluxes at larger scales. The formulation of the closure is organized following this logic.
In section 3.1, we first adapt a minimum TKE dissipation closure proposed for LES subgrid models (Abkar & Moin, 2017) to the EDMF framework. Given the diffusive closure (3) and the eddy diffusivity (9), the minimum dissipation assumption can be used to construct a mixing length closure. This mixing length closure is shown to be equivalent to other proposed closures (e.g., Grisogono, 2010) for stable stratification, but additional entrainment terms appear in the general case. Section 3.2 highlights the limitations of this closure for climate modeling and weather prediction purposes when a prognostic TKE equation is used. Section 3.3 then introduces a modified mixing length closure, which builds on the minimum dissipation model and corrects its shortcomings by introducing additional mechanisms of net TKE dissipation.

Minimum Dissipation of Environmental TKE
As in Verstappen (2011) and Abkar and Moin (2017), we assume that at the small scales represented by the environment in the EDMF scheme, TKE is dissipated at least at the rate at which it is produced. This condition translates into an inequality for the production and dissipation terms in the environmental TKE budget (8): Here, the terms involving TKE injection from entrained air are also taken to be locally balanced by dissipation, consistent with the assumption that entrainment events occur over timescales much shorter than the eddy turnover time K ψ; i =ē i . Note that the net dissipation condition (13) does not include redistribution terms, such as the turbulent transport or the velocity-pressure gradient correlation P. Moreover, the inequality (13) represents a local condition for the environment, and it does not preclude net subgrid-scale energy production due to processes such as convection, represented by plumes. Denoting the difference between the right-hand side and the left-hand side of (13) as the net environmental TKE dissipation γ 0 , the prognostic environmental TKE Equation 8 reduces to Here, P captures the effect of plumes on the environmental TKE. The evolution of the grid mean TKE that follows from decomposition (5) and the simplified prognostic Equation 14 is 10.1029/2020MS002161

Journal of Advances in Modeling Earth Systems
A detailed derivation of this equation and the subgrid-scale kinetic energy pathways in the extended EDMF scheme is described in Appendix B. Under the net dissipation closure (13), grid mean TKE production occurs through the first two terms on the right-hand side of (15): the convective buoyancy flux and the subdomain-scale shear production.
The net dissipation condition (13) can be written in terms of the mixing length by introducing the closures described in section 2. Using Taylor's dissipation surrogate (12) and downgradient closures for the shear and buoyancy terms of the form the inequality (13) leads to a condition for the maximum value of the mixing length l at which the net dissipation γ 0 is still positive semidefinite: Here, the environmental buoyancy gradient is computed following Tan et al. (2018), taking into account possible phase change effects. In (16) and (17), x k and u k, 0 represent the k-th coordinate and k-th velocity component in the environment, respectively. From the inequality (17), an expression for the mixing length that minimizes turbulent dissipation can be obtained by solving for l. This is equivalent to setting γ 0 ¼ 0 in (14) and (15). For the resulting value of the mixing length, production and dissipation of TKE are locally balanced: Here, Δ is the discriminant and the different terms are given by In (18), the product ðS þ BÞD is independent of the mixing length, so l tke can be readily evaluated. Although the term ðS þ BÞ is sign indefinite, the discriminant (18) can be shown to remain positive semidefinite even when the shear and buoyancy terms result in TKE destruction, provided that the inequality (13) holds. This is because the minimum dissipation balance requires so that the expression for the discriminant Δ is of the form The mixing length l tke depends on local characteristics of the environment and on the vertical velocity difference between subdomains, which enters the injection term I in (19). Hence, convection modifies the environmental diffusive transport directly through entrainment processes. In addition, convection also regulates the time evolution of turbulent fluxes through its effect on the prognostic environmental TKE Equation 14, captured by P.
This approach can also be applied to turbulence models that retain covariance terms w ′ i ψ ′ i for other subdomains and not only for the environment. In this case, the minimum dissipation condition may be used to obtain a characteristic mixing length l tke, i for each subdomain. However, variance within plumes can also be accounted for by variance among plumes when the number of subdomains is increased.
In stably stratified boundary layers, where convection is inhibited, pressure work and entrainment fluxes in (6) act to homogenize the different subdomains, such that ψ * i →0 for any variable ψ and a 0 → 1 (i.e., there are no convective plumes). Under these conditions, the minimum dissipation mixing length (18) reduces to the expression proposed by Grisogono (2010) for steady-state SBL flow: The balance between shear production, destruction due to stratification, and dissipation, which arises when using this mixing length, is a well-known leading-order state in neutral (Spalart, 1988) and moderately stable boundary layer flows (Li et al., 2016).

Limitations of the Minimum-Dissipation Closure
Expression (18) for the mixing length l tke captures the leading-order balance in the environmental TKE budget at small scales. However, a model with a diffusive closure based on l tke cannot fully describe the dynamics of the boundary layer at the coarse resolutions typical of GCMs, on the order of 10 4 m in the horizontal and 10-100 m in the vertical. At these scales, the resolved horizontal gradients are weak, and the environmental TKE Equation 14 that results from using l tke can be simplified using the boundary layer approximation (neglecting horizontal relative to vertical derivatives): Note that we set γ 0 ¼ 0 to obtain (23), since we are considering the case where l = l tke and production locally balances dissipation. In stable conditions (P ¼ 0), integrating the conservative form of (23) from the surface layer (z s ) to the free troposphere above (z i ) yields the evolution equation for the vertically integrated environmental TKE: In stable conditions, a 0 ≈ 1 and ψ * i ≈ 0 for any variable ψ. In addition, the absence of plumes implies that detrainment and entrainment processes are negligible. From (24), it follows that the evolution of the vertically integrated TKE under the minimum dissipation closure only depends on the flux from the unresolved surface layer in stable conditions. But unbalanced TKE dissipation has been observed to become increasingly important as stratification develops in field studies of the atmospheric boundary layer (Li et al., 2016), and it can be expected to play a role in conditions of strong surface cooling. The budget (24) cannot capture unbalanced TKE destruction within the boundary layer due to stratification. Furthermore, the minimum dissipation mixing length l tke leads to enhanced eddy diffusion with increasing stratification, as deduced from (22). This is contrary to the evidence of turbulent mixing being inhibited in strong stratification, such as near strong inversions.
The limitations of a minimum dissipation model also become apparent in convectively unstable boundary layers. Integrating the TKE Equation 23 in the vertical, the evolution of the vertically integrated environmental TKE in convective conditions reads Here, the last term only accounts for changes in environmental area fraction and does not result in a source or sink of ē 0 (Tan et al., 2018). A major difference between the SBL budget (24) and the convective budget (25) is the contribution of the velocity-pressure gradient correlation. From the velocity-pressure gradient relation (11), pressure work captures the important energization of turbulence in the environment owing to ascending or descending plumes (Schumann & Moeng, 1991). At the grid scale, the source of this subgrid-scale energy term is the convective buoyancy flux in (15), which accelerates the plumes in convective conditions.
The TKE balance (25) shows that, in convective conditions, the source of environmental TKE from updrafts or downdrafts can only be compensated by the flux from the unresolved surface layer. This is often a source term rather than a sink term, because shear production is surface intensified. Thus, the TKE balance (25) suggests an unbalanced growth of TKE in convective boundary layers. This continuous TKE increase in convective conditions is inconsistent with LES results showing quasi-stationary TKE levels in convective boundary layers (Nieuwstadt et al., 1993).
The TKE balances (24) and (25) highlight the shortcomings of the minimum dissipation balance (18) as a general closure for diffusive mixing in the boundary layer in stable and convective conditions. The lack of net dissipation mechanisms in the vertically integrated TKE balance hinders the correct representation of important processes, such as the shallowing of the boundary layer in the late afternoon or the sharp mixing inhibition near inversions. Moreover, it precludes reaching a quasi-stationary state in statically unstable boundary layers. Nevertheless, the limitations of the minimum dissipation model can be used to inform the construction of a generalized master length scale based on the TKE production-dissipation inequality (13).
The limitations of the minimum dissipation balance showcased in this section are not necessarily applicable to other turbulence models. For example, He et al. (2019) use the production-dissipation condition to diagnose TKE and eddy diffusivity from a mixing length l. This allows an instantaneous adjustment of TKE to a new balanced state, at the cost of representing convection with an empirical parameterization that has no subgrid interaction with turbulent diffusion.

Constrained Minimization of TKE Dissipation
A master length scale that corrects the limitations of the minimum-dissipation model can be constructed by taking dissipation to be higher than production under certain circumstances. Using closures of the form (9) and (16) for the turbulent fluxes and (12) for the dissipation, it follows from the production-dissipation inequality (13) that excess dissipation occurs for l < l tke . Hence, unbalanced TKE dissipation arises naturally in regions of the boundary layer where the characteristic size of environmental eddies is constrained to be smaller than l tke . A general mixing length capturing this condition can be written as where l j (j ¼ 1; 2; …; N) are candidate mixing lengths based on flow constraints and s min ðxÞ is a smooth minimum function defined in Appendix A. The TKE production-dissipation inequality (17) with the closures substituted implies that the minimum length scale (26) provides maximum TKE dissipation among the candidate length scales. Thus, the use of the minimum length scale (26) is equivalent to the minimization of TKE dissipation in (13) subject to the constraint that dissipation exceeds the candidate dissipation rates, where Dj l¼lj is the candidate dissipation rate evaluated at length scale l j .
Our suggestion for choosing a general mixing length as a smooth minimum of various candidates contrasts with the common practice (e.g., Han & Bretherton, 2019;He et al., 2019) to use the expression suggested by Blackadar (1962), for a master length scale l h . This length scale l h , proportional to the harmonic mean of the candidates l 1 and l 2 , is smaller than both l 1 and l 2 . If closures similar to (9) and (12) are used in a prognostic equation for TKE, the mixing length (28) results in an unrealistic intensification of TKE dissipation in regions where the candidate length scales l 1 and l 2 are similar. This undesirable characteristic is avoided by using the smooth minimum (26).
We consider two limiting factors for the characteristic length scale of turbulent motion in boundary layer flows: stable stratification and the distance to solid boundaries.

Stratification Constraints
Environmental stratification constrains the size of turbulent eddies by inhibiting the vertical displacement of air masses. Stably stratified turbulence is known to show high vertical variability and reorganization into layered structures, with most mixing occurring within the layers (Waite, 2011). The thickness of these layers is determined by the vertical scale at which the governing dynamic equations become self-similar (Augier et al., 2012;Billant & Chomaz, 2001), known as the buoyancy scale l b . For a flow with an imposed stratification given by the Brunt-Väisälä frequency N e , this length scale is where c b is an empirical coefficient. It is important to note that imposing l b as an upper bound for the size of eddies is similar to doing so by the Ozmidov scale l o ∼ ffiffiffiffiffiffiffiffiffiffiffiffi ffi D=N 3 e q only if turbulent motions at the scale in question are assumed to be in the inertial subrange, such that (12) holds. In this case, an expression equivalent to (29) for the Ozmidov scale is However, recent experimental studies suggest that under strong stratification, turbulence may not display an inertial subrange (Grachev et al., 2013). In that case, Expression (12) and the Ozmidov scale (30) may not be applicable (Li et al., 2016), whereas the buoyancy scale (29) still holds.
The buoyancy frequency of moist air depends on the latent heat release and evaporative cooling associated with the vertical displacement of air parcels. In general, the effective static stability N e lies between the dry and the moist adiabatic limits. In the same spirit as O'Gorman (2011), we use an effective static stability of the form where θ v is the virtual potential temperature, λ represents the area fraction of environmental air undergoing phase change, and θ vl is the liquid water virtual potential temperature, defined below. In Expression (31), we have used the definition Expression (31) differs from that presented in O'Gorman (2011) in the diagnosis of λ and the use of θ vl; 0 instead of the saturated equivalent potential temperature. These differences arise from the much smaller vertical scale considered here. At scales of 10-100 m, the occurrence of phase changes is not necessarily correlated with the sign of the vertical velocity of air parcels (O'Gorman, 2011;O'Gorman & Schneider, 2006). Thus, λ cannot be diagnosed from vertical velocity statistics. In the nonprecipitating cases considered here, λ is given by the environmental cloud fraction f c, 0 . Cloud fraction diagnosis is cloud-type dependent in many current GCMs (Collins et al., 2004). In our EDMF scheme, we use a regime-independent probabilistic cloud scheme (see Appendix C).
The liquid water virtual potential temperature θ vl appearing in the effective static stability measures the buoyancy of cloudy air parcels when moist adiabatically returned to clear conditions (Grenier & Bretherton, 2001;Marquet, 2011), We use θ vl instead of the saturated equivalent potential temperature because θ vl converges to θ v in the dry limit, while also including the effects of latent heat release. Here, η ¼ R v =R d − 1, L v is the latent heat of vaporization, c p is the specific heat of air, q t and q l are the total and liquid water specific humidities, θ l is the liquid water potential temperature, T is the temperature, and R v , R d are the gas constants for water vapor and dry air, respectively. Expression (33) can be used to evaluate ∂θ v; 0 =∂θ vl; 0 in (31) and (32).
Note that the effective static stability (31) converges to the dry limit when q l → 0 for all values of λ; it reduces to N 2 e ¼ ð1 − λÞN 2 , with dry buoyancy frequency N, in conditions that are well mixed in θ l and q t .

Wall Constraints
The presence of boundaries also imposes an upper limit on the size of eddies near them. Following Monin and Obukhov (1954), the eddy diffusivity in the surface layer has the form where ξ = z/L, ϕ ψ (ξ) is an empirical stability function, κ is the von Kármán constant, L is the Obukhov length, and u * is the friction velocity. The upper bound for the mixing length near the surface is obtained by matching this eddy diffusivity with the Expression (9) for the eddy diffusivity: Here, κ * ¼ ē 1=2 0 =u * is the ratio of rms turbulent velocity to the friction velocity in the surface layer. The friction velocity in our model is diagnosed using the flux-profile relationships of Byun (1990), except in free convective conditions. When the conditions for free convection are satisfied, the diagnostic of u * , which is a function of the horizontal wind speed at the lowest model level, is modified following Beljaars (1995).
The choice of a common master length for momentum and tracer diffusion implies c h ϕ h = c m ϕ m , such that ϕ h ¼ Pr t ϕ m . In our formulation, the turbulent Prandtl number is taken to be a function of the gradient Richardson number Ri, based on a simplified cospectral budget of momentum and heat transport (Katul et al., 2013;Li, 2019): 10.1029/2020MS002161

Journal of Advances in Modeling Earth Systems
Here, ω 2 ¼ 53=13 is a phenomenological constant, and Pr t, 0 is the turbulent Prandtl number in neutral conditions. The stability function ϕ m is often written in the form (Businger et al., 1971;Nakanishi, 2001) where H(·) is the Heaviside function and a − i ; a þ i are empirical functions. The values of a − i are taken as negative definite to reflect the convective elongation of eddies in unstable conditions. In stable conditions, self-similarity of the flow requires that a þ 2 ¼ 1 and a þ 1 > 0, such that under strong stratification, the mixing length (35) becomes independent of ξ. As shown by Monin and Obukhov (1954), the asymptotic limit of ϕ m under strong stratification also requires that a þ 1 ¼ Pr t ðRi st Þ=Ri st : Here, Ri st is the asymptotic Richardson number at ξ ≫ 1=a þ 1 in the surface layer. The empirical function (37) has been shown to become increasingly inaccurate with stability for ξ > 0.5 (Optis et al., 2016;Sorbjan & Grachev, 2010). Moreover, extending the use of the limiting scale l w above the surface layer precludes the use of a þ 1 ≠ 0 in stable conditions, since the Obukhov length characterizes stratification only in the constant flux layer near the surface. Although the use of l w in Expression (26) mandates a þ 1 ¼ 0, the effect of stability in eddy diffusion is still captured by l b . In the constant flux layer, the limiting length l b is equivalent to the use of the empirical function (37) in the strongly stable limit, with Under weaker stratification, turbulence in the surface layer can reach a quasi-steady state (Spalart, 1988). In this case, the limiting scale l w should converge to l tke . Assuming that entrainment processes are limited to dynamical entrainment by the plumes in the surface layer, the ratio of the two length scales can be written as which is constant under neutral stratification and is slowly varying with Ri due to the opposing effect of the Prandtl number (36). From (39), the convergence of l tke to l w in the surface layer is satisfied for ðc d c m Þ 1=2 κ 2 * ≈ 1. The use of a soft minimum function for the mixing length (26) allows for a smooth transition from Monin-Obukhov similarity theory near the surface to a local turbulent closure farther away from it, where the use of Monin-Obukhov scaling may be inaccurate (Optis et al., 2016). In addition, Expressions (38) and (39) show that this transition is asymptotically consistent.

Master Mixing Length
Finally, the smooth minimum of the three candidate length scales determines the mixing length, where and The proposed diffusive closure is implemented using Equations 40-43, as well as the prognostic environmental TKE Equation 8 in flux form (see Equation B1). The mixing length (40) depends on a group of nondimensional parameters C that must be obtained empirically:

Journal of Advances in Modeling Earth Systems
Values for these parameters are reported in studies of boundary layer turbulence, obtained from field observations (Businger et al., 1971) or LES (Nakanishi, 2001). However, the direct use of some of these values in the EDMF scheme is not justified due to the decomposition of the subgrid-scale flow into different subdomains. Because of the large size of the parameter space C and the presence of other parameters in the EDMF scheme, we limit the parameter optimization process to C * ¼ fc m ; c d ; c b g in this study. C * contains the parameters that appear in the closures that are most strongly affected by the domain decomposition. All other parameters in C, except Pr t, 0 , arise from similarity theory arguments for the unresolved surface layer. Here, it is assumed that similarity arguments are valid outside convective updrafts, and all values are taken from Nakanishi (2001). For the simulations reported in the next section, the parameter space used is shown in Table 1. The rest of parameters used in the EDMF scheme, which do not appear explicitly in the formulation of the mixing length closure, are reported in Cohen et al. (2020).

Results for Single-Column Simulations
Here we focus on case studies targeting the simulation of the Arctic SBL, stratocumulus clouds, and dry convection. The performance of the extended EDMF scheme in moist convective cases is explored in Cohen et al. (2020), using the same set of closures and parameters. The extended EDMF scheme is tested for horizontal resolutions typical of GCMs. Invoking the boundary layer approximation (neglecting horizontal derivatives), we perform simulations in a single-column model (SCM). The SCM is a one-dimensional vertical model that aims to represent a single atmospheric column within a GCM. Results from single-column simulations using the extended EDMF scheme are then compared to horizontal averages obtained from LES over the same domain. LES are set up by further discretizing the atmospheric column horizontally and using horizontal doubly periodic boundary conditions. The EDMF scheme used here differs from the one described in Tan et al. (2018) in the parameterizations of the eddy diffusivity K ψ , the vertical pressure anomaly gradients in (6) and (10), entrainment and detrainment, and the addition of turbulent entrainmentÊ ij . The parameterization of the eddy diffusivity follows (9) and (40)-(43). The entrainment parameterization is described in Cohen et al. (2020), and the treatment of the pressure anomaly term is shown in (11). In addition, although the theoretical framework presented here allows for the use of downdrafts, the implementation used in this section decomposes the domain solely into one updraft and its turbulent environment.
LES are performed using PyCLES, an anelastic fluid solver in which the subgrid-scale fluxes are treated implicitly by the weighted essentially non-oscillatory (WENO) scheme used to discretize the prognostic equations (Pressel et al., 2015). Implicit LES using WENO numerics have been shown to result in higher effective resolution than other combinations of numerics and explicit SGS closures . Finally, LES results from previous model intercomparison projects are also reported where available.

Stable Boundary Layer
Statically stable conditions in the boundary layer inhibit convection, reducing the EDMF scheme to a diffusive closure. In the implementation of the scheme, this translates to conditioning the surface updraft area fraction on the sign of the surface buoyancy flux, such that it becomes zero in conditions of surface cooling. With no updrafts or downdrafts, the only contribution to the subgrid-scale flux (2) comes from the environmental downgradient turbulent flux (3). This leads to a high sensitivity of SCM simulations to changes in the mixing length formulation. Here we focus on the GEWEX Atmospheric Boundary Layer Study (GABLS), discussed in Beare et al. (2006).

Simulation Setup
The initial and boundary conditions of the simulation are adapted from observations during the Beaufort and Arctic Storms Experiment (Curry et al., 1997) and follow Beare et al. (2006). The velocity field is initialized as ð⟨u⟩; ⟨v⟩Þ ¼ ðu g ; 0Þ, where the geostrophic velocity is u g ¼ 8 m s −1 . The initial temperature sounding

Journal of Advances in Modeling Earth Systems
is given by a mixed layer with potential temperature θ ¼ 265 K up to 100 m, overlain by an inversion with a potential temperature gradient of 10 K km −1 . The surface boundary condition is given by constant cooling, For both the SCM and LES, the domain height is 400 m. In the LES configuration, the domain spans 400 m in both horizontal directions as well. The LES data are generated using an isotropic mesh with Δx i ¼ 3:125 m resolution, which translates into 2 × 10 6 degrees of freedom. The full range of LES results from Beare et al. (2006), using the same resolution, is also included for reference. The SCM simulations are performed at vertical resolutions of Δz ¼ 3:125 m, 12.5 m, and 50 m (128, 32, and 8 degrees of freedom, respectively). This range characterizes the performance of the EDMF scheme both at high resolution and for coarser resolutions typical of regional and global climate models in the lower troposphere. The time steps used in the SCM simulations, in order of increasing Δz, are Δt ¼ 5, 15, and 60 s. Figure 1 shows vertical profiles of ⟨θ⟩, ⟨u⟩, and ⟨v⟩ time averaged over the ninth hour of simulation. The EDMF scheme captures well the boundary layer height and the intensity of the low-level jet, with little resolution dependence of the mean profiles up to Δz ¼ 12:5 m. At 50 m resolution, the SCM predicts a slightly deeper boundary layer. The EDMF-simulated TKE follows closely the LES data, as shown in Figure 2. The time series show periods of TKE growth due to the subgrid momentum flux from the surface layer and periods of decay due to the increasing stratification. These changes in vertically integrated TKE are much smaller than the integrated TKE production and dissipation terms, as shown in Figure 3. The domain-mean TKE budget, which coincides with the environmental budget for stable conditions, is shown in Figure 3.

Results
The two main causes of grid sensitivity at 50 m resolution are the inability to capture the region of maximum shear production close to the surface and the deterioration of the friction velocity diagnosis. The effect of the former can be observed in Figure 3. Even if the budget is correctly captured above 50 m, the absence of grid cells at the lower levels results in a significant reduction of the vertically integrated production and dissipation. In addition, the diagnosis of u * based on Byun (1990) overestimates the friction velocity at coarser resolutions. This can be observed by comparing the normalized TKE profile to the vertically integrated time series in Figure 2.
The dominant mixing length throughout the simulation is shown in Figure 2 for all heights. Initially, the wall-limited mixing length l w is dominant below the inversion, due to the absence of mean shear and stratification. As the shear and stratification develop, the dominant mixing length profile attains a three-layered structure. Closest to the bottom boundary, the distance to the wall constrains the size of eddies. Farther away from the surface, the mixing length is determined by the local TKE balance. As stratification

Journal of Advances in Modeling Earth Systems
increases with height, the stratification-limited mixing length l b becomes dominant, depleting TKE and limiting turbulent mixing. The eddy diffusivity, shown in Figure 2, is maximum near the transition from l tke to l b , where the mixing length is largest. Again, the overestimation of the friction velocity and the absence of grid points in the lower layers result in an overestimation of the eddy diffusivity at coarse resolutions.
Both the LES and EDMF budgets show the quasi-balance of TKE sources and sinks throughout the boundary layer, even in regions where l tke is not dominant. The downgradient parameterization of shear production S, buoyant production B, and the turbulent transport T results in profiles that follow closely the LES data, particularly at higher resolution. This validates the assumptions used to model the second-order moments in the extended EDMF scheme under stable stratification.

Stratocumulus-Topped Boundary Layer
The ability of the extended EDMF scheme to represent the dynamics of the STBL is tested using as a baseline the second Dynamics and Chemistry of Marine Stratocumulus (DYCOMS-II) field study (Stevens et al., 2003), performed near the coast of San Diego, California. In particular, the conditions observed during the first research flight (RF01) are considered, for which precipitation was not observed.

Simulation Setup
The simulation setup for DYCOMS-II RF01 is reported in Stevens et al. (2005). The base state is initialized with a two-layer structure in θ l and q t , separated by a strong inversion at z i ¼ 840 m. The bottom layer is well mixed in both conserved variables, with saturation and cloud formation occurring above 600 m. The cloud top, located at z i , is characterized by Δθ l ¼ 8:5 K and Δq t ¼ −7:5 g kg −1 . The free troposphere is warmer and drier than the mixed layer, with a θ l lapse rate of (1/3)(z − z i ) −2/3 K m −1 and constant q t . The surface sensible and latent heat fluxes are set to 15 and 115 W m −2 , respectively. The vertical

Journal of Advances in Modeling Earth Systems
water distribution induces radiative cloud base warming and radiative cooling at cloud top and in the free troposphere.
The domain height is 1.5 km. In the LES, the horizontal domain extent is set to 3.36 km. The resolution used for the LES is Δz ¼ 5 m in the vertical and Δx ¼ 35 m in the horizontal. This corresponds to 2.76 × 10 6 degrees of freedom. The SCM simulations are performed with vertical resolutions of Δz ¼ 5, 20, 50, and 75 m or 300, 75, 30, and 20 degrees of freedom, respectively. In the SCM simulations, the time step is diagnosed from a Courant-Friedrichs-Lewy (CFL) condition based on the maximum updraft velocity in the domain, using a Courant number of 0.9. This results in average time steps of 3, 14, 39, and 63 s, respectively.

Results
The mean profiles obtained with the extended EDMF scheme display very little resolution sensitivity compared to the spread of results from LES, as shown in Figure 4. LES of STBLs are strongly dependent on the discretization numerics and the treatment of subgrid-scale fluxes . Overly diffusive LES models result in excessive cloud top mixing, reducing the water content of the cloud layer and transitioning to decoupled cumulus-like conditions.
Similarly, the ability of SCM simulations to capture the stratocumulus cloud layer is contingent upon the cloud top mixing not being too strong. With large gradients in q t and θ l across the inversion, the mixing length is the main limiter of cloud top diffusive mixing. As shown in Figure 5a, the buoyancy scale (29) is crucial to limit the cloud top eddy diffusivity and maintain a sharp inversion over the mixed layer (see Appendix D for details). It is important to note that in our formulation, the mixing length may be smaller than Δz. This allows to maintain a coupled cloud layer even at coarse vertical resolution.
How the dominant mixing length varies with height in the STBL is shown in Figure 5a. Throughout most of the boundary layer, environmental mixing is determined by the minimum dissipation balance. Mixing is constrained by stratification at cloud top and in the lower part of the cloud, where the environmental cloud fraction f c, 0 is less than unity. The vertically integrated TKE obtained in the SCM simulations is similar across resolutions and follows closely the WENO-based LES statistics, as shown in Figure 5b. Again, the variation of TKE with resolution in the SCM simulations is significantly lower than the spread of values obtained with different LES, not all of which successfully simulate the presence of a stratocumulus cloud layer.
The liquid water path (LWP) time series from the SCM simulations are in agreement with the LES results. At coarse resolution, cloud top entrainment of dry air is too low, which leads to an overestimation of q l and LWP, as shown in Figures 4c and 5c. However, even at this resolution, the water content bias obtained with the EDMF scheme is significantly lower than the dry bias of some of the LES models. The vertical heat and moisture fluxes, as well as the contributions from the turbulent flux (eddy diffusivity) and subdomain-mean terms (mass flux), are shown in Figure 4. The SCM simulations slightly overestimate the heat flux in the cloud layer and underestimate the moisture flux throughout the boundary layer. These biases compensate each other to some extent, leading to a small bias in the buoyancy flux. Similar biases are reported for models using the EDMF scheme and different parameterizations (Wu et al., 2020).
In the extended EDMF scheme, the environmental turbulent flux is the leading contributor to the buoyancy flux. The context of this decomposition should be considered when comparing these results to LES studies of the dynamics governing the STBL (e.g., Davini et al., 2017). Since we do not consider downdrafts in our SCM simulations, the environment contains all dynamic structures of the flow except updrafts. Therefore, the turbulent flux here also represents the transport due to downdrafts. Although LES studies emphasize the importance of convective transport due to downdrafts in stratocumulus clouds (Davini et al., 2017), we find that their implementation is not necessary to reproduce the STBL using the extended EDMF scheme. This is in agreement with Wu et al. (2020), where the authors show that the implementation of downdrafts in an EDMF scheme does not significantly improve simulations of the STBL.
Vertical profiles of TKE and eddy diffusivity are shown in Figure 6. The magnitude of TKE is underestimated by the SCM in the cloud layer, while the SCM maintains similar values to LES in the subcloud layer. The eddy diffusivity is also diagnosed from LES, with the Prandtl number Expression (36) used in our model imposed as a constraint. From the closures (16), and since the simulations are performed in a single column (i.e., the horizontal derivatives in (16) are zero), we can estimate the eddy diffusivity from the environmental shear and buoyancy production as In the diagnosis of (45), updraft and environment identification is performed using the methodology proposed in Couvreux et al. (2010). The eddy diffusivity in the SCM simulations follows a similar profile to K les m below the cloud layer, underestimating K les m at coarser resolution but converging toward K les m as the resolution is refined. The peak in K les m within the cloud layer is due to a positive w ′ 0 b ′ 0 under a vanishing b 0 gradient. This environmental buoyancy production may be attributed to convective downdrafts, which are considered to be part of the environment in our analysis.

Dry Convection
The dry convective boundary layer differs from the previous cases in that the mass flux term is the leading-order contribution to the subgrid-scale vertical transport throughout most of the boundary layer. However, an accurate parameterization of the eddy diffusivity contribution is still necessary for a correct simulation of the dry convective boundary layer.

Simulation Setup
The simulation setup follows Nieuwstadt et al. (1993). The base state is initialized as a mixed layer with potential temperature θ ¼ 300 K up to z 1 ¼ 1,350 m, above which potential temperature increases at a rate of 3 K km −1 . The flow, which is initialized with a horizontal velocity of 1 cm s −1 , is driven by a constant surface heat flux of ⟨w * θ * ⟩ ¼ 6 K cm s −1 .
The simulation is performed in a domain spanning 3.75 km in the vertical. For the LES, the horizontal cross-sectional area is 6.4 × 6.4 km 2 , and the resolution is Δz ¼ 25 m in the vertical and Δx ¼ 50 m in the horizontal. The SCM simulations are performed with vertical resolutions of 25, 50, and 150 m. As for the DYCOMS-II simulations, the time step is diagnosed from a CFL condition. The average time step for these simulations is 14, 30, and 100 s, respectively.

Results
Time-averaged profiles of potential temperature and vertical buoyancy flux are shown in Figure 7 for the fifth hour of simulation. The potential temperature mixed layer and its associated vertical heat flux are well captured for all resolutions considered, with little resolution sensitivity. The convective heat flux is roughly constant throughout the boundary layer, while the diffusive flux decreases with height.
All simulations show a small cold bias throughout the boundary layer and a warm bias below the inversion. The latter is due to the absence of plume overshooting in the SCM simulations, as shown in Figure 7b. The evolution of the boundary layer depth, diagnosed as the height of minimum buoyancy flux (Stevens, 2007), is shown in Figure 7d. The boundary layer growth in the SCM simulations is slower than in LES, with the bias decreasing as the vertical resolution is refined.
Reducing these biases with the extended EDMF scheme is possible, albeit with a different set of parameters controlling the pressure closure (11). These results are not shown here, since the goal of the model is to simulate all boundary layer regimes with a single set of parameters. Learning a set of parameters that minimizes these and other biases in the results shown here and in Cohen et al. (2020) is left for future work.
Finally, Figure 8 shows vertical profiles of TKE and eddy diffusivity. All SCM simulations display a TKE bias above 1,500 m due to the absence of plume overshooting, which leads to zero convective TKE following Equation 5. Eddy diffusion, which is a function of environmental TKE in our model, is not affected by this convective bias. Indeed, the diffusive closure (9) leads to an accurate prediction of the depth of the diffusive layer, as shown in Figure 8b.

Summary and Discussion
The mixing length formulation proposed in this study provides a regime-independent closure of turbulent fluxes for EDMF schemes. The results for the stable boundary layer, stratocumulus-topped boundary layer, and dry convection demonstrate the ability of EDMF schemes with this mixing length closure to accurately describe the structure of the boundary layer in regimes where existing parameterizations currently in use in climate models fail or are inaccurate.
In the stable boundary layer, where convection and the subdomain decomposition in the EDMF scheme do not play a role, the proposed closure is able to reproduce the vertical structure and time evolution of turbulence over a range of vertical resolutions, down to O(10 m). In the stratocumulus-topped boundary layer, where convective fluxes do play a role, the transport owing to environmental diffusion still provides the leading-order contribution to the subgrid-scale vertical fluxes in our EDMF scheme. The way in which environmental stratification limits the mixing length seems to be the crucial feature that allows our EDMF scheme to reproduce the sharp inversion at the stratocumulus cloud top, even at relatively coarse vertical resolution.
Several characteristics differentiate this closure from others proposed in the literature. First, choosing the smooth minimum (40) of various candidate mixing lengths is consistent with the idea that estimates of the mixing length arising from different physical arguments should converge to a similar master length scale if they are simultaneously valid. For widely used expressions such as the harmonic mean (28), this does not hold, leading to unrealistic reductions in mixing. Second, our formulation explicitly links the eddy diffusivity to the effect of convective cells on the environment, leading to a consistently closed TKE balance. This results, for example, in the TKE injection term I appearing in the length scale (18), for which TKE production and dissipation are in balance. Third, the mixing length does not depend on integral quantities such as the boundary layer thickness or Deardorff's convective scale. The inclusion of these terms in other models often leads to regime-dependent nonlocal terms that are noncausal and hence difficult to justify in general. Finally, the closure smoothly connects with Monin-Obukhov similarity theory near the surface with no assumptions about the height at which the transition occurs. This is particularly relevant for climate models with low vertical resolution, for which the use of similarity theory even in the first model level may be inaccurate.
A similar approach to the one shown here may be used to develop increasingly complex closures for high-order turbulence models. As an example, the net dissipation argument used in the TKE production-dissipation inequality (13) can also be applied to the temperature variance budget to diagnose the turbulent Prandtl number. The same could be done for other second-moment budgets in models with additional second-order prognostic equations, to obtain independent diffusivities for different tracers.
Finally, the optimization of the full parameter space was beyond the scope of this study and is left for future work. The access to LES data for a wider range of atmospheric conditions is necessary to enable a more comprehensive optimization of the parameter space in the EDMF scheme.

Journal of Advances in Modeling Earth Systems
The prognostic equation for the convective kinetic energy in subdomain i can be obtained as 10.1029/2020MS002161

Journal of Advances in Modeling Earth Systems
Summing over all subdomains, we obtain the subgrid-scale convective TKE balance where, under the EDMF assumptions, all terms involving within-subdomain covariances are only nonzero in the environment (i ¼ 0). The first and second terms on the right-hand side are turbulent transport terms. The third term represents shear production of convective energy. The fourth and fifth terms yield shear production of subdomain TKE by the convective flow, representing an advective sink in the balance (B3). The sixth and seventh terms are the convective components of the buoyant production and velocity-pressure gradient terms. Finally, the dynamical and turbulent entrainment terms act to transfer subgrid kinetic energy from the plumes to within-subdomain variance. Note that the velocity-pressure gradient term can be rewritten as This yields the definition of the pressure work on the plumes used in Expression (10).
Some of the terms in Budgets (B1) and (B3) transfer subgrid energy among the environment and plumes, resulting in a null contribution to the global budget. The grid-mean TKE prognostic equation that results from their sum is where there is no contribution from pressure-velocity correlations in our model. The evolution of the grid-mean TKE under the net dissipation closure (13) can be obtained by subtracting (13) from (B5): where γ 0 is the net environmental dissipation. According to (B6), grid-mean TKE is generated through convective buoyant production B * and the vertical convergence term S * . Both dynamical and turbulent entrainment act as a transfer term from subgrid-scale convective kinetic energy to environmental TKE, resulting in a grid-mean TKE sink under the net dissipation closure. A schematic of the energetic pathways between budgets (B1) and (B3) and the overall evolution of grid-mean TKE under the mixing length closure presented here is shown in Figure B1.

Appendix C: Probabilistic Model for Cloud Fraction
We consider θ l and q t to be log-normally distributed with expected values θ l; 0 and q t; 0 , variances σ 2 θl and σ 2 q t , and covariance σ 2 q t ; θl . The log-normal distribution is preferred over the commonly used Gaussian distribution (e.g., Sommeria & Deardorff, 1977) for two reasons: Both θ l and q t remain non-negative, and positive skewness is allowed. Under the Gaussian assumption, negative values of q t may be drawn from the distribution if σ 2 q t =q 2 t; 0 is not sufficiently small (Mellor, 1977). In addition, distributions with positive skewness have been shown to capture the development of cumulus convection better (Bougeault, 1981).
The expected value of cloud fraction f c, 0 can be computed as (Mellor, 1977) where H(·) is the Heaviside function and p(θ l , q t ) is the probability density function (PDF) of the log-normal bivariate distribution with marginal PDFs given by q t ∼ LNðμ q t ; s 2 q t Þ; and θ l ∼ LNðμ θl ; s 2 θl Þ; μ θl ¼ ln θ where w j and w i are the Gauss-Hermite weights corresponding to evaluation points θ l, j and q t, i , respectively. The evaluation points (θ l, j , q t, i ) of the log-normal distributions (C2) and (C4) are related to the Gauss-Hermite mass points (ξ j , χ i ) through the normal distributions x and y with same parameters: and q t; i ¼ e y i ; y i ¼ μ q t þ ffiffi ffi 2 p s q t χ i ; y ∼ Nðμ q t ; s 2 q t Þ: Note that the evaluation points θ l, j are drawn from the conditional PDF (C4). In (C6), the liquid water specific humidity q l is obtained as q l = q t, i − q s (θ l, j , q t, i ) for q t, i > q s (θ l, j , q t, i ), where q s is the equilibrium saturation specific humidity. Thus, supersaturation is not considered and all excess water vapor is immediately converted to liquid water condensate. The equilibrium saturation specific humidity is found iteratively using a saturation adjustment procedure (see Bryan & Fritsch, 2002, for details). Consistent with this approach, the environmental liquid water specific humidity q l; 0 is computed as q l; 0 ¼ 1 π ∑ ni i w i ∑ nj j w j ½q t; i − q s ðθ l; j ; q t; i ÞHðq t; i − q s ðθ l; j ; q t; i ÞÞ: In this study, the probabilistic cloud model is implemented using n i ¼ n j ¼ 3.

Journal of Advances in Modeling Earth Systems
The effect of the entrainment term I in the formulation of l tke is noticeable in the buoyancy flux and LWP of the STBL, both of which show a positive bias with respect to LES in Figure D2. The effect of the entrainment term I is more significant in the dry convective case, where its absence leads to very low diffusive mixing in the boundary layer and low TKE compared to LES. Finally, the results for NLW show that the implied balance in l tke can approximate the wall constraints relatively well, leading only to small biases in the SBL simulation.

Data Availability Statement
(a) (b) (c) (d) Figure D3. Profiles of (a) potential temperature, (b) buoyancy flux, (c) normalized TKE, and (d) eddy diffusivity for the fifth hour of the dry convection simulation. Results shown for LES and for the SCM (Δz ¼ 150 m) with alternative mixing length formulations.

Journal of Advances in Modeling Earth Systems
LOPEZ-GOMEZ ET AL. 26 of 28