## Abstract

Lagrangian particle tracking is used in a large-eddy simulation to study an individual cumulus congestus. This allows for the direct measurement of the convective entrainment rate and of the residence times of entrained parcels within the cloud. The entrainment rate obtained by Lagrangian direct measurement is found to be higher than that obtained using the recently introduced method of Eulerian direct measurement. This discrepancy is explained by the fast recirculation of air in and out of cloudy updrafts, which Eulerian direct measurement is unable to resolve. By filtering these fast recirculations, the Lagrangian calculation produces a result in very good agreement with the Eulerian calculation.

The Lagrangian method can also quantify some aspects of entrainment that cannot be probed with Eulerian methods. For instance, it is found that more than half of the air that is entrained by the cloud during its lifetime is air that was previously detrained by the cloud. Nevertheless, the cloud is highly diluted by entrained air: for cloudy air above 2 km, its mean height of origin is well above the cloud base. This paints a picture of a cloud that rapidly entrains both environmental air and its own detritus.

## 1. Introduction

A recent large-eddy simulation (LES) of deep atmospheric convection has revealed fractional entrainment rates at the astonishingly high rate of about 1 km^{−1} throughout the depth of the troposphere (Romps 2010). Those rates were obtained using a technique referred to as “direct measurement” of entrainment. Although the method has been validated in the case of shallow cumulus (Dawe and Austin 2011a), no such independent validation has yet been made for deeper convection. Here, we validate the previous results in an LES of a single cumulus congestus using a new and independent method: direct measurement of entrainment with massless Lagrangian particles.

Recent evidence has also indicated that convective clouds entrain a substantial amount of air whose equivalent potential temperature *θ _{e}* is substantially larger than in the clear-air environment. This was argued to contribute to the discrepancy between entrainment rates obtained by direct measurement and the bulk-plume equations (Romps 2010). Such a discrepancy was demonstrated for the case of total water and vertical velocity by Dawe and Austin (2011b). Here, we use the Lagrangian particles to explore the difference in

*θ*between entraining air and the surrounding environment.

_{e}The paper is organized as follows. The numerical methods and the simulation are briefly summarized in sections 2 and 3, respectively. Section 4 describes the details of the Lagrangian direct measurement of entrainment. The in-cloud residence times of entrained air parcels are studied in section 5. In section 6, we discuss the height of origin for entrained air parcels and their thermodynamic properties. Finally, section 7 gives some concluding remarks.

## 2. Numerical methods

The LES is performed using the fully compressible, nonhydrostatic, cloud-resolving Das Atmosphärische Modell (DAM) (Romps 2008). The governing equations are discretized using a fifth-order weighted essentially nonoscillatory (WENO) method (Shu 1997), and time is advanced using a total-variation-diminishing third-order Runge–Kutta scheme (TVD-RK3) (Shu and Osher 1988). Microphysics is modeled by the six-class Lin–Lord–Krueger scheme (Lin et al. 1983; Lord et al. 1984; Krueger et al. 1995a). In the present study, a subgrid-scale model is not included because it was previously found to lead to poor distributions of lower-tropospheric cloud (Romps and Kuang 2010a) and violations of the second law of thermodynamics (Romps 2008).

For the Lagrangian particle tracking, the velocity **V**(*t*) and thermodynamic variables at the particle location **X**(*t*) are computed by an interpolation from the Eulerian data. This approach has been used previously in LES to study the convective boundary layer (Weil et al. 2004; Kim et al. 2005; Dosio et al. 2005; Gopalakrishnan and Avissar 2000), stratocumulus and cumulus (Krueger et al. 1995b; Heus et al. 2008; Yamaguchi and Randall 2012), and congestus and cumulonimbus (Lin and Arakawa 1997a,b; Carpenter et al. 1998; Böing et al. 2012). We have tested several interpolation schemes. It is found that a simple linear interpolation, as used by Lin and Arakawa (1997a,b), is not sufficient to resolve the motion of a Lagrangian particle near a cloud boundary, where the spatial gradient of an Eulerian field becomes very large. To account for such a discontinuity in the numerical solution, the Eulerian fields are mapped onto the Lagrangian particles using a third-order WENO approximation. Consistent with the LES solver, the equation of particle motion, *d***X**/*dt* = **V**, is integrated with TVD-RK3. Like the thermodynamic fields, the velocity **V** is obtained by interpolation of the resolved Eulerian velocity to the location of the particle. To ensure that the particles are tracking the resolved flow, this velocity is not supplemented with any subgrid-scale turbulent velocity. This same choice was also made by Gopalakrishnan and Avissar (2000), Dosio et al. (2005), and Böing et al. (2012).

## 3. Simulation

In the present study, convective entrainment is investigated by analyzing the Lagrangian statistics of air parcels in a three-dimensional LES of a single cumulus congestus. To initiate this cumulus congestus, a Gaussian-distributed temperature disturbance is added to a background state. The background temperature profile *T*_{0}(*z*) and specific humidity profile *q _{υ}*

_{0}(

*z*) are taken from an LES of radiative–convective equilibrium over a 300-K ocean. The initial temperature field is then given by

where Δ*T* = 1 K is the magnitude of the temperature perturbation; |**r**| is the radial distance from the center of the computational domain; *r*_{ref} = 1000 m and *z*_{ref} = 500 m are the length scales of the warm thermal in the horizontal (*x*–*y*) and vertical (*z*) directions, respectively; and *ξ* is a random variable uniformly distributed in (−1, 1). The specific humidity is then increased to match the relative humidity at that height before the temperature disturbance was added,

where an asterisk denotes the saturation humidity. The initial density is then specified to give hydrostatic balance in each grid column.

The computational domain has dimensions of (*L _{x}* ×

*L*) = (20 km × 20 km) in the horizontal plane and 27.5 km in the vertical direction. The computational grid spacing is (

_{y}*δx*,

*δy*,

*δz*) = (50 m, 50 m, 50 m). The time-step size of

*δt*= 1 or 2 s is specified adaptively based on the Courant–Friedrichs–Lewy (CFL) condition. Periodic boundary conditions are used in the horizontal directions and the lower boundary is no-slip. Since this is a transient simulation of a single cloud, neither surface fluxes nor radiation is used. For the Lagrangian particles, a bounce-back boundary condition is imposed at

*z*

_{low}= 0 and

*z*

_{upp}= 20 km.

Figure 1 shows the three-dimensional cloud boundary at 15, 20, 25, and 30 min. The cloud boundary in this figure is defined by the 10^{−5} kg kg^{−1} isosurface of the mass fraction of nonprecipitating cloud condensates *q _{c}*. The base of the cloud is located at roughly

*z*= 500 m, and the cloud top reaches just over 8 km. The radius of the cloud ranges from about 500 m at the cloud base to 1000 m near the cloud top. Figure 2 shows the trajectories of five particles over

*t*= 0–1 h. These five trajectories are chosen among the particles whose location is below the cloud base at

*t*= 0 and around the cloud top at

*t*= 1 h. The colors on the trajectories denote

*q*. The particle trajectories exhibit rotational motions because of turbulent eddies, and the changes in the color contours indicate the high spatial variability of

_{c}*q*within the cloud. After detraining, the particles travel a few kilometers in the horizontal directions due to divergence.

_{c}Here, we use 4 × 10^{7} Lagrangian particles as a conjugate solution of the Eulerian LES. As such, a subgrid-scale fluctuation is not added to the equation of motion of the Lagrangian particles. This allows the Lagrangian particles to follow the resolved flow. In this way, changes in a particle’s conserved variables, such as equivalent potential temperature, are attributable to unresolved mixing processes, such as subgrid or numerical diffusion. This choice also allows the massless particles to be associated with a constant effective mass of dry air, as described below.

Following the volume-averaging procedure of Jackson (1997), the governing equation of the number density of the Lagrangian particles is given by the conservation equation,

in which *d*/*dt* is the total derivative, *n* is the number density of the Lagrangian particles (*m*^{−3}), **u** is the velocity of air, and *ρ _{a}* is the density of dry air. If the fluid is incompressible,

**∇**·

**u**= 0, then Eq. (2) implies that the number density field is invariant in time. In a compressible flow, however, the volume of an air parcel is not conserved and, hence, neither is

*n*. If we define the specific number density (kg

^{−1}) as

*n** =

*n*/

*ρ*(i.e., the number of particles in a unit mass) then it is trivial to show that the specific number density field is conserved in time,

_{a}*dn**/

*dt*= 0.

Figure 3a shows the horizontally averaged number density at each 50-m grid level. For the initial condition, the particles are distributed randomly in the computational domain with a constraint that the horizontally averaged *n** is constant in the vertical direction. As *ρ _{a}* is a decreasing function of height,

*n*is also largest near the ground and decreases with the height. The horizontally averaged

*n** on each grid level is shown in Fig. 3b at the beginning of the simulation (black line) and after 1 h (red circles). The equations for

*n*and

*n** indicate that, if the continuity equation and the particle velocity are computed with a high accuracy,

*n*and

*n** should be invariant in time, respectively, in incompressible and compressible flows. In the present simulation, as expected, the horizontally averaged

*n** remains virtually unchanged in time. Heus et al. (2008) found significant spurious clumping of particles at the lower boundary when the subgrid velocity fluctuation was omitted (see their Fig. A1), but we find that the difference between the smallest and largest

*n** in the first vertical grid cell remains less than 1% throughout the computational period. We have tested a few interpolation schemes and found that using a linear interpolation may result in a larger fluctuation of the specific number density near the lower boundary.

In the present Lagrangian framework, each particle represents a fraction of the total mass of dry air in the computational domain from *z* = 0 to *L _{z}* = 20 km (recall that the bounce-back level is 20 km). The mass per particle is simply

where *N _{p}* = 4 × 10

^{7}is the total number of the particles. The motion of each Lagrangian particle represents the transport of dry air of

*δm*≃ 100 tons, which corresponds to a cubic mass of near-surface air with a width of about 45 m.

It is important to note that the Eulerian thermodynamic fields are subject to numerical diffusion in a way that the Lagrangian particles are not. Recall that the thermodynamic properties associated with a Lagrangian particle are given by the interpolation of the Eulerian fields to the position of the particle. Therefore, as particles rise up along with the cloud, one might imagine that water would diffuse out of the cloud numerically, leading the in-cloud particles to have a smaller mean equivalent potential temperature than the mean value with which they started. Here, the equivalent potential temperature is a conserved variable for all adiabatic and reversible transformations (Romps and Kuang 2010a). To test for this effect, we compute the mean *θ _{e}* of the cloud using two different methods: once using

*θ*of the in-cloud particles at the time of sampling and once using the initial

_{e}*θ*(i.e., at

_{e}*t*= 0) of the in-cloud particles. First, we define an activity operator to identify a cloudy region. In particular, we focus here on cloudy updrafts, as they are responsible for the majority of latent heat release and are the focus of most convective parameterizations. Cloudy updrafts are defined as those regions where is positive. The activity operator used here is

where *w* is the vertical velocity. Following Romps (2010), we use *q*_{threshold} = 10^{−5} kg kg^{−1} and *w*_{threshold} = 1 m s^{−1}. Using this activity operator, we define

in which *T* is the total simulation time (1 h) and *Z _{i}*(

*t*) and are the vertical location and the equivalent potential temperature of the

*i*th particle at time

*t*, respectively. The term

_{I(z)}[

*Z*] is an indicator function of a small interval centered at

*z*,

*I*(

*z*) ≡ [

*z*−

*δz*/2,

*z*+

*δz*/2], which is defined as

_{I(z)}[

*Z*] = 1 if

*Z*∈

*I*(

*z*) and

_{I(z)}= 0, otherwise. The time-averaged cloud mass 〈

*m*〉(

*z*) is computed from

The results are shown in Fig. 4. Below 9 km, the agreement between 〈*θ _{e}*〉 and is satisfactory. The difference between two computations is less than about 0.5 K. It is not straightfoward to pin down the source of this difference. It may stem from numerical diffusion, numerical error involved in the interpolation of Eulerian variables, or from the fact that

*θ*is not exactly conserved in a precipitating cloud. Above 9 km, shows a significant fluctuation, which seems to be related to the small sample size in the region. The total number of particles sampled in each 100-m-wide bin below 9 km over the simulation period of 1 h is about

_{e}*O*(10

^{4}–10

^{5}), while above 9 km the number of particles reduces dramatically to less than 2000. We conclude from this figure that numerical diffusion likely affects the

*θ*of individual particles (as illustrated by the noise above 9 km, where particle numbers are small), but it does not introduce any significant mean bias (as illustrated by the excellent match below 9 km, where particle numbers are large).

_{e}## 4. Entrainment rate

The activity operator, defined in Eq. (3), can be applied to either the *i*th Lagrangian particle [i.e., _{i}(*t*)] or the Eulerian location **x** [i.e., (**x**, *t*)]. Let us consider a thin horizontal slice of the atmosphere defined by the height interval *I*(*z*) ≡ [*z* − *δz*/2, *z* + *δz*/2]. In the Eulerian framework, the dry air within a small volume of size *dxdyδz* centered on **x** ≡ (*x*, *y*, *z*) contributes the following amount to the vertical momentum of cloudy updrafts in *I*(*z*) at time *t*:

Averaging this over *x* and *y* and dividing by *δz* gives the mean mass flux *M* as calculated in the Eulerian framework,

In the Lagrangian framework, the dry air associated with particle *i* contributes the following amount to the vertical momentum of cloudy updrafts in *I*(*z*) at time *t*:

Here, *W _{i}* is the vertical velocity of particle

*i*. Summing this over all particles and dividing by the control volume

*L*gives

_{x}L_{y}δz*M*as calculated in the Lagrangian framework,

Figure 5a plots *M* as calculated in the Eulerian framework (solid line) and in the Lagrangian framework (circles) averaged over the first hour of the simulation. The agreement between the Lagrangian and Eulerian calculations is excellent.

In the Eulerian framework, the entrainment rate is calculated using the direct measurement method of Romps (2010). In particular, we diagnose the sources of an Eulerian continuity equation for active air in each grid cell. This source is averaged over the time it takes the boundary of active air to pass over the grid cell. If this averaged source is positive, then it is counted as entrainment. If it is negative, then it is counted as detrainment. As a consequence of this averaging procedure, Eulerian direct measurement cannot resolve the detailed partitioning between entrainment and detrainment on time scales less than the time it takes the updraft boundary to move 2*δx* = 100 m. For typical speeds of large, turbulent eddies of about 5 m s^{−1}, this implies a time scale of about 20 s. For example, if there is actually 10 kg of entrainment and 3 kg of detrainment in that 20 s, then Eulerian direct measurement will report 10 − 3 = 7 kg of entrainment and 0 kg of detrainment during that time. This effect will tend to bias low the entrainment rates obtained from Eulerian direct measurement.

In the Lagrangian framework, the local entrainment rate *e*(**x**, *t*) is the number of particles that switch from inactive to active in each grid cell over some averaging time. Here, the entrainment rate is computed at each time step, *δt* = 2 s. The horizontally averaged entrainment rate is then calculated in the Lagrangian framework as

where *H* is the Heaviside unit step function: *H*(*x*) = 1, if *x* > 0, and *H*(*x*) = 0, otherwise. We refer to this as Lagrangian direct measurement.

Figure 5b plots the rate of entrainment averaged over the first hour using Eulerian direct measurement (solid line) and Lagrangian direct measurement (circles). The Eulerian and Lagrangian entrainment rates are qualitatively similar. Both Eulerian and Lagrangian measurements show strong cloud-base entrainment at *z* ≃ 500 m, which is followed by a local minimum at *z* ≃ 1000 m. For 1000 < *z* < 8000 m, the Lagrangian entrainment rate profile is very similar to the Eulerian entrainment rate, but it is consistently larger in the free troposphere by about 20%.

In Fig. 5c, the net entrainment rate, defined by the entrainment rate subtracted by the detrainment rate, is plotted to test the consistency between the Lagrangian and Eulerian measurements. The Lagrangian detrainment rate is computed similar to Eq. (8) but with *H*[_{i}(*t* − *δt*) − _{i}(*t*)]. It is found that the net entrainment rate shows an excellent agreement between the Lagrangian and the Eulerian measurement methods. This result indicates that, although the mass of air entrained into a cloud computed by the Lagrangian method is higher than that of the Eulerian measurement, the Lagrangian estimate also predicts a higher detrainment rate so that the mass flux of the cloud estimated by the Lagrangian method is consistent with that of the Eulerian method.

## 5. Residence times

When we browse through the particle trajectories that pass through the cloud, we find that the activity of a typical Lagrangian particle changes several times on a fairly short time scale, from tens of seconds to a few minutes. These fluctuations occurring on a time scale of less than about 20 s may contribute to the higher entrainment rate measured by the Lagrangian method. For demonstration, we randomly choose a particle from a pool of the Lagrangian particles that have an initial location below the cloud base, *Z*(0) < 500 m; a final destination near the cloud top, *Z*(1 h) > 7000 m; and at least four reentry events (detrainment followed by a subsequent entrainment). The trajectory of this randomly selected particle is shown in Fig. 6. This figure displays its time series of vertical location (Fig. 6a), activity operator (Fig. 6b), vertical velocity (Fig. 6c), and equivalent potential temperature (Fig. 6d). As can be seen from Fig. 6b, the particle exhibits multiple entrainment and detrainment events on a short time scale. Comparing Figs. 6b,c, it is clear that the rapid detrainment and reentry events for this particle are due to the oscillations in its vertical velocity. For this particular trajectory, only the first entrainment event at *t* ≃ 7 min is associated with condensation. After that initial entrainment, *q _{c}* remains well above

*q*

_{threshold}with

*q*≃ 10

_{c}^{−3}–10

^{−2}. Figure 6d shows the changes of

*θ*, which is conserved for all adiabatic and reversible transformations. The equivalent potential temperature of the particle remains constant until the particle reaches

_{e}*z*≃ 1500 m, suggesting that this particle remains in a protected core up to 1500 m. This will be discussed further in section 6. A rapid decrease of

*θ*is observed after 10 min due to the mixing with environmental air entrained laterally. From 10 min,

_{e}*θ*drops about 8 K in less than 2 min, suggesting that the time scale of mixing is relatively small. A similar rapid decrease of the potential temperature has been observed in the bubble simulation of Lasher-Trapp et al. (2005).

_{e}An important time scale in the entrainment process is the time a parcel remains in the cloud once it has entrained. We denote the length of time between a particle’s entrainment and subsequent detrainment by ^{e}; we will refer to this as the “residence time.” Similarly, we define the length of time between a particle’s detrainment and subsequent entrainment by ^{d}; we will refer to this as the “reentry time.” Figure 7a shows the probability density function (PDF) of particle entrainment events *p _{e}* as a function of

^{e}and the PDF of particle detrainment events

*p*as a function of

_{d}^{d}. It is found that the PDFs are strongly peaked around zero time. Both

*p*(

_{e}^{e}) and

*p*(

_{d}^{d}) decrease by about three orders of magnitude between = 0.5 and = 20 min. The behavior of the PDFs at small suggests that air parcels are moving in and out of the cloud (defined by = 1) on a very short time scale due to small-scale turbulent eddies. As we will see shortly, the discrepancy in Fig. 5b can be explained by this high-frequency oscillation, which is not captured by Eulerian direct measurement.

Although Fig. 7a shows that the most common entrainment events are those that bring air into the cloud for less than 20 s, this does not imply that the cloud is ventilated on this time scale. Instead, these fast entrainment/detrainment events account for a small fraction of the cloud mass precisely because these parcels reside in the cloud for such a short time. The blue circles in Fig. 7b plot the PDF of cloud updraft mass Ψ as a function of residence time (as opposed to the PDF of particle entrainment events as a function of residence time). At every time step, each active particle (or, rather, its corresponding mass *δm*) is binned according to the residence time of that particle within the cloud (i.e., the time of its next detrainment event minus the time of its previous entrainment event). In other words, Ψ is related to *p _{e}* and the cloud mass as

in which *T* is the total simulation time, *m*(*t*) is the mass of the cloud at time *t*, and Λ(*t*) is the total rate of entrainment at *t*, as shown:

We see from Fig. 7b that the most common type of parcel in the cloud has a residence time of ^{e} ≃ 1.5 min. Also, Ψ(^{e}) decreases slowly for ^{e} > 1.5 min, compared to *p*(^{e}), which exhibits a rapid exponential decay with ^{e}. The cumulative distribution function (CDF) of Ψ(^{e}), the dashed black line in Fig. 7b, shows that 90% of the cloudy-updraft mass comprises parcels with residence times greater than 1 min.

As argued above, Eulerian direct measurement is likely to be biased low because of the way it must average over the time it takes the active boundary to move two grid steps (a time estimated to be about 20 s). It is also possible that the Lagrangian method is biased high due to numerical artifacts leading to spurious oscillations in the time series of a particle’s *w* or *q _{c}*, which would lead to spurious oscillations in . It is not possible at this point to say which of either Eulerian or Lagrangian direct measurement is more accurate. Confidence in both methods is buoyed by the fact that they give quantitatively similar estimates of the entrainment rate, as shown by Fig. 5b.

To confirm that the difference between the Eulerian and Lagrangian methods lies in the fast recycling events, we recompute the Lagrangian estimate of entrainment rate by discarding particles’ entrainment and detrainment couplets that are separated in time by less than a chosen filter *τ*. In Fig. 7c, the Lagrangian entrainment rate with *τ* = 0, 20, 60, and ∞ s are displayed together with the Eulerian measurement; *τ* = ∞ indicates that an air parcel is not counted as reentrained after its first entrainment event. As expected, the Lagrangian entrainment rate is a decreasing function of *τ* except for the cloud-base entrainment. When the fast reentry events associated with the sharp peak of *p _{e}* and

*p*for < 20 s are filtered, the Lagrangian entrainment rate is in excellent agreement with the Eulerian measurement. The decrease of the Lagrangian entrainment rate at larger

_{d}*τ*indicates that a large amount of the entrained air is not directly from the environment but from the cloud itself. Among all 368 702 entrainment events in the simulation, 61% (223 833) are from reentry events. In contrast, the cloud-base entrainment is mainly due to the saturation of warm air as it ascends adiabatically from the ground to the lifting condensation level and, hence,

*τ*has no effect on the cloud-base entrainment. The entrainment rate estimated with

*τ*= ∞ sets the lower bound of the Lagrangian entrainment rate. This lower bound of Lagrangian entrainment rate is about half of the Eulerian direct measurement, which is closer to a bulk-plume estimate (Romps 2010). Figure 7d shows the fractional entrainment rate calculated by the ratio of the entrainment rate to the mass flux of active air: that is,

*e*/

*M*. The fractional entrainment rate is diagnosed to have a value around 1–2 km

^{−1}, and the lower bound of the Lagrangian fractional entrainment rate is about 0.5 km

^{−1}.

## 6. Heights of origin

In this section, we focus on the origin of the cloud air as estimated by the Lagrangian particles. First, we show snapshots of the most recent entrainment height *z*_{ent} of the cloudy-updraft or active air at 10 (developing stage), 20 (fully developed stage), and 30 min (dissipating stage) in Figs. 8a–c. The height of entrainment is computed by an azimuthal average and shown as a function of the radial distance from the center of the domain *r* and the altitude *z*. In the developing stage (Fig. 8a), most of the active air in the cloud has been entrained near the cloud base. In Fig. 8b, it is shown that there is a cloud core for *z* < 2 km, where most of the air in that region comes from the cloud-base entrainment. In the middle of the cloud, 4 < *z* < 7 km, the cloud air has most recently entrained, on average, from about 1 to 2 km below. In Figs. 8a,b, there is a narrow band near the cloud top, where *z*_{ent} ≃ *z*. This sharp front of the newly entrained air indicates cloud-top entrainment. Finally, in the dissipating stage (Fig. 8c), the core region has disappeared and the air remaining in the cloud is mostly entrained near the current height of the air parcels, that is, *z*_{ent} ≃ *z*. The empty holes of the active air observed near the center of the cloud are due to the rain-induced downdrafts.

In section 5, we have shown that an air parcel may experience multiple detrainment and reentrainment events. As a consequence, the entrainment height *z*_{ent} as it is defined here may not be a good measure of the origin of entrained air. Instead, the original height of active air *z*_{org} (i.e., the height of the Lagrangian particles at *t* = 0) is shown in Figs. 8d–f. In general, *z*_{org} is less than *z*_{ent}. For example, in Figs. 8b,e, the average *z*_{org} for the cloud in 5 < *z* < 7 km is around *z*_{org} ∈ (2, 4) km, while *z*_{ent} in that region indicates the cloud air is from 4 < *z* < 6 km. This difference in *z*_{org} and *z*_{ent} can be explained by reentrainment events due to energetic turbulent eddies.

In convective parameterizations, it is standard to assume that air entrained at *z* has thermodynamic properties equal to that of the environment at *z*. To assess this assumption, we plot the mean equivalent potential temperature of air when it entrains as a function of where it entrains (*z*_{ent}) for the set of particles that remains in the cloud at *t* = 15 min (Fig. 9a) and *t* = 20 min (Fig. 9b). The mean *θ _{e}* in the cloud and the mean

*θ*in the environment are also plotted for reference. This shows that the equivalent potential temperature of entraining air is much higher than . This is, in large part, because a significant fraction of the entrained air is previously detrained by the cloud and, therefore, has a higher

_{e}*θ*(i.e., is moister) than the environment. Another contribution comes from numerical diffusion. The red circles in Fig. 9 show the mean equivalent potential temperature of air when it first entrains as a function of where it first entrains (

_{e}*z*

_{ent,0}) for the same set of particles in the cloud at

*t*= 15 (Fig. 9a) and

*t*= 20 min (Fig. 9b). As expected, is lower than . Nevertheless, is still about 1–2 K higher than

*θ*

_{env}. Below 6 km, this is likely partially due to the lifting of parcels (from lower heights with higher

*θ*) before they entrain. There is likely also a contribution to the

_{e}*θ*of entrainment parcels from moistening by numerical diffusion in the vicinity of the cloud.

_{e}Figure 9 also shows that, as the top of the cloud is approached, the cloud *θ _{e}* decreases rapidly and eventually approaches . This is consistent with the entrainment of environmental air at the cloud top as observed in Fig. 8. Further investigation has shown that an air parcel near the cloud top ascends, on average, a few hundred meters before being entrained. The detailed mechanism of this cloud-top entrainment is beyond the scope of the present study.

## 7. Concluding remarks

In this study, we introduce a Lagrangian framework to diagnose cloud processes in large-eddy simulations, and we use that framework to study a single cumulus congestus. Although it is standard practice to include a subgrid-scale (SGS) model in the equation for particle motion (e.g., Heus et al. 2008), such an SGS model is intentionally omitted here. Another choice that is made deliberately in this study is to initialize the particles with a random spatial distribution subject to the constraint that the specific number density is spatially invariant in a statistical sense. These two choices, along with a careful choice of numerics, guarantee that the constant specific number density remains invariant in space and time (see Fig. 3b). This allows the Lagrangian histories to serve as a conjugate solution for the resolved Eulerian flow. This approach has some obvious benefits over the Eulerian tracer method of Romps and Kuang (2010a,b), which can track only one real number per additional tracer. Here, a single Lagrangian particle can carry a wealth of history information with little computational expense.

Given the kinematic and thermodynamic time series for each of the Lagrangian particles, it is straightforward to calculate the convective entrainment rate [see Eq. (8)]. We refer to this procedure as Lagrangian direct measurement of entrainment. As shown in Fig. 7c, the Lagrangian direct measurement gives an entrainment rate (solid blue circles) that is very similar to the rate obtained from Eulerian direct measurement (black line; Romps 2010), but it is consistently higher by about 20% in the free troposphere. It is shown here that this is due to the recycling of particles into and out of the cloud on short time scales (less than 20 s) that are not resolved by the Eulerian method. By filtering out those fast recycling events, Lagrangian direct measurement produces an entrainment rate (open green circles in Fig. 7c) that is in excellent agreement with the Eulerian result.

Studying the kinematic and thermodynamic histories of the Lagrangian particles has revealed additional information on the entrainment process. Kinematically, it is found that 61% of all entrainment events are actually reentrainment of air previously detrained by the cloud. Nevertheless, about 90% of the total mass of the cloud is composed of air parcels whose in-cloud residence time is larger than 1 min. In agreement with the observed fractional entrainment rate of 1–2 km^{−1}, the mean entry height of cloud updraft air ranges from 0 to 1 km below the height of observation (see Figs. 8a–c). In agreement with the observation that the majority of entrainment events are actually reentrainment events, the mean original height of cloud-updraft air ranges from 0 to 2 km below the height of observation (see Figs. 8d–f). Contrary to Heus et al. (2008), who concluded that there is not significant cloud-top entrainment in shallow cumulus, some clear evidence is found for cloud-top entrainment by the cumulus congestus in this study. Further work will be needed to quantify the relative magnitudes of lateral entrainment and cloud-top entrainment.

Thermodynamically, the entraining parcels of air are found to have an equivalent potential temperature *θ _{e}* (blue circles in Fig. 9) that exceeds that of the environment (dashed line). This is primarily a consequence of the large number of reentrainment events. Calculating

*θ*of entraining air for parcels’ first entrainment event only (red circles),

_{e}*θ*is found to be much closer to the environment.

_{e}## Acknowledgments

This work was supported by the Laboratory Directed Research and Development Program of Lawrence Berkeley National Laboratory under U.S. Department of Energy Contract DE-AC02-05CH11231. This research was also supported in part by the director, Office of Science, Office of Biological and Environmental Research of the U.S. Department of Energy under Contract DE-AC02-05CH11231 as part of their Atmospheric System Research Program. This research used computing resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract DE-AC02-05CH11231, and computing resources of the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant OCI-1053575. Feedback from three anonymous reviewers helped improve this manuscript.

### APPENDIX

#### Sensitivity to Activity Operator

The entrainment rate and other statistics will depend on the definition of cloud or cloud updrafts as defined by the activity operator . In this study, we have used = *H*(*q _{c}* −

*q*

_{threshold})

*H*(

*W*−

*w*

_{threshold}), where

*w*

_{threshold}= 1 m s

^{−1}and

*q*

_{threshold}= 10

^{−5}kg kg

^{−1}. To show the sensitivity of the Lagrangian statistics to , the residence time and reentry time PDFs (

*p*and

_{e}*p*) are recomputed using =

_{d}*H*(

*q*−

_{c}*q*

_{threshold}) with

*q*

_{threshold}= 10

^{−5}kg kg

^{−1}(Fig. A1). Here, the original activity operator is denoted as

_{w,q}and the new operator is

_{q}. It is shown that all of the probability density functions are peaked near the origin and exhibit a rapid exponential decay at short times, suggesting the short-time-scale recycling is a robust feature. As expected, while

*p*for both

_{e}_{w,q}and

_{q}exhibit exponential decays,

*p*for

_{e}_{q}decreases much more slowly compared to that of

_{w,q}. More interestingly, it is found that the reentry probability density functions for

_{w,q}and

_{q}collapse onto one curve for

^{d}> 5 min.

Figure A2 shows the trajectory of a sample particle. This particle is selected randomly among the set of particles whose initial location is below the cloud base, whose final location is above 7 km, and that experiences at least six reentry events based on _{q}. In Fig. A2a, it is shown that the particles’ activity oscillates on a short time scale. Unlike the high-frequency activity oscillations seen in Fig. 6b, the high-frequency oscillations seen here are related to the oscillation of *q _{c}* around

*q*

_{threshold}(Fig. A2c).

## REFERENCES

*J. Atmos. Sci.,*

**55,**3440–3455.

*J. Atmos. Sci.,*

**62,**1175–1191.

*J. Atmos. Sci.,*

**52,**2851–2868.

*J. Atmos. Sci.,*

**69,**1118–1136.