## Abstract

The vertical velocities of convective clouds are of great practical interest because of their influence on many phenomena, including severe weather and stratospheric moistening. However, the magnitudes of forces giving rise to these vertical velocities are poorly understood, and the dominant balance is in dispute. Here, an algorithm is used to extract thousands of cloud thermals from a large-eddy simulation of deep and tropical maritime convection. Using a streamfunction to define natural boundaries for these thermals, the dominant balance in the vertical momentum equation is revealed. Cloud thermals rise with a nearly constant speed determined by their buoyancy and the standard drag law with a drag coefficient of 0.6. Contrary to suggestions that cloud thermals might be slippery, with a dominant balance between buoyancy and acceleration, cloud thermals are found here to be sticky, with a dominant balance between buoyancy and drag.

## 1. Introduction

The vertical velocities of cloud updrafts strongly affect aerosol activation rates (Abdul-Razzak et al. 1998), formation of hail (Danielsen et al. 1972), clear-air turbulence (Lane et al. 2012), aircraft hazard (Lane et al. 2003), lightning flash rates (Romps et al. 2014), tornado occurrence (Davies-Jones 1984), gravity wave generation (Fovell et al. 1992), the depth of convective overshooting (Wang 2007), and the convective moistening of the stratosphere (Grosvenor et al. 2007). Despite the importance of vertical velocities, the balance of forces giving rise to those motions is poorly understood. To make some headway on elucidating this balance of forces, we focus here on the quantum of moist convection: the cloud thermal. Our goal is to answer the following question: what is the dominant balance in the vertical momentum budget of mature cloud thermals?

One hypothesis is that the dominant balance in the vertical momentum budget of a mature cloud thermal is between buoyancy and acceleration . In this hypothesis, drag plays no significant role. Since cloud thermals take the form of a quasi-spherical vortex ring (Levine 1959), evidence in favor of this no-drag picture can be found in the experimental literature on vortex rings. Reynolds (1876) was the first to study the momentum budget of vortex rings in a laboratory setting, and he concluded that “these rings do move without any appreciable resistance.” This conclusion was bolstered by the laboratory work of Maxworthy (1974), who estimated the drag coefficient to be much less than one at a value of about . More recently, Sherwood et al. (2013) argue that cloud thermals should have no form drag, just like Hill’s vortex (Hill 1894). Because of this presumed lack of drag, Sherwood et al. (2013) refer to cloud thermals as “slippery.” Based on this assumption, the parcel model in that study neglects form drag and wave drag by setting the drag coefficient to zero.

The alternate (and opposite) hypothesis is that the dominant balance is between buoyancy and drag . We will refer to the thermals in this scenario as “sticky.” The evidence in favor of this hypothesis comes from recent cloud-resolving studies of moist convection. By conditionally sampling grid cells in a large-eddy simulation (LES) of shallow convection, de Roode et al. (2012) find that “the pressure gradient is the dominant term balancing the buoyancy,” and this conclusion is bolstered by the similar analysis of Wang and Zhang (2014). In fact, even the numerical results of Sherwood et al. (2013) exhibit a pressure perturbation acceleration comparable in magnitude to the buoyancy. Since form drag and wave drag manifest as a net pressure perturbation gradient within the thermal, these results support the hypothesis of sticky thermals. Here, we extend the studies of de Roode et al. (2012), Sherwood et al. (2013), and Wang and Zhang (2014) by studying thousands of snapshots of cloud thermals in an LES of deep convection with the goal of providing a definitive answer to the slippery-or-sticky question.

To put this issue in a mathematical context, let us imagine that we have some objective rule for identifying the region of space occupied by a thermal. Let be the subset of space identified as the thermal at time *t*. For a general function of space and time , the rate of change of the integral of *f* over can be written as

where is the two-dimensional boundary of , is the normal unit vector on pointing outside the thermal, and is a three-dimensional vector on and normal to that is defined to give the movement of . In particular, if **x** lies on , then lies on for infinitesimal . Let us define the entrainment velocity as , where is the fluid velocity. Like , is defined only on . Substituting into Eq. (1) and using Gauss’s theorem, we get

Setting in Eq. (2), using the momentum equation, and then subtracting a hydrostatic base state from the pressure and density, we get

where is the pressure perturbation (i.e., the deviation from the hydrostatic base state) and the buoyancy *b* equals , where is the density perturbation. The left-hand side is the rate of change of the thermal’s momentum. This is equal to the sum of the three terms on the right-hand side: the pressure perturbation gradient force, the buoyancy, and the entrainment and detrainment of momentum. We have neglected molecular viscosity here as its direct effect on the momentum budget of cloud thermals is wholly insignificant.

The slippery-thermal hypothesis is that the dominant balance in Eq. (3) is between terms 1 and 3, as indicated. The sticky-thermal hypothesis is that the dominant balance is between terms 2 and 3, as indicated. Evidence from large-eddy simulations indicates that the fourth term—acceleration by entrainment and detrainment—is far smaller than originally thought (Dawe and Austin 2011) and may be practically negligible (de Roode et al. 2012; Sherwood et al. 2013). Therefore, the question of whether mature cloud thermals are slippery or sticky boils down to a question of the magnitude of form drag and wave drag, as manifested in the pressure perturbation gradient force (term 2 above).

As a brief aside, it is worth emphasizing a simple, but important, point. The retarding force on a blob of fluid moving in the direction is equal to the integral of over the surface of the blob, where is the pressure perturbation, points normal to the surface outward, and both and are unit vectors. Writing this as the surface integral of and using Gauss’s theorem, we see that this can be written as the integral of the over the volume of the blob. Therefore, term 2 in Eq. (3) is the drag force, as claimed. One could argue that a piece of force should be thought of as contributing a virtual mass to an accelerating blob of air, giving rise to an effective inertia equal to times the mass of the blob in the case of a sphere (Batchelor 2005; Romps and Kuang 2010) and an infinite effective inertia in the case of an infinite, horizontal slab (giving rise to hydrostatic balance even if the slab is buoyant). This piece of the pressure perturbation force, which exists when the blob of fluid is accelerating even if it is momentarily motionless, cannot be described by a drag law of the form . Nevertheless, we will not attempt any semantic distinction between effective inertia and the other sources of drag. Instead, we will simply refer to the volume integral of as the drag; this terminology is particularly appropriate in this study since the cloud thermals are observed to rise without acceleration.

To assess the magnitude of the drag on cloud thermals, we study output from a high-resolution (100-m grid spacing) large-eddy simulation of radiative–convective equilibrium (RCE) over a 300-K ocean surface (section 2). An objective tracking algorithm is used to track cloud tops, and the thermals underneath those tops are identified from their azimuthally averaged streamfunctions (section 3). This produces 4852 snapshots of cloud thermals (of 715 unique cloud thermals), which reveal the dominant role for drag, supporting the sticky-thermal hypothesis (section 4). The results are then briefly summarized (section 5).

## 2. The large-eddy simulation

Simulations were run to RCE with Das Atmosphärische Modell (DAM; Romps 2008) on a square doubly periodic domain with a width of 32 km and a model top at 30 km. The time step adjusts automatically to satisfy the Courant–Friedrichs–Lewy (CFL) condition. The lower boundary was specified to be an ocean surface with a fixed temperature of 300 K, and surface fluxes were calculated using a bulk formula. Both shortwave and longwave radiation were calculated interactively using the Rapid Radiative Transfer Model (Clough et al. 2005; Iacono et al. 2008), and the top-of-atmosphere insolation was specified to be a constant diurnal average for the equator on 1 January.

A simulation with a 500-m horizontal spacing (and a vertical spacing of 50 m below 600 m and 500 m above 5 km) was started from an RCE sounding and run for 2 days to guarantee a fully convecting RCE state. The state at the end of that simulation was then interpolated to a grid with the same domain size (36 km × 36 km × 30 km) but with a 100-m horizontal spacing (and a vertical spacing of 50 m below 600 m and 100 m between 1100 m and 17 km). This was run for 20 h to ensure that the RCE state had adjusted to the new grid spacing. The last hour was run with a time step of 0.25 s (explained below) and with snapshots saved every minute (for a total of 60 snapshots).

This setup is similar to that used by Romps and Kuang (2010) and Romps (2011), and it produces a similar state of deep convective RCE. Figure 1 shows various profiles from the simulation, starting with the virtual potential temperature . Relative humidity is calculated with respect to liquid for and with respect to ice for . The cloud fraction is defined as the fraction of area occupied by a mass fraction of cloud liquid plus cloud ice greater than . The total convective mass flux is calculated from the three-dimensional snapshots as the average of , in which wherever vertical velocity exceeds 1 m s^{−1} and the mass fraction of cloud liquid plus cloud ice exceeds , and elsewhere. The convective mass flux of “good” thermals is discussed in the next section.

## 3. Tracking thermals

For each three-dimensional snapshot, cloud tops are identified as a grid point at the local top (highest point within ±500 m in *x* and *y*) of ascending cloudy regions. Cloud tops are then connected in time to build up sequences of cloud-top trajectories. Although this procedure may be straightforward to implement by eye, we do not attempt that here. Instead, we have built a set of objective criteria and then written a piece of software to automatically identify cloud tops and connect them in time. This algorithm, described in the appendix, identifies 4852 cloud tops (715 unique cloud sequences) from the 1 h of LES output. On average, a cloud sequence follows a cloud top for 6 min. For each cloud top, an azimuthally averaged map is made for each prognostic variable around a vertical axis that passes through the cloud top. These two-dimensional *r*–*z* maps, along with information on how they are connected in time, constitute the data used in this study.

Of course, the cloud top is just the tip of the iceberg, so to speak. Underneath lies the cloud thermal that drives the cloud-top ascent. The procedure for identifying the thermal’s outline is depicted in Fig. 2. First, the cloud-top vertical velocity is calculated using adjacent cloud-top positions in the cloud sequence. Next, a streamfunction is defined in the *r*–*z* coordinates using the azimuthal average of ,

Since the anelastic continuity equation is

where *u* is the radial component of velocity, is related to *u* and *w* by

where is the mass flux (kg m^{−1} s^{−1}) in the *z* direction (in the thermal’s reference frame) and is the mass flux (kg m^{−1} s^{−1}) in the *r* direction. Finally, the streamline is chosen as the boundary of the thermal. This choice is motivated by two facts: 1) there is no net mass flux across that boundary (net mass flux is zero across a streamline), and 2) the net vertical mass flux within the contour at each height is zero in the cloud top’s reference frame (or, in other words, the mean vertical velocity in the Earth’s reference frame is equal to ). This streamline definition of the boundary may be thought of as a generalization of the spherical approach taken by Sherwood et al. (2013). For Hill’s vortex, which is spherical, both approaches would identify the correct boundary, but our streamline approach has the advantage of being able to identify the boundaries of vortex rings that are not spherical.

We refer to the region that is contiguous with the cloud top as the thermal’s mask. The mask is deemed “bad” and not used in subsequent analyses if the top of the mask slopes upward with increasing *r*, if the mask forks with increasing radius (i.e., if, at constant radius, flips sign with height more than two times), if the mask reaches higher than the cloud top, if the mask reaches to a radius greater than 2 km, or if the mask extends more than 3 km below the cloud top. Of the 4852 unique cloud tops that are elements of cloud sequences, one quarter of those (1224 cloud tops) have well-behaved (i.e., good) streamfunction masks by these criteria.

The right panel in Fig. 1 shows the mass flux of the good thermals. For a height interval , this is calculated as the sum of vertical momenta (of all thermals instantaneously observed in the height interval) divided by , the number of snapshots, and the area of the simulation domain. Given the fact that we are tracking cloud thermals only when their tops are rising through clear air, and using several other restrictive criteria to define good thermals to simplify the analysis, it is no surprise that the mass flux of good thermals is a small fraction of the total convective mass flux (to be precise, the height-integrated good-thermal mass flux is 3% of the height-integrated total mass flux). Nevertheless, the profile of good-thermal mass flux shows that the cloud-tracking algorithm is sampling convection throughout the depth of the troposphere (aside from an intentional avoidance of very shallow thermals).

## 4. Results

First, we perform a sanity check of the data by comparing the mean Eulerian *w* inside a thermal with its Lagrangian . Whenever we refer to the “mean” of some quantity *X* in a thermal, this is calculated as

where the factor *r* in the two integrals makes this equal to the volume-weighted average in three dimensions. By definition of the thermal mask, the mean *w* should be nearly identical to . The only deviations should come from the coarseness of the grid (the uncertainty in the location of the streamline is roughly equal to the grid spacing) and the small variations in over the thermal. Figure 3 plots the mean *w* against for all thermals with good masks. The points fall on the 1-to-1 line with a coefficient of determination of 0.98.

With confidence in the procedure for identifying thermals, we can now proceed to examine the properties of these thermals in detail. The following subsections discuss the thermal structure, the momentum budget, and the internal circulation of the cloud thermals.

### a. Thermal structure

Because of the discrete nature of the 100-m grid, 138 of the 1224 thermals with good masks have a volume exactly equal to the median volume. These thermals have a volume equal to a sphere with a diameter of 560 m. Figure 4 shows various mean properties of these median-volumed thermals (cloud condensate, vertical velocity, streamfunction, pressure perturbation, pressure perturbation gradient acceleration, and buoyancy).

The cloud condensate is located just below the top of the cloud top, as expected. The vertical velocity has a maximum underneath the cloud top and weakly negative values at the periphery of the thermal, consistent with previous large-eddy simulations (Heus and Jonker 2008; Glenn and Krueger 2014). For the streamfunction, only positive values are plotted here; the negative values outside the thermal would saturate the color bar if plotted. By the procedure outlined in section 3, the thermal is identified as the region with a positive streamfunction. Therefore, the shape of the mean thermal can be identified here as the region with any red or pink color.

The pressure perturbation has a dipole structure with high pressure at and above the cloud top and a larger region of low pressure underneath the cloud top. If the perturbations of pressure and density were in hydrostatic balance, then the pressure perturbation gradient acceleration and the buoyancy would be mirror images of each other. Of course, cloud thermals are not hydrostatic, so this is not the case. The core of the thermal is positively buoyant (with the most buoyant air just below the cloud top), and that buoyant core is surrounded by an egg-shaped shell of negative buoyancy. In contrast, the pressure perturbation gradient acceleration resembles a sandwich with a wide and flat region of downward acceleration capped above and below by upward acceleration.

Returning now to all 1224 thermals with good masks, Fig. 5 plots the thermal-mean virtual temperature anomaly against the thermal-mean temperature anomaly. Here, virtual temperature is defined to include the loading by condensates. The virtual temperature anomaly is related to the buoyancy *b*, the density perturbation , the temperature anomaly , the water vapor mass fraction anomaly , and the condensate mass fraction anomaly by

The average temperature anomaly is 0.6 K and the average virtual temperature anomaly (accounting for the virtual effect of both water vapor and condensates) is 0.4 K. For thermals with positive buoyancy, the relationship is approximated fairly well by . Many in situ observations of cloud temperatures are reported as an anomaly of unloaded virtual temperature . To compare with those studies, we note that the mean unloaded virtual temperature anomaly (not shown in Fig. 5) is 0.8 K. Therefore, the mean buoyancy of thermals (+0.4 K) stems from the following combination of temperature anomaly, water vapor anomaly, and condensate anomaly:

These numbers match well with the buoyancies calculated from in situ observations of temperature and condensate loading within maritime tropical updrafts. Mean unloaded virtual temperature anomalies for convective updrafts have been reported as 0.55, 0.5–0.7, and 0.4–0.6 K over the western Pacific Ocean during the Taiwan Area Mesoscale Experiment (TAMEX; Jorgensen and LeMone 1989), the Equatorial Mesoscale Experiment (EMEX; Lucas et al. 1994), and the Tropical Ocean and Global Atmosphere Coupled Ocean–Atmosphere Response Experiment (TOGA COARE; Wei et al. 1998), respectively. Mean buoyancies (reported as loaded virtual temperature anomalies) in maritime updrafts have been found to be approximately 0.5 K during the GARP Atlantic Tropical Experiment (GATE; Lawson and Cooper 1990) and 0.2–0.3 K during TOGA COARE (Wei et al. 1998).

### b. Momentum budget

Before we attempt to diagnose the momentum budget, it is useful to first look at the cloud-top trajectories. For each cloud sequence, the time series of cloud-top displacement is plotted as a thin black line in Fig. 6. When the streamfunction mask is bad, the trajectory is plotted with the cloud-top heights at those times omitted; this choice is made to be as consistent as possible with the parcel trajectories described below, which can only use data from good masks. Although they are difficult to see on this busy plot, the trajectories resemble a splay of straight lines. The average of all of these time series is shown as the thick black line, which is plotted up to the time when the number of thin black lines drops below 5. It, too, is a straight line, indicating the absence of any significant acceleration or deceleration of these thermals. Recall that none of the conditions used to select cloud tops had anything to do with the stage of the cloud’s life cycle; therefore, the linearity of these trajectories is truly remarkable. This linearity strongly suggests a balance of forces, which would cause the thermals to ascend with a terminal velocity.

To check this balance, we calculate theoretical time series of cloud-top heights using the thermodynamic variables inside the thermal masks. In particular, we integrate Eq. (3) twice to give height, but using only the second term on the right-hand side (buoyancy only; blue lines), or only the first term on the right-hand side (pressure perturbation gradient acceleration only; red lines), or the sum of the first two terms on the right-hand side (buoyancy plus pressure perturbation gradient acceleration; purple lines). The thin blue, red, and purple lines show the trajectories calculated for individual cloud sequences. The thick blue, red, and purple lines show the respective means up until when the number of thin lines drops below 5. Clearly, buoyancy alone (blue lines) would cause the thermals to accelerate upward much more quickly than observed. Similarly, the pressure perturbation alone (red lines) would cause the thermals to rapidly decelerate and descend, which is in stark contrast to the observed ascent. But, by combining these two terms (purple lines), we closely replicate the observed ascent. This demonstrates that the dominant balance in the momentum budget is between buoyancy and drag; cloud thermals are sticky. Note that this success is achieved despite our neglect of the entrainment–detrainment term. This agrees with previous findings that the entrainment drag is weak (Dawe and Austin 2011; de Roode et al. 2012; Sherwood et al. 2013).

To demonstrate this balance in another way, Fig. 7 plots the thermal-mean pressure perturbation gradient acceleration against the thermal-mean buoyancy. Here, the mean *b* values from the 1224 thermals with good masks are divided into 20 quantiles, each containing 56 thermal snapshots. For each of those quantiles, a circle is plotted at the average *b* and average . To good approximation, these points fall on the 1-to-1 line, indicating that the buoyancy is largely canceled by the pressure perturbation gradient force; in other words, buoyancy and drag cancel. The drag deviates from the 1-to-1 line most significantly for values of *b* near zero. This is caused by a bias in our sampling: we are collecting only rising thermals, and doing so regardless of their buoyancy. Thermals will have a drag acting opposite their motion, which means that all of these rising thermals will have a positive even though some of their buoyancies may be negative.

How does this observed drag relate to the usual drag law? To find out, we plot the pressure perturbation drag against , where the thermal’s area *A* is calculated from its volume *V* assuming a spherical shape; that is, , where (Fig. 8). As in Fig. 7, the scatter is reduced by averaging over 20 quantiles of . The resulting data fall close to a straight line with a slope of 0.6. This tells us that the standard drag law applies,

with a drag coefficient of . This value of is not far from the value measured for solid spheres moving through a fluid. The highest Reynolds number (Re) for which the drag on a solid sphere has been measured is 10^{7}. At that Reynolds number, the drag coefficient for a sphere is 0.2 and rising with Re (Blevins 1992). Bear in mind that a cloud thermal is rising with a Reynolds number on the order of 10^{8}–10^{9}, although the effective Reynolds number of this LES is significantly less. And, unlike the typical sphere-in-flow laboratory experiment, the thermals simulated here are rising through a stratified fluid, which leads to wave drag. Therefore, the value of is certainly plausible. Note that 0.6 is larger than the value of 0.2 assumed in the parcel model of Romps and Kuang (2010). It is, of course, also larger than the value of zero used by Sherwood et al. (2013) and others (e.g., Nie and Kuang 2012).

It is important to note that care must be taken to extract the correct pressure field from the large-eddy simulation. DAM uses a split-time scheme in which acoustic modes are integrated with a small, but inexpensive, acoustic time step, and the rest of the dynamics is integrated with a large time step that is advanced using a third-order Runge–Kutta scheme (Romps 2008). The pressure field that is saved to the output is the pressure at the end of the final Runge–Kutta step. For a large time step that is too large, there are many intervening acoustic steps and the pressure field at the end of the large time step may not be representative of the mean pressure field experienced during the large time step. This difference manifests itself in a dependence of the inferred on the length of the large time step. This is shown in Fig. 9, where the calculated from the LES has not converged until a time step of ≲1 s is used. As a result, all of the results reported here are from the simulation using a 0.25-s time step.

### c. Internal circulation

The sticky behavior found here (i.e., ) is quite dissimilar from Hill’s vortex, which experiences zero drag. We can see the source of this difference by comparing a composite thermal for the LES side by side with Hill’s vortex. For this purpose, we average together the *r*–*z* maps of all thermals with good masks that have volumes in between the 90th and 95th percentiles; this leads to an averaging of 65 thermals with a mean volume equal to a sphere with a diameter of 920 m. These high percentiles of volume are chosen to give a relatively high-resolution view of the circulation. Figure 10 plots the mean of these thermals on the left-hand side next to Hill’s vortex (scaled to the same size and velocity) on the right-hand side. Contours represent streamlines ranging from −50 to 50 Mg s^{−1} in increments of 10 Mg s^{−1}. The contours are the largest closed streamlines plus the streamline at . Colors represent the perturbation pressure; the same scale is used for both panels.

The streamfunction for the Hill’s vortex is

where *a* is the radius of the contour. This definition of the streamfunction differs from the most commonly reported versions by a factor of to ensure that this satisfies Eqs. (4) and (5). The pressure perturbation for Hill’s vortex is

The pointy noses on the contours at (on the top of the cloud thermal and on the bottom of Hill’s vortex) are artifacts of plotting contours on this coarse grid.

Note that the circulations are quite similar between the two. The contours in the cloud thermal are expected to be a bit washed out as a result of having averaged over 65 different thermals, all with slightly different shapes and sizes. Therefore, we should not read too much into the differences in the contours between the mean cloud thermal and Hill’s vortex; to the contrary, the agreement between the two circulations is quite remarkable.

With its broader features, the pressure field is less affected by the averaging. Note that the thermal has an obvious fore–aft pressure gradient, which Hill’s vortex does not have. This dipole distribution of pressure in the cloud thermal is real and visible in individual cloud thermals. It differs markedly from the tripole distribution of pressure—with fore–aft symmetry—in Hill’s vortex. It is this lack of fore–aft symmetry that causes the drag on the cloud thermals.

What breaks the fore–aft symmetry in cloud thermals? One could certainly point a finger at turbulence: cloud thermals trail a turbulent, cloudy wake, which is distinctly different in character from the laminar, clear air above it. But, even in the absence of turbulence, it is impossible for a buoyant cloud thermal to have both a symmetric circulation and a symmetric pressure distribution like Hill’s vortex. The reason is buoyancy. The vertical acceleration is given by , where is the pressure perturbation and *b* is the buoyancy. For both *w* and to have fore–aft symmetry, *b* would need to have fore–aft antisymmetry, but this is not possible for a cloud thermal with net positive buoyancy.

Finally, let us compare the ratio of cloud-top speed to cloud-core speed for cloud thermals and Hill’s vortex. Figure 11 plots against for all 1224 thermals with good masks. Here, is calculated as the maximum *w* within the thermal at . The best-fit slope is calculated using total least squares regression with a zero intercept. We see that . This will be a small underestimate of the maximum velocity in the thermal as a result of the maximum sometimes not being directly underneath the top of the cloud and, therefore, getting averaged down in the azimuthal average. We can conclude that a rough rule of thumb is that the core velocity is twice the cloud-top velocity, which is not far from the ratio of 2.5 in Hill’s vortex.

## 5. Summary

In summary, we have run an hour-long large-eddy simulation of maritime, tropical, deep convection and have automatically tracked 4852 cloud tops in that simulation. Azimuthally averaging around a vertical axis through the cloud top, we use the zero streamline in the cloud top’s reference frame to define the boundary of the cloud thermal (Fig. 2). These cloud thermals have unloaded virtual temperature anomalies and net buoyancies (Fig. 5) that are consistent with in situ measurements of convection over tropical oceans. The cloud tops are found to ascend with nearly constant vertical velocity, suggestive of a terminal rise speed. This conclusion is bolstered by the finding that the dominant balance in the thermals’ vertical momentum equation is between buoyancy and drag (Figs. 6 and 7). The standard drag law predicts this drag force well using a drag coefficient of (Fig. 8). The zero used by Sherwood et al. (2013), who argued that cloud thermals experience no drag and, therefore, should be thought of as slippery, is inconsistent with these results. Instead, cloud thermals are sticky, ascending with a nearly constant terminal rise speed set by the balance between buoyancy and drag. A comparison of cloud thermals to Hill’s vortex shows that Hill’s vortex is drag free thanks to its fore–aft symmetry—a symmetry that cloud thermals do not have (Fig. 10).

## Acknowledgments

This work was supported by the U.S. Department of Energy’s Atmospheric System Research, an Office of Science, Office of Biological and Environmental Research program under Contract DE-AC02-05CH11231, and by a grant from the Undergraduate Research Apprentice Program (URAP) at the University of California, Berkeley.

### APPENDIX

#### Cloud-Tracking Algorithm

The first step is to take the three-dimensional snapshots, which are defined on the stretched vertical grid that is native to the simulations, and linearly interpolate them to a uniform 100-m vertical spacing. The stretched grid of the simulations has a spacing of 50 m below 600 m, a spacing of 100 m above 1100 m, and a smoothly transitioned spacing in between. The data on this stretched grid are interpolated to a uniform 100-m grid that exactly coincides with the 100-m levels above 1 km, so as to leave undistorted the data in the majority of the troposphere.

The data are on an Arakawa C grid with density, pressure, temperature, and mass fractions at the center of the grid boxes and velocity components on the faces. Let us denote the center of a grid box at time slice *n* by a 4-tuple of integers , with *i*, *j*, and *k* representing the position in number of grid cells east, north, and up, respectively. Let us denote the centers of faces by fractional positions—for example, is the location of *u* on the eastern face, is the location of *υ* on the northern face, and is the location of *w* on the upper face. We define a grid box to be a “cloud-top candidate” if and only if the Eulerian vertical velocity *w* at is greater than 1 m s^{−1} and the hydrometeor mass fraction (consisting of cloud liquid, ice, rain, snow, and graupel) at is greater than 10^{−5} kg kg^{−1} and less than or equal to 10^{−5} kg kg^{−1} at . A cloud-top candidate is defined to be a “cloud top” if and only if there is no other cloud-top candidate at a larger *k* within a square column of 1-km width centered on . A “cloud” is a set of cloud tops defined such that two cloud tops belong to the same cloud if one is within a square column of width 1 km centered on the other. Note that, by definition of a cloud top, all cloud tops within a single cloud will have the same *k*.

For a given cloud, the “cloud peak” is defined as the cloud top that has the largest Eulerian vertical velocity *w* at . In an attempt to attain a finer precision on the cloud peak’s “true” height, we linearly interpolate the hydrometeor mass fractions between and and define the cloud peak’s height to be the height at which the hydrometeor mass fraction equals 10^{−5} kg kg^{−1}.

The next step is to track cloud peaks through time. For each cloud peak at time *n*, we consider a 1-km-wide square column centered on that cloud peak and search for higher cloud peaks at time that fall within that column. If there is more than one that satisfies this criterion, the lowest one (in height) is taken, as it is deemed more likely to be connected with the peak in the previous time step. Conversely, if there are multiple peaks whose square columns share a common peak in the subsequent step, the highest one is taken, using a similar reasoning. This process is then repeated for all adjacent pairs of time steps, and sets of four or more cloud peaks linked in this way are called a cloud sequence. In other words, a “cloud sequence” is defined as a set of cloud peaks such that

the number

*N*of cloud peaks in the set is greater than or equal to 4;for some integer

*m*, the set contains exactly one cloud peak for each time ;for each , the cloud peak in the set from time is the lowest cloud peak at that time that is located within a 1-km-wide square column centered on the cloud peak in the set from time

*n*;for each , the cloud peak in the set from time

*n*is the highest cloud peak at that time that is located within a 1-km-wide square column centered on the cloud peak in the set from time ; andthis is the largest such set containing these cloud peaks.

To eliminate outlier clouds that seem to be ascending too rapidly, two final criteria are imposed. We define a “cloud-top velocity” to be a cloud’s cloud-peak ascent speed averaged over a 3-min interval (i.e., total vertical displacement during that interval, as calculated using the interpolated heights, divided by 3 min). Next, we define the “max velocity” associated with a cloud top located at as the maximum Eulerian vertical velocity *w* from all points such that

is within a square column of width 1 km centered on the cloud top’s location ,

, and

every grid box between and inclusive has a hydrometeor mass fraction greater than 10

^{−5}kg kg^{−1}.

Then, we discard all cloud sequences that have either

a cloud-top velocity greater than the maximum of the cloud-core velocities associated with each of the cloud peaks during the associated 3-min interval or

a coefficient of variation of three adjacent vertical displacements that exceeds 0.5.

## REFERENCES

*J. Adv. Model. Earth Syst.,*

**6,**1–8, doi:.