## Abstract

The entrainment and detrainment rates are important quantities characterizing airmass exchange between clouds and the environment. One of the challenges in calculating the rates is the need to know the velocity vector of the cloud interface in relation to that of the cloud-free air; however, the interface is not well resolved in most cloud model simulations and so the precise value of the vector is not known. Here a new method is described to approximately calculate mass fluxes across the cloud surface in well-resolved simulations of cumulus convection. The method does away with the need to calculate a cloud interface velocity and instead uses gradients of a defined cloud scalar across the cloud interface. As a result, the entrainment and detrainment rates are expressed as an integration over a small region around the cloud interface. The integrand is composed of the total derivative of the cloud scalar. The new method is applied to large-eddy simulations (LES) of a shallow cumulus case and a deep convection case. Compared to a previous method, the approach described here gives 1.5–2 times smaller exchange rates and shows less noise. The smaller exchange rates are explained as the result of differences in how the two methods correct for the advective contribution to variations of cloud volume. Derived two-dimensional distributions of the exchange rates agree well for both methods. Spatial correlation coefficients are about 0.69–0.88 for entrainment and 0.55–0.78 for detrainment.

## 1. Introduction

Cumulus convection is an important process in the atmosphere system because of its transport of heat, moisture, momentum and pollutants, as well as its influences on solar and terrestrial radiation. Due to these important effects it is necessary to correctly represent cumulus processes in atmospheric models. In mesoscale and large-scale climate/weather models a parameterization approach is applied to represent cumulus convection effects. Uncertainties in the parameterization of cumulus convection result in large uncertainties in the estimation of climate sensitivity and the cloud radiative response to changing climate (Bony and Dufresne 2005; Williams and Webb 2009).

The most common scheme for parameterization of convection is the mass-flux method in which cumulus convection is simplified as one or more updrafts. Thermodynamic properties of the updraft change during its ascending process through exchanges between the updraft and its environment. The rates at which air is mixed in (entrainment) and lost (detrainment) are key factors that determine the skill of the parameterization scheme. Several methods exist for parameterizing the entrainment and detrainment rates in atmospheric models (e.g., fixed values or inversely proportional to the updraft radius; Turner 1963); empirical expressions including factors like the environmental humidity, cloud depth, and the updraft buoyancy excess (Bechtold et al. 2008; de Rooy and Siebesma 2008); and the more physics-based buoyancy sorting approach (Kain and Fritsch 1990; Kain 2004). However, all of these methods have deficiencies and further efforts are required to investigate the entrainment and detrainment processes (de Rooy et al. 2013). The cumulus convection processes are usually investigated by observations and cloud-resolving model (CRM) or large-eddy simulation (LES). As a starting point we need accurate estimates of these lateral mixing rates.

The bulk plume method (Siebesma and Cuijpers 1995) and the direct measure (Siebesma 1998) are two approaches for calculating entrainment and detrainment rates. In the bulk plume method, the atmosphere is separated into a cloudy region that moves upward quickly and a more calm environment. Air properties are uniform within each region and airmass exchange occurs between them. Properties of the air moving upward change with altitude due to the exchange. The entrainment/detrainment rates can be inferred from vertical variations of some conserved quantities, like total water, within the plume. In the direct approach airmass flow across the interface between cloud and environment is estimated. In this case air velocities relative to the interface should be known. The bulk plume method is easy to implement in practice. However, the direct measure is difficult to realize, and since the work of Siebesma (1998) only two algorithms have been described, each several years later (Romps 2010; Dawe and Austin 2011). The reason for the difficulty is that the velocity of the cloud interface is hard to determine in the coarse grid LES model. In Romps (2010) (referred to as the Romps method hereafter) entrainment/detrainment is represented as the sources/sinks of air activity (equal to unity in the cloud and zero otherwise) and then does not require the interface velocity. Dawe and Austin (2011, hereafter DA11) calculate the sources/sinks of air activity by using a complicated interpolation algorithm. Applications to LES of cumulus convection shows the two algorithms become closer to each other as the LES resolution decreases from 100 to 25 m, but the Romps method requires a long time average to obtain robust results (at least a 1-min average for a 1.5- s sampling frequency) (DA11).

This work describes an alternative method for providing a direct measure of entrainment and detrainment in cumulus convection. The method is directly based on the formula given in Siebesma (1998) which uses the interface velocity explicitly. The Romps (2010) method is based on the sources/sinks of the air activity and bypasses the determination of the interface velocity, as does the DA11 method. The Romps method shows about 50%–100% noise (virtually random) over one model time step (e.g., 1 and 2 s) in times series of the estimated entrainment/detrainment rates as revealed in DA11. Averaging over long time (e.g., 1 min) is required to reduce the noise. Compared to the Romps method, the DA11 method gives about 1.5 times lower entrainment/detrainment rates at 25-m grid spacing, about 2 times lower rates at 50-m grid spacing and 4 times lower rates at 100-m grid spacing. So it is not well suited to the LES model with grid spacing larger than 50 m. By developing a different kind of approach for calculating entrainment/detrainment rates, the study here aims to obtain more insight on entraining/detraining processes and look for possibility of overcoming the deficiencies encountered by existing methods. In section 2 the alternative method is derived and comparisons to the existing methods are discussed. Applications to LES models are presented in section 3. Section 4 presents discussion and in section 5 conclusions are drawn.

## 2. Method

### a. Method description

According to Siebesma (1998) the entrainment (*E*) and detrainment (*D*) rates are defined as

where *A* is the area of the model domain in the horizontal (m^{2}), *ρ* is the air density (kg m^{−3}), **n** is the unit cloud interface vector pointing outside, and **u** and **u**_{i} are the three-dimensional airflow velocity and cloud interface velocity (m s^{−1}). The integration is along the cloud boundary that is an intersection between the cloud interface and a horizontal plane. *E* and *D*, in kg m^{−3} s^{−1}, are measures of the airmass flux flowing into and out from the cloud across its interface. The cloud interface velocity has to be determined in order to apply Eq. (1) directly.

Consider a three-dimensional field of some scalar *ξ*(*x*, *y*, *z*, *t*) at some time *t*. The moving velocity of isosurfaces of this scalar in the direction perpendicular to the isosurfaces can be written as

Then the airmass flux across the isosurface is

Both sides of Eq. (3) are evaluated at the isosurface as well.

If we define some cloud scalar *ξ* with its value larger than a critical value *ξ*_{c} as the definition of the cloud, then the interface velocity in Eq. (1) can be explicitly evaluated with Eq. (2). Movement of the cloud interface is caused by advection associated with airflow and microphysical (e.g., condensation/evaporation, conversion into precipitation) and diffusive processes. Accordingly, the local tendency of the cloud scalar in Eq. (2) includes contributions by advection, microphysics and diffusion. The cloud interface movement that can be attributed to processes other than advection is referred as the entraining/detraining velocity.

Accordingly, the entraining/detraining airmass flux can be estimated using Eq. (3). We obtain an alternative expression for the integrand in Eq. (1). In a comparison with Eq. (2) the total derivative *dξ*/*dt* replaces the local tendency of the cloud scalar. It is well known that the total derivative represents variabilities of a quantity after correcting for advective transport and referred as source/sink. This is consistent with the fact that entrainment/detrainment (airflow across the cloud interface) arises only from physical processes instead of advection. In a case where the cloud is transported by airflow only, there is not a source/sink for the airmass enclosed by the cloud interface and therefore no entrainment/detrainment occurs. In cumulus convection studies, the cloud is usually defined as upwardly moving cloudy air parcels. Therefore all physical processes that influence the vertical velocity and cloud condensates can contribute to entrainment/detrainment, such as buoyancy, precipitation, condensation/evaporation and diffusion. Diffusion strength depends on how the resolved motion is defined. In the real atmosphere, air motion occurs at spatial scales down to the Kolmogorov scale (~1 mm), below which the effects of viscosity and diffusion are important. In a discrete numerical model, the resolved motion is air motion resolvable by grid spacing and diffusion is parameterized as subgrid turbulent mixing. The two processes are numerically simulated following different procedures. Contribution by the resolved motion in a numerical model is completely removed by using the total derivative. The parameterized subgrid turbulent mixing contributes to entraining/detraining processes.

Replacing the integrand with the right side of Eq. (3) and taking *E* as an example, we obtain

The term **n** ⋅ ∇*ξ* represents the cloud scalar gradient along the interface direction and can be rewritten as ∂*ξ*/∂*l*_{i}; *l*_{i} is a distance perpendicular to the interface. Thus Eq. (4) can be written as

where *δξ* = (∂*ξ*/∂*l*_{i})*dl*_{i} is the infinitesimal variation of *ξ* within the infinitesimal distance *dl*_{i} along the interface direction. If the projection of *dl*_{i} on the horizontal plane along the interface is *dl*_{x} and the angle between **n** and the horizontal plane is *β* (see Fig. 1), then there is

where *dS* is the area of an infinitesimal region bounded by isolines of *ξ*_{c} and *ξ*_{c} + *δξ* with a length of *dl* along the isolines of the cloud scalar on the horizontal plane (see Fig. 1). Inserting Eq. (6) into Eq. (5) leads to

where the integration is performed in the region bounded by the isolines of *ξ*_{c} and *ξ*_{c} + *δξ*. From a numerical viewpoint areal integration is much easier to implement than line integration and this is why Eq. (4) is converted into Eq. (7). In practice a small value of *δξ* is used and since the cloud is defined as the region with *ξ* > *ξ*_{c}, the quantity *dS* has a negative sign. The cloud interface direction is the opposite of the gradient of *ξ*, so cos*β* is calculated as

Similarly, the detrainment *D* can be expressed as

In Eqs. (7) and (9) the integrations are performed in a small region at the cloud boundary, which requires that the cloud and its boundary are resolved well. However, the model grid spacing does not always satisfy this requirement. For example, in shallow cumulus the diameter of the cloud plume can be as small as a few tens of meters. In this case a model grid spacing of 25 m (as frequently applied grid spacing in shallow cumulus simulations) cannot resolve all the plumes reasonably. In order for the proposed method to be applied in the coarse-grid model, the cloud scalar and its total derivative must be remapped onto finer grids. In this study an interpolation method as follows is applied. Assuming we have an interpolation operator *F*, then $\xi i=F\u2061(\xi )$, where the superscript *i* denotes the quantity after applying interpolation and *ξ* is some scalar. The total derivative of the interpolated scalar can be estimated as

where Δ*t* is the model time step. The interpolation operator applied here is the piecewise parabolic (PPM) interpolation with monotonicity constraints (Lin and Rood 1996). The PPM method conserves mass and sharp gradients and is well suited to the cloud scalar that usually shows a significant gradient around the cloud boundary. The one-dimensional PPM interpolation algorithm is applied sequentially to *x*, *y*, and *z* coordinates to obtain model variables and their total derivatives on finer grids.

In previous studies the cloud updraft is defined as the region with positive cloud mixing ratios and positive vertical velocities, and additionally with the cloud core having positive buoyancy. To be consistent with the previous definition, the cloud scalar applied here is a product of the cloud mixing ratio *r*_{c} [in kg (cloud) per kg (dry air)] and the vertical velocity *w* (in m s^{−1}). When mapped onto finer grids the cloud scalar and its total derivative can be written as

### b. Comparison with previous methods

Previous methods for directly calculating the *E* and *D* include the method described in Romps (2010) and the algorithm developed in DA11. In the Romps method a quantity *A* termed as the activity of the air is defined as unity within the cloud and zero otherwise. The activity source ∂*ρA*/∂*t* + ∇ ⋅ *ρ***u***A* is evaluated with the first order upwind algorithm and summed during the whole period when a grid cell is adjacent to the boundary between the cloud and the environment. Once the grid cell is away from the boundary, the summation of the activity source is attributed to *E* if it is positive and *D* if negative. The meaning of the procedure is straightforward because when the activity source is positive, the cloud extends into a grid cell after correcting for the contribution by advection. There is some similarity between the air activity *A* and the cloud scalar defined here. The two methods differ in the ways of removing the contribution by advection. The method described here removes the advective contribution by applying the modeled total derivative of the cloud scalar, in contrast to the Romps method of calculating the total derivative of *A* itself. In DA11,*E* and *D* are estimated based on the source of a cloud volume enclosed by the grid cell walls lying in the cloud and a cloud surface intersecting with the grid cell. The cloud volume source is $\rho \u2061(dV/dt)+\u222bW\rho u\u22c5dW$, where *V* is the cloud volume contained in the grid cell, **W** is the portion of the grid cell walls within the cloud, and the air density is assumed to be time independent. As in the Romps method, the expression for the source includes two parts; the first term is the total tendency of the volume and the second is minus of the volume variation due to advection. The source is attributed to *E* if it is positive and *D* if negative. This cloud volume source is basically an integral expression of the activity source introduced in Romps (2010) except for an interpolation is applied in Dawe and Austin (2011).

The Romps and DA11 methods can handle entrainment/detrainment along both lateral and vertical cloud interfaces. But the described method here includes those along the lateral interface only and does not work when the cloud interface is nearly horizontal (e.g., at the bottom and top of cylindrically shaped interface). This is because that Eq. (1) as a starting point of the described method does not account for interfaces with vertical normals. In summary, the new method mainly differs from the previous ones in the way it corrects for the advective contribution to the change in the cloud volume.

## 3. Application to LES of cumulus convection

In this section the method developed in section 2 is applied to large-eddy simulations of both a nonprecipitating shallow cumulus convection case and a precipitating deep convection case. The shallow cumulus case was observed during the Barbados Oceanographic and Meteorological Experiment (BOMEX; Holland and Rasmusson 1973). This is a trade wind cumulus case with steady behavior and without apparent precipitation. Initial profiles and forcings including the subsidence and prescribed cooling and drying representing large-scale processes are detailed in Siebesma et al. (2003). The LES was performed with the anelastic nonhydrostatic mesoscale model Meso-NH (Lafore et al. 1998). The turbulence parameterization was a 1.5-order closure (Cuxart et al. 2000). The microphysical scheme was based on Pinty and Jabouille (1998). Fluxes at the surface were dynamically calculated in the model. The rain processes were excluded in the simulation. The simulation was conducted in a domain of 5 × 5 km^{2} in the horizontal and 3 km in the vertical. The model grid spacing was set to 25 m (or 50 m) in three directions. Periodic lateral boundary conditions weer used and a sponge layer was specified from 2.5 to 3 km. The model was run for 5.52 h and the outputs during the last 16 min were retained for every time step (1 s for 25-m grid spacing and 2 s for 50-m grid spacing) for later analysis.

The selected deep convection case was a 100-km-long squall line observed during the Tropical Ocean and Global Atmosphere Coupled Ocean–Atmosphere Response Experiment (TOGA COARE) on 22 February 1993. The initial sounding used by the model is from Redelsperger et al. (2000). The turbulence parameterization, microphysical scheme and surface fluxes calculation were the same as in the shallow cumulus case except that precipitation was allowed. The simulation was performed in a domain of 24 × 24 km^{2} in the horizontal and 32 km in the vertical. The model grid spacing was 100 m (or 200 m) in the three directions. Similarly, Periodic lateral boundary conditions were used and a sponge layer was specified above 22 km. Instead of applying large-scale forcings, nudging to the initial conditions of potential temperature, water vapor and horizontal winds was applied with a nudging time scale of 16 h. The model was run for 22.25 h for the model with 100-m grid spacing (19.25 h for the model with 200-m grid spacing) and the outputs during the last 15 min were retained for every time step (2 s for the 100-m grid spacing and 4 s for the 200-m grid spacing). Different times were used to keep outputs because cumulus convection is completely different between models with 100- and 200-m grid spacing. Around 22 h there was no convective activity in the model with 200-m grid spacing.

As described in section 2, performing an integration in a small region at the cloud boundary is required to calculate the *E* and *D.* First the cloud scalar and its total derivative are interpolated to 3 times refined grids for 25- and 100-m grid spacing (5 times refined for 50- and 200-m grid spacing). The cloud boundary is defined as *ξ*_{c} = 5.0 × 10^{−6} m kg s^{−1} kg^{−1}; *δξ* at each altitude level is determined by decreasing a large value continuously until either of the two following conditions is broken. One condition is that area bounded by the contour lines *ξ*_{c} and *ξ*_{c} + *δξ* is not smaller than 0.25 times the cloud area. The other is that the number of the refined grid points located in the bounded region is not smaller than 1.3 times the number of the points just adjacent to the environment. These conditions were applied to make sure enough points were included and those far away from the cloud boundary were excluded. One example is shown in Fig. 2. As a comparison the *E* and *D* were calculated as well with the method developed in Romps (2010) as described in section 2b.

### a. Shallow cumulus

The entrainment and detrainment rates are given in Fig. 3 for the shallow cumulus case simulated by the model with 25- and 50-m grid spacing. These rates are calculated based on Eqs. (7) and (9) for the method described in this study (represent horizontal averages by definition), and are horizontally averaged values for the Romps method. The estimated *E* and *D* should agree with the cloud continuity equation, ∂*ρa*/∂*t* + ∂*M*/∂*z* = *E* − *D*, where *a* is the cloud fraction and $M=\u222b\xi >\xi c\rho w\u2009dx\u2009dy/A$ is the fractional vertical mass flux of the cumulus convection. Both sides of the continuity equation are shown in Fig. 3. For the shallow cumulus case the cloud amount (the cloud fraction and the updraft mass flux) grows below about 600 m and continuously declines above this level (not shown). The dominance of detrainment over entrainment is a common characteristic for the shallow cumulus convection. The *E* and *D* estimated with the method developed in section 2 are consistent with this picture; *D* is larger than *E* at most altitudes. The airmass exchange between the cloud and the environment occurs through turbulence induced mixing processes, condensation/evaporation of cloud water, buoyant forcing, and pressure gradient force-induced vertical accelerations.

Compared to the *E* and *D* calculated with the Romps method, the approach described here does not agree as well with the cloud continuity equation and gives lower magnitudes for the exchange rates. As has been stated the new method requires that each convective plume is well resolved and so *E* and *D* are calculated with the interpolated cloud scalar and its total derivative. The interpolated variables could not be dynamically consistent since the modeled cloud evolves on the original model grids. On the other hand, it is not surprising that the Romps method agrees well with the cloud continuity equation. Performing an integration for the activity source over the *i*th model layer, $\u222b\u2009dx\u2009dy\u222bzi\u22121/2zi+1/2\u2061(\u2202\rho A/\u2202t+\u2207\u22c5\rho uA)\u2009dz$, we arrive at the cloud continuity equation naturally. So the Romps method agrees with the cloud continuity equation by design. The amplitudes of exchange rates from the Romps method are about 1.5 times greater than those from the method here.

### b. Precipitating deep convection

The entrainment and detrainment rates for the precipitating deep convection case simulated with 100- and 200-m grid spacing are displayed in Fig. 4. These rates are horizontally averaged values as in section 3a. As in the shallow cumulus case the *E* and *D* calculated with the method described here show reasonable agreement with the cloud continuity equation and the Romps method achieves a better agreement. One exception occurs at altitudes below around 3000 m, where *E* − *D* from the method here shows consistent negative biases compared to the cloud continuity equation. The exchange rates are 2 times lower than those estimated by the Romps method. It should be noted that the model with 100-m grid spacing simulates a different dynamics than the one with 200-m grid spacing. However, the Romps method gives largely different *E* and *D* profiles between the 100- and 200-m grid spacing models. In a comparison the difference is much smaller for the method described. This indicates less sensitivity to grid spacing using the new method compared to the Romps method.

### c. Temporal stability

To investigate the temporal resolving ability of the new method, Fig. 5 presents time series of the *E* and *D* at 1000 m for the shallow cumulus case modeled by the 25-m grid spacing model. The results show that the new method is more stable than the Romps method. With a temporal resolution of 1 s the estimated *E* and *D* show some noise. After applying a moving average to the time series with a window width of 3 s, the estimated *E* and *D* show much improved stability. The results here are similar to those presented in DA11 where the Romps method showed much higher noise than their method. As explained in section 2b, the method described in DA11 is an integral version of the Romps method except that an interpolation is applied. Without the interpolation, the cloud volume can only change by integral multiples of the grid cell volume per model time step. Fractional changes are allowed when the interpolation is applied, which reduce the noise encountered in the Romps method.

### d. Spatial variability

Similar to the Romps and DA11 methods, the method described here is able to resolve spatial variation of the exchange rates. Figure 6 shows the two-dimensional distributions (height versus, horizontal distance) of the *E* and *D* for the shallow cumulus simulated by the model with 25 m grid spacing. There are similarities between the spatial distributions of the exchanging rates from the Romps method and the method used here. Correlation coefficients are 0.77/0.55 (coefficients of *E*/*D*), 0.88/0.70, 0.70/0.78 and 0.69/0.78 for cumulus convection simulated by the models of grid spacing 25, 50, 100 and 200 m. *E* and *D* used for calculating the correlation coefficients are averaged over 1, 2, 2, and 4 min, respectively. Different averaging times are used because these models have different model time steps.

## 4. Discussion

Above we have described a new method to directly calculate entrainment and detrainment rates, applied it to case studies of shallow cumulus and deep convection events, and compared results with the previous method of Romps (2010). Here we discuss some of the underlying factors that affect the new method and can help explain differences with prior approaches.

### a. Influences of the choice of variables to define cloud scalar

As discussed in section 2a, the entraining/detraining rates expressed by Eqs. (7) and (9) depend on the sources and sinks of the cloud scalar as encapsulated by *dξ*/*dt*. As a result, the particular choice of variables to define *ξ* has important influences on estimated rates. The influences are investigated in this section to explain characteristics of the described method and reveal dependence of the exchange rates on how the cloud is defined. In section 3, a definition *ξ* = *wr*_{c} was applied, where *r*_{c} is the cloud mixing ratio, and the sources/sinks include buoyancy, precipitation, condensation/evaporation and subgrid mixing. Here we test another definition *ξ* = *wr*_{t}, where *r*_{t} is total water mixing ratio. In the revised definition the condensation/evaporation no longer contribute to entrainment/detrainment processes. Figure 7 shows exchange rates calculated based on these two definitions in the shallow cumulus case simulated by the model with 25-m grid spacing. For total water updraft, *ξ*_{c} = 8.0 × 10^{−3} m kg s^{−1} kg^{−1} is selected to make sure that cloud fractions are similar between the two definitions. The results are an average over a 1-min period. As expected the total water updraft definition presents generally smaller exchange rates than the cloud updraft definition for most altitudes. These differences are consistent with different response of total water and cloud condensates to vertical motion.

### b. Differences between the Romps method and the new method

It is shown in section 3 that the new method gives rates that are about 1.5–2 times smaller than the Romps method. One possible explanation is that the differences come from different approaches to correcting for the advective contribution to the cloud interface velocity. The Romps method represents the source of the cloud volume as a residual between the local tendency and advection-induced variation of the air activity. Then it has to specify a transport scheme to estimate the advective transport of the air activity. The transport scheme used for calculating advective transport of the air activity has important impacts on results. In DA11 two algorithms are applied to estimate the advective contribution, the first order upwind approach as proposed in Romps (2010) and the transport scheme of the LES model they used. The latter algorithm leads to lower *E* and *D* compared to the first order upwind algorithm. However, the method described here uses the total derivative of the cloud scalar which is an output of the model and does not include the advective contribution at all. The Romps method could be used with any transport scheme since the cloud activity is not a model-simulated scalar. One could expect that a scheme consistent with numerics of the LES model (PPM scheme for the model applied here) is more suitable to the Romps method. Figure 8 presents exchange rates calculated with the new and the Romps method using the first order upwind and PPM transport schemes for the shallow cumulus case simulated by the model with 25-m grid spacing. The results indicate that applying the PPM scheme gives exchange rates that are lower than applying the first order upwind scheme but still 1.5–2 times higher than those from the new method. This is because air activity is not a modeled scalar and even applying a numerically consistent transport scheme does not ensure a numerically consistent advective contribution to the interface velocity. A way to achieve complete numerical consistency is modifying the sources and sinks of air activity in the Romps method as ∂*ρA*/∂*t* − (∂*ρA*/∂*t*)_{adv}, where

*H* is a function of *ξ* which is unity when *ξ* ≥ *ξ*_{c} and zero otherwise. −**u** ⋅ ∇*ξ* represents advective source of the cloud scalar. This source can be derived from outputs by the LES model since the cloud scalar consists of modeled quantities. So in the modified Romps method specifying a transport scheme for calculating advection of *A* is not required anymore. The advective contribution to the cloud volume is estimated in a way consistent with the model numerics and then the new method. However, the modified Romps method cannot be applied to the cloud definition used in section 3. In this definition the advective source of the cloud scalar could be estimated as

or

where the subscript adv means advective source outputted by the model. In the model, advection and subgrid turbulent mixing are calculated at first for cloud condensates and other thermodynamic quantities. Then a saturation adjustment is applied to ensure a consistency of the model state with thermodynamic laws. The cloud state that solely undergoes advection, *ξ* − **u** ⋅ ∇*ξ*Δ*t*, could not exist. This is because the state *r*_{c} + *r*_{c,adv}Δ*t* can be inconsistent with the thermodynamic laws and is not a physical state. To test the modified Romps method, a pollutant-liked tracer with a lifetime of 1 h and unit mixing ratio at the surface is transported in the model. The cloud scalar is taken as the tracer mixing ratio and clouds are defined as regions where the mixing ratio is larger than 0.2 kg kg^{−1}. The advective source of the cloud scalar, −**u** ⋅ ∇*ξ*, is just taken from outputs by the LES model. In Fig. 9, exchange rates calculated with the new method, and the Romps method and the modified Romps method are presented for the shallow cumulus case simulated by the model with 25-m grid spacing. Results show that the modified Romps method gives the same entraining/detraining rates as the new method. It is reasonable to speculate that differences between the Romps method and the new method are caused by their methods of correcting for the advective contribution to the interface velocity, even though an explicit demonstration given above is achieved only for the case where the cloud is defined by the tracer plume.

### c. Influences of δξ

As described in section 3, the small value *δξ* was determined by making sure enough number of the refined grid points are bounded by *ξ*_{c} and *ξ*_{c} + *δξ*. Strictly speaking, Eqs. (7) and (9) make sense in the limit where *δξ* approaches zero. However, as *δξ* becomes small, the number of times the interpolation is refined should increase to make sure there are enough grid points included by isolines *ξ*_{c} and *ξ*_{c} + *δξ*; as a result, the numerical calculation becomes difficult. Figure 10 shows the dependence of entrainment and detrainment rates on values of *δξ* for the shallow cumulus case simulated by the model with 25-m grid spacing. The values of *δξ* are determined by making sure that the number of the bounded refined grid points is about 1.3 times that of the grid points just adjacent to the environment; the sum of the area of the bounded refined grid cell is about 0.25, 0.10, and 0.05 times the area of the cloud. An area fraction limit of 0.25 is chosen and applied in section 3 for easing numerical work and is sufficient for illustrating the methodology.

## 5. Conclusions

This study describes a method to directly calculate the entrainment and detrainment rates in cumulus convection. The method represents the velocity of the cloud interface explicitly through use of a cloud scalar, and the airmass flux across the cloud interface is then calculated directly. The method is applied to a shallow cumulus case simulated with 25- and 50-m grid spacing models and a deep convection case with 100- and 200-m grid spacing models. Compared to the method of Romps (2010), the new method agrees less accurately with the cloud continuity equation and gives 1.5–2 times smaller exchange rates. The reason could be that the two methods use different ways of correcting for advective contributions to cloud volume variabilities. The new method shows less noise than the Romps method. An examination of the spatial distribution of exchange rates calculated with the two methods reveals a high similarity between them.

In the previous work (Siebesma and Cuijpers 1995; Romps 2010) a combination of several conditions is used to recognize the cloud region, such as the cloud core (positive buoyancy relative to the horizontal mean in addition to the cloud updraft conditions). However, the cloud core defined in these studies cannot be applied in this study. It is not straightforward to define a cloud scalar that recognizes fluid volumes with a condensate mixing ratio larger than 1.0 × 10^{−5} kg kg^{−1}, a positive vertical velocity and positive buoyancy. One distinct feature of the new method is that the exchange rates are represented based on the total derivative of the cloud scalar. This could make it possible to separate the exchanging rates into contributions by different physical processes as conceptually proposed by de Rooy et al. (2013). The problem of extending the method to include the cloud core definition and resolve different physical contributions is a topic for future research.

## Acknowledgments

Thanks are due to the editor and anonymous reviewers for suggestions that improved the manuscript. This work is supported by the Fundamental Research Funds for the Central Universities lzujbky-2018-5.

## REFERENCES

*Quart. J. Roy. Meteor. Soc.*

*Geophys. Res. Lett.*

*Quart. J. Roy. Meteor. Soc.*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*Quart. J. Roy. Meteor. Soc.*

*Mon. Wea. Rev.*

*J. Appl. Meteor.*

*J. Atmos. Sci.*

*Ann. Geophys.*

*Mon. Wea. Rev.*

*Quart. J. Roy. Meteor. Soc.*

*J. Atmos. Sci.*

*Buoyant Convection in Geophysical Flows*, E. J. Plate, Ed., Kluwer Academic Publishers, 441–486

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Fluid Mech.*

*Climate Dyn.*

Conf. on Cloud Physics, Everett, WA, Amer. Meteor. Soc., 217–220