A method is introduced for directly measuring convective entrainment and detrainment in a cloud-resolving simulation. This technique is used to quantify the errors in the entrainment and detrainment estimates obtained using the standard bulk-plume method. The bulk-plume method diagnoses these rates from the convective flux of some conserved tracer, such as total water in nonprecipitating convection. By not accounting for the variability of this tracer in clouds and in the environment, it is argued that the bulk-plume equations systematically underestimate entrainment. Using tracers with different vertical profiles, it is also shown that the bulk-plume estimates are tracer dependent and, in some cases, unphysical. The new direct-measurement technique diagnoses entrainment and detrainment at the gridcell level without any recourse to conserved tracers. Using this method in large-eddy simulations of shallow and deep convection, it is found that the bulk-plume method underestimates entrainment by roughly a factor of 2. The directly measured entrainment rates are then compared to cloud height and cloud buoyancy. Contrary to existing theories, fractional entrainment is not found to scale like the inverse of height, the cloud buoyancy, or the gradient of cloud buoyancy. On the other hand, fractional detrainment is found to scale linearly with cloud buoyancy. Finally, direct measurement is used to diagnose the spatial distribution of entrainment and detrainment during the evolution of an individual deep cumulonimbus.
Since the 1970s, bulk-plume equations have been used to diagnose convective entrainment and detrainment rates from observations of the large-scale budgets of deep and shallow convection (Yanai et al. 1973; Esbensen 1978). More recently, bulk-plume equations have been used to diagnose the fractional rates of entrainment ε and detrainment δ from large-eddy simulations (LES) of shallow convection (Siebesma and Cuijpers 1995; Siebesma et al. 2003). But the bulk-plume model invokes two approximations that may not be valid. First, it is assumed that the properties of the detrained air are the same as the average properties in the plume. Second, the properties of entrained air are assumed to be the same as the average properties of environmental air. Since both the plume and the environment are represented by their bulk properties in this approach, the term “bulk-plume” is a bit of a misnomer, but it is certainly less cumbersome than “bulk-plume-and-bulk-environment.” Therefore, we will refer to the fractional entrainment and detrainment rates obtained from these approximations as the bulk-plume ε and δ.
The bulk-plume ε and δ are diagnosed using the budgets for the convective fluxes of two quantities: dry air and some conserved tracer. In shallow, nonprecipitating convection, there are two different variables that have been used for that conserved tracer: total water qt and liquid-water potential temperature θl. It has been argued that the accuracy of the bulk-plume ε and δ can be gauged by calculating ε and δ twice—once with qt and once with θl—and then comparing the two sets of estimates (Siebesma 1996). Indeed, it is found that the two sets of estimates agree (Siebesma and Cuijpers 1995). But qt and θl are very highly anticorrelated both within clouds (see Fig. 2 of Romps and Kuang 2010b) and in the environment (see Fig. 1 of Siebesma et al. 2003), so using θl is almost identical to using −qt. This means that the accuracy of the bulk-plume ε and δ has not yet been independently tested.
There are many reasons why it is important to check the validity of the bulk-plume ε and δ. For example, if these values are dependent on the tracer distribution, then the bulk-plume ε and δ obtained using one tracer will give the wrong answer when used to predict the transport of a tracer with a different vertical distribution. Or, if we take the bulk-plume ε and δ at face value, when they are actually too low, then we may underestimate the sensitivity of convection to changes in free-tropospheric humidity. Other problems will arise if we try to use the bulk-plume ε and δ in a convective parameterization that does not use the bulk-plume approximations, such as a buoyancy-sorting scheme (Kain and Fritsch 1990) or a stochastic parcel model (Romps and Kuang 2010b). Even using the bulk-plume ε and δ in a bulk-plume model may lead to an incorrect estimate of the processing of trace chemicals by aqueous reactions. Finally, and perhaps most importantly, we may miss an opportunity to learn about the true entrainment and detrainment rates, which may be more fundamental than the effective rates estimated by the bulk-plume equations.
This paper has two main objectives: 1) to assess the validity of the bulk-plume entrainment and detrainment rates and 2) to measure entrainment and detrainment directly. Section 2 reviews the standard bulk-plume method and, in so doing, presents an extension of the method to deep, precipitating convection. In that section, it is argued from first principles that the bulk-plume ε and δ systematically underestimate the true ε and δ. Section 3 describes the large-eddy simulations of shallow and deep convection that are used in this study. In section 4, those simulations are used to perform sanity tests on the bulk-plume method, which performs poorly. To remedy this situation, a new method is introduced in section 5 that allows for the direct measurement of entrainment and detrainment at the gridcell level. In section 6, the directly measured entrainment and detrainment rates are subjected to sanity tests, which yield positive results. The direct-measurement and bulk-plume methods are compared quantitatively in section 7, confirming a large underestimation from the bulk-plume method. Section 8 explores the implications of the new results for popular models of entrainment and detrainment, such as ε ∼ 1/z and ε ∼ db/dz, where b is the cloud buoyancy. Since direct measurement allows entrainment and detrainment to be viewed as functions of space and time, section 9 presents the distributions of entrainment and detrainment during the development of an individual deep cumulonimbus. Section 10 summarizes the findings.
2. Bulk-plume method—Theory
Before we define entrainment, we must answer the question “Entrainment into what?” In particular, we must define some criteria for dividing the atmosphere into two categories between which mass is exchanged. Referring to the two categories as “active” and “inactive,” we say that a Lagrangian parcel “entrains” when it flips from the inactive category to the active category. Likewise, a parcel of air “detrains” when it flips from active to inactive. To maintain full generality in the equations that follow, let us define an “activity operator” 𝒜 that is equal to one at the location of active air and zero at the location of inactive air. Later, we will choose to define active air as that which has a condensate mixing ratio and a vertical velocity above some thresholds. In that case, 𝒜 would take the form
where qc is the mixing ratio of condensates, and w is the vertical velocity. Denoting the horizontal average by angled brackets, 〈𝒜〉 is the fractional area covered by active air as a function of height, which is often denoted by σ. We define the entrainment rate as the local rate at which air flips from inactive to active. This quantity has units of kg m−3 s−1. Similarly, we define the detrainment rate as the local rate at which air flips from active to inactive. By definition, e and d are nonnegative quantities and ed is identically zero, since no point can be simultaneously entraining and detraining.
The bulk-plume method for diagnosing e and d relies on the continuity equations for three densities: dry air ρ, active dry air 𝒜ρ, and a tracer ϕρ, where ϕ is the mixing ratio. Since dry air is conserved, it has no sources or sinks. For active air, the sources and sinks are, by definition, entrainment and detrainment. For the ϕ tracer, we will only require that its sources and sinks Sϕ be zero where 𝒜 = 1. For example, in nonprecipitating convection, we can set ϕ equal to the total-water mixing ratio qt. More generally, we can use any artificial tracer whose sources and sinks are restricted to inactive air. This will give us the freedom to specify the horizontally averaged vertical profile of such a tracer.
Keeping the ϕ tracer fully general for the time being, the three continuity equations are
Here, we have used the condition that 𝒜Sϕ = 0. Finally, using Eq. (2) to pull 𝒜, ρ, and inside the derivatives, we are left with the continuity equation for active ϕ tracer,
Here, we have assumed that the horizontal average is performed over an area that is either doubly periodic or has u = υ = 0 m s−1 at the lateral boundaries. Finally, by defining ϕe and ϕd as the average ϕ of entraining and detraining parcels
these equations may be written as
Ultimately, what we want to diagnose from these equations are the horizontally averaged rates of entrainment 〈e〉 and detrainment 〈d〉. In a cloud-resolving simulation, only the terms on the left-hand side of Eqs. (7) and (8) are easily measured. Therefore, we have only two equations but four unknowns: 〈e〉, 〈d〉, ϕe, and ϕd. To solve this problem, the bulk-plume method approximates the average ϕ of entraining air, ϕe, as equal to the average ϕ of inactive air and approximates the average ϕ of detraining air, ϕd, as equal to the mass-flux-weighted average of active air:
where the subscript ϕ has been added to e and d to denote the fact that these estimates were obtained using the mixing ratio ϕ. Equations (11) and (12), and variations on them, have been used to diagnose fractional entrainment and detrainment profiles with ϕ = qt in many simulations of nonprecipitating convection (e.g., Tiedtke 1989; Schumann and Moeng 1991; Siebesma and Cuijpers 1995; Siebesma 1996).
In deep convection, which is precipitating, there is no naturally conserved tracer. Moist static energy h and equivalent potential temperature θe are only approximately conserved. The hydrostatic approximation used to derive h is violated by nonhydrostatic convection. Equivalent potential temperature is not conserved under diffusion, and neither h nor θe is strictly conserved when there is precipitation. Therefore, we introduce a passive tracer that will be explicitly conserved within active air, so it can serve as the conserved quantity in the mass-flux budgets. In particular, we use the “purity” tracer, which is advected in the same way as dry air and has sources and sinks only some distance away from active air (Romps and Kuang 2010a,b). Those sources and sinks are specified in the following way. First, the purity ϕ is reset to one everywhere below the cloud base at every time step. Above the cloud base, the purity ϕ is set to zero everywhere at every time step except within the vicinity of active air. A grid cell is defined to be in the vicinity of active air if a cube—seven grid spaces wide and centered on the grid cell—has any active air within it. This scheme guarantees that, above the cloud base, inactive air has an average purity very close to zero. The reason it is called the purity tracer is that its mixing ratio within active air is very nearly equal to the fraction of air that has come directly from below the cloud base. When ϕ is the purity-tracer mixing ratio, 〈ϕρ(1 − 𝒜)〉 is very nearly zero. Equations (11) and (12) then simplify dramatically for steady-state convection. Defining fractional entrainment and detrainment as ε ≡ 〈e〉/〈ρw𝒜〉 and δ ≡ 〈d〉/〈ρw𝒜〉, Eqs. (11) and (12) can be approximated as
The expressions for eϕ and dϕ are only as good as the bulk-plume approximations for ϕe and ϕd. For general ϕe and ϕd, the expressions for eϕ and dϕ are
To assess the effect that different choices for ϕe and ϕd have on the estimates of entrainment and detrainment, we calculate the partial derivatives of εϕ and δϕ with respect to ϕe and ϕd:
Since the average purity (total water) of active air is greater than the average purity (total water) of inactive air, ϕd − ϕe is positive. And, since entrainment and detrainment are, by definition, nonnegative quantities, one would certainly hope that Eqs. (15) and (16) would give positive values for εϕ and δϕ. Therefore, we see from these derivatives that εϕ and δϕ increase when ϕe increases or ϕd decreases.
With these derivatives in hand, we can now see why the bulk-plume method systematically underestimates ε and δ. In Eq. (10), ϕd is set equal to the mass-flux-weighted average of ϕ within active air. This approximation misses an important point. Convecting clouds are highly heterogeneous and, as shown by Romps and Kuang (2010b), that heterogeneity is due to variations in the amount of entrained air, as measured by the purity tracer. Furthermore, the purity, total water, and buoyancy of cloudy air parcels are all highly correlated with each other. In a heterogeneous cloud with a large variance of buoyancy, we expect detraining parcels to be those parcels that are less buoyant than average. Given the high correlations of buoyancy with purity and total water, we also expect detraining parcels to have mixing ratios of total water and purity that are below average. Therefore, for ϕ equal to total water or purity, the bulk-plume equation [Eq. (10)] systematically overestimates ϕd.
What about ϕe? Equation (9) approximates ϕe as the average values within inactive air. This might be a suitable approximation if convecting clouds or their constituent plumes were initiated randomly in the horizontal plane—but they are not. A single deep-convective cloud can contain many plumes bubbling up within it, and on larger scales clouds develop preferentially in clusters. Therefore, a cloud tends to entrain air that was recently detrained from itself or from other clouds. Since that detrained air has higher total water and higher purity than the average inactive air, the bulk-plume equation [Eq. (9)] systematically underestimates ϕe when ϕ is the total water or purity.
For ϕ equal to total water or purity, the traditional method overestimates ϕd and underestimates ϕe. From the derivatives of εϕ and δϕ, we see that the effects of too high a ϕd and too low a ϕe work in the same direction: they give εϕ and δϕ that are below their true values. Thus, the bulk-plume approximations lead to a systematic bias that causes fractional entrainment and detrainment rates to be underestimated. The obvious follow-up question is “By how much?”.
3. Large-eddy simulations
We will quantify the amount of underestimation using simulations of shallow and deep convection in the cloud-resolving model Das Atmosphärische Modell (DAM; Romps 2008). DAM is a three-dimensional, fully compressible, nonhydrostatic model of the atmosphere. For microphysics, the model uses the six-class Lin–Lord–Krueger scheme (Lin et al. 1983; Lord et al. 1984; Krueger et al. 1995). Since the debut of DAM, the shortwave and longwave radiation schemes have been upgraded to the Rapid Radiative Transfer Model (RRTM) (Clough et al. 2005; Iacono et al. 2008). In addition, the finite-volume advection scheme has been upgraded to use the three-dimensional Uniformly Third-Order Polynomial Interpolation Algorithm (Leonard et al. 1993) combined with the three-dimensional, monotonic flux limiter of Thuburn (1996).
For the simulation of shallow, nonprecipitating, marine convection, we use the initial conditions and large-scale forcings described in the intercomparison study of Siebesma et al. (2003). Those initial conditions and forcings were designed to replicate the undisturbed period during the third phase of the Barbados Oceanographic and Meteorological Experiment (BOMEX) (Holland and Rasmusson 1973). In this simulation, the RRTM radiation scheme is replaced by the radiative cooling profile diagnosed during that period of the BOMEX campaign. Furthermore, in the microphysics scheme, the autoconversion of cloud water to rain is turned off to ensure nonprecipitating convection. The simulation is performed in a domain that is 12.8 km wide in the horizontal directions and that has a model top at 3 km. The grid spacing is set to 50 m in all directions. After allowing the simulation to equilibrate for 3 h, statistics are collected for the next 5 h. For this simulation, active air is defined as having a liquid-water mixing ratio greater than 10−5 kg kg−1 and a vertical velocity greater than 0.5 m s−1.
For the simulation of deep, precipitating convection over the ocean, we fix the sea surface temperature to 300 K. The magnitude and direction of the shortwave radiation flux are set to constant values that provide the average daily insolation at the equinoctial equator using a zenith angle whose cosine is equal to the daily insolation-weighted average. Aside from the influence of the sea surface fluxes and the RRTM radiation, there are no forcings applied to the atmosphere. The simulation is run in a cubic domain with a 25.6-km width and periodic boundary conditions in the x and y directions. The grid spacing used here is 200 m in all directions, and the simulation has been run to a steady state for several weeks of model time. Statistics are collected over a three-week period of radiative–convective equilibrium. For this simulation, active air is defined as having a condensate mixing ratio greater than 10−5 kg kg−1 and a vertical velocity greater than 1 m s−1.
4. Bulk-plume method—Sanity tests
The left panel of Fig. 1 shows the fractional entrainment rate (solid) and fractional detrainment rate (dashed) as diagnosed from the total-water budget (green) and from the purity budget (blue) in the simulation of shallow convection. At this scale, the two sets of curves overlap each other so closely that they are difficult to distinguish from one another. Given the large correlation between purity and total water found previously by Romps and Kuang (2010b), this overlap is not surprising. What are surprising are the large excursions of the fractional detrainment going to the unphysical value of −0.04 m−1 near the cloud base and to the very large value of 0.04 m−1 in the inversion. This aspect may be understood by recognizing that where the mass flux is changing rapidly, Eq. (14) may be approximated as
Here, we have used the fact that the fractional change in the mass flux is much larger than the fractional change in the purity of active air. At the extremal height ranges of the active air, the mass flux m changes from one level to the next by multiplicative factors much greater than one. This implies that the length scale m/(dm/dz) is much less than the grid spacing. Therefore, the use of a finite-difference approximation to the derivatives effectively regularizes this length scale. In particular, −(dm/dz)/m is calculated as
where Δz is the grid spacing. When m(z + Δz) ≫ m(z), as is the case near the cloud base, this value is nearly equal to −2/Δz, which for Δz = 50 m is −0.04 m−1. At the upper reach of the mass flux, m(z + Δz) ≪ m(z), so δϕ is equal to 2/Δz = 0.04 m−1. In addition, where the fractional changes in the mass flux are large, Eq. (13) may be approximated as
so εϕ goes to zero at the extremities. This explains why εϕ is spuriously small at the cloud base when, in fact, we know the true ε is quite large there.
The right panel of Fig. 1 shows the fractional entrainment rate (solid) and fractional detrainment rate (dashed) as determined using the mass-flux budget of the purity tracer in the simulation of deep convection. Again, we see the large excursions of δϕ at the extremities of the mass-flux height range. Since the grid spacing is 200 m, δϕ reaches the unphysical value of −1/(200/2) = −0.01 m−1 near the cloud base and 0.01 m−1 near the tropopause. The other unusual feature in these plots is the excursion of εϕ to negative values between 11.4 and 13.6 km, which is the height range over which the average purity of active air 〈ϕρw𝒜〉/〈ρw𝒜〉 actually increases with height. In reality, changes in the active-air purity are determined by a balance between two effects: the detrainment of air whose purity is less than the active-air average and the entrainment of air whose purity is also less than the active-air average. Detrainment tends to increase the purity of active air, while entrainment tends to decrease it. In the upper reaches of convection, there is much more detrainment than entrainment, so the former effect wins out. But, since Eq. (10) assumes that the detrained air has the same purity as the active-air average, detrainment cannot affect the purity of active air in the bulk-plume budget. Therefore, in the bulk-plume equations, the job falls to entrainment to match the observed increase in active-air purity. Because inactive air has a lower purity than active air, only a negative entrainment can produce the observed increase in purity, which illustrates, rather strikingly, one of the inadequacies of the bulk-plume method.
According to the arguments given here, we expect the bulk-plume ε and δ to be lower than the true values. One reason for this is that detraining parcels have a different purity than the cloud average, and that difference is comparable to the difference between the mean-cloud and mean-environment purities. One way to mitigate this problem is to ramp up the cloud-to-environment difference so as to minimize the effect of the in-cloud variance. We can accomplish this with a tracer whose mixing ratio increases exponentially with height in the environment. This “exponential” tracer is advected along with dry air, but sources outside the vicinity of active air set its mixing ratio to ez/λ with some length scale λ. We choose this length scale such that the tracer increases by an order of magnitude once every 500 m in the shallow simulation and once every 2 km in the deep simulation (i.e., once every 10 grid spacings). The bulk-plume estimates of ε and δ obtained from this exponential tracer are shown in Fig. 2 alongside the corresponding estimates from the purity tracer. As expected, the use of the exponential tracer increases the estimates of ε and δ. Furthermore, the large magnitude of the increase illustrates how sensitive the bulk-plume method can be to the distribution of the conserved tracer. In other words, the bulk-plume estimates of entrainment are tracer dependent.
This tracer dependence bedevils the very reason for using the bulk-plume method. The bulk-plume approach to convection is supposed to work as follows. First, the convective (i.e., active) flux of some tracer ϕ is measured or observed. This flux is then used in Eqs. (11) and (12) to calculate εϕ and δϕ. Then, for any conserved tracer ζ whose clear-air (i.e., inactive) profile is known, we can calculate the convective flux of ζ by integrating the following equation:
Figure 3 compares the fluxes calculated from Eq. (17) (red curves) with the actual fluxes measured in the deep-convective LES (black curves); the fluxes are plotted for the middle and upper troposphere (5–15 km). For the red curves, Eq. (17) is initialized with the actual flux at one kilometer and is then integrated upward. In the top left panel, the red curve is the purity flux calculated from Eq. (17) using εpurity and δpurity. The agreement with the true flux in black is excellent, which is to be expected, because Eq. (17) simply inverts the procedure used to calculate εpurity and δpurity from the actual purity flux. The exponential flux calculated using εexponential and δexponential, shown in the bottom right panel, also produces excellent agreement. On the other hand, when εpurity and δpurity are used to calculate the flux of the exponential tracer, as shown by the red curve in the bottom left panel, the estimated flux deviates wildly from the true flux. In fact, the negative εpurity causes the estimated flux to become very large and negative when in fact the true flux is very large and positive. When εexponential and δexponential are used to calculate the purity flux, as shown in the top right panel, Eq. (17) underestimates the true purity flux by an order of magnitude over the majority of the troposphere. We see, therefore, that the discrepancy in Fig. 2 is not merely an academic matter; instead, it has grave consequences for the estimates of convective fluxes.
5. Direct measurement—Theory
What we would really like is a way to measure the entrainment and detrainment rates directly in a cloud-resolving model. Given how complicated the interface is between active and inactive air in turbulent convection, direct measurement has been deemed impractical (Siebesma 1996); however, there is a straightforward way to define the local entrainment and detrainment rates. Recall that e is the local source of active air and d is the local sink of active air. By definition, both of these are nonnegative quantities. Also, since a point in the atmosphere entrains when it switches from 𝒜 = 0 to 𝒜 = 1 and detrains when it switches from 𝒜 = 1 to 𝒜 = 0, e and d are never positive at the same point and time. Therefore, we see from Eq. (2) that e and d can be diagnosed as follows:
Where the expression is positive, there is a source of active air (i.e., entrainment). Where it is negative, there is a sink of active air (i.e., detrainment). Therefore, we will refer to as the “activity source.”
It is important to note that Eqs. (18) and (19) are exactly the same definitions of entrainment and detrainment used in the bulk-plume approach: subtracting Eq. (19) from Eq. (18) produces Eq. (2). What differs here from the bulk-plume method is the proposal that we diagnose Eqs. (18) and (19) directly from the simulation without recourse to tracers or horizontal averaging. We will refer to the direct evaluation of these two equations as “direct measurement.”
To get comfortable with these expressions, let us work through three very simple examples. The most trivial example is the case where no infinitesimal Lagrangian parcel ever flips from one state of activity to the other. In that case, ρ𝒜 is a conserved density, so the activity source is zero and, by Eqs. (18) and (19), e and d are zero. Next, consider a motionless fluid that, at time t = 0, suddenly qualifies as active. An activity operator that would give this result is 𝒜 = ℋ(t), where ℋ is the Heaviside step function. Since the fluid is motionless, ∂ρ/∂t = 0 and , so the activity source reduces to ρ∂/∂tℋ(t) = ρδ(t), where δ is the Dirac delta function. Therefore, e = ρδ(t) and d = 0, which confirms that there is a Delta-function burst of entrainment at t = 0. Finally, consider the case of constant-density air moving in the x direction with velocity u. If we define active air as air that is located at positive x, then 𝒜 = ℋ(x). Therefore, if u > 0, then e = ρuδ(x) and d = 0; in other words, there is a source of active air at the origin. If u < 0, then e = 0, d = −ρuδ(x), and there is a sink of active air at the origin. For details on the implementation of these equations in a numerical model, refer to the appendix.
6. Direct measurement—Sanity tests
Before we use direct measurement—Eqs. (18) and (19)—to answer substantive questions of science, we first want to evaluate any possible errors in the directly measured e and d. In particular, we want to investigate any errors of the form
From Eq. (7), we know that 〈e − d〉actual = ∂/∂t〈ρ𝒜〉 + ∂/∂z〈ρw𝒜〉 at each height. Therefore, our first test is to check that 〈e − d〉measured also equals this quantity. We calculate ∂/∂t〈ρ𝒜〉 + ∂/∂z〈ρw𝒜〉 by keeping track of the active-air mass flux at each vertical level in the simulation. For the simulation of shallow convection, this is plotted as the dashed line in the left panel of Fig. 4. Next, we calculate 〈e − d〉measured by averaging the directly measured e − d at each vertical level; the result is plotted as the solid line. As expected, there is very good agreement between these two curves. Although the peak values differ near the cloud base, the integrals of the curves over their positive ranges near the cloud base are in excellent agreement. For the simulation of shallow convection, these values differ by less than 1%. The right panel of Fig. 4 plots the same two curves for the simulation of deep convection. Again, there is very good agreement; the integrals of the positive values near the cloud base agree to within 2%. Therefore, we conclude that 〈e − d〉measured = 〈e − d〉actual, which implies that errore = errord.
To put an upper bound on errore and errord, we use the fact that 〈e〉measured, 〈e〉actual, 〈d〉measured, and 〈d〉actual are nonnegative. This provides us with the following inequalities:
Using the fact that errore = errord, we conclude that
We can then put an upper bound on errore and errord by simulating a turbulent flow with 𝒜 that would be expected to give either a small 〈e〉 or a small 〈d〉. For this purpose, we simulate a warm bubble that rises out of the boundary layer leading to a deep cumulonimbus. This flow is simulated in a 25.6-km wide cubic domain with periodic boundary conditions and an isotropic 100-m grid. A passive tracer is initialized to one below 1 km and to zero above, and this tracer is given no sources or sinks during the simulation. We define active air as air that has a mixing ratio of this passive tracer greater than or equal to 0.99 kg kg−1. Above 1 km, we expect detrainment to dominate over entrainment because the mixing of active and inactive parcels is likely to produce mixing ratios less than 0.99 kg kg−1.
The directly measured profiles of 〈e〉 and 〈d〉 from this bubble-release experiment are shown in the left panel of Fig. 5 plotted on a log axis. We can see from this plot that detrainment is everywhere greater than entrainment. To quantify the relative amounts of entrainment and detrainment, the right panel plots 〈e〉/〈d〉. Broadly speaking, the curve may be described as having values between 0.5 and 0.8 in the first 500 m, a drop to a local minimum of 0.03 at 1 km, and values in the vicinity of 0.2 throughout most of the troposphere.
How can we understand the relatively large values of 〈e〉/〈d〉 near the surface? At the end of the simulation, the fraction of air that is active ranges from zero at 1 km to about 60% at the surface. In other words, there is still a large reservoir of active air near the surface. Furthermore, during the simulation, the inactive air in the lowest 500 m had an average purity of 95%. Thus, the air in the lower 500 m is split fairly evenly between active air with a mixing ratio greater than 99% and inactive air with an average mixing ratio of 95%. This is a decent environment for generating entrainment through the diffusion of tracer from air with a mixing ratio of 1 to air whose mixing ratio is just below 0.99 kg kg−1. In contrast, at a height of 1 km, which is the original boundary between tracer values of 1 and 0, there is intense detrainment as active air rising from below 1 km abruptly encounters air with a mixing ratio of 0 kg kg−1. This large amount of detrainment explains the low 〈e〉/〈d〉 there.
Throughout the rest of the troposphere, where there is fully developed turbulence, the ratio takes values around 0.2. The question, then, is how to interpret this. Does the actual flow really lead to 20% as much entrainment as detrainment? Or is this entrainment an artifact of numerical error in the calculation of 〈e〉 and 〈d〉? Or did the algorithm for calculating 〈e〉 and 〈d〉 accurately measure entrainment and detrainment that the numerical simulation generated spuriously? Unfortunately, it is not clear how to distinguish between these possibilities. On the bright side, however, this simulation allows us to put a rather stringent upper bound on the errors:
In other words, the worst-case scenario is that 0.4/1.2 = 33% of 〈e + d〉measured is spurious.
Having performed well on these tests, we can now apply the direct-measurement method to the simulations of shallow and deep convection. The left panel of Fig. 6 plots the resulting ε and δ for the simulation of shallow, nonprecipitating convection. The most striking thing about this plot is the large ranges of ε and δ. The range for ε is 1.6–29 km−1, and the range for δ is 2.5–120 km−1. Another interesting feature is that entrainment dominates over detrainment below 700 m, but above that height detrainment exceeds entrainment. Beginning at the bottom of the inversion at 1500 m, the detrainment increases rapidly, as expected.
The plot of ε and δ for the simulation of deep convection is shown in the right panel of Fig. 6. Here, we see that many of the qualitative features of shallow convection are present in deep convection as well. For example, ε and δ span very large ranges of values: from 0.18 to 11 km−1 for fractional entrainment and from 0.73 to 15 km−1 for fractional detrainment. Furthermore, entrainment dominates over detrainment in the first kilometer—where convecting clouds are born—and detrainment increases rapidly between 10 and 15 km, where convection terminates. Between 1 and 10 km, the detrainment rate typically exceeds the entrainment rate, but there are a couple of altitude ranges (near 4 and 8 km) where entrainment exceeds detrainment.
7. Direct measurement versus bulk plume
As we have seen, the bulk-plume method gives unphysical values for ε and δ near the bottom and top of the cloud layer. For that reason, we will compare the bulk-plume and direct-measurement methods where the bulk-plume method has at least some chance of giving the correct answer. For shallow convection, we will focus on the main part of the cloud layer, starting 100 m above the cloud base (725 m) and ending 100 m below the inversion (1400 m). The left panel of Fig. 7 shows ε and δ together with εϕ and δϕ for two tracers—total water and purity—zoomed in on this region of interest. We see that the bulk-plume method underestimates the directly measured entrainment and detrainment. In this altitude range, the directly measured fractional entrainment ranges from 2.2 to 2.8 km−1, but the bulk-plume fractional entrainment ranges from only 0.8 to 1.7 km−1. For detrainment, the directly measured values range from 3.5 to 4.1 km−1 and the bulk-plume values range from 2.2 to 2.9 km−1.
For the simulation of deep convection, we will focus on the main part of the cloud layer, starting 1 km above the cloud base (1.5 km) and ending 1 km below the cold-point tropopause (13.7 km). The right panel of Fig. 7 shows ε and δ, as well as εϕ and δϕ calculated using the purity tracer, zoomed in on this region of interest. As in shallow convection, εϕ and δϕ are smaller than the directly measured values. In this altitude range, the directly measured ε ranges from 0.5 to 1.3 km−1, but εϕ ranges from −0.6 to 0.6 km−1. For detrainment, δ ranges from 0.7 to 3.7 km−1, and the bulk-plume δϕ ranges from 0.2 to 2.9 km−1.
It is interesting to ask how these results would change given a different activity operator. For example, what if eddies or gravity waves are causing parcels to oscillate in and out of the active category because of the w ≥ 1 m s−1 criterion? Exactly how this phenomenon would affect the diagnosed entrainment rates is far from clear: it would contribute to the directly measured entrainment (as the same parcel flips to active over and over again), it would contribute to ∂/∂z〈ρw𝒜〉 and ∂/∂z〈ϕρw𝒜〉 in the numerator of Eq. (15), and it would increase the overall mass flux 〈ρw𝒜〉, which enters the equation for fractional entrainment through ε = 〈e〉/〈ρw𝒜〉. Given all of this complexity, it is at least conceivable that the net effect of the w ≥ 1 m s−1 criterion would be to increase the difference between the directly measured and bulk-plume estimates of ε. To test this, we can diagnose both simulations using 𝒜 = ℋ(qc − 10−5) so as to eliminate any contribution from parcels flipping in and out of the active category because of oscillations in their vertical velocity. The results, shown in Fig. 8 indicate that the gap between the two estimates has not closed. In fact, with this definition of active air, the bulk-plume estimates are more erroneous in both absolute and relative terms. Above 13 km, 〈ρw𝒜〉 is negative because of overshooting clouds having a positive e − d before they sink back down. Since ε = 〈e〉/〈ρw𝒜〉 and ϕ = 〈d〉/〈ρw𝒜〉, ε and δ become negative. This is not a failure of direct measurement but rather a failure of the notion of a “rate per distance,” which has its conceptual origin in plume models with strictly nonnegative vertical velocities.
We see from these comparisons that the bulk-plume method produces estimates of entrainment and detrainment that are significantly lower than the directly measured values. Although this error varies by height and by the type of convection, the bulk-plume equations can be grossly characterized as underestimating entrainment by at least a factor of 2. These differences are greater than the upper bound on the error of the directly measured values. Therefore, we conclude with confidence that the bulk-plume entrainment and detrainment rates are systematically biased low.
8. Theories for ε and δ
Now that we have reliable estimates of ε and δ from two types of convection, we can test various theories for ε and δ. One such theory is that ε scales like the inverse of the largest eddy sizes in the cloud, which, as the argument goes, scale with the height of the cloud; in other words, ε ∼ 1/z (e.g., Siebesma 1996). It has been argued that this theory is supported by the general tendency of ε to decrease with height as diagnosed using the bulk-plume method in shallow convection (de Roode et al. 2000; Siebesma 1996; Siebesma et al. 2003). As a result, the proposed ε ∼ 1/z relationship has been used in convective parameterizations (Jakob and Siebesma 2003; Bretherton et al. 2004); however, it is clear from Fig. 7 that neither the bulk-plume entrainment rates nor the directly measured entrainment rates bear even a passing resemblance to a 1/z relationship. Furthermore, with ε and δ diagnosed for deep convection, we can put the proposed 1/z relationship through an even more basic test. Since the shallow convection reaches a height of 1.5 km and the deep convection reaches a height of 15 km, a 1/z scaling predicts that ε would be one-tenth as large in deep convection as it is in shallow convection. In the main part of the cloud layer in shallow convection, ε ranges from 2.2 to 2.8 km−1, so we would expect ε in the main cloud layer of deep convection to range from approximately 0.2 to 0.3 km−1. Instead, we find that ε ranges from 0.5 to 1.3 km−1 in deep convection. Furthermore, it is clear from the profile of ε in the right panel of Fig. 7 that ε does not even decrease monotonically with height in the main part of the cloud layer. In fact, it increases almost monotonically from approximately 5 to 11 km. Using the less restrictive activity operator, Fig. 8 shows that the fractional entrainment rate increases over the majority of the cloud layer in shallow convection. Therefore, fractional entrainment does not scale as 1/z.
Another approach has been to relate the entrainment and detrainment rates to b or the buoyancy gradient db/dz. One suggestion is that ε decreases with increasing buoyancy (Lin 1999). Another proposal is that entrainment increases with increasing db/dz (Bretherton and Smolarkiewicz 1989), an idea that has been used in convective parameterizations (Emanuel and Živković-Rothman 1999; von Salzen and McFarlane 2002). To test these ideas in the context of the shallow and deep simulations, Fig. 9 plots the profiles of (left) ε and (right) δ against the profiles of (top) b and (bottom) db/dz. The shallow data are shown in black, and the deep data are shown in blue. These plots disprove any simple relationships between ε and b, ε and db/dz, or δ and db/dz. The only glimmer of hope is the relationship between δ and b. In both the shallow and deep simulations, decreases in buoyancy lead to a roughly linear increase in fractional detrainment. Of course, the slope and intercept of the linear fits to the deep and shallow data are different, but this is to be expected: deep convective clouds are, typically, much wider than shallow clouds, so their smaller area-to-volume ratio leads to smaller fractional detrainment rates.
The linear relationship between buoyancy and δ can be illustrated quite nicely around the melting line in the simulation of deep convection. The domain-wide heating caused by the release of the latent heat of fusion is shown in the top curve of Fig. 10. The domain-averaged temperature profile crosses 273.16 K at a height of 4.6 km, denoted by a thin vertical line in the figure. (Another line at 4 km is also added to aid the eye in lining up features from the various curves.) As that top curve makes clear, the melting of snow and graupel below the melting line provides a strong diabatic cooling in a 1-km layer below 4.6 km. To see the effect of this cooling, we define ΔT as the domain-wide horizontally averaged temperature minus the running mean of that temperature over a 2-km interval. The second curve of Fig. 10 shows that this cooling produces a temperature excursion in ΔT of about −0.15 K at 4 km. The third curve shows the effect of this temperature excursion on the average buoyancy of active air. The buoyancy has a local maximum at 4.0 km, where ΔT is most negative, and a local minimum at 4.6 km, where ΔT has recovered. Between 4.0 and 4.6 km, the 0.006 m s−2 change in the average buoyancy of active air roughly agrees with the 0.2-K change in ΔT: 0.2 K/273 K × 9.81 m s−2 = 0.007 m s−2.
The variations in ε and δ at these heights, shown in the bottom two curves of Fig. 10, can be plausibly explained as a boom-and-bust cycle for the convecting clouds. As a cloud enters the vicinity of 4.0 km, the negative excursion in ΔT leads to an increase in the buoyancy of all parcels within that cloud. Parcels that were close to detraining experience a boost in buoyancy that temporarily spares them that fate, which reduces the overall fractional detrainment of the cloud. Meanwhile, since the whole cloud is experiencing a larger buoyant acceleration, it is able to activate surrounding air that would not have otherwise reached a vertical velocity of 1 m s−1. But when ΔT recovers around 4.6 km, much of that extra baggage—those parcels that would have detrained earlier or would not have entrained in the first place—detrains, causing a rise in δ. This sudden reduction in buoyancy means the cloud is unable to accelerate as much air to the 1 m s−1 threshold, which lowers ε. In other words, the boom goes bust.
9. Spatial and temporal distributions of e and d
Of course, there is only so much that can be accomplished by studying horizontally averaged profiles. Fortunately, direct measurement provides access to the full four-dimensional distribution of e and d. This paves the way to correlating e and d with the local properties of the flow to reveal the mechanisms responsible for entrainment and detrainment. To illustrate the level of detail that is possible, Fig. 11 displays the azimuthally averaged (top) active air, (middle) entrainment, and (bottom) detrainment during the release of a deep convective bubble simulated with an isotropic 100-m grid spacing. Each panel shows data averaged over 20 min, with averaging intervals ending (left) 40, (middle) 60, and (right) 80 min into the simulation.
Between 20 and 40 min, lateral entrainment takes place predominantly at radii between r = 1 and 2 km and heights between z = 1 and 5 km. There is also a smooth, striated pattern of entrainment at an altitude of 7.5–10 km that is suggestive of pileus caps. At 60 min, there are still heavy amounts of entrainment between radii of 1 and 2 km and extending now up to an altitude of 10 km. There is also a ring of entrainment in the outflow at a radius of 4 km and a height of about 12 or 13 km. Finally, at 80 min, the pattern of entrainment has settled into a mushroom shape, resembling something of a question mark in this r–z plot. The highest concentrations of entrainment are confined between r = 0 and 2 km and z = 1 and 7 km.
At both 40 and 60 min, the location of maximum detrainment is confined to the upper central core of the cloud. At 40 min, the detrainment maximum occurs at r = 1–2 km and a height of 5 km. At 60 min, the region of maximum detrainment has moved to r = 0–3 km and z = 12–15 km. At 80 min, the pattern has largely disintegrated, but the maximum detrainment still appears to be located in the upper central region of the cloud.
As shown analytically, the standard bulk-plume method is systematically biased toward diagnosing low values of entrainment and detrainment. In previous work, the bulk-plume analysis has relied on naturally occurring conserved tracers, such as total water in nonprecipitating convection. By introducing artificial tracers, the standard analysis can be extended to precipitating convection. The use of these tracers in large-eddy simulations of shallow and deep convection reveals more problems with the bulk-plume method: it produces unphysical negative values of entrainment and it produces estimates of entrainment that depend on the distribution of the tracer being used.
Direct measurement of entrainment and detrainment is proposed here as an alternative to the bulk-plume method of diagnosing these rates. Direct measurement diagnoses these quantities at the gridcell level by analyzing the local sources and sinks of convecting air. Comparing the results of direct measurement to the results of the bulk-plume method in both shallow and deep convection, the bulk-plume estimates are found to underestimate entrainment and detrainment by roughly a factor of 2.
It is not being proposed here that the directly measured ε and δ should be used in the bulk-plume equations. As we have seen, there is no one set of entrainment/detrainment profiles that will make the bulk-plume method work for all tracers. This deficiency has to do with the inability of the bulk-plume method to accurately estimate the mixing ratios of entraining and detraining air ϕe and ϕd. Fortunately, there are other models of convection, such as the stochastic parcel model of Romps and Kuang (2010b), that intrinsically simulate cloud heterogeneity and can therefore predict ϕd. Coupled with a parameterization of convective organization, such models might be able to accurately predict ϕe as well. For the validation of such models, it will be critical that we have accurate measurements of ε and δ from large-eddy simulations.
The profiles of ε measured here for shallow and deep convection have already allowed us to evaluate some existing theories of entrainment. For example, fractional entrainment is found to not scale like the inverse of height 1/z. Also, fractional entrainment scales neither like the cloud buoyancy b nor like the gradient of cloud buoyancy db/dz. Instead, evidence is found that the fractional detrainment scales linearly with −b. Fortunately, direct measurement raises the possibility that future work could tease out the properties of the local flow that lead to entrainment.
Thanks are due to Zhiming Kuang and Chris Bretherton for helpful conversations, and to two anonymous reviewers for suggestions that improved the manuscript. The simulations were performed on the Odyssey cluster supported by Harvard’s FAS Research Computing Group.
Implementation of Direct Measurement
In a numerical model, “active air” is not a prognostic field that is explicitly advected around the domain. Instead, the activity operator 𝒜 is applied at the end of every time step to determine which grid cells are active. To say that “active air is being advected” is really to say that the quantities upon which 𝒜 depends are being advected along with the flow. For example, we could define active air such that some mixing ratio q is greater than some threshold, say, ½ kg kg−1. Let us suppose that we begin at t = 0 with grid cells that have either q = 0 kg kg−1 or q = 1 kg kg−1. With a monotonic and diffusive advection scheme, grid cells at later times will take values in the range of 0 to 1 kg kg−1; however, the activity will still only take binary values of 0 or 1 depending on the sign of q − ½. As the air with positive q − ½ advects along the grid, the active region will advect with it, albeit in a more halting fashion.
Since active air is not a quantity that is moved around the domain by the model’s advection scheme, there is no explicit calculation of the fluxes in the model. Therefore, we must provide our own interpretation of this quantity. In a finite-volume model, we must decide how to interpolate 𝒜 onto the faces of grid cells. In the calculation of , we use the facial values of ρ and that are used for the advection of dry-air mass. For 𝒜, we use first-order upwind interpolation onto faces for the following reasons: it is monotonic, its small stencil provides the highest spatial resolution of e and d, and it facilitates the use of fractional time steps to be discussed shortly.
In Eqs. (18) and (19), the contributions to e and d at generally occur only at the precise moments when the boundary of active air passes over . In a numerical grid, the location of the boundary is only known to within the grid spacing, so the moment of entrainment or detrainment cannot be known as precisely. Therefore, we must add up the numerical approximation to over the full time it takes the boundary to pass through a grid cell of width Δx. In particular, the activity source at a grid cell is added up for the entire time that the grid cell is adjacent to a cell of the opposite activity (i.e., the period of adjacency).
To see what goes wrong if we do not sum over the period of adjacency, let us consider the one-dimensional advection of active air with a Courant number of ½ as depicted in Fig. A1. The definition of 𝒜 is not important for this gedanken experiment so long as the mixing ratios on which it depends are all being advected along with the air. Since 𝒜 takes only binary values, it shifts by one grid cell once every two time steps. Figure A1 shows the values of 𝒜 at five different times for three grid cells oriented along the direction of the flow. Here, shaded cells are active, and blank cells are inactive. On the right are depictions of the contributions to the activity source for the middle grid cell during each time step. All values have been normalized to a unit square cell, a unit time step, and a unit density. The numbers next to the arrows indicate the magnitude of the flux through the respective faces (i.e., ρu𝒜), and the numbers in the center of the middle cell denote the local tendency of active air [i.e., ∂/∂t(ρ𝒜)].
Let us now walk through this example one step at a time. At time t = 0, the grid cell just to the left of the three depicted here is active. During the first time step, the upwind values of 𝒜 on the middle cell’s two faces are both zero, so there is no contribution to the activity source in the middle grid cell. At t = 1, the mixing ratios of the left grid cell have changed sufficiently to qualify that grid cell as active. At this point, the middle cell is adjacent to the activity boundary, so we begin summing the activity source from this point on. During the second time step from t = 1 to t = 2, the upwind value of activity on the first face is 1, so the flux through that face is equal to one times the Courant number, or ½. Therefore, during the time step from t = 1 to t = 2, the activity source for the middle cell is −½. From t = 2 to t = 3, the flux through the first face is again equal to ½, contributing an additional −½ to the activity source. But at time t = 3, the middle cell qualifies as active, so the switch from inactive to active contributes 1 to the activity source. During the next time step, the upwind value of 𝒜 on both faces is 1, so there is no longer a convergence of active air into the middle cell. By time t = 5 (not shown), the right cell will qualify as active, and the middle cell will therefore no longer be adjacent to the activity boundary. At that point, we sum the activity sources for the middle cell and assign the result to either e or d. As we would expect for advection, the sum of sources is −½ − ½ + 1 = 0. Had we assigned the sources at every time step, then we would have concluded, erroneously, that the middle cell had detrained max[0, −(−½)] = ½ units of active air during the second time step and had entrained max(0, −½ + 1) = ½ units of active air during the third time step. Only by summing the activity source over the period of adjacency do we properly account for the finite grid spacing.
We must also account for the finite size of the time step. In general, cells will flip between active and inactive at some point in the middle of a time step. This complication is avoided if the Courant number is much less than one or, in the case of 1D advection, is equal to 1/n for some integer n. In the case depicted in Fig. A1, the Courant number was ½, so we did not run into any difficulty with regards to the finite time step. But let us consider the case depicted in Fig. A2 in which air with q = 0 kg kg−1 and q = 1 kg kg−1 is advected to the right with a Courant number of 0.7. Consider an activity operator defined as 𝒜 = ℋ(q − ½). At t = 0, the average q in the first cell is 0.2 kg kg−1, so it is inactive. Therefore, the upwind value of activity on the middle cell’s left face is zero, as shown by the flux numbers labeled “W/o averaging.” At the end of the first time step, the first cell has an average q = 0.9 kg kg−1 > ½ kg kg−1, so it is active. Therefore, the flux through the left face during the second time step is 0.7. During that second time step, the activity of the middle cell also switches to one, so the total activity source for the middle cell is 1 − 0.7 = 0.3. During subsequent time steps, the convergence of active air into the middle cell is zero, and, by t = 4, the right cell becomes active and the middle cell is no longer adjacent. Therefore, we would conclude that the middle cell entrained 0.3 units of air.
Since the active air is simply advecting, we should not have measured any entrainment or detrainment. So what has gone wrong here? When calculating the fluxes during the first time step, we have not accounted for the activity of the first cell during four-sevenths of the time step. By using the value of the first cell’s activity at t = 0 throughout the entire time step, we have underestimated the flux through that left face. Instead, we should use the average upwind value during the time step, which was . The average flux was, therefore, . Similarly, during the second time step, the average value of the middle cell’s activity was one-seventh, so the average flux on the second face was . By averaging activity during time steps when calculating the fluxes, we find that the activity source in the middle cell was −0.4 − 0.7 + 0.1 + 1 = 0.
To simplify the algorithm for averaging activity in a numerical model, we restrict attention to 𝒜 that can be defined as
for some set of n gridcell properties ξi. For example, if n = 2 with ξ1 = qc − 10−5 and ξ2 = w − 1, then air is active when it has a condensate mixing ratio greater than 10−5 kg kg−1 and a vertical velocity greater than 1 m s−1. Let us denote the value of ξi before the time step by a superscript 1 and the value afterward with a superscript 2. For a ξi that changes sign during the time step, the time at which it crosses zero is |ξi1|/|ξi2 − ξi1|, which can also be written as |ξi1|/(|ξi1| + |ξi2|).
Now, consider a grid cell that becomes active during a time step. This means that there are some negative ξi1 but no negative ξi2. The crossing time for a ξi that switches from negative to positive can be written as −ξi1/(|ξi1| + |ξi2|). Because of the logical “ands” in the definition of 𝒜, the cell becomes active when the last ξi becomes positive. Therefore, we want the largest crossing time for all the ξi that cross zero during the time step. Since the ξi that do not cross zero all have positive ξi1, the switch from inactive to active takes place at the maximum of −ξi1/(|ξi1| + |ξi2|) for all i. Similarly, if the cell flips from active to inactive, the switch takes place at the minimum of ξi1/(|ξi1| + |ξi2|) for all i such that ξi2 < 0. Therefore, 𝒜 may be summarized as
Corresponding author address: David M. Romps, Harvard University, 416 Geological Museum, 24 Oxford St., Cambridge, MA 02138. Email: firstname.lastname@example.org