Buoyancy-driven entrainment in dry thermals

Over 50 years ago it was proposed that dry thermals entrain because of buoyancy (via a constraint which requires an increase in the radius a ). However, this runs counter to the scaling arguments commonly used to derive the entrainment rate, which rely on either the self-similarity or a turbulent entrainment hypothesis. The assumption of turbulence-driven entrainment has been investigated and it has been found that the entrainment efficiency e varies by less than 20 % between laminar ( Re = 630) and turbulent ( Re = 6300) thermals. This motivated us to utilize the argument of buoyancy-controlled entrainment in addition to the ther-mal's vertical momentum equation to build a model for thermal dynamics which does not invoke turbulence or self-similarity. We derive simple expressions for the thermals' kinematic properties and their fractional entrainment rate 𝜖 and find close quantitative agreement with the values in direct numerical simulations. In particular, our expression for entrainment rate is consistent with the parametrization 𝜖 ∼ B ∕ w 2 , for Archimedean buoyancy B and vertical velocity w . We also directly validate the role of buoyancy-driven entrainment by running simulations where gravity is turned off midway through a thermal's rise. The entrainment efficiency e is observed to drop to less than one third of its original value in both the laminar and turbulent cases when g = 0, affirming the central role of buoyancy in entrainment for dry thermals.


INTRODUCTION
Since the introduction of the entrainment hypothesis by Morton et al. (1956), the physical mechanism of entrainment has been thought to originate from turbulent eddies. The wide acceptance of this theory is due to its success in analyzing flows such as jets and plumes in a large variety of physical contexts (Turner, 1986). The entrainment hypothesis states that the average inflow velocity across the edge of a turbulent flow is assumed to be proportional to a characteristic velocity. However, there are some flows where entrainment persists even in the laminar regime. For instance, "thermals", or regions of isolated buoyant fluid thought to be the basic unit of convection (Yano, 2014;Romps and Charn, 2015), were recently found to vary in entrainment by less than 20% in laminar (Reynolds number Re = 630) and turbulent (Re = 6300) simulations in a dry atmosphere (Lecoanet and Jeevanjee, 2019).
The observed expansion rates in thermals is much larger than that of plumes and is sensitive to the initial conditions (Turner, 1973;1986). According to Turner (1957), these properties may be better understood if the thermal is modelled as a buoyant vortex ring. Turner (1957) shows why a buoyant ring must expand with time due to momentum conservation, and an illuminating, complementary perspective was provided by Zhao et al. (2013) which shows that a ring's expansion is related to its baroclinicity. We refer the reader to Zhao et al. (2013) for a more in-depth explanation.
Here, we combine the buoyant vortex ring argument of Turner (1957) with the thermal's vertical momentum equation by regarding the thermal as the region of fluid that moves coherently with the vortex ring (i.e., as the ellipsoidal region of fluid whose average velocity equals the gross translational velocity of the vortex ring itself). This aligns with the notion of a vortex ring "bubble" (Shariff and Leonard, 1992;Akhmetov, 2009), and we assume that the radius of the thermal a is proportional to the radius of the vortex R, i.e., a = R. Utilizing this connection, we derive analytical expressions for the density, vertical velocity, height, and radius as functions of time. We then relate the radius to the height to get an expression for the fractional entrainment rate, ≡ d log V∕dz = e∕a, where e is a proportionality constant called the "entrainment efficiency" and a is the radius of the thermal 1 . This is significant because our predictions require no tuning parameters to fit the data, unlike previous works (Morton et al. 1956;Levine, 1959;Turner, 1962;1964;Simpson and Wiggert, 1969;Escudier and Maxworthy, 1973;Simpson, 1983;Turner, 1986;Johari, 1992). We validate this model by comparing the predictions of the radius, height, density, vertical velocity, and fractional entrainment to direct numerical simulations of dry thermals.
Our model does not invoke turbulence, contrary to the substantial amount of literature which attributes entrainment in convection to turbulent eddies (Morton et al. 1956(Morton et al. , 1956Turner, 1964;Escudier and Maxworthy, 1973;Baker et al. 1984;Ferrier and Houze, 1989;Craig and Dörnbrack, 2008;de Rooy and Siebesma, 2010;Romps and Kuang, 2010;de Rooy et al. 2013;Sherwood et al. 2013;Yano, 2014). In addition to validating our simple non-turbulent model against simulation data, we directly test the role of buoyancy by performing numerical simulations where gravity is removed midway through a thermal's rise. The thermals' observed entrainment rates decrease significantly. 1 This is equivalent to the mass flux M formulation of entrainment ≡ d log M∕dz when the fluid is Boussinesq and detrainment is assumed negligible (Lecoanet and Jeevanjee, 2019).
Our expressions and simulations seem to indicate that entrainment in dry thermals is predominantly a laminar process, consistent with Sànchez et al. (1989) and contravening the original entrainment hypothesis. Furthermore, buoyancy-driven entrainment is distinct from the "dynamic" entrainment proposed in Houghton and Cramer (1951) and used in Asai andKasahara (1967), Ferrier andHouze (1989), de Rooy and Siebesma (2010), and Morrison (2017), since dynamic entrainment is a consequence of positive vertical acceleration, whereas our theory and simulations show negative vertical acceleration. We conclude with a discussion of the implications of these findings on future studies of entrainment, given the importance of this process in understanding cumulus convection and climate sensitivity (Klocke et al. 2011;Zhao, 2014).
We begin by reviewing previous work on buoyant vortex rings in Section 2 and discuss the relationship between a vortex ring and a thermal. In Section 3 we introduce a new model which allows us to derive the characteristics of a thermal and its entrainment. In Section 4 we verify the predictions made by the model and then affirm the role of buoyancy in entrainment. In Section 5 we discuss second-order contributions to entrainment and conclude by looking at future research directions.

Fundamental fluid mechanics
In our analysis of thermals we will use the Boussinesq approximation to the momentum equation, where density differences ′ in the flow are small compared to the constant background reference density . We decompose the density as = + ′ . The pressure can be decomposed similarly as p = p+p ′ where p is in hydrostatic balance with , ie. ∇p = g, where g = −gẑ is the gravitational acceleration. Formally, the Boussinesq approximation is valid when vertical length-scales in the problem are smaller than the scale height. This approximation is helpful, as it eliminates the effect of adiabatic expansion of a fluid parcel as it rises. Therefore, the observed expansion of thermals will be purely due to the mechanical effect of entrainment. The integrated buoyancy force in the Boussinesq approximation is F = − ∫ ′ g dV, where the integral is taken over the region of interest. The momentum equation is where u is the velocity field and is the kinematic viscosity. When a spherical thermal is released from rest, it accrues vorticity from its baroclinicity and generates a vortex ring structure, where density perturbations ′ are largely confined to the vortex core that develops ( Figure  1). Buoyancy is present only where there are density perturbations, therefore, the buoyancy resides mostly within the vortex core. In our model, we will assume that the buoyancy resides entirely within the vortex core. This assumption will be tested in Section 4. We introduce some basic aspects of vortex rings below.

Vortex ring dynamics
A vortex ring can be characterized by its kinematic properties such as circulation and impulse. The circulation of a vortex ring can be calculated along a circuit  passing through the centre of the thermal's vortex ring and then returning through the ambient fluid, as in Figure 1. The integral can also be computed as an area integral of vorticity over the region bounded by the circuit, . For vortex rings, which are azimuthally symmetric, the vorticity is =̂and this integral simplifies in cylindrical coordinates (r, , z), where the origin is located in the centre of the vortex ring, The impulse for incompressible fluids is defined as (Lamb and Caflisch, 1993;Batchelor, 2000;Akhmetov, 2009;Shivamoggi, 2010), where x is the position vector and Vol is the entire domain. It can be interpreted as the time and volume integral of external forces that must be applied to a flow in order to generate the observed fluid motion from rest (Lamb and Caflisch, 1993). Therefore, internal forces such as pressure or viscosity do not come into play. Equation (2) simplifies if we assume the flow is Boussinesq, that the radius of the core is much smaller than the radius of the ring R, and use the second part of Equation (1), We only get an impulse in the z direction, which is a consequence of the azimuthal symmetry of the flow. We turn to the circulation equation to consider how circulation evolves in time, F I G U R E 1 2D vertical slices at the thermal midpoint of the density perturbation (top) and vorticity (bottom) of two thermals with different Reynolds numbers at = 2.5. The left (right) panel shows the Re = 630 (Re = 6300) simulations. The same colour scale is used for both Re. The thermals originally start as a buoyant sphere (with noise added to break the symmetries of the problem), but then develop into a vortex ring that induces a flow structure that moves coherently with it. We identify this region as the thermal, and details of how it is computed is given in section 3 of Lecoanet and Jeevanjee (2019). The contour of the thermal is plotted in black. Both the density perturbations and the vorticity become primarily concentrated within the core region of the vortex in the laminar thermal. This is also true in the turbulent thermal, but to a lesser extent due to vorticity fluctuations at many scales, characteristic of a turbulent fluid. The dashed blue line is an example of a circuit  which passes through no buoyant fluid and can be used for Equations (1) and (3) In our current analysis we will assume the Reynolds number is sufficiently high that the viscous effects can be neglected (though Section 5.1 gives more discussion of viscous torques). If the circuit  passes through no buoyant fluid, then the circulation is constant with time (Fohl, 1967). This circuit will be possible once the vortex ring has formed, so we expect the circulation to initially increase with time as the thermal "spins up" but then reach a constant value. We will verify this property of the circulation in Section 4.

2.3
Expansion by conservation of momentum Turner (1957)'s argument applies after the vortex ring reaches a constant circulation and can be summarized succinctly: internal forces such as pressure or viscosity do not affect the impulse of the ring, so the change in impulse depends only on buoyancy.
F and are constants and therefore R must increase with time, i.e. the thermal must entrain. Anders et al. (2019) give a generalization of this argument to the density-stratified case where is not constant.

Expansion by baroclinicity
Although the argument in Equation (5) shows why the radius of thermal increases with time, it is helpful to look at the vorticity equation to explain how. We follow Zhao et al. (2013) below. The vorticity equation in the Boussinesq approximation is Once again ignoring viscous effects, in cylindrical coordinates we have The first term represents the intensification of vorticity due to stretching of vortex lines (Thorne and Blandford, 2017), but should disappear when integrated over the entire thermal because of the symmetry of the field (Figures 1 and 2). Equation (7) then tells us that the vorticity evolution is dictated by the second term, baroclinicity, which depends on gradients in density perturbations. The gradients are always set up such that vorticity is constantly F I G U R E 2 Azimuthally averaged vertical slices of the change in vorticity ⟨D ∕Dt⟩ at = 2.5. Angle brackets denote an azimuthal average has been taken, which is the reason for the symmetry of the figures across x = 0. Radial gradients in the density field produce a torque (baroclinicity) that drives the continual destruction of vorticity in the centre of the vortex and the continual creation vorticity in the outer region. This causes the centre of rotation of the vortex to continually move outward (Zhao et al. 2013) and is the mechanical description of how dry thermals expand and therefore entrain as they rise being created in the exterior of the ring and destroyed in the interior, as shown in Figure 2, and thus the core region of the ring appears to continually move radially outward. As the ring expands, it entrains more fluid. This will be shown in Section 4 to be the dominant mechanism of entrainment.

A NEW MODEL FOR THERMALS
Given the above picture for buoyancy-driven entrainment in thermals, we endeavour to build an analytical model for thermals which captures this picture, inspired by the formulation of Escudier and Maxworthy (1973) who invoke the vertical momentum equation. In particular, we assume: 1. The flow is Boussinesq and begins as a buoyant sphere at t = 0, but then fully spins up into a vortex ring at t = t o with radius R o , height z o , and vertical velocity w o . 2. The impulse of a buoyant vortex ring is given by Equation (3), with circulation Γ and background density being constant.
3. The only external force imparting an impulse is the integrated buoyancy F, which is constant because detrainment is assumed to be negligible (Lecoanet and Jeevanjee, 2019). The gross entrainment does not affect F because the surrounding fluid has no mass anomaly. While ⟨ ′ ⟩ will change and V will change, their combination ⟨ ′ ⟩V will remain constant. For a thermal with volume V = ma 3 , where ⟨ ′ ⟩ is the average density anomaly of the thermal and m is a constant for the volume since the thermal is no longer a sphere. m is verified to be roughly constant in Anders et al. (2019). 4. The thermal is the region of fluid that moves coherently with the vortex ring and so the thermal's radius a is proportional to the vortex radius R, i.e., a = R. 5. The impulse of a thermal with volume ma 3 and vertical velocity w can be written as (Akhmetov, 2009), where we parametrize virtual mass effects via a virtual mass coefficient (1 + k) (Batchelor, 2000). For all ellipsoidal fluid geometries, the value of k is the same as the solid body value of the same geometry Tarshish et al. (2018).
These assumptions will be discussed in more detail through the rest of the paper and will be tested along with the model predictions in Section 4.
We begin our analysis with the vortex momentum Equation (3) which is equal to the impulse imparted by the integrated buoyancy, Ft. Assuming that F is constant, we evaluate this equality at t and t o and take their ratio to get the radius of the vortex ring as a function of time, We utilize the proportionality of a and R and introduce a nondimensional time, = t∕t o to get the thermal's radius. a Now consider the equation for integrated buoyancy conservation, Equation (8). If we plug in Equation (11) for a, we can solve for ⟨ ′ ⟩ explicitly as a function of time, where ⟨ ′ o ⟩ is the average density anomaly of the thermal at t o . Now we turn to the Boussinesq momentum equation for a spherical thermal, Equation (9). In this formulation, a thermal is the region of fluid moving together with the vortex, so the impulse will be equal to the sum of the vortex momentum and the momentum created by virtual mass (Akhmetov, 2009). Evaluating Equation (9) at t and t o and taking their ratio gives the velocity of the thermal as a function of time: Now we integrate w to get z.
Assuming that detrainment is negligible implies that integrated buoyancy is conserved. Therefore the fractional entrainment rate is also the fractional change in thermal volume with respect to height (Lecoanet and Jeevanjee, 2019). Straightforward algebra shows Expressed as a function of the thermal's radius, This is one of the main results of this paper, along with having an analytical model for all the thermal's variables (a, w, ⟨ ′ ⟩) which does not invoke similarity or turbulence. Unlike Escudier and Maxworthy (1973), our approach does not require any additional empirical constants to determine the thermal's behavior. This is because Escudier and Maxworthy (1973) use the turbulent entrainment assumption, d( a 3 )∕dt = 3 a 2 |w|, ( is an unspecified constant tuned to match experiments), instead of Equation (3). They now have four unknowns (a, w, ⟨ ′ ⟩, ) and only three equations. In our case however, we employ Equation (3) instead of an entrainment assumption to fully solve the system for a, w, and ⟨ ′ ⟩.
Our approach also gives a simple expression for the thermal's entrainment efficiency, We can also use Equations (3) and (9), and I z = Ft to rewrite the efficiency in terms of the thermal's integrated buoyancy and circulation. This yields another key result, namely that e is a constant dictated by the initial spin-up of the thermal: Thermals with different initial conditions will thus have different efficiencies, which perhaps explains the variance of efficiencies with respect to Re found in Lecoanet and Jeevanjee (2019). Equation (17) suggests in particular that differences in e between the laminar and turbulent cases in that paper may stem from differences in Γ (Figure 3). Furthermore, variations in initial aspect ratio, which Lai et al. (2015) showed could explain the large range of thermal entrainment rates found in the literature, might also be understood via Equation (17) in terms of their effect on F and Γ.
Another important aspect of Equation (17) is that it embodies the well-known scaling e ∼ F∕ Γ 2 (e.g., Turner, 1957;Fohl, 1967;Bond and Johari, 2010;Lai et al. 2015). If we apply this scaling to the entrainment rate and write it using Γ ∼ wa and F ∼ Ba 3 (where B = g⟨ ′ ⟩∕ is the average Archimedean buoyancy), we find Though this scaling is not often applied to atmospheric convection, it was recently employed (on largely dimensional grounds) in the parametrization of Tan et al. (2018), and was also used in Gregory (2001). Equation (17) provides a precise foundation for Equation (18), at least in the idealized dry case.

Simulation set-up
We analyze the direct numerical simulations of thermals found in Lecoanet and Jeevanjee (2019) 2 . We outline the simulation set-up below, but for more details (including the scalings used for non-dimensionalization)are given in section 2 of Lecoanet and Jeevanjee (2019).
Simulations are run with Dedalus, an open source pseudo-spectral framework (Burns et al. 2019). We solve the non-dimensionalized Boussinesq equations: After non-dimensionalization, the simulations are entirely described by the initial condition, the Reynolds number 2 Although they ran an ensemble of simulations, we restrict our analysis to the simulations with the fifth initial condition.

F I G U R E 3 Verification of constant circulation Γ for laminar
(yellow) and turbulent (purple) simulations. For both cases, Γ is roughly constant, which validates the approximation used to derive Equation (11). The circulation also shows the thermals' spin-up periods. We chose our time-scale such that = 1 is approximately the spin-up time and the Prandtl number, Pr. We take Pr = 1 for all simulations, Re = 630 for laminar runs, and Re = 6300 for turbulent runs. We initialize the thermal with a spherical density anomaly ′ of magnitude negative 1 and diameter L th , plus a noise field to break symmetries in the problem. The spherical anomaly is placed near the bottom of a 3D domain with height 20L th and horizontal lengths 10L th . The simulation time ranges from t ∈ [0, 63.2], long enough for the thermal to approach or hit the top boundary. We utilize the thermal tracking algorithm described in section 3 of Lecoanet and Jeevanjee (2019). This defines the thermal to be the axisymmetric volume whose averaged vertical velocity matches the velocity of the thermal's top, where the thermal top is defined by a ′ threshold. Utilizing this tracking method, we can calculate the height, radius, velocity, and other dynamic quantities of the thermal. These measurements will serve as the test for the equations developed in Section 3. All comparisons use a spin-up time of t o = 4 √ 10 ≈ 12.6, which is determined by looking at when the thermal's circulation reaches a constant value (Figure 3). All of the code used to analyze the simulations in this work can be found online in a Github repository (https://github.com/mckimb/buoyant_entrainment; accessed 1 November 2019).

Comparison to simulations
The equation for the radius of the vortex ring (Equation (10)) depended on an assumption of constant circulation Γ within the thermal. While the constancy of circulation in Boussinesq thermals was confirmed in tank experiments in Zhao et al. (2013), we thought it worthwhile to also confirm this property in our F I G U R E 4 Verification of impulse expressions. The laminar thermal's impulse is shown in yellow and the turbulent's impulse in purple. Equation (2) was derived assuming the radius of the core region of a vortex is small compared to its overall radius. We test the validity of this assumption by explicitly comparing it to the impulse computed by Equation (2). We also assume integrated buoyancy is constant and therefore the impulse generated by it, F , can be calculated too. We find that in both laminar and turbulent thermals, F tracks Equation (2) closely.
simulations. We obtain the circulation by first calculating the azimuthal vorticity field, selectively choosing only the data points which lie within the thermal boundary, and then azimuthally averaging the data. We then take an area integral over the thermal's domain which is equivalent to the integral in the second part of Equation (1). The resulting Γ is fairly constant in the laminar and turbulent simulations (Figure 3). The laminar thermal's circulation does increase slightly over time, which may be due to small baroclinic contributions in Equation (4). Variations in the turbulent case are larger, which can be attributed to buoyant fluid not contained entirely within the core of the vortex ring and small detrainment events which occur throughout the thermal's rise. Simplifying Equation (2) to Equation (3) made relating the impulse to the thermal's dynamic variables analytically tractable. In doing so, we assume the integrated buoyancy is constant and imparts an impulse F , and that the radius of the core region of the vortex is much smaller than the radius of the entire vortex. Visually inspecting Figure 1 casts doubt on the validity of this assumption as the core region is clearly finite. Despite this, we find that the approximations closely follow the actual impulse for both laminar and turbulent simulations (Figure 4), verifying their validity. Deviations in the laminar simulation grow for > 3.5 because the thermal begins to interact with the top boundary.
We now test the predictions for the thermals' characteristics, a( ), ⟨ ′ ⟩( ), and z( ) made in Section 3 and plotted in Figure 5. All plots are normalized by their spin-up value at = 1. We find that Equations (11), (12), and (14) agree quite closely with the simulations. Deviations occur for  (14) are the dashed lines and simulation values are the solid lines for the laminar (yellow) and turbulent (purple) cases. The deviations in the laminar data for > 3.5 is due to the thermal interacting with the top boundary. We find that both thermals exhibit power law behaviour, consistent with Lecoanet and Jeevanjee (2019) the laminar simulation at > 3.5, because the thermal starts to interact with the top boundary.
We verify the prediction Equation (15) of the fractional entrainment rate, (a) in Figure 6. The predictions match the directly measured entrainment closely. Deviations in the laminar case occur for > 3.5 because the thermal begins to interact with the top boundary of the simulation. We find that the turbulent prediction overestimates the actual entrainment, but the reasons for this are unclear, given how well the approximations in the theory hold up.

Mechanism denial of buoyancy-driven entrainment
To verify that buoyancy (not turbulence) is the dominant source of entrainment, we ran both laminar and turbulent F I G U R E 6 Verification of entrainment theory for Re = 630 and 6300. Plotted is the fractional entrainment rate as a function of thermal radius a. Predicted values from Equation (15) are dashed lines and simulation data are solid lines. Our model accurately predicts the entrainment rate in both laminar and turbulent simulations, although deviations are slightly larger in the turbulent case. Note the similarity in between laminar and turbulent simulations. Deviations occur for the laminar simulation at a > 2.7, because the thermal starts to interact with the top boundary simulations in which the gravitational force is set to zero part way through the simulation. We turn off gravity at = 1.5, after the thermal has already spun up into a vortex ring. The sudden removal of gravity in the simulation skews the balance of forces in the fluid, causing the thermal to assume a new trajectory which is no longer self-similar (Figure 7). With entrainment absent, the vortex ring should rise with constant velocity. In the laminar case, the height indeed increases linearly with time. However, in the turbulent case, we find the height still increases as the square-root of time, similar to buoyant thermals. To calculate the thermal volume, we respectively use linear and square-root time dependence to approximate the height (Lecoanet and Jeevanjee, 2019 give more details about the thermal tracker). In both cases, the thermal no longer grows in size appreciably compared to the unmodified simulations, as shown in Figure 7.
Calculating entrainment accordingly, we see that in both the laminar and turbulent cases, e drops significantly (Figure 8). Without buoyancy, the observed entrainment efficiency drops dramatically, giving the most direct evidence of the central role of buoyancy in entrainment. This nicely complements the mechanism affirmation experiments of Lecoanet and Jeevanjee (2019) who kept F I G U R E 7 2D vertical slices at the thermal midpoint of the density perturbation at = 0, = 1.5, and = 4. The left panel shows the unaltered Re = 6300 simulation, where the thermal's evolution is self-similar. The right panel shows the mechanism denial simulation which starts with the same initial conditions, but then is altered by removing gravity when > 1.5. The thermal continues to rise because of its impulse, but it no longer expands appreciably because the buoyancy-driven mechanism of entrainment has been removed. Note that the density perturbation field reduces to a passive scalar when g = 0 because it is no longer present in the governing equations buoyancy and instead removed turbulence to show that the Reynolds number had little effect on the measured entrainment rate.

Second-order effects
Buoyancy is essential and is the leading-order effect in entrainment in dry thermals. However, even after removing gravity from the simulations, both the laminar and turbulent thermals continue to entrain a small amount (Figure 8). Note that there are no baroclinic torques without gravity. Equation (3) hence implies R ∼ 1∕ √ Γ so an increase in R suggests the circulation is decreasing with time. The residual entrainment appears to depend on the level of turbulence of the thermal, as the entrainment in the Re = 6300 simulation is roughly double the entrainment in the Re = 630 simulation. Furthermore, we find the laminar thermal's height increases roughly linearly with time, whereas the turbulent thermal's height follows a square-root time dependence.

F I G U R E 8 Confirmation of buoyancy-driven entrainment,
showing the entrainment efficiency e as a function of time. We directly test the idea of buoyancy-driven entrainment by setting g = 0 at = 1.5, sufficiently after the thermals have spun up. The original simulations are dashed lines, and the simulations with gravity removed midway through are solid lines. We find in both laminar (yellow) and turbulent (purple) thermals the entrainment efficiency sharply drops off to less than one third of their original values, giving the most direct evidence of the central role of buoyancy in entrainment. The residual entrainment that remains can be attributed to viscous effects (Section 5.1). Deviations occur for the laminar simulations at > 3.5, because the thermal starts to interact with the top boundary In the absence of gravity, the circulation evolves according to We hypothesize the entrainment seen in Figure 8 is viscous entrainment. It may then be puzzling that the turbulent thermal (with lower viscosity) entrains more than the laminar thermal (with higher viscosity). To understand this, we will now estimate the size of the integrand in Equation (20) in both laminar and turbulent cases.
In the laminar thermal, the vorticity is largest on the length-scale of the vortex radius R, so we can estimate ∼ w∕R and ∇ ∼ R −1 . Then the integrand scales like However, in the turbulent case, we see the vorticity is largest on small scales (Figure 1). For turbulent flows, the enstrophy is typically maximum near the viscous scale, L v ∼ R Re −3∕4 (Thorne and Blandford, 2017). We then estimate ∼ w∕L v and ∇ ∼ L −1 v , so the integrand scales like This illustrates that the viscous torques may change the circulation more in turbulent thermals than laminar thermals! This is consistent with the substantially larger entrainment in turbulent thermals than laminar thermals when there is no gravity (Figure 8). Hence, with gravity absent, we find there is a small amount of residual entrainment, which is turbulent in origin. These effects are small for buoyant thermals in the presence of gravity.
One subtlety in the argument above is that the integrand in the laminar case is everywhere negative, whereas in the turbulent case, it has both negative and positive contributions, leading to substantial cancelations. We expect that cancelations in the turbulent case will make the net change in circulation due to viscous torques independent of Reynolds number, rather than increasing with Reynolds number. We also note that the slight increase in circulation for the laminar simulation including gravity (Figure 3) is likely due to baroclinic torques, given that the viscous torques should decrease the circulation. For a more complete discussion of viscous torques, see Nikulin (2014), which builds a similar analytical model for buoyant vortex rings and includes the effects of turbulent dissipation on the circulation.

Conclusion
In this paper, we have developed a simple theory for thermals based on the vertical momentum Equation (9) and the kinematic constraint Equation (3). This allowed us to derive a complete theory for dry thermals without any additional, empirical constants. We also demonstrated that dry thermals entrain primarily through a buoyancy-driven processes by setting g = 0 midway through the simulation runs. The entrainment is seen to drop off sharply and the small, residual entrainment that remains can be attributed to viscous dissipation that depends on the Reynolds number of the thermal. This is in contrast to many other explanations of entrainment in convection which rely on a primarily turbulent or stochastic entrainment hypothesis (Morton et al. 1956;Escudier and Maxworthy, 1973;Romps, 2010). Our theory also resulted in the explicit expression (Equation (17)) for the entrainment efficiency, which captures the well-known F∕Γ 2 scaling for e and is consistent with parametrizations of gross entrainment as proportional to B∕w 2 , as in Gregory (2001) and Tan et al. (2018). The verification of the B∕w 2 scaling here and the 1∕a scaling in Lecoanet and Jeevanjee (2019) might also explain why these scalings outperformed others in the moist convection study of Hernandez-Deckers and Sherwood (2018).
However, to truly bridge the gap between our idealized dry thermals and real-world cumulus convection, future work should develop a new suite of simulations for moist thermals, perhaps with a simplified model of condensation as in Vallis et al. (2018). The current theory will likely need modification because the assumption of constant integrated buoyancy breaks down in moist thermals, which can generate buoyancy in their centre due to latent heat release (Romps and Charn, 2015;Morrison and Peters, 2018). As for our Boussinesq approximation, the framework presented in this paper has already been generalized to anelastic and compressible thermals in Anders et al. (2019). They study thermals in the context of solar convection, and the excellent agreement between theory and simulation in their more general paradigm suggests that our framework might also be amenable to moist thermals. A simple model for the moist case would provide a foundation on which a thorough understanding of cumulus convection may be built, bridging the gap between simulation and theory in climate and atmospheric modelling (Held, 2005;Jeevanjee et al. 2017).