Evaluation of the Bulk Mass Flux Formulation Using Large-Eddy Simulations

Inthis study,bulk massﬂux formulationsforturbulentﬂuxesareevaluatedforshallowand deepconvection using large-eddy simulation data. The bulk mass ﬂux approximation neglects two sources of variability: the interobject variability due to differences between the average properties of different cloud objects, and the intraobject variability due to perturbations within each cloud object. Using a simple cloud–environment decomposition, the interobject and intraobject contributions to the heat ﬂux are comparable in magnitude with thatfrom thebulk mass ﬂux approximation, but do not share a similarverticaldistribution,and so cannot be parameterized with a rescaling method. A downgradient assumption is also not appropriate to parame- terize the neglected ﬂux contributions because a nonnegligible part is associated with nonlocal buoyant structures. A spectral analysis further suggests the presence of ﬁne structures within the clouds. These points motivate investigations in which the vertical transports are decomposed based on the distribution of vertical velocity. As a result, a ‘‘core-cloak’’ conceptual model is proposed to improve the representation of total vertical ﬂuxes, composed of a strong and a weak draft for both the updrafts and downdrafts. It is shown that the core-cloak representation can well capture the magnitude and vertical distribution of heat and moisture ﬂuxes for both shallow and deep convection.


Introduction
The representation of moist convection in general circulation and numerical weather prediction models plays a central role in understanding the multiscale processes of the atmosphere and also the climate sensitivity (Arakawa and Schubert 1974;Randall et al. 2003; Denotes content that is immediately available upon publication as open access. Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JAS-D-19-0224.s1. Arakawa 2004;Bony et al. 2015). The major task of convection parameterization is to represent the subgrid vertical transports due to an ensemble of unresolved convective elements, and specifically their effects on the resolved-scale variables. The majority of current convection parameterizations are based on the bulk mass flux formulation, which approximates the subgrid vertical flux of a scalar quantity as being the product of the convective mass flux with the departure from the gridbox average of the transported quantity (for mathematical details see section 2b). This formulation is based upon a decomposition of the flow field using a top-hat assumption (Randall et al. 1992) or the segmentally constant approximation (Yano et al. 2010). It is also common to assume the model grid spacing to be large enough for grid boxes to contain a large number of clouds and to assume that the area fraction of convection is much less than unity.
As the grid spacing of many global weather and climate models will be reduced to the order of 10 km or even finer, convection can be partly resolved, and this point has motivated reconsiderations and reassessments of these and other convective parameterization assumptions. The bulk mass flux approximation has been evaluated using cloudresolving models (Guichard et al. 1997;Yano et al. 2004) for deep convection and large-eddy simulations (Siebesma and Cuijpers 1995) for shallow convection. These studies found that the bulk mass flux approximation can substantially underestimate the vertical fluxes by 30%-50%, depending on the variable considered and the horizontal resolution. As a result, a parameterization of the neglected contributions to the vertical flux would appear to be necessary. How might this be achieved without sacrificing the computational efficiency, which is arguably the main attraction of the bulk mass flux approach? A drawback of these earlier studies, however, is their relatively coarse resolution by modern standards, so that some of the fine or coherent structures (e.g., cloud-top overturning structures, thin subsiding shells around the cloud, downdrafts within the stratocumulus-topped boundary layer) may not have been well resolved. Such coherent structures have been shown to be important for the vertical transport in recent works (Heus and Jonker 2008;Glenn and Krueger 2014;Park et al. 2016;Davini et al. 2017;Brient et al. 2019). Zhu (2015) investigated the mass flux representation using high-resolution simulations, but did not consider the role of fine structures. It is thus worthwhile to revisit the analysis of the bulk mass flux approximation, and to ask whether the approximation is able to provide an adequate representation of the ensemblemean effect of these fine and coherent structures.
Efforts have been made to parameterize the neglected subplume fluxes. Lappen and Randall (2001), for example, attempted to do so as a downgradient effect in a unified parameterization of boundary layer and moist convection. This basically assumes that these subplume fluxes result from small eddies, which is not necessarily the case since inhomogeneity within the plumes could arise from more substantial internal motions. Moeng (2014) relates the total subgrid turbulent flux to the horizontal gradients of resolvable variables for deep convection. Generally, subplume variability consists of two parts: the interobject variability due to the differences among the average properties of different coherent cloud objects, and the intraobject variability due to the inhomogeneity within each cloud object (see details in section 2b). An assessment of these aspects of variability, including their vertical distributions and any relationships with bulk mass flux terms, is desirable for their parameterization but has not been addressed in previous studies. Here, a thorough analysis of the bulk mass flux formulation and interobject and intraobject variabilities for deep and shallow convection will be performed by using large-eddy simulations. The analysis is designed to investigate several questions: 1) Can the bulk mass flux approximation represent the ensemble effect from fine structures of clouds on the vertical transport? 2) What are the characteristics of interobject and intraobject variability that constitute the subplume fluxes? 3) What are the key elements that need to be considered in convection parameterization in order to provide an efficient and accurate representation of the vertical fluxes of both heat and moisture using a mass flux approach?
The paper is organized as follows. Section 2 describes the large-eddy simulations and introduces the bulk mass flux approximation alongside formulations for the neglected inter-and intraobject variability. The algorithms to identify coherent cloud objects are presented in section 3. Section 4 provides an evaluation of bulk mass flux approximation, the features of interobject and intraobject variability and spectral representation of them, and points out the necessity of understanding the fine structures of clouds. Section 5 investigates the key elements that are responsible for vertical transport and a core-cloak conceptual model is proposed to improve the mass flux approximation. A discussion and a summary are provided in sections 6 and 7, respectively.  Brown et al. 2015Brown et al. , 2018 model is used for the large-eddy simulations of both shallow and deep convection. The simulation of shallow convection is based on the Barbados Oceanographic and Meteorological Experiment (BOMEX), and the model configuration follows that of Siebesma et al. (2003). The grid spacing is 25 m in all directions and the domain size is 15 km 3 15 km. The 3D Smagorinsky-Lilly scheme (Smagorinsky 1963;Lilly 1962) is used for the parameterization of subgrid turbulence. A simple saturation adjustment cloud scheme is used to represent the conversion between water vapor and cloud liquid water as this is a nonprecipitating case without ice water.
The evaluation of deep convection is based on a radiative-convective equilibrium (RCE) simulation. The simulation has a horizontal resolution of 200 m and domain size of 132 km 3 132 km. The model top is at 40 km, using 99 stretched vertical levels. Sea surface temperature is held fixed at 300 K, and surface pressure is 1000 hPa. The simulation is initialized with horizontally homogeneous tropical profiles of potential temperature and water vapor. Constant initial horizontal wind profiles are specified (U 5 5 m s 21 , V 5 0 m s 21 ), and the domain-mean wind fields are relaxed toward these values with a time scale of 6 h. A prescribed cooling profile is applied in order to destabilize the atmosphere, and is 1.5 K day 21 from the surface to 12 km, with a linear decay to zero at 16 km. The microphysics is parameterized using Cloud Aerosol Interaction Microphysics (CASIM; Grosvenor et al. 2017;Miltenberger et al. 2018) in double-moment configuration. The subgrid turbulence is parameterized through the 3D Smagorinsky-Lilly scheme (Smagorinsky 1963;Lilly 1962).

SIMULATIONS
To further evaluate the ''core-cloak'' representation proposed in section 5c, we also consider RCE simulations using interactive radiation, with and without selfaggregation. An aggregated simulation is performed over a 100 km 3 100 km domain at 1-km horizontal resolution with 300-K sea surface temperature. This simulation is a part of the Met Office Unified Model (UM; Davies et al. 2005) contribution to the Radiative-Convective Equilibrium Model Intercomparison Project (RCEMIP; Wing et al. 2018), which is designed to investigate cloud and climate sensitivity, quantify the dependence of the degree of convective aggregation on temperature, and to assess robustness across a spectrum of models. Details of the simulation design are available in Wing et al. (2018). A nonaggregated simulation uses the same configuration except that it homogenizes the radiative tendencies at each time step. As the interaction between radiation and water vapor or cloud plays a key role in self-aggregation (Bretherton et al. 2005;Muller and Held 2012;Wing and Emanuel 2014;Muller and Bony 2015), the organization of convection is inhibited in this second simulation. Comparison between these two simulations is conducted to assess the robustness of a core-cloak representation in organized convection.
All of the calculations in this study are taken from periods when the simulations have achieved an equilibrium state. For the BOMEX simulation, we take data at 10-min intervals from hour 5 to hour 6. For the RCE and RCEMIP simulations, our evaluation data is sampled every 6 h for the last 5 days of simulation. The RCE and RCEMIP simulations last for 54 and 125 days, respectively.

b. Decomposition of total resolved vertical turbulent transport
In this study, we will only consider the resolved vertical fluxes of scalars. The subgrid turbulent fluxes have been checked and are small compared to the resolved fluxes in these large-eddy simulations (not shown). We have also applied the analysis to simulations with different resolutions (1 km, 400m, and 200 m for RCE simulation; 100, 50, and 25 m for BOMEX simulation) and the conclusions do not change.
At each vertical level and time, multiple convective objects can be identified using certain criteria and these are scattered across the domain (see details in section 3). The remaining part of the domain is considered as the environment. Each object is composed of a coherent cluster of contiguous grid points that are identified as updraft or downdraft and is denoted with a subscript i. For convenience of presentation, the environment is also considered as an extra object denoted by i 5 0. An atmospheric quantity within the object is f i , the average of this quantity over the object is f i , and the perturbation from the average over the object is The domain average is denoted as hfi, and the departure from the domain average is denoted as f i * 5 f i 2 hfi. The difference between the average of an object and the average of the full domain is denoted as f i * , and follows from the definition of f i * ; that is, f i * 5 f i 2 hfi. The area fraction of each object is denoted by a i . By definition, the domain average can be computed by hfi 5 å n i50 a i f i , where n is the number of identified objects. If we also apply the same definition to the vertical velocity, the total vertical turbulent flux of f can be represented as JUNE 2020 G U E T A L .
hw*f * i 5 hwfi 2 hwihfi . (1) The last step uses the identity å n i50 a i f i * 5 0 according to the definition of f i * . The domain-average vertical flux of quantity f can be divided into two terms. Term (1.1) is due to the difference between the average of each object and the domain average and here is called the mass flux term. The reader should keep in mind that this term is different from the ''mass flux'' in conventional convection parameterizations, which would include a factor of density and refers to the vertical transport of air mass. Term (1.2) is due to the perturbations within each object [term (1.2a)] and within the environment [term (1.2b)].
Instead of considering each object explicitly, simplifications could be made by parameterizing the vertical fluxes for selected objects under certain conditions. For example, in a conventional convection parameterization, the bulk plume is an ensemble of all the updrafts. Thus, it is equivalent to a collection of grid points in the LES within a particular category (updrafts) and these grid points do not necessarily need to be physically connected. To simplify the representation in this manner, we define the average of f over all updraft objects as f p 5 The superscript p means that all of the identified updraft objects have been collected together as a single draft and the areaweighted average is taken over all such objects. The domain averages hfi and hwi can now be expressed as hfi 5 a 0 f 0 1 (1 2 a 0 )f p and hwi 5 a 0 w 0 1 (1 2 a 0 )w p , with the downdrafts considered here to be part of the environment. Term (1.1) on the right-hand side (rhs) of Eq. (1) can then be decomposed as follows: (2) Term (2.2) in Eq.
(2) can be further simplified as Substituting for terms (2.1), (2.2), and (2.3) into Eq. (1), the total vertical flux can be written as On the rhs of Eq. (4), term (4.1) represents the vertical flux due to the difference between the bulk average and domain mean. It has contributions from the environment and from a bulk plume composed of all updraft objects. Using the definitions of hfi and hwi, term (4.1) can be manipulated as Term (4:1) 5 a 0 (1 2 a 0 )(w p 2 w 0 )(f p 2 f 0 ) . (5) In conventional convection schemes, the area fraction of updrafts (1 2 a 0 ) is assumed to be much less than 1 within a GCM grid box so that term (4.1) can be expressed as (1 2 a 0 )(w p 2 hw 0 i)(f p 2 hf 0 i). Term (4.2) represents the contribution due to the difference between the average of each updraft object and the average of all the updraft objects and is called the interobject variability in this study. Term (4.2) would vanish if we were to assume that all of the objects composing the bulk updraft were the same. Term (4.3) results from the fluctuations within each object and is called the intraobject variability. It would vanish if we adopt the top-hat assumption. Approximating the vertical flux using term (4.1) only is called the bulk mass flux approximation and has been widely used in convection parameterization. Equation (4) only accounts for updrafts and an environment. However, contributions from downdrafts can also be important, and a simple generalization of the derivation from Eqs. (2)-(4) leads to where c u 5 and c may represent vertical velocity w or a transported quantity f. The superscripts u and d indicate an area-weighted average over all updraft and downdraft objects, respectively, while a u and a d represent the total area fraction of the updrafts and downdrafts. This decomposition will be assessed in section 5b. Terms (6.1)-(6.3) in turn generalize terms (4.1)-(4.3). Note that the interobject variability, term (6.2), is no longer due to the difference between the average of each object and the bulk plume average, but results from the difference between the average of each object and the bulk updraft or downdraft.
The mass flux term, term (6.1), and the intraobject variability, term (6.3), are now divided into three contributions from updrafts, downdrafts and the environment.

Definition of objects and drafts
To evaluate the bulk mass flux approximation, we first need to define the objects under consideration. As described in section 2b, an object is a collection of spatially adjacent grid points each of which satisfies certain criteria. There are various ways to define the cloud objects such as using cloud water, perturbation of virtual potential temperature and vertical velocity, individually or a combination of these. We first apply a traditional sampling method, that is, small thresholds of cloud liquid water, q l . 10 25 kg kg 21 in BOMEX or liquid water and ice q l 1 q i . 10 25 kg kg 21 in RCE, to label grid points as cloudy. Contiguous labeled grid points are identified as an individual object by checking the neighboring grid points (south, north, west, and east) around the cloudy points until no more cloudy points are found. We will use this algorithm to investigate the interobject and intraobject variability in section 4. We also combine thresholds of cloud water and positive buoyancy in section 4 to examine the subplume fluxes contributed by the cloud core.
We have also investigated the application of criteria based on labeling grid points using percentile thresholds of vertical velocity. Different types of updrafts and downdrafts can be further defined based on different percentile ranges. For example, in the BOMEX simulation we investigated a three-draft partition (weak, medium, strong). At each vertical level, we produced distributions of vertical velocities for upward and downward motions. Grid points exceeding the top 0.1% of upward vertical velocity were identified as strong updrafts, those within the top 0.5%-0.1% of upward motions were identified as medium updrafts, and those within the top 5%-0.5% of the upward motions were identified as weak updraft. The same percentile ranges were also applied to downward motions to identify the weak, medium and strong downdrafts, and the rest of the domain is considered to be the environment (Fig. 1a). Each type of draft is therefore an ensemble of grid points within velocity space. This algorithm will be used to evaluate a multidraft model in section 5c, and in particular a core-cloak representation, in which the core represents the strong drafts and the cloak represents weak drafts. Compared to a multiobject algorithm, this algorithm continues to identify fine structures in the cloud objects but it merges similar parts of the objects together as abstract drafts rather than dealing with individual objects explicitly.
The use of a fixed percentile of vertical velocity is somewhat different from previous studies where a fixed value has been used to identify convective clouds, sometimes with an additional cloud liquid water threshold. Recent studies Efstathiou et al. 2020) identified the coherent structures by optimizing the vertical transport of scalars (e.g., total water and liquid water potential temperature). Such methods can characterize structures contributing the most to the vertical transport and covering the smallest area fraction possible. But the identified structure may be different depending on what flux the algorithm aims to optimize. This is because the distributions of different scalars (e.g., cloud water and potential temperature) differ from each other (see section 4b). The percentile method is taken to be preferable here, in part because we wish to treat the shallow and deep cases on the same basis, and it would be difficult to choose a suitable value threshold for different types of convection at different heights. The use of percentile thresholds, calculated separately at each time and each level, to detect the objects and drafts ensures that only the grid points on the tail of the distribution are chosen. Another advantage is that dry drafts are detected in the subcloud layer, so that we can extend the assessment of the bulk mass flux formulation there also, as has been adopted in eddy diffusivity-mass flux (EDMF) parameterizations (e.g., Siebesma and Teixeira 2000).
As an example illustrating the draft and object definitions, a snapshot from the BOMEX case is shown with horizontal and vertical cross sections in Figs. 1b and 1c, respectively. A snapshot from the RCE simulation shows similar features (see Fig. S1 in the online supplemental material). Most of the strong updrafts (top 0.1%) are collocated with cloud and form the core of individual cloud objects (e.g., cloud A in Fig. 1c). They are surrounded by medium and weak updrafts. Some clouds have downdrafts on their periphery, indicating a shell structure (Heus and Jonker 2008). Other clouds do not have detected updrafts but do have strong downdrafts in their vicinity (e.g., cloud B in Fig. 1c). Such clouds are in the decaying stage of their life cycles, when the upward vertical velocities within the clouds are no longer on the tail of the distribution. There are also some updrafts that can be seen, but without any cloud liquid water. Some of these updrafts are associated with gravity waves propagating away from the convection (e.g., the updraft signals above the cloud A). Others are in their developing stage and clouds have not formed yet. In addition, our decomposition also identifies clouds (specifically clouds C and D in Fig. 1c) that have just begun to form and so have low cloud tops and are still connected with their dry precursors in the subcloud layer. This means that it would be possible to study the life cycles of convection throughout the vertical range extending from subcloud layer to cloud top if the decomposition were to be combined with suitable 3D object tracking. We do not pursue that here, but simply observe that our draft decomposition can capture the gross features of clouds from cloud base to cloud top, even though no conditions on cloud liquid water has been applied. Figure 2 shows vertical profiles of vertical velocities corresponding to different percentiles and also the averaged vertical velocity of cloud and buoyant cloud (defined as q l 1 q i . 10 25 kg kg 21 and u 0 y . 0, where u 0 y 5 u y 2 hu y i). In the BOMEX case, the cloud-mean vertical velocity is close to the top 5% threshold near cloud base and cloud top and close to the top 1% threshold in the rest of the cloud layer. The mean vertical velocity of the buoyant cloud core increases with height to exceed the top 0.5% threshold (Fig. 2a) above 1 km. In the RCE case, the distribution of vertical velocity is more skewed toward extreme positive values (Fig. 2b). The mean vertical velocity of cloud is close to the top 0.5% threshold between 3 and 6 km and close to top 5% threshold below 1 km and above 8 km. The vertical velocity of buoyant cloud is 1-2 m s 21 larger than the cloud-mean value. Both in-cloud profiles have a maximum near 6 km. In this study, to keep consistency for both shallow and deep convection, unless otherwise noted, the percentile thresholds of top 0.5% and top 5%-0.5% bins are taken as indicative of the updraft cores and weak updrafts, respectively.

Cloud-environment decomposition a. Inter-and intraobject variability
Section 2 showed that subplume turbulent fluxes consist of two contributions, due to interobject and intraobject variability. Understanding their features is necessary to examine whether a downgradient assumption or a rescaling method may be reasonable to parameterize them.  The bulk mass flux approximation works well for the vertical transport of total water mixing ratio q t , and liquid water potential temperature u l (Figs. 3c and 3d; note that u l 5 u[1 2 (L y q l /c p T)]; Betts 1973) since the cloud objects are defined using cloud liquid water. It captures about 80% of the total fluxes and its vertical distribution is similar to the total fluxes. The interobject and intraobject variability within the clouds are small and share similar shapes. The environmental variability is comparable in magnitude but has opposite sign to these through most of the cloud layer. While the vertical fluxes of q t and u l may be enough for nonprecipitating shallow convection, the vertical heat fluxes also need to be evaluated since they are typically used in most numerical models and are important for the parameterizations that predict the turbulent kinetic energy using the buoyancy flux as an important source term.
However, the bulk mass flux approximation provides a rather poor representation for the vertical fluxes of u and u y (Figs. 3a,b). It is negative throughout the cloud layer for hw*u * i while the total flux is positive from 800 to 1800 m (Fig. 3a). For the vertical buoyancy flux hw*u y * i, the bulk mass flux approximation has the opposite sign to total flux in the inversion layer (Fig. 3b). The inter-and intraobject variability within the cloud are comparable with bulk mass flux approximation and with the environmental variability. Most importantly, these terms do not share similar vertical profiles with each other, nor with the bulk mass flux approximation (Figs. 3a,b). As a result, the subplume fluxes cannot be parameterized by rescaling the bulk mass flux contribution. This is because the vertical gradients of interand intraobject variability have opposite signs at some levels (e.g., from 1000 to 1500 m). On the other hand, a large part of subplume fluxes is associated with the buoyant cloud (gray lines in Fig. 3) instead of small-scale eddies. In addition, the total subplume fluxes do not share similar shapes with vertical gradient of mean u y profile (not shown) and therefore the downgradient assumption is also not sufficient to reproduce all the subplume fluxes.
For the RCE simulation, the bulk mass flux approximation based on the traditional cloudy sampling cannot well capture the vertical fluxes of both heat and moisture (not shown), especially at high levels. This is because the anvil clouds in the upper troposphere cover a large area but have small vertical velocities. The top-hat assumption gives a small mean vertical velocity over the cloudy regions and thus results in significant underestimation.

b. Spectral distribution of vertical fluxes
The different performance of the bulk mass flux approximation for the vertical heat and water fluxes under a cloud-environment decomposition indicates that the internal distributions of temperature and cloud water within the cloud are different. This implies that a bulk cloud is unable to represent well both the temperature and cloud water variability. One way to reduce the subplume fluxes would be to deal with each cloud object explicitly in Eq. (1), which would eliminate interobject variability. Although treating each object explicitly is impractical, we might hope that a spectral parameterization of convection would be able to reduce interobject variability substantially, under the assumption that objects with similar sizes differ much less than those with different sizes.
To explore this idea, the resolved turbulent flux, mass flux term, and intraobject variability are calculated separately for each cloud object, and the statistics are collected with respect to cloud size for BOMEX in Fig. 4. The size of each cloud object is defined as the equivalent size (square root of area coverage). For the turbulent heat flux, intraobject variability (Fig. 4c) is positive and dominates the total flux (Fig. 4a)  near cloud top, where the total heat flux is negative. The mass flux term is negative for small-sized clouds (,200 m) but is weakly positive for medium and large clouds below the inversion layer (Fig. 4b). From 1500 m and above, the mass flux term is negative across almost the whole cloud spectrum and makes an important contribution to the total heat flux (Fig. 4b).
Turning to the buoyancy flux, we find that the mass flux term (Fig. 4e) is the major contributor to the total flux (Fig. 4d). It has a maximum (or minimum) for mediumsize clouds of 200-300 m throughout the cloud layer. While the intraobject variability (Fig. 4f) is relatively small at cloud top, it is comparable with the mass flux term at 1600 m and about half of the mass flux term below the inversion layer. Typically, the intraobject variability for medium size clouds is about 1/3 of the total turbulent flux and thus is nonnegligible. Our results therefore indicate that a spectral method is not enough to provide a good representation for turbulent fluxes by just using the mass flux approximation. This is because there are finerscale structures responsible for vertical transport within each cloud object. This will be the focus of next section.

Key elements for vertical transport
While the cloud condensate is the most intuitive criteria for cloud object identification, it may not be the best choice for an efficient description of vertical fluxes produced by finer structures: for example, overturning circulations near cloud top. In this section, we examine the key elements for describing the vertical fluxes step by step by using the decomposition based on the vertical velocity distribution as described in section 3.

a. Bulk updraft and environment
We begin with the simplest possibility and decompose the domain into two parts: updrafts and environment. For both the BOMEX and RCE simulations, the updrafts are identified as the vertical velocity exceeding the top 0.5% of upward motions. This threshold was found to be most suitable for maximizing the contribution of the bulk mass flux term to the turbulent fluxes of heat, which is significantly underestimated by traditional cloud sampling and it also approximately captures the core of the updrafts (section 3). Once the updrafts are identified, the ensemble of them is considered as a bulk updraft. Figure 5 shows the total resolved turbulent fluxes in BOMEX and the contributions from the bulk mass flux approximation [term (4.1)], interobject [term (4.2)], and intraobject [term (4.3)] variability. In the cloud layer, the bulk mass flux approximation can capture the gross feature of the total fluxes. The interobject and intraobject variabilities within the updrafts are very small, presumably in part because of the small area fraction we set for the decomposition. The variability in the environment dominates the total fluxes in the subcloud layer. In other words, the largest vertical motions do not play a major role in the subcloud fluxes. The environmental variability has two peaks above cloud base: one in the lower part of the cloud layer and one in the inversion layer, and at those heights it has a similar importance to the bulk mass flux term.
For the vertical fluxes of heat hw*u * i (Fig. 5a) and buoyancy hw*u y * i (Fig. 5b), the bulk mass flux term accounts for most of the total fluxes from cloud base above to just below the inversion layer. However, within the inversion layer, the bulk mass flux term makes a strong negative contribution while the total flux is positive or near zero. This indicates the presence of overshooting updrafts with negative buoyancy. The positive contribution from environmental variability in the inversion layer might arise from negatively buoyant downdrafts associated with the overshooting updrafts. We return to this point in section 5b. For the fluxes of total water hw*q t *i (Fig. 5c) and liquid water potential temperature hw*u l * i (Fig. 5d), the bulk mass flux term captures 50% or less of the total fluxes in the lower part of the cloud layer, where the environmental variability plays an important role. This is worse than the bulk mass flux approximation using traditional cloudy sampling as some cloudy points have been considered as environment by the decomposition based on vertical velocity.
In the RCE simulation, the bulk mass flux term is a major component of the total fluxes in the free troposphere (Fig. 6). The interobject variability is comparable with intraobject variability within the updrafts, both being small throughout the troposphere. These terms do not have similar shapes compared with bulk mass flux term. The environmental variability has a similar magnitude to the interobject variability for the heat (Fig. 6a) and buoyancy (Fig. 6b) fluxes. For the fluxes of q t (Fig. 6c) and liquid-ice potential temperature u li ( Fig. 6d; note that u li 5 u[1 2 (L y q l /c p T) 2 (L s q i /c p T)]; Tripoli and Cotton 1981), the environmental variability is comparable with bulk mass flux term and thus is nonnegligible. It has two maxima: one in the lower troposphere, and another in the upper troposphere (near 11 km). The anvil structures emerging from deep convection could explain the maximum of environmental variability in the upper troposphere. In the anvil cloud at high levels, the vertical velocities are small and thus are defined as environment in our decomposition even though they are responsible for part of the vertical transport. Similar to the BOMEX simulation, the environmental variability dominates in the subcloud layer. This indicates that the vertical fluxes in the lowest part of the atmosphere are mainly contributed by drafts with less extreme vertical velocities.

b. Updrafts, downdrafts, and environment
To see the role of downdrafts in the turbulent fluxes, we use a 0.5% threshold to pick up the strong downdrafts, consistent with the threshold for updrafts. The various contributions to the turbulent fluxes are now calculated according to Eq. (6). There is no significant improvement of the bulk mass flux approximation in the RCE simulation (not shown). However, in the BOMEX simulation, the mass flux approximation [term (6.1)] is improved near cloud top for all fluxes considered (cf. Figs. 5 and 7). This is due to the reduction of environmental variability because extreme downward motions

J O U R N A L O F T H E A T M O S P H E R I C S C I E N C E S VOLUME 77
near cloud top have been identified as separate downdraft contributions instead of as the environment. The improvement emphasizes the importance of overturning structures near cloud top. These structures entrain dry air, initiate downdrafts to penetrate around the cloud edge through evaporative cooling and form a shell structure (Blyth et al. 1988;Heus and Jonker 2008). Therefore, the calculation suggests a model that includes a representation of downdrafts near cloud top would be beneficial to better represent the turbulent fluxes. Despite the improvement of the bulk mass flux approximation near cloud top, the intraobject variability in the environment still explains a nonnegligible portion of total fluxes in the lower part of the cloud layer. However, this term can be made negligible if the percentile threshold for updrafts is relaxed to cover the top 5% (see in section 5c) and in the case the mass flux term accounts for most of the total fluxes. This point indicates the potentially important role of less extreme updrafts in maintaining vertical flux in the lower part of cloud layer. The dominance of environmental variability within the subcloud layer for both shallow and deep convection suggests that the vertical fluxes in the boundary layer are controlled by relatively weak vertical motions (predominantly the top 30%-40%, see in section 6). In summary, the above analysis indicates that the turbulent fluxes are composed from drafts with a range of magnitudes, and that representing the total fluxes with a bulk model (with traditional sampling or vertical velocity sampling) results in underestimation. The plume model used in a convection parameterization needs to include at least strong and weak drafts.

c. Improving the mass flux approximation-Core-cloak representation
Basing a parameterization on the bulk mass flux approximations of terms (4.1) or (6.1) above, neglects the contribution from interobject and intraobject variability. As we have found, interobject variability may be important if a system has a broad spectrum of cloud size, and intraobject variability may be important if cloud objects have complicated spatial distributions of different quantities due to the complexity of internal updraft dynamics and their interaction with the environment. Furthermore, the intraobject or interobject variability may become more important considerations as the grid spacing of GCMs decreases to O(10) km or less because of the much more limited sampling of cloud objects (Plant and Craig 2008;Sakradzija et al. 2016). The complexity of parameterizing interobject and intraobject variability using terms (4.2) and (4.3) or terms (6.2) and (6.3) is that physically coherent objects need to be considered explicitly.
In section 4b we considered simplifying the problem using a spectrum of cloud sizes. Here, we consider a possible simplification by collecting together similar parts of the flow as different types of drafts. For example, we might categorize the updrafts or downdrafts into three types: strong, medium and weak. As discussed in section 2b, each type of draft would be composed of multiple disconnected objects that have a similar range of vertical velocity. While the definitions of interobject and intraobject variability in section 2b use the concept of physically coherent objects, the mathematical derivation does not need this constraint and is easily extended to abstract drafts. In this case, the total vertical flux can be written as Note that the subscript j does not label the coherent objects as did i in Eqs. (4) and (6) but rather it labels the different types of drafts. Term (7.1) is the mass flux term and is analogous to terms (1.1) and (6.1). The interobject variability from section 2b is now absorbed within term (7.1). As a result, only the intraobject variability is retained in term (7.2) and may need parameterization. This is referred to as intradraft variability hereafter, because there are no longer explicit objects but abstract drafts. As shown in section 4, a major reason that the bulk plume model fails to approximate the total turbulent fluxes with the bulk mass flux approximation is that the bulk model only describes the mean property of the ensemble of drafts while the vertical transport is actually controlled by the combination of drafts with different values of vertical velocity. The idea is that Eq. (7) could form the basis of a computational cheap multidraft model that includes the major components responsible for the full fluxes vertical transport. More specifically, the hope is that the intradraft variability may be smaller than the intraobject variability in the bulk plume model because part of the intraobject variability has been captured by the different draft types. Equation (7) is written for three types of draft but a simpler starting point for evaluating the idea is to take a two-draft representation, composed of strong and weak updrafts and downdrafts plus the environment. Figure 8 demonstrates the decomposition of the buoyancy and total water fluxes for a weak-strong draft representation in BOMEX. The percentile bins of top 5%-0.5% and top 0.5% are used to pick up the weak and strong drafts, since the top 5% value is close to the cloud-mean vertical motion, and the top 0.5% should capture the core of the clouds. In comparison with the bulk model from Fig. 5, which is shown in Figs. 8a and 8d as the blue dashdotted lines, and the bulk mass flux approximation based on traditional cloud sampling (the gray dashdotted lines in Figs. 8a,d), the mass flux approximation of term (7.1) has been improved to better match the total FIG. 7. As in Fig. 5, except that the contributions associated with downdrafts are now included. We use the top 0.5% threshold to identify updrafts and downdrafts. The red line represents the total flux, the blue line the mass flux approximation [term (6.1)], the magenta line the interobject variability [term (6.2)], the green solid line the intraobject variability within the updrafts and downdrafts [term (6.3a)], and the green dashed line the intraobject variability within the environment [term (6.3b)].

J O U R N A L O F T H E A T M O S P H E R I C S C I E N C E S VOLUME 77
buoyancy flux (Fig. 8a). The improvement mostly comes from reduction of intraobject (or intradraft) variability in the lower part of the cloud and at the cloud top (cf. Figs. 8a and 6b). The mass flux term is controlled by the strong updrafts throughout much of the cloud layer (magenta line in Fig. 8b). The weak updraft plays an important role in successfully capturing the flux from cloud base to 1000-m height (yellow line in Fig. 8b). As a result, the mass flux approximation is improved by more than 40% in the lower part of the cloud layer. The strong downdraft controls a large portion of vertical buoyancy flux near cloud top (1600-2000 m, Fig. 8b). The vertical structures of contributions to intradraft variability within the cloud layer from the different drafts are also consistent with the mass flux term, with strong updraft dominating most parts of the cloud, weak updrafts contributing to the lower part of the cloud and strong downdrafts controlling the values around cloud top ( Fig. 8c). In the subcloud layer, the environmental variability accounts for most of the total fluxes. The decomposition of the turbulent flux of total water for a weak-strong draft model (Figs. 8d-f) has broadly similar characteristics to that of the buoyancy flux. The mass flux term is significantly improved over that in the bulk plume model by up to 50% below 1000 m (Fig. 8d) and the improvement mainly comes from the contribution of weak updrafts (Fig. 8e). One difference is that the weak downdrafts also contribute negatively to the total water flux and have similarly sized contributions to the strong downdrafts throughout the cloud layer except near to the cloud top, where the strong downdrafts dominate (Fig. 8e). This illustrates the role of shell structures surrounding the cloud in transporting moist air downward. Weak updrafts are important below 1000 m and improve the mass flux approximation in the lower part of the cloud layer. Another point of , and (f), the contributions are shown for weak updrafts (wu, yellow), strong updrafts (su, magenta), weak downdrafts (wd, black), strong downdrafts (sd, blue), and the environment (env, cyan). difference is that the weak updrafts make a nonnegligible contribution within the cloud layer (Fig. 8e) whereas they contributed to the buoyancy flux only below 1000 m (Fig. 8b). This point serves to exclude the possibility that the weak updrafts identified in our decomposition are mostly signals of gravity waves outside the clouds, because in that case the corresponding contribution to the flux of total water would be very small throughout the cloud layer. The vanishing buoyancy flux by weak updrafts above 1000 m suggests rather that the weak updrafts cover a transition zone where the buoyancy changes from positive to negative due to the turbulent mixing between the updraft core and the environment, and hence overall vertical transport of buoyancy is near zero.
The above analyses are consistent with a picture of drafts that originate from the subcloud layer and ultimately make their way to the inversion layer. Within the subcloud layer, further experimentation with percentile thresholds (see in section 6) reveals that the top 30%-40% of updrafts transport moisture and heat upward. Only updrafts in the top 5% with positive buoyancy then survive to make important contributions to fluxes within the cloud layer. Ultimately, only the more extreme drafts within the top 0.5% are able to penetrate throughout the full depth of the cloud layer and end within the inversion layer. Cloud-top overturning initiates strong downdrafts that also make a nonnegligible contribution to the total fluxes near the inversion layer.
The same decomposition is applied to deep convection in the RCE simulation (Fig. 9). The bulk mass flux approximation based on traditional cloud sampling significantly underestimates the vertical fluxes (gray dashdotted lines in Figs. 9a,d). Compared to the bulk plume representation (blue lines in Figs. 6c,d and also blue dash-dotted lines in Figs. 9a,d), the two-draft model improves the mass flux approximation for the fluxes of total cloud water (Figs. 9a-c) and liquid ice potential temperature (Figs. 9d-f). However, there is little improvement for the heat and buoyancy fluxes (not shown), perhaps because the deep convective core is more collocated with positive buoyancy. This further indicates the different spatial distributions of variables within the drafts. For the vertical flux of cloud water, the mass flux approximation is improved by about 30% between 6 and 12 km and by up to 50% between 1 and 2 km (Fig. 9a). FIG. 9. As in Fig. 8, but for the results from the RCE simulation for the vertical profiles of the time-(last 5 days) and domain-averaged vertical fluxes of (a)-(c) total water hw*q t *i and (d)-(f) liquid ice potential temperature hw*u li * i. The bulk mass flux approximation from term (4.1) based on updraft sampling (blue dash-dotted line) and the bulk mass flux approximation based on cloud sampling (gray dashdotted line) is also plotted for comparison.
The intradraft variability is reduced to about 8% of the resolved flux below 6 km and about 16% above 8 km. The strong updraft dominates the mass flux term (Fig. 9b) with the weak updraft important in describing the two peaks of intradraft variability that occur at low levels and above 8 km in the bulk plume representation (cf. Figs. 9c and 6c). The contributions of weak and strong downdrafts to the cloud water transport are comparable but are relatively small throughout the whole troposphere with maxima near the cloud top (Fig. 9b). For the vertical flux of liquid ice potential temperature, the main improvement to the mass flux approximation is by about 50% at upper levels (above 8 km, Fig. 9d), where strong and weak updrafts contribute comparably and strong and weak downdrafts also make nonnegligible positive contributions (Fig. 9e).
An important aspect of the improvements obtained from a two-draft representation compared to a bulk model is that better shapes are produced for the vertical profiles of fluxes (e.g., the peaks occur at similar heights as the total resolved fluxes). This is true for both deep and shallow convection and is important because the tendency of a variable within a convection parameterization is determined by the vertical gradient of the vertical fluxes and is essential for vertical distributions of heat, moisture, and hydrometeor (Wong et al. 2015).
We have also extended the two-draft model to a threedraft model, as in Eq. (7), with weak, medium, and strong drafts for both updrafts and downdrafts. For example, one way to do so would be to further split the strong draft in the two-draft model into separate medium and strong drafts in order to account for more intradraft variability in the cloud layer. However, the improvement was found to be minor for both shallow and deep convection (not shown). Therefore, a twodraft model seems to be an attractive approach for the free-tropospheric fluxes considering that the intradraft variability is much reduced (Figs. 8c,f and 9c,f) although no doubt further efforts could be made to refine the definitions to further improve its formulation.
Our results suggest a possible extension of the bulk plume model that is applied in many current convection parameterizations in GCMs. We call this two-draft conceptual model a core-cloak representation of convection and a schematic is shown in Fig. 10. The collection of strong updrafts is depicted as the core while the cloak corresponds to the collection of weak updrafts. This core-cloak structure is also applied to the downdrafts. Parameterization of this core-cloak model would need careful treatment of exchanges between the different types of drafts. As per the schematic, we would anticipate a treatment in which the strong drafts are only able to entrain (detrain) air from (to) the weak drafts, while the weak drafts would behave as a buffer region that can entrain (detrain) air from (to) both the environment and the strong drafts. The updraft and downdraft can be coupled through cloud-top overturning structures. Given that the intradraft variability of the strong and weak downdrafts is very small in both shallow and deep convection (Figs. 8c,f and 9c,f), the core-cloak representation could credibly also be simplified by allowing weak and strong updrafts but only one type of downdraft.
One may question that our sampling method based on the vertical velocity would pick up signals associated with gravity waves or isolated motions as the ''cloak'' part in our conceptual model. These weak drafts may contribute to the total mass flux but do not contribute to the vertical fluxes of scalars. To investigate this further, we have performed an additional analysis to identify the objects that have both core and cloak structures and are also cloudy. At each vertical level, we first identify the objects using the top 5% percentile threshold for upward and downward motions. However, only the objects that include the grid points with vertical velocity exceeding top 0.5% threshold and also have cloudy points are retained to calculate vertical fluxes, and are named to have core-cloak structure. With this sampling, these objects are most probably not associated with gravity waves. Figure 11e shows that the core-cloak objects only occupy around less than 10% of all the identified objects (objects with both strong and weak drafts and objects with only weak drafts). The ratio of core-cloak updrafts has the maximum near cloud base while that of corecloak downdraft maximizes near cloud top. This indicates that the core-cloak structures originate from the FIG. 10. A schematic diagram of the core-cloak representation of convection. Both updrafts and downdrafts are represented as the combination of a strong core (su, sd) at the center and weak cloak (wu, wd) around the center.

JUNE 2020
G U E T A L .
subcloud layer for updraft and from cloud top for downdrafts. But these convective cloudy objects with core-cloak structure contribute most of the vertical transport of heat and moisture and the vertical fluxes associated with them are very close to the mass flux contribution from the two-draft calculation that comprises all the isolated weak drafts (blue dashed lines in Figs. 11a-d). The core-cloak updrafts (gray lines) dominate the transport throughout most of the cloud layer (except in the lower part of the cloud layer) while the core-cloak downdrafts (gray dashed lines) highlight their importance near cloud top. The dominance of core-cloak drafts on vertical transport has also been confirmed in the RCE simulation (not shown). These results suggest that our core-cloak conceptual model is a true realization of the convective objects that are responsible for the vertical transport of scalars, not only for a bulk description of the weak and strong drafts, but also for individual convective elements.

Discussion and further tests
There is weak or no convective organization in the BOMEX and RCE simulations, and one may ask whether the proposed core-cloak representation could also provide a good description of the fluxes in a situation of organized convection. Becker et al. (2018) showed that a self-aggregated state can result in enhanced horizontal turbulent mixing, and plausibly this may affect the level of inhomogeneity within cloud and hence the intraobject variability. We have therefore extended our analysis to two RCE simulations as described in section 2a: one has interactive radiation and produces self-aggregation and the other has homogenized radiation and does not. The turbulent flux profiles were different in the two cases and the core-cloak representation was able to successfully account for those differences. Figure 12 shows that the bulk mass flux approximation based on traditional cloudy sampling significantly underestimates the vertical fluxes of heat within the whole troposphere and the vertical fluxes of moisture from mid-to high levels. Our core-cloak representation can well capture the vertical transport of heat and moisture both in magnitude and vertical distribution throughout the free troposphere.
In section 5c, we showed that the core-cloak representation is able to account for the turbulent fluxes in the cloud layer. In contrast, within the subcloud layer, the environmental variability is dominant (Figs. 6c,f and 7c,f). FIG. 11. Vertical profiles of time-(last 1 h) and domain-averaged vertical fluxes in BOMEX simulation for (a) the potential temperature flux hw*u * i, (b) the buoyancy flux hw*u y * i, (c) the total water flux hw*q t *i, and (d) the liquid water potential temperature flux hw*u l * i. The red line represents the total resolved flux. The blue solid line represents the mass flux approximation using two-draft representation. The top 0.5% of vertical velocities are used for strong drafts and the bin from 0.5% to 5% for weak drafts. The gray solid line represents the vertical fluxes associated with updrafts that have both core and cloak structure and are also cloudy. The gray dash-dotted line represents the vertical fluxes associated with downdrafts that have both core and cloak structure and are also cloudy. The blue dash-dotted line is the sum of gray solid line and dash-dotted line. (e) Percentage of cloud objects that have both core and cloak structures with respect to all the identified objects in BOMEX simulation: updraft (black solid line) and downdraft (black dashed line).
It would therefore be natural to envisage the use of a boundary layer parameterization within the subcloud layer alongside a core-cloak convection parameterization. Another possibility would be to make an extension to a three-draft model that also treats the nonlocal transport in the subcloud layer. Such a possibility is demonstrated in Fig. 13. If an additional plume type is included to cover the top 5%-40% of vertical velocities, then the mass flux approximation can represent well the resolved buoyancy and total water fluxes in the BOMEX case, not only in the cloud layer but also in the subcloud layer (Figs. 13a,c). The improvement in the subcloud layer is due to the mass flux contribution from both the weak updrafts and downdrafts (Figs. 13b,e). This suggests that a single updraft and downdraft may be enough for the transport in the subcloud layer. The fact that the strongest drafts do not play a major role in subcloud vertical fluxes is due to the less skewed distribution of vertical velocity in the subcloud layer than in the cloud layer for both shallow and deep convection (see Fig. S2 in supplemental materials).
Our results provide some support for extensions of EDMF schemes. While the original formulation of EDMF uses a single updraft (Siebesma and Teixeira 2000;Soares et al. 2004;Siebesma et al. 2007), it could also include multiple draft types. This idea has been tested using dual (Neggers et al. 2009;Neggers 2009) or multiple mass flux schemes (Cheinet 2003(Cheinet , 2004Su selj et al. 2012). Our study further emphasizes the important contribution from cloud-top downdrafts to the heat fluxes in the inversion layer. This has also been confirmed in a detailed study on the cloud life cycle (Zhao and Austin 2005) and the coherent structures (Park et al. 2016) in shallow cumulus clouds and also in stratocumulus clouds (Davini et al. 2017;Brient et al. 2019). Knowing the key physical processes for vertical FIG. 12. Vertical profiles of time-(last 5 days) and domain-averaged resolved vertical turbulent fluxes in the selfaggregation simulation: (a) hw*u * i, (b) hw*u y * i, (c) hw*q t * i, and (d) hw*u l * i. The red line represents the total flux, the blue solid line represents the mass flux approximation using core-cloak decomposition [term (7.1)], the green solid line represents the intradraft variability [term (7.2)], and the blue dashed line represents the bulk mass flux approximation [term (4.1)] using bulk cloud-environment decomposition based on traditional cloud sampling. transport throughout the cloud layer will provide valuable guide for future development of EDMF schemes, given that previous studies did not consider the downdrafts. Following this idea, the EDMF scheme could also be extended to include deep convection. Some recent studies have also suggested more general extensions of EDMF with multiple drafts (Thuburn et al. 2018;Tan et al. 2018).
Other recent studies have also argued that the description of convective clouds with only a bulk updraft or downdraft is inadequate (Heus and Jonker 2008;Hannah 2017). The core-cloak representation of convection in this study shares some similarity with other proposals but also differs from them in important ways. The three-layer model of Heus and Jonker (2008) divides the flow into the cloud core with positive velocity and buoyancy, the subsiding shell structure wrapping around the core with negative vertical velocity and buoyancy, and the environment. The buffered-core model of Hannah (2017) proposes a core in the center, the environment in the outmost region and a buffered region composed of a mixture of detrained core air and entrained environmental air. Our core-cloak representation treats both the updrafts and downdrafts as having a core of strong draft surrounded by a weak draft and does not require a particular sign for the buoyancy. We should stress, however, that the core-cloak model provides simply a possible decomposition of the flow that gives an accurate and efficient description of turbulent fluxes using a mass flux approximation. To implement our model as a full parameterization scheme would of course require the development of descriptions of triggering, closure and the exchange terms between weak and strong drafts and the environment.

Summary
The bulk mass flux formulation has been evaluated for both shallow and deep convection using large-eddy simulation data. It is found that the bulk mass flux approximation cannot capture the right magnitudes and vertical distributions of turbulent heat and water fluxes at FIG. 13. As in Fig. 8), but for a three-draft representation (and without the gray lines for the bulk mass flux approximation based on traditional cloud sampling). The top 0.5% of vertical velocities are used for strong drafts, a bin from 0.5% to 5% for medium drafts, and a bin from 5% to 40% for weak drafts. In (b), (c), (e), and (f), the contributions are shown for weak updrafts (wu, yellow), medium updrafts (mu, red), strong updrafts (su, magenta), weak downdrafts (wd, black), medium downdrafts (md, green), strong downdrafts (sd, blue), and the environment (env, cyan). the same time using a cloud-environment decomposition. A bulk mass flux approximation neglects contributions that arise from inter-and intraobject variability. The interand intraobject variabilities of the turbulent heat flux are comparable in magnitude to the estimate from the bulk mass flux approximation and do not share similar shapes. Hence, they cannot be parameterized through a rescaling method. In addition, a large part of the subplume fluxes is associated with the buoyant core of clouds and therefore cannot be represented through a downgradient assumption as applied in Lappen and Randall (2001). A spectral analysis emphasizes the comparable contribution of intraobject variability and the mass flux term to the total fluxes across the whole spectrum of cloud size, although interobject variability can be much reduced in such a representation. The above results show that there are nonnegligible contributions to the fluxes from fine structures within and outside the cloud, which are ignored by the bulk mass flux approximation.
To understand the key elements of cloud circulations responsible for the turbulent transport, a decomposition based on the distribution of vertical velocity was used to construct different types of drafts. The decomposition using a single bulk updraft and its environment substantially underestimate the fluxes of thermodynamic quantities using the bulk mass flux approximation, consistent with previous studies. With a single downdraft also included, the bulk mass flux approximation improves near cloud top in shallow convection but still underestimates the fluxes in the lower part of the cloud and in the subcloud layer. The downdraft motions produced in response to the overshooting updrafts near cloud top do contribute substantially to the vertical fluxes in the inversion layer and should be parameterized in the shallow convection scheme. There are important contributions to fluxes in the lower part of the cloud layer from the less extreme updrafts. This indicates that the vertical transport is controlled by a combination of drafts of different strengths. Accordingly, we proposed a ''core-cloak'' conceptual model for both updrafts and downdrafts. It is found that such a core-cloak representation can well capture the vertical fluxes with a mass flux approximation in terms of both the magnitudes and vertical distributions. It improves the mass flux approximation of both heat and water fluxes significantly (50% at some levels) for both shallow and deep convection. Therefore, this study shows that a simple minimal extension of the bulk mass flux framework would be sufficient to correct the underestimation of turbulent fluxes without the need for more complicated parameterizations of intraobject variability. We intend to pursue the practical implications of this conceptual model within the future development of a convection parameterization.