Composited structure of non-precipitating shallow cumulus clouds

The normalized distributions of thermodynamic and dynamical variables both within and outside shallow clouds are investigated through a composite algorithm using large-eddy simulations of oceanic and continental cases. The normalized magnitude is maximum near the cloud centre and decreases outwards. While relative humidity (RH) and cloud liquid water ( q l ) decrease smoothly to match the environment, the vertical velocity, virtual potential temperature ( 𝜃 v ), and potential temperature ( 𝜃 ) perturbations have more complicated behaviour towards the cloud boundary. Below the inversion layer, 𝜃 ′ v becomes negative before the vertical velocity has turned from an updraft to a subsiding shell outside the cloud, indicating the presence of a transition zone where the updraft is negatively buoyant. Due to the downdraft outside the cloud and enhanced horizontal turbulent mixing across the edge, the normalized turbulent kinetic energy (TKE) and horizontal turbulent kinetic energy (HTKE) decrease more slowly from the cloud centre outwards than the thermodynamic variables. The distributions all present asymmetric structures in response to the vertical wind shear, with more negatively buoyant air, stronger downdrafts, and larger TKE on the downshear side. We discuss several implications


INTRODUCTION
Atmospheric convective mixing was found to account for half the spread of climate sensitivity in a study of CMIP5 global climate models (Sherwood et al., 2014). Bulk plume models have been widely used within convection parameterization schemes to estimate the vertical transport of heat, moisture, and momentum produced by a large ensemble of clouds that cannot be resolved explicitly within a climate model grid box. Such plume models apply the top-hat assumption (Randall et al., 1992;Yano et al., 2010), which neglects the variability of vertical velocity and transported variables, and which can lead to underestimation of vertical transport (Siebesma and Cuijpers, 1995;Guichard et al., 1997;Yano et al., 2004). A more complete description may be provided by considering multiple plumes with different sizes, entrainment and detrainment rates etc. Such spectral plume models (e.g., Arakawa and Schubert, 1974;Zhang and MacFarlane, 1995) attempt to account for interobject variability due to differences in mean properties between different types of convective cloud. However, the approach still neglects the inhomogeneity within each cloud, that is, the intraobject variability, which has been shown to contribute a non-negligible proportion of vertical heat flux (Gu et al., 2020a). An appropriate representation of intraobject variability depends on a detailed understanding of the internal structures within clouds. Previous studies of aircraft observations and large-eddy simulations have taken the first steps towards a fuller understanding of the cloud margin and the subsiding shells around shallow convective cloud. The near-cloud shells are mostly driven by evaporative cooling due to turbulent mixing between the cloud and the environment, rather than mechanical forcing (Rodts et al., 2003;Heus and Jonker, 2008;Wang et al., 2009;Wang and Geerts, 2010). They also compensate for a large part of the upwards mass flux within the cloud Heus et al., 2009;Glenn and Krueger, 2014). While these studies have demonstrated that variables are inhomogeneously distributed within the cloud rather than following a top-hat distribution, there is no consensus on alternative forms of the distributions. Some observations using aircraft and radar wind profilers suggest that the vertical velocity has a "triangular" shape (Zipser and LeMone, 1980;Wang et al., 2020), while a recent convection scheme has used a parabolic distribution (Leger et al., 2019). It is also important to determine the differences and similarities between the distributions of different variables, since these will affect the fluxes and may give further insights into the dynamics and thermodynamics within the cloud. Knowing the differences and similarities can also provide useful information for parameterizations that use assumed distributions. Composite structures obtained from observations have focused on larger cloud sizes (e.g., > 400 m), partly to ensure sufficient sampling within each cloud. However, the structures of medium-and small-sized clouds are also important, because the medium-sized clouds dominate the total vertical transport (Gu et al., 2020a), whereas the structure of small clouds might provide useful information on how they are inhibited from developing and also improve our understanding of cloud life cycles. In addition, small cumulus clouds can also help to sustain larger clouds by detraining heat and moisture above the cloud base and thus maintaining cloud-subcloud coupling (Neggers, 2015). Another caveat is that intersections taken during flights are not always through the cloud centre, which can cause difficulties in determining size distributions (e.g., Barron et al., 2020) and may bias the construction of composites.
In this study, we undertake a detailed investigation of the spatial distributions of thermodynamic and dynamical variables within and near the edge of numerically simulated shallow cumulus clouds. The large-eddy simulation and the algorithm for compositing the cloud structure are described in Sections 2.1 and 2.2, respectively. The composites themselves are discussed in Section 3, focusing on the distributions of key thermodynamic and dynamical variables (Section 3.1), their dependence on cloud size and variations in the vertical (Section 3.2), and possible analytic forms of the distributions that can be used within parameterization schemes (Section 3.3). These distributions for individual variables are then combined in Section 4 to examine whether products of the distributions can be used to reproduce the vertical fluxes of moisture and heat. Discussion and a summary are given in Section 5.

Large-eddy simulation of shallow cumulus clouds
The Met Office-NERC (Natural Environment Research Council) Cloud (MONC: Brown et al., 2015; model is used to perform a large-eddy simulation of oceanic shallow convection based on the Barbados Oceanographic and Meteorological Experiment (BOMEX). Most of the model configuration follows that of Siebesma et al. (2003) with slight changes (Gu et al., 2020a). The grid spacing used is 25 m in all directions and the domain size (L × W × H) is (15 × 15 × 3) km 3 . 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 (a) (b) F I G U R E 1 Schematic diagram of the composite algorithm. The grey patch in (a) denotes a horizontal cross-section of a sample cloud to be analyzed. The black dashed lines represent the simulation grid. The red dot is the cloud-liquid-water weighted centroid. Four neighbouring grid lines highlighted in blue and cyan are taken as sampling slices for the composite. The thick cyan line is taken as an example slice to illustrate the normalization in (b). Midpoint and cloud edge points of the slice are denoted with green and yellow dots, respectively. The distance x along the slice is measured from the midpoint and normalized by the length between the two cloud edge points, 2R. The composited variable is interpolated on to equally spaced points along the slice (eight points within the cloud and four points outside) and is normalized by its maximum value m within the full cloud object (the maximum value may not necessarily lie on the slice) [Colour figure can be viewed at wileyonlinelibrary.com] between water vapour and cloud liquid water, as this is a nonprecipitating case without ice water. More details of the simulation can be found in Gu et al. (2020a). Our analyses cover a period in the equilibrium state at hour 5 of the simulation, with 10-min output frequency. The composited cloud structure in the BOMEX case was also checked with another large eddy model, the CM1 model (Bryan and Fritsch, 2002) developed at the National Center for Atmospheric Science (NCAR), using the same resolution, and the conclusions do not change.
Another case was also studied in order to test the robustness of the composited structures of nonprecipitating shallow cumulus clouds. That simulation is based on observations at the Southern Great Plains (SGP) site of the Atmospheric Radiation Measurement (ARM) Program on June 21, 1997, and follows the configuration described by Brown et al. (2002). The model grid spacing is the same as in the BOMEX simulation and the domain size is (6.4 × 6.4 × 4) km 3 . The subgrid turbulence was again parameterized with the Smagorinsky-Lilly scheme, while the microphysics is parameterized using Cloud AeroSol Interaction Microphysics (CASIM: Grosvenor et al., 2017;Miltenberger et al., 2018) in a double-moment configuration. The ARM simulation runs for 24 hr with 15-min output frequency and our analysis focuses on a 6-hr period centred at 1900 UTC, when the surface sensible and latent heat fluxes achieve their maximum values of 140 and 500 W⋅m −2 , respectively. The general features of composited structures in these oceanic and continental nonprecipitating shallow cumulus clouds are consistent from cloud base to cloud top (not shown) and therefore only the results from the BOMEX case are presented here.

2.2
Composite algorithm Figure 1 shows a schematic diagram of the algorithm for compositing the cloud structure. At each vertical level, all cloudy points are first identified with the cloud liquid water criterion q l > 10 −5 kg⋅kg −1 . Contiguous labelled points are combined to form an individual cloud object by checking the eight neighbouring grid points around the cloudy points until no more cloudy points are found. For each cloud object, the liquid water weighted centroid is then determined as the cloud centre. Four neighbouring model lines (two along the west-east direction and another two along the south-north direction) nearest to the centroid are taken as sample slices for the cloud structure composite. Any slices that cut across multiple boundaries of cloud objects are not used in the composite. This may neglect some organized clouds with multiple cores. However, these clouds occupy a small portion (about 5%) and thus we can obtain good sampling for compositing the cloud structures. In each slice, the midpoint of the intersection between the slice and the cloud edge is taken as the centre for the composite. The slice is then divided into 12 equal-length bins, of which eight bins are fully within the cloud and a further two on each side represent its immediate environment. Extra tests on the number of bins indicate consistent conclusions. Distances to the centre are normalized by the length of intersection, 2R. The values of a particular variable for each bin are linearly interpolated from the model output and are normalized by its maximum absolute value within the cloud object. Therefore, normalized values are usually less than 1, as the maximum value may not always be located in the slice. After repeating the above process for all cloud objects at a particular level, the data from every slice are composited to obtain the normalized distribution of different variables from the cloud centre outwards and across the cloud boundaries.
Through this compositing algorithm, we are able to determine how different regions (e.g., cloud core, transition zone, subsiding shell) are configured relative to the cloud edge. We remark that the cloud edge here is defined as the outermost cloudy boundary of each cloud object, rather than the innermost cloud-free boundary outside the cloud, as defined in Heus and Jonker (2008). While there is no established "self-similarity" theory for the scaling of cumulus cloud properties, our rescaling of spatial coordinates and normalizing by maximum perturbations has a practical advantage to aggregate some similar features of distributions within clouds with different sizes and geometries, and is convenient for a direct comparison of the distributions of different thermodynamic and dynamical variables. The way that we take the slices also has the advantage of capturing asymmetric structures with respect to the vertical wind shear (see Section 3), which is approximately from west-southwest to east-northeast within the cloud layer during the equilibrium state of our simulation.
To confirm the robustness of the composited structures, we have also tested some alternative algorithms: for example, by starting from the cloud edge and shrinking the cloud boundary inward or extending the cloud boundary outward in single grid-box steps to assess the internal or near-environment structures, respectively. We first identify the boundary of each cloud object, and then move the boundary outward for one grid box each time. The mean properties of each boundary can be calculated and composited to get the distribution outside the cloud. For the distribution within the cloud, we move the cloud boundary inward for one grid box each time until no more grid boxes are within the boundary. If there is an isolated grid box within the innermost boundary, that grid box will be considered as the cloud centre. If not, the innermost boundary will be viewed as the cloud centre. The mean properties of each boundary are calculated and interpolated further on to the normalized distance to the cloud centre. The resulting composited structures are consistent with results to be presented below from the final algorithm chosen (not shown). Sensitivity tests on the model simulation grid length (10, 25, 50, and 100 m) also give consistent results for the main features of the composited structures (not shown).

3.1
Overall structure Figure 2 shows the composited structure of both thermodynamic (Figure 2a-c) and dynamical (Figure 2c-e) variables for all cloud objects at vertical levels of 600, 1,000, and 1,800 m. The simulation has a well-mixed layer below 520 m followed by a conditionally stable layer up to 1,500 m, and an inversion layer from 1,500-2,000 m, where the cloud numbers gradually cease. The chosen levels are near the cloud base, around the middle of the cloud layer, and at the cloud top.
The normalized values of all variables except the relative humidity are less than 1.0, indicating that the maximum values are not always along the cloud slices. The broad features of the mean normalized distributions for different variables and levels are consistent, with peak amplitudes occurring near the cloud centre and decreasing outwards, indicating that the cloud core is less affected by mixing with the environment. However, different variables have different distribution shapes in the vicinity of the cloud edges.
As expected, the relative humidity is almost homogeneously distributed within the cloud and normalized values are 1.0 in most parts. As a result, the variability between slices of the relative humidity is very small within the cloud. Outside the cloud, the variability increases with distance away from the cloud, indicating large fluctuations in the near-cloud environment. The averaged distribution of normalized liquid water mixing ratio decreases quickly to zero outside the cloud, while the mean relative humidity decreases more slowly to match the environment, forming a moist buffering region that likely helps prevent the cloud from being subject to direct mixing with dry environmental air. The region itself is likely to result in part from mixing between the cloud and the environment. The size of this moist patch is not much larger than the size of the cloud at the midlevel of the cloud layer ( Figure 2b) and near the cloud top ( Figure 2c). Near the cloud base, however, the size of this region is larger than the cloud size, indicating that cloud triggering within large moist patches is favoured (Figure 2a: Park et al., 2016).
The presence of a buffering region mediates the exchange between cloud and environment in the near-cloud region. Bulk entraining plume models describe the lateral mixing in terms of entrained air bringing the mean property of the environment and detrained air F I G U R E 2 Composited normalized distributions of (a-c) thermodynamic and (d-f) dynamical variables at heights of (a,d) 600 m, (b,e) 1,000 m, and (c,f) 1,800 m. The thermodynamic variables are relative humidity perturbation (RH ′ , blue line), virtual potential temperature perturbation ( ′ v , red line), potential temperature perturbation ( ′ , brown line), and cloud liquid water perturbation (q ′ l , grey line). We show the perturbation of relative humidity because the environmental relative humidity varies at different vertical levels. The dynamical variables are perturbations of total turbulent kinetic energy (TKE, purple line), horizontal turbulent kinetic energy (HTKE, yellow line), and vertical velocity (w ′ , black line). The normalized distributions of cloud liquid water q ′ l and vertical velocity w ′ are included on each panel for ease of comparison. Solid lines represent the mean distributions and the light shadings with corresponding colours show one standard deviation around the mean (+∕− one standard deviation). The number of slices used for the compositing is indicated in the header of each panel. On the x-axis, "C" represents the centre of the composite slice and "L" denotes the length between cloud edge points (2R in Figure 1) [Colour figure can be viewed at wileyonlinelibrary.com] taking the mean property of the bulk plume. However, different direct measures of entrainment and detrainment show that the air mixed during the exchanges takes the values of variables within the buffering region, rather than in the further environment or cloud core (Romps, 2010;Dawe and Austin, 2011a;2011b). Therefore, the difference between the properties of entrained and detrained air is smaller than that between the mean properties of cloud and the far environment, resulting in entrainment and detrainment rates larger by roughly a factor of 2 than in the bulk plume approach.
The strongest upward motion occurs near the cloud centre, but has smaller normalized vertical velocity than the liquid water content. At the cloud base ( Figure 2a) and in the middle of the cloud layer ( Figure 2b) the cloud edge exhibits some weak upward motions. Outside the cloud, the vertical velocity becomes negative, indicating the presence of a subsiding shell. However, the variability of vertical velocity outside the cloud suggests that not all clouds have distinct shell structures. Near the cloud top, the cloud edge marks the approximate boundary between upward and downward motions (Figure 2c). The downdrafts outside the cloud are most pronounced near the cloud top and also cover the largest normalized size, compared with those in the middle of the cloud layer and at the cloud base. These features suggest that the subsiding shell may originate near the cloud top and that it weakens, in the normalized sense, during descent.
The buoyancy near the cloud base ( Figure 2a) and in the middle of the cloud layer ( Figure 2b) decreases from its peak value near the centre and changes sign some distance before the cloud edge. This decay is more rapid than that of the vertical velocity, which means that there is a transition zone where there is upward motion with negative buoyancy. Beyond that lies the region of subsiding shells outside the cloud, which is negatively buoyant, presumably due to evaporative cooling induced by turbulent mixing as found in previous studies (Rodts et al., 2003;Heus and Jonker, 2008;Wang et al., 2009;Wang and Geerts, 2010).
The buoyancy composite does lend some support to the concept of buoyancy-sorting mechanisms used in some parameterizations (e.g., Kain and Fritsch, 1990). The central idea is that air with negative buoyancy is detrained outside the cloud, while air with positive buoyancy is entrained into the cloud. However, the existence of a transition zone with negative buoyancy but upward motion indicates that not all air with negative buoyancy should be detrained out of the updraft and a fraction is retained within the cloud. This finding is consistent with previous observational and simulation studies (Taylor and Baker, 1991;Siebesma and Cuijpers, 1995;Zhao and Austin, 2005), and has also been applied in the shallow convection scheme of Bretherton et al. (2004).
The normalized ′ has a similar distribution to the buoyancy, but its amplitude is quite different. ′ in the centre of the cloud is mostly positive in the middle of the cloud layer (also see Section 3.2), but negative near the cloud base and cloud top (Figure 2a-c). Near the cloud base, the opposite signs of buoyancy and ′ emphasize the importance of water-vapour anomalies in cloud triggering (Fabry, 2006;Madaus and Hakim, 2016) and contrast with the use of only a temperature perturbation in the triggering tests for some parameterization schemes (e.g., Gregory and Rowntree, 1990;Kain, 2004). At the cloud top, there is no positive buoyancy within the cloud. The normalized ′ at the midlevel and cloud top has distribution and amplitude very similar to the buoyancy, despite some variability between slices, presumably because the virtual effect of anomalous water vapour is cancelled by that of liquid water. Both ′ and ′ v have a local maximum near the cloud centre at the cloud top, despite their negative amplitude. This type of distribution supports continuous generation of horizontal vorticity with rising thermals and results in an increasing dynamical pressure drag with height for clouds, as is shown in Gu et al. (2020b).
The normalized distributions of turbulent kinetic energy (TKE) and its horizontal component (HTKE) decay outward more slowly than those of the thermodynamic variables. Within the cloud, the TKE is dominated by its vertical component, which accounts for more than 80% of the total near the cloud base ( Figure 2d) and in the middle of the cloud layer (Figure 2e), and more than 60% near the cloud top ( Figure 2f). However, the HTKE becomes important outside the cloud. The slow decay of TKE across the cloud edge while the vertical motion decreases more rapidly indicates enhanced horizontal turbulent mixing between the cloud and the environment.
The distributions of all the variables have a somewhat asymmetric structure with respect to the vertical wind shear (which is directed from left to right in Figure 2). The negative buoyancy is slightly stronger on the downshear side, resulting in a stronger subsiding shell. As a result, the distributions of TKE and HTKE also have higher TKE and HTKE on the downshear side over the range of heights for which there is large-scale vertical wind shear (above 700 m). Note that near the cloud top (Figure 2c,f) the buoyancy and ′ distributions have rather weak asymmetry, but the downdraft is nonetheless stronger on the downshear side. This may suggest that mechanical forcing plays some role in favouring downdraft on the downshear side near cloud top.
It should be noted that the variability of the normalized distributions between slices is larger for vertical velocity, buoyancy, and potential temperature than for cloud liquid water and relative humidity. This indicates that the vertical heat fluxes may be more difficult to estimate than the vertical moisture fluxes if only the averaged distributions of thermodynamic and dynamical variables are used, without taking into account the variations between cloud slices. This will be demonstrated in further detail in Section 4.

Dependence on cloud size
The general pattern of the normalized distribution for clouds with different sizes does not change a great deal compared with the composite of all clouds, a point that provides some a posterori justification of the normalizations used. Nonetheless, there are some interesting differences in some of the details, as shown in Figure 3. For example, the moist buffering region denoted by the slowly decaying relative humidity has a larger normalized size for smaller clouds, indicating that the absolute size of buffer region may not scale linearly with the cloud size in a simple way. The feature that the buffering region is larger at the cloud base (Figure 3a,d,g) than above (Figure 3b,c,e,f,h,i) is consistent for clouds with different sizes. For small clouds, the normalized positive buoyancy (Figure 3a,b) within a cloud is relatively weak compared with larger clouds (Figure 3d,e,g,h). However, the normalized ′ remains negative in all of the small cloud composites, showing the important role of water vapour in maintaining the buoyancy of these clouds. For medium and large clouds, the normalized ′ is much closer to the buoyancy in the middle of the cloud layer, indicating that the buoyancy of these clouds is largely latent-heating F I G U R E 4 As in Figure 3 but for dynamical variables, presented in the same way as in Figure 2d-f [Colour figure can be viewed at wileyonlinelibrary.com] driven. At the cloud top, the buoyancy and ′ are negative, but with somewhat less buoyant air near the cloud centre for all cloud sizes.
For dynamical variables (Figure 4), the normalized downdraft outside the cloud is strongest for small clouds. The strong downdraft results in secondary maxima at the cloud top in the TKE distribution (Figure 4c). Also note that the transition zone with negatively buoyant updraft has a larger normalized size in smaller clouds (Figure 3). Within the cloud, the HTKE accounts for almost 50% of the total TKE near the cloud top for small clouds (Figure 4c), but remains below 20% for medium and large clouds at all F I G U R E 5 Scatter plots of mean vertical velocity (w, y-axis) and the maximum vertical velocity (w max , x-axis) within each cloud object at (a) 600 m, (b) 1,000 m, and (c) 1,800 m. Red solid dots are data points from smaller clouds (0-50% of the size distribution) and blue open circles are data points from larger clouds (50-100% of the size distribution). Red and blue dashed lines represent the least-squares regression for smaller and larger sizes, respectively. The equations of the least-squares fits are given on the legend with corresponding colours. The black solid line represents the relationship w = 0.5w max , which has been assumed in some previous theoretical studies [Colour figure can be viewed at wileyonlinelibrary.com] levels ( Figure 4d-i). Asymmetric distributions with respect to the weak shear exist for all cloud sizes, with relatively stronger downdrafts, TKE, and more negative buoyancy on the downshear side. However, the asymmetry is most pronounced from the small clouds. The enhanced strength of the downdraft with height also holds for all cloud sizes, suggesting a potentially important role for the downdraft in vertical transport, which will be discussed in Section 4.

Possible power-law distributions
It has been shown in Sections 3.1 and 3.2 that the normalized distributions of buoyancy and ′ have some dependence on cloud size and height in the vertical, although both variables consistently produce a local maximum near the cloud centre. Relative humidity is almost homogeneously distributed within the cloud, due to the saturation, and therefore may not be able to provide useful information on in-cloud dynamics. However, the distributions of cloud liquid water and vertical velocity demonstrate generally consistent features regardless of cloud size and vertical level. Therefore, it is worth considering whether one can describe their inhomogeneity within the clouds using an explicit form of the normalized distribution function. Such a distribution would have benefits for both scientific understanding of cloud dynamics and the development of novel parameterizations of convection. For example, a radial dependence of vertical velocity within a cloud has been taken into account in the parameterization of Leger et al. (2019), motivated by the representation of pressure drag. A recent theoretical study (Yano, 2020) also suggested that the inhomogeneity of vertical velocity must be introduced in plume models, otherwise the buoyancy and entrainment effect will perfectly cancel out with the counterbalancing force from the pressure perturbation, leaving a pure drag force, which prevents a steady-state solution. In this subsection, we focus on the normalized distribution of vertical velocity, as it is the key to coupling cloud dynamics, microphysical processes, and other processes (e.g., aerosol loading and radiative effects). Inspection of the composited structure in Figure 2 suggests that a simple power-law distribution might be suitable within the cloud, as in Leger et al. (2019). We denote the normalized distribution as f (r∕R), where r is the distance from the cloud centre and R the cloud radius. A power-law distribution takes the form f (r∕R) = a 0 + a 1 (r∕R) m , where a 0 and a 1 are constants that can be chosen to satisfy the conditions at r = 0 and r = R, respectively, and m is the power of the distribution. We need an extra constraint to determine the m parameter, and a natural choice is to set the first-order moment, which is the mean vertical velocity. Figure 5 shows the mean and the maximum vertical velocity within each cloud object at the cloud base (Figure 5a), the middle of the cloud layer (Figure 5b), and the cloud top ( Figure 5c). The data points lie approximately along a line on which the mean vertical velocity is half the maximum value. The ratio between the mean vertical velocity w and maximum vertical velocity w m is larger than 0.5 for small clouds (red solid dots) and lower than 0.5 for large clouds (blue open circles) below the cloud top (Figure 5a,b). This relationship also holds for the shallow clouds of the ARM simulation (not shown) and in a smaller-domain higher-resolution BOMEX simulation (10-m grid spacing; not shown). Note that the definition of "small" and "large" clouds is slightly different from that in Sections 3.1 and 3.2 (See captions of Figure 5). The ratio is about 0.6 at the cloud top, where it does not differ significantly between large and small clouds. Such a linear relationship between w and w m has been used as a prior assumption in previous theoretical studies on cloud dynamics (Kuo and Raymond, 1980;Morrison, 2016) and in a new parameterization scheme (Peters et al., 2021), but has not been thoroughly validated. A recent observational study (Wang et al., 2020) found a similar relationship with a ratio between 0.5 and 0.8 that decreases systematically with updraft size. Our study provides further evidence to support this linear hypothesis using high-resolution simulations of shallow cumulus clouds. Despite the spread among clouds with different sizes, the 0.5 ratio between w and w m provides a reasonable approximation for the cloud ensemble.
We are now able to determine a suitable power for the normalized distribution. The cloud object is assumed to have a symmetric round shape and thus the vertical velocity w at any location within the cloud can be written as w(r) = w m f (r∕R). For simplicity, we assume that the radial distribution f (r∕R) is independent of height. Figures 2-4 do indicate some height variations, but it may nevertheless be a reasonable first approximation. Therefore, the mean vertical velocity of the cloud object can be obtained as follows.
(1) For a 2D symmetric cloud, Using to denote the ratio of w to w m , we have (2) Similarly, for a 3D axisymmetric cloud, at each vertical level, we have Substituting the power-law normalized distribution into Equations 2 and 3, we can relate and m in a 2D cloud, and in a 3D cloud, If we set a 0 = 1 and a 1 = −1 to make sure that f (0) = 1.0 at the cloud centre and f (R) = 0 at the cloud edge, the F I G U R E 6 The relationship between the power m of the assumed normalized distribution and the ratio between the mean and maximum vertical velocity within the cloud ( = w∕w max ) for axisymmetric 2D (black) and 3D (blue) cloud using Equations 6 and 7 [Colour figure can be viewed at wileyonlinelibrary.com] equations can be written as and = m m + 2 (3D cloud). (7) Figure 6 gives the ratio in terms of the power m. For large m, approaches 1 slowly and the distribution gradually tends to the top-hat distribution in the limit of m→∞. According to Figure 6, the observational value of between 0.5 and 0.8 implies that the power of the distribution should be around 1-4. If we choose = 0.5 as in the theoretical studies, then the power m = 1 for the 2D cloud and m = 2 for the 3D cloud. The linear m = 1 distribution is consistent with a "triangular" shape of vertical velocity, which has been observed in aircraft transects or time-height sections through deep clouds (e.g., Zipser and LeMone, 1980;Wang et al., 2020). Therefore the 2D-like cloud structure might arise as a result of observational strategy or different cloud types. Our result also supports a recent convection parameterization that takes into account internal structures of 3D axisymmetric updrafts and uses a quadratic function to represent the vertical velocity distribution (Leger et al., 2019).

A SIMPLE APPLICATION OF THE NORMALIZED DISTRIBUTIONS
In this section, we attempt to use the composited normalized distributions of vertical velocity and thermodynamic variables from Section 3 to estimate the vertical fluxes of heat and moisture. Two key elements need to be considered: the maxima and the normalized distributions for different variables. Rather than discussing the form of normalized distribution, as in Section 3.3, we are concerned here about the extent to which we can reproduce the vertical fluxes using composited distributions.

Formulation
A straightforward way to explore whether the normalized distributions of quantities within clouds could provide an accurate representation of vertical fluxes is to calculate the fluxes for each cloud object separately and then sum over all cloud objects. Using the normalized distributions in this way, an assumption is needed that the shape of each cloud object can be approximated as round, thereby enabling the product of distributions to be straightforwardly integrated over the cloud object. The analysis requires the maximum perturbations within each cloud object in order to construct the summation, as well as each of the normalized distributions. Hence, it would not be appropriate as a practical basis for a parameterization approach. However, this estimation provides a baseline for investigation. If it can provide a good representation, then it is possible to seek various simplifications. Comparison with the simplifications can also help to indicate those aspects of an assumed-distribution formulation that are critical for an adequate representation of vertical fluxes. We first introduce some notations for convenience of discussion. An atmospheric quantity within a cloud object i is i , and the maximum perturbations of vertical velocity w and transported variable within this object are w ′ mi and ′ mi , respectively. The normalized distributions for vertical velocity and within each cloud object are f w i (r∕r i ) and f i (r∕r i ), where r i is the equivalent radius of cloud object i and is calculated as √ S i ∕ , where S i is the area coverage of cloud object i. The contribution from cloud object i to the domain-averaged vertical flux of can then be calculated as where S tot is the total area of the domain. Summing over all cloud objects across the domain gives the total domain-averaged vertical fluxes: Thus, any contributions from the environment are neglected.
As noted above, further simplification would be necessary for practical use. First, instead of using the radial distribution for each cloud object, we could use the averaged normalized distribution over all of the cloud objects at each vertical level. The flux estimation would then become where f w i (r∕r i ) and f i (r∕r i ) are the averaged normalized distributions of all cloud objects, as shown in Figure 2 and discussed in Section 3.
Alternatively, we might try to use the averages over all cloud objects of the maximum perturbations w ′ m and ′ m at each level, leading to (11) We can further use both the averaged maximum perturbations and the averaged normalized distribution to estimate the vertical fluxes as Figure 7 shows the results of different estimations of vertical heat and moisture fluxes, together with the true vertical fluxes (black dash-dotted line). Using the maximum perturbations and normalized distributions of each cloud object (Equation 9, black solid line) gives a remarkably accurate estimation for all of the in-cloud vertical fluxes (black dash-dotted line), in terms of not only the magnitude but also the vertical distribution. This implies that the round shape assumption may be appropriate for shallow cumulus clouds. It also demonstrates that a good representation of vertical transport by shallow cumulus clouds can be obtained if we know how the variables are distributed radially from the cloud centre. However, the above representation would be difficult to implement for practical use.

Representation of vertical heat and moisture fluxes
One simplification is that of Equation 10, replacing the normalized distribution of each cloud with the averaged normalized distribution of all clouds. However, for the F I G U R E 7 Domain-averaged resolved vertical turbulent fluxes of (a) v , (b) , (c) l , and (d) q t . The total fluxes are marked by grey dash-dotted lines, while the fluxes contributed from all of the clouds are represented with a black dash-dotted line ("cloud" in the legend). The estimations using Equations 9-12 are denoted with black, red, blue, and green solid lines, respectively. . The vertical buoyancy flux is also significantly underestimated, ranging from less than 25% at 800 m to above 50% at 1,600 m. The underestimation of vertical heat fluxes may arise from complicated mixing between the cloud and the environment near the cloud edge, where negative heat fluxes occur. Taking the averaged distribution may result in stronger negative buoyancy or ′ within the transition zone and thus significantly cancel out the positive fluxes contributed by the buoyant cloud core, so that the total heat fluxes are underestimated. In contrast, though, the formulation of Equation 10 does give encouraging results that the vertical fluxes of q ′ t and ′ l are reproduced reasonably well. These estimations are close to those of Equation 9 within most of the cloud layer, albeit with a small underestimation at the cloud top. Unlike the distributions of ′ v and ′ , cloud liquid water q t and liquid water potential temperature ′ l (q t and ′ l are anticorrelated) do not have negative perturbations in the transition zone. As a result, there is less uncertainty in water distribution near the cloud edge, which therefore leads to the good performance of Equation 10 for water fluxes. The results here are a manifestation of the "evaporation-aware" temperature being more distributed in a more complex fashion within the cloud than water.
We can also calculate the vertical fluxes using individual normalized distributions for each cloud, but with the averaged maximum perturbations of all clouds, as in Equation 11 (blue solid line in Figure 7). The result is generally worse than when using Equation 10, since the vertical fluxes of both heat and moisture variables are underestimated. For vertical fluxes of q t and l , the underestimation is about 10% near the cloud base and gradually increases with height to a maximum of about 50% at 1,600 m. The vertical fluxes of v and are both underestimated below the inversion layer. Unlike the estimation with Equation 10, it reproduces the correct sign of the flux, although the magnitude is weak within the cloud layer. This is because the normalized distribution near the cloud edge is dealt with explicitly for each cloud object, confirming the suggestion that a careful treatment of the transition zone is important for parameterizing the vertical heat fluxes. Within the inversion layer, however, the heat and buoyancy fluxes (Figure 7a,b) are positive using this form of estimation, but the actual heat and buoyancy fluxes are negative from 1,700 m and above because of the negative buoyancy and ′ in the updraft (Figure 2c). Using the averaged maximum perturbation results in relatively weak downward motion and thus produces slightly positive heat fluxes. Therefore, the diversity of maximum perturbation values across cloud objects also results in large uncertainty, leading to underestimation of heat and moisture fluxes. This suggests that, in order to estimate the maxima in a parameterization with an assumed distribution, one may require a spectral model, using a range of different cloud sizes and entrainment and detrainment rates.
We can further simplify the estimation of vertical fluxes using Equation 12, with both the maximum perturbations and normalized distributions replaced by their averaged values over all clouds (green solid line in Figure 7). This gives similar results to Equation 11 for vertical fluxes of q t and l and similar results to Equation 10 for the fluxes of v and , confirming that both the normalized distribution and maximum perturbation of each cloud are important. However, this does not mean that a reasonable estimation of vertical fluxes requires knowledge of the exact distribution and extreme values within each cloud, which is impractical. Comparison between the results of Equations 10 and 11 suggests a better performance of Equation 10, because the water fluxes are captured quite well. The underestimation of heat and buoyancy fluxes comes from the complicated distribution within a transition zone containing weak updraft motions, and therefore could be improved by considering detailed structures within the cloud. In particular, the cloud could be divided into different parts to take into account such inhomogeneity due to the radial dependence. Gu et al. (2020a) proposed a "core-cloak" conceptual model to parameterize the strong and weak updrafts separately, and showed that such a representation can improve the parameterization of vertical heat fluxes significantly.
Another important issue is that the vertical transport averaged across the domain is not necessarily well approximated, even if the vertical fluxes within the clouds can be well represented given sufficient understanding of internal cloud properties. For vertical water fluxes (Figure 7c,d), the contributions from clouds produce an underestimation near the cloud base but a slight overestimation at the cloud top, compared with the total flux across the full domain (grey dash-dotted line in Figure 7). For vertical heat fluxes (Figure 7a,d), the clouds produce less negative flux but slightly more positive buoyancy flux than the total fluxes near the cloud base. Most importantly, the clouds only contribute negative heat fluxes in the inversion layer, but the true total heat fluxes remain positive up to 1,800 m. The discrepancy in the inversion layer between the total fluxes and that from clouds comes from the contribution of downdrafts outside the cloud. The downdrafts lead to positive heat fluxes in the inversion layer because of the evaporative cooling associated with the mixing. However, the cloud-top downdrafts associated with a single cloud may be discontinuous around the cloud and also have asymmetry relative to the ambient vertical shear (Figures 2-4), so that they cannot be well represented with a distribution as a function of r∕R. Thus, an alternative representation of cloud-top downdrafts would be necessary for future convection parameterizations, as suggested by some previous studies from various perspectives (Zhao and Austin, 2005;Glenn and Krueger, 2014;Park et al., 2016;Davini et al., 2017;Brient et al., 2019;Gu et al., 2020a).

SUMMARY
In this study, normalized distributions for nonprecipitating shallow cumulus clouds at different vertical levels are constructed using a compositing algorithm applied to large-eddy simulation data of the BOMEX and ARM cases. Despite the variability between cloud slices, the distributions of different variables share some similar features, with the normalized magnitude having a maximum near the cloud centre and decreasing outward. The relative humidity is almost homogeneously distributed within the cloud due to saturation. It decreases slowly outward from the cloud edge, forming a moist buffering region to protect the cloud from interacting directly with the environment beyond.
The virtual potential temperature ( v ) and potential temperature ( ) mean distributions have more complicated behaviour in the vicinity of the cloud boundary. Below the inversion layer, moving from the centre of the cloud outwards, v decreases rapidly from positive to negative before the vertical velocity turns from updraft to downdraft near the cloud edge, indicating the presence of a transition zone. At the cloud base, it is shown that moisture plays an important role in cloud triggering. In the inversion layer, the normalized buoyancy and have very similar magnitudes and are both negative, showing that the evaporative cooling associated with cloud-top mixing is important there. The vertical component of turbulent kinetic energy (TKE) dominates the total TKE within the cloud at all vertical levels, but there is an increase in the contribution from the horizontal part within the inversion layer. The fact that normalized TKE and horizontal TKE decrease more slowly with increasing distance from the cloud centre also indicates the presence of a downdraft/shell and enhanced horizontal turbulent mixing.
We find that the mean vertical velocity is about half the maximum vertical velocity within the clouds at different levels. This finding is consistent with observations of deep convection (Wang et al., 2020) and also provides evidence for previous theoretical studies and parameterization schemes that took this linear relationship as an assumption without clear justification (Kuo and Raymond, 1980;Morrison, 2016). We use this relationship to discuss possible power-law formulations for the vertical velocity distributions within the cloud and find that a linear distribution is suitable for a 2D updraft and a quadratic distribution for 3D updrafts. This provides a reasonable explanation for the "triangular" shape distribution of vertical velocity captured by aircraft observations penetrating through the convection (Zipser and LeMone, 1980). It also provides some support for a recent convection parameterization that used a prescribed quadratic distribution of vertical velocity within the updraft to take into account the pressure perturbation effect on the vertical momentum budget (Leger et al., 2019).
A simple application of the normalized distributions of different variables is made, in an attempt to represent the vertical heat and moisture fluxes. The fluxes can be captured well if the maximum perturbations and normalized distributions for each cloud object are used, along with the assumption that all cloud objects have round shapes even if their geometric shapes are complicated due to turbulent mixing. Simplified calculations with averaged maximum perturbations, or averaged normalized distributions, or both, degrade both heat and moisture fluxes, suggesting that the variability of maximum perturbations between clouds and the variability in the transition region near the cloud edge need to be carefully treated in parameterizations that use assumed distributions. We also find that downdrafts outside the cloud make important contributions to the total heat and moisture fluxes in the inversion layer and therefore also need to be represented reasonably in future parameterizations. An alternative way to consider the distribution near the cloud edge would be to adopt a "core-cloak" representation, in which the cloud core region and the transition zone region are treated separately within a bulk mass flux approach (Gu et al., 2020a).
It should be recognized that the composited structures in this study are only representative for nonprecipitating shallow cumulus clouds. It is not clear to what extent our conclusions can be applied to precipitating shallow convection, stratocumulus clouds, deep convection, and organized convection. For example, organized clouds may require an extension of our current algorithm to consider the preference in organization direction. The distribution of dynamical and thermodynamic variables within the cloud can also be affected by microphysical processes (Endo et al., 2019) and subgrid unresolved processes. While we have found that the general features of composited structures do not change significantly with model grid spacings from 10-100 m, the detailed structures near the cloud edge might present richer information in even higher resolution simulations, because the interaction between cloud and environment cannot necessarily be resolved fully even at 10 m. Whether any such details would impact the estimation of vertical fluxes must be left for future studies.