Sub-cloud turbulence explains cloud-base updrafts for shallow cumulus ensembles: First observational evidence

Sub-cloud turbulent kinetic energy has been used to parameterize the cloud-base updraft velocity ( w b ) in cumulus parameter-izations. The validity of this idea has never been proved in observations. Instead, it was challenged by recent Doppler lidar observations showing a poor correlation between the two. We argue that the low correlation is likely caused by the diﬃculty of a ﬁxed-point lidar to measure ensemble properties of cumulus ﬁelds. Taking advantage of the stationarity and ergodicity of early-afternoon convection, we developed a lidar sampling methodology to measure w b of a shallow cumulus (ShCu) ensemble (not a single ShCu). By analyzing 128 ShCu ensembles over the Southern Great Plains, we show that the ensemble properties of sub-cloud turbulence explain nearly half of the variability in ensemble-mean w b , demonstrating the ability of sub-cloud turbulence to dictate w b . The derived empirical formulas will be useful for developing cumulus parameterizations and satellite inference of w b .


Introduction
Cloud-base updraft velocity (wb) is a crucially important variable as it influences various aspects of cumulus clouds (Rogers and Yau, 1996).The wb modulates the aerosol cloud-mediated effect by governing the supersaturation near cloud bases (Twomey, 1959;Rosenfeld, 2014).In polluted conditions, cloud droplet size and number concentration are more sensitive to wb than aerosol concentration and size (Reutter et al., 2009).Moreover, wb dictates lateral entrainment of cumulus that remains an unresolved bottleneck for climate modeling (Donner et al., 2016).
Despite its importance, current cumulus parameterization schemes rarely express wb explicitly (Donner et al., 2016).Most schemes parameterize the cloud-base mass flux (Mb) without specifying the wb.For example, Arakawa and Schubert (1974) determine the Mb by adjusting the cloud work function towards a value maintaining an equilibrium between the large-scale forcing and the convection.Krishnamurti et al. (1983) determine Mb under the assumption that convection must balance the column integrated vertical advection of moisture.Kain and Fritsch (1993) and Grell (1993) parameterize Mb by requesting the convection to remove the large-scale instability over the convective time scale.
The earliest effort that explicitly represents the wb in Mb closure is Brown (1979) who approximates the wb using the environmental vertical velocity from the surrounding nine points at lower tropospheric levels.This scheme is physically flawed by the fact that the air masses that initiate cumulus clouds are convective in nature.This issue is addressed by Neggers et al. (2009) and Fletcher and Bretherton (2010) (FB10) who argued that the wb could be dictated by the subcloud turbulent intensity.FB10 used a set of cloud-resolving simulations to empirically derive the following formula to represent the wb: wb = 0.28×TKEML 1/2 + 0.64, (1) in which the TKEML is the turbulent kinetic energy averaged horizontally and vertically in the subcloud mixed layer.FB10 shows that such a boundary-layer-based mass flux closure scheme outperforms several commonly used schemes for three cumulus cases.
Still lacking is observational evidence of the ability of TKEML to explain the wb.As quoted by Donner et al. (2016): "… parameterizations that do provide vertical velocities have been subject to limited evaluation against what have until recently been scant observations."The only observational pursuit to evaluate the Eq. ( 1) is from Lareau et al. (2018) who analyzed Doppler lidar observations of ~1500 individual shallow cumulus (ShCu) over the Southern Great Plains (SGP), finding that sub-cloud vertical velocity variance (a proxy for TKEML) explains only a few percent of the wb variability.This led them to cast doubt upon the relationship.They argue that sub-cloud updrafts must work against negative buoyancy near the top of the mixed layer to generate wb, and such a penetrative nature of the convection deteriorates their correlations.
Given the contrasting results, it is imperative to answer the question of whether or not subcloud turbulence explains the wb.This is not only important for cumulus parameterizations but also crucial for advancing other pursuits in the field of cumulus dynamics.First, theoretical inquiries of cumulus dynamics often rely on the assumption of a tight coupling between the subcloud turbulence and wb.For example, in one-dimensional bulk models of boundary layer clouds, a key variable is the Deardoff velocity scale, w * , which dictates the sub-cloud turbulence intensity (Betts, 1973;Neggers et al., 2006;Stevens, 2006;Zheng, 2019).Linking the w* with the wb is the basis for several important coupling processes between the cloud and sub-cloud layers (Neggers et al., 2006;van Stratum et al., 2014;Zheng et al., 2020).Second, recently emerging new satellite remote sensing methodologies of retrieving wb (Zheng and Rosenfeld, 2015;Zheng et al., 2015Zheng et al., , 2016) ) have offered great insights into the aerosol indirect effect and climate change (Rosenfeld et al., 2016;Seinfeld et al., 2016;Li et al., 2017;Grosvenor et al., 2018;Rosenfeld et al., 2019).
These studies infer the wb via quantifying the TKEML or its equivalents.Evaluating if the TKEML explains the wb is essential to evaluate the physical validity of these techniques.
To that end, this study examines the relationship between the wb and sub-cloud turbulence for ShCu using DL observations over the SGP.We focus on wb of ShCu ensembles, not single ShCu, because the former is more relevant to cumulus parameterization.We show that ensembleaveraged wb and sub-cloud turbulence are highly correlated with statistical significance (correlation coefficient greater than 0.7).Evaluating the relationship on ensembles but not on individual ShCu might explain the disparities with the previous finding (Lareau et al., 2018).The next session discusses the difference between the ensemble-mean wb and the wb of single cumuli.
It lays the foundation for developing the sampling strategy of ShCu ensembles.Section 3 introduces the observational data and methodology.Section 4 shows the results, followed by a summary.

wb of cumulus ensembles
Distinguishing between the ensemble and individual ShCu is necessary.The concept of cumulus ensemble is a fundamental building block for all cumulus parameterizations (Arakawa and Schubert, 1974).A cumulus ensemble on spatial scales of several tens of kilometers is composed of individual cumulus with a wide range of distributions in size and age.Since the individual cumulus clouds are at different stages of their lifetime, their physical properties differ considerably even if the surface and large-scale forcing are uniform.
The difference could be illustrated by Figure 1  distinctive pattern with strong updrafts within clouds surrounding by shells of downdrafts (Fig. 1a).We can see a rough correspondence between the vertical velocity field at the cloud-base level (Fig. 1a) and the TKEML (Fig. 1b): regions with larger TKEML typically have stronger updrafts near cloud bases.Such a correspondence, however, breaks down on the length scale of a single ShCu.For example, the vertical velocity field shows strong updrafts within individual clouds surrounding by shells of downdrafts whereas the TKEML variability across the cloud edges is considerably more uniform.This is not surprising since both updrafts and downdrafts contribute to the vertical mixing, jointly regulating the TKEML.As a result, their covariation on the length scale of individual ShCu tends to be noisy, which is confirmed by Figure 1c that compares the two quantities averaged over individual ShCu.The degree of scattering is likely to increase substantially when the synoptic and surface forcings are allowed to change.uncertainties.Fortunately, a convective boundary layer often experiences a quasi-steady state (Moeng, 1984;Lensky and Rosenfeld, 2006;Stull, 2012).In atmospheric science, whether a dynamical system can be considered quasi-steady depends on the difference between the characteristic time scale of the system and the time scale of external forcing.For a typical convective boundary layer over the continent, the surface forcing time scale is on the order of a few hours (defined as half of the period when the surface heat fluxes remain positive) whereas the time scale for shallow convective circulations is several tens of minutes (i.e. the convective time scale) (Fig. S1a).Such a time scale separation allows the mixed layer to remain in a quasi-steady state in which changes in turbulent properties are negligible compared with the turbulence production and dissipation terms (Stull, 2012).This quasi-steady assumption is particularly valid in the early afternoon when the surface fluxes reach their plateau and their time derivatives minimize (Fig. S1b).As such, focusing on early-afternoon ShCu can reduce the uncertainty of sampling due to temporal evolution.
In summary, to measure the wb of ShCu ensembles from surface-mounted DL, the sampling window must be at least a few hours to sample enough amount of individual ShCu.Moreover, an ideal sampling period is the early afternoon when the boundary layer is close to stationarity.

Data and Methodology
We use observations from the Department of Energy's Atmospheric Radiation Measurements (ARM) SGP observatory.The key instrument used in this study is the DL.The DL measures vertical velocity with ~ 1 s temporal and 30 m vertical resolution.The transmitted wavelength is 1.5 µm.In addition to DL, we also use data from radiosondes, a ceilometer, a Kaband cloud radar (KAZR), and ARM instruments measuring surface meteorological variables routinely.

An example case
To illustrate the sampling principle of ShCu ensembles, Figure 2a shows a MODIS satellite imagery of a ShCu field over the SGP at 20:30 UTC on June 10, 2012.The wind is southeasterly at a speed of ~ 9 m/s, corresponding to a horizontal distance of ~ 70 km over the two hours (the red solid line in Fig. 2a).One can see a few dozens of single cumuli drifting over the SGP site along the wind direction.Figure 2b shows a time-height plot of the DL from 19 to 21 UTC, corresponding to 13 ~ 15 local standard time (LST).Black dots mark the cloud-base heights (zb) measured by the ceilometer.To count how many individual cumuli are sampled during this period, we use the DL reflectivity to identify single cumuli.Figure 2c shows the zoomed-in window near cloud bases during the 19:48 ~ 20:00 UTC.The navy contours encompass pixels with DL reflectivity greater than 10 -4.6 m -1 sr -1 , a threshold that defines cloudy pixels (Lareau et al., 2018).
Based on the reflectivity threshold, a total of 84 individual clouds are identified during the 2-h period.The majority of them have a duration shorter than 4 s, which seems too short to constitute a single cloud.Thus, we congregate clouds with gaps < 20 s, reducing the cloud population to 29, with 12 of them lasting longer than 30 s.

Computing the wb
We select "cloud-base" DL pixels through two steps.First, to exclude the decoupled cloud elements and elevated cloud sides, pixels with cloud bases higher than 30% of lifting condensation level (LCL) are removed.Second, for the remaining coupled clouds, we select pixels within three gates below the cloud base (~ 100 m) and cloudy pixels above the cloud base.These pixels are defined as "cloud-base" pixels.Because of the strong signal attenuation, the DL only penetrates < 100 m into the clouds.Therefore, the cloudy pixels are mostly concentrated near several tens of meters above the cloud base.Figure S2 shows a comparison of the vertical velocity probability density function (PDF) between the two sub-groups of "cloud-base" pixels.Their PDF distributions are overall similar, suggesting that it is tenable to combine them as "cloud-base" pixels.
To compute the ensemble-mean wb, we average the selected vertical velocities in two ways.
The first is to simply average the vertical velocities above a threshold:  ̅ = ∑     ∑   ⁄ , in which the Ni represents the frequency of occurrence of positive vertical velocity wi that is greater than a critical value (wcrit).This is the common way for cloud-base mass fluxes study.The second way of averaging is weighted by volume:  ̅ vol = ∑     2 ∑   ⁄   .The volume-averaged updraft speed has been considered as more relevant to the understanding of aerosol cloud-mediated effects because it gives more weight to the larger vertical velocities that generate clouds with greater volume (Rosenfeld et al., 2014;Zheng et al., 2015;Rosenfeld et al., 2016).

Other quantities
Ideally, the TKEML should be computed as 0.5( ′ 2 +  ′ 2 +  ′ 2 ) averaged below the cloud base.However, the DL can only measure the vertical component, 0.5′ 2 , denoted as TKE w ML.In this study, we use the TKE w ML to approximate the TKEML, motivated by the fact that TKE w ML dominates the TKEML in typical convective boundary layers (Stull, 2012).The potential contributions from horizontal components of TKEML will be taken into account in our analyses in section 3.
We used the surface-measured temperature and moisture to compute the LCL using the exact analytical formula of Romps (2014).As described in the example case, we used the threshold of DL reflectivity to identify single cumuli.To compute the chord length of individual cumuli, we used the DL product of horizontal wind speed near cloud-base, which is derived from a velocity azimuth display algorithm (Teschke and Lehmann, 2017).The multiplication of cloud-base horizontal wind speed and cloud duration yields the cloud chord length.

Case selection
A total of 128 ShCu days were selected between 2011 ~ 2014.The selection criterion is in principle similar to previous studies (Zhang and Klein, 2013;Lareau et al., 2018), which involves both objective and subjective criteria.The objective criteria include three steps: (1) the cloud-base height (defined as the mean of the lowest quartile within the 2-h period) has to be within 30% of LCL to ensure coupling, (2) the KAZR reflectivity cannot exceed 0 dBZ between the surface and cloud base to ensure no considerable precipitation, and (3) the cloud duration cannot exceed 30 min to exclude stratiform clouds.Besides, we examine KAZR imageries to ensure ShCu-like characteristics.This is the best we can do since a completely objective method for selecting ShCu remains missing, although the emerging new technique of machine learning is promising to address this issue in the near future (Rasp et al., 2019).
Based on these criteria, we obtain 32 ShCu days per year, similar to the 28 ShCu days per year in Zhang and Klein (2013) and Lareau et al. (2018), suggesting that there is no marked sampling difference between this study and previous ones.Fig S3 shows the statistics of these selected ShCu ensembles.On average, each ensemble is composed of ~ 20 individual ShCu, with half lasting longer than 30 secs.The majority of the ensembles have the maximum cloud chord length shorter than 5 km, consistent with prior knowledge. 1/2 for different wcrit.

Sub-cloud turbulence explains cloud-base updrafts
Overall, the (TKE w M) 1/2 is a good predictor of cloud-base updrafts, explaining ~ 50% of their variabilities.Note that the degree of scattering is still noticeable, but given the instrument error of the DL (~ 0.1 m/s) and potential sampling errors due to the assumption of stationarity, such degrees of correlation are good enough for demonstrating the physical validness.To our knowledge, this is the first observational evidence supporting the ability of the sub-cloud turbulence to dictate cloud-base updrafts that was only found in high-resolution models (Grant and Brown, 1999;Fletcher and Bretherton, 2010;van Stratum et al., 2014).Such good correlations suggest a continuity of vertical momentum between the sub-cloud layer and cloud base, despite the inbetween weakly stable layer (i.e.cloud-base transition layer) (Neggers et al., 2007;Stevens, 2007).
Indeed, the stability of the transition layer interacts with the convective circulation, a manifestation of the dynamical coupling between the sub-cloud and cloud layers, to reach an equilibrium that maintains the mass conservation (Neggers et al., 2006;Fletcher and Bretherton, 2010).In this regard, the transition layer property should not be considered an external forcing that alters the coupling between the sub-cloud and cloud-base dynamics, but an internal parameter that responds to the circulation.
Both   ̅̅̅̅ and    ̅̅̅̅̅̅ increase with the wcrit, but the    ̅̅̅̅̅̅ shows much weaker sensitivity primarily because the    ̅̅̅̅̅̅ gives more weight to the larger vertical velocities.The intercepts also increase with wcrit, which is an artificial consequence of using non-zero wcrit.Physically speaking, a zero TKE w M should lead to zero cloud-base updraft speed.Therefore, we will focus our subsequent discussions on the slopes that bear more physical meaning than intercepts.
To compare our results with that from FB10, we visualize the Eq.(1) in Figure 3a (light blue curve).FB10 uses the wcrit of 0.5 m/s.Our empirical estimate (the red line) shows a stronger sensitivity of   ̅̅̅̅ to the sub-cloud turbulence than FB10 by more than a factor of 3. What causes the difference?One possible reason is that we used the TKE w M that does not include the horizontal components of the TKE, leading to smaller values of TKE and, thus, a steeper slope.Another more likely reason is that the horizontal resolutions of the model used by FB10 are too coarse (1 km) to accurately simulate the vertical velocities.For instance, modeled vertical velocities decrease with the model resolution by a power law of -2/3 (Rauscher et al., 2016;Donner et al., 2016).The underestimated   ̅̅̅̅ due to low resolution may flatten the slope of   ̅̅̅̅ versus (TKEML) 1/2 in FB10.
To understand which factor is responsible, we use the LES data of 18 ShCu days from the LASSO project (Text S1).The LASSO horizontal resolution is 100 m, 10 times finer than that used in FB10.With the model output of three-dimensional winds, we are able to diagnose the full components of TKEML so that we can conduct an "apple-to-apple" comparison between the LASSO and FB10.As shown by the green lines in Fig. 3a, LASSO models (WRF and System for Atmospheric Modeling, SAM) show slopes steeper than the FB10 by more than a factor of 3 (see Fig. S4 for their scatter plots with statistical details).This confirms that the flatter slope of FB10 is likely caused by the coarse model resolution.The comparison between the LASSO and DL, which is not the focus of this study, is discussed in the supplementary material (Text S2).
We have tabulated the empirical formulas for   ̅̅̅̅ and    ̅̅̅̅̅̅ for different wcrit (Table S1) so that readers can use what suits their research interests.

4.2.Diurnal dependence
Given that all cases are in the early afternoon, one may ask how the observed relationship is representative of the other times of a diurnal cycle.To address this question, we use the LAASO data to examine its diurnal dependence.We chose the wcrit = 0 m/s for determining the   ̅̅̅̅ because, as noted above, using an ad-hoc wcrit, say 0.5 m/s, leads to a markedly positive   ̅̅̅̅ for zero 1/2 .By using wcri = 0 m/s, we can force the best-fit line through the origin through the least-square algorithm, freeing us from the unphysical meaning of positive intercepts.Figure 4a and b show the scatterplots of the   ̅̅̅̅ versus (TKE w M) 1/2 in different local times simulated by WRF and SAM, respectively.Both models show notably significant correlations between the two quantities in different phases of a diurnal cycle, confirming the ability of (TKE w M) 1/2 to explain the variability of   ̅̅̅̅.More importantly, the slope of the relationship varies little with local time, except in the early morning and late afternoon (Fig. 4c and d).In the early morning, the stronger capping inversion weakens the speeds of rising thermals when they penetrating into the inversion, leading to smaller   ̅̅̅̅ for given sub-cloud turbulence (Fig. S1c).Such a stabilization effect becomes less influential as the convection kicks up, which lessens the inversion strength.In the late afternoon, as the solar insolation weakens, the surface fluxes decrease considerably whereas the boundary layer remains deep (Fig. S1d).This leads to a decoupling between the ShCu and the surface (Stull, 2012), which may explain the flatter slope between   ̅̅̅̅ and (TKE w M) 1/2 in the late afternoon.
In summary, the diurnal dependence of the coupling between the wb and sub-cloud turbulence is small, except in the early morning and late afternoon when the strong capping inversion and cloud-surface decoupling may lead to flatter slopes, respectively.

Conclusion
This study examines the relationship between the sub-cloud turbulence and cloud base updrafts using Doppler lidar (DL) observations of 128 shallow cumulus (ShCu) ensembles over the Southern Great Plains.We proposed a new DL sampling method that allows measuring the cloud-base updrafts for an ensemble, instead of individual, ShCu.Specifically, we take advantage of the stationarity and ergodicity of ShCu-topped boundary layers in the early afternoon when the temporal change in the surface forcing is minimum.For each ShCu case, we selected a 2-hour window of DL that includes an average amount of ~ 20 individual cumuli with varying sizes, constituting an ensemble.This allows us to compute the ensemble-averaged quantities from DL measurements made at a fixed point.By analyzing the 128 ShCu ensembles, we found that the vertical velocity variance explains ~ 50% variability of ensemble-mean cloud-base updrafts, thus supporting the widely-held hypothesis and practice of using the sub-cloud turbulent kinetic energy to parameterize the cloud-base updrafts in some state-of-the-art mass flux closure schemes of convection parameterization (Bretherton et al., 2004;Neggers et al., 2009;Fletcher and Bretherton, 2010).To our knowledge, this is the first observational evidence that demonstrates the ability of showing a ShCu ensemble simulated by the Weather Research and Forecasting (WRF) in the Large-Eddy Simulation (LES) Atmospheric Radiation Measurements (ARM) Symbiotic Simulation and Observation (LASSO) project (Text S1).The surface fluxes and large-scale forcing are uniform over the 14.4 × 14.4 km domain with a horizontal grid size of 100 m.The vertical velocity field at the cloud-base level shows a

Figure 1 :
Figure 1: Examples of the different length scales of spatial variability of wb and TKEML using WRF-simulated ShCu on 21 UTC, June 6, 2015.(a) Spatial distribution of vertical velocity at the cloud-base level with maximum cloud coverage.Black contours mark the cloudy regions with liquid water content greater than 0.01 g/m 3 .(b) The same scene but the color shading is the TKEML.(c) Scatter plot of cloud-base vertical velocity versus TKEML, with each point representing mean over individual cumuli.The size of a point is proportional to the size of cumuli.The data are obtained from the first phase of LASSO project.The TKEML is computed as 0.5( ′ 2 +  ′ 2 +  ′ 2 ) averaged below the cloud base.

Figure 2 :
Figure 2: An example case of the shallow cumulus field on Jun 10, 2012, over the SGP.(a) MODIS image centered on the SGP site (red star) at ~20:30 UTC.The red solid line marks the rough direction and travel distance of the mean horizontal wind during the 19 ~ 21 UTC.(b) Height-time plot of Doppler lidar image of vertical velocity during a twohour window from 19 to 21 UTC.The black dots mark the cloud-base heights measured by a ceilometer.The blue rectangle marks a smaller window shown in the (c).Navy contours mark the cloudy regions defined as groups of pixels with reflectivity greater than 10 -4.6 m -1 sr -1 .

Figure
Figure 3 shows the scatter plots of   ̅̅̅̅ (a) and    ̅̅̅̅̅̅ (b) versus (TKE w M)

Figure 4 :
Figure 4: Scatterplots of   ̅̅̅̅ (wcrit = 0 m/s) versus the (TKE w M) 1/2 grouped by the local standard time, simulated by WRF (a) and SAM (b).Each group of points corresponds to a bestfit linear regression line forced through zero.The slopes of the best-fit lines are plotted in (c) and (d) for WRF and SAM, respectively.
sub-cloud turbulence intensity to dictate the cloud-base updrafts.With the observational data, we derived empirical relationships between the square-root of sub-cloud turbulent kinetic energy and ensemble-mean cloud-base updraft speeds that are computed for different thresholds of vertical velocity and by different averaging schemes.Although all the 128 cases were sampled in the early afternoon, the diurnal variation of the relationship is weak (except in the early morning and late afternoon), as shown by the LES simulations of 18 ShCu cases over the SGP.These empirical formulas are useful for the developments of cumulus parameterizations, theoretical studies of ShCu dynamics, and satellitebased inference of cloud-base updrafts.

Figure S1 :
Figure S1: WRF-simulated composite diurnal variations of t * (a), TKEML (b), and height-time plots of Brunt-Vaisala frequency (c), and vertical velocity variance (d).In (c) and (d), the black lines mark the diagnosed mixed-layer height (h).All plotted are the composite means of the 18 ShCu cases from the 1 st phase of the LAASO project.

Figure S2 :
Figure S2: Probability density functions of vertical velocity for pixels at 100 m below (black) and above (red) the cloud bases.

Figure S3 :
Figure S3: Statistical distribution of key quantities for the 128 ShCu cases including (a) the number of individual cumuli (the red marks those that last longer than 30 secs), (b) maximum cloud duration, (c) horizontal wind speed near cloud base, and (d) maximum cloud chord length.

Table S1 :
Statistics of the relationships between the ensemble-mean cloud-base updrafts and