We formulate a model for the dynamics of respiratory droplets and use it to study their airborne lifetime in turbulent air representative of indoor settings. This lifetime is a common metric to assess the risk of respiratory transmission of infectious diseases, with a longer lifetime correlating with higher risk. We consider a simple momentum balance to calculate the droplets' spread, accounting for their size evolution as they undergo vaporization via mass and energy balances. The model shows how an increase in the relative humidity leads to higher droplet settling velocity, which shortens the lifetime of droplets and can, therefore, reduce the risk of transmission. Emulating indoor air turbulence using a stochastic process, we numerically calculate probability distributions for the lifetime of droplets, showing how an increase in the air turbulent velocity significantly enhances the range of lifetimes. The distributions reveal non-negligible probabilities for very long lifetimes, which potentially increase the risk of transmission.

Various viral diseases, including SARS-CoV-2, spread through respiratory transmission.1 A leading precautionary measure to mitigate respiratory transmission is maintaining a “safe distance,” reflecting a priori knowledge on the spread of droplets and their time airborne. This distance, typically 2 m, was calculated based on the model by Wells,2 who employed Stokes' solution to set the droplet velocity as proportional to its surface area, which decreases at a constant rate due to vaporization. This model, however, is too simplistic to represent the dynamics of respiratory droplets, as was demonstrated by various recent works (e.g., Chong et al.3 and Wang et al.4). These works along with others5–7 conducted direct numerical simulations on the dynamics of droplets in expiratory events, showing how small droplets can travel long distances and remain airborne longer than was considered previously. In the present work, we formulate a particularly simple model to study the airborne lifetime of respiratory droplets, which allows us to simulate a large number of droplets, and, thus, quantify the probability of anomalously long droplet airborne lifetime, which can have a disproportionate effect on virus transmission.

Droplets produced by expiratory events vaporize according to the vapor–pressure balance with the air that surrounds them, possibly after a short period of growth.8 Saliva droplets do not completely vaporize due to nonvolatile substances, resulting in droplets saturating to a small finite size as they reach equilibrium with the surrounding air.9,10 Here, we formulate a model for the vaporization of saliva droplets to capture this unique size evolution and calibrate it with experimental measurements.11 This model is then coupled with a simple momentum balance to describe the dynamics of a single, representative droplet. To retain our model simplicity, we consider a constant relative humidity in the air, thus not accounting for the humid turbulent puff exhaled with the droplets.12 This prevents our analysis from describing the dynamics of small droplets (with initial radius R010μm), which initially travel within a puff of high humidity that slows their vaporization and increases their airborne lifetime.3,4,13

Droplets expelled from the human body span a wide range of sizes, from less than a micrometer to hundreds of micrometers in diameter.14,15 The vast majority of expelled droplets travel at sufficiently low velocity16 to satisfy Re=|vu|R/νa1, where Re is the Reynolds number, v and u are the droplet and air velocities (bold letters denote vector quantities), respectively, R is the droplet radius, and νa is the air kinematic viscosity. Taking advantage of the small Reynolds number, we describe the dynamics of airborne droplets using the aerosol limit of the Maxey–Riley model17 

dxdt=v,
(1)
mdvdt=6πRμaC(vu)mg,
(2)

where x is the droplet position, m=4πR3ρ/3 is its mass, with ρ the droplet density, μa is the air dynamic viscosity, g is the gravitational acceleration, and C is the slip correction for drag force over small spheres18 

C=1+lR[A+Bexp(CRl)],
(3)

where l is the air mean free path and A,B, and C are constants determined by fitting (3) to experimental results. In this work, we use the values calculated by Davies19 (see Table I). Equation (2) is obtained from the complete Maxey–Riley model20 by neglecting all terms proportional to the air-to-liquid density ratio ϑρ=ρa/ρ1. Thus, we treat droplets as liquid particles, small enough to be considered rigid (and, therefore, spherical) and heavy enough that ϑρ1031.

TABLE I.

A list of all dimensionless parameters in our model, along with representative values calculated for R0=30μm and T=Ta=20°C. Reference values for properties of air and water were taken from Poling et al.25 

NotationInterpretationValue @ 20°C
Pr Air Prandtl number 0.7 
Sc Air–water Schmidt number 0.62 
λ Approximation for Pr and Sc 0.66 
g Scaled gravity 0.35 
ϑρ Air–water density ratio 1.2×103 
ϑc Air–water heat capacity ratio 0.24 
ϑM Water–air molecular weight ratio 0.62 
Hd H/(cpTa) 1.8 
Hv HMw/(Ta) 16 
Tb Scaled boiling temperature 1.27 
η Scaled surface tension 0.009 
β Saliva compound constant 8×104 
F Convective mass/heat 0.707 
 transfer coefficient  
Kn Droplet Knudsen number 0.002 
A,B,C Eq. (3) coefficients19  1.257, 0.4, 1.1 
NotationInterpretationValue @ 20°C
Pr Air Prandtl number 0.7 
Sc Air–water Schmidt number 0.62 
λ Approximation for Pr and Sc 0.66 
g Scaled gravity 0.35 
ϑρ Air–water density ratio 1.2×103 
ϑc Air–water heat capacity ratio 0.24 
ϑM Water–air molecular weight ratio 0.62 
Hd H/(cpTa) 1.8 
Hv HMw/(Ta) 16 
Tb Scaled boiling temperature 1.27 
η Scaled surface tension 0.009 
β Saliva compound constant 8×104 
F Convective mass/heat 0.707 
 transfer coefficient  
Kn Droplet Knudsen number 0.002 
A,B,C Eq. (3) coefficients19  1.257, 0.4, 1.1 

Volatile droplets lose mass while vaporizing, and hence, their radius varies. To account for this variation, we include a mass balance equation which depends on the droplet temperature. This temperature is, in turn, governed by energy balance, leading to the system

dmdt=h¯mAMw(pTpaTa),
(4)
mcpdTdt=Hdmdth¯hA(TTa),
(5)

where A=4πR2 is the droplet surface area, p and T are the droplet vapor pressure and temperature, pa and Ta are their counterparts for air, is the universal gas constant, Mw is water's molecular weight, cp is the droplet isobaric heat capacity, H is the heat of vaporization, and h¯m and h¯h are the convective mass and heat transfer coefficients, respectively.21,22 In (4) and (5), we assume that convection dominates over diffusion between the droplet and air. Next, we non-dimensionalize variables by letting

t=τt̂,v=Uv̂,u=Uû,x=τUx̂,R=R0R̂,g=Uτĝ,T=TaT̂,andp=Pp̂,
(6)

where the hat denotes dimensionless quantities, U is a characteristic velocity, R0R(t=0) is the droplet initial radius, τ=2R02/(9νaϑρ) is the Stokes timescale, and P is the atmospheric pressure. To simplify our model, we assume that Pr=Sc=λ, where Pr=μacp,a/ka is the Prandtl number with cp,a and ka the air isobaric heat capacity and thermal conductivity and Sc=νa/D is the Schmidt number with D the molecular diffusion coefficient for air–water vapor mixture. This assumption is not strictly necessary for the numerical calculations that follow; however, it greatly assists in drawing physical insight from the model with only a marginal sacrifice in accuracy (see Table I for typical values of Pr and Sc). Building upon the scaling theory for convection over a sphere,21 the mass and heat convection dimensionless parameters—the Sherwood and Nusselt numbers, Sh(Re,Sc)=h¯mR/D and Nu(Re,Pr)=h¯hR/ka, respectively—are identical for Pr = Sc, and therefore, we write Sh=Nu=F. Finally, we write the model in a dimensionless form, omitting all hats for convenience and recalling that dm=ρAdR to obtain

dxdt=v,
(7a)
dvdt=(vuCS+g),
(7b)
dSdt=4ϑMΦ(S,θ)9λF(Re,λ),
(7c)
dθdt=2[HdϑMΦ(S,θ)+ϑcθ]3λSF(Re,λ),
(7d)

where SR2,θT1 is the dimensionless temperature difference between the droplet and air, ϑM=Mw/Ma and ϑc=cp,a/cp,d are the ratios between water and air molecular mass and heat capacity, respectively, Hd=H/(cpTa) is the droplet latent heat parameter, C(S;Kn) is the dimensionless form of (3) with Kn=l/R0 the Knudsen number, and Φ is the scaled vapor pressure difference between the droplet and air. The latter is expressed through the Clausius–Clapeyron relation, in the form

Φ=exp[Hv(1θ+11Tb+βS3/2)+ηS1/2]/(θ+1)ϕexp[Hv(11Tb)],
(8)

where ϕ is the air relative humidity, which is assumed to be constant, Hv=HMw/(Ta) is the vapor latent heat parameter, and Tb is the water scaled boiling temperature. The term η/S1/2 is the addition of capillary evaporation,23 stemming from the droplet curvature, where η=2γcϑρϑM/(PR0) is the scaled surface tension with γc the water–air surface tension. The term βS3/2, with β a constant, is the correction to the droplet boiling temperature due to the increasing concentration of nonvolatile substances in its composition.24 Here, we assume that saliva may be modeled as a dilute solution throughout its vaporization, such that the loss of mass from the droplet is expressed only through a change in volume—which is proportional to S3/2—while the droplet density remains that of water, independent of the (small) concentration of nonvolatiles. The dimensionless constant β depends on the saliva reference composition; its value is determined below using experimental measurements of the saliva droplet vaporization.

To make the system (7a)–(7d) complete, an explicit form must be assigned to F(Re,λ). In order to recover the well-known D2 law of vaporization,26 i.e., that S decreases linearly with t, F must be constant. This decouples (7c) and (7d) from the droplet dynamics (7a) and (7b) and forms a closed system. To determine the value of F, one could invoke Ranz and Marshall's theoretical result,22 giving F=2 for convection over a single sphere in the limit Re0. Instead, we fit values to F and β according to experimental measurements of the saliva droplet vaporization.

The vaporization process of saliva droplets may be divided into three distinct stages: initial cooling, D2 vaporization, and saturation to equilibrium. Figure 1 shows the size evolution of saliva droplets vaporizing in air, where the two last stages are clearly depicted—the linear decrease marking the D2 vaporization stage and the saturation to a constant value corresponding to the last stage. The initial cooling stage typically lasts less than a second and is, therefore, difficult to observe in Fig. 1. Below, we describe the process in each of these stages and derive analytic approximations that are key for the analysis that follows.

FIG. 1.

Time evolution of saliva droplets vaporizing in air, expressed as the dimensionless radius squared S(t)R(t)2. The dots mark experimental measurements from Lieber et al.,11 solid lines are the numerical solution to (7c) and (7d), dashed and dashed dotted lines are the D2 linear laws with slopes α and the equilibrium size Seq, calculated through (9) and (10), respectively. The experimental results for ϕ=0.534 (blue) were used to fit the model constants—F=0.707 and β=8×104—and hence, theory and experiment match by construction. These values were then used to predict the time evolution at ϕ=0.067 (green), showing good quantitative agreement. The results indicate that the droplet size evolution may be approximated as a piecewise function with a linear decay followed by an equilibrium value Seq.

FIG. 1.

Time evolution of saliva droplets vaporizing in air, expressed as the dimensionless radius squared S(t)R(t)2. The dots mark experimental measurements from Lieber et al.,11 solid lines are the numerical solution to (7c) and (7d), dashed and dashed dotted lines are the D2 linear laws with slopes α and the equilibrium size Seq, calculated through (9) and (10), respectively. The experimental results for ϕ=0.534 (blue) were used to fit the model constants—F=0.707 and β=8×104—and hence, theory and experiment match by construction. These values were then used to predict the time evolution at ϕ=0.067 (green), showing good quantitative agreement. The results indicate that the droplet size evolution may be approximated as a piecewise function with a linear decay followed by an equilibrium value Seq.

Close modal

1. Stage 1: Initial cooling

Droplets expelled from a human body are typically warmer than the surrounding air, i.e., θ(t=0)>0. The vaporization process begins with a rapid decrease in θ as both vaporization and convection—HdϑMΦ and ϑcθ in (7d), respectively—cool the droplet, while its size remains nearly unchanged. As θ falls below zero, convection heats the droplet until a balance with vaporization is reached and its temperature stabilizes. The temperature then remains nearly constant throughout the second vaporization stage, and therefore, we denote it by θD2. The value of θD2 may be computed numerically by setting the right-hand side (RHS) of (7d) to zero and neglecting the terms βS3/2 and ηS1/2 in (8) since these only become significant for S1, whereas here S1.

2. Stage 2: D2 vaporization

Once a droplet temperature stabilizes at θD2, the RHS of (7c) becomes nearly constant, giving rise to a linear decrease of S with time—the well-known D2 law.26 The slope,

α=4ϑMFΦ(θD2)9λ,
(9)

can be estimated by fitting a straight line to experimental measurements of S vs t in the vaporizing stage. Eq. (9) then yields an estimate for the value of the constant F. We use the results of Lieber et al.,11 who recorded the vaporization of levitating saliva droplets, to evaluate F by least squares fitting to the data with ϕ=0.534 in Fig. 1 (blue), yielding F=0.707. This value is then substituted back to (7c) and (7d) and used to predict the vaporization at ϕ=0.067 (green), showing a good quantitative agreement. Encouraged by this agreement, we use the model to describe the vaporization throughout the calculations that follow.

3. Stage 3: Saturation to equilibrium

As water vaporizes from a droplet and its size diminishes, the concentration of nonvolatiles increases, leading to an increase in the droplet boiling temperature. This process is accounted for by the term βS3/2 in (8), recalling that S3/2 is proportional to the droplet volume. Following the assumption of the constant relative humidity, the droplet–air vapor pressure gradient decreases with a decrease in S until, at a finite size Seq > 0, the droplet and air reach a state of thermodynamic equilibrium and the vaporization terminates. This equilibrium corresponds to the vanishing of the right-hand side of (7c) and (7d), resulting in θ=Φ=0. This readily provides the means for calculating Seq, by solving Φ(Seq,θeq=0)=0. The limit η0, which neglects the minor effect of surface curvature, gives the simple expression

Seq=[β(Hv+Tblogϕ)Tb2logϕ]2/3.
(10)

By comparing (10) to the droplet equilibrium size measured experimentally, we fit the constant β=8×104. The dependence of Seq on temperature is weak as both Hv,TbTa1, yielding SeqTa2/3, which varies very little in the indoor temperature range 285–305K. The relative humidity, on the other hand, strongly affects the equilibrium size, with Seq increasing with ϕ as shown in Fig. 2 (green curve and right vertical axis). As a result, droplets at a high relative humidity are larger and, therefore, fall faster to the ground.

FIG. 2.

Left: airborne lifetime of droplets released 1.5m above the ground in quiescent air for two initial radii, R0=20 and 30μm, as a function of the relative humidity ϕ. The solid curves are calculated through (13) and the dashed curves are obtained by solving (7a)–(7d) numerically with u =0. The analytic and numeric results are in very good agreement, ratifying that (11) is a good quantitative approximation for S(t). The monotonically decreasing trend of the curves emphasizes that an increase in ϕ decreases the lifetime of respiratory saliva droplets. Right: droplet equilibrium size expressed as the scaled radius squared, Seq, as a function of ϕ. An increase in ϕ sets the equilibrium between the saliva and air at a lower concentration of nonvolatiles, which increases the droplet equilibrium size.

FIG. 2.

Left: airborne lifetime of droplets released 1.5m above the ground in quiescent air for two initial radii, R0=20 and 30μm, as a function of the relative humidity ϕ. The solid curves are calculated through (13) and the dashed curves are obtained by solving (7a)–(7d) numerically with u =0. The analytic and numeric results are in very good agreement, ratifying that (11) is a good quantitative approximation for S(t). The monotonically decreasing trend of the curves emphasizes that an increase in ϕ decreases the lifetime of respiratory saliva droplets. Right: droplet equilibrium size expressed as the scaled radius squared, Seq, as a function of ϕ. An increase in ϕ sets the equilibrium between the saliva and air at a lower concentration of nonvolatiles, which increases the droplet equilibrium size.

Close modal

The results in Fig. 1 demonstrate that a saliva droplet size evolution may be approximated by the piecewise form

S(t){1|α|t,t<teq,Seq,t>teq,
(11)

in which teq=(1Seq)/|α|, with only a marginal sacrifice in accuracy. This approximation is used to obtain analytic results for the droplet airborne lifetime in Sec. II B.

We now turn our focus to the dynamics of respiratory saliva droplets, governed by (7a)–(7b) with S(t) given by the solution to (7c)–(7d). We concentrate our analysis on the vertical motion of droplets as it determines the lifetime—the time for a droplet released from a height z0>0 to reach the ground—which is a key metric for assessing the risk of respiratory transmission. We note that care needs to be exercised in inferring viral transmission risks from droplet airborne lifetimes because aerosolized viruses can be deactivated at varying environmental conditions. Herinafter, z,v,u, and g are scalar quantities representing vertical components, i.e., aligned with gravity. We begin our analysis by considering the simplest case of u =0, corresponding to a room of completely quiescent air.

1. Quiescent air

The case u =0 provides a benchmark result for free-falling droplets, from which valuable physical insight can be drawn. By employing the piecewise form (11) of S, we derive an analytic solution for z(t). The explicit solution to z(t) is tedious and difficult to interpret (see  Appendix A). However, using the smallness of Kn for droplets larger than R05μm (see Table I), we simplify this expression by focusing on the motion after the vaporization terminates and taking a series expansion about Kn =0, yielding

z(t)z0+v0(1+AKn)g6|α|(3+4AKn)g(Seq+AKnSeq)(tteq),t>teq,
(12)

with z0 and v0 the droplet initial height and vertical velocity, respectively. In deriving (12), we considered the case |α|1 and Seq3/21, which holds throughout the realizable range for Ta and ϕ (see  Appendix A for the complete derivation). The first and second terms on the RHS of (12), z0+v0(1+AKn), denote the droplet initial conditions, where |v0|z0 is typically obtained since zτ1 in our scaling (6). For simplicity, all the results for quiescent air are calculated with v0=0. The third term represents the altitude decrease by time teq, inversely proportional to the D2 vaporization rate, |α|, which decreases nearly linearly with ϕ (see  Appendix B). This indicates that an increase in the relative humidity increases the altitude drop during the vaporization stage, which shortens the droplet airborne lifetime. The fourth term recovers the expected linear descent at the terminal velocity vt=g(Seq+AKnSeq), involving Seq that increases with an increase in ϕ. This demonstrates how an increase in ϕ translates to higher droplet settling velocities that lead to shorter droplet airborne lifetimes. All the terms proportional to Kn, which stem from velocity slip on the droplet surface, naturally increase the droplet settling velocity since the drag force is reduced. By setting z(t)=0 in (12) and substituting teq=(1Seq)/|α|, one easily derives an analytic approximation for a droplet airborne lifetime in quiescent air with v0=0,

tfz0(1AKnSeq1/2)gSeq3(12Seq)AKn(43Seq1/2)6|α|Seq.
(13)

Figure 2 (blue curve, left vertical axis) shows the lifetime of two respiratory droplets (R0=20,30μm) in quiescent air as a function of the relative humidity, calculated both analytically through (13) and numerically by solving (7a)–(7d) with u =0 (solid and dashed curves, respectively). The curves clearly show that (i) a droplet airborne lifetime monotonically decreases with an increase in ϕ and that (ii) the analytic approximation (13) closely follows the numerical result, verifying the use of (11) to approximate S(t) as well as the asymptotic approximations |α|1 and Seq3/21.

2. Indoor turbulence

The air motion in indoor settings is turbulent, driven by natural and forced ventilation, people motion, and other factors. To retain the model simplicity, we represent indoor turbulence by an Ornstein–Uhlenbeck process

du=γudt+σdW,
(14)

in which W is a Wiener process, γ=0.58s1 is fitted using experimental measurements of the indoor air velocity,27 and σ=2γU with U the root mean square (RMS) velocity. Any effect of mean air velocity on the droplet motion can be incorporated in future work by accounting for air flow in specific settings and conditions.

We calculate the airborne lifetime of droplets, defined as the time for a droplet released at z=1.5m to reach the ground, by solving (7a)–(7d) and (14) numerically. In what follows, we only present results for droplets in the range R0=2050μm; calculation with smaller droplets replicated the qualitative behavior of R0=20μm, and larger droplets reach the ground within seconds regardless of the relative humidity and air velocity. Figure 3 shows trajectories of five droplets with R0=20μm at the relative humidity ϕ=0.6 (cyan to dark blue colors) and 0.8 (yellow to red). The random trajectories obtained in turbulent air with RMS velocities (a) U =0.1 and (b) 0.3m/s are compared with the deterministic trajectories obtained for quiescent air (black dashed and dashed-dotted curves). The distinction between the ϕ=0.6,0.8 bundles of trajectories in Fig. 3(a) demonstrates that the decrease in droplet airborne lifetime as ϕ is increased, predicted analytically for quiescent air, also holds for u0. An increase in air velocity results in greater variability of droplet airborne lifetime—as clearly seen in Fig. 3(b)—which stems from the enhanced entrainment of droplets with the turbulent air flow.

FIG. 3.

Representative trajectories of the vertical motion of saliva droplets at relative humidities ϕ=0.6 (cyan to dark blue) and ϕ=0.8 (yellow to red), for air root mean square velocity (a) U =0.1 and (b) U=0.3m/s. The black dashed (ϕ=0.6) and dashed-dotted (ϕ=0.8) curves are the deterministic trajectories for a droplet in quiescent air, calculated from (12).

FIG. 3.

Representative trajectories of the vertical motion of saliva droplets at relative humidities ϕ=0.6 (cyan to dark blue) and ϕ=0.8 (yellow to red), for air root mean square velocity (a) U =0.1 and (b) U=0.3m/s. The black dashed (ϕ=0.6) and dashed-dotted (ϕ=0.8) curves are the deterministic trajectories for a droplet in quiescent air, calculated from (12).

Close modal

The findings above can clearly be noted by observing random droplet trajectories, such as the ones depicted in Fig. 3. In the context of disease transmission, however, it is essential to quantify the probability for an anomalously long lifetimes which can dramatically increase the rate of transmission. Accordingly, we characterize the entire range of lifetimes statistically by calculating 5000 droplet trajectories and collating their lifetimes into a probability density function (PDF). Figure 4 shows these PDFs for saliva droplets with initial radii R0=20,35 and 50μm, at the relative humidity of 60% (red), 70% (blue), and 80% (green), for (a) U =0.1 and (b) U=0.3m/s. The ambient indoor temperature is 20°C. We emphasize that the horizontal axis is logarithmic, demonstrating the extensive variability in droplet airborne lifetime. The black solid and dashed curves are distributions fitted to the data (discussed below), and the markers on the horizontal axis give the lifetimes for u =0 (colors matching the histograms), which nearly coincide with the mean values of the PDFs throughout our calculations. This agreement was anticipated since the mean value of u for an Ornstein–Uhlenbeck process decays exponentially, and so the process mean velocity is u¯=0.

FIG. 4.

Probability density function for the airborne lifetime of respiratory droplets released 1.5m above ground with initial radii of 20 (right), 35 (center), and 50μm (left). The red, blue, and green histograms correspond to relative humidity of ϕ=0.6,0.7, and 0.8, respectively. Solid and dashed black curves are lognormal and log-lognormal distributions fitted to the data, respectively. The air root mean square velocity is (a) U =0.1 and (b) U=0.3m/s. Markers on the horizontal, logarithmic axis mark the analytic result for quiescent air (colors matching the histograms). As expected, smaller droplets remain airborne longer. An increase in relative humidity slows the vaporization rate and increases droplets' equilibrium size, thus accelerating their descent to the ground and shortening their lifetime. An increase in air velocity increases the lifetime variability.

FIG. 4.

Probability density function for the airborne lifetime of respiratory droplets released 1.5m above ground with initial radii of 20 (right), 35 (center), and 50μm (left). The red, blue, and green histograms correspond to relative humidity of ϕ=0.6,0.7, and 0.8, respectively. Solid and dashed black curves are lognormal and log-lognormal distributions fitted to the data, respectively. The air root mean square velocity is (a) U =0.1 and (b) U=0.3m/s. Markers on the horizontal, logarithmic axis mark the analytic result for quiescent air (colors matching the histograms). As expected, smaller droplets remain airborne longer. An increase in relative humidity slows the vaporization rate and increases droplets' equilibrium size, thus accelerating their descent to the ground and shortening their lifetime. An increase in air velocity increases the lifetime variability.

Close modal

All the histograms in Fig. 4(a) (U=0.1m/s) are reasonably well fitted by the black solid curves, which represent lognormal distributions. This fit demonstrates the data's heavy “tails,” indicating that the probability for anomalously long lifetime is much larger compared with normal distributions with the same mean and variance. As the air RMS velocity increases [Fig. 4(b)], the lifetime variability is significantly enhanced. Recalling that lognormal distributions appear as Gaussians on a logarithmic scale, we note that the PDFs for R0=35μm in Fig. 4(b) deviate strongly from a lognormal distribution, displaying substantially fatter tails. To showcase this deviation, we fit these PDFs with log-lognormal distributions, denoted by dashed black curves.

To understand the remarkable statistics in Fig. 4, we separate the discussion on the lifetime of large (R0=50μm) and small (R0=20μm) droplets. The large droplets PDFs at varying ϕ are grouped together, indicating that vaporization only weakly affects their dynamics. These droplets reach the ground well before the vaporization terminates, and their motion is approximately ballistic regardless of the relative humidity. The PDFs, in this case, are skewed due to the absorbing boundary condition at z =0, inducing an asymmetry in the effect of air velocity—which fluctuates about a zero mean—on the droplet motion. Indeed, increasing the initial height above the ground allows more time for the air velocity to change direction through the droplet airborne lifetime and entrain it more symmetrically, resulting in convergence toward normal statistics.

Small droplets, on the other hand, are less affected by gravity and do not fall a significant distance, on the average, during their vaporization. The fat tails in their PDFs derive from the nonlinear interplay between vaporization and the drag force acting to entrain droplets to the air flow. At early times, when a droplet vaporizes and S(t) decreases, the drag force (uv)/CS(t) in (7b) increases non-linearly. The symmetry in u about a zero mean, imposed by the Ornstein–Uhlenbeck process, then leads to an asymmetric effect on the lifetime.

Each of the two effects described above—finite domain for large droplets and vaporization for small droplets—leads to a departure from the normal distribution. For intermediate size droplets (R0=35μm), sufficiently large air velocity can trigger a combined effect that dramatically increases the probability for anomalously long lifetimes, as manifested by the ϕ=0.7 PDF in Fig. 4(b). The average lifetime for these relatively large droplets is 45 s; however, 11% of all droplets are predicted to remain airborne more than 90 s. For comparison, less than 1% of droplets for the equivalent PDF in Fig. 4(a) remain airborne after 90 s. As such large droplets can potentially carry a significant viral load, the non-negligible probabilities for an anomalously long lifetime can have an effect on virus transmission.

The results throughout show that increasing the relative humidity shortens the airborne lifetime of droplets, which can potentially lower the risk of respiratory transmission. This is, indeed, the case provided that such an increase does not elongate the pathogen viability, which increases the risk.28 The statistical analysis, in which a stochastic process was used to represent indoor turbulence, proposes that even modest intensity of turbulence significantly increases the probability for a very long droplet airborne lifetime. This can be a cause for concern if droplets are entrained in ventilation-induced vortices, to be balanced by the clear need for ventilation to discharge droplets with a mean flow of air.

This research was supported by EPSRC Programme Grant No. EP/R045046/1: Probing Multiscale Complex Multiphase Flows with Positrons for Engineering and Biomedical Applications (PI: Professor M. Barigou, University of Birmingham).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

We derive an analytic solution for the descent of free-falling droplets by solving the system (7) with the following simplifications: setting the air velocity to u =0 (quiescent air) and approximating the droplet size evolution S(t) using the piecewise form (11). The system then simplifies to

dzdt=v,
(A1)
dvdt=vCS(t)g,
(A2)

with

S(t)={1|α|t,t<teq,Seq,t>teq,
(A3)

subject to the initial conditions z(t=0)=z0 and v(t=0)=v0. The equations are solved successively, with the solution to (A2) substituted into (A1) and integrated, to obtain a closed-form expression. We simplify the expression, which is tedious and is, therefore, not written explicitly, by employing several asymptotic approximations: based on experimental measurements of the saliva droplet vaporization,11 we note that |α|104102 and Seq0.020.2 for the realizable range of air temperature and humidity. Further, for R014μm, we have Kn0.05. Accordingly, we consider the asymptotic limits |α|1,Seq3/21, and Kn21, to finally obtain

z(t)={z0+[v0(1+AKn)+g(1+2AKn)]×(1SL1/2[SL1/2+AKn1+AKn]2/|α|)g6|α|(3[1SL2]+4AKn[1SL3/2]),t<teq,z0+v0(1+AKn)g(3+4AKn)6|α|g(Seq+AKnSeq)(tteq),t>teq,
(A4)

where teq=(1Seq)/|α| and SL(t)=1|α|t.

We begin by rewriting the definition of the D2 vaporization rate α,

α=4ϑMFΦ(θD2)9λ.
(B1)

Recalling that ϑM,F, and λ are constants, α varies only according to Φ, which is generally a function of both θ and S and is given as Eq. (8). However, during the first two vaporization stages, the droplet is large enough so that both effects involving S are negligible—namely, surface tension ηS1/2 and composition-driven increase in boiling temperature βS3/2—and hence

Φ(θD2)=exp[Hv(1θD2+11Tb)]/(θD2+1)ϕexp[Hv(11Tb)].
(B2)

The temperature θD20 reflects a balance between vaporization and convection, acting to decrease and increase the droplet temperature, respectively. In practice, droplets cool to within several degrees below the air temperature, which in scaled form translates to very small variations from zero. Accordingly, we expand (B2) as a series near θD2=0,

Φ(θD2)=(1ϕ)exp[Hv(11Tb)]+O(θD2),
(B3)

showing that Φ(θD2), to leading order, is linearly dependent on the relative humidity. In turn, we find that |α| decreases nearly linearly with ϕ.

1.
L.
Bourouiba
, “
The fluid dynamics of disease transmission
,”
Annu. Rev. Fluid Mech.
53
,
473
508
(
2021
).
2.
W.
Wells
, “
On air-borne infection: Study II. Droplets and droplet nuclei
,”
Am. J. Epidemiol.
20
,
611
618
(
1934
).
3.
K.
Chong
,
C.
Ng
,
N.
Hori
,
R.
Yang
,
R.
Verzicco
, and
D.
Lohse
, “
Extended lifetime of respiratory droplets in a turbulent vapor puff and its implications on airborne disease transmission
,”
Phys. Rev. Lett.
126
,
034502
(
2021
).
4.
J.
Wang
,
M.
Alipour
,
G.
Soligo
,
A.
Roccon
,
M.
De Paoli
,
F.
Picano
, and
A.
Soldati
, “
Short-range exposure to airborne virus transmission and current guidelines
,”
Proc. Natl. Acad. Sci. U. S. A.
118
,
e2105279118
(
2021
).
5.
M.
Abuhegazy
,
K.
Talaat
,
O.
Anderoglu
, and
S.
Poroseva
, “
Numerical investigation of aerosol transport in a classroom with relevance to COVID-19
,”
Phys. Fluids
32
,
103311
(
2020
).
6.
A.
Khosronejad
,
S.
Kang
,
F.
Wermelinger
,
P.
Koumoutsakos
, and
F.
Sotiropoulos
, “
A computational study of expiratory particle transport and vortex dynamics during breathing with and without face masks
,”
Phys. Fluids
33
,
066605
(
2021
).
7.
S.
Olivieri
,
M.
Cavaiola
,
A.
Mazzino
, and
M.
Rosti
, “
Transport and evaporation of virus-containing droplets exhaled by men and women in typical cough events
,”
Meccanica
57
,
567
575
(
2022
).
8.
C. S.
Ng
,
K. L.
Chong
,
R.
Yang
,
M.
Li
,
R.
Verzicco
, and
D.
Lohse
, “
Growth of respiratory droplets in cold and humid air
,”
Phys. Rev. Fluids
6
,
054303
(
2021
).
9.
L.
Liu
,
J.
Wei
,
Y.
Li
, and
A.
Ooi
, “
Evaporation and dispersion of respiratory droplets from coughing
,”
Indoor Air
27
,
179
190
(
2017
).
10.
R.
Dhand
and
J.
Li
, “
Coughs and sneezes: Their role in transmission of respiratory viral infections, including SARS-CoV-2
,”
Am. J. Respir. Crit. Care Med.
202
,
651
659
(
2020
).
11.
C.
Lieber
,
S.
Melekidis
,
R.
Koch
, and
H.
Bauer
, “
Insights into the evaporation characteristics of saliva droplets and aerosols: Levitation experiments and numerical modeling
,”
J. Aerosol Sci.
154
,
105760
(
2021
).
12.
L.
Bourouiba
,
E.
Dehandschoewercker
, and
J.
Bush
, “
Violent expiratory events: On coughing and sneezing
,”
J. Fluid Mech.
745
,
537
563
(
2014
).
13.
L.
Bourouiba
, “
Turbulent gas clouds and respiratory pathogen emissions potential implications for reducing transmission of COVID-19
,”
JAMA
323
,
1837
1838
(
2020
).
14.
R.
Papineni
and
F.
Rosenthal
, “
The size distribution of droplets in the exhaled breath of healthy human subjects
,”
J. Aerosol Med.
10
,
105
116
(
1997
).
15.
S.
Yang
,
G.
Lee
,
C.
Chen
,
C.
Wu
, and
K.
Yu
, “
The size and concentration of droplets generated by coughing in human subjects
,”
J. Aerosol Med.
20
,
484
494
(
2007
).
16.
P.
Bahl
,
C.
Silva
,
C.
Macintyre
,
S.
Bhattacharjee
,
A.
Chughtai
, and
C.
Doolan
, “
Flow dynamics of droplets expelled during sneezing
,”
Phys. Fluids
33
,
111901
(
2021
).
17.
M.
Maxey
, “
The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields
,”
J. Fluid Mech.
174
,
441
465
(
1987
).
18.
M.
Knudsen
and
S.
Weber
, “
Luftwiderstand gegen die langsame Bewegung kleiner Kugeln
,”
Ann. Phys.
341
,
981
994
(
1911
).
19.
C.
Davies
, “
Definitive equations for the fluid resistance of spheres
,”
Proc. Phys. Soc.
57
,
259
270
(
1945
).
20.
M.
Maxey
and
J.
Riley
, “
Equation of motion for a small rigid sphere in a nonuniform flow
,”
Phys. Fluids
26
,
883
889
(
1983
).
21.
T. L.
Bergman
,
A. S.
Lavine
,
F. P.
Incropera
, and
D. P.
DeWitt
,
Fundamental of Heat and Mass Transfer
(
Wiley
,
2010
).
22.
W.
Ranz
and
W.
Marshall
, “
Evaporation from drops
,”
Chem. Eng. Prog.
48
,
141
146
(
1952
).
23.
S.
Thomson
, “
On the equilibrium of vapour at a curved surface of liquid
,”
Philos. Mag.
42
,
448
452
(
1871
).
24.
P.
Atkins
,
Physical Chemistry
(
Oxford University Press
,
1990
).
25.
B.
Poling
,
J.
Prausnitz
, and
J.
O'Connell
,
The Properties of Gases and Liquids
(
McGraw-Hill
,
2001
).
26.
D.
Spalding
, “
Combustion of liquid fuels
,”
Nature
165
,
160
(
1950
).
27.
M.
Loomans
,
The Measurement and Simulation of Indoor Air Flow
(
Technische Universiteit Eindhoven
,
1998
).
28.
E.
Huynh
,
A.
Olinger
,
D.
Woolley
,
R.
Kohli
,
J.
Choczynski
,
J.
Davies
,
K.
Lin
,
L.
Marr
, and
R.
Davis
, “
Evidence for a semisolid phase state of aerosols and droplets relevant to the airborne and surface survival of pathogens
,”
Proc. Natl. Acad. Sci. U. S. A.
119
,
e2109750119
(
2022
).