## Abstract

A method is introduced for diagnosing a transilient matrix for moist convection. This transilient matrix quantifies the nonlocal transport of air by convective eddies: for every height *z*, it gives the distribution of starting heights *z*′ for the eddies that arrive at *z*. In a cloud-resolving simulation of deep convection, the transilient matrix shows that two-thirds of the subcloud air convecting into the free troposphere originates from within 100 m of the surface. This finding clarifies which initial height to use when calculating convective available potential energy from soundings of the tropical troposphere.

## 1. Introduction

A transilient matrix for a convecting fluid describes the vertical transport of mass by eddies. Each column of the matrix corresponds to an initial height and each row corresponds to a final height: the element *b _{ij}* describes the amount of mass transported from

*z*to

_{j}*z*. If this matrix can be diagnosed for a convecting fluid, it can provide valuable information on the transilient (i.e., nonlocal) transport by eddies.

_{i}### a. Methods of diagnosing

A transilient matrix was first defined for convection by Stull (1984) [see also the review by Stull (1993)]. As originally defined, the element *b _{ij}* is proportional to the mass of air transported from

*z*at

_{j}*t*= 0 to

*z*at

_{i}*t*= Δ

*t*(per infinitesimal intervals around

*z*and

_{i}*z*). As such, is a function of the time interval Δ

_{j}*t*. In a numerical simulation,

*b*is easily diagnosed using a “set-and-go” method with tracers: initialize a horizontally uniform tracer

_{ij}*q*at

^{j}*t*= 0 in the layer containing

*z*, run the simulation until

_{j}*t*= Δ

*t*while advecting the tracer passively, and then set

*b*proportional to the average concentration of

_{ij}*q*in the layer containing

^{j}*z*. By using a tracer for every level, the entire matrix can be diagnosed in one simulation.

_{i}Although the set-and-go method is straightforward, it is not well suited to discerning the origins and destinations of convective eddies. If Δ*t* is chosen to be comparable to or shorter than the lifetime of a convective eddy, then *b _{ij}* will have significant contributions from eddies that were only passing through

*z*at

_{j}*t*= 0 or passing through

*z*at

_{i}*t*= Δ

*t*. On the other hand, if Δ

*t*is chosen to be comparable to or larger than the lifetime of a convective eddy, then layers other than

*z*will become significantly polluted by

_{j}*q*, which will cause

^{j}*b*to have contributions from eddies arriving at

_{ij}*z*but coming from layers other than

_{i}*z*.

_{j}This paper introduces a way to diagnose that avoids these problems. The method uses radioactive tracers that are continuously injected at their respective heights and that decay away with some time scale *τ*; we may refer to this method as “inject-and-decay.” The steady-state profiles of the tracers are arranged in a matrix that is inverted to obtain the transilient matrix. This transilient matrix can then be used to answer several questions about convection. In a simulation of moist convection, it can tell us where in the boundary layer cumuli originate.

### b. Importance of cloud origin

The calculation of convective available potential energy (CAPE) often depends sensitively on the height of origin used for the lifted parcel (Williams and Renno 1993; Mapes and Houze 1992; Smith 1997; Kingsmill and Houze 1999). This sensitivity is caused by the fact that the equivalent potential temperature *θ _{e}* need not be well mixed within the mixed layer even though the virtual potential temperature

*θ*is. To deal with the ambiguity of initial height, a common approach has been to calculate CAPE using a parcel with the average mixed-layer potential temperature and specific humidity (e.g., Bluestein and Jain 1985; Kingsmill and Houze 1999). If, however, cumuli draw their mass preferentially from a surface layer where

_{υ}*θ*is higher than the mixed-layer average, then these approaches will underestimate the “true” CAPE.

_{e}The details of cloud origin are also relevant to theories of the triggering of moist convection. Current lines of research have begun to parameterize the dry thermals that become convecting clouds (e.g., Cheinet 2003; Siebesma et al. 2007; Neggers et al. 2009). An important input into these models is how much air the thermals entrain as they rise from the surface to the lifting condensation level (LCL). The thermal’s entrainment rate affects the probability that a newborn cumulus will have the kinetic energy and buoyancy to reach its level of free convection.

### c. Evidence on cloud origin

It is well known that cumuli are formed by thermals fed from plumes rising off the surface. In a cumulus-topped boundary layer, glider pilots, hang gliders, paragliders, and pilots of model gliders are often able to ride an individual updraft from near the surface to the base of a newborn cumulus (Wallington 1961; Reichmann 1978; Pagen 1995; Rennó and Williams 1995). What is less certain is what fraction of a newborn cumulus originates from the surface as opposed to having been entrained by the thermal from other heights in the subcloud layer. This leads us to our main research question: what fraction of an average newborn cumulus is composed of air transported to the cloud base directly from a surface layer?

In observations (Crum et al. 1987; Raymond 1995) and in cloud-resolving models (CRMs) (Kuang and Bretherton 2006; Romps and Kuang 2010a), the thermodynamic properties of thermals and cumuli at the top of the boundary layer are very similar to the thermodynamic properties of surface air. The most likely explanation for this is that thermals transport surface air almost undiluted to the top of the boundary layer. An alternative explanation, however, is that thermals leave the surface with higher-than-average heat and moisture and then entrain enough air during their journey to reduce their thermodynamic properties to the surface-air mean. Another line of evidence—comparing observed cloud-base heights to the theoretical LCL of surface air—has produced mixed results. In some studies, the cloud base was found to agree with the LCL of undiluted surface air (Wilde et al. 1985; Stull and Eloranta 1985), while another study found the cloud base to agree with the LCL for air taken equally from all heights in the subcloud layer (Craven et al. 2002).

The use of Lagrangian particles in cloud-resolving simulations has also yielded mixed results. In a 3D simulation of shallow trade-wind cumuli, Heus et al. (2008) studied the vertical distribution of particles half an hour prior to being observed in a cloud at 1000 m (third panel of their Fig. 6). Of the particles that were below the 500-m cloud base half an hour before, 30% were below 100 m; this suggests an important, but not dominant, role for near-surface air. In another study using Lagrangian particles, Lin (1999) defined the origination height for a particle as the height where the particle’s vertical velocity was smaller than 0.0125 m s^{−1} for eight consecutive minutes. In a 2D simulation of shallow trade-wind cumuli, it was found that, of the particles in convecting clouds that originate from below the cloud base, about 60% of them originate from below 100 m (Lin 1999, his Fig. 4d). But, when the same methodology was applied to a 2D simulation of deep convection, no preference was found for subcloud particles to originate from below 100 m (Lin 1999, his Fig. 7d). One possibility for this discrepancy is that the 0.0125 m s^{−1} threshold is too sensitive for deep convection, in which downdrafts and cold pools frequently disturb the boundary layer.

### d. Outline of paper

The theory of the transilient matrix and the inject-and-decay method for diagnosing it are introduced in section 2. In section 3, low-resolution simulations are used to assess various properties of the transilient matrix, including the effect of the advection scheme and the matrix’s invariance with respect to the time scale for tracer decay. Section 4 presents high-resolution simulations, which reveal that the majority of subcloud air in cumuli originates from within 100 m of the surface. The transilient matrix used to obtain these results is compared to the transilient matrix of Stull (1984) in section 5. Section 6 summarizes the results and discusses applications.

## 2. Theory

We wish to describe fluid flow with a transilient function *b*(*z*, *z*′) that gives the mean advective mass flux from *z*′ to *z*. This function is useful for two important reasons, among others. First, it reveals where large eddies begin and end, which is relevant for the construction of convective parameterizations. Second, it can predict the convective redistribution of tracers without the need for explicit simulation.

We begin by defining as the flux of mass from *z*′ to *z* due to individual large eddies, which is a nonnegative quantity. This has dimensions of mass per time, per horizontal area, per height interval around *z*′, and per height interval around *z* (kg m^{−4} s^{−1}). Here, an “eddy” is any fluid motion whose advective time scale is short compared to the time scales of interest. A “large” eddy is any eddy whose vertical extent is large enough that it acts nondiffusively on the tracers of interest; we could refer to this synonymously as a transilient eddy. We refer to all other eddies as “small” eddies. The effect of a small eddy on the relevant tracers is simply to diffuse the tracers in the vertical; we could also refer to this as a diffusive eddy. The vertical scale that distinguishes between large and small eddies is the vertical scale of features whose transport we wish to study. Later, when we discretize the transilient matrix, this vertical scale is implicitly chosen to be the discretization length.

In the moist atmosphere, we are often interested in the transport of chemical tracers whose sources and sinks act on time scales much longer than the advective time scales of moist convective updrafts, moist convective downdrafts, and eddies in clear-air turbulence. In addition, atmospheric tracers typically have mean profiles that vary vertically on length scales that are large compared to the depth of clear-air eddies, which are confined vertically by static stability. Therefore, we can identify moist convective updrafts and downdrafts as large (i.e., transilient) eddies and the motions in clear-air turbulence as small (i.e., diffusive) eddies.

Now, consider a passive scalar with a mixing ratio *q* in a convecting fluid. Let us assume that three conditions are satisfied:

the fractional area of large eddies is much less than one,

the large eddies begin with the average

*q*from their height of origin, andthe sources and sinks of

*q*act on time scales much longer than the large-eddy advective time scales.

As shown in the appendix, these conditions allow us to write the vertical convergence and advection of *q* in terms of :

where an overbar denotes an average over *x*, *y*, and *t*. The compensating mass flux *M* is defined as

because the first term on the right is the eddy flux from above *z* to below *z* and the second term is the eddy flux from below *z* to above *z*. The diffusivity *k* is the turbulent diffusivity effected by the small eddies. Equation (1) has a straightforward interpretation. The first term on the right-hand side describes the displacement of with the brought in by large eddies starting at *z*′ and ending at *z*. The second term describes the advection of *q* by compensating subsidence. The third term describes the diffusion of *q* by small eddies.

The transilient function is defined as

where *δ* is the Dirac delta function. The terms on the right-hand side describe the net source of mass at *z* due to four effects: individual large eddies arriving at *z* from *z*′, large eddies leaving *z* for all possible *z*″, compensating subsidence, and turbulent diffusion by small eddies. As shown in the appendix, the vertical fluxes of *q* are related to *b* and by

It is also shown in the appendix that *b* conserves mass. In particular, the net sink of mass at *z*′ (due to motions going to all possible *z* destinations) is zero:

This guarantees that the integral over *z* of the right-hand side of Eq. (4) is zero, just like the left-hand side. Also, the net source of mass at *z* (due to motions leaving from all possible *z*′ origins) is zero:

This guarantees that, if is replaced by a constant, the right-hand side of Eq. (4) is zero, just like the left-hand side.

Using a set of radioactive tracers in a cloud-resolving simulation, we can diagnose a discretized version of the transilient function (i.e., the transilient matrix). Consider a set of *n* radioactive tracers with mixing ratios . Each tracer will be given its own time-independent source *S ^{k}*(

*z*) and will decay with a time scale

*τ*. The governing equation for the

*k*th tracer is

Averaging over *x*, *y*, and time, this equation reduces to

where the average is denoted by an overbar. If the three previously listed conditions are satisfied, then we can use Eq. (4) to write this as

Let us now choose a vertical grid {*z _{i}*;

*i*= 1, … ,

*n*} and an associated set of grid spacings {Δ

*;*

_{i}*i*= 1, … ,

*n*}. Defining , , , and

*S*=

_{ik}*S*

*(*

^{k}*z*), we can approximate Eq. (7) as

_{i}Solving for *b _{ij}* gives

where is the matrix inverse of .

The and matrices have straightforward interpretations, as shown in Fig. 1. Row *i* of the matrix contains the mean tracer values at level *i*. Column *j* is the mean profile of tracer *j*. Roughly speaking, each row of tells us where air at a specific level originated, and each column tells us where air from a specific level is destined to go. More precisely, row *i* gives the sources (due to large eddies, small eddies, and compensating subsidence) at level *i* binned by origin. Similarly, column *j* gives the sinks at level *j* binned by destination.

This definition of the transilient matrix does not suffer from the “convective structure memory” (Ebert et al. 1989; Stull 1993; Sobel 1999; Larson 1999) that hampers the definition used by Stull (1984). If we define the matrix , then is the matrix operator that evolves a horizontally averaged profile of a conserved tracer forward in time by Δ*t*. In other words, the vector (i.e., discretized profile) of tracer values *q** _{i}*(

*t*+ Δ

*t*) is related to

*q*

*(*

_{i}*t*) by

If *q** _{i}* has sources

*Q*(

_{i}*t*), this generalizes to

where, to avoid clutter, the indices and summation are implicit. If **Q** is time independent, then this becomes

where is the matrix inverse of . See section 6 for further comparison of this transilient matrix and the one defined by Stull (1984).

## 3. Low-resolution runs

Simulations are performed using the fully compressible cloud-resolving model called Das Atmosphärische Modell (DAM) (Romps 2008). The finite-volume advection scheme uses the three-dimensional Uniformly Third-Order Polynomial Interpolation Algorithm (UTOPIA) (Leonard et al. 1993) combined with the three-dimensional, monotonic flux limiter of Thuburn (1996). All variables are stored as double precision, and all prognostic variables are conserved by the finite-volume scheme to within roundoff error. For microphysics, DAM uses the six-class Lin–Lord–Krueger scheme (Lin et al. 1983; Lord et al. 1984; Krueger et al. 1995) and, for radiation, it uses the fully interactive Rapid Radiative Transfer Model (RRTM) (Clough et al. 2005; Iacono et al. 2008). The top-of-atmosphere downwelling shortwave flux is set to 572 W m^{−2} incident at a zenith angle of 43°. This provides the average 1 January equatorial insolation using a zenith angle whose cosine is equal to its insolation-weighted average on that day. The computational domain is three-dimensional with a doubly periodic boundary condition in the horizontal. The surface fluxes are parameterized with a bulk-flux scheme. Since there is no mean advection in these simulations, a 5 m s^{−1} eastward velocity is added to the local wind speed when calculating the bulk fluxes in order to simulate the flux enhancement from large-scale flow. The sea surface temperature is set to 300 K.

Several simulations are performed, all of which use the setup just described. The simulations differ in only three respects: grid size, duration, and tracer decay time. The first simulation, designed as a proof of concept, has a cubic domain that is 32 km wide. The grid spacing is relatively coarse, with a 2-km horizontal grid spacing and 64 vertical levels with a spacing ranging from 75 m near the surface to 500 m in the upper troposphere. Higher-resolution simulations are described in section 4. Since there are 64 vertical levels in this simulation, there are 64 radioactive tracers, each injected at its respective level. Each tracer is given a decay time of *τ* = 24 h. The simulation is run for 1 month and the first 10 days are discarded as spinup time. The tracer profiles are averaged over the remaining days and plugged into Eq. (9) to calculate *b _{ij}*.

### a. Near-diagonal terms

In both the *q _{ij}* and

*b*matrices, the largest elements are clumped in a band around the diagonal. This band is shaped by small eddies, but also by gravity waves and numerical diffusion. To check that our main results are insensitive to the numerics, all simulations have been run twice: once advecting the tracers with a first-order upwind scheme, and once advecting the tracers using the third-order UTOPIA scheme with the Thuburn flux limiter. Only the advection scheme for the tracers differs between the two runs; physical variables are always advected using the UTOPIA/Thuburn scheme.

_{ij}We might guess that the diagonal bands in and are formed primarily by vertical diffusion. Ignoring large eddies and compensating subsidence, the effective 1D equation for a diffusive, radioactive tracer injected at *z*_{0} is

which can be written in terms of a transilient matrix:

The discretized version of the transilient term is

so the rows of have the form

The steady-state solution to Eq. (10) is

Given *z*_{0}, *s*, *ρ*, and *τ*, there is only one free parameter: the vertical diffusivity *k*. In fact, by comparing the and matrices diagnosed from the CRM to the theoretical predictions in Eqs. (11) and (12), we can read off the vertical diffusivity in the CRM.

Figure 2 compares columns of the matrix to Eq. (12) and rows of the matrix to Eq. (11). The solid lines in the top panels show the column of that is equal to the profile of the tracer injected at 6.4 km. The solid lines in the bottom panels show the corresponding row of the transilient matrix, which equals the sources of air at 6.4 km binned by origin. To make visible the very small oscillations in , the ordinate comprises two log axes with a linear axis in between. The results in the left (right) panels were generated with the first-order (third-order) advection scheme.

There are three ways to estimate the vertical diffusivity from these profiles. We can use the peak value of the mixing ratio *q*, the three middle elements of the corresponding *b _{i}_{•}*, or the slope of the log of

*q*:

For the first-order advection, the value of *k* calculated in these three ways ranges from 14 to 19 m^{2} s^{−1}. For the third-order advection, *k* = 7–9 m^{2} s^{−1}. These are within the wide range of estimated diffusivities for the troposphere (e.g., Table 1 in Wilson 2004).

The dashed lines in Fig. 2 plot Eqs. (11) and (12) with *k* chosen as the mean of the three estimates. With first-order advection, we see that the profile of closely matches the expected exponential decay and the row of closely resembles the diffusion operator. With third-order advection, the tracer decays more rapidly than exponential. A consequence of this nonexponential decay is the more complicated structure of the transilient matrix, which corresponds to an operator with higher-order derivatives in addition to second-order diffusion. We also see that is higher below the injection point than above, which is caused by compensating subsidence.

### b. Off-diagonal terms

The parts of the transilient matrix that are of most interest to us are the off-diagonal terms corresponding to large eddies (i.e., convective clouds). Figure 3 shows the base-10 logarithm of the elements of (kg m^{−4} s^{−1}). The abscissa corresponds to the origin of the eddies and the ordinate corresponds to the destination of the eddies. For example, we can read off that the mass flux in eddies originating between 1000 and 1001 m and ending between 7000 and 7001 m is about 10^{−10} kg m^{−2} s^{−1}.

The most striking feature of Fig. 3 is the leftmost column, which corresponds to convective eddies that begin in the first layer (0–75 m) and end between 2 and 16 km. This column contains the largest elements away from the diagonal. The relative statistical error in Fig. 3 is about 10% or less, so this feature is statistically robust. In fact, the flux from that first layer is so large that it dominates the transport out of the subcloud layer. Of course, this simulation uses a very coarse vertical grid, so this result is explored at much higher resolution in section 4.

### c. Independence of τ

As stated in section 1 and shown in the appendix, Eq. (9) gives the transilient matrix when *τ* of the tracers is large compared to the lifetime of eddies. In moist convection, the lifetimes of clouds are on the order of 1 h, so *b _{ij}* should be independent of

*τ*so long as

*τ*is much greater than an hour. When

*τ*is comparable to an hour, Eq. (9) should give a matrix that is dependent on

*τ*and no longer equal to the true transilient matrix.

To test this, a 2-month simulation is performed using *τ* = 48 h. Month-long simulations are also run using *τ* = 24, 12, 6, 3, 1.5, and 0.75 h. Figure 4 shows the fractional differences between the *b _{ij}* calculated from these runs and the 48-h run. We see that

*b*changes very little over the first few halvings of

_{ij}*τ*, but it diverges significantly for

*τ*= 1.5 and 0.75 h. This confirms the invariance of

*b*for

_{ij}*τ*much larger than 1 h and the divergence of

*b*from the true transilient matrix for

_{ij}*τ*comparable to 1 h.

## 4. High-resolution runs

The most striking result from the previous section is that the air convecting out of the boundary layer comes predominantly from the layer adjacent to the surface. We check this phenomenon now with a high-resolution simulation. This simulation uses a 200-m grid spacing in the horizontal and 175 levels in the vertical with a spacing ranging from 5 m near the surface to 200 m in the free troposphere. The domain is a 25.6-km-wide square with a model top at 30 km. To expedite the approach to equilibrium while keeping *τ* much greater than the eddy lifetimes, *τ* is set to 12 h. The simulation is run twice, advecting the radioactive tracers once with first-order advection and once with third-order advection. Recall that the physical variables (momentum, water, etc.) are always advected with third-order advection, so the convection is exactly the same in both runs. The simulation is run for at least one week with the matrix and its tendency obtained from an average over the last day.

### a. Height of origin

The transilient matrices obtained from this simulation are shown in Fig. 5. With both advection schemes, there is a column of high values along the left edge of the matrix. This band is spread out over 10–20 vertical levels, so it is very well resolved. The matrix obtained using third-order advection has a strip of negative values centered on the column corresponding to 200 m. Because the elements of are, by definition, nonnegative, the transilient matrix must be nonnegative away from the diagonal. Clearly, the matrix obtained from tracers advected with the third-order scheme violates this criterion. The matrix obtained with first-order advection does not suffer from this problem, which inspires more confidence. With either scheme, however, the band of high values along the leftmost edge of the matrix is robust, persistent, and well resolved.

The band of high values extends from 0 to 100 m on the abscissa. In other words, the function *b*(*z*, *z*′) returns large values for *z*′ < 100, which means that individual large eddies are pumping air into the free troposphere from below 100 m. To quantify this effect, we can calculate the fraction of air passing upward through the cloud base that originated from below 100 m. For the high-resolution simulation, the cloud base is at 567 m. Therefore, for each height *z* in the free troposphere, we calculate

which is the fraction of subcloud air detrained at *z* that originated from below *h*. This is a cumulative distribution as a function of *h*.

Figure 6 plots this function of *h* for a variety of height ranges. The negative values in *b _{ij}* for third-order advection show up as negative slopes in the cumulative distributions. With either advection scheme, the majority of air transported out of the subcloud layer originates from below 100 m. With first-order (third-order) advection, the fraction is 67% (81%). This result concurs qualitatively with the laboratory experiments of Sánchez et al. (1989), who found that incipient thermals have low entrainment rates.

### b. Check with purity tracers

Recall that three conditions were used to identify the **b** in Eq. (4) as the transilient matrix. The third condition—that the tracer sources and sinks have long time scales—is easily satisfied by choosing a *τ* that is sufficiently large. The other two—that the fractional area of eddies is much less than one and that eddies begin with the average tracers at their height of origin—are not under our control. In the subcloud layer, the circulation comprises dry thermals that occupy a substantial fraction of the domain, which violates the first condition and may also lead to violations of the second condition. It is incumbent upon us, then, to check the results with an independent method.

To verify the conclusion that air from 0 to 100 m dominates the subcloud air in cumuli, we use the purity tracer that was first introduced by Romps and Kuang (2010a). The purity tracer is initialized to a mixing ratio of 1 below some height *h* and to 0 above *h*. As eddies rise up past *h*, they carry with them an initial purity of 1 that is subsequently reduced by the entrainment of air from above *h*. The mixing ratio of the purity tracer tells us what fraction of the air in the eddy originated from below *h*.

To analyze steady-state convection, the purity tracer is reset to 1 below *h* at every time step. Above *h*, the purity is set to 0 well outside the eddies at every time step. In Romps and Kuang (2010b), a grid cell was defined as “well outside” if the surrounding 7 × 7 × 7 box of grid cells contained only noneddying air. Here, we want to err on the side of too small a purity, so we use a box that is 3 × 3 × 3 grid cells. This provides a buffer of one grid cell in any direction from an eddying grid cell; a smaller box (i.e., 1 × 1 × 1) would cause the purity tracer to be zeroed as an eddy advects from one grid cell to the next. Since there is no cloud in the subcloud layer and updraft velocities are small, a grid box is defined as eddying if its vertical velocity is positive. At any given moment, roughly half of the boundary layer satisfies this condition, but only those parcels that maintain a continuously positive velocity avoid having their purity set to 0.

We introduce two purity tracers into the high-resolution simulation: “low” and “high.” The low purity tracer is set to 1 below *h* = 100 m and 0 above. The high purity tracer is set to 1 below the cloud base (*h* = 567 m) and 0 above. For a parcel above the cloud base, the ratio of the low to high tracers tells us the fraction of subcloud air that originated from below 100 m.

Figure 7 shows the horizontally averaged low tracer divided by the horizontally averaged high tracer. Below 100 m, both tracers are set to 1, so their ratio is 1. As we move up in the boundary layer, the high purity is still 1, but the average low purity is equal to the fractional area of eddies times the average low purity of the eddies. We see that the average low purity dips to as low as 0.1 just below the cloud base.

Above the cloud base, however, the ratio rises to approximately 70% (~80%) for first-order (third-order) advection. This matches the result from the transilient matrix, which gave a ratio of 67% (81%). This confirms that the majority of subcloud air in cumuli originates from below 100 m. Compare this conclusion with the null hypothesis that new cumuli are composed of air taken equally from all heights below cloud base. If the null hypothesis were true, then the ratio of low to high purities above the cloud base would be 100/567 = 18%, which is denoted by the dashed line in Fig. 7.

### c. Implications of near-surface origin

One implication of a predominantly near-surface origin is that CAPE is more accurately calculated by lifting parcels from near the surface. In Fig. 8, CAPE is calculated by adiabatically lifting (without the latent heat of fusion) a parcel initialized to the average properties of air at the height on the abscissa. The left panel shows the results of this calculation using 4096 soundings from snapshots of the high-resolution simulation. The right panel shows the same calculation using 1934 soundings collected at the Atmospheric Radiation Measurement Program (ARM) site in Nauru. The observations exhibit much more variability than the simulation due to the presence of diurnal, mesoscale, and synoptic variability, but the qualitative picture is the same. For both the CRM and ARM data, CAPE decreases almost monotonically with the height of origin. Therefore, calculating CAPE for an average subcloud parcel underestimates the true CAPE available to cumuli, which mostly originate from the first 100 m.

What is so special about 100 m? In dry boundary layers or boundary layers topped by nonprecipitating clouds, a surface layer can often be identified in the bottom 10% of the boundary layer by larger mean values of *θ _{υ}*. Beneath precipitating cumulonimbi, however, the boundary layer is strongly affected by cold pools (e.g., Moeng et al. 2009), which introduce horizontal heterogeneity that complicates the interpretation of the mean

*θ*profile. Nevertheless, a surface layer can be identified from the largest percentiles of

_{υ}*θ*. This is shown in Fig. 9, which plots the 95th percentile of CAPE for parcels taken from the initial height on the abscissa and lifted to the cloud base. In this calculation, we use profiles from ARM soundings at Nauru and soundings from snapshots of the cloud-resolving simulation. Again, the observations exhibit much more variability, which causes the 95th percentile to be much larger for Nauru than for the CRM. Nevertheless, in both the observations and the model, the CAPE value decreases rapidly within the first 100 m. This identifies the first 100 m as a surface layer in which some parcels have enough buoyancy to reach their level of free convection.

_{υ}## 5. Comparison to Stull

To compare the transilient matrix introduced here with the one defined by Stull (1984, 1993), we use the thought experiment depicted in Fig. 10. For simplicity, we work in nondimensional units. Consider an incompressible fluid of unit density in a rectangular 2D domain. Let us conceptually divide the domain into five horizontal layers of unit depth. Imagine that fluid rises from the first layer to the fifth layer in a fraction *ε* ≪ 1 of the domain at a speed of 1 (i.e., one layer per unit time). The remaining 1 − *ε* of the domain sinks at a speed of *ε*/(1 − *ε*) ≈ *ε*. A transilient matrix can be diagnosed in two different ways: the set-and-go method of Stull (1984) and the inject-and-decay method described in section 2.

In Stull’s set-and-go method, five tracers are initialized to a value of one in their respective layers and zero elsewhere. Using the notation of section 2, the matrix is defined with columns equal the vertical profiles of the tracers. At time *t* = 0, is equal to the identity matrix,

Following convention, the matrix is shown upside-down, which makes the diagonal rise from bottom left to upper right. Advecting the tracers for time Δ*t*, the equation

is solved for by discretizing in time:

At Δ*t* = 1, the ascending column has risen by one layer while the rest of the fluid has descended a distance *ε*. Since the ascending branch covers a fraction *ε* of the domain and the descending branch covers a fraction ~1, is

In the first column, we see that the ascending branch has removed *ε* of fluid from the first layer and deposited it in the second layer. In the second column, we see two effects: the ascending branch has removed *ε* of fluid from the second layer and deposited it in the third layer, and the descending branch has removed *ε* of fluid from the second layer and deposited it in the first layer. Analogous explanations apply to the other columns.

Now, let us wait a time equal to the eddy turnover time, which is Δ*t* = 4 in this case. This gives us the following:

The 1s in the top row tell us that fluid has been transported to the fifth layer from each of the first four layers. This corresponds to the flushing out of the initial tracers in the ascending column. The 1s in the leftmost column tell us that fluid was transported from the first layer to each of the other layers; this corresponds to the ascending column filling itself with the tracer from the first layer.

There are two main problems with the matrices in Eqs. (13) and (14). First, neither matrix conveys the simple story of convective transport from the bottom to the top with compensating subsidence in between. The matrix in Eq. (13) is simply the diffusion operator: each row has the {… , 1, −2, 1, …} pattern of the discretized ∂^{2}/∂*z*^{2} operator. Therefore, it contains no information about the large eddies, which defeats the purpose of a transilient matrix. The matrix in Eq. (14), on the other hand, has nonzero elements corresponding to large eddies that are not present in the flow. The second problem is that of convective structure memory (Ebert et al. 1989). If the two matrices in Eqs. (13) and (14) were consistent, then would equal 4 times , but it does not. This leads to the question of which, if either, of these is the correct description of the fluid flow.

If we are concerned with the redistribution of tracers on time scales long compared to the eddy turnover time, the answer is that neither is correct. We get the correct matrix using the inject-and-decay method, whereby we continuously inject the five tracers into their respective levels and let the tracers decay on a time scale that is large compared to the eddy’s advective time scale (i.e., *τ* ≫ 4). For simplicity of illustration, let *τ* be small compared to 1/*ε*, which is the time it takes the compensating subsidence to move one unit of height; this assumption of *τ* ≪ 1/*ε* is not necessary, but it simplifies our thought experiment. In a steady state, the factor of *q _{ik}*/

*τ*−

*S*in Eq. (8) is equal to the source of tracer

_{ik}*k*at height

*i*due to the fluid motion. The ascending branch generates a sink

*ε*of tracer 1 in layer 1 and a source

*ε*of tracer 1 in layer 5. Since we have assumed that the decay time

*τ*is much smaller than the subsidence time scale 1/

*ε*, the surplus of tracer 1 decays in layer 5 before it is exported by subsidence. Therefore,

*q*

_{51}/

*τ*−

*S*

_{51}=

*ε*. Similarly, the export of tracer 1 from layer 1 by the ascending eddy is matched locally by

*q*

_{11}/

*τ*−

*S*

_{11}= −

*ε*. By the same token, the compensating subsidence, which also occurs at a rate of

*ε*, must be matched by for

*i*> 1. There is also nonzero decay of tracer 1 in the ascending column, but this is of magnitude

*ε*/

*τ*, and we have required that

*τ*≫ 4 (see condition 3 of section 2), so

*ε*/

*τ*is much smaller than the other elements of . With our simplifying assumption of

*τ*≪ 1/

*ε*, is equal to the identity matrix plus

*O*(

*ε*) terms. Therefore, the nondimensionalized version of Eq. (9) yields

This transilient matrix describes the flow accurately: the first column describes the transport of fluid from the bottom layer to the top layer by the ascending branch, and the other columns describe the compensating subsidence. Since this matrix contains all the information about the flow, it can be safely exponentiated to any time interval, large or small.

The relationship between this inject-and-decay transilient matrix and Stull’s set-and-go matrix bears some similarity to the relationship between the two approaches of constructing linear response functions of cumulus ensembles to large-scale temperature and moisture anomalies (Kuang 2010). Like the set-and-go transilient matrix, a matrix that describes the response of the atmosphere some time after an initial temperature or moisture perturbation suffers from complications due to the finite lifetime of convective eddies. On the other hand, introducing temperature or moisture sources, building a matrix of the resulting temperature and moisture anomalies, and then performing a matrix inversion, like the inject-and-decay method introduced here, is appropriate for flows with time scales longer than the convective time scale and results in a linear response functions matrix that can be exponentiated [see the appendix of Kuang (2010)].

## 6. Conclusions and discussion

We have introduced a transilient matrix that describes the rearrangement of mass by convective eddies in moist convection. It is shown how this matrix can be diagnosed using a set of radioactively decaying tracers in cloud-resolving simulations. The resulting transilient matrix tells us that the majority of air in newborn cumuli originates from within 100 m of the surface. This conclusion is confirmed independently using purity tracers.

This conclusion has important implications for how CAPE is calculated. There are many choices that must be made when calculating CAPE for a lifted parcel, including the entrainment rate, how condensates fall out, and the effective temperature of condensate freezing. One of the most significant choices to be made is the initial height from which to lift the parcel. The influence of the initial height is so large that operational forecasts distinguish CAPE calculations by this choice: surface-based CAPE (SBCAPE) uses a parcel from the surface, mean-layer CAPE (MLCAPE) initializes the parcel to the mean potential temperature and specific humidity in a layer typically 50–100 mb deep, and maximum unstable CAPE (MUCAPE) picks the initial height that gives the largest value of CAPE (Doswell and Rasmussen 1994; Bunkers et al. 2002). The results obtained here suggest that CAPE should be calculated using a parcel with the mean potential temperature and specific humidity of the air below 100 m.

For a study of boundary layer processes, a 200-m horizontal grid spacing is far from ideal. Unfortunately, this is near the limit of computational feasibility for simulations with the following requirements: a run time much greater than an hour, a domain sufficient to contain deep convection, and a passive tracer for every vertical layer. There are two obvious concerns about this coarse spacing, and they would have opposite effects. One concern is that the simulated boundary layer thermals may be wider, and therefore more protected from lateral entrainment, than real thermals. This would tend to overestimate the fraction of air originating from below 100 m. The other concern is that the 200-m horizontal spacing may imprint itself on the vertical scale of eddies, which would cause a spuriously large amount of air to be scooped up from above 100 m. This would tend to underestimate the fraction of air originating from below 100 m. Since both of these effects would naively scale with the grid spacing, and since the 2-km and 200-m simulations agree quite well, it appears that these effects are not significant.

Although the simulations presented here are in a steady state, there is no requirement of a steady state in the derivation of . All that is needed is for the time derivative in Eq. (9) to be averaged over a sufficiently large domain or a sufficiently long time to capture the ensemble effects of convection. This makes possible a wide range of applications. One such application is in a superparameterized general circulation model (GCM), in which a cloud-resolving model is embedded within each GCM grid cell. In this case, the CRM’s large-scale features (including tracer profiles) are time-dependent. By diagnosing in this time-dependent convection, and then exponentiating to the full GCM time step, we can faithfully model the effect of convection on the GCM’s passive tracers. If the rank of were small compared to the number of GCM tracers, this approach would be less computationally expensive than simulating all of the GCM’s tracers explicitly in the CRM.

Another potential application of the transilient matrix is as a way to improve the transport of chemical tracers by standard convective parameterizations. If the transilient matrix could be parameterized in terms of bulk convective mass fluxes, then convective parameterizations could be augmented to include the transilient transport of chemicals in a GCM. There is no guarantee that such a parameterization is feasible, but the possibility merits investigation.

## Acknowledgments

This research was supported by the Office of Biological and Environmental Research of the U.S. Department of Energy under Grant DE-FG02-08ER64556 as part of the Atmospheric Radiation Measurement Program, by the Director, Office of Science, of the U.S. Department of Energy under Contract DE-AC02-05CH11231, and by NSF Grant ATM-0754332.

### APPENDIX

#### Deriving the Transilient Matrix

For convenience, we will denote the four-vector (*x*, *y*, *z*, *t*) by *x ^{μ}*. Note that

*x*is a nonunique label for the parcel of air located at (

^{μ}*x*,

*y*,

*z*) at time

*t*. If

*x*is within a large eddy, then we define

^{μ}*z*(

_{i}*x*) as the initial height of parcel

^{μ}*x*when it enters the large eddy and

^{μ}*z*(

_{f}*x*) as the final height of the parcel when it exits the large eddy. We then define as

^{μ}Note that the double integral of is

To reduce clutter, integrals are assumed to be taken over their variables’ full range (i.e., from 0 to ∞ above) unless specified otherwise. When the integrals are taken over finite ranges, we get

Furthermore, if we pick *z* between *z*′ and *z*″, multiply the above expression by *ρw*(*x ^{μ}*), and integrate

*x*and

*y*over the domain area

*A*, then we get the net upward mass flux eddying from [

*z*

_{1},

*z*

_{1}+ Δ

*z*] to [

*z*

_{2},

*z*

_{2}+ Δ

*z*] through height

*z*at time

*t*. If we then integrate over time, we get the net mass transported from [

*z*

_{1},

*z*

_{1}+ Δ

*z*] to [

*z*

_{2},

*z*

_{2}+ Δ

*z*] by large eddies through height

*z*during that time period.

This leads naturally to a definition for . For *T* much larger than eddy turnover times, is defined as

We can use any value of *z* between *z*′ and *z*″ because a parcel eddying from *z*′ to *z*″ must pass through all intermediate heights. Even if the parcel doubles back through *z* on its convoluted path from *z*′ to *z*″, the fact that *w* will be negative during downward passes through *z* guarantees that the parcel’s net contribution to the integral will be its mass and no more. If we were to choose a *z* outside *z*′ and *z*″, then the integral on the right-hand side would be zero because parcels eddying from *z*′ to *z*″ will pass as many times downward through *z* as they will upward through *z*. Therefore, we can also write

We can use Eq. (A3) to see that

is the flux of mass into [*z*_{2}, *z*_{2} + Δ*z*] due to individual large eddies originating in [*z*_{1}, *z*_{1} + Δ*z*].

Let us define *q*(*z*″, *z*′; *z*) as the mass-weighted average of *q* in large eddies at *z* traveling from *z*′ to *z*″:

Denoting the average over *x*, *y*, and *t* by an overbar, is the net vertical mass flux of *q* at height *z*. Using Eq. (A2), this can be separated into an integral over large eddies and the region outside large eddies:

The last term in Eq. (A7) is integrated over the region of fluid that contains only compensating mass flux and small eddies. We can write the flux of *q* due to the compensating mass flux *M* (i.e., the average *ρw* in the region outside large eddies) as *M* times the horizontal average of *q* in the noneddying region. The compensating mass flux is given by

because the first term on the right is the eddy flux from above *z* to below *z* and the second term is the eddy flux from below *z* to above *z*. The small eddies act diffusively by assumption, so their effect can be described by a vertical diffusivity *k*(*z*) acting on the horizontal average of *q* in the noneddying region. Assuming that the fractional area of large eddies is small, we can approximate the noneddying *q* by the horizontal average of *q* over the whole domain. Therefore, we can write as

We will now make two approximations. First, we will assume that the sources and sinks of *q* act on time scales long compared to the large-eddy advective time scales. This implies that tracers are relatively unaffected by sources and sinks as they are transported by eddies, which allows us to approximate *q*(*z*″, *z*′; *z*) as the mean mixing ratio that parcels have as they enter a large eddy at *z*′. Second, we will assume that the average large eddy starting at *z*′ has the same *q* as the horizontal average at *z*′. These two approximations reduce this equation to

Taking minus the vertical derivative and using Eq. (A8), this becomes

In other words, the vertical convergence and advection of *q* into *z* (left side) is equal to the flux of large eddies to *z* times their *q* surplus (first term on right), plus the advection of *q* by the compensating subsidence (second term on right), plus the diffusive convergence of *q* (third term on right). We define the transilient function *b*(*z*, *z*′) as

The physical interpretation of this equation is discussed in section 2.

It follows immediately from Eq. (A10) that the integral of *b* over its first argument is zero. It can also be shown that the integral over its second argument is zero:

These two identities reflect the fact that *b* conserves mass. In particular, the net source of mass originating from *z*′ is zero,

and the net source of mass arriving at *z* is zero,

Since the action of *b* on can be rearranged as

Eq. (A9) can be written in terms of *b* as

To arrive at Eq. (A13), we required three conditions: the fractional area of large eddies is much less than one; the large eddies begin with the average *q* from their height of origin; and the sources and sinks of *q* act on time scales much longer than the large-eddy advective time scales.