## Abstract

Entrainment in cumulus convection remains ill understood and difficult to quantify. For instance, entrainment is widely believed to be a fundamentally turbulent process, even though Turner pointed out in 1957 that dry thermals entrain primarily because of buoyancy (via a dynamical constraint requiring an increase in radius *r*). Furthermore, entrainment has been postulated to obey a 1/*r* scaling, but this scaling has not been firmly established. Here, we study the classic case of dry thermals in a neutrally stratified environment using fully resolved direct numerical simulation. We combine this with a thermal tracking algorithm that defines a control volume for the thermal at each time, allowing us to directly measure entrainment. We vary the Reynolds number (Re) of our thermals between laminar (Re ≈ 600) and turbulent (Re ≈ 6000) regimes, finding only a 20% variation in entrainment rate *ε*, supporting the claim that turbulence is not necessary for entrainment. We also directly verify the postulated *ε* ~ 1/*r* scaling law.

## 1. Introduction

The rate at which cumulus clouds mix with their environment, or entrain, has long been known to be central to their dynamics (Simpson 1983a; Cotton 1975; Simpson et al. 1965; Stommel 1947). This led to a large number of studies, particularly in the early days, focused on the entrainment and dynamics of discrete, transient, convecting “thermals,” believed to be the fundamental unit of convection [see the review by Yano (2014), as well as references given below]. With the advent of large-scale numerical modeling and convective parameterization schemes, however, attention shifted to describing the average entrainment of an *ensemble* of convecting clouds, often conveniently modeled as one or more continuous, steady-state, entraining “plumes” (Yano 2014; de Rooy et al. 2013). While such plume convection schemes are ubiquitous in global climate models, they lie in tension with the wealth of evidence that cumulus clouds are actually composed of discrete thermals (Romps and Charn 2015; Sherwood et al. 2013; Heus et al. 2009; Damiani et al. 2006; Blyth et al. 2005; Zhao and Austin 2005; Carpenter et al. 1998; Miller et al. 1983; Saunders 1961; Malkus and Scorer 1955; Scorer and Ludlam 1953). Furthermore, single-plume schemes suffer from an “entrainment paradox” in which no optimal entrainment rate exists (e.g., Sherwood et al. 2013; Mapes and Neale 2011). This leads to uncertainties in the parameterization of entrainment, which turn out to be some of the largest contributors to uncertainties in climate sensitivity (Zhao 2014; Klocke et al. 2011; Murphy et al. 2004).

Given the uncertainties in plume convection schemes and their tenuous connection to cumulus phenomenology, it seems worthwhile to turn back to thermals as a basis for understanding clouds and building parameterization schemes, as suggested by Sherwood et al. (2013) [see also Morrison (2017) for a recent effort to reconcile the thermal and plume pictures]. As with plume models, however, entrainment rates are key. The parameter of interest is typically the fractional gross entrainment rate *ε*, defined to be the fraction of a parcel’s volume (or mass, assuming small horizontal variations in density) that it entrains per unit vertical distance traveled, in units of inverse meters. A long-standing, widely used assumption regarding *ε*, sometimes known as the *entrainment assumption*, is that for a thermal of radius *r*,

(Johari 1992; Turner 1986; Simpson 1983b; Simpson and Wiggert 1969; Turner 1962; Levine 1959; Morton et al. 1956, hereafter MTT56). The number *e* is a constant which we refer to here as the *entrainment efficiency*. The entrainment assumption is also widely applied in plume models, where it is sometimes formulated in terms of an “inflow velocity” (e.g., Turner 1986; MTT56), but this is equivalent to (1). Importantly, in virtually all these cases, entrainment is thought to arise from turbulent mixing at the interface between buoyant fluid and the environment.

The 1/*r* form of the entrainment assumption (1) is a plausible inference from basic dimensional analysis, and has been indirectly verified in laboratory experiments. In particular, “dry” thermals propagating through a neutrally stratified (neutral) ambient with no phase change or buoyancy source were studied by Richards (1961, hereafter R61), Scorer (1957, hereafter S57), and MTT56; these authors wrote down theories based on (1), and found their experiments in broad agreement with those theories.

This work on idealized dry thermals in a neutral ambient provides a foundation on which a thorough and fundamental understanding of cumulus thermals might be built, presumably in the hierarchical manner often adopted in climate modeling (Jeevanjee et al. 2017; Held 2005). Despite the simplicity of this case, however, questions remain. Perhaps most surprisingly, the well-cited work of Turner (1957) suggests that entrainment in thermals arises from a relationship between the thermal’s impulse *P*, circulation *K*, and radius *r*:

Here the circulation *K* is computed along a roughly semicircular circuit passing through the center of the thermal’s vortex ring and then returning through the ambient fluid. If the density anomaly is negligible along this circuit, then the circulation should be constant (Fohl 1967). This conclusion is supported by experiments (Zhao et al. 2013).

As for the impulse *P*, it can be regarded as a generalized momentum of the fluid which can only be changed by external forces, such as gravity, and which is unaffected by internal forces such as pressure or viscosity (e.g., Shivamoggi 2010; Batchelor 2000). In the presence of positive buoyancy and no other external forces, we must then have *dP*/*dt* > 0. Following Turner (1957), we can then combine this with constant *K* and (2) to deduce that

that is, *buoyant thermals must entrain*. (Similarly, for a vortex ring of descending, negatively buoyant fluid, both *dP*/*dt* < 0 and *K* < 0, so the radius also increases.) Although this argument implicitly neglects the effects of viscosity, it does not rely on turbulent fluctuations. Thus, it predicts that laminar and turbulent thermals should entrain in the same way. Physically, this is because a vortex ring of a specified size and circulation can only translate at a given velocity; if the velocity changes, then the radius must change. This implication is quite surprising in light of the vast literature attributing entrainment in thermals to turbulence, but so far has received little attention.

Another question is the veracity of the 1/*r* scaling in (1). The laboratory measurements of R61, S57, MTT56, and others are consistent with such a scaling but have not verified it directly, and previous numerical investigations seeking a direct confirmation have yielded mixed results (e.g., Hernandez-Deckers and Sherwood 2018; Dawe and Austin 2013; Stirling and Stratton 2012). Finally, S57 and others studying the neutrally stratified case claim that detrainment is negligible, though this is never quantified and a detrained wake is clearly visible in photographs of the experiments. Thus, there are still several open questions:

Does entrainment depend strongly or weakly on the Reynolds number (Re)?

Does directly measured entrainment indeed obey the 1/

*r*scaling in (1)?Is detrainment indeed negligible for dry thermals?

To answer these questions, we perform fully resolved direct numerical simulation (DNS) of three-dimensional dry thermals rising through a neutral environment, closely analogous to the laboratory experiments of S57. While thermals have certainly been simulated before, using LES (e.g., Morrison and Peters 2018; Morrison 2017; Moser and Lasher-Trapp 2017; Langhans et al. 2015; Yeo and Romps 2013; Romps 2010; Craig and Dörnbrack 2008; Stiller and Craig 2001; Grabowski 1995, 1993; Grabowski et al. 1993a,b, 1991; Ogura 1962; Lilly 1962), none of these works have studied the dry case in three dimensions, and more importantly questions 1–3 above remain open. Here we take advantage of the fact that DNS simulations have a well-defined Reynolds number, and vary it to address question 1 above. We also combine our simulations with a thermal tracking algorithm inspired by that of Romps and Charn (2015), which allows us to precisely define the thermal’s control volume as a function of time. Differentiating this volume with respect to height then yields the thermal’s (net) entrainment rate, allowing us to address question 2. Question 3 will be addressed using our thermal tracking, passive tracers, as well as a density anomaly budget for the environment, which is straightforward to calculate from simulation output.

We begin by briefly reviewing in the following subsection the analytic similarity theory of S57 for dry thermals in a neutral environment. Then in section 2 we describe our simulation setup and the broad properties of the simulated thermals. In section 3 we describe and apply a thermal tracking algorithm to find the volume of each thermal as a function of time. In section 4 we differentiate this function to find the entrainment, verifying the 1/*r* scaling in (1), and also verifying that detrainment in dry thermals is negligible. In section 5 we calculate the entrainment efficiency *e* in our simulations, exploring its Re dependence and comparing to previous work. We conclude in section 6.

### Analytic theory of thermals

In this subsection we review the classic similarity theory of S57. This theory does not invoke the relation (2), but instead assumes a geometric similarity between states of the thermal at different times, leading to the similarity constraints (4) and (7) below. These similarity constraints contain dimensionless constants which must be determined empirically. Though these constants can themselves be constrained by invoking (2) or the vertical momentum equation (e.g., Escudier and Maxworthy 1973; Turner 1964a), we nonetheless use the S57 theory to describe our thermals, for the sake of connecting with previous literature. We discuss possible elaborations to the S57 theory in the conclusions.

The S57 theory describes a “self-similar” thermal whose shape at any given time must be geometrically similar to its shape at any other time.^{1} From this self-similarity condition S57 deduces that the thermal’s height *z* must be linearly related to its radius *r* as

for constant *n* and some suitably chosen “virtual origin” *z*_{0}. This equation just describes the cone traced out by the flanks of the thermal as it rises.

S57 further argues, again by similarity, that the thermals’ volume *V* = *mr*^{3}, for constant *m*. S57 also assumes that detrainment of parcel fluid into the environment is negligible, in which case the gross fractional entrainment rate *ε* is also equal to the *net* fractional entrainment rate, which can be written as the fractional change in parcel volume with respect to height:

By measuring *r*(*z*) from photographs, S57 infers from (4) an average value of *n* ≡ (*dr*/*dz*)^{−1} ≈ 4, and hence *e* = 3/*n* ≈ 0.75.

S57 also argues that the Froude number *C*^{2} is constant, and hence that

where *B* is the thermal’s average Archimedean buoyancy (m s^{−2}). Additionally, in the absence of detrainment the thermal’s “mass anomaly” (i.e., integrated density anomaly, in kg) is conserved. This is proportional to *BV* and hence

## 2. Simulation setup and qualitative description

### a. Simulation setup

We perform direct numerical simulations of dry thermals in the Boussinesq approximation, in which the background density is assumed to be nearly constant. A buoyant fluid parcel which does not entrain or detrain thus maintains constant volume as it rises. Our simulations thus neglect effects from atmospheric density stratification which may be included in, for example, the anelastic approximation.^{2} The Boussinesq equations are

All dimensional variables are indicated with a tilde $\u2061(.~)$. Below we will nondimensionalize the equations. The fluid velocity $u\u02dc$ is assumed to be nearly incompressible, $p\u02dc$ is the pressure, $g\u02dc$ is the gravitational acceleration, and **e**_{z} is the unit vector in the vertical direction (parallel to gravity). We expand the density $\rho \u02dc$ as $\rho \u02dc=\rho \u02dc0+\rho \u02dc\u2032$ with $\rho \u02dc0$ constant, and assume $\rho \u02dc\u2032\u226a\rho \u02dc0$ such that density fluctuations are only important in the buoyancy term (the Boussinesq approximation). The density changes due to temperature fluctuations $T\u02dc\u2032$ according to $\rho \u02dc\u2032=\u2212\alpha \u02dcT\u02dc\u2032$, where $\alpha \u02dc$ is the (constant) coefficient of thermal expansion. The viscous and thermal diffusivity are $\nu \u02dc$ and $\kappa \u02dc$, respectively.

The thermal is initialized as a spherical temperature perturbation with diameter $L\u02dcth$ and temperature $\Delta T\u02dc$ (but with zero velocity). With these, we can define *dimensionless* variables, which do not have a tilde,

where [consistent with the energetics argument of (7)],

With this nondimensionalization, the equations become

Now the problem is entirely characterized by only two dimensionless numbers, the Reynolds number and the Prandtl number,

as well as the choice of noise we add to the temperature field. In this paper, we fix Pr = 1, close to the atmospheric value of Pr = 0.7.

One of our main goals is to determine how the thermal evolution varies with Re (question 1 above). We run a series of simulations with $Re=(2/10)\xd7104\u22486300$, as well as a series of simulations with $Re=(2/10)\xd7103\u2248630$. We refer to the higher-Reynolds-number simulations as *turbulent*, as they exhibit a range of spatial scales, chaos, etc. (see Fig. 2). These simulations are intended to be representative of the turbulent entrainment found in atmospheric convection (although their Reynolds number is still much smaller than in the atmosphere). In contrast, the lower-Reynolds-number simulations, which we term *laminar*, are not chaotic and only have a single spatial scale; these simulations may be more similar to the thermals described in simulations of ensembles of convective clouds such as Hernandez-Deckers and Sherwood (2016), Romps and Charn (2015), and Sherwood et al. (2013), in which each thermal is only resolved by a handful of grid cells.^{3} We also ran several simulations with Re ≈ 1800, but found them difficult to interpret as they are on the verge of turbulence, and thus exhibit large fluctuations. Thus, we will not discuss them further.

Reproducibility is important for science (Oishi et al. 2018). To aid readers to reproduce our results or to perform their own analysis of our data, we have created a github repository containing all the simulation scripts used to generate, analyze, and plot the data used in this paper. Although we could not include all the raw data from the simulations (which constitutes around 10 TB of data), we do include the reduced data necessary to reproduce the plots in this paper. The repository is available at https://github.com/lecoanet/thermals_entrainment.

We solve (10) using the open-source Dedalus^{4} pseudospectral code (Burns et al. 2016, 2019). We discretize the problem by expanding each quantity in a certain number of sine or cosine modes in *x*, *y*, and *z*. The domain extends from $\u22125L\u02dcth$ to $5L\u02dcth$ in the *x* and *y* directions, and from 0 to $20L\u02dcth$ in the *z* direction. For each direction, the normal velocity is expanded as a sine series, and the perpendicular velocities are expanded as cosine series. The pressure is expanded in cosine series in all directions and the temperature perturbation is expanded in cosine series in the horizontal directions and a sine series in the vertical direction. This corresponds to stress-free boundaries, and no-flux boundaries in the horizontal direction and isothermal boundaries in the vertical direction. For laminar simulations, we use 256 modes in the horizontal directions and 512 modes in the vertical direction. For turbulent simulations, we use 512 modes in the horizontal directions and 1024 modes in the vertical direction. We use the 3/2 padding rule when evaluating nonlinear terms to prevent aliasing errors. To time step the problem, we use a third-order, four-stage implicit–explicit Runge–Kutta timestepper [Ascher et al. (1997), where linear terms are treated implicitly, and nonlinear terms are treated explicitly], with the time step size set by a Courant–Friedrichs–Lewy condition with prefactor 0.7. We run the simulations for ≈60 time units, which is enough time for the thermals to approach the top boundary.

To construct the initial condition, we first specify a spherical density perturbation,

where

and erf is the error function (e.g., Abramowitz and Stegun 1968), *x*_{0} = *y*_{0} = 0, *z*_{0} = 1.5, and *r*_{0} = 0.5 is the initial radius. We take Δ*r* = 0.1 as a smoothing length. To break the symmetries of the problem, we use the initial density perturbation $\rho sph\u2032\xd7\u2061[1+N\u2061(x,y,z)]$, where *N* is a noise field specified in terms of the amplitude and phase of its sine and cosine modes. For each mode **k** = (*k*_{x}, *k*_{y}, *k*_{z}) with *k*_{x}, *k*_{y}, *k*_{z} < 128 × 2*π*/10, we set *N*_{k} by *A*(1 + *k*^{2})^{−1/6}*ξ*sin(*ϕ*), where *A* is an amplitude, *ξ* is a normally distributed random variable and *ϕ* is a random variable uniformly distributed from 0 to 2*π*. We pick *A* such that the root-mean-square density perturbations are 0.21.^{5} At both Reynolds numbers, we run five simulations with different choices of the initial random seed, and thus different noise fields *N*. This allows us to compute an ensemble average over our simulations.

As the thermal evolves, we expect its radius *r* and typical velocity *u*_{rms} to evolve such that *ru*_{rms} stays approximately constant (since $r~t$ and $u~w~1/t$; see the “Analytic theory of thermals” subsection above). This means that the Reynolds number of the thermal is approximately constant over the simulation. At early times, the thermal takes up a small part of the domain, and is thus harder to resolve. The resolution of our turbulent simulations is insufficient to fully resolve the flow at the early stages of the simulations, leading to low-amplitude Gibbs’ ringing, especially in the density field. In these early stages, the viscous dissipation scale is about twice the grid spacing. We thus limit our analysis to later stages of the simulation (e.g., *t* > 10), where the thermal is larger and well resolved (i.e., the viscous dissipation scale is larger than the grid spacing).

### b. Qualitative analysis

In Fig. 1 we show 2D vertical slices of the time evolution of a laminar, Re = 630 thermal and a turbulent, Re = 6300 thermal (left and right panels of each pair, respectively), drawn from our ensemble of simulations and initialized with the same initial noise field. ( Appendix A discusses the effect of changing the initial noise field on the thermal.) Movies showing the evolution of these thermals are included in the online supplemental material (https://doi.org/10.1175/JAS-D-18-0320.s1). Both thermals take the form of buoyant vortex rings (see Fig. 2) and entrain ambient fluid as they rise, at seemingly comparable rates; this is consistent with Turner’s claim that (3) drives entrainment, and provides a preliminary answer to question 1 above. The thermals grow significantly larger with time, and because of the associated dilution exhibit a decrease in vertical velocity as well as density perturbation, as evident in the changing color scales. This decrease in vertical velocity despite increasing impulse is due to “mixing drag,” that is, the fact that newly entrained, neutrally buoyant mass must accelerate by absorbing momentum from the existing thermal. (This growth and deceleration of the thermal with time contrasts the behavior of moist thermals in more realistic simulations, which we discuss further in the conclusions.) Both thermals also detrain, as there is a trail of fluid left below each thermal, but it is unclear from this figure how significant this is.

In Fig. 2 we plot the azimuthal vorticity (*ω*_{ϕ} = ∂_{z}*u*_{r} − ∂_{r}*u*_{z}, where **u** is represented in cylindrical coordinates) of these two thermals at *t* = 30. Although the turbulent simulation has substantial vorticity at small scales, both simulations show strong vorticity maxima in the core of the vortex ring. Unlike, say, Hill’s spherical vortex which has only positive vorticity throughout, these thermals also have some slight negative vorticity near their centers; thus the vertical velocity is not maximal at the center of the thermal.

While our laminar and turbulent thermals have much in common, differences do exist between the two cases. The laminar thermal forms a well-defined (though somewhat asymmetric) vortex ring, whereas in the turbulent case the vortex ring structure is narrower in vorticity, but more diffuse in density. Furthermore, the turbulent thermal clearly entrains slightly more than the laminar one, leading to lower heights, vertical velocities, and smaller density anomalies at a given time. The rest of this paper will focus on measuring this entrainment, as well as detrainment, to yield more quantitative answers to questions 1–3 above.

## 3. Thermal volume

In Fig. 1 we also plot the boundary of the thermal, whose definition and calculation we now describe. This boundary will be important for distinguishing the thermal from ambient fluid, as well as from its trail of detrained fluid, and will thus be critical for measuring both entrainment and detrainment.

Heuristically, we think of a thermal as a region of fluid that moves coherently. Following Romps and Charn (2015), we thus define the thermal to be the axisymmetric volume whose averaged vertical velocity matches the velocity of the thermal’s top.^{6} This definition agrees with the well-known notion of the vortex ring “bubble” or “atmosphere” (Akhmetov 2009; Shariff and Leonard 1992, their section 3.2). This definition is based on fundamental fluid dynamics and treats thermals as coherent dynamic entities, and as such is arguably preferable (at least in this context) to conventional definitions of cloud control volumes based on arbitrary thresholds of cloud condensate, vertical velocity, and/or buoyancy (e.g., de Rooy et al. 2013; Dawe and Austin 2011; Romps 2010). Similar tracking algorithms to ours have been used in the recent studies of Hernandez-Deckers and Sherwood (2018), Morrison and Peters (2018), Hernandez-Deckers and Sherwood (2016), and Sherwood et al. (2013). Our algorithm seems to work well for our simulations (see movies in the online supplemental material), but is of course not perfect, and sometimes misclassifies fluid as being outside the thermal, even though it continues to rise coherently with the rest of the thermal volume. This complicates our detrainment calculation in section 4.

We now briefly describe the algorithm for identifying this boundary; full details are given in appendix B. The first step is to calculate the position of the top of the thermal *z*_{t}, which we do using a density threshold. To calculate the thermal-top velocity *w*_{top}, we first perform a fit to *z*_{t}, and then take the derivative analytically.

In Fig. 3, we plot the thermal-top height as a function of time for each of our simulations. The curves have negative concavity, again indicating deceleration. The turbulent thermals have lower heights than the laminar thermals at the same time (consistent with Figs. 1 and A1). For each curve, we also plot the fit we use to calculate *w*_{top} [see (B1)]. The fit is very good except at early times *t* < 10. This is expected, as it takes some time for the thermal to spin up and reach the self-similar regime. Furthermore, as described earlier, our turbulent simulations are not very well resolved for *t* < 10, although this is not expected to lead to large differences in bulk properties like the thermal-top height. The good fit to the $t$ time dependence gives a first quantitative check of the classical theory of S57 (see the “Analytic theory of thermals” subsection). Because the thermal boundary is calculated based on *z*_{t} and our fit is not good for *t* < 10, we expect the thermal properties will not be accurately calculated during these early stages of the simulations.

To find the axisymmetric volume translating with this velocity, we calculate the Stokes streamfunction *ψ* associated with the azimuthally averaged velocities, **u**_{axi} = [*u*_{axi}(*r*, *z*), *w*_{axi}(*r*, *z*)]. The Stokes streamfunction satisfies (Batchelor 2000)

where **e**_{ϕ} is the unit vector in the azimuthal direction. This is the streamfunction associated with the thermal velocity, in the frame moving upward at the thermal-top velocity. The thermal boundary is the line of constant *ψ* starting from the point (*r*, *z*) = (0, *z*_{top}) where the axisymmetric vertical velocity matches the thermal-top velocity, that is, *w*_{axi}(0, *z*_{top}) = *w*_{top}. Thus, by definition, the thermal has no net mass flux in the frame moving with the thermal-top velocity, and its average vertical velocity in the rest frame is equal to the thermal-top velocity. More details about this calculation can be found in appendix B.

Figure 1 shows our calculation of the thermal boundary in our simulations. This automated procedure seems to do a good job of identifying the thermal, and seems to match the region one might pick out “by eye.” When visualizing the simulation data, it is easiest to identify the thermal when looking at the vertical velocity (rather than the density), because it is smoother. The vertical velocities are strong and coherent within the thermal volume, as opposed to the density field which is highly variable and exhibits pockets of completely unmixed ambient fluid within the thermal (see, e.g., the *t* = 30 panel of the Re = 6300 simulation in Fig. 1). This makes it more difficult to determine the thermal volume by only looking at the density field. Note that these pockets of unmixed environmental fluid seem to result from engulfment at the rear of the thermal by the large-scale vortex ring circulation, and are only later mixed/diffused in, as also noted by previous authors (Moser and Lasher-Trapp 2017; Johari 1992; Sànchez et al. 1989; Turner 1964b; Woodward 1959; S57). This is consistent with entrainment originating from (3), rather than turbulence.

The boundary of the thermal tells us the radius of the thermal as a function of height *r*(*z*). We can then calculate the thermal volume by simply summing *πr*(*z*)^{2}*dz* over *z*. This thermal volume is shown in Fig. 4. The classical S57 theory predicts that *z* ~ *r* and hence that volume should increase like $V~zt3$, which we also find in our simulations. The turbulent thermals have larger volumes than the laminar thermals at the same height, again suggesting somewhat larger entrainment rates for the turbulent thermals. There is also a nonintuitive decrease in the volume in our lower-Reynolds-number simulations when *z*_{t} becomes large, but this is an artifact of the thermal interacting with the top boundary (see Fig. 3).

## 4. Net entrainment and detrainment

While the increasing volume of our thermals with time must be due to entrainment of environmental fluid, it can also in principle be offset by detrainment of fluid from the thermal to the environment. This possibility is not so far-fetched, as we clearly see in the simulations a trail of detrained fluid below the thermal (Fig. 1). Accordingly, and analogously to *ε*, we define the fractional *detrainment rate δ* to be the fraction of the thermal’s volume that it detrains per unit vertical distance traveled, again in units of m^{−1}. We may then generalize (5) by allowing for detrainment and properly equating *d*ln*V*/*dz* to the net fractional entrainment rate *ε*_{net}, giving

Since *ε*_{net} can be calculated from the volume only, we focus on it first; estimating *δ* will require additional machinery, to be discussed later.

One goal of our paper is to check the *ε*_{net} ~ 1/*r* scaling suggested by (6). To do so, we must calculate *r* and *d*ln*V*/*dz* in our simulations. We define the thermal radius *r*_{th} to be the maximum radius of the thermal boundary. To calculate the entrainment rate, we calculate *dV*/*dt* with a second-order central differencing formula (as implemented in numpy’s gradient function). Then we calculate *ε*_{net} by dividing *dV*/*dt* by *Vw*_{top}. Taking the derivative via finite differences introduces noise in the entrainment rate. Although we could have reduced this noise by fitting the data and taking a derivative analytically, or by smoothing, we decided to minimize the data processing to eliminate any biases so introduced. On the other hand, it seems acceptable to fit the thermal-top heights via (B1) because it is in an early stage of the data analysis.

Figure 5 shows the net entrainment rate *ε*_{net} as a function of thermal radius *r*_{th} on log–log scales. A slope of −1 is evident, providing a direct verification of the 1/*r* scaling of *ε*_{net}.^{7} We find this dependence for both individual simulations (thin lines), and the mean (thick lines; calculated by binning *r*_{th} in steps of 0.1, and then calculating the average *ε*_{net} and average *r*_{th} over simulations with a given Reynolds number in each bin). Note that the turbulent thermals indeed have slightly larger net entrainment rates than the laminar thermals. There is significantly more scatter in entrainment rates among the turbulent thermals than the laminar thermals, likely due to turbulent chaos.

Next we turn to detrainment. A simple approach to measuring detrainment would be to diagnose changes in the thermal’s mass anomaly *BV* with time. This approach, however, assumes that *B* is uniform throughout the thermal, so that a change in *BV* can be confidently associated with a change in *V*. Such an assumption seems reasonable for our turbulent thermals, but not for the laminar ones (Fig. 1). Also, if some mass anomaly leaves our idealized, axisymmetric thermal volume at some time but then reenters it soon afterward and continues to move upward, this would still be measured as detrainment, which seems undesirable.

To deal with these complications, we take the following approach. In the turbulent case we integrate *ρ*′ in the *environment below the thermal* as a measure of detrained mass anomaly *M*_{d}, and track the rate at which *M*_{d} changes. We only track detrained mass anomaly below the thermal because the truly detrained fluid rises less quickly than the thermal, and thus falls below it. This method for measuring detrainment neglects the significant mass anomaly exiting the thermal volume’s sides and top which subsequently reenters the thermal from the sides (in contrast to entrained fluid, which enters from below). We think of this fluid as always contained within some “platonic” thermal, and never as having really been detrained. The fact that it leaves our axisymmetric thermal volume is an imperfection of our idealized definition of the thermal boundary.

As for the laminar case, because *B* is not uniform in this case, we felt compelled to use a different method entirely. For this case we use a passive tracer initialized within the thermals at several times within each simulation, and track the rate at which the tracer is detrained into the environment. Full details of these calculations are given in appendix C.

With these methodologies in place, we calculate detrainment rates *δ* and directly compare to net entrainment rates *ε*_{net}. Ensemble averages of these quantities are shown in Fig. 6 (the averages are taken in the same way as for Fig. 5, but with *z*_{ct} bins of size 0.5). Both the laminar and turbulent detrainment rates are an order of magnitude smaller than the entrainment rates, justifying the approximation of S57 and others in neglecting detrainment. Furthermore, the detrainment rates, despite being calculated using different algorithms, are similar for the two sets of calculations, and thus also only marginally sensitive to Re. However, it does seem that the detrainment process itself depends on Re—it is steady and primarily through the bottom for the laminar thermals, but intermittent and primarily through the sides for the turbulent thermals. In contrast, in both cases the entrainment is predominantly from the bottom of the thermal. Note also that the detrainment rate for the laminar thermal decreases like *z*^{−1}, and since *z*_{t} ~ *r*_{th}, this suggests that *δ* ~ 1/*r*, just like the net entrainment rate. The detrainment rate for the turbulent thermals is too noisy to determine a temporal variation, but we hypothesize that it may also decrease like *r*^{−1}. Figures 5 and 6 thus show that *ε*_{net} ~ 1/*r* in all cases and *δ* ~ 1/*r* in the laminar case, which in conjunction with *ε* = *ε*_{net} + *δ* provides strong support for the entrainment assumption (1).

## 5. Entrainment efficiency

Since Fig. 6 shows that *δ* ≪ *ε*_{net} and hence that *ε* ≈ *ε*_{net}, we may estimate the entrainment efficiency *e* in (1) by multiplying *ε*_{net} by *r*_{th}. This should be approximately constant over the course of a simulation. We plot *e* as a function of time in Fig. 7. We find that the averaged entrainment efficiency (dark thick lines) is indeed roughly constant with time, except at the very beginning of the simulations or near the end when the thermals begin to interact with the top boundary. Furthermore, the turbulent thermals on average have a larger entrainment efficiency than the laminar thermals. However, there is substantial spread between the different simulations (especially for the turbulent thermals), and an individual turbulent thermal may have a smaller entrainment efficiency at a given time than an individual laminar thermal.

The average entrainment efficiency drops quickly for heights *z*_{t} > 16 for our laminar simulations and for heights *z*_{t} > 18 for our turbulent simulations. As described above, this is because the thermals are approaching the top boundary of the simulation domain. However, there is a slower, more secular decrease in the entrainment efficiency with height, for both Reynolds numbers. This is not explained by the S57 theory (see the “Analytic theory of thermals” subsection). Nevertheless, it does seem that the lowest-order description of the problem is that the entrainment efficiency is roughly constant.

To more quantitatively compare the entrainment efficiencies of the simulations with different Reynolds numbers, we also take a time average of the entrainment efficiency between *z*_{t} = 6 and *z*_{t} = 16. This ensures we are not influenced by initial transient effects at early stages of the simulations, nor by interactions with the boundaries at late stages. The values are reported in Table 1. We find that the laminar thermals have a lower entrainment efficiency than the turbulent thermals by only about 20%, yielding a quantitative answer to question 1. The entrainment efficiencies of our turbulent simulations are similar to those found in laboratory experiments with an initial aspect ratio ≥1, but smaller than those of R61, S57, and MTT56 which had an initial aspect ratio of 0.5. This is quantitatively consistent with the dependence on aspect ratio found in Lai et al. (2015). Also note that entrainment efficiencies of order slightly less than 1 are broadly consistent with typical fractional entrainment rates of *O*[(1 km)^{−1}] measured in simulations and observations (e.g., de Rooy et al. 2013; Malkus 1954), since a typical radius for clouds is *r* ≈ *O*(1 km).

## 6. Discussion and conclusions

In this paper we present a suite of direct numerical simulations of resolved, dry thermals. We initialize the thermals with a sphere of buoyant fluid released at the bottom of our simulation domain, with random perturbations to break symmetries. We run an ensemble of simulations with different random perturbations (but with the same statistical properties), and with Reynolds numbers of either Re ≈ 630 (laminar) or Re ≈ 6300 (turbulent). To quantitatively measure the entrainment and detrainment associated with each thermal, we implement a thermal tracking algorithm inspired by Romps and Charn (2015), where the resulting thermal boundary is shown with a black line in Fig. 1.

Our entrainment calculation allows us to directly verify the long-standing entrainment assumption *ε* = *e*/*r* for both laminar and turbulent thermals. Furthermore, we find that both cases entrain comparably, with entrainment efficiencies *e* that differ by only 20%. This is consistent with the idea of entrainment governed by (3) for buoyant vortex rings (Turner 1957). We also found that detrainment rate is roughly an order of magnitude smaller than the net entrainment rate, suggesting that net entrainment rates are approximately equal to *gross* entrainment rates.

The finding that entrainment has a strong laminar component should be contrasted with the plume picture of convection, in which entrainment is thought of as purely turbulent (e.g., Turner 1986; Squires and Turner 1962; Kuo 1962; MTT56). Also note that entrainment governed by (3) is distinct from the “dynamic” entrainment proposed by Houghton and Cramer (1951) and incorporated in later studies (e.g., Morrison 2017; de Rooy and Siebesma 2010; Ferrier and Houze 1989; Asai and Kasahara 1967), since dynamic entrainment is a consequence of *positive* vertical acceleration, whereas we see negative vertical acceleration here.

While the results presented here help paint a picture of entrainment in dry convection, a key question is how these results extend to atmospheric moist convection. As mentioned in section 2b, simulated moist atmospheric thermals (without resolved turbulent mixing) exhibit a fairly constant size and vertical velocity (Hernandez-Deckers and Sherwood 2016; Romps and Charn 2015; Sherwood et al. 2013), in contrast to our dry thermals. These moist thermals also appear to have a qualitatively different momentum budget than dry thermals, with little to no mixing drag but a significant perturbation pressure drag (Morrison and Peters 2018; Romps and Charn 2015; Sherwood et al. 2013; de Roode et al. 2012); this contrasts with the dry thermal’s momentum budget, which can be completely characterized by mixing drag and virtual mass effects (Turner 1964a).

What might be the sources of these differences? Regarding entrainment/detrainment, the fact that moist thermals exhibit a constant size implies that *ε*_{net} ≈ 0 and hence *δ* ≈ *ε*. Since our gross entrainment rates seem to be similar to those measured in more realistic simulations, the discrepancy is likely due to detrainment. There are two likely candidates for a sharp increase in detrainment in more realistic cases. One is the stable stratification of the real atmosphere, which reduces the buoyancy of unsaturated parcels as they rise. The other candidate is the mixing of cloudy and clear air, which causes evaporation of condensate and again a decrease in buoyancy. These processes ultimately generate negatively buoyant air which is likely to detrain, a process formalized in the “buoyancy-sorting” scheme of Raymond and Blyth (1986) and still invoked in contemporary convection schemes (e.g., Romps 2016).

As for differences in dynamics, the argument of Turner (1957) relies on a constant circulation *K*, which itself relies on the existence of a circuit through the center of the vortex ring along which the density anomaly is negligible. Moist thermals, however, can exhibit self-generated buoyancy in their center due to latent heat release (e.g., Morrison and Peters 2018; Romps and Charn 2015), so such a circuit may not exist. To the contrary, such buoyancy in the thermal’s center may very well *increase* the circulation *K*, which by (2) would help *r* to increase more slowly or even stay constant in the face of increasing *P*, as is indeed observed. Future work could assess how the constraint (2) operates in the moist case, and also assess how a stable stratification and condensate evaporation contribute to buoyancy sorting and enhanced detrainment.

In closing, we note that while this paper underwent review, several additional studies have appeared (employing these or similar simulations) which round out the analysis given here. The “virtual mass” or “effective buoyancy” of these thermals (e.g., Batchelor 2000; Jeevanjee and Romps 2016) is studied in Tarshish et al. (2018). A theory for thermals based not on the similarity assumptions (4) and (7), but rather on the relation (2) and the vertical momentum equation, is given in McKim et al. (2019). This theory requires no additional, empirical constants, and further suggests that the 20% increase in *e* in the turbulent case stems largely from changes in the large-scale characteristics of the flow, such as the circulation *K* in (3). This theory is extended from the Boussinesq case to the anelastic, density-stratified case in Anders et al. (2019).

## Acknowledgments

DL is supported by a PCTS fellowship and a Lyman Spitzer Jr. fellowship. NJ was supported by the Visiting Scientist Program in Princeton’s Atmosphere and Ocean Science program, as well as a Harry Hess fellowship from the Princeton Geoscience Department. Computations were conducted with support by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center on Pleiades with allocations GID s1647 and s1439.

### APPENDIX A

#### Ensemble Characteristics

We find that different initial noise fields can lead to different thermal evolution, even though the statistical properties of the noise (e.g., power spectrum, root-mean-square fluctuations) are identical. Figure A1 shows 2D vertical slices of the density in all 10 of our simulations at *t* = 30. Each vertical pair represents simulations with the same initial condition; the top plots show laminar thermals and the bottom plots show turbulent thermals. The rightmost plots are from the same simulations as shown in Fig. 1.

Although there is substantial variation between the different simulations with the same Reynolds number, we can still clearly identify some trends with the Reynolds number. In particular, as we described for the single choice of initial noise field in Fig. 1, the laminar thermals systematically rise higher and have larger density anomalies than the turbulent thermals. Thus, this observation is robust to the exact choice of initial noise.

### APPENDIX B

#### Thermal Tracking Algorithm

Here we will describe the details of our thermal tracking algorithm. The first step is to determine the thermal-top height *z*_{t}(*t*). At each time (we have full volume outputs every $1/10\u22480.36$ time units), we first calculate the horizontal average of *ρ*′ at each height, ⟨*ρ*′⟩_{x,y}. We define a cutoff value of 0.1max_{z}⟨*ρ*′⟩_{x,y}. The height is very insensitive to the cutoff value because the density field is very sharp. The thermal-top height is the highest point at which ⟨*ρ*′⟩_{x,y} is greater than this cutoff value in absolute value.

We find that the volume of the thermal shows sensitivity to the value of *w*_{top}, so care must be taken when calculating *dz*_{t}/*dt*. We find the best way to calculate this derivative is to first fit

based on (9), using a nonlinear least squares fit, and then take the derivative analytically to determine *w*_{top}.

Next we calculate the horizontal midpoint of the thermal. We define the midpoint using the vertical velocity. At each height, we define *x*_{m} by

and similarly for *y*_{m}. We then pick the vertical height which maximizes $\u2211w>0w$. The horizontal midpoint of the thermal is then (*x*_{m}, *y*_{m}) at this height. Using a mass weighting yields a similar midpoint.

We then calculate the azimuthal average of the vertical velocity at every vertical height. To do this, we define a radial grid which has the same grid spacing as the *x* or *y* grids, but only extends from *r* = 0 to *r* = 5. To calculate *w*_{axi}(*r*, *z*), we take the average of *w* in azimuthal rings of width Δ*r* centered around *r* + Δ*r*/2.

Once we have *w*_{axi}(*r*, *z*), we can calculate the Stokes streamfunction *ψ* satisfying (11). We do this using Dedalus (Burns et al. 2016). First, we (spectrally) interpolate *w*_{axi}(*r*, *z*) onto a Chebyshev grid in the radial direction. We then solve the linear boundary value problem,

with the boundary condition *ψ*(*r* = 0) = 0.

The boundary of the thermal is then the contour *ψ* = 0. To find this contour, at every height, we find the maximum of *ψ*. If the maximum occurs at radius larger than 0.18 and is positive, then we use root finding to find a zero of *ψ* at larger radii than the maximum. This defines the radius of the thermal at each height. The lower cutoff of 0.18 is used to ensure the thermal is one simply connected volume, and to prevent the thermal tracking from thinking that pockets of detrained fluid near the midpoint are part of the thermal. In one simulation, there were several times at which such pockets of fluid were identified as being part of the thermal (even though they were far below and disconnected from the real thermal). In this case, we manually removed these spurious thermal components (i.e., set the radius of the thermal at these heights to zero).

### APPENDIX C

#### Detrainment

Here we detail our detrainment calculations. These are somewhat involved, as the detrainment rate is low, which makes it difficult to calculate. Furthermore, the physical process differs between laminar and turbulent thermals, leading us to different methodologies for the two cases.

We conceptually think of a thermal as a volume of fluid rising coherently at a certain velocity. Then for fluid to be detrained, it must be rising slower than the thermal, and thus must ultimately end up below the thermal. In the turbulent case, however, we find that often fluid initially leaves the thermal volume on the tops and sides, and that some of this fluid is subsequently reincorporated into the thermal. Although one might think of this process as a detrainment followed by subsequent entrainment event, we instead think of this fluid as being part of the thermal the whole time, even if it was not correctly identified as such by our thermal tracking algorithm. This is a limitation of our assumption of axisymmetry when defining the thermal volume.

Thus, to diagnose detrainment in turbulent thermals, we note that the tail of detrained fluid (e.g., as seen in Fig. 1) has similar density to the average density of the thermal (this does not seem to be the case for laminar thermals). If we assume the average density of the detrained fluid matches the average density of the thermal, then the volume of detrained fluid can be measured by the mass anomaly of fluid detrained (integral of *ρ*′ over that volume). Since the background *ρ* is not even defined in our simulations due to the nondimensionalization, we will refer in this appendix to such mass anomalies simply as “mass.”

Because the environmental fluid has *ρ*′ = 0 (from the Boussinesq approximation), entrainment of environmental fluid cannot change the mass of the thermal $Mth\u2261\rho th\u2032V$, where $\rho th\u2032$ is the average density perturbation of the thermal. Entrainment changes $\rho th\u2032$ and *V* individually but not their product. Thus, if we know the total “mass” of detrained fluid *M*_{d} throughout the simulation, we may approximate the detrainment rate as

Defining *M*_{d} requires care, however. If we simply integrate *ρ*′ over the whole domain excluding the thermal, we integrate not only over detrained fluid below the thermal but also fluid above or on the sides as well (e.g., above the thermal at *t* = 50), much of which will subsequently be re-entrained. With this in mind, we define the detrained mass *M*_{d} to be only the mass below *z*_{bot}, the lowest point in the thermal volume. Once a detrained fluid element is below the thermal, it will very rarely be re-entrained because it is rising less quickly than the thermal. (Recall that in our Boussinesq model, the background has *ρ*′ = 0, so does not contribute to *M*_{d}.) We also calculated the mass below *z*_{bot} − *r*_{th}/4. Although there is less mass below that point than below *z*_{bot} at any given time, it leads to similar detrainment rates. This suggests that almost all of the mass that falls below *z*_{bot} will soon fall below *z*_{bot} − *r*_{th}/4 (i.e., it has truly been detrained).

Figure C1 shows in thin lines *M*_{d} for all of our simulations, normalized by the total initial mass, $M0=\u2212(4\pi /3)r03$. In the turbulent simulations, the signal is noisy; there can be abrupt increases in the detrained mass, as detrainment events are often discrete, and the detrained mass occasionally decreases as the thermal rises, due to some of the mass that falls below *z*_{bot} becoming re-entrained. To alleviate some of these issues, we also plot the ensemble-average detrained mass (dark thick lines in Fig. C1, computed as in Fig. 5 but with *z*_{ct} bins of size 0.5).

In its early evolution (*z*_{t} < 5, not shown), the thermal detrains a nonnegligible amount of fluid, roughly 5%–10% of *M*_{0}. After this initial transient, however, there is on average a relatively slow growth of detrained mass. Indeed, on average the detrained mass only grows by a factor of 4–5 from *z*_{ct} = 5 through the end of the simulation for turbulent thermals. In contrast, the volume changes by almost a factor of 100 (Fig. 4), suggesting detrainment is indeed much weaker than entrainment.

Figure C1 also shows the detrained mass for laminar thermals, which increases continuously and smoothly as the thermal rises. The laminar thermals detrain much less mass than the turbulent thermals. However, the detrained fluid in laminar thermals (see the fluid immediately below the thermal in Fig. 1) is much less dense than the average density of the thermal. This is because laminar thermals’ density is much more concentrated in a vortex core than turbulent thermals’ density. Thus, we cannot approximate the detrained *volume* to be proportional to the detrained “mass” for laminar thermals.

To measure the detrainment rate for laminar thermals, then, we instead run a set of dye simulations which track the movement of a passive scalar (“dye”) both into and out of the thermal. We solve the equation

for the passive scalar field *s*. We use Schmidt number of unity (*D* = *ν*). The field *s* is represented as a cosine series in the *x* and *y* directions, and a sine series in the *z* direction. The field *s* denotes the density of the passive tracer; we will use *c* to denote a volume integral of *s*, that is, the amount of the passive tracer.

Every $10$ time units, we reinitialize *s* to be nearly unity within the thermal, and nearly zero outside, as described below. To define the thermal, we use the thermal radius as a function of *z* calculated with our volume tracker (section 3). We use this to define the radius of the thermal as a function of angle from the vertical from the center of the thermal—the vertical center of the thermal is defined as the average of the highest and lowest points in the thermal. Then for each grid point in the domain, we calculate the angle to the thermal center relative to the vertical, and then calculate the radius of the thermal for that angle *r*_{th}(*θ*) using linear interpolation. Then the dye field is initialized to

where Δ*r* = 0.1 is a smoothing length.

After each reinitialization, the dye field is advected along with the thermal. Almost all of the dye stays within the thermal, but is diluted due to the entrainment of environmental fluid. If we calculate the entrainment rate by measuring the amount of environmental fluid that enters the volume, we find entrainment consistent with our measurements purely in terms of the thermal volume. This is because entrainment is much stronger than detrainment, so the entrainment rate is almost entirely determined by the change in volume. There is a small amount of dye which is initialized within the thermal that leaves the thermal. Some dye leaves the thermal volume due to diffusion, but some is advected out, that is, is detrained. As above, we consider dye detrained if it is below the lowest point of the thermal. We denote amount of dye below the lowest point of the thermal by *c*_{d}.

To quantify detrainment, we calculate how *c*_{d} changes relative to the total amount of dye, *c*_{tot}. The total dye is approximately equal to the volume of the thermal at the previous reinitialization time. We plot the fraction of detrained dye, *c*_{d}/*c*_{tot}, in Fig. C2. The fraction of detrained dye at reinitialization is marked with a cross and the last fraction recorded before the following reinitialization is marked with a circle. There is a small amount of detrained dye immediately after reinitialization because the dye field is not set to exactly zero outside the thermal [see (A6)]. Sometimes the dye near the boundaries is advected into the thermal, causing the amount of detrained dye to decrease in time. Nevertheless, the amount of detrained dye mostly increases until the dye field is reinitialized.

We then define the detrainment rate for laminar thermals to be

where Δ*c*_{d} and Δ*z*_{t} are the change in the detrained dye and height over the course of one reinitialization step (of time $10$). Thus, this derivative represents the slope of the line between the cross and circle of each line segment in Fig. C2. We use this “course-grained” derivative because there are transients associated with the reinitialization which makes it difficult to interpret derivatives on finer scales. Since *δ* is a function of *z*_{t}, we interpret our calculated value as corresponding to the height *z*_{t,0} + Δ*z*_{t}/2.

Figure C2 also shows the fraction of detrained dye for a turbulent simulation (corresponding to the rightmost simulation in Fig. A1). Strikingly, this calculation shows almost no detrainment for *z*_{t} > 13. This appears to contradict Fig. 1, which appears to show detrained fluid below the turbulent thermal. Turbulent thermals have very little detrained dye because it takes longer than $10$ time units for dye which leaves the thermal (typically from the top or the sides) to fall below the lowest point of the thermal. The detrained fluid is then reinitialized to having *c* = 0 before it can fall below the thermal. This suggests the calculation needs more “memory” in order to correctly calculate detrainment. However, this would require longer times between reinitialization, which decreases the time resolution of the detrainment rate. For these reasons, we appeal to the assumption that the density of detrained fluid is equal to the average density of the thermal to calculate detrainment for turbulent thermals.

## REFERENCES

## Footnotes

Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JAS-D-18-0320.s1.

^{a}

Current affiliation: Program in Atmosphere and Ocean Sciences, Princeton University, Princeton, New Jersey.

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

This self-similar solution is supposed to hold until the thermal interacts with an environmental boundary, discontinuity, or additional length scale, e.g., the thermal size becomes comparable to a scale height, or the thermal approaches the tropopause in the atmosphere or a boundary in experiments or simulations.

^{2}

The results of this paper are extended to include the effects of density stratification in Anders et al. (2019).

^{3}

Although our thermals could entrain via diffusion, the diffusion time scale is $~rth2/\kappa $, whereas the entrainment time scale is ~(*wε*_{net})^{−1}; taking their ratio we find diffusive entrainment is less efficient by a factor of ≈Re and thus is negligible for our simulations.

^{4}

More information is available at dedalus-project.org.

^{5}

We also found similar results for larger *A*; smaller values of *A* led to artificially symmetric flows.

^{6}

It is relatively easy to identify the top of the thermal, as there is a sharp change in density at the top of the thermal.

^{7}

Although the thermal is still spinning up for *r* ≲ 1, the 1/*r* scaling appears to roughly hold even then.