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,

 
ε=e/r
(1)

(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:

 
P=πρKr2.
(2)

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

 
d(r2)dt=1πρKdPdt>0,
(3)

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:

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

  2. Does directly measured entrainment indeed obey the 1/r scaling in (1)?

  3. 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

 
zz0=nr
(4)

for constant n and some suitably chosen “virtual origin” z0. 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 = mr3, 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:

 
ε=dlnVdz(no detrainment).
(5)

Equation (5) along with (4) yields a 1/r entrainment law:

 
ε=dlnVdz=3dlnrdz=3n1r.
(6)

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 C2 is constant, and hence that

 
w2=C2Br,
(7)

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

 
B=B0r03/r3,
(8)

where subscript “0” denotes initial values. Substituting (8) and (4) into (7) and integrating with respect to time yields

 
zz0=att0
(9)

for some suitably chosen t0, where a is a function of n, C2, and m (see S57 for details). We will verify and utilize the thermal trajectory, (9), in the course of constructing our thermal tracking algorithm.

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

 
t˜u˜+ρ˜01˜p˜ν˜˜2u˜+g˜ρ˜ρ˜0ez=u˜˜u˜,t˜T˜κ˜˜2T˜=u˜˜T˜,˜u˜=0.

All dimensional variables are indicated with a tilde (.~). Below we will nondimensionalize the equations. The fluid velocity u˜ is assumed to be nearly incompressible, p˜ is the pressure, g˜ is the gravitational acceleration, and ez is the unit vector in the vertical direction (parallel to gravity). We expand the density ρ˜ as ρ˜=ρ˜0+ρ˜ with ρ˜0 constant, and assume ρ˜ρ˜0 such that density fluctuations are only important in the buoyancy term (the Boussinesq approximation). The density changes due to temperature fluctuations T˜ according to ρ˜=α˜T˜, where α˜ is the (constant) coefficient of thermal expansion. The viscous and thermal diffusivity are ν˜ and κ˜, respectively.

The thermal is initialized as a spherical temperature perturbation with diameter L˜th and temperature ΔT˜ (but with zero velocity). With these, we can define dimensionless variables, which do not have a tilde,

 
˜(L˜th1),T˜(ΔT˜)T,ρ˜(α˜ΔT˜)ρ,u˜(u˜th)u,p˜(ρ˜0u˜th2)p,t˜(u˜th/L˜th)t,

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

 
u˜th2g˜L˜thα˜ΔT˜ρ˜0.

With this nondimensionalization, the equations become

 
tu+pRe12u+ρez=uu,
(10a)
 
tρPr1Re12ρ=uρ,
(10b)
 
u=0.
(10c)

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

 
Re=u˜thL˜thν,
 
Pr=ν˜κ˜,

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)×1046300, as well as a series of simulations with Re=(2/10)×103630. 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 Dedalus4 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 5L˜th to 5L˜th in the x and y directions, and from 0 to 20L˜th 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,

 
ρsph=12[erf(r0rICΔr)1],

where

 
rIC(xx0)2+(yy0)2+(zz0)2,

and erf is the error function (e.g., Abramowitz and Stegun 1968), x0 = y0 = 0, z0 = 1.5, and r0 = 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 ρsph×[1+N(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 = (kx, ky, kz) with kx, ky, kz < 128 × 2π/10, we set Nk by A(1 + k2)−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 urms to evolve such that rurms 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.

Fig. 1.

2D vertical slices at y = 0 of the (top) vertical velocity and (bottom) density of two thermals with different Reynolds numbers. The two thermals are initialized with the same noise field. At each time, the left panel shows the laminar thermal (Re = 630), and the right panel shows the turbulent thermal (Re = 6300). We use the same color scale for the two thermals at the same time; we change the color scale at different times as the thermal becomes dilute and slows down over its evolution. We also plot the boundary of the thermal as computed in section 3. The laminar and turbulent thermals exhibit comparable entrainment, with the turbulent case entraining slightly more. The supplemental material includes movies of these two thermal simulations.

Fig. 1.

2D vertical slices at y = 0 of the (top) vertical velocity and (bottom) density of two thermals with different Reynolds numbers. The two thermals are initialized with the same noise field. At each time, the left panel shows the laminar thermal (Re = 630), and the right panel shows the turbulent thermal (Re = 6300). We use the same color scale for the two thermals at the same time; we change the color scale at different times as the thermal becomes dilute and slows down over its evolution. We also plot the boundary of the thermal as computed in section 3. The laminar and turbulent thermals exhibit comparable entrainment, with the turbulent case entraining slightly more. The supplemental material includes movies of these two thermal simulations.

Fig. 2.

2D vertical slices of ωϕ at the y thermal midpoint (defined in section 3) of two thermals with different Reynolds numbers, at t = 30. These slices show positive vorticity cores of the thermals, which coincide with the buoyancy anomaly. Note the much smaller scales of ωϕ in the turbulent case.

Fig. 2.

2D vertical slices of ωϕ at the y thermal midpoint (defined in section 3) of two thermals with different Reynolds numbers, at t = 30. These slices show positive vorticity cores of the thermals, which coincide with the buoyancy anomaly. Note the much smaller scales of ωϕ in the turbulent case.

In Fig. 2 we plot the azimuthal vorticity (ωϕ = ∂zur − ∂ruz, 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 zt, which we do using a density threshold. To calculate the thermal-top velocity wtop, we first perform a fit to zt, 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 wtop [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 zt 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.

Fig. 3.

Thermal-top height as a function of time in our 10 thermal simulations. Lower Reynolds numbers are shown in yellow, and higher Reynolds numbers are shown in blue. The dashed lines represent the t fit given in (A1). The fit is very good for times greater than ≈10.

Fig. 3.

Thermal-top height as a function of time in our 10 thermal simulations. Lower Reynolds numbers are shown in yellow, and higher Reynolds numbers are shown in blue. The dashed lines represent the t fit given in (A1). The fit is very good for times greater than ≈10.

To find the axisymmetric volume translating with this velocity, we calculate the Stokes streamfunction ψ associated with the azimuthally averaged velocities, uaxi = [uaxi(r, z), waxi(r, z)]. The Stokes streamfunction satisfies (Batchelor 2000)

 
×(ψeϕ)=(2πr)(uaxiwtopez),
(11)

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, ztop) where the axisymmetric vertical velocity matches the thermal-top velocity, that is, waxi(0, ztop) = wtop. 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)2dz 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 zt becomes large, but this is an artifact of the thermal interacting with the top boundary (see Fig. 3).

Fig. 4.

Thermal volume as a function of thermal-top height in our 10 thermal simulations. Laminar simulations are shown in yellow, and turbulent simulations are shown in blue. Note the use of (left) a linear scale and (right) a log scale (to illustrate that the volume increases like zt3). The decreases in volume at late times in our Laminar simulations are due to the thermal interacting with the top boundary.

Fig. 4.

Thermal volume as a function of thermal-top height in our 10 thermal simulations. Laminar simulations are shown in yellow, and turbulent simulations are shown in blue. Note the use of (left) a linear scale and (right) a log scale (to illustrate that the volume increases like zt3). The decreases in volume at late times in our Laminar simulations are due to the thermal interacting with the top boundary.

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 dlnV/dz to the net fractional entrainment rate εnet, giving

 
εnetdlnVdz=εδ.
(12)

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 dlnV/dz in our simulations. We define the thermal radius rth 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 Vwtop. 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 rth 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 rth in steps of 0.1, and then calculating the average εnet and average rth 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.

Fig. 5.

Net entrainment rate as a function of radius. The mean net entrainment rate across the thermals with a given Reynolds number is plotted with the dark thick line, and each individual simulation is plotted with a thin line. The range between simulations of a given Reynolds number is shaded. The entrainment rate decreases as r−1 (dashed line) over the course of the simulation.

Fig. 5.

Net entrainment rate as a function of radius. The mean net entrainment rate across the thermals with a given Reynolds number is plotted with the dark thick line, and each individual simulation is plotted with a thin line. The range between simulations of a given Reynolds number is shaded. The entrainment rate decreases as r−1 (dashed line) over the course of the simulation.

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 Md, and track the rate at which Md 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 zct 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 zt ~ rth, 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).

Fig. 6.

Average net entrainment and detrainment rates as a function of the height of the thermal. The solid lines show the net entrainment rate, and the dotted lines show the detrainment rate. Because the detrainment rate is typically at least an order of magnitude smaller than the net entrainment rate, we have εnetε. The dashed black line shows a z−1 power law.

Fig. 6.

Average net entrainment and detrainment rates as a function of the height of the thermal. The solid lines show the net entrainment rate, and the dotted lines show the detrainment rate. Because the detrainment rate is typically at least an order of magnitude smaller than the net entrainment rate, we have εnetε. The dashed black line shows a z−1 power law.

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 rth. 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.

Fig. 7.

Entrainment efficiency e [see (1)] as a function of thermal-top height. The mean efficiency across the thermals with a given Reynolds number is plotted with the dark thick line, and each individual simulation is plotted with a thin line. The range between simulations of a given Reynolds number is shaded. Although there is a wide spread in the entrainment efficiency in any individual simulation (especially for higher Reynolds number), the turbulent thermals typically have higher entrainment efficiency than the laminar thermals.

Fig. 7.

Entrainment efficiency e [see (1)] as a function of thermal-top height. The mean efficiency across the thermals with a given Reynolds number is plotted with the dark thick line, and each individual simulation is plotted with a thin line. The range between simulations of a given Reynolds number is shaded. Although there is a wide spread in the entrainment efficiency in any individual simulation (especially for higher Reynolds number), the turbulent thermals typically have higher entrainment efficiency than the laminar thermals.

The average entrainment efficiency drops quickly for heights zt > 16 for our laminar simulations and for heights zt > 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 zt = 6 and zt = 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 rO(1 km).

Table 1.

The average entrainment efficiency e [see (1)] in our simulations at two Reynolds numbers, and as reported in laboratory experiments. We also provide the equivalent n parameter described in S57 (see the “Analytic theory of thermals” subsection), as well as the initial aspect ratio.

The average entrainment efficiency e [see (1)] in our simulations at two Reynolds numbers, and as reported in laboratory experiments. We also provide the equivalent n parameter described in S57 (see the “Analytic theory of thermals” subsection), as well as the initial aspect ratio.
The average entrainment efficiency e [see (1)] in our simulations at two Reynolds numbers, and as reported in laboratory experiments. We also provide the equivalent n parameter described in S57 (see the “Analytic theory of thermals” subsection), as well as the initial aspect ratio.

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.

Fig. A1.

2D vertical slices at y = 0 of density perturbation of thermals at t = 30 for (top) laminar and (bottom) turbulent simulations. Each vertical pair of simulations has identical initial noise fields, but this noise field varies in the simulations arranged horizontally, leading to substantial variation in thermal evolution across simulations at a given Reynolds number. However, the laminar thermals are systematically higher and have larger density anomalies at this time than the turbulent thermals. We also plot the boundary of the thermal as computed in section 3.

Fig. A1.

2D vertical slices at y = 0 of density perturbation of thermals at t = 30 for (top) laminar and (bottom) turbulent simulations. Each vertical pair of simulations has identical initial noise fields, but this noise field varies in the simulations arranged horizontally, leading to substantial variation in thermal evolution across simulations at a given Reynolds number. However, the laminar thermals are systematically higher and have larger density anomalies at this time than the turbulent thermals. We also plot the boundary of the thermal as computed in section 3.

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 zt(t). At each time (we have full volume outputs every 1/100.36 time units), we first calculate the horizontal average of ρ′ at each height, ⟨ρ′⟩x,y. We define a cutoff value of 0.1maxzρ′⟩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 wtop, so care must be taken when calculating dzt/dt. We find the best way to calculate this derivative is to first fit

 
ztat+z0,
(B1)

based on (9), using a nonlinear least squares fit, and then take the derivative analytically to determine wtop.

Next we calculate the horizontal midpoint of the thermal. We define the midpoint using the vertical velocity. At each height, we define xm by

 
xm=w>0xww>0w,
(B2)

and similarly for ym. We then pick the vertical height which maximizes w>0w. The horizontal midpoint of the thermal is then (xm, ym) 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 waxi(r, z), we take the average of w in azimuthal rings of width Δr centered around r + Δr/2.

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

 
rψ=2πr(waxiwtop),
(B3)

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ρthV, where ρth is the average density perturbation of the thermal. Entrainment changes ρth and V individually but not their product. Thus, if we know the total “mass” of detrained fluid Md throughout the simulation, we may approximate the detrainment rate as

 
δ1MthdMddz.
(C1)

Defining Md 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 Md to be only the mass below zbot, 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 Md.) We also calculated the mass below zbotrth/4. Although there is less mass below that point than below zbot at any given time, it leads to similar detrainment rates. This suggests that almost all of the mass that falls below zbot will soon fall below zbotrth/4 (i.e., it has truly been detrained).

Figure C1 shows in thin lines Md for all of our simulations, normalized by the total initial mass, M0=(4π/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 zbot 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 zct bins of size 0.5).

Fig. C1.

The fraction of detrained mass, defined as mass below zbot, the lowest point of the thermal, as a function of thermal-top height. The blue curves show turbulent thermals, and the yellow curves show laminar thermals. The dark thick lines show the detrained mass averaged over our ensemble of turbulent or laminar thermals.

Fig. C1.

The fraction of detrained mass, defined as mass below zbot, the lowest point of the thermal, as a function of thermal-top height. The blue curves show turbulent thermals, and the yellow curves show laminar thermals. The dark thick lines show the detrained mass averaged over our ensemble of turbulent or laminar thermals.

In its early evolution (zt < 5, not shown), the thermal detrains a nonnegligible amount of fluid, roughly 5%–10% of M0. 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 zct = 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

 
tsD2s=us
(C2)

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 rth(θ) using linear interpolation. Then the dye field is initialized to

 
sinit(r,θ)=12{1erf[rrth(θ)Δr]},
(C3)

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 cd.

To quantify detrainment, we calculate how cd changes relative to the total amount of dye, ctot. The total dye is approximately equal to the volume of the thermal at the previous reinitialization time. We plot the fraction of detrained dye, cd/ctot, 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.

Fig. C2.

The fraction of detrained dye, defined as the dye below zbot, the lowest point of the thermal, as a function of thermal-top height. The dye field is reinitialized to be almost entirely within the thermal every 10 time units, denoted with a cross. The final detrained dye fraction before the next reinitialization is denoted with a circle. The yellow curves show the fraction for laminar thermals, and the blue curve for one turbulent thermal. The slope of the line connecting the cross and circle for each curve is the detrainment rate. We are unable to measure the detrainment rate for turbulent thermals using this method because it takes more than 10 time units for the detrained fluid to fall below zbot.

Fig. C2.

The fraction of detrained dye, defined as the dye below zbot, the lowest point of the thermal, as a function of thermal-top height. The dye field is reinitialized to be almost entirely within the thermal every 10 time units, denoted with a cross. The final detrained dye fraction before the next reinitialization is denoted with a circle. The yellow curves show the fraction for laminar thermals, and the blue curve for one turbulent thermal. The slope of the line connecting the cross and circle for each curve is the detrainment rate. We are unable to measure the detrainment rate for turbulent thermals using this method because it takes more than 10 time units for the detrained fluid to fall below zbot.

We then define the detrainment rate for laminar thermals to be

 
δ=1ctotΔcdΔzt,

where Δcd and Δzt 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 zt, we interpret our calculated value as corresponding to the height zt,0 + Δzt/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 zt > 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

REFERENCES
Abramowitz
,
M.
, and
I. A.
Stegun
,
1968
:
Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables. Dover Publications, 1046 pp
.
Akhmetov
,
D. G.
,
2009
:
Vortex Rings. Springer, 164 pp., https://doi.org/10.1007/978-3-642-05016-9
.
Anders
,
E. H.
,
D.
Lecoanet
, and
B. P.
Brown
,
2019
:
Entropy rain: Dilution and compression of thermals in stratified domains
.
Astrophys. J.
,
884
,
65
, https://doi.org/10.3847/1538-4357/ab3644.
Asai
,
T.
, and
A.
Kasahara
,
1967
:
A theoretical study of the compensating downward motions associated with cumulus clouds
.
J. Atmos. Sci.
,
24
,
487
496
, https://doi.org/10.1175/1520-0469(1967)024<0487:ATSOTC>2.0.CO;2.
Ascher
,
U. M.
,
S. J.
Ruuth
, and
R. J.
Spiteri
,
1997
:
Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations
.
Appl. Numer. Math.
,
25
,
151
167
, https://doi.org/10.1016/S0168-9274(97)00056-1.
Batchelor
,
G. K.
,
2000
:
An Introduction to Fluid Dynamics. Cambridge University Press, 658 pp
.
Blyth
,
A. M.
,
S. G.
Lasher-Trapp
, and
W. A.
Cooper
,
2005
:
A study of thermals in cumulus clouds
.
Quart. J. Roy. Meteor. Soc.
,
131
,
1171
1190
, https://doi.org/10.1256/qj.03.180.
Bond
,
D.
, and
H.
Johari
,
2010
:
Impact of buoyancy on vortex ring development in the near field
.
Exp. Fluids
,
48
,
737
745
, https://doi.org/10.1007/s00348-009-0761-z.
Burns
,
K. J.
,
G. M.
Vasil
,
J. S.
Oishi
,
D.
Lecoanet
, and
B.
Brown
,
2016
:
Dedalus: Flexible framework for spectrally solving differential equations. Astrophysics Source Code Library, http://www.ascl.net/1603.015
.
Burns
,
K. J.
,
G. M.
Vasil
,
J. S.
Oishi
,
D.
Lecoanet
, and
B.
Brown
,
2019
:
Dedalus: A flexible framework for numerical simulations with spectral methods. arXiv, https://arxiv.org/abs/1905.10388
.
Carpenter
,
R. L.
,
K. K.
Droegemeier
, and
A. M.
Blyth
,
1998
:
Entrainment and detrainment in numerically simulated cumulus congestus clouds. Part II: Cloud budgets
.
J. Atmos. Sci.
,
55
,
3433
3439
, https://doi.org/10.1175/1520-0469(1998)055<3433:EADINS>2.0.CO;2.
Cotton
,
W. R.
,
1975
:
Theoretical cumulus dynamics
.
Rev. Geophys.
,
13
,
419
448
, https://doi.org/10.1029/RG013i002p00419.
Craig
,
G. C.
, and
A.
Dörnbrack
,
2008
:
Entrainment in cumulus clouds: What resolution is cloud-resolving?
J. Atmos. Sci.
,
65
,
3978
3988
, https://doi.org/10.1175/2008JAS2613.1.
Damiani
,
R.
,
G.
Vali
, and
S.
Haimov
,
2006
:
The structure of thermals in cumulus from airborne dual-Doppler radar observations
.
J. Atmos. Sci.
,
63
,
1432
1450
, https://doi.org/10.1175/JAS3701.1.
Dawe
,
J. T.
, and
P. H.
Austin
,
2011
:
Interpolation of LES cloud surfaces for use in direct calculations of entrainment and detrainment
.
Mon. Wea. Rev.
,
139
,
444
456
, https://doi.org/10.1175/2010MWR3473.1.
Dawe
,
J. T.
, and
P. H.
Austin
,
2013
:
Direct entrainment and detrainment rate distributions of individual shallow cumulus clouds in an LES
.
Atmos. Chem. Phys.
,
13
,
7795
7811
, https://doi.org/10.5194/acp-13-7795-2013.
de Roode
,
S. R.
,
A. P.
Siebesma
,
H. J.
Jonker
, and
Y.
de Voogd
,
2012
:
Parameterization of the vertical velocity equation for shallow cumulus clouds
.
Mon. Wea. Rev.
,
140
,
2424
2436
, https://doi.org/10.1175/MWR-D-11-00277.1.
de Rooy
,
W. C.
, and
A. P.
Siebesma
,
2010
:
Analytical expressions for entrainment and detrainment in cumulus convection
.
Quart. J. Roy. Meteor. Soc.
,
136
,
1216
1227
, https://doi.org/10.1002/QJ.640.
de Rooy
,
W. C.
, and Coauthors
,
2013
:
Entrainment and detrainment in cumulus convection: An overview
.
Quart. J. Roy. Meteor. Soc.
,
139
,
1
19
, https://doi.org/10.1002/qj.1959.
Escudier
,
M. P.
, and
T.
Maxworthy
,
1973
:
On the motion of turbulent thermals
.
J. Fluid Mech.
,
61
,
541
552
, https://doi.org/10.1017/S0022112073000856.
Ferrier
,
B. S.
, and
R. A.
Houze
,
1989
:
One-dimensional time-dependent modeling of GATE cumulonimbus convection
.
J. Atmos. Sci.
,
46
,
330
352
, https://doi.org/10.1175/1520-0469(1989)046<0330:ODTDMO>2.0.CO;2.
Fohl
,
T.
,
1967
:
Turbulent effects in the formation of buoyant vortex rings
.
J. Appl. Phys.
,
38
,
4097
4098
, https://doi.org/10.1063/1.1709085.
Grabowski
,
W. W.
,
1993
:
Cumulus entrainment, fine-scale mixing, and buoyancy reversal
.
Quart. J. Roy. Meteor. Soc.
,
119
,
935
956
, https://doi.org/10.1002/qj.49711951305.
Grabowski
,
W. W.
,
1995
:
Entrainment and mixing in buoyancy-reversing convection with applications to cloud-top entrainment instability
.
Quart. J. Roy. Meteor. Soc.
,
121
,
231
253
, https://doi.org/10.1002/qj.49712152202.
Grabowski
,
W. W.
,
T. L.
Clark
,
W. W.
Grabowski
, and
T. L.
Clark
,
1991
:
Cloud-environment interface instability: Rising thermal calculations in two spatial dimensions
.
J. Atmos. Sci.
,
48
,
527
546
, https://doi.org/10.1175/1520-0469(1991)048<0527:CIIRTC>2.0.CO;2.
Grabowski
,
W. W.
,
T. L.
Clark
,
W. W.
Grabowski
, and
T. L.
Clark
,
1993a
:
Cloud-environment interface instability: Part II: Extension to three spatial dimensions
.
J. Atmos. Sci.
,
50
,
555
573
, https://doi.org/10.1175/1520-0469(1993)050<0555:CEIIPI>2.0.CO;2.
Grabowski
,
W. W.
,
T. L.
Clark
,
W. W.
Grabowski
, and
T. L.
Clark
,
1993b
:
Cloud-environment interface instability. Part III: Direct influence of environmental shear
.
J. Atmos. Sci.
,
50
,
3821
3828
, https://doi.org/10.1175/1520-0469(1993)050<3821:CEIIPI>2.0.CO;2.
Held
,
I. M.
,
2005
:
The gap between simulation and understanding in climate modeling
.
Bull. Amer. Meteor. Soc.
,
86
,
1609
1614
, https://doi.org/10.1175/BAMS-86-11-1609.
Hernandez-Deckers
,
D.
, and
S. C.
Sherwood
,
2016
:
A numerical investigation of cumulus thermals
.
J. Atmos. Sci.
,
73
,
4117
4136
, https://doi.org/10.1175/JAS-D-15-0385.1.
Hernandez-Deckers
,
D.
, and
S. C.
Sherwood
,
2018
:
On the role of entrainment in the fate of cumulus thermals
.
J. Atmos. Sci.
,
75
,
3911
3924
, https://doi.org/10.1175/JAS-D-18-0077.1.
Heus
,
T.
,
H. J. J.
Jonker
,
H. E. A.
Van den Akker
,
E. J.
Griffith
,
M.
Koutek
, and
F. H.
Post
,
2009
:
A statistical approach to the life cycle analysis of cumulus clouds selected in a virtual reality environment
.
J. Geophys. Res.
,
114
,
D06208
, https://doi.org/10.1029/2008JD010917.
Houghton
,
H. G.
, and
H. E.
Cramer
,
1951
:
A theory of entrainment in convective currents
.
J. Meteor.
,
8
,
95
102
, https://doi.org/10.1175/1520-0469(1951)008<0095:ATOEIC>2.0.CO;2.
Jeevanjee
,
N.
, and
D. M.
Romps
,
2016
:
Effective buoyancy at the surface and aloft
.
Quart. J. Roy. Meteor. Soc.
,
142
,
811
820
, https://doi.org/10.1002/qj.2683.
Jeevanjee
,
N.
,
P.
Hassanzadeh
,
S.
Hill
, and
A.
Sheshadri
,
2017
:
A perspective on climate model hierarchies
.
J. Adv. Model. Earth Syst.
,
9
,
1760
1771
, https://doi.org/10.1002/2017MS001038.
Johari
,
H.
,
1992
:
Mixing in thermals with and without buoyancy reversal
.
J. Atmos. Sci.
,
49
,
1412
1426
, https://doi.org/10.1175/1520-0469(1992)049<1412:MITWAW>2.0.CO;2.
Klocke
,
D.
,
R.
Pincus
, and
J.
Quaas
,
2011
:
On constraining estimates of climate sensitivity with present-day observations through model weighting
.
J. Climate
,
24
,
6092
6099
, https://doi.org/10.1175/2011JCLI4193.1.
Kuo
,
H.
,
1962
:
On the controlling influences of eddy diffusion on thermal convection
.
J. Atmos. Sci.
,
19
,
236
243
, https://doi.org/10.1175/1520-0469(1962)019<0236:OTCIOE>2.0.CO;2.
Lai
,
A. C. H.
,
B.
Zhao
,
A. W. K.
Law
, and
E. E.
Adams
,
2015
:
A numerical and analytical study of the effect of aspect ratio on the behavior of a round thermal
.
Environ. Fluid Mech.
,
15
,
85
108
, https://doi.org/10.1007/s10652-014-9362-3.
Langhans
,
W.
,
K.
Yeo
, and
D. M.
Romps
,
2015
:
Lagrangian investigation of the precipitation efficiency of convective clouds
.
J. Atmos. Sci.
,
72
,
1045
1062
, https://doi.org/10.1175/JAS-D-14-0159.1.
Levine
,
J.
,
1959
:
Spherical vortex theory of bubble-like motion in cumulus clouds
.
J. Meteor.
,
16
,
653
662
, https://doi.org/10.1175/1520-0469(1959)016<0653:SVTOBL>2.0.CO;2.
Lilly
,
D. K.
,
1962
:
On the numerical simulation of buoyant convection
.
Tellus
,
14
,
148
172
, https://doi.org/10.3402/tellusa.v14i2.9537.
Malkus
,
J. S.
,
1954
:
Some results of a trade-cumulus cloud investigation
.
J. Meteor.
,
11
,
220
237
, https://doi.org/10.1175/1520-0469(1954)011<0220:SROATC>2.0.CO;2.
Malkus
,
J. S.
, and
R. S.
Scorer
,
1955
:
The erosion of cumulus towers
.
J. Atmos. Sci.
,
12
,
43
57
, https://doi.org/10.1175/1520-0469(1955)012<0000:TEOCT>2.0.CO;2.
Mapes
,
B.
, and
R.
Neale
,
2011
:
Parameterizing convective organization to escape the entrainment dilemma
.
J. Adv. Model. Earth Syst.
,
3
,
M06004
, https://doi.org/10.1029/2011MS000042.
McKim
,
B.
,
N.
Jeevanjee
, and
D.
Lecoanet
,
2019
:
Buoyancy-driven entrainment in dry thermals. arXiv.org, https://arxiv.org/abs/1906.07224
.
Miller
,
L. J.
,
J. E.
Dye
, and
B. E.
Martner
,
1983
:
Dynamical-microphysical evolution of a convective storm in a weakly-sheared environment. Part II: Airflow and precipitation trajectories from Doppler radar observations
.
J. Atmos. Sci.
,
40
,
2097
2109
, https://doi.org/10.1175/1520-0469(1983)040<2097:DMEOAC>2.0.CO;2.
Morrison
,
H.
,
2017
:
An analytic description of the structure and evolution of growing deep cumulus updrafts
.
J. Atmos. Sci.
,
74
,
809
834
, https://doi.org/10.1175/JAS-D-16-0234.1.
Morrison
,
H.
, and
J. M.
Peters
,
2018
:
Theoretical expressions for the ascent rate of moist deep convective thermals
.
J. Atmos. Sci.
,
75
,
1699
1719
, https://doi.org/10.1175/JAS-D-17-0295.1.
Morton
,
B. R.
,
G.
Taylor
, and
J. S.
Turner
,
1956
:
Turbulent gravitational convection from maintained and instantaneous sources
.
Proc. Royal Soc. London
,
234A
,
1
23
, https://doi.org/10.1098/rspa.1956.0011.
Moser
,
D. H.
, and
S.
Lasher-Trapp
,
2017
:
The influence of successive thermals on entrainment and dilution in a simulated cumulus congestus
.
J. Atmos. Sci.
,
74
,
375
392
, https://doi.org/10.1175/JAS-D-16-0144.1.
Murphy
,
J. M.
, and Coauthors
,
2004
:
Quantification of modelling uncertainties in a large ensemble of climate change simulations
.
Nature
,
430
,
768
772
, https://doi.org/10.1038/nature02771.
Ogura
,
Y.
,
1962
:
Convection of isolated masses of a buoyant fluid: A numerical calculation
.
J. Atmos. Sci.
,
19
,
492
502
, https://doi.org/10.1175/1520-0469(1962)019<0492:COIMOA>2.0.CO;2.
Oishi
,
J. S.
,
B. P.
Brown
,
K. J.
Burns
,
D.
Lecoanet
, and
G. M.
Vasil
,
2018
:
Perspectives on reproducibility and sustainability of open-source scientific software from seven years of the Dedalus project. arXiv, https://arxiv.org/abs/1801.08200
.
Raymond
,
D.
, and
A.
Blyth
,
1986
:
A stochastic mixing model for nonprecipitating cumulus clouds
.
J. Atmos. Sci.
,
43
,
2708
2718
, https://doi.org/10.1175/1520-0469(1986)043<2708:ASMMFN>2.0.CO;2.
Richards
,
J. M.
,
1961
:
Experiments on the penetration of an interface by buoyant thermals
.
J. Fluid Mech.
,
11
,
369
384
, https://doi.org/10.1017/S0022112061000585.
Romps
,
D. M.
,
2010
:
A direct measure of entrainment
.
J. Atmos. Sci.
,
67
,
1908
1927
, https://doi.org/10.1175/2010JAS3371.1.
Romps
,
D. M.
,
2016
:
The Stochastic Parcel Model: A deterministic parameterization of stochastically entraining convection
.
J. Adv. Model. Earth Syst.
,
8
,
319
344
, https://doi.org/10.1002/2015MS000537.
Romps
,
D. M.
, and
A. B.
Charn
,
2015
:
Sticky thermals: Evidence for a dominant balance between buoyancy and drag in cloud updrafts
.
J. Atmos. Sci.
,
72
,
2890
2901
, https://doi.org/10.1175/JAS-D-15-0042.1.
Sànchez
,
O.
,
D. J.
Raymond
,
L.
Libersky
, and
A. G.
Petschek
,
1989
:
The development of thermals from rest
.
J. Atmos. Sci.
,
46
,
2280
2292
, https://doi.org/10.1175/1520-0469(1989)046<2280:TDOTFR>2.0.CO;2.
Saunders
,
P. M.
,
1961
:
An observational study of cumulus
.
J. Meteor.
,
18
,
451
467
, https://doi.org/10.1175/1520-0469(1961)018<0451:AOSOC>2.0.CO;2.
Scorer
,
R. S.
,
1957
:
Experiments on convection of isolated masses of buoyant fluid
.
J. Fluid Mech.
,
2
,
583
594
, https://doi.org/10.1017/S0022112057000397.
Scorer
,
R. S.
, and
F. H.
Ludlam
,
1953
:
Bubble theory of penetrative convection
.
Quart. J. Roy. Meteor. Soc.
,
79
,
94
103
, https://doi.org/10.1002/qj.49707933908.
Shariff
,
K.
, and
A.
Leonard
,
1992
:
Vortex rings
.
Annu. Rev. Fluid Mech.
,
24
,
235
279
, https://doi.org/10.1146/annurev.fl.24.010192.001315.
Sherwood
,
S. C.
,
D.
Hernández-Deckers
,
M.
Colin
, and
F.
Robinson
,
2013
:
Slippery thermals and the cumulus entrainment paradox
.
J. Atmos. Sci.
,
70
,
2426
2442
, https://doi.org/10.1175/JAS-D-12-0220.1.
Shivamoggi
,
B. K.
,
2010
:
Hydrodynamic impulse in a compressible fluid
.
Phys. Lett.
,
374A
,
4736
4740
, https://doi.org/10.1016/j.physleta.2010.09.062.
Simpson
,
J.
,
1983a
:
Cumulus clouds: Early aircraft observations and entrainment hypotheses. Mesoscale Meteorology—Theories, Observations and Models, D. K. Lilly and T. Gal-Chen, Eds., Springer, 355–373
.
Simpson
,
J.
,
1983b
:
Cumulus clouds: Interactions between laboratory experiments and observations as foundations for models. Mesoscale Meteorology—Theories, Observations and Models, D. K. Lilly and T. Gal-Chen, Eds., Springer, 399–412
.
Simpson
,
J.
, and
V.
Wiggert
,
1969
:
Models of precipitating cumulus towers
.
Mon. Wea. Rev.
,
97
,
471
489
, https://doi.org/10.1175/1520-0493(1969)097<0471:MOPCT>2.3.CO;2.
Simpson
,
J.
,
R. H.
Simpson
,
D. A.
Andrews
, and
M. A.
Eaton
,
1965
:
Experimental cumulus dynamics
.
Rev. Geophys.
,
3
,
387
431
, https://doi.org/10.1029/RG003i003p00387.
Squires
,
P.
, and
J. S.
Turner
,
1962
:
An entraining jet model for cumulo-nimbus updraughts
.
Tellus
,
14
,
422
434
, https://doi.org/10.3402/tellusa.v14i4.9569.
Stiller
,
O.
, and
G. C.
Craig
,
2001
:
A scaling hypothesis for moist convective updraughts
.
Quart. J. Roy. Meteor. Soc.
,
127
,
1551
1570
, https://doi.org/10.1002/qj.49712757505.
Stirling
,
A. J.
, and
R. A.
Stratton
,
2012
:
Entrainment processes in the diurnal cycle of deep convection over land
.
Quart. J. Roy. Meteor. Soc.
,
138
,
1135
1149
, https://doi.org/10.1002/qj.1868.
Stommel
,
H. M.
,
1947
:
Entrainment of air into a cumulus cloud
.
J. Meteor.
,
4
,
91
94
, https://doi.org/10.1175/1520-0469(1947)004<0091:EOAIAC>2.0.CO;2.
Tarshish
,
N.
,
N.
Jeevanjee
, and
D.
Lecoanet
,
2018
:
Buoyant motion of a turbulent thermal
.
J. Atmos. Sci.
,
75
,
3233
3244
, https://doi.org/10.1175/JAS-D-17-0371.1.
Turner
,
J. S.
,
1957
:
Buoyant vortex rings
.
Proc. Roy. Soc. London
,
239A
,
61
75
, https://doi.org/10.1098/rspa.1957.0022.
Turner
,
J. S.
,
1962
:
The ‘starting plume’ in neutral surroundings
.
J. Fluid Mech.
,
13
,
356
368
, https://doi.org/10.1017/S0022112062000762.
Turner
,
J. S.
,
1964a
:
The dynamics of spheroidal masses of buoyant fluid
.
J. Fluid Mech.
,
19
,
481
490
, https://doi.org/10.1017/S0022112064000854.
Turner
,
J. S.
,
1964b
:
The flow into an expanding spherical vortex
.
J. Fluid Mech.
,
18
,
195
208
, https://doi.org/10.1017/S0022112064000155.
Turner
,
J. S.
,
1986
:
Turbulent entrainment: the development of the entrainment assumption, and its application to geophysical flows
.
J. Fluid Mech.
,
173
,
431
471
, https://doi.org/10.1017/S0022112086001222.
Woodward
,
B.
,
1959
:
The motion in and around isolated thermals
.
Quart. J. Roy. Meteor. Soc.
,
85
,
144
151
, https://doi.org/10.1002/qj.49708536407.
Yano
,
J. I.
,
2014
:
Basic convective element: Bubble or plume? A historical review
.
Atmos. Chem. Phys.
,
14
,
7019
7030
, https://doi.org/10.5194/acp-14-7019-2014.
Yeo
,
K.
, and
D. M.
Romps
,
2013
:
Measurement of convective entrainment using Lagrangian particles
.
J. Atmos. Sci.
,
70
,
266
277
, https://doi.org/10.1175/JAS-D-12-0144.1.
Zhao
,
B.
,
A. W.
Law
,
A. C.
Lai
, and
E. E.
Adams
,
2013
:
On the internal vorticity and density structures of miscible thermals
.
J. Fluid Mech.
,
722
,
1
12
, https://doi.org/10.1017/jfm.2013.158.
Zhao
,
M.
,
2014
:
An investigation of the connections among convection, clouds, and climate sensitivity in a global climate model
.
J. Climate
,
27
,
1845
1862
, https://doi.org/10.1175/JCLI-D-13-00145.1.
Zhao
,
M.
, and
P. H.
Austin
,
2005
:
Life cycle of numerically simulated shallow cumulus clouds. Part II: Mixing dynamics
.
J. Atmos. Sci.
,
62
,
1291
1310
, https://doi.org/10.1175/JAS3415.1.

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/κ, whereas the entrainment time scale is ~(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.

Supplemental Material