Airborne lifetime of respiratory droplets Open Access

Various viral diseases, including SARS-CoV-2, spread through respiratory transmission.¹ 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,² 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.³ and Wang et al.⁴). These works along with others^5–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.⁸ 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.¹¹ 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.¹² This prevents our analysis from describing the dynamics of small droplets (with initial radius $R0≲10 μ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 velocity¹⁶ to satisfy $Re=|v−u|R/νa≪1$⁠, 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 model¹⁷ 

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 spheres¹⁸ 

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 Davies¹⁹ (see Table I). Equation (2) is obtained from the complete Maxey–Riley model²⁰ 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 $ϑρ≈10−3≪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

where $A=4πR2$ is the droplet surface area, p and T are the droplet vapor pressure and temperature, p_a and T_a are their counterparts for air, $ℜ$ is the universal gas constant, M_w is water's molecular weight, c_p 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

where the hat denotes dimensionless quantities, U is a characteristic velocity, $R0≡R(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 k_a 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,²¹ 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

where $S≡R2, θ≡T−1$ 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

where $ϕ$ is the air relative humidity, which is assumed to be constant, $Hv=HMw/(ℜ Ta)$ is the vapor latent heat parameter, and T_b is the water scaled boiling temperature. The term $η/S1/2$ is the addition of capillary evaporation,²³ stemming from the droplet curvature, where $η=2γcϑρϑM/(PR0)$ is the scaled surface tension with γ_c the water–air surface tension. The term $βS−3/2$⁠, with β a constant, is the correction to the droplet boiling temperature due to the increasing concentration of nonvolatile substances in its composition.²⁴ 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 D² law of vaporization,²⁶ 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,²² giving $F=2$ for convection over a single sphere in the limit $Re→0$⁠. 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, D² 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 D² 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.

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 $βS−3/2$ and $ηS−1/2$ in (8) since these only become significant for $S≪1$⁠, whereas here $S≈1$⁠.

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 D² law.²⁶ The slope,

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.,¹¹ 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.

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 $βS−3/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 S_eq > 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 S_eq, by solving $Φ(Seq,θeq=0)=0$⁠. The limit $η→0$⁠, which neglects the minor effect of surface curvature, gives the simple expression

By comparing (10) to the droplet equilibrium size measured experimentally, we fit the constant $β=8×10−4$⁠. The dependence of S_eq on temperature is weak as both $Hv,Tb∝Ta−1$⁠, yielding $Seq∝Ta2/3$⁠, which varies very little in the indoor temperature range 285–$305 K$⁠. The relative humidity, on the other hand, strongly affects the equilibrium size, with S_eq 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.

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

in which $teq=(1−Seq)/|α|$⁠, 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.

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 $R0≳5 μ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

with z₀ and v₀ the droplet initial height and vertical velocity, respectively. In deriving (12), we considered the case $|α|≪1$ and $Seq3/2≪1$⁠, which holds throughout the realizable range for T_a 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 t_eq, inversely proportional to the D² 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 S_eq 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=(1−Seq)/|α|$⁠, one easily derives an analytic approximation for a droplet airborne lifetime in quiescent air with $v0=0$⁠,

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/2≪1$⁠.

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

in which W is a Wiener process, $γ=0.58 s−1$ is fitted using experimental measurements of the indoor air velocity,²⁷ 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.5 m$ to reach the ground, by solving (7a)–(7d) and (14) numerically. In what follows, we only present results for droplets in the range $R0=20–50 μ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.3 m/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 $u≠0$⁠. 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.

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.3 m/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$⁠.

All the histograms in Fig. 4(a) (⁠$U=0.1 m/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 $(u−v)/C S(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.²⁸ 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

with

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,¹¹ we note that $|α|∼10−4−10−2$ and $Seq∼0.02−0.2$ for the realizable range of air temperature and humidity. Further, for $R0≳14 μm$⁠, we have $Kn≲0.05$⁠. Accordingly, we consider the asymptotic limits $|α|≪1, Seq3/2≪1$⁠, and $Kn2≪1$⁠, to finally obtain

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

We begin by rewriting the definition of the D² vaporization rate α,

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 $ηS−1/2$ and composition-driven increase in boiling temperature $βS−3/2$—and hence

The temperature $θD2≤0$ 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$⁠,

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