An Extended Eddy-Diffusivity Mass-Flux Scheme for Unified Representation of Subgrid-Scale Turbulence and Convection.

Large-scale weather forecasting and climate models are beginning to reach horizontal resolutions of kilometers, at which common assumptions made in existing parameterization schemes of subgrid-scale turbulence and convection-such as that they adjust instantaneously to changes in resolved-scale dynamics-cease to be justifiable. Additionally, the common practice of representing boundary-layer turbulence, shallow convection, and deep convection by discontinuously different parameterizations schemes, each with its own set of parameters, has contributed to the proliferation of adjustable parameters in large-scale models. Here we lay the theoretical foundations for an extended eddy-diffusivity mass-flux (EDMF) scheme that has explicit time-dependence and memory of subgrid-scale variables and is designed to represent all subgrid-scale turbulence and convection, from boundary layer dynamics to deep convection, in a unified manner. Coherent up and downdrafts in the scheme are represented as prognostic plumes that interact with their environment and potentially with each other through entrainment and detrainment. The more isotropic turbulence in their environment is represented through diffusive fluxes, with diffusivities obtained from a turbulence kinetic energy budget that consistently partitions turbulence kinetic energy between plumes and environment. The cross-sectional area of up and downdrafts satisfies a prognostic continuity equation, which allows the plumes to cover variable and arbitrarily large fractions of a large-scale grid box and to have life cycles governed by their own internal dynamics. Relatively simple preliminary proposals for closure parameters are presented and are shown to lead to a successful simulation of shallow convection, including a time-dependent life cycle.

Introduction

The representation of clouds and convection in weather forecasting and climate models continues to be at the core of biases and uncertainties in weather forecasts and climate change projections (e.g., Bauer et al., 2015; Bechtold et al., 2004; Bony & Dufresne, 2005; Bony et al., 2015; Brient et al., 2016; Brient & Schneider, 2016; Cess et al., 1990, 1996; Dufresne & Bony, 2008; Klocke & Rodwell, 2014; Stephens et al., 2010; Vial et al., 2013; Webb et al., 2006). The dynamics controlling clouds and convection occur at scales smaller than the resolutions of current general circulation models (GCMs). Although GCM resolutions are being refined at a rate that will make explicitly resolving deep convection routine within the next decade (Kajikawa et al., 2016; Miyamoto et al., 2013; Ohno et al., 2016; Palmer, 2014; Prein et al., 2015), shallow convection will still need to be parameterized for the foreseeable future (Schneider et al., 2017). Thus, we will continue to need subgrid‐scale (SGS) schemes for clouds and convection, and ideally SGS schemes that are resolution‐adaptive, so that they can adjust as GCM resolutions increase and the need for parameterization of deep convection diminishes (Arakawa & Jung, 2011; Wyngaard, 2004).

However, currently used SGS schemes make a seamless representation of the dynamics controlling clouds and convection at different scales difficult. GCMs usually contain several disparate SGS schemes, for example, separate schemes for boundary layer turbulence, shallow convection, and deep convection (e.g., Donner et al., 2011; Neale et al., 2010). Turbulent and convective dynamics that in reality form a continuous spectrum are represented discontinuously by the separate SGS schemes. As a result, convergence of the schemes as the resolution increases (e.g., to SGS closures for large‐eddy simulations) usually does not occur. The different schemes also converge to different physical limits, for example, as the latent heat of vaporization approaches zero. In that latter case, all schemes should represent dry turbulence, but they do so in very different ways, for example, often through diffusive closures in boundary layer turbulence schemes and mass flux closures in convection schemes (McFarlane, 2011; Zhang et al., 2013). This leads, in practice, to competing effects of the different schemes on the grid‐scale (GS) dynamics (Brient et al., 2016; Vial et al., 2016, 2017; Zhang et al., 2013). Because the schemes only interact through their GS effects, they usually do not adequately capture interactions and transitions among different cloud and convective regimes (e.g., D'Andrea et al., 2014; Hohenegger & Bretherton, 2011; Kuang & Bretherton, 2006; Lappen & Randall, 2001a). Moreover, the plethora of different SGS schemes affecting clouds and convection contributes to the proliferation of adjustable parameters in GCMs, increasing model uncertainty and making it a challenging problem to tune models (e.g., Flato et al., 2013; Golaz et al., 2013; Hourdin et al., 2013, 2017; Mauritsen et al., 2012; Randall & Wielicki, 1997; Ruiz et al., 2013; Schirber et al., 2013).

Unified representations of all convective and cloud dynamics have the potential to ameliorate some of these problems, by representing all SGS dynamics—be they boundary layer turbulence, shallow convection, or deep convection—in a unified framework without artificial discontinuities. Several unified SGS schemes have been proposed. Lappen and Randall (2001a, 2001b, 2001c) present a scheme that combines into a single unified scheme ideas from mass flux closures, commonly used for cumulus convection, and from higher‐order closures, sometimes used for boundary layer turbulence. The unified scheme is based on assuming joint probability density functions (PDFs) for the relevant state variables and on deriving equations for the statistics characterizing the joint PDFs. Larson et al. (2002), Golaz et al. (2002), and Larson and Golaz (2005) developed these ideas further into the Cloud Layers Unified by Binormals (CLUBB) scheme, which assumes binormal PDFs jointly for vertical velocities and conserved scalars. First and second‐order moment equations with closure assumptions are solved for the dynamics of the statistics characterizing the PDFs. This scheme has been successfully used across scales and across low‐cloud regimes (Guo et al., 2015; Larson et al., 2012). However, the PDFs are local in the vertical; vertical correlations such as those that arise in vertically coherent updrafts and downdrafts are not taken explicitly into account. Therefore, CLUBB has usually not been used, or has been used with only limited success (Guo et al., 2015), to parameterize deep convection, for which vertically coherent updrafts and downdrafts are crucial. By contrast, the Unified Convection Scheme (UNICON) by Park (2014a) explicitly parameterizes the dynamics of updrafts, downdrafts, mesoscale organized flows, and their interactions. It helps to reduce some of the convection‐related biases in GCMs (Park, 2014b). But it still requires a separate scheme for boundary layer turbulence. Its diagnostic formulation of updrafts and downdrafts makes its application at high resolutions questionable, when GS dynamics evolve on similar timescales as SGS convection and so cannot be presumed to be instantaneously in equilibrium with the GS dynamics. Gentine et al. (2013a, 2013b) and D'Andrea et al. (2014) developed a probabilistic plume model that unified the representation of dry, shallow, and deep convection in a simplified bulk model, which can also be regarded as a prototypical unified convection scheme.

Here we lay the theoretical foundations for an alternative unified SGS scheme, which closes first and second‐order moment equations for all turbulent and convective motions and which can explicitly represent nonlocal updrafts and downdrafts. The scheme is based on the eddy‐diffusivity mass‐flux (EDMF) approach (Angevine, 2005; Angevine et al., 2010; Siebesma et al., 2007; Siebesma and Teixeira, 2000; Soares et al., 2004; Suselj et al., 2012; Neggers, 2009; Neggers et al., 2009). The EDMF approach decomposes all SGS motions into two components: an approximately isotropic turbulent environment, and vertically coherent convective structures. The approximately isotropic turbulent environment is represented by down‐gradient SGS fluxes, with an eddy diffusivity (ED) that can depend, for example, on the environmental turbulence kinetic energy (TKE). The vertically coherent convective structures, such as dry thermals or cumulus updrafts, are represented as entraining and detraining plumes that satisfy a mass flux (MF) equation. A single bulk plume is usually used to represent the effect of an ensemble of updrafts, but multiple plumes have also been considered (Cheinet, 2003, 2004; Neggers, 2012; Suselj et al., 2013). Such schemes can be coupled to a probabilistic cloud scheme (e.g., Sommeria & Deardorff, 1977) to form a complete SGS closure. EDMF schemes have been successfully implemented to represent boundary layer turbulence and shallow convection in several numerical weather prediction models (Han et al., 2016; Köhler, 2005; Köhler et al., 2011; Suselj et al., 2014).

We extend the EDMF approach by:

Making plumes prognostic to represent the life cycle of updrafts and downdrafts;

Extending the plume formulation to include downdrafts;

Allowing cross‐sectional areas of plumes to be variable and cover arbitrary fractions of a grid box footprint; and

Consistently partitioning second‐order moments (e.g., TKE and scalar variances) between plumes and environment.

These extensions make the new scheme more flexible for representing convection at different resolutions and lead to a unified and physically consistent scheme. While this paper focuses on the theoretical framework, we also present preliminary proposals for the closure parameters and parametric functions that arise. We show that the resulting extended EDMF scheme captures both the average condition and the transient evolution of shallow convective clouds.

Section 2 presents the concepts of the environment‐plume decomposition and the key equations of the extended EDMF scheme. Section 3 lists a set of preliminary closure formulations (e.g., entrainment and detrainment rates) for the EDMF equations, which we use for the first tests. Section 4 presents an evaluation of the extended EDMF scheme against large‐eddy simulations (LES) of low clouds. Section 5 summarizes the results and presents an outlook of further developments we envision for the closure formulations.

Extended EDMF Scheme

The EDMF scheme is based on a decomposition of the horizontal domain into n + 1 subdomains, indexed by (Yano, 2014a). The total domain can be thought of as the horizontal plane in a grid box of a large‐scale model. Among the subdomains, the environment (i = 0) is unique in that we assume it to contain quasi‐isotropic turbulence, so that turbulent dispersion therein can be well represented by local eddy diffusion. The other subdomains ( ) represent vertically coherent plumes, which can be updrafts or downdrafts. They are assumed to consist of nonlocal mass fluxes that interact with the environment and potentially with each other through entrainment and detrainment. However, individual plumes are assumed to be horizontally homogeneous. Therefore, from a probabilistic point of view, the EDMF scheme models PDFs of scalars as superpositions of an environmental PDF (approximately Gaussian) and delta distributions for the plumes (Lappen & Randall, 2001a).

We begin by deriving the dynamic equations for each subdomain in general terms, without assuming any specific form of PDFs. The distinction between the environment (i = 0) and plumes ( ) will be made explicit thereafter, when we apply the EDMF assumptions to reduce the number of terms that require closure. The appendix summarizes the notations and symbols.

Domain Decomposition

Each subdomain is characterized by the area fraction of the horizontal domain it occupies at a given time t and level z, so that

We denote by the horizontal mean of variable over subdomain i, and by fluctuations about the subdomain mean. We denote by the horizontal mean of over the total domain, and by fluctuations about the total mean. The difference between subdomain and domain means is denoted as . Here, the flow field can be a scalar such as liquid water potential temperature θ or total water specific humidity q, or a velocity component (u, v, w).

Mean fields and covariances over the total domain can be decomposed as where ψ is a second flow field. The second equation (3) is a decomposition of the total covariance into covariances among subdomain means and (first term on the right‐hand side) and covariances within subdomains (second term). From the large‐scale model perspective, represents the resolved GS mean, and represents the SGS fluxes and (co‐)variances of scalars that need to be parameterized.

Dynamic Equations for Subdomains

In deriving dynamic equations for mean fields and covariances in the subdomains, we make the following simplifying assumptions:

Horizontal variations of density ρ are neglected, except in the calculation of vertical accelerations. This makes the EDMF scheme similar to a subdomain‐averaged anelastic system, and area‐weighted averages over subdomains as in equations (2) and (3) are equivalent to mass‐weighted averages.

Horizontal variations of SGS statistics (mean fields and covariances) are neglected, so that only derivatives with respect to time t and height z appear (boundary‐layer approximation).

Mean horizontal velocities in any subdomain are taken to be equal to the domain‐mean values , so that only advection by domain‐mean horizontal velocities contributes to SGS horizontal fluxes.

Fluid masses exchanged between any two subdomains by entrainment or detrainment carry with them the mean properties of the subdomains (mean‐field approximation). This also applies to exchange of covariances among subdomains: they are entrained or detrained like other fluid properties.

With these assumptions, the continuity equation for the area fraction a becomes the equation for the subdomain mean becomes and the equation for the within‐subdomain covariance becomes

Here, is the del operator in the horizontal plane. The ‐terms are included to allow for the horizontal advection of SGS properties across grid cells. The fractional entrainment rate gives the rate of entrainment into subdomain i from subdomain j, defined so that , where E is the mass entrained per unit time into subdomain i from j (normalized by the area of the entire domain). The fractional detrainment rate δ gives the rate of detrainment from subdomain i into all other subdomains, defined so that , where Δ is the mass detrained from subdomain i. (Into which subdomain the mass is detrained does not matter for the subdomain i from which it is detrained. Hence, the subscript j only appears in the entrainment rate for subdomain i, because the properties of the air entrained from subdomain j matter for i.) By mass conservation, any mass detrained from subdomain j must be entrained by other subdomains (or re‐entrained by j), so that , and thus

Exact definitions of entrainment and detrainment rates have been given, e.g., by de Rooy et al. (2013) and Yano (2014a). They are reproduced with slight modifications in Appendix A for reference. A detailed derivation of the covariance equation (6) is given in Appendix B.

The three terms on the right‐hand side of the mean equation (5) represent the vertical turbulent transport within a subdomain, entrainment, and detrainment to and from a subdomain, and sources and sinks within a subdomain. The source term may include pressure gradients for the velocity components (u, v, w), buoyancy acceleration for the vertical velocity w, and precipitation and evaporation for thermodynamic variables and scalars such as liquid water potential temperature θ or total water specific humidity q. It may also include radiative heating or cooling for thermodynamic variables.

The right‐hand side of the covariance equation (6) includes terms representing covariance generation through down‐gradient turbulent transport within a subdomain (or covariance destruction through upgradient turbulent transport), entrainment and detrainment of covariance to and from a subdomain and conversion of mean differences between subdomains to covariances within a subdomain, vertical transport of covariance by turbulence within a subdomain (nonlinear triad interactions), and sources and sinks (including molecular diffusion) that covary with and ψ.

The dynamic equations for the total domain follow by summing the subdomain equations (4), (5), (6) over all indices i and eliminating terms with the domain decomposition equations (1), (2), (3) and (7); details are given in Appendix C.

Vertical Velocity and TKE Equations

There are two special cases of the mean equation (5) and of the covariance equation (6) that we will use extensively and that are worth stating explicitly. The first is the mean equation for the vertical velocity , which is

We further define a domain‐mean hydrostatic pressure . By adding the relation to the right hand side of equation (8), we get a more familiar form of the vertical velocity equation:

Here, is the buoyancy anomaly with respect to the domain mean, where θ is the virtual potential temperature; is the pressure anomaly with respect to . Equation (9) holds generally, even for fully compressible GS equations. Note that by definition, but if the GS motion is nonhydrostatic, the domain‐mean pressure acceleration generally is nonzero. In the special case when the GS equations are anelastic, we may replace b with the anelastic buoyancy and replace with the anelastic dynamic pressure, and equation (9) still holds.

The second special case is the covariance equation (6) for TKE. The domain‐mean TKE can be decomposed as where , and we have assumed that the subdomain‐means of horizontal velocities are the same for all subdomains (assumption 3 in section 2.2), so that and for all i. After decomposing the source/sink term into three components representing the effects of buoyancy, pressure, and viscous dissipation, the TKE equation for subdomain i becomes

The TKE equation (11) for each subdomain represents a consistent partitioning of TKE across the subdomains, which is important, for example, for using TKE in diffusive closures of turbulent fluxes within a subdomain, as is done in the EDMF scheme we are developing here.

The pressure terms in both equations include contributions from the normal and form stresses, while the horizontal stresses and the dilatation effect also contribute to the pressure term in the TKE equation. Detailed derivations are presented in Appendix D.

EDMF Assumptions

The dynamic equations (4), (5), (6), (9), and (11) hold quite generally, under mild assumptions that allowed us to carry out a domain decomposition and drop horizontal derivatives of subdomain statistics. However, the dynamic equations contain many terms that require closure, from entrainment rates among all subdomains, detrainment rates for all subdomains, to turbulent fluxes within subdomains. The EDMF scheme makes additional assumptions that reduce the number of terms requiring closure (Neggers, 2009; Neggers et al., 2009; Siebesma & Teixeira, 2000; Siebesma et al., 2007; Soares et al., 2004; Suselj et al., 2012):

There is one distinguished subdomain, the environment (i = 0), which has turbulent fluctuations , etc. In all other subdomains, the plumes ( ), turbulent fluctuations , etc. are assumed to be negligible.

Turbulent fluxes in the environment, such as , are assumed to be diffusive, so that, for example,

with an eddy diffusivity K. Turbulent fluxes of (co‐)variances such as are also closed diffusively in an analogous fashion.

In addition, EDMF closures usually assume that entrainment and detrainment only occur between plumes ( ) and the environment (i = 0) (e.g., Yano, 2014a). However, we do not make this assumption here, to allow, for example, for updrafts to detrain air into neighboring downdrafts.

With these assumptions, we are left with five principal equations in addition to the continuity equation (4) for the area fraction a: equations for vertical velocity and scalars in the plumes, and equations for scalars , covariances , and TKE in the environment. The equation for the vertical velocity in the plumes ( ) becomes and that for the scalar in the plumes,

The equation for the vertical velocity in the environment becomes that for the environmental scalar becomes and that for environmental covariances becomes

Finally, the environmental TKE equation is

Note that the diffusive closure is not applied to the buoyancy production term in the TKE equation: the discontinuities of and at saturation make it necessary to compute the buoyancy production separately for cloudy (saturated) and clear (unsaturated) regions, as we will discuss further in section 3.6.

The pressure terms in the environmental TKE equation (19) are related to the pressure terms in the vertical velocity equations (13) and (15): the TKE production and destruction by pressure perturbations should sum to zero over the domain, except for a term involving the normal stress, , which is usually absorbed into the eddy diffusion of TKE, and a dilatation term (see Appendix D). This means, for example, that kinetic energy lost in plumes by pressure drag exerted by their environment should appear as kinetic energy in the environment, to the extent it is not dissipated by diffusion. Detailed closure formulations are discussed in section 3.6.

The equation for the total domain are recovered by summing over all indices i, as discussed in Appendix C. Specifically, the equation for a scalar is

The first and second terms on the right hand side represent the MF and ED components of vertical SGS fluxes, respectively, with mass flux . We assume that mass fluxes in plumes and in the environment compensate each other within a grid box, so that . This assumption of compensating subsidence occurring in the same grid box as the updrafts it compensates may need to be relaxed at high GS resolution, when compensating subsidence may occur in neighboring grid boxes.

For horizontal velocity , we only need to solve for its domain mean, because of the assumption of horizontal homogeneity across subdomains. The prognostic equation for in the approximation of the atmosphere as a thin shell is where f is the Coriolis parameter, and k is the upward pointing vertical unit vector. Note that the convective momentum transport term (i.e., the momentum transport by plumes) vanishes because by assumption. To consistently represent convective momentum transport, it would be necessary to relax the homogeneity assumption for horizontal velocities.

Alternatively, the EDMF equations in flux form (13)–(18) can also be written in advective form using the continuity equation (4). For example, the advective form of the scalar equation (14) for the plumes is

The detrainment rate no longer appears explicitly because only entrainment of fluid with scalar values different from leads to material tendencies of the scalar in subdomain i; detrainment affects mass fluxes in subdomain i but not material scalar tendencies (e.g., Randall, 2015, chapter 6). For and with the source term S including buoyancy and pressure gradient accelerations, the scalar equation (21) becomes an advective equation for the vertical velocity. When mean SGS fields adjust essentially instantaneously to GS conditions, so that they are in quasi‐equilibrium with the GS (Arakawa & Schubert, 1974), the explicit time derivative can be neglected, and the scalar equation (21) reduces to the familiar equation for the rate of change of the scalar with height z used in many convection closures (e.g., Arakawa & Schubert, 1974; Bretherton et al., 2004; Gregory, 1997; Yano, 2014a). However, for numerical discretization, the equations in flux form are often preferable because conservative discretization schemes can be more easily designed on their basis.

Summary of Extensions to Standard EDMF Schemes

The EDMF equations (13), (14), (15), (16), (17), (18) differ from standard EDMF equations (Neggers, 2009; Neggers et al., 2009; Siebesma et al., 2007; Siebesma & Teixeira, 2000; Soares et al., 2004), principally in five ways:

Explicit time derivatives appear in the equations, which enables SGS memory and explicit life cycles of convective plumes. Such memory terms have generally been neglected in parameterization schemes (Yano, 2014a), although there is evidence for their importance, e.g., for deep convection (e.g., Mapes, 1997). There are also fundamental reasons, rooted in response theory (Ruelle, 2009), why memory terms should arise in parameterizations if the timescale separation between parameterized and resolved processes is finite, as it is for convective clouds, especially at high GS resolutions (Lucarini et al., 2014; Wouters et al., 2016; Wouters & Lucarini, 2013).

The area fractions a are not necessarily small and vary with height z and time t, rather than being constant and small for plumes. Variable area fractions that can become large are more appropriate at high GS resolutions (Arakawa & Jung, 2011) and allow greater flexibility in representing entrainment and detrainment while considering the buoyancy accelerations in plumes.

All vertical advection and entrainment/detrainment terms are written consistently in flux form, which facilitates the design of conservative numerical schemes for them.

Entrainment and detrainment are allowed among plumes, not just between plumes and the environment, which makes it possible to include downdrafts that entrain air from adjacent updrafts.

The environmental TKE equation (18) couples the environmental TKE in a consistent fashion with the dynamics of plumes, making it possible, for examples, to consider transport of TKE in plumes and subsequent detrainment into the environment. This enables nonlocal vertical redistribution of TKE, which is important, for example, in convective boundary layers and deep convection (e.g., Khairoutdinov & Randall, 2002; Shin et al., 2013; Witek et al., 2011a).

Requirements for Convergence at Increasing Resolution

Resolutions of atmospheric models are continuously increasing. Resolutions in the gray zone for deep convection (Wyngaard, 2004), with horizontal scales of order kilometers, at which deep convective systems begin to be resolved, will become routine in the coming years (e.g., Ban et al., 2015; Miyamoto et al., 2013; Palmer, 2014; Schneider et al., 2017). SGS parameterizations should be able to adapt to this increasing resolution. As turbulent eddies and convective plumes are increasingly resolved by the GS model, the SGS fluxes, and covariances parameterized by the EDMF scheme should gradually diminish. At LES resolutions, the EDMF scheme should act like an SGS closure for an LES (if its eddy diffusion is extended to be isotropic). This can be achieved if the following conditions are satisfied:

At sufficiently high resolution, the MF component should disappear because, as the convective plumes are resolved on the GS, the SGS motion becomes more isotropic and coherent SGS plumes contribute less to vertical fluxes. This is satisfied if (i) scalar and buoyancy anomalies in SGS plumes become small, so that SGS plume velocities become small, or (ii) mixing rates (entrainment and detrainment) become large, so that SGS plumes are locally mixed out. These two conditions are related: a larger mixing rate leads to a smaller buoyancy anomaly and thus a smaller vertical velocity; this conversely lengthens the time over which air masses in plumes can mix and hence increases the mixing rate (e.g., Romps, 2016; Romps & Kuang, 2010a, 2010b). To achieve this, the formulations of initial and boundary conditions for scalars and buoyancy in subdomains, mixing rates, and/or the pressure drag for vertical velocity, should all depend on resolution.

At intermediate resolution, the number of plumes in a grid box can fluctuate strongly over time because of small‐sample variability. Therefore, some stochasticity should be included in the MF equations (Berner et al., 2017), for example, stochasticity in the number of plumes in the domain, in their initial or boundary conditions (e.g., Plant & Craig, 2008), or in entrainment and detrainment rates (e.g., Nie & Kuang, 2012; Romps, 2016; Romps & Kuang, 2010a; Suselj et al., 2013).

The eddy diffusivity should decrease as resolution increases, and it should asymptote to an SGS diffusivity used in LES. For a TKE‐based closure for the diffusivity, , the mixing length l should not exceed the grid size h as h becomes small (Deardorff, 1980; Moeng, 1984). However, one remaining difference between the EDMF scheme at high resolution and common SGS closures in LES is that the EDMF scheme is not isotropic: eddy diffusion occurs only in the vertical direction. An isotropic extension of the eddy diffusion in EDMF and fully three‐dimensional extensions of the TKE budget are thus required for convergence at LES‐like resolution.

At LES‐like resolution, the SGS TKE and scalar covariances should scale with the grid size h according to similarity theory for the inertial subrange. For three‐dimensional homogeneous and isotropic turbulence, and (based on an integration over the Kolmogorov spectrum of TKE). Since the SGS TKE and covariance dissipation rates are constrained by their generation rates at resolved scales, which do not necessarily change with resolution, these scaling laws should hold as long as the TKE and covariance dissipation rates scale as and , with a mixing length (see section 3.6).

Closure Formulations

To close the EDMF equations (4) and (13), (14), (15), (16), (17), (18), we need to parameterize the entrainment rates among subdomains, the detrainment rates δ, and various source terms . For the environment, we additionally need to parameterize the TKE dissipation term and the eddy diffusivity K. We also need to formulate the boundary conditions for the plumes and the environment. Here we discuss relatively simple, preliminary closures, and their numerical implementation, for now focusing on a single bulk updraft (with index 1). Multiple updrafts and downdrafts, which can represent variations among plume properties, closures for microphysics, radiation, and environmental scalar covariances, as well as possible extensions to represent convective momentum transport, will be discussed in a later paper. Table 1 summarizes the closure parameters that arise in what follows.

Table 1

Closure Parameters

Symbol	Description	Value
a_1,ll	Updraft fraction at the lowest atmospheric level	0.10
(al, bl)	Constants for the surface‐based mixing length formula	(−100,0.2) when Λ<0
		(2.7,−1) when Λ≥0
αb	Scaling constant for virtual mass term	1/3
αd	Scaling constant for drag term	0.375
c_ϵ	Scaling constant for entrainment rate	0.12
c_δ,0	Scaling constant for detrainment rate	0.12
ce	Scaling constant for TKE dissipation	2.0
cK	Scaling constant for eddy diffusivity and eddy viscosity	0.1
δB	Cloud layer background detrainment rate	0.004  m^−1
rd	Length scale of the horizontal spacing between plumes	500 m

Boundary Conditions

No‐penetration conditions (i.e., ) are enforced at what we assume to be a flat surface in equations (14), (15), (16), (17), (18), (19), (20). (In the presence of topography, the vertical velocities at the surface would not generally vanish, but the domain‐mean vertical velocity at the surface, , would need to be partitioned across the subdomains.) At the top of the model domain, we assume that the atmospheric dynamics are predominantly large‐scale and resolved on the GS. Therefore, the turbulent fluxes, TKE, and scalar variances are all set to zero, and no‐penetration conditions are enforced.

The domain‐mean turbulent fluxes at the surface (subscript sfc) are usually either prescribed (as in our test cases in section 4) or can be computed from bulk aerodynamic formulas and Monin‐Obukhov similarity theory (Businger et al., 1971; Byun, 1990).

The domain‐mean TKE and scalar covariances at the lowest atmospheric level (subscript ll) are prescribed as and where is the height of the lowest level above the surface, Λ is the Monin‐Obukhov length, is the friction velocity scale, (or ) is the scalar scale satisfying , and is the convective velocity scale for unstable boundary layers. These length and velocity scales are computed from the surface fluxes based on Monin‐Obukhov similarity theory, with depending on the dry convective boundary layer depth , following Wyngaard and Coté (1974) and Wyngaard (1975). The choice of coefficients in (22) and (23) is based on field measurements (see also Potty et al., 1997 and Witek et al., 2011b). The closures (23) and (24) may need to be modified at high GS resolution, when TKE and scalar covariances become resolved at the grid scale.

The TKE and scalar covariances serve as lower boundary conditions for the environment equations (17) and (18), after applying the subdomain decomposition given by equation (3) to exclude the covariances caused by differences between plume and environmental means. We further assume that the turbulent fluxes at the surface go into the environment, so that they provide the lower boundary conditions for the ED term in equation (16). This is consistent with the no‐penetration condition of the updraft ( ), which implies that the updraft does not contribute to surface turbulent fluxes.

At the lowest atmospheric level, we for now assume that the updraft fraction is fixed and that the updraft air carries the mean properties of air in the top fraction of the buoyancy distribution. If buoyancy fluctuations near the surface are linearly related to fluctuations in θ and q, and given that for positive sensible and latent heat fluxes at the surface, θ and q are perfectly positively correlated according to the surface covariance condition (23), the updraft air also consists of the top fraction of the θ and q distributions at the lowest level. For Gaussian distributions, the mean values over this top fraction of the distributions can be written as where represents the mean of the top a fraction of a standard normal distribution (Neggers et al., 2009), with the analytical form where is the inverse cumulative distribution function of the standard Gaussian distribution evaluated at the ‐th quantile. The prognostic EDMF is able to run with as large as 0.99. For shallow cumulus (e.g., the BOMEX case in section 4), we find that the vertical fluxes are insensitive to the choice of because changes in updraft area fraction and vertical velocity largely compensate; however, the cumulus cloud fraction is somewhat sensitive to and increases with . For all following tests, we set , following previous EDMF studies (e.g., Neggers, 2009; Soares et al., 2004; Witek et al., 2011a), so that .

Entrainment and Detrainment

The entrainment and detrainment rates in the plume continuity equation (4) are assumed to consist of two distinct processes (de Rooy et al., 2013): turbulent entrainment and detrainment, primarily across the lateral edges of plumes, and dynamical entrainment and detrainment, primarily at the top and bottom boundaries of plumes.

Lateral turbulent entrainment and detrainment involve disorganized inflows and outflows across plume boundaries. For now, we prescribe the fractional entrainment and detrainment rates through simple relations that depend on updraft buoyancy and vertical velocity, adapting an expression from Gregory (2001) (see also de Roode et al., 2012). The fractional entrainment and detrainment rates take the form and where is the Heaviside step function, and the “background” cloud layer detrainment rate is , which is consistent with the value diagnosed by Siebesma and Cuijpers (1995). A motivation for the functional form of the entrainment and detrainment rate dependence on updraft buoyancy and velocity is that is of dimension 1/length, i.e., of the dimension of entrainment and detrainment rates, making and nondimensional coefficients or functions. It is also physical for slower and more buoyant updrafts to entrain more per unit length. We choose based on LES of shallow convection, which is roughly consistent with the value of 0.0833 used by Gregory (2001); a smaller may be necessary for deep convection. We further choose , with and with a function that enforces and through dynamical detrainment. More sophisticated formulations based on buoyancy sorting (Bretherton et al., 2004; Kain & Fritsch, 1990) and with stochasticity (e.g., Nie & Kuang, 2012; Romps, 2016; Romps & Kuang, 2010a; Suselj et al., 2013) may prove to be preferable in the future.

At the top and bottom of plumes, buoyant acceleration and deceleration produce dynamical entrainment and detrainment through massive inflows and outflows. For updrafts, the dynamical entrainment is assumed to occur only at the lowest atmospheric level, and it is implicitly realized by resampling the updraft boundary conditions from the domain distribution at every time step (section 3.1). The dynamical detrainment is represented through the function that modulates the coefficient in (27), to ensure that the updraft area fraction does not exceed one. In our test cases in this paper, this does not occur, even without an explicit dynamical detrainment to prevent it. So the function Γ is taken to be zero in our test cases. To limit updraft area fractions more generally, one may take Γ to be a function that rapidly increases as the environmental area fraction a 0 becomes small, e.g., a logistic function with suitable coefficients k > 1 and . This would ensure that dynamical detrainment is active when a 0 becomes small, and it would detrain decelerating updrafts completely as because of the scaling of the detrainment rate.

The above discussions in principle also apply to precipitative downdrafts, with two major differences. First, precipitative downdrafts are caused by precipitate falling into and evaporating in an unsaturated environment, generating negative buoyancy. Thus, the dynamical entrainment for downdrafts is closely coupled to microphysics, involves updraft‐downdraft interactions, and may occur over a range of heights. Second, downdrafts must terminate at the surface even if they are still negatively buoyant. Thus, the dynamical detrainment is concentrated near the surface, where it has strong impacts on updrafts and the environment. Note that with our current definitions, the entrainment and detrainment rates for downdrafts are negative because mass is entrained or detrained over negative vertical displacements. The parameterization of downdrafts is left for future study.

Eddy Diffusivity

As in common diffusive closures, the environmental eddy diffusivity K is assumed to depend on the environmental TKE through where is a scaling constant, and l is the eddy mixing length (e.g., Witek et al., 2011b). The eddy mixing length is a Lagrangian decorrelation length, and it generally depends on the quantity being mixed (Galperin et al., 1988; Grenier & Bretherton, 2001; Mellor & Yamada, 1974, 1982). That is, eddy viscosities and eddy diffusivities for different quantities usually differ. For simplicity, we use one mixing length for all quantities, including scalars, their covariances, and TKE, prescribing it following Witek et al. (2011b) as

Here, the length scale l captures the limitation of the vertical extent of eddies by their distance from the surface and by surface‐layer stability, where Λ is again the Monin‐Obukhov length and is the von Kármán constant. The coefficients are for unstable surface layers ( ) and for neutral or stable surface layers ( ) (Nakanishi, 2001). The other length scale l relates the eddy size to the eddy turnover timescale τ by Taylor's hypothesis, with (Teixeira & Cheinet, 2004). Our choice of mixing length l is broadly consistent with previous studies (e.g., Grenier & Bretherton, 2001; Witek et al., 2011b), with the important distinction that our ED closure is restricted to the environment, at the exclusion of convective plumes. We additionally limit the mixing length l to be no less than the vertical mesh size h.

Momentum Source Terms

The horizontal momentum equation (20) contains three source terms: the first term is closed with the ED closure, and the two other terms (GS pressure gradient and Coriolis acceleration) are computed by the GS model. Therefore, no additional closure is needed in this equation. However, in the vertical momentum equation (9) for the plumes and the environment, the pressure term needs to be parameterized. We can separate it into the GS pressure gradient that is computed by the GS model (and contributes to GS acceleration), and a SGS pressure perturbation term.

The SGS pressure perturbation term for a plume is usually formulated as the sum of a virtual mass term and a drag term (e.g., Bretherton et al., 2004), both of which generally reduce the buoyant acceleration of the plume. The virtual mass term represents the reduction of a plume's effective buoyancy due to its mechanical pushing on and pulling of the environment (Simpson & Wiggert, 1969). As a rough approximation, we assume that it is proportional to the buoyancy acceleration, as it is in simple geometries such as a buoyant sphere (Odar & Hamilton, 1964; Romps & Kuang, 2010b; Turner, 1963). We further assume that the drag only occurs between the plume and the environment (no plume‐plume interactions), and that it is quadratic in updraft velocity (Romps & Charn, 2015; Simpson & Wiggert, 1969). For dimensional reasons, it also needs to be inversely proportional to a length scale, for which the plume radius is natural. This leads to where α and α are positive dimensionless coefficients, and r is a characteristic horizontal spacing between plumes, so that represents the plume radius. For our single updraft, we choose based on Gregory (2001), and , which is broadly consistent with the drag coefficient of 0.225 given by Romps and Charn (2015) when the value given there ( ) is multiplied by a factor of 3/4 arising from assuming spherical geometry of the rising air parcel (Simpson & Wiggert, 1969). We further choose based on single‐column tests and comparison with LES of shallow convection (details in section 4). These parameter choices can likely be improved by making them explicitly dependent, for example, on plume geometry (e.g., Morrison, 2016a, 2016b)

The SGS pressure perturbation term for the environment is assumed to be equal and opposite to the sum of the SGS pressure perturbation terms of all plumes, so that the sum of the SGS pressure terms over the total domain is zero. While this is a conservation requirement, it is not respected by all closures that have been published; some closures introduce additional nonconservative body forces acting on updrafts through the pressure term.

Scalar Source Terms

Various source terms for the conserved scalars θ and q also need to be specified. In our preliminary tests (section 4), precipitation and cloud cover in the environment (as opposed to in plumes) are negligible, while the radiative cooling, horizontal advection, and subsidence tendencies are all horizontally homogeneous. Therefore, no microphysics scheme is needed, and the other tendencies are simply added to the scalar equations (14) and (16).

More general microphysics and radiation schemes can be incorporated consistently in the EDMF scheme, including representations of direct interactions between these processes and the SGS dynamics (instead of only interactions on the grid scale as in most parameterizations). Here, we briefly outline the approach and discuss some of its advantages, while detailed implementations are left for future work.

Microphysical processes (including cloud and precipitation processes) are tightly coupled to turbulent and convective motions. Within the EDMF framework, the corresponding source terms are computed separately in the plumes and in the environment, using the respective subdomain variables. We may further assume a joint probability distribution of the conserved scalar variables in the environment, e.g., a normal distribution with covariances provided by equation (17). This allows for the implementation of probabilistic cloud schemes (e.g., Sommeria & Deardorff, 1977), and, more generally, sophisticated microphysics schemes through a numerical quadrature over the probability distribution. It permits a direct coupling of SGS dynamics to microphysical processes.

Radiative processes are usually assumed to be less strongly coupled to SGS dynamics, and are often computed outside the closures and applied homogeneously across the domain. However, this may be inaccurate, for example, for stratocumulus clouds that only partially cover the domain, because radiative cooling at their cloud tops drives their dynamics (Wood, 2012). Including radiative tendencies directly in the source terms for plumes and the environment, given cloud variables provided by a cloud scheme, may lead to more accurate parameterizations, for example, of stratocumulus. It likewise permits a direct coupling of SGS dynamics to radiation.

TKE Source Terms

In the TKE equation (18), the shear production and SGS transport terms are closed with the EDMF closure, while the buoyancy, pressure, and dissipation terms require additional closures.

The buoyancy flux is calculated separately for the unsaturated (clear, subscript d for dry) and saturated (cloudy, subscript s) parts of the environment, under the assumption that the turbulent fluxes are homogeneous across both parts. We linearize the relation around the environmental mean, i.e., for the clear part, and for the cloudy part. With the ED closure, the buoyancy flux is then parameterized as where is the environmental cloud fraction. The derivatives are given by (Deardorff, 1976) where R and R are gas constants of dry air and water vapor, and c is the isobaric specific heat of air. As a simple approximation, , and are calculated directly from the environmental averages , and . Calculating them separately from the conditional mean values of , and over the cloudy and clear portions of the environment would be more accurate.

In our cumulus test cases (section 4), clouds are concentrated in updrafts. Therefore, we assume that , and thus the buoyancy flux is simply the clear buoyancy flux. For cases with environmental clouds (e.g., stratus, stratocumulus, or anvil‐topped cumulonimbus), the EDMF scheme can be coupled to a probabilistic cloud scheme (e.g., Sommeria & Deardorff, 1977) to obtain nonzero environmental cloud fractions to use in (31).

As in Mellor (1973) and Witek et al. (2011b), we assume that pressure work does not provide a net source/sink of TKE for the total domain. Therefore, we parameterize the pressure term in the environmental TKE equation (18) to be equal and opposite to the sum over all subdomains of the analogous pressure terms involving the subdomain‐mean vertical velocities , obtained by multiplying the vertical velocity equations (13) and (15) for the plumes and environment by . This gives

Separating out the environmental term on the right‐hand side leads to a formulation that relates the pressure term in the environmental TKE equation to the SGS pressure perturbation term (30) for a plume:

The dissipation term is parameterized as: where l is the mixing length (29) and is a constant scaling parameter.

Numerical Implementation

The extended EDMF scheme described above has been implemented in a newly developed single column model. It uses an anelastic reference state as described in Pauluis (2008) and Pressel et al. (2015), and the model equations are discretized on the one‐dimensional analog of an Arakawa‐C grid, with vertical velocities defined on grid cell edges, and all other variables defined at the cell centers.

Each time step of the model begins with an update of the surface flux scheme, as the EDMF scheme requires knowledge of the surface sensible and latent heat fluxes and momentum fluxes. With these quantities known, the lower boundary conditions for the updraft scalars ( ) and updraft area fraction a can be set. All prognostic updraft equations are solved explicitly in time, with adaptive substepping based on a Courant number computed with the maximum updraft velocity. A simple upwind discretization is used for the advective terms, and first‐order operators are used to interpolate values from cell centers to cell edges, or vice versa. The area fraction and vertical velocity equations are integrated upward together, level by level. This approach allows both equations to use consistent values of detrainment rates, incorporating the dynamical detrainment component as described in section 3.2. The updraft vertical velocity equation uses the most current updraft buoyancy value as obtained at the end of the last substep. After completing the update of the area fraction and vertical velocity, the updraft scalar equations are advanced together. For numerical reasons, the current upwind scheme stops integrating whenever updraft velocities become negative ( ), leading to instant detrainment of the entire updraft. This affects the model performance when convection is intermittent (section 4.4). In the future, we will implement an improved numerical scheme that permits downdrafts to eliminate this shortcoming.

Once the values of updraft quantities have been updated, mass‐flux tendencies are applied to the domain‐mean prognostic variables and explicitly. The tendency due to turbulent transport in the environment, parameterized using the eddy diffusivity closure of section 3.3, is treated implicitly.

Finally, equation (18) for the environmental TKE is solved using an upwind, implicit discretization, except for the source terms due to shear, buoyancy fluxes, and entrainment of TKE from updrafts, which are computed explicitly. This completes the update of the EDMF scheme. Any additional tendencies of the domain‐mean variables, such as large‐scale forcings or radiation, are applied explicitly, and the simulation is advanced to the next time step.

This implementation is sufficient for the numerical experiments on which we focus here, which are designed to test the properties of the scheme, not to maximize its computational efficiency. Computationally more efficient implementations (e.g., semi‐Lagrangian or semi‐implicit schemes) are possible and will be explored in the future.

Single‐Column Tests

As a first test, we illustrate the ability of the extended EDMF scheme to realistically simulate the dynamics of boundary layer turbulence and shallow convection by comparison with LES. We compare the extended EDMF scheme (“prognostic EDMF” hereafter) to a more traditional variant (“diagnostic EDMF” hereafter) that is in all respects the same except that it solves steady‐state updraft equations. That is, the diagnostic EDMF scheme eliminates the terms in equations (4), (13), and (14). However, the updraft area fraction still varies with altitude, and the diagnostic EDMF scheme also still solves a prognostic TKE equation, as in Witek et al. (2011b). In the life cycle experiments (section 4.4), we have to further restrict the updraft fraction of the diagnostic EDMF scheme to be no more than 0.7 to ensure its robustness, but otherwise it uses the same parameters as the prognostic EDMF scheme. In addition, numerical stability of the diagnostic form of the EDMF continuity equation requires us to limit the entrainment rate to less than 1/Δz (we choose 0.9/Δz).

Case Description

We focus on a case designed by Siebesma and Cuijpers (1995) and Siebesma et al. (2003), which represents an undisturbed period of the Barbados Oceanographic and Meteorological Experiment (BOMEX) (Holland & Rasmusson, 1973). This BOMEX case features a cumulus‐topped marine boundary layer, in which the interaction between the subcloud turbulence and the cumulus updraft is the essential dynamical process. Cloud radiative effects and microphysics are not as important, because the cumulus cover is relatively small and precipitation does not occur. Therefore, the BOMEX case is an appropriate test for the dynamical aspects of the extended EDMF scheme, without the complications that would arise from processes such as microphysics.

The BOMEX case specification follows Siebesma et al. (2003). The initial profiles are composed of a nearly well‐mixed subcloud layer below , a conditionally unstable cumulus layer up to 1,500 m, and an inversion layer and a free troposphere above. The boundary layer is destabilized by surface sensible and latent heat fluxes and radiative cooling; it is stabilized by large‐scale subsidence. Horizontal advection causes additional drying in the subcloud layer. All of these forcing terms are prescribed and held fixed over the simulation period of 6 h, at the end of which a near‐steady state is established.

We consider two cases: (i) The approximately statistically steady state that develops after approximately 2 h of spin‐up with a uniform and constant surface sensible heat flux of 9.5 W m−2 and latent heat flux of 147 W m−2. (ii) Periodically repeated transient life cycles of convective clouds triggered by surface fluxes that are modulated in time according to where W(t) is a nondimensional weighting function that multiplies the mean surface fluxes of the steady case under (i). This yields fluxes that vary between 1% and 100% of those in the steady case. This life cycle case is designed to test the ability of the prognostic EDMF scheme to capture the time variability of shallow convection, which is not well simulated in many current GCMs (Nuijens et al., 2015a, 2015b). It is a straightforward extension of the widely studied steady BOMEX case, and it may prove to be a generally revealing test of parameterization schemes. In additional tests in which the time‐mean surface fluxes in the life cycle simulations are chosen to match the constant surface fluxes in the statistically steady simulation, we found similar time‐mean states in LES in either case, indicating there is no substantial nonlinear rectification of the oscillations in this setup.

LES and SCM Setup

The LES results are generated with PyCLES, which solves an energetically consistent form of the anelastic equations of motion using total water‐specific humidity and moist‐specific entropy as prognostic thermodynamic variables (Pressel et al., 2015). The model uses a second‐order Runge‐Kutta time stepping scheme (Shu & Osher, 1988) with adaptive time step. Following the implicit LES approach of Pressel et al. (2017), advection of momentum and scalar variables is discretized with a nominally fifth‐order weighted essentially nonoscillatory (WENO) scheme, and no explicit SGS closures are applied outside the near‐surface layer. The domain height is 3,000 m, and the horizontal domain is square and doubly periodic. For the standard BOMEX case with fixed surface fluxes, the horizontal domain size is , with an isotropic mesh size of . A three‐dimensional rendering of this simulation is shown in Figure 1. The transient life cycle case is performed over a much larger horizontal domain, extending on each side, to help ensure the robustness of its time‐varying statistics. The horizontal and vertical grid spacings are reduced to and , respectively.

Figure 1

Visualization of cloud field in standard BOMEX case simulated with PyCLES. Volume rendering of the q field is shown in white, and the lowest‐level q is contoured as indicated by the colorbar.

The single‐column model (SCM) domain likewise is 3,000 m deep, with a vertical mesh size of . With no other adjustments to the model parameters, the results only show modest sensitivity to the choice of grid spacing between 10 and , and this sensitivity may be reduced further by modifying the preliminary closures discussed in section 3 to incorporate dependence on the vertical resolution. The SCM is integrated for 6 hours with a time step of , and adaptive substepping is implemented for the prognostic updraft equations. Large‐scale subsidence tendencies are computed from the prescribed subsidence velocity using an upwind scheme.

BOMEX Results

Figure 2 shows the profiles of domain‐mean liquid water potential temperature and total water vapor specific humidity, and . The prognostic EDMF scheme reproduces the mean profiles simulated by LES reasonably well, albeit with some deficiencies. The prognostic EDMF scheme's subcloud layer is deeper, drier, and warmer than that in the LES; its cloud layer is cooler and moister. The height of the inversion layer is generally consistent with the LES. The diagnostic EDMF scheme produces similar mean profiles as the prognostic EDMF scheme.

Figure 2

Domain‐mean liquid water potential (a) temperature and (b) humidity profiles in BOMEX simulation. (a) Liquid water potential temperature . (b) Total water specific humidity . Both are averages of BOMEX simulations averaged over the sixth hour. Results are obtained with the prognostic EDMF scheme (blue), the diagnostic EDMF scheme (orange), and LES (black). The gray shading shows the plus/minus one standard deviation range of LES models that participated in the intercomparison by Siebesma et al. (2003).

Consistent with the similarities in the domain‐mean profiles of and , the prognostic and diagnostic EDMF schemes produce similar profiles of domain‐mean liquid water specific humidity, (Figure 3a). Both the amount of liquid water and the extent of the cloud layer are similar. Environmental cloudiness is neglected by the EDMF schemes, so their similar predictions of reflect similar updraft area fractions and similar updraft thermodynamic properties (in particular, the extent to which entrainment causes deviations from an adiabatic liquid water profile). The discrepancy between the implicit LES with high‐order WENO schemes and the Siebesma et al. (2003) model intercomparison ensemble in Figure 3a arises because the combined numerical and SGS model dissipation is weaker with implicit LES than with the explicit SGS schemes used in the intercomparison ensemble (Pressel et al., 2017). We have verified that simulated by PyCLES can be brought within the intercomparison range by including a Smagorinsky SGS closure (Pressel et al., 2015).

Figure 3

Domain‐mean liquid water‐specific humidity and TKE profiles in BOMEX simulation. (a) Liquid water specific humidity . (b) Total TKE . Plotting conventions as in Figure 2.

The updraft properties are also similar between the prognostic and diagnostic EDMF schemes. Compared to LES, the updraft of the EDMF schemes has a lower area fraction (Figure 4a), a higher updraft velocity (Figure 4b), a larger updraft mass flux (Figure 4c), and a much higher in‐cloud liquid water specific humidity (Figure 4d). However, the EDMF updraft penetrates to a similar height as that of the LES. These biases can likely be reduced with further optimization of the EDMF parameters. Above , a prominent cloud‐top peak in occurs in the EDMF scheme (Figure 3a), which results from an increase in the updraft area fraction as the updraft rapidly decelerates where it approaches its termination height; the profile in the LES is much smoother and extends higher. This suggests that the closure for the dynamical detrainment in the EDMF scheme needs to be improved. Also, this portion of the cloud layer in the LES is dominated by a few large clouds, which are difficult to capture with a single bulk updraft; multiple updrafts in the EDMF scheme may ameliorate this problem.

Figure 4

Updraft profiles in BOMEX simulation. (a) Updraft area fraction a 1. (b) Updraft velocity . (c) Updraft mass flux m 1. (d) Updraft liquid water‐specific humidity . Plotting conventions as in Figure 2. LES updraft profiles are obtained using cloud conditional sampling as in Siebesma et al. (2003). Cloudy and clear portions of the EDMF‐predicted updraft profiles are indicated by solid and dashed lines, respectively.

To compare the TKE in the EDMF schemes with that in the LES, we consider the domain‐mean TKE , which in the EDMF schemes consists of the environmental TKE and the TKE associated with the mean vertical velocities of the updraft and environment. Both EDMF schemes produce TKE that is too low in the subcloud layer (Figure 3b), because of excessive dissipation. The strong TKE dissipation is needed to limit the diffusivity in the environment. This may be remedied by using a smaller c in equation (28) to compute the eddy diffusivity for TKE. The EDMF schemes do not reproduce well the upper portion of the TKE profile in the LES. However, the TKE in the LES above approximately 2 km altitude is not associated with turbulent fluxes of liquid potential temperature or total water specific humidity (Figure 5), and hence does not appear to lead to irreversible mixing. This suggests that much of the TKE in the LES in this stably stratified region is associated with gravity waves, rather than turbulence, which cannot be expected to be represented by the EDMF schemes, where the TKE only appears as a measure of irreversible mixing.

Figure 5

Domain‐mean vertical fluxes of conserved variables in BOMEX simulation. (a) Flux of liquid water potential temperature θ. (b) Flux of total water‐specific humidity q. Plotting conventions as in Figure 2.

The vertical turbulent fluxes produced by the EDMF schemes agree broadly with the LES results (Figure 5) but are weaker in the upper cumulus layer (above ) where the bulk updraft assumption is less accurate.

Life Cycle Results

Some of the new capabilities of the prognostic EDMF scheme are illustrated by the transient life cycle simulations, in which birth and decay cycles of convective clouds in the LES are triggered on a timescale of an hour (Figure 6a). On such short timescales, the internal memory and time‐dependence of cloud dynamics becomes important and can, to some degree, be captured by the prognostic EDMF scheme. The prognostic EDMF scheme reproduces the birth and decay of clouds (Figure 6b) and the resulting fluctuation in the liquid water path (Figure 6d). It approximately reproduces their phase relation to the forcing variations (Figure 6e), including their growth over time after the surface fluxes peak, and the delayed growth of the first cloud in the simulation. It also roughly reproduces some features of the cloud structure such as their height and their liquid water content. In contrast, the diagnostic EDMF scheme does not capture the cloud life cycle as well: the fluctuations in liquid water path are in phase with the surface forcing (Figures 6c and 6d). However, there are also evident deficiencies in the prognostic EDMF scheme, such as the too early onset of the first two clouds, and the high liquid water specific humidity near the top of the growing clouds, which leads to a time‐mean bias in (right plot of Figure 6b). We have obtained better reproduction of the LES results with stronger dynamical detrainment at the cloud top, which again suggests that refinements of the dynamical detrainment formulation are worth exploring further. The abrupt termination of the clouds (Figure 6b) is a limitation of the plume vertical advection scheme that does not allow for downdrafts (section 3.7), which is likewise in need of further development.

Figure 6

Time evolution of clouds and forcing in life cycle simulations. Domain‐mean liquid water specific humidity in (a) LES, (b) prognostic EDMF scheme, and (c) diagnostic EDMF scheme. Contours are plotted at 0.0005, 0.001, 0.003, and then in intervals of 0.003 up to . The plots on the right show 6 h mean profiles of the liquid water‐specific humidity. The dashed black lines in Figures 6b and 6c show the LES profile of (a) for comparison. (d) Liquid water path (LWP) as a function of time in LES (black), prognostic EDMF scheme (blue), and diagnostic EDMF scheme (orange). (e) Sensible (red) and latent (green) heat fluxes at the surface that are imposed as forcing.

We find results similar to those for the cloud liquid water for the cloud mass flux (Figure 7) and cloud fraction (Figure 8).

Figure 7

Cloud mass flux (i.e., mass flux where ) in life cycle simulations. (a) LES; (b) prognostic EDMF scheme; (c) diagnostic EDMF scheme. Contours are plotted at 0.005, 0.01, and then in intervals of 0.01 up to . As in Figure 6, plots on the right show 6 h mean profiles.

Figure 8

Cloud fraction in life cycle simulations. (a) LES; (b) prognostic EDMF scheme; and (c) diagnostic EDMF scheme. Contours are plotted at 0.005, 0.01, 0.03, and then in intervals of 0.03 up to 0.33. Plots on the right are again 6 h averages.

The success of the prognostic EDMF scheme in capturing the phase relation between the rapid birth and decay of clouds and the oscillatory forcing in the life cycle simulation bodes well for further developments and applications in large‐scale models with relatively high resolution, where capabilities to resolve such life cycles will become crucial.

Discussions and Conclusions

We have presented the theoretical foundations for an extended EDMF scheme that can serve as a unified closure of a broad class of SGS turbulent motions, from boundary layer turbulence to deep convection. The scheme builds on well‐established ideas of modeling convective transport through mass fluxes and more isotropic turbulent transport through diffusive fluxes, unifying them in an EDMF framework as originally proposed by Siebesma and Teixeira (2000) and Soares et al. (2004). However, the new scheme extends the original EDMF approach in several ways. For example, the extended EDMF scheme makes plumes (updrafts or downdrafts) prognostic, to be able to represent convective life cycles with their own memory on the subgrid scale, and it allows for plume cross‐sectional areas to vary with height and time and potentially to cover a large portion of the host model grid box. It also consistently averages the dynamical equations over the subdomains occupied by plumes and the environment, similar to Yano (2014a). This leads to a consistent partitioning of TKE and other second‐order moments between plumes and the environment, which interact through entrainment, detrainment, and pressure drag. These relaxations of common restrictions in more traditional parameterization schemes are crucial for GCMs with horizontal mesh sizes of tens of kilometers to kilometers, which will become routine in the coming years.

While the focus of this paper was on the theoretical foundations, we have also proposed preliminary, physically informed closures for parameters and parametric functions appearing in the extended EDMF scheme, such as an eddy diffusivity that depends on TKE, and entrainment and detrainment rates that depend on updraft vertical velocity and buoyancy. Tests with the preliminary closures show that the extended EDMF scheme is able to simulate an approximately steady state of boundary layer turbulence and shallow convective clouds under constant forcing. It is also able to simulate aspects of the birth and decay of shallow convective clouds in novel transient life cycle experiments we conducted, in which forcing by surface fluxes oscillates. The memory terms in the plume equations are crucial to reproduce the phase relation between the birth and decay of clouds and the oscillatory forcing in the life cycle experiments. At the same time, these tests make clear that further improvements of the closures (e.g., for dynamical detrainment) and likely a more flexible structure with multiple plumes are necessary to improve the goodness‐of‐fit of the extended EDMF scheme to LES statistics.

That the extended EDMF scheme is able to capture many aspects of the steady state and transient dynamics of convective clouds bodes well for further developments. For example, its subgrid‐scale memory may help in the reproduction of the diurnal cycle of convection, which has been a persistent problem for parameterization schemes (e.g., Bechtold et al., 2004; Betts & Jakob, 2002; Dai, 2006; Dai & Trenberth, 2004; Dirmeyer et al., 2012). It may also help with the representation of out‐of‐equilibrium thermodynamic processes, such as the formation of supercooled liquid, which forms because timescales for updrafts and for microphysical processes such as ice fallout and vapor deposition differ (e.g., Korolev & Field, 2008; Morrison et al., 2012; Rauber & Tokay, 1991; Shupe et al., 2008). The time‐dependence of the SGS dynamics in the extended EDMF scheme allows explicit simulation of such multitimescale dynamics. It also makes the traditional distinction between “plume” and “bubble” representations of convection unnecessary (e.g., Yano, 2014b): A time‐dependent updraft plume, which can decay from the bottom up, can represent the dynamics of a rising bubble.

While the theoretical framework of the extended EDMF scheme is designed to be applicable to situations from boundary layer turbulence to deep convection, several issues remain to be addressed to realize successful simulations, for example, of deep convection. These include:

A representation of microphysical processes, including cloud and precipitation processes, that is consistent with the EDMF assumptions needs to be coupled to the scheme.

Cloud‐radiative effects that are important, for example, for the dynamics of stratocumulus clouds need to be incorporated into the EDMF dynamics.

Multiple plumes representing an ensemble of convection types need to be incorporated in the EDMF scheme.

How precipitative downdrafts are to be initialized and how they interact with updrafts needs to be addressed.

Stochastic elements (e.g., in entrainment rates or the number of plumes in the domain) will likely need to be included to improve goodness‐of‐fit of the scheme.

Improvements may be needed in the numerical methods used to integrate the governing equations, particularly the prognostic updraft equations.

We are currently addressing these and other issues, to further develop the extended EDMF scheme to the point that it can serve as a unified parameterization scheme for a broad range of SGS dynamics and GS resolutions in GCMs.