Open Access

Abstract. Probability distribution functions of shallow cumulus cloud core entrainment and detrainment rates are calculated using 4362 individual cumulus clouds isolated from LES (large eddy simulation) using a cloud tracking algorithm. Calculation of the mutual information between fractional entrainment/detrainment and a variety of mean cloud core properties suggests that fractional entrainment rate is best predicted by the mean cloud buoyancy B and the environmental buoyancy lapse rate d θ ρ /d z at that level, while fractional detrainment is best predicted by the mean vertical velocity w and the critical mixing fraction χ c . Fractional entrainment and detrainment rates are relatively insensitive to cloud core horizontal area, and the perimeter of horizontal cloud core sections display an a 0.73 dependence. This implies that cloud core mass entrainment flux E is proportional to cloud core cross-sectional area instead of cloud core surface area, as is generally assumed. Empirical best-fit relations for e( B , d θ ρ /d z and δ( w , χ c ) are found for both individual shallow cumulus clouds and cloud ensembles. It is found that clouds with high buoyancy in strong stratification experience low entrainment rates, while clouds with high vertical velocities and critical mixing fractions experience low detrainment rates.


Introduction
Shallow cumulus clouds, sometimes referred to as tradewind cumulus, occur in the tropics as a transitional state between stratus decks, which occur in strongly stratified downwelling regions, and deep cumulus clouds, which occur in weakly stratified upwelling regions.Shallow cumulus reach heights of 2-3 km, transporting heat and moisture upward which erodes the inversion stratification and preconditions the atmosphere for deep convection.Biases in the parameterization of shallow cumulus in general circulation models (GCMs) have impacts on the distribution and intensity of deep convection, which can result in poor representations of the Hadley and Walker circulations (Stevens, 2005).
Additionally, shallow cumulus serve as a test for GCM cloud parameterizations, which in general have been developed for stratus or deep cumulus.Because of this, several shallow cumulus test cases, based upon field campaigns, have been created by the Global Energy and Water Cycle Experiment (GEWEX) Cloud System Studies (GCSS; Randall et al., 2003) boundary layer cloud group, suitable for modelling via large eddy simulation (LES; Stevens et al., 2001;Brown et al., 2002;Siebesma et al., 2003;vanZanten et al., 2011).Much of this work has focused on the entrainment and detrainment rates of shallow cumulus, which strongly affect shallow cumulus properties and constitute one of the largest sources of uncertainty in GCMs (Sanderson et al., 2008;Klocke et al., 2011).
Entrainment and detrainment of mass is defined as the rate at which mass crosses into (entrainment) or out of (detrainment) some region in a fluid, such as the region containing condensed liquid water (i.e. a cloud).The entrainment and detrainment rates of a cloud at a given height can be formally defined as (Siebesma, 1998) where E and D are the entrainment and detrainment rates (kg m −1 s −1 ), ρ is the density of air (kg m −3 ), u is the velocity of the air (m s −1 ), u i is the velocity of the cloud surface (m s −1 ), n is a unit vector directed out of the cloud surface, and the path integral is taken around the cloud surface at a constant vertical level.However, mass entrainment and detrainment are more often represented with the fractional mass entrainment and detrainment rates = E/M and δ = E/M (both m −1 ), where M = ρwa is the vertical mass flux (kg s −1 ), w is the vertical velocity (m s −1 ), and a is the cross-sectional area (m 2 ) of the cloud.These can be thought of as the fraction of the cloud mass that is being entrained and detrained per metre of rise through the cloud.
Many parameterizations of cumulus entrainment and detrainment rates have been proposed and tested against LES output (de Rooy et al., 2012).Turner (1963) proposed a simple scaling for entrainment as being proportional to the cloud vertical velocity times the perimeter of a cross section.This results in the fractional entrainment at a given height being inversely proportional to the cloud radius (assuming the cloud cross section is roughly circular).This has served as the basis of many parameterizations (Arakawa and Schubert, 1974;Tiedtke, 1989;Kain and Fritsch, 1990;Wagner and Graf, 2010), some of which make the further assumption that variations in the effective radius of the cloud field are negligible and so and δ can be treated as constants (Tiedtke, 1989;Bretherton and Park, 2008).Others parameterize the effective cloud radius as proportional to the height of cloud top (Bretherton et al., 2004), or simply allow to be inversely proportional to height (de Rooy and Siebesma, 2008).
Buoyancy sorting schemes allow entrainment and detrainment to depend on the properties of cloud and environment by assuming that cloud parcels experience a range of mixing rates, and the parcels which become negatively buoyant as a result of this mixing detrain from the cloud plume (Kain and Fritsch, 1990;Bretherton et al., 2004;de Rooy and Siebesma, 2008).The critical mixing fraction χ c -the fraction of environmental air in a mixture of cloudy and environmental air needed to make the mixture neutrally buoyant -is the primary control on entrainment and detrainment in these parameterizations, with larger χ c resulting in larger and smaller δ.In a similar spirit, Bechtold et al. (2008) and Stirling and Stratton (2012) allow entrainment to depend directly upon the atmospheric specific humidity.
Several parameterizations use various arguments to link entrainment and detrainment to the dynamic variables of the clouds.Neggers et al. (2002) proposed an inverse relation-ship between and vertical velocity w.Using arguments concerning the rate turbulent kinetic energy is produced in the cloud, Gregory (2001) proposed ∝ B/w 2 , where B is the buoyancy of the cloud (m s −2 ).von Salzen and McFarlane (2002) use ∝ dB/dz, while de Rooy and Siebesma (2010) present relations for and δ dependent on B/w 2 , w −1 dw/dz, and a −1 da/dz.Finally, Romps and Kuang (2010) proposed that entrainment is essentially random, and that entrainment rate should be parameterized as a stochastic process with a set probability of a discrete mixing event occurring for every L metres a parcel rises.
The wide range of parameterization forms present in the literature for the entrainment and detrainment rates suggests the modelling community has not yet reached agreement on which variables are the best predictors of these processes.It is therefore important to develop better ways to test these hypotheses over a wide range of cloud and environmental conditions.
Traditionally entrainment and detrainment rates are diagnosed in LESs using mean cloud field tracer budgets (Siebesma and Cuijpers, 1995).Recently, Romps (2010) and Dawe and Austin (2011a) have developed methods to calculate these rates directly from model velocity, humidity, and temperature fields.These directly calculated entrainment/detrainment rates are ≈ 3 times larger than those calculated via tracer budgets due to the presence of a shell of recirculated air surrounding the clouds which biases the tracer budget calculations (Dawe and Austin, 2011b).Unlike tracer budget calculations, these new direct calculation methods allow us to easily localize entrainment and detrainment to individual clouds, and provide us with a new way to study the dependence of entrainment and detrainment rates on cloud properties.
Since and δ of the cloud ensemble are the result of the entrainment and detrainment of the individual clouds in the ensemble, studying the entrainment and detrainment of the individual clouds should give some insight into the behaviour of the ensemble.Since a single LES simulates hundreds or thousands of clouds, analysis of individual clouds will produce several orders of magnitude more statistical samples of , δ, and other cloud properties from an LES than simply analysing the mean cloud field properties.To this end, this study uses the direct entrainment/detrainment rate calculation method detailed in Dawe and Austin (2011a) and the cloud tracking algorithm detailed in Dawe and Austin (2012) to estimate joint probability distribution functions of fractional entrainment and detrainment rates with a variety of cloud properties for individual LES shallow cumulus clouds.Using measures of the mutual information shared between cloud properties and the fractional entrainment and detrainment rates, we develop a parameterization to predict the mean fractional entrainment and detrainment rates of individual shallow cumulus clouds, and extend this to the prediction of the bulk entrainment and detrainment rates of the cloud ensemble.

Model description and output data sets
All LES calculations in this paper were made using the System for Atmospheric Modelling (SAM version 6.8.2; Khairoutdinov and Randall, 2003).Two model runs were performed, configured as standard GCSS cases: a Barbados Oceanographic and Meteorological Experiment (BOMEX; Siebesma et al., 2003) run, and an Atmospheric Radiation Measurement study (ARM; Brown et al., 2002) run.The BOMEX run was performed on a 6.4 km × 6.4 km horizontal × 3.2 km vertical domain for 6 h, and the first 3 h of simulation were discarded.The ARM run was performed on a 6.4 km × 6.4 km × 4.5 km domain and 8.5 h of output between hour 4.5 and 13 were saved.Both models were run with a 25 m grid size in all directions and a time step of 1 s.Precipitation was disabled in both runs.
Instantaneous model fields were output each minute, generating 180 snapshots for the BOMEX run and 510 snapshots for the ARM run.Individual cloud histories were then identified from the model outputs using the cloud tracking algorithm detailed in Dawe and Austin (2012).The algorithm identified 2838 individual clouds in the BOMEX run and 1524 clouds in the ARM run.
Note that few of the calculations performed in this paper rely on the time histories of individual clouds, and could have been performed equally well by identifying connected cloudy regions in the model snapshots.Using the cloud tracking algorithm allows us to connect detritus from a dissipating cloud to its parent cloud, reducing the effective number of small clouds identified in the simulation.Nevertheless, we do not expect our use of the cloud tracking algorithm to significantly alter our results relative to using clouds identified from snapshots of model output.
Cloud core properties of each cloud as a function of height were calculated at each saved time, where cloud core was defined as grid points with condensed liquid water, upward velocity, and positive buoyancy.Cloud core total specific moisture q t (units of kg H 2 O per kg moist air), specific condensed liquid water q l (kg H 2 O per kg moist air), liquid-water potential temperature θ l (K), density potential temperature θ ρ (K), and vertical velocity w profiles were calculated using conditionally sampled horizontal means.Cloud core horizontal area a was found by summing the horizontal area of cloud core grid cells at each height, and cloud core surface area S (m 2 ) was determined by summing the areas of cloud core grid cell faces adjacent to non-core grid cells at each height.Mean horizontal properties for the entire model slab were also recorded to generate cloud anomalies relative to the background mean and mean environmental stratification.
Cloud core buoyancy was calculated as where g (m s −2 ) is the acceleration due to gravity, and the bar denotes the horizontal mean over the entire model domain.
Additionally, for each cloud height we calculate the critical mixing fraction χ c via de Rooy and Siebesma (2008): where θ ρ = θ ρ − θ ρ , θ l = θ l,c − θ l,e and q t = q t,c − q t,e are the mean cloud-core properties minus the mean properties of the environment, c p (J kg −1 K −1 ) is the specific heat capacity of dry air at constant pressure, π = T /θ is the Exner function, the mean temperature T (K) of the cloud divided by the mean potential temperature θ (K), and α and β are constants with values of α ≈ 0.12 and β ≈ 0.4.Finally, the direct entrainment/detrainment estimation method of Dawe and Austin (2011a) was used to calculate vertical profiles of and δ.These calculations were done by horizontally summing the instantaneous mass entrainment E and detrainment D over a region including the cloud core plus all points immediately outside the cloud core.These extra points were included because the tetrahedral interpolation scheme used by Dawe and Austin (2011a) to track the motion of the cloud core surface occasionally locates the surface outside of the cloud core grid cells, which results in entrainment and detrainment occurring outside of the cloud core.This misplacement of the entrainment locations reduces the accuracy of the direct entrainment calculation (which itself is low-biased ≈ 20 % by the interpolation used to generate the cloud core surface) as some mass entrainment and detrainment is displaced vertically.However, the error introduced by this will be random and should not alter the dependence of the entrainment and detrainment rates upon the cloud core properties.The summed E and D values are then divided by the cloud core vertical mass flux M calculated using horizontal cloud core areas calculated by the tetrahedral surface interpolation algorithm to generate self-consistent and δ values (Fig. 1).
This results in 147 060 samples of cloud core properties at various heights and times for the BOMEX output, and 134 949 samples for the ARM output.Instantaneous cloud samples at a given height consisting of less than 16 grid cells (cross-sectional area 10 000 m 2 ) were then filtered from the sample set, as they were subject to large amounts of gridscale noise.This mainly has the effect of removing small clouds and the tops and bottoms of larger clouds.Excluding these small area cloud samples removes nearly half of the cloud samples at a given height (Fig. 2a); however, the total cloud fraction (Fig. 2b) and vertical mass flux (Fig. 2c) of the cloud field is only reduced by ≈ 5 %.After filtering, 65 303 samples remain for the BOMEX output and 87 327 samples remain for the ARM output.

Cloud core property PDFs
Here we examine probability density functions (PDFs) of cloud core properties in the BOMEX output.Since the BOMEX case forcing does not vary in time, we amalgamate all three hours of model output into a single data set.This results in over 1000 cloud property samples at each height.However, since the decorrelation timescale for individual cloud properties is ≈ 15 min, only ≈ 100 of these samples are actually independent.
Visual inspection of the distribution of cross-sectional area a of the cloud samples suggests it is consistent with the power law distribution found in studies of shallow cumulus clouds, as the distribution decreases monotonically with size, while total specific water q t , liquid-water potential temperature θ l and vertical velocity w appear normally distributed (Fig. 3).The mean values of the q t , θ l and w PDFs coincide with the overall horizontal mean values conditionally sampled on the cloud core.The variance of a is relatively constant with height, while the variances of q t , θ l and w steadily increase from cloud base to the start of the inversion at 1500 m.Once in the inversion, the variance of a, q t , and θ l rapidly decreases with height, while the variance of w remains high.
Next we examine some derived cloud core properties: buoyancy, critical mixing fraction χ c , and fractional entrainment and detrainment rates (Fig. 4).Buoyancy displays a small positive skewness, and combined with the requirement that B is positive in the cloud core this suggests B is best modelled with a log-normal distribution.Critical mixing fraction shows a normal distribution, while and δ show strong log-normal distributions.The mean values of the cloud core B and χ c PDFs again agree with the horizontal mean value of the conditionally sampled core, while the mean of the cloud core log 10 ( ) and log 10 (δ) distributions agree with the the log 10 of the net cloud core ensemble and δ.The variance of B and χ c increases through the cloud layer then decreases rapidly in the inversion, while the variances of log 10 ( ) and log 10 (δ) are essentially constant with height.

Mutual information analysis
In this section we analyse merged output from both the ARM and BOMEX cases to determine which cloud core properties are the strongest predictors of the cloud core mass entrainment and detrainment rates.This analysis is complicated by strong correlations between cloud properties (Dawe and Austin, 2012) and non-linear relationships between cloud core variables and entrainment/detrainment rates, which make it difficult to unambiguously link variability in and δ with a single cloud property.In order to overcome these problems, we quantify the strength of dependencies between entrainment/detrainment and cloud properties using the mutual information (MI) shared between them (Shannon and Weaver, 1949).
MI is defined as I (X; Y ) = P (x, y) ln P (x, y) P (x)P (y) dxdy, (5) where P (x), P (y), and P (x, y) are the marginal and joint probability density functions for the variables X and Y .Similar to the Pearson correlation coefficient, a high MI between two variables implies a strong functional relationship between those variables, but unlike correlation, MI measures non-linear as well as linear relationships.Additional details on the MI calculation are provided in Appendix A.
We estimate the joint PDFs between variables using histograms.PDF estimates generated via histogram are dependant on proper bin choice: too few bins results in a poorlyresolved PDF, while too many bins results in each bin containing too few samples for a reliable PDF estimate.To determine appropriate bin spacing we performed our calculations for a range of bin choices.We restricted the data range so that the majority of bins contained more than 10 samples and found that between 20 and 30 bins generated similar PDFs and MI estimates.All PDFs we present here were calculated with 20 equal-width bins spread across the data range, with the exception of the and δ PDFs which, due to their log-normal distribution, were log-transformed before histogramming.(Repeating our calculations on the un-transformed and δ values gave similar results.)Data limits and bin widths are summarized in Table 1.
We note here that MI provides a purely statistical analysis of the relationships between variables, without reference to the dynamics of the clouds.The relationships the MI analysis finds have no physical basis and may actually result from indirect correlations between the variables we examine and the true underlying dynamics of the system.The relationships we find may be best considered a kind of null hypothesis: a useful physically-based parameterization of entrainment and detrainment should outperform this statistical analysis.In light of this, we refrain from attempting to interpret our results in dynamical terms until the discussion in Sect. 5.

Entrainment
In this section we examine the dependence of the fractional mass entrainment rate on a variety of cloud variables.The literature provides several examples of entrainment parameterizations using a variety of variable combinations (Turner, 1963;Tiedtke, 1989;Kain and Fritsch, 1990;Neggers et al., 2002;de Rooy andSiebesma, 2008, 2010), but we have  chosen to focus on the basic cloud properties in our analysis for several reasons.First, if the parameterizations have predictive power the MI analysis should pick out the parameterization variables automatically.Second, we perform our calculations using directly measured mass entrainment rates, which differ from the modified rates used in entrainment parameterizations which must account for the influence of the moist cloud shell (Dawe and Austin, 2011b).Third, when we calculated the MI between log 10 ( ) and several parameterizations, they generally showed MI values smaller than the cloud variables we present here.
We estimate the joint PDFs between log 10 ( ) and the following cloud core properties: vertical velocity w, cloud core horizontal area a, buoyancy B, critical mixing fraction χ c ,  the lapse rate of environmental density potential temperature dθ ρ /dz (K m −1 ), and the height z (m).We consider the joint PDF of log 10 ( ) and height z as a null hypothesis, as there is little reason the absolute height above ground should, by itself, affect the entrainment rate.The resulting joint PDFs display remarkably similar behaviour for all variables, with larger variable values associated with smaller log 10 ( ) (Fig. 5).This is not surprising in light of the strong correlations present between shallow cumulus cloud properties (Dawe and Austin, 2012).
MI values for log 10 ( ) are given in Table 2. Buoyancy B shows the largest MI value with log 10 ( ), with a value nearly double the next largest, I (log 10 ( ); χ c ).All variables show MI values larger than the maximum value generated by calculating the MI between log 10 ( ) and 100 random permutations of each variable, which we use as a measurement of statistical significance.However, I (log 10 ( ); a) is smaller than I (log 10 ( ); z), suggesting a has little influence on the entrainment rate.
Cross-sectional area a shows the smallest MI value with log 10 ( ), and the relative lack of dependence of the mean value of log 10 ( ) on a is readily apparent in the PDF (Fig. 5c).This is surprising, as we would expect entrainment rate to be related to the surface area of the core surface, which in turn should be related to the area occupied by the clouds.However, the variance in log 10 ( ) is strongly dependent on a, with the largest and smallest values of log 10 ( ) only occurring for the smallest area clouds.We take this to indicate a strong patchiness and spatial localization in the distribution of entrainment.Small clouds may be subject to flow structures that drive large or small amounts of entrainment, but large clouds average over these flow structures, mitigating the variability in entrainment they experience.Nevertheless, even at the largest cloud sizes, there is still nearly an order of magnitude range in the variability of .
One possible cause of the relative constancy of log 10 ( ) versus a is the existence of correlations between a and other cloud core properties.For example, cloud core area is positively correlated with buoyancy (Dawe and Austin, 2012).If larger area clouds tended to have reduced log 10 ( ) this would be offset in the joint PDFs by the tendency for high buoyancies to increase log 10 ( ); the true dependence of log 10 ( ) on a would be masked by the covariance of a and B. We can separate out the effects of these correlations by generating joint PDFs of log 10 ( ) with two variables simultaneously.To do this we calculate three dimensional histograms to estimate P (log 10 ( ), y, ζ ), where y and ζ are various combinations of the cloud core properties z, a, B, w, and dθ ρ /dz.(We use ζ for the third PDF variable as we have already designated z to represent height.)The resulting histograms show little change when calculated using 20 and 30 bins along each dimension, so we maintain the same bin widths as used in generating the two dimensional histograms.
These three dimensional joint PDFs reveal a great deal of information about the behaviour of , but can be difficult to visualize in two dimensions.Visual inspection of the distributions of log 10 ( ) at various points in the (y, ζ ) space show reasonably Gaussian distributions.If is log-normally distributed at all points in the (y, ζ ) space, then half of the distribution should have larger values than the mean of log 10 ( ) and half should have smaller values.Thus, to visualize these PDFs we plot the mean of log 10 ( ) over the (y, ζ ) space to show how the distributions change with variables.
The easiest plot to interpret is probably the mean of the joint PDF of log 10 ( ), height, and area P (log 10 ( ), z, a) (Fig. 6, row 1, column 1), which clearly shows the vertical variation in with height and the slight decrease in as cloud area increases.Since we do not expect height to directly influence , apparent variations in with height actually arise due to changes in the mean cloud properties.This is apparent comparing P (log 10 ( ), z, a) to P (log 10 ( ), z, B) (Fig. 6, row 1, column 2), in which nearly all the variation in collapses onto changes in B. In fact, at nearly every height the mean of is better correlated with B, χ c , w and dθ ρ /dz than z (Fig. 6, row 1).Similarly, the mean of the joint PDFs of log 10 ( ) and a (Fig. 6, column 1) show the apparent variation of with a is better explained by correlations between a and other cloud properties.
The remaining plots are less clear-cut, with buoyancy, critical mixing fraction, and vertical velocity all displaying strong independent covariability with .Buoyancy shows the strongest covariance with (Fig. 6, column 2), in agreement with the MI calculations, but it is difficult to judge which variable is the second most important.To quantify which variable provides the most information about that is independent of B, we calculate the conditional mutual information (CMI) for each PDF: (We have designated samples of the random variable Z with ζ to avoid confusion with the height z.)By conditioning the PDF of Y on the value of Z, CMI removes the mutual information between X and Z, revealing the MI between X and Y .
We calculate CMI between log 10 ( ) and the cloud core properties conditioned on B to determine which variable provides the most information that is not already provided by B (Table 2).dθ ρ /dz shows the largest CMI with log 10 ( ) when conditioned on B, despite the small MI between log 10 ( ) and dθ ρ /dz.Note that the CMI between log 10 ( ) Calculating the remaining CMI of log 10 ( ) conditioned on both B and dθ ρ /dz shows values close to the noise level of the calculation (Table 2).This indicates that nearly all the information about recoverable from the cloud core state can be found using only B and dθ ρ /dz.Note that this may be an artifact of an insufficient number of samples to properly resolve the full multi-dimensional histograms.

Detrainment
In this section we repeat the previous analysis to examine the dependence of the fractional mass detrainment rate δ on the cloud core properties.The joint PDFs of log 10 (δ) with the cloud properties (Fig. 7) are a little more complex than the log 10 ( ) PDFs.As with the entrainment, larger w, a, B, and χ c values are associated with smaller log 10 (δ).Unlike the entrainment, log 10 (δ) appears to increase with stronger stratification when dθ ρ /dz ≥ 3 K km −1 .Detrainment shows slightly more dependence on a than entrainment, though this covariance is still small relative to the other variables.Finally, log 10 (δ) decreases with w between 0-3 m s −1 but increases between 3-6 m s −1 .Fig. 6.Mean values of log 10 ( ) for each bin of joint probability density functions of various cloud core properties for individual clouds in the combined BOMEX and ARM LES output.The y axis of each row shows height, environmental stratification, critical mixing fraction, buoyancy, and vertical velocity (from top to bottom), and the x axis of each column shows cross-sectional area, vertical velocity, buoyancy, critical mixing fraction, and environmental stratification (from left to right).
Values of MI between log 10 (δ) and the cloud properties are given in Table 3.The largest MI value results from I (log 10 (δ); χ c ), implying that the critical mixing fraction is the best predictor of detrainment rate.This result is in broad agreement with a range of previous work on parameterization of cloud core detrainment (Kain and Fritsch, 1990;de Rooy and Siebesma, 2008;Bretherton and Park, 2009).
Mean values of joint PDFs for log 10 (δ) are presented in Fig. 8.These clearly display the strong relationship between δ and χ c .The relatively strong dependence of log 10 (δ) on buoyancy disappears completely when P (log 10 (δ), χ c , B) (Fig. 8, row 4, column 2) is examined.Area a shows a moderate effect on log 10 (δ) that is independent of χ c , but the vertical velocity w shows the largest effect on log 10 (δ) independent of χ c (Fig. 8, row 3, column 3).The largest detrainment rates occur when both χ c and w are small.This is confirmed by calculating CMI values between log 10 (δ) and the cloud core properties conditioned on χ c (Table 3).The vertical velocity shows the largest CMI with log 10 (δ), over twice the CMI of a with log 10 (δ) conditioned on χ c .This strong inverse relationship between log 10 (δ) and w is reminiscent of the parameterization of Neggers et al. (2002).Neggers et al. proposed a w −1 behaviour for , not δ, but it is not implausible that turbulent entrainment and detrainment would follow the same behaviour.Furthermore, since w and B are correlated (Dawe and Austin, 2012), a dependence of on B would also cause a correlation between and w.As our analysis is purely statistical we are unable to unambiguously attribute the behaviour of and δ to dependence on B or w, but either way, these results support the constant timescale w −1 behaviour observed by Neggers et al.
Finally, we calculate CMI values of log 10 (δ) conditioned on both χ c and w (Table 3).The largest CMI value in this case results from z, and is only roughly twice the statistical noise level, so we conclude there is little meaningful information remaining.

Cloud perimeter vs. area
An interesting result of the previous analysis is the apparent independence of and cloud cross-sectional area a.Many entrainment parameterizations follow the assumption made by Turner (1963) that entrainment follows the scaling where k is a dimensionless constant, C is the perimeter of the cloud cross section (m), and R is the cloud radius (m) (Arakawa and Schubert, 1974;Kain and Fritsch, 1990).(The second form of the equation is derived by assuming the cloud is cylindrical, so C = 2πR and a = πR 2 .)This assumption appears to be incorrect; fractional entrainment rate is almost independent of cloud area, at least for shallow cumulus.This may help explain the efficacy of the assumption made by some parameterizations (Tiedtke, 1989;Bretherton and Park, 2009) that R is constant, which implies C ∝ a, E ∝ ρwa and thus = E/M is independent of area.
Real clouds, of course, are not cylindrical.If the perimeter of the fractal cloud surface were to scale linearly with a, this would explain the relatively constant value of with a. Siebesma and Jonker (2000) found a C ∝ a 0.66 relationship in a BOMEX LES where a and C were calculated from a two-dimensional projection of the cloud area to mimic satellite observations, significantly different from the C ∝ a 0.5 relationship one would expect for a cylindrical cloud.However, entrainment will depend on the area-perimeter relationship of horizontal cross sections through the cloud, which may differ from the two-dimensional projection used in their study.
In light of this, we calculate the perimeter-area relationship for horizontal cloud slices in our LES output.Since the LES is a discrete model, we calculate a pseudo-perimeter for each cloud by taking the cloud surface area at a given height and dividing it by the LES vertical grid spacing dz (25 m).We calculate a fit of the curve C = ka n to the data by performing a linear least-squares best fit between log a and log C to find log C = n log a + log(k), which results in n = 0.73 and k = 1.50 (C = 1.50a 0.73 , correlation 0.950 ± 0.001, root mean square (RMS) error 2158 m, Fig. 9).This relationship shows a significantly larger correlation than either a linear (C = 0.043a, correlation 0.934 ± 0.001, RMS error 3008 m) or a square root (C = 26.0√ a, correlation 0.945 ± 0.001, RMS error 2695 m) relationship between C and a, where we have constrained these fits so that C(0) = 0.If C ∝ a 0.73 , then Eq. ( 7) implies ∝ a −0.27 .However, our analysis also shows that is independent of a when other variables are held constant.This contradiction implies that the basic concept underlying Eq. ( 7) -that mass entrainment flux E is proportional to the cloud surface area -is not true for these simulated clouds.
We can check this using the LES output by fitting power law relationships between the entrainment and detrainment fluxes and the cloud core area and pseudo-perimeter.Doing so shows E ∝ a 0.97±0.01, E ∝ C 1.29±0.01, D ∝ a 0.87±0.01, and D ∝ C 1.08±0.01(Fig. 9).Thus, it appears that E is indeed proportional to cross-sectional area, while D is not obviously proportional to either a or C.
Why this surprising result should be the case is not readily apparent.We find approximately the same results when we filter cloud heights within 100 m of cloud base, where E might reasonably be expected to be proportional to crosssectional area due to the condensation-produced buoyancy of the rising thermals.In any case, the linear dependence between E and a is clearly fortuitous for the purposes of cloud parameterizations.

Parameterization of entrainment and detrainment rates
While mutual information provides us with a way to measure the dependencies between variables in a data set, it says nothing about the functional form of those dependencies.In this section we attempt to construct a parameterization for the and δ of individual shallow cumulus clouds by curve fitting simple power law relationships for (B, dθ ρ /dz) and δ(w, χ c ).Additionally, we attempt to extend these fits to parameterize and δ values for the overall cloud ensemble.We wish to emphasize that these parameterizations are purely statistical in nature, with little reference to the underlying dynamics of the system.For example, they do not produce relationships with units of m −1 and require constant multipliers with units that correct for dimensional consistency, unlike most published parameterizations.We do not advocate that statistical fits be used as parameterizations without an understanding of the dynamics of the system, but instead we suggest they be used as a null hypothesis for the behaviour of and δ.In other, words, parameterizations of shallow cumulus mass entrainment and detrainment should at minimum display higher correlation and lower RMS error when compared with statistical power law fits to be considered valid.However, since most currently published parameterizations predict tracer, not mass, entrainment rates, we are not able to directly compare our results.

Entrainment
In this section we examine the dependence of fractional entrainment on buoyancy B and stratification dθ ρ /dz.We fit power law relationships between the variables by performing linear least-squares best fits between log 10 ( ), log 10 (B), and log 10 (dθ ρ /dz) to find relationships of the form = 10 b x m , where m and b are the slope and intercept of the line fit log 10 ( ) = m log 10 (x) + b. (We perform these fits using logarithms in base 10 instead of natural logs simply for easier interpretation of the resulting plots.)The joint PDF of log 10 ( ) and log 10 (B) shows a relationship with a great deal of variance (Fig. 11a).Nevertheless, the correlation between the variables of −0.78 is fairly high, and a linear fit with a slope of −1.29 does a reasonable job matching the data (RMS error of 0.46) except at small and large values of B where is systematically underestimated.
The joint PDF of log 10 ( ) and log 10 (dθ ρ /dz) shows a less robust linear relationship (correlation of −0.48, Fig. 11b).There appear to be two regimes, one between log 10 (dθ ρ /dz) values of −2.75 and −2, and a second between −3 and −2.75.The linear fit to this data (slope −1.35, RMS error of 0.64) does reasonably well for larger stratification, but at low stratification is significantly underestimated.
Since the relationships between and B and between and dθ ρ /dz have the same sense -stronger buoyancy and stratification mean weaker entrainment -we also try a fit to the product of the two variables.The joint PDF of log 10 ( ) and log 10 (Bdθ ρ /dz) shows a stronger linear relationship than either variable individually (correlation of −0.83, Fig. 11c).The resulting curve fit (slope −1.06, RMS error 0.35) still underestimates the entrainment rate at low and high values of buoyancy and stratification, but many of the extremely low entrainment values present at low buoyancy are raised by the addition of the stratification.Additionally, the fit is tantalizingly close to indicating a simple inverse relationship between and Bdθ ρ /dz.We therefore conclude that a power law fit between and Bdθ ρ /dz provides a simple but skillful estimate of the entrainment rate of individual shallow cumulus clouds.

Detrainment
Again we repeat the previous analysis to examine the dependence of fractional detrainment δ on vertical velocity w and critical mixing fraction χ c .The joint PDF of log 10 (δ) and log 10 (w) shows a relatively poor linear relationship (corre-−0.52,Fig.However, the fit to this data shows a good inverse relationship between the variables −0.92, RMS error 0.34).However, unlike the other the sense of the relationship between the variables has two regimes: δ decreases with increasing w between log 10 (w) values of −0.5 and 0.5, but increases with w above a log 10 (w) of 0.5.
The joint PDF of log 10 (δ) and log 10 (χ c ) shows a much more linear relationship than log 10 (w) (correlation −0.71, Fig. 12b), but the curve fit between these variables is weaker (slope −1.34, RMS error 0.62) and overestimates δ at small values of χ c .
However, as with , a fit between δ and the product wχ c does a better job than either alone.The joint PDF displays a stronger linear relationship than either variable alone (correlation −0.76, Fig. 12c), and the curve fit matches the mean PDF values well (slope −0.86, RMS error 0.32), also displaying a nearly inverse relationship between δ and wχ c .

Resulting best fits
Figure 13 shows the resulting best-fit power-law curves for and δ.We find the relationships = 5.19 × 10 −8 (Bdθ ρ /dz) −1.07 (8) and δ = 2.72 × 10 −3 (wχ c ) −0.89 (9) provide reasonable fits to the individual cloud data.(The power law exponents for these fits are slightly different than those reported in the previous section, as we have expanded the data range over which we are performing these fits to include some of the more extreme model values.)As noted above, while suffering from a lack of dimensional consistency with and δ, these parameterizations do have the intriguing property of nearly being a simple inverse relationship.In any case, these relationships provide a reasonable first-order estimate of the magnitude of and δ for individual cumulus clouds in the BOMEX and ARM cases.However, large-scale models require parameterizations not of individual clouds, but of cloud ensembles.Translating these individual curve fits into equations usable for whole cloud fields is a problematic task.Formally, the ensemble entrainment and detrainment rates at a given height can be written in terms of (B, dθ ρ /dz) and δ(w, χ c ) as Instead of performing this complicated transformation, we simply refit (B, dθ ρ /dz) and δ(w, χ c ) using the ensemble values of , δ, and the mean cloud core properties.Figure 14  shows the resulting best-fit power-law curves for : and for δ: Surprisingly, both fits display a ≈ −0.7 power law, a coincidence for which we have no explanation.

Discussion
As far as we are aware, this study represents the first time PDFs of cloud core entrainment and detrainment with various other cloud core properties have been calculated for individual clouds in LES.This allows us to easily examine the behaviour of shallow cumulus ensembles and find novel results, such as the negligible role in vertical transport played by smaller cumulus, or the C ∝ a 0.73 relationship between perimeter and area for clouds.Applying this technique to problems such as examining the cloud PDF changes that occur during the transition from shallow to deep convection would no doubt provide equally novel results.
Our analysis implies that δ is roughly inversely proportional to w, which is highly reminiscent of the multiparcel entrainment model of Neggers et al. (2002).While w is not the strongest predictor of in our results, larger w is undoubtedly associated with smaller (Fig. 5b), adding support to the Neggers et al. model. Romps and Kuang (2010) criticized the Neggers et al. model on the basis that cloud base properties are very uniform, so a dependence of entrainment rate on cloud properties would not produce the wide variance observed in cloud properties, proposing that stochastic entrainment events instead represent multiparcel entrainment better.Our results suggest that both models are partially correct, as and δ show strong dependence on cloud properties but also display a great deal of randomness (Fig. 13).For example, calculating the standard deviation of log 10 for moderate values of B and dθ ρ /dz gives log 10 = (−2.0± 0.3), equivalent to an range of ≈ (0.005-0.2) m −1 .
This creates the following picture of shallow cumulus dynamics: all clouds start with uniform properties at cloud base.Randomness in the mixing events experienced by each cloud produces property differences as the clouds rise, and these property differences then feed back upon the entrainment and detrainment rates.Entrainment events tend to decrease B, w and χ c , and reductions in B, w and χ c all imply higher rates of entrainment and detrainment.A positive feedback thus appears to exist, so that clouds which experience large entrainment events early in their evolution tend to be more vulnerable to further dilution.
Our results have all been calculated for direct and δ values, which differ substantially from the and δ values needed for GCM parameterizations due to the influence of the moist cloud shell, Dawe and Austin (2011b).This makes it difficult to compare our results directly with previously published parameterizations, as we cannot be sure that differences between the parameterizations and our results are due to real differences or due to neglect of the modifying effects of the cloud shell.We hope to address deficiency in a future paper.
However, some preliminary indications of the validity of cloud parameterizations can be drawn from our results.We find little dependence between and cloud area, contradicting parameterizations that vary entrainment rate with cloud radius (Arakawa and Schubert, 1974;Kain and Fritsch, 1990).This disagrees with the results of Stirling and Stratton (2012), who find a clear relationship with both and δ decreasing with mean area per cloud in simulations of the diurnal cycle of deep convection.However, Stirling and Stratton examine much larger clouds than our shallow cumulus, are looking at the behavior of the bulk cloud field rather than individual clouds, and measure entrainment with tracer budgets rather than direct calculations, which may help explain the discrepancy.
As assumed by buoyancy sorting parameterizations (Tiedtke, 1989;Bretherton and Park, 2009), critical mixing fraction χ c has strong effects on δ, with large χ c suppressing detrainment.However, large χ c is also associated with reduced , while buoyancy sorting parameterizations predict enhanced .Other parameterizations relate to the quantity or some simplification thereof (de Rooy and Siebesma, 2010;Gregory, 2001).The inverse dependence of and B do not rule out these types of relations, as w and B are correlated, and it is possible that a B/w 2 relationship could appear as a 1/B due to these correlations.However, none of these relationships take into account the effect of background stratification, which appears to be important in our results and has also been incorporated into an entrainment parameterization by Stirling and Stratton (2012).Finally, as mentioned previously, we find support for the w −1 parameterization of Neggers et al. (2002).Finally, we acknowledge the awkwardness inherent in our parameterization relationships not being dimensionally consistent with the m −1 of and δ.However, it is not clear that a complex turbulent phenomena like these should necessarily display dimensional consistency; the C ∝ a 0.73 relationship we find for the dependence of cross-sectional area and perimeter certainly does not, for example.Either way, we believe the main value of these results is as a guide to developing better parameterizations, acting as a statistical null hypothesis that a mass entrainment and detrainment scheme should better in terms of larger correlations and smaller RMS error.

Conclusions
Joint probability density functions for fractional cloud core mass entrainment/detrainment rates and horizontal mean cloud core properties were calculated for individual clouds isolated from BOMEX and ARM LES with a cloud tracking algorithm.Clouds with cross-sectional cloud core area a less than 10 000 m 2 were found to have negligible effects on the vertical mass and property transports of the cloud field despite occupying ≈ 50 % of the total cloud core area, and were excluded from the analysis.
PDFs of cloud core properties showed most properties having normal distributions, with the exception of cloud core area, which was exponentially distributed, and fractional entrainment and detrainment rates, which displayed log-normal distributions.Joint PDFs between and δ with a showed little dependency of fractional entrainment or detrainment rates on cloud core area.Examination of the relationship of cloud core cross-sectional perimeter and area showed C ∝ a 0.73 .
Dependence between , δ and various cloud core properties was quantified using mutual information.was found to have the highest MI with B and dθ ρ /dz, and δ with χ c and w.Overall, and δ appear to be primarily governed by buoyancy, either directly or through the critical mixing fraction.Highly buoyant clouds experience less fractional entrainment and detrainment than less buoyant clouds.Similarly, highly stratified environments reduce cloud core entrainment and large upward velocities reduce cloud core detrainment.Power law fits of the form (B, dθ ρ /dz) and δ(w, χ c ) were found to provide reasonable predictions of entrainment/detrainment rates both for individual clouds and the overall cloud ensemble.However, this study has only examined two shallow cumulus regimes; a bigger parameter space is needed to validate these results.
We have not directly compared the results presented here with entrainment and detrainment rate parameterizations in the literature, primarily because these parameterizations have been tuned to values derived from tracer budget calculations, which are more applicable to ensemble cloud rate calculations needed for GCMs.Performing this calculation requires correcting the direct entrainment/detrainment rates for the effect the moist cloud shell has on tracer fluxes.Evaluating these equivalent tracer budget rates and comparing them to cloud parameterizations currently in use would help understand how applicable these results are to GCM parameterizations.

Fig. 2 .
Fig. 2. Vertical profiles of (a) number of cloud core samples, (b) cloud core cross-sectional area, and (c) cloud core vertical mass flux summed over all cloud samples in the entire BOMEX LES run (black line) and all cloud samples with instantaneous cross-sectional area larger than 10 000 m 2 (red line).

Fig. 3 .
Fig. 3. Probabilities in range x (P (x) x) at each model height (top row) and at 1 km height (bottom row) for cloud core (a) cross-sectional area, (b) total specific humidity, (c) liquid-water potential temperature, and (d) vertical velocity in the BOMEX LES output.White lines indicate the horizontal mean of each variable conditionally sampled on the cloud core over the entire model domain.

Fig. 4 .
Fig. 4. Probabilities in range x (P (x) x) at each model height (top row) and at 1 km height (bottom row) for cloud core (a) buoyancy, (b) critical mixing fraction, (c) log 10 of the fractional mass entrainment, and (d) log 10 of the fractional mass detrainment in the BOMEX LES output.White lines indicate the horizontal mean of each variable conditionally sampled on the cloud core over the entire model domain.

Fig. 5 .
Fig. 5. Joint probability density functions multiplied by bin area (P (x, y) x y) for individual clouds in the combined BOMEX and ARM LES output of log 10 ( ) versus (a) height, (b) vertical velocity, (c) cross-sectional area, (d) buoyancy, (e) critical mixing fraction, and (f) vertical gradient of environmental density potential temperature.PDFs are plotted using a logarithmic scale.White lines indicate the mean log 10 ( ) value as a function of the x axis variable.

Fig. 7 .
Fig. 7. Joint probability density functions multiplied by bin area (P (x, y) x y) for individual clouds in the combined BOMEX and ARM LES output of log 10 (δ) versus (a) height, (b) vertical velocity, (c) cross-sectional area, (d) buoyancy, (e) critical mixing fraction, and (f) vertical gradient of environmental density potential temperature.PDFs are plotted using a logarithmic scale.White lines indicate the mean log 10 (δ) value as a function of the x axis variable.

Fig. 8 .
Fig.8.Mean values of log 10 (δ) for each bin of joint probability density functions for various cloud core properties of individual clouds in the combined BOMEX and ARM LES output.The y axis of each row shows height, environmental stratification, critical mixing fraction, buoyancy, and vertical velocity (from top to bottom), and the x axis of each column shows cross-sectional area, vertical velocity, buoyancy, critical mixing fraction, and environmental stratification (from left to right).

Fig. 9 .
Fig. 9. Joint probability density function multiplied by bin area of cloud core perimeter C versus cloud core cross-sectional area a (P (a, C) a C) for individual clouds in the combined BOMEX and ARM LES output.PDF is plotted using a logarithmic scale.White line shows the mean perimeter of clouds as a function of cross-sectional area.Black dotted, dashed, and solid lines show best-fit lines for linear, square root, and an arbitrary power law relationships, respectively.

Fig. 10 .
Fig. 10.Joint probability density function multiplied by bin area (P (x, y) x y) for individual clouds in the combined BOMEX and ARM LES output of log 10 (E) versus (a) log 10 of crosssectional area, (b) log 10 of cloud core perimeter, and log 10 (D) versus (c) log 10 of cross-sectional area, (d) log 10 of cloud core perimeter.PDFs are plotted using a logarithmic scale.White lines indicate the mean log 10 (E)/ log 10 (D) value as a function of the x axis variable, and black lines show linear least-square best fits of log 10 (E)/ log 10 (D) versus the x axis variable.

Fig. 11 .
Fig. 11.Joint probability density functions multiplied by bin area (P (x, y) x y) for individual clouds in the combined BOMEX and ARM LES output of log 10 ( ) versus (a) log 10 of buoyancy, (b) log 10 of environmental stratification, and (c) log 10 of buoyancy times environmental stratification.PDFs are plotted using a logarithmic scale.White lines indicate the mean log 10 ( ) value as a function of the x axis variable, and black lines show linear least-square best fits of log 10 ( ) versus the x axis variable.

Fig. 12 .
Fig. 12. Joint probability density functions multiplied by bin area (P (x, y) x y) for individual clouds in the combined BOMEX and ARM LES output of log 10 (δ) versus (a) log 10 of vertical velocity, (b) log 10 of critical mixing fraction, and (c) log 10 of vertical velocity times critical mixing fraction.PDFs are plotted using a logarithmic scale.White lines indicate the mean log 10 (δ) value as a function of the x axis variable, and black lines linear least-square best fits of log 10 (δ) versus the x axis variable.

Fig. 13 .
Fig. 13.Joint probability density functions multiplied by bin area (P (x, y) x y) for individual clouds in the combined BOMEX and ARM LES output between fractional mass entrainment/detrainment rates and best-fit entrainment/detrainment rates as predicted by cloud core properties.(a) log 10 ( ) versus the best-fit (B, dθ ρ /dz) relationship.(b) log 10 (δ) versus the best-fit δ(w, χ c ) relationship.PDFs are plotted using a logarithmic scale.White lines indicate the mean log 10 ( ) or log 10 (δ) values as functions of the best-fit relationship, and black lines show the log 10 ( ) or log 10 (δ) values predicted by the best-fit relationship.

Fig. 14 .
Fig. 14.Joint probability density functions multiplied by bin area (P (x, y) x y) for the horizontal mean cloud ensembles of the BOMEX and ARM LES output between fractional mass entrainment/detrainment rates and best-fit entrainment/detrainment rates as predicted by cloud core properties.(a) log 10 ( ) versus the best-fit (B, dθ ρ /dz) relationship.(b) log 10 (δ) versus the best-fit δ(w, χ c ) relationship.PDFs are plotted using a logarithmic scale.White lines indicate the mean log 10 ( ) or log 10 (δ) values as functions of the best-fit relationship, and black lines show the log 10 ( ) or log 10 (δ) values predicted by the best-fit relationship.

Table 1 .
Data limits and bin widths used to calculate histograms.

Table 2 .
Mutual information between log 10 ( ) and various cloud core properties for individual clouds in the combined BOMEX and ARM LES output.Noise level is found by taking the maximum of 100 Monte Carlo trials of mutual information between log 10 ( ) and a random permutation of each variable.Maximum mutual information values in each comparison category are in bold.

Table 3 .
Mutual information between log 10 (δ) and various cloud core properties for individual clouds in the combined BOMEX and ARM LES output.Noise level is found by taking the maximum of 100 Monte Carlo trials of mutual information between log 10 (δ) and a random permutation of each variable.Maximum mutual information values in each comparison category are in bold.
ρ /dz conditioned on B is higher than the MI between log 10 ( ) and dθ ρ /dz, indicating that correlations between B and dθ ρ /dz were obscuring the true strength of the dependence of on dθ ρ /dz.(Clouds in strong stratification tend to have low buoyancy, which increases .)Examination of P (log 10 ( ), dθ ρ /dz, B) (Fig.6, row 2, column 2) shows the largest values are present at low buoyancy and stratification.