Lei ZHU ,Chunsong LU* ,Xiaoqi XU ,Xin HE ,Junjun LI ,Shi LUO ,Yuan WANG,and Fan WANG

1Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters (CIC-FEMD) and Key Laboratory for Aerosol-Cloud-Precipitation of China Meteorological Administration, Nanjing University of Information Science and Technology, Nanjing 210044, China

2Nanjing Joint Institute for Atmospheric Sciences, Nanjing 210019, China

3College of Aviation Meteorology, Civil Aviation Flight University of China, Guanghan 618307, China

4Collaborative Innovation Center for Western Ecological Safety, Lanzhou University, Lanzhou 730000, China

5Department of Geography, Hong Kong Baptist University, Hong Kong SAR 999077, China

ABSTRACT The process of entrainment-mixing between cumulus clouds and the ambient air is important for the development of cumulus clouds.Accurately obtaining the entrainment rate (λ) is particularly important for its parameterization within the overall cumulus parameterization scheme.In this study,an improved bulk-plume method is proposed by solving the equations of two conserved variables simultaneously to calculate λ of cumulus clouds in a large-eddy simulation.The results demonstrate that the improved bulk-plume method is more reliable than the traditional bulk-plume method,because λ,as calculated from the improved method,falls within the range of λ values obtained from the traditional method using different conserved variables.The probability density functions of λ for all data,different times,and different heights can be well-fitted by a log-normal distribution,which supports the assumed stochastic entrainment process in previous studies.Further analysis demonstrate that the relationship between λ and the vertical velocity is better than other thermodynamic/dynamical properties;thus,the vertical velocity is recommended as the primary influencing factor for the parameterization of λ in the future.The results of this study enhance the theoretical understanding of λ and its influencing factors and shed new light on the development of λ parameterization.

Key words: large-eddy simulation,cumulus clouds,entrainment rate,probability density functions,spatial and temporal distribution

1.Introduction

Cumulus clouds have important effects on the vertical transport of heat,moisture,and momentum and play an important role in affecting the energy budget of the Earth–atmosphere system and climate change (Guo et al.,2014,2015b;Donner et al.,2016;Wang et al.,2018;Yang et al.,2019;Jeyaratnam et al.,2021;Sheng et al.,2022).In weather and climate models,because model grid spacing is insufficient to resolve convective processes,the parameterization of cumulus clouds must be achieved using empirical hypotheses and physical variables at grid points (Zhang and McFarlane,1995;Lin et al.,2000;Donner et al.,2001;Wu,2012;Xie and Liu,2013).The quality of the cumulus parameterization scheme is crucial for simulating precipitation,monsoons,Madden-Julian oscillations,and tropical cyclones (Zou and Zhou,2011;Del Genio et al.,2012;Yang et al.,2015;Wang et al.,2017;Zhao et al.,2018a).Entrainment-mixing between cumulus clouds and ambient air is a key process that affects cloud development and the formation of precipitation,by causing cloud droplet evaporation to alter temperature,humidity,buoyancy (B),and vertical velocity (w) in clouds (Jonas,1990;Houze,1993;Telford,1996;de Rooy et al.,2013;Yeom et al.,2017;Wang,2020;Luo et al.,2022;Wang et al.,2023).In many cumulus parameterization schemes,the lateral entrainment rate (λ) is a physical variable characterizing the strength of entrainment,which has a significant influence on climate sensitivity,precipitation,and monsoons (Klocke et al.,2011;Yang et al.,2013,2021;Zhao,2014;Lu and Ren,2016;Hanf and Annamalai,2020).

To better understand the process of entrainment-mixing and improve cumulus parameterization schemes in models,an accurate calculation ofλis essential.Stommel (1947)first calculatedλusing the temperature and specific humidity inside and outside the cloud.Based on this,Yanai et al.(1973) and Betts (1975) proposed the bulk-plume method to calculateλaccording to conserved physical variables (i.e.,total moisture,liquid water potential temperature (θl),and moist static energy) in clouds and ambient air.The bulkplume method has been used in many studies.For example,Esbensen (1978) used this method to calculate the observedλin large-scale shallow cumulus clouds,and this was also applied to calculateλfrom aircraft observations (Neggers et al.,2003;Gerber et al.,2008;Lu et al.,2018) and satellite data (Luo et al.,2010;Takahashi and Luo,2012;Li et al.,2022;Takahashi et al.,2023).Romps (2010) proposed a method to estimateλby directly calculating the amount of air entrained into the cloud in large-eddy simulation (LES)experiments and found that the results were approximately twice those from the bulk-plume method.Dawe and Austin(2013) noted that the reason for the difference may be related to the presence of a moist cloud shell around the cloud core that is more humid than the ambient air.Nevertheless,accurately calculatingλremains a challenge,hindering the improvement of cumulus parameterization schemes.

The parameterization ofλis an important way to describe the process of entrainment-mixing in models (de Rooy et al.,2013;Lu et al.,2016;Zhang et al.,2016;Villalba-Pradas and Tapiador,2022).However,the process of entrainment-mixing is also often considered to be a stochastic process.Romps and Kuang (2010) found that a parcel model with a constant or continuousλcould not reproduce the observed variability of clouds,whereas a stochastic parcel model could,in which entrainment behaves like a stochastic Poisson process.Observations have shown that the probability density functions (PDFs) ofλfor shallow and deep convection follow a log-normal distribution (Lu et al.,2012a;Guo et al.,2015a),supporting the stochastic process of entrainment-mixing.Romps and Kuang (2010) and Romps (2016)presented Lagrangian and Eulerian implementations of the stochastic parcel model,respectively,to treat the process of entrainment-mixing in cumulus clouds as a stochastic process.Böing et al.(2014) studied the process of entrainmentmixing using two parcel models that describeλin different ways.The stochastic mixing model withλsatisfying a gamma distribution better captured the scatter of the incloud thermodynamic properties.Yang et al.(2021) linked the mass fluxes between shallow and deep convection by assuming thatλsatisfies a log-normal distribution,and a modified deep convection scheme improved precipitation simulations in mean state and variability at various timescales.However,a simply assumed log-normal distribution ofλmay cause uncertainty since this assumption only considers statistical characteristics.Therefore,analyzing the PDF ofλand its influencing factors provides an important reference for numerical models to improve cumulus parameterization schemes.

Despite significant progress,further improvement of the calculation accuracy ofλand determination of the PDF ofλand its influencing factors remain worthy of further study.This study developed an improved bulk-plume method and applied it to a simulated cumulus cloud case using an LES.First,the accuracy of the method for calculatingλwas evaluated.Second,the PDFs ofλwere fitted to a lognormal distribution,then the variation in the fitting parameters over time and the influencing factors were analyzed.Finally,the variation in the PDFs ofλwith height was analyzed.The LES was chosen because it is difficult to obtain the PDFs ofλat different heights and times owing to limitations to aircraft observations.An LES can provide the threedimensional structure of cumulus convection and its temporal evolution (Neggers et al.,2003;Endo et al.,2015),which provides important data for studyingλ(Romps,2010;Dawe and Austin,2013;Drueke et al.,2019).

The remainder of this paper is organized as follows.Section 2 introduces the model setting and simulation results.Section 3 describes and verifies the improved method for calculatingλ.Section 4 analyzes the spatial and temporal distribution ofλ,the variation in PDFs ofλwith time and height,and its influencing factors.The conclusions and discussions are presented in section 5.

2.Simulation results

The LES model was used to simulate cumulus clouds over the Southern Great Plains (SGP) during the LES Atmospheric Radiation Measurement Symbiotic Simulation and Observation (LASSO) experiment (Gustafson et al.,2020)on 11 June 2016 [Coordinated Universal Time (UTC) and Local Standard Time (LST,LST=UTC– 6)] (Xu et al.,2022).The model adds time-varying large-scale forcing(Endo et al.,2015) to a Weather Research and Forecasting model tailored for solar irradiance forecasting (WRF-Solar),which can better handle radiation-related processes,as well as cloud-aerosol-radiation feedback (Haupt et al.,2016;Jimenez et al.,2016).The large-scale forcing data and model settings were the same as those used by Xu et al.(2022).The horizontal grid spacing of the model was 100 m × 100 m,the number of grids was 144 × 144,and the vertical grid spacing was 30 m with 225 layers.The simulation started at 1200 UTC on 11 June 2016,and ended at 0300 UTC on 12 June 2016.The output of simulation results occurred every 10 minutes.The quasi-steady state in an LES may be reached within a few large-eddy turnover times(Moeng and Sullivan,1994;Nakanish,2001).Here,the spinup time was set to 3 h (Xu et al.,2022),which was more than eight large-eddy turnover times.The chosen microphysical scheme was the Thompson aerosol-aware scheme (Thompson and Eidhammer,2014) with an improved parameterization of the entrainment-mixing mechanism (Luo et al.,2020;Xu et al.,2022).

For this cumulus cloud case,the diurnal variations of the simulated cloud fraction (CF) and liquid water path(LWP) in Fig.1a illustrate a first peak around 1600–1700 UTC and a secondary peak around 2000–2100 UTC.Since LWP is a vertically integrated quantity and CF is a horizontal coverage value,the peak timings of the two variables do not always coincide (Shin et al.,2021).The mean thermodynamic profiles (Fig.1b) show that the mean characteristics of the convective planetary boundary layer are as follows: the surface layer from 0 m to about 330 m,the mixed layer from about 330 m to 1600 m,and the inversion layer from about 1600 m to 2800 m.The maximum updraft in the simulation is 9.40 m s–1,and the mean of the maximum updraft for each moment is 6.19 m s–1.The simulated cumulus clouds formed at approximately 1500 UTC and dissipated completely at approximately 2400 UTC.The simulations generally captured the diurnal variation of cumulus clouds (Xu et al.,2022).

Fig.1.(a) Time series of the cloud fraction (CF) and liquid water path (LWP) from the simulation on 11 June 2016.(b) The simulated mean thermodynamic profiles of virtual potential temperature (θv) and total water vapor mixing ratio (qt).(c) Distribution of the grid numbers of cloud grid points with time and height.The black,dash-dotted lines divide the three stages of cloud ensemble evolution: formation (1500–1900 UTC),maintenance (1900–2200 UTC),and dissipation (2200–0020 UTC the next day).The blank values in (c) indicate that no grid point of the area meets the criteria for a non-precipitating cloud.

Cloud grid points without precipitation were selected based on the following criteria: their liquid water mixing ratios (ql) must exceed 0.01 g kg–1,Bmust be greater than 0 m s–2,and rain-water mixing ratios must be less than 0.005 g kg–1(Lu et al.,2012b).Figure 1c shows the timeheight distribution of the number of cloud grid points in the simulated region.For the case of cumulus cloud,the height of the cloud base continued to rise from 1500 to 1900 UTC,and the cloud ensemble vigorously developed within a large number of cloud grids,accompanied by the decay of small cumulus clouds and the appearance of new clouds;after 1900 UTC the number of cloud grid points declined and after 2200 UTC the cumulus cloud ensemble gradually dissipated.Based on the temporal and spatial distribution characteristics,the simulated cumulus clouds were divided into three stages: formation (1500–1900 UTC),with continuous increases in cloud base and top heights,maintenance(1900–2200 UTC),with relatively consistent cloud base and top heights,and dissipation (2200–0020 UTC the next day),with decreasing cloud base and top heights.method (Betts,1975;Neggers et al.,2003;Gerber et al.,2008):

where ϕ is a conserved physical variable during cloud ascent,the subscripts c and e represent the cloud and ambient air,respectively,andzis height.Substituting the total water mixing ratio (qt) andθlas the conserved variables into Eq.(1) yields:

The total water vapor mixing ratio is defined as:

whereqvis the water vapor mixing ratio.θlis defined as(Betts,1973):

whereθ=γTis the potential temperature withγ=(p0/p)0.286;p0andpare standard atmospheric pressure and air pressure,respectively;Tis temperature;Lvis latent heat;andcpis the specific heat capacity at constant pressure.

Substituting Eq.(4) and Eq.(5) into Eq.(2) and Eq.(3)yields:

3.An improved method for calculating λ

3.1.Development of the improved method

The parameterλcan be calculated using the bulk-plume whereqvcin the cloud is the saturated water vapor mixing ratio calculated according to the Clausius-Clapeyron equation(Wallace and Hobbs,2006).Theqlecancels and does not appear in Eq.(6) considering the nature of ambient air.Combining Eqs.(6) and (7) withqlcin the cloud at different heights and meteorological properties of ambient air,bothTcandλin the cloud can be calculated simultaneously by iteration.Note that the input variables of this method areqlcin the cloud andqve,Te,andpeof ambient air;the output variables areqvc,Tc,andλ;and the remaining variables are intermediate variables.

The improved method does not depend on which conserved variable is used to calculateλ.On the contrary,the iteration of the two equations takes into account the simultaneous conservation of the two physical variables,making the calculation ofλmore accurate.For cumulus clouds with multilayer data (in this study referring to LES simulations),λis calculated layer-by-layer from the cloud base to the cloud top,allowing for the vertical distribution ofλto be obtained.

3.2.Validation of the improved method

The variables required to calculateλfrom the model results were determined as follows: First,the domain-averagedqve,Te,andpewere taken as the properties of the ambient air involved in the process of entrainment-mixing (Lu et al.,2012b).Second,an averageqlcof the cloud grid points with positiveB(the virtual potential temperature in the cloud exceeds that of the environment) in each layer of the domain was considered as the cloud properties of that layer.Finally,the cloud base height at each moment was the lowest height of the cloud grid points in the domain at that moment.

Before analyzingλ,as calculated using the improved method,it is necessary to examine the accuracy of the improved method with respect to the bulk-plume method.Figure 2 compares the PDF ofλcalculated by the improved method described in section 3.1,with that calculated by the traditional bulk-plume method [Eqs.(2) and (3)].The results show thatλcalculated by the improved method was relatively large in the lower layers,decreased gradually with increasing height,and increased gradually in the upper layers(Fig.2b).Theλcalculated from the traditional bulk-plume method usingqtwas smaller than that calculated using the same method but withθl.Interestingly,λfrom the improved method was between those of the traditional method usingqtorθl.The mean value ofλcalculated by the improved method was 1.93 km–1,and the mean value ofλcalculated fromqtandθlwas 1.34 km–1and 2.16 km–1,respectively.Therefore,the improved method can reduce the uncertainty associated with traditional methods.

Fig.2.Probability density function (PDF) of entrainment rate (λ) as a function of height calculated using (a) the bulk-plume method using total water mixing ratio (qt),(b) the improved method,and (c) the bulk-plume method using liquid water potential temperature (θl).The mean values of λ calculated by the three methods are shown in the figure.

4.Spatial and temporal distribution of λ and its influencing factors

4.1.Overall characteristics of λ

Figure 3 shows the spatial and temporal distributions ofλafter examining the accuracy of the improved method for calculatingλ.In the formation and dissipation stages,λexceeded that in the maintenance stage,that is,λdecreased at first and then increased with time.In the vertical direction,λwas generally large in the lower and upper layers,but smaller in the middle layers.The average profile ofλalso demonstrates thatλfirst decreased and then increased along with the increase in height (Fig.2b) and the overall vertical distribution ofλthat decreased with height was consistent with that obtained in previous studies (Dawe and Austin,2013;Xu et al.,2021).However,λexhibited some differences at different stages of the cumulus ensemble life cycle.During formation,λwas relatively large near the cloud base and initially decreased before increasing with increased height.During the maintenance stage,the vertical distribution ofλwas relatively consistent,and the overall trend featured a decrease ofλwith increasing height and then increasing,whereasλin the middle layer (～2000–2200 m) changed only slightly with height.During dissipation,λincreased with increasing height.

In the formation stage,the cloud ensembles experienced significantλover a large depth due to the rapid development accompanied by the decay of small cumulus clouds and the appearance of new clouds,some of which were initial cumulus clouds;the main reason is that the cloud ensemble in the formation stage possessed a smallerwthan the subsequent maintenance stage [Fig.S1b in the electronic supplementary material (ESM)] sinceλfrequently increases with decreasingw(Neggers et al.,2002;de Rooy et al.,2013;Lu et al.,2016;Xu et al.,2021) andλhas the strongest relationship withwamong other thermodynamic/dynamical properties(see sections 4.2.2 and 4.3.2 for detailed analysis).In contrast,the largerwin the maintenance stage (Fig.S1b in the ESM) leads to the development of smallerλin the lower part of the cloud (Fig.3).The cloud dissipation was responsible for the increasingλnear cloud tops and toward the end of the simulation,as the dissipating clouds share similar physical properties with the ambient air [resulting in small ϕe–ϕcin Eq.(1)].It is also worth noting that smallerwoccurred in the cloud ensembles near the cloud tops and toward the end of the simulation (Fig.S1b;in the ESM),resulting in largerλ(Fig.3).The increasingλtowards the cloud top has been found in many model simulations (Romps,2010;Dawe and Austin,2011,2013;Xu et al.,2021).

Fig.3.Distribution of entrainment rate (λ) with time and height.The black dash-dotted lines divide the three stages of cloud ensemble evolution: formation (1500–1900 UTC),maintenance (1900–2200 UTC),and dissipation (2200–0020 UTC the next day).The blank values indicate that no grid point of the area meets the criteria for a non-precipitating cloud.

Dawe and Austin (2013) diagnosed the spatial and temporal distributions ofλfor each individual shallow cumulus cloud using a cloud-tracking algorithm in an LES.The results for one of the cumulus clouds showed thatλwas larger in the lower levels and gradually decreased with height,and there were several large values in the upper layers.Note thatλinitially decreased and then increased with time.The spatial and temporal distributions of the cumulus ensemble in the present study represented the ensemble characteristics of all individual cumulus clouds and were very similar to the results of individual cumulus clouds in Dawe and Austin (2013).This indicates that the spatial and temporal distributions of the cumulus ensembleλsimulated in this study can also represent the characteristics of individual cumulus clouds.In addition,the variation with height in PDF ofλacross all cumulus clouds simulated by Dawe and Austin(2013) demonstrated that ensembleλinitially decreased and then increased with height,consistent with the vertical variation of ensembleλsimulated in the present study.

4.2.Variation of the PDF of λ over time and its influencing factors

4.2.1.Variation of the PDF over time

Figure 4 shows the PDF of 1293 domain-averagedλvalues at all heights and times.Due to the assumptions described in section 3.2,we note that the calculatedλis the cloud ensembleλ.The entrainment characteristics of cloud ensembles (e.g.,cloud ensembleλ) are required in current general circulation models (GCMs) (Dawe and Austin,2013).The PDF ofλwas well fitted using a log-normal distribution(coefficient of determinationR2=0.99);similar results were obtained by Lu et al.(2012a).The log-normal distribution was chosen as the fitting function because previous studies have demonstrated that the PDF ofλconforms to a log-normal distribution (Lu et al.,2012a;Dawe and Austin,2013;Guo et al.,2015a).In climate models,the grid spacing is coarse and cumulus convection is often described through parameterization.The different effects of entrainment on cumulus convection can be obtained by parameterizingλrather than by assuming a constantλ(Gregory,2001;Neggers et al.,2002;Wu,2012).Therefore,a wide distribution ofλis reasonable.It is also important to treat the process of entrainment-mixing as a stochastic process (Romps and Kuang,2010;Böing et al.,2014;Romps,2016;Yang et al.,2021).The PDF ofλis well described by the log-normal distribution (Fig.4),which supports stochastic entrainment and has important implications for treating convection in models.

Fig.4.Probability density function of the entrainment rate (λ)at all times and heights.The number of samples,the coefficient of determination (R2),mean (μ),and standard deviation (σ) of ln(λ) for the log-normal fit (red line) of all data are provided.

Fig.5.Panels (a–i) represent the probability density function of the entrainment rate (λ) in hourly intervals from 1500 to 2400 UTC.The number of samples,coefficient of determination (R2),mean (μ),and standard deviation (σ) of ln(λ) for the log-normal fit (red line) in each time interval are provided.

The PDF of cloudλfor all times and heights in Fig.4 can present the characteristics of ensembleλduring the entire life cycle of the cumulus ensemble.It is also interesting to see whether there are differences in the PDFs ofλat different times,which could further improve the description of the entrainment process in the cumulus parameterization scheme.Figure 5 shows the PDF for all cloudλfor each hour from 1500 to 2400 UTC,along with the parameters of the log-normal fit andR2.For the consistency of the time interval,the data of the two moments after 2400 UTC (13 samples) were not involved in the analysis.TheR2in each time interval exceeded 0.90 (more than half attain 0.95),consistent with the results for all times and heights (Fig.4).Furthermore,λin each time interval was also described well by a log-normal distribution,further confirming the stochastic distribution ofλ.

Differences still existed between the distributions ofλfor each time interval.The two parameters of the log-normal fit ofλwere the mean (μ) and standard deviation (σ) of ln(λ).Theμinitially decreased and then increased with time,andσinitially increased with time and decreased later.Thus,in the early stage of cloud ensemble development,λwas relatively large and the distribution was relatively concentrated.Over time,λinitially decreased and then increased,and its distribution gradually widened.At the end of the cloud life cycle,λbecame relatively large and concentrated.

4.2.2.Influencing factors

An accurate description of the PDF ofλin the model is of value for improving the simulation,which requires investigation of the factors affecting the PDF ofλ.Previous studies have shown thatw,B,and the environmental relative humidity(RHe) are important factors,and these variables are often used to parameterizeλ(Stirling and Stratton,2012;de Rooy et al.,2013;Lu et al.,2016,2018;Zhang et al.,2016;Bera and Prabha,2019;Xu et al.,2021).Therefore,the effects ofw,B,and RHeon the PDF ofλare analyzed in this section.Note thatw,B,and RHerefer to the mean values of the cloud grid points as mentioned in section 2.

Similar to the analysis of the PDF ofλ,PDFs were also calculated forw,B,and RHeat each moment and fitted using a log-normal distribution.The hourly PDFs ofwfitted well to a log-normal distribution (Fig.S2;in the ESM),and the results for the other variables were similar (figures not shown).This facilitated the analysis of the relationship between these variables andλusing the fitting parameters of the log-normal fit.

To determine the factors influencing the PDF ofλ,the temporal variation of the parameters in the log-normal fit of the hourly PDF ofλ,w,B,and RHewas examined(Fig.6).Consistent with the analysis in section 4.2.1,the mean value ofλinitially decreased and then increased with time (Fig.6a).In contrast,winitially increased and then decreased with time (Fig.6b);thus establishing a negative correlation withλ.Bgenerally decreased with time (Fig.6c)and was not correlated withλ.The RHeinitially decreased and then increased with time (Fig.6d) and was positively correlated withλ,although the timing of inflection points,from decreasing to increasing values,differed forλand RHe.

Fig.6.Mean (μ) for the log-normal fit of the hourly probability density functions of (a) entrainment rate (λ),(b)vertical velocity (w),(c) buoyancy (B),and (d) environmental relative humidity (RHe),as a function of time.Error bars represent the standard deviation (σ).

In summary,λwas negatively correlated withw,consistent with the results of many previous studies (de Rooy et al.,2013;Lu et al.,2016).Note that there was no significant correlation betweenλandB,although previous studies have found a negative correlation between them (Lin,1999;von Salzen and McFarlane,2002).However,a lack of a significant correlation betweenλandBwas also found by Romps(2010).Gregory (2001) usedB/w2to parameterizeλ,a scheme that has been applied to climate models (Kim and Kang,2012;Song and Zhang,2018),showing thatλis positively correlated withBand negatively correlated withw.In previous studies,λandBwere found to be positively correlated,negatively correlated,or not correlated,suggesting that the effect ofBonλmight be indirect.In contrast,whas a more direct effect onλbecause whenwis smaller,the cloud is allotted more time to mix with ambient air during ascent,which results in a largerλ(Neggers et al.,2002;Lu et al.,2016;Zhang et al.,2016;Xu et al.,2021).Lu et al.(2016) demonstrated that a largerλreduces the temperature in the cloud when the ambient air and the cloud mix,thereby reducingB,while the decrease inBdrives a decrease inw,which leads to an increase inλ.They also noted thatwis the optimum choice when using a single variable to parameterizeλ.Previous studies have shown that the relationship betweenλand RHecan be positive (Axelsen,2005;Lu et al.,2018;Bera and Prabha,2019;Stanfield et al.,2019;Zhu et al.,2021) or negative (Bechtold et al.,2008;Zhao et al.,2018b);and the results of this study tend to suggest a positive correlation.In conclusion,a relationship betweenwandλis optimal.

4.3.Variation of the PDF of λ with height and its influencing factors

4.3.1.Variation of the PDF with height

Section 4.2 discussed variations in the PDF ofλover time.Because clouds andλexhibit different distribution characteristics at different heights (Figs.1–3),it is necessary to examine the variation of the PDF ofλwith height.Similar to section 4.2.1,PDFs were calculated forλper 200 m in the range 1000 m to 2600 m above sea level and fitted using a log-normal distribution.The PDFs ofλwere well fitted in each height range (Fig.7;R2exceeded 0.96,except for the PDF of 2400–2600 m),which further confirms that entrainment is a stochastic process.In the lower layer (Figs.7a,b),μwas relatively large andσwas relatively small,indicating thatλof the lower layer was relatively large and the distribution was relatively concentrated.As height increased,μinitially decreased and then increased,whereasσincreased gradually.Correspondingly,λinitially decreased and then increased with height,and the distribution became increasingly dispersed.Interestingly,the PDF was narrow at lower altitudes (mainly due to the identical properties of most cloud ensembles in the formation stage) and became wider as height increased (mainly due to cloud variability brought forth by the coexistence of cloud ensembles in multiple stages).The PDFs ofλat all heights in an LES,as presented by Dawe and Austin (2013),showed thatλstrongly followed a log-normal distribution and the mean PDF of all individual cumulus cloudλwas consistent with the cloud ensembleλ.Furthermore,λinitially decreased with height and then remained unchanged,and increased at high levels,similar to the vertical distribution obtained in this study.

Fig.7.Panels (a–h) represent the probability density function of the entrainment rate (λ),in 200-m increments,over the range of 1000 to 2600 m.The number of samples,coefficient of determination (R2),mean (μ),and standard deviation (σ) of ln(λ) of the log-normal fit (red line) in each 200 m are provided.

4.3.2.Influencing factors

The factors influencing the PDF ofλwere analyzed in section 4.2.2.In this section,we discuss whether these same conclusions apply to the PDFs ofλat different heights.Similar to section 4.2.2,log-normal distributions were fitted to the PDFs ofw,B,and RHefor each 200-m range.The PDFs ofwwithin each 200-m range were well fitted by a log-normal distribution (Fig.S3;in the ESM) and the results for the other variables were similar (figures not shown).Figure 8 shows variations in the fitting parameters of the PDFs ofλand its influencing factors,as a function of height,to analyze the relationship betweenλand influencing factors at different heights.

The mean value ofλinitially decreased and then increased with height (Fig.8).In contrast,winitially increased and then decreased with increasing height,andλandwwere negatively correlated.Binitially increased and then decreased with increasing height and was negatively correlated withλ.However,there was a small value in 1800–2000 m,which was inconsistent with the vertical change inλ.The vertical profiles ofλandBhave a better correspondence compared with the time-series plots ofλandBin Fig.6.The RHeinitially decreased and then increased with height,consistent with the vertical distribution ofλ;the corresponding relationship between these two variables was also better than that apparent in Fig.6.In conclusion,the relationship betweenλand each influencing factor at different heights is similar to that at different times (section 4.2.2).Furthermore,similar to Fig.8,the vertical distributions of the fitting parameters for the PDFs ofλand its influencing factors were plotted for the three stages of the cloud ensembles (Fig.S1;in the ESM).Comparisons of the results for the three stages and that for all times show similar vertical distributions,andwhas the best relationship withλin the three stages.Combining the results of different heights and times,the relationship betweenλandwis the strongest.

Fig.8.Mean (μ) of the log-normal fit of the probability density functions of (a) entrainment rate (λ),(b) vertical velocity (w),(c) buoyancy (B),and (d) environmental relative humidity (RHe) in each 200-m layer as a function of height.Error bars represent the standard deviation (σ).

5.Conclusion and discussions

The entrainment rate (λ) is an important physical variable in cumulus parameterization schemes.The accurate calculation ofλfrom observations or high-resolution simulations is key to improvingλparameterization.The probability density function (PDF) ofλcan be very useful for treating cumulus convection in models.Therefore,this study applied a largeeddy simulation of cumulus clouds to calculateλbased on the improved bulk-plume method and analyzed the spatial and temporal distributions ofλ,including variations of the PDF ofλwith time and height and its influencing factors.The main conclusions are as follows:

First,an improved method for calculatingλwas developed and its accuracy was validated.The improved method,which solves the conservation equations of two variables simultaneously,is an improvement over the traditional bulkplume method.The calculatedλ,based on the improved method,numerically falls within the range ofλvalues calculated by the traditional bulk-plume method using two different conserved variables,indicating the reliability of the improved method.

Second,the spatial and temporal distributions ofλwere examined.During the entire life cycle of the cumulus ensemble,λinitially decreased and then increased;that is,λin the formation and dissipation stages was larger than that during the maintenance stage.In terms of its vertical distribution,λgenerally exhibited an initial decline and then an increase with increasing height.Regardless of its overall characteristics or variations with time and height,the PDF ofλwas well fitted by a log-normal distribution,demonstrating that it is reasonable to treatλas a stochastic process in numerical simulations.This study provides the PDF ofλto facilitate the randomization of the entrainment process in models.

Finally,the main factors affecting the spatial and temporal distributions ofλwere determined.For the variation of PDF with time and height,λwas negatively correlated with the vertical velocity (w) and positively correlated with environmental relative humidity (RHe),consistent with previous studies.There is a poor correlation betweenλand buoyancy(B) with respect to time,but the relationship betweenλandBwas generally negative as a function of vertical height.Overall,the relationship betweenλandwwas the strongest.These findings offer insight and a new point of reference for handling stochastic entrainment processes in a cumulus parameterization scheme (i.e.,assuming a log-normal PDF).

Two points are noteworthy.First,in future parameterizations ofλ,it is recommended thatwis used as the primary factor.In addition,the influences of RHe,B,and other factors cannot be ignored.It is necessary to comprehensively analyze the effects of these factors onλby combining more observations and high-resolution simulations so that the model convection scheme can describe entrainment more accurately.Moreover,as pointed out by Böing et al.(2012),the detrainment rate is even more important thanλin determining the vertical distribution of convective cloud mass.Thus,accurately obtaining the detrainment rate in cumulus clouds is fairly important to improve cumulus parameterization schemes.

Second,the improved method proposed here can be applied to aircraft observations.Because it is difficult to obtain multiple observations of the same shallow cumulus cloud by aircraft,there is usually only one horizontal penetration for each cumulus cloud.For aircraft observations,the detection height of the aircraft and the height of the cloud base can be taken as two height levels to obtain a large number of observationalλvalues of individual cumulus clouds.To increase the accuracy of the measurement ofλin future aircraft observations,it is highly recommended to measure multiple layers of the same cumulus cloud and the interval between each layer should be as small as feasibly possible.Additionally,the improved method can also be applied to remote sensing data such as enhanced ground-based radar observations,which can obtain the vertical distribution ofλin a large horizontal range.

Acknowledgements.The authors would like to thank Jianfeng GU at Nanjing University for the helpful discussions.This research is supported by the National Natural Science Foundation of China (Grant Nos.42175099,42027804,42075073) and the Innovative Project of Postgraduates in Jiangsu Province in 2023 (Grant No.KYCX23_1319).Shi LUO is supported by the National Natural Science Foundation of China (Grant No.42205080),the Natural Science Foundation of Sichuan (Grant No.2023YFS0442),and the Research Fund of Civil Aviation Flight University of China (Grant No.J2022-037).This research is also supported by the National Key Scientific and Technological Infrastructure project “Earth System Science Numerical Simulator Facility”(EarthLab).The numerical calculations in this paper were conducted on the supercomputing system in the Supercomputing Center of Nanjing University of Information Science &Technology.

Data availabilityThe large-scale forcing data used in this study can be downloaded from the Atmospheric Radiation Measurement (ARM) user facility (https://doi.org/10.5439/1647300 (Tao and Xie,2004);https://doi.org/10.5439/1647174 (Tao and Xie,2012)).LASSO data can be downloaded from https://doi.org/10.5439/1342961 (Gustafson et al.,2017).

Electronic supplementary material:Supplementary material is available in the online version of this article at https://doi.org/10.1007/s00376-023-2357-6.