Airborne lifetime of respiratory droplets

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


I. INTRODUCTION
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 meters, 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 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 8 . Saliva droplets do not completely vaporize due to non-volatile 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 R 0 10 µm), which initially a) Electronic mail: avshalom.offner@ed.ac.uk travel within a puff of high humidity that slows their vaporization and increases their airborne lifetime 3,4,13 .

II. MODEL
Droplets expelled from the human body span a wide range of sizes, from less than a micron to hundreds of microns in diameter 14,15 . The vast majority of expelled droplets travel at sufficiently low velocity 16 to satisfy Re = |v − u|R/ν a ≪ 1, where Re is the Reynolds number, v and u the droplet and air velocities (bold letters denote vector quantities), respectively, R the droplet radius, and ν a 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 17 where x is the droplet position, m = 4πR 3 ρ/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 18 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 19 (see table I).
Eq. (2) is obtained from the complete Maxey-Riley model 20 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πR 2 is the droplet surface area, p and T are the droplet vapor pressure and temperature, p a and T a their counterparts for air, R 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)-(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, R 0 ≡ R (t = 0) is the droplet initial radius, τ = 2R 2 0 / 9ν a ϑ ρ is the Stokes time scale, and P is the atmospheric pressure. To simplify our model, we assume that Pr = Sc = λ , where Pr = µ a c p,a /k a is the Prandtl number with c p,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 21 , the mass and heat convection dimensionless parameters -the Sherwood and Nusselt numbers, Sh (Re, Sc) = h m R/D and Nu (Re, Pr) = h h R/k a , respectively -are identical for Pr = Sc, and therefore we write Sh = Nu = F . Finally, we write the model in dimensionless form, omitting all hats for convenience and recalling that dm = ρAdR to obtain where S ≡ R 2 , θ ≡ T − 1 is the dimensionless temperature difference between the droplet and air, ϑ M = M w /M a and ϑ c = c p,a /c p,d are the ratios between water and air molecular mass and heat capacity, respectively, H d = H/ (c p T a ) is the droplet latent heat parameter, C (S; Kn) is the dimensionless form of (3) with Kn = l/R 0 the Knudsen number, and Φ is the scaled vapor pressure difference between the droplet and air. The   25 .
latter is expressed through the Clausius-Clapeyron relation, in the form where φ is the air relative humidity, which is assumed to be constant, is the vapor latent heat parameter, and T b is the water scaled boiling temperature. The term η/S 1/2 is the addition of capillary evaporation 23 , stemming from the droplet curvature, where η = 2γ c ϑ ρ ϑ M / (PR 0 ) 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 the increasing concentration of non-volatile substances in its composition 24 . Here we assume that saliva may be modelled as a dilute solution throughout its vaporization, such that the loss of mass from the droplet is expressed only through a change in volumewhich is proportional to S 3/2 -while the droplet density remains that of water, independent of the (small) concentration of non volatiles. The dimensionless constant β depends on the saliva reference composition; its value is determined below using experimental measurements of saliva droplet vaporization. To make the system (7a)-(7d) complete, an explicit form must be assigned to F (Re, λ ). In order to recover the wellknown D 2 law of vaporization 26 , i.e. that S decreases linearly with t, F must be constant. This decouples (7c)-(7d) from the droplet dynamics (7a)-(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 Re → 0. Instead, we fit values to F and β according to experimental measurements of saliva droplet vaporization.
The dots mark experimental measurements from Lieber et al. 11 , solid lines are the numerical solution to (7c)-(7d), dashed and dashed dotted lines are the D 2 linear laws with slopes α and the equilibrium size S eq , 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 · 10 −4 -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, S eq .

A. Droplet vaporization
The vaporization process of saliva droplets may be divided into three distinct stages: initial cooling, D 2 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 2 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 figure 1. Below we describe the process in each of these stages, and derive analytic approximations that are key for the analysis that follows.

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 vaporiza-tion and convection -H d ϑ 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.

Stage 2: D 2 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 D 2 law 26 . The slope, can be estimated by fitting a straight line to experimental measurements of S vs. t in the vaporizing stage. (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 figure 1 (blue), yielding F = 0.707. This value is then substituted back to (7c)-(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.

Stage 3: saturation to equilibrium
As water vaporizes from a droplet and its size diminishes, the concentration of non-volatiles 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 S 3/2 is proportional to the droplet volume. Following the assumption of 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 Φ (S eq , θ 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 H v , T b ∝ T −1 a , yielding S eq ∝ T  (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, S eq , as a function of φ . An increase in φ sets the equilibrium between the saliva and air at a lower concentration of non-volatiles, which increases the droplet equilibrium size. strongly affects the equilibrium size, with S eq increasing with φ as shown in figure 2 (green curve and right vertical axis). As a result, droplets at high relative humidity are larger and therefore fall faster to the ground.
The results in figure 1 demonstrate that a saliva droplet size evolution may be approximated by the piecewise form in which t eq = (1 − S eq ) / |α|, with only a marginal sacrifice in accuracy. This approximation is used to obtain analytic results for the droplet airborne lifetime in the next section.

B. Droplet dynamics
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 z 0 > 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.

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 R 0 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 0 and v 0 the droplet initial height and vertical velocity, respectively. In deriving (12) we considered the case |α| ≪ 1 and S 3/2 eq ≪ 1, which holds throughout the realisable range for T a and φ (see appendix A for the complete derivation). The first and second terms on the RHS of (12), z 0 + v 0 (1 + A Kn), denote the droplet initial conditions, where |v 0 | ≪ z 0 is typically obtained since z ∝ τ −1 in our scaling (6). For simplicity, all the results for quiescent air are calculated with v 0 = 0. The third term represents the altitude decrease by time t eq , inversely proportional to the D 2 vaporization rate, |α|, which decreases nearly linearly with φ (see appendix B). This indicates that an increase in 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 terminal velocity v t = −g S eq + A Kn S eq , 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 t eq = (1 − S eq ) / |α|, one easily derives an analytic approximation for a droplet airborne lifetime in quiescent air with v 0 = 0, Figure 2 (blue curve, left vertical axis) shows the lifetime of two respiratory droplets (R 0 = 20, 30 µm) in quiescent air as a function of relative humidity, calculated both analytically through (13) as well as 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 S 2 eq ≪ 1.

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 in which W is a Wiener process, γ = 0.58 s −1 is fitted using experimental measurements of 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.

III. RESULTS AND DISCUSSION
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 R 0 = 20 − 50 µm; calculation with smaller droplets replicated the qualitative behavior of R 0 = 20 µm, and larger droplets reach the ground within seconds regardless of the relative humidity and air velocity. Figure 3 show trajectories of 5 droplets with R 0 = 20 µm at 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 figure 3a 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 figure 3b -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 figure 3. In the context of disease transmission, however, it is essential to quantify the probability for anomalously long lifetimes which can dramatically increase the rate of transmission. Accordingly, we characterize the entire range of lifetimes statistically by calculating 5,000 droplet trajectories and collating their lifetimes into a probability density function (PDF). Figure 4 show these PDFs for saliva droplets with initial radii R 0 = 20, 35 and 50 µm, at relative humidity of 60% (red), 70% (blue), and 80% (green), for (a) U = 0.1 and 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 figure 4a (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 (figure 4b), the lifetime variability is significantly enhanced. Recalling that lognormal distributions appear as Gaussians on a logarithmic scale, we note that the PDFs for R 0 = 35 µm in figure 4b deviate strongly from 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 figure 4, we separate the discussion on the lifetime of large (R 0 = 50 µm) and small (R 0 = 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 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 towards 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 normal distribution. For intermediate size droplets (R 0 = 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 figure 4b. The average lifetime for these relatively large droplets is 45 seconds, however 11% of all droplets are predicted to remain airborne more than 90 seconds. For comparison, less than 1% of droplets for the equivalent PDF in figure 4a remain airborne after 90 seconds. As such large droplets can potentially carry a significant viral load, the non-negligible probabilities for anomalously long lifetime can have an effect on virus transmission.

IV. CONCLUDING REMARKS
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 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 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.
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 S (t) = 1 − |α|t, t < t eq S eq , t > t eq , subject to the initial conditions z (t = 0) = z 0 and v (t = 0) = v 0 . 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 saliva droplet vaporization 11 , we note that |α| ∼ 10 −4 − 10 −2 and S eq ∼ 0.02 − 0.2 for the realisable range of air temperature and humidity. Further, for R 0 14 µm we have Kn 0.05. Accordingly, we consider the asymptotic limits |α| ≪ 1, S 3/2 eq ≪ 1, and Kn 2 ≪ 1, to finally obtain − g S eq + A Kn S eq (t − t eq ) , t > t eq (A4) where t eq = (1 − S eq ) / |α| and S L (t) = 1 − |α|t.