Modeling the role of respiratory droplets in Covid-19 type pandemics.

In this paper, we develop a first principles model that connects respiratory droplet physics with the evolution of a pandemic such as the ongoing Covid-19. The model has two parts. First, we model the growth rate of the infected population based on a reaction mechanism. The advantage of modeling the pandemic using the reaction mechanism is that the rate constants have sound physical interpretation. The infection rate constant is derived using collision rate theory and shown to be a function of the respiratory droplet lifetime. In the second part, we have emulated the respiratory droplets responsible for disease transmission as salt solution droplets and computed their evaporation time, accounting for droplet cooling, heat and mass transfer, and finally, crystallization of the dissolved salt. The model output favourably compares with the experimentally obtained evaporation characteristics of levitated droplets of pure water and salt solution, respectively, ensuring fidelity of the model. The droplet evaporation/desiccation time is, indeed, dependent on ambient temperature and is also a strong function of relative humidity. The multi-scale model thus developed and the firm theoretical underpinning that connects the two scales-macro-scale pandemic dynamics and micro-scale droplet physics-thus could emerge as a powerful tool in elucidating the role of environmental factors on infection spread through respiratory droplets.

INTRODUCTION

It has been well established that the SARS-CoV-2 virus responsible for the Covid-19 pandemic transmits via respiratory droplets that are exhaled during talking, coughing, or sneezing. Each act of expiration corresponds to different droplet sizes and myriad trajectories for the droplets embedded in the corresponding jets. Wells was the first to investigate the role of respiratory droplets in respiratory disease transmission. Expelled respiratory droplets from an average human being contain dissolved salt with a mass fraction of about 0.01 as well as various proteins and pathogens in varying concentrations. In this paper, to model the outbreaks, we extensively use the evaporation and settling dynamics of NaCl–water droplets as a surrogate model of the infectious droplets. Stilianakis and Drossinos included respiratory droplets in their epidemiological models. However, they neglected the droplet evaporation dynamics and assumed that characteristic post-evaporation droplet diameters are half of the pre-evaporation droplet diameters based on Nicas et al. In the context of the present Covid-19 pandemic, while the role of droplet nuclei and corresponding “aerosol transmission” route are not clear, it is widely accepted that respiratory droplets are definitely a dominant vector in transmitting the SARS-CoV-2 virus. This merits a detailed investigation of the evaporation dynamics of respiratory droplets and development of a pandemic model that is explicitly dependent on the respiratory droplet characteristics. As such, the evaporation mechanism of respiratory droplets are laced with complexities stemming from droplet aerodynamics, initial droplet cooling, heat transfer, mass transfer of the solvent and solute, respectively, and finally, crystallization of the solute—a phenomenon known as efflorescence. All these are strongly affected by ambient conditions in which the droplet evaporates. These urgently necessitate a model based on first principles, which connects the detailed evaporation dynamics of respiratory droplets with the pandemic evolution equations. In this paper, a model for the infection rate constant based on collision theory incorporates the evaporation physics of respiratory droplets, ab initio.

The droplet evaporation model thus developed is first validated with new experimental results obtained from droplets observed to evaporate in an acoustic levitator. While very interesting insights can be obtained from sessile droplet evaporation, after an expiratory event, the floating droplet evaporates in the absence of surface contact. Thus, the levitated droplets are similar to the droplets in atmosphere compared to their sessile counterpart. Furthermore, the desiccation dynamics necessitates a contact-less environment for the droplet. Alongside a droplet evaporation model, a chemical kinetics based reaction mechanism model is developed with final rate equations similar to that yielded by the SIR (Susceptible, Infectious, Recovered) model. In general, the resemblance of the equations modeling kinetics to those of population dynamics is well known. However, the rigorous framework (analytical as well as computational) of chemical reaction mechanisms that can at present handle few thousands of species and tens of thousands of elementary reactions seems particularly attractive. This could be utilized toward adding further granularity in the pandemic model, if required large mechanisms can be reduced systematically with mechanism reduction techniques. Furthermore, it can be integrated into advection-diffusion-reaction equations, and their moments could be solved using appropriate moment-closure methods. However, for any reaction mechanism, the key inputs are the parameters for the reaction rate constant. In our case, one rate constant is shown to be a strong function of the droplet lifetime. Therefore, next, the droplet lifetime is evaluated over a wide range of conditions relevant to the ongoing Covid-19 pandemic, and the growth rate exponents (eigenvalues) are presented. The results do not suggest that factors not considered in this paper play a secondary role in determining the outbreak spread. Rather, this paper aims to establish the mathematical connection between the pandemic and the respiratory droplet dynamics using a well defined framework rooted in physical sciences. This paper is arranged as follows: first, we provide details of the experiments used to obtain the evaporation characteristics of the water and salt solution droplets. This is followed by the reaction mechanism model that yields the equations for the growth rate and the infection rate constant of the outbreaks. This infection rate constant provides the connection and motivation for modeling the droplet evaporation time scales. Next, to evaluate the rate constant, detailed modeling of the droplet evaporation is presented. This is followed by results and discussions. Finally, we summarize the approach and findings in Sec. VII.

EXPERIMENTS

The experiments with isolated evaporating droplets were conducted in a contact-less environment of an ultrasonic levitator (tec5) to discount boundary effects, generally present in suspended, pendant, or sessile droplet setups. The experimental setup with the diagnostics is shown in Fig. 1. A droplet was generated and positioned near one of the stable nodes of the levitator by using a micropipette. The levitated droplet was allowed to evaporate in the ambient condition of the laboratory at 30 °C and at about 50% relative humidity (RH). The transient dynamics of evaporation and precipitation of the evaporating droplet was captured with the shadowgraphy technique using a combination of a CCD camera (NR3S1, IDT Vision) fitted with a Navitar zoom lens assembly (6.5× lens and 5× extension tube) and a backlit-illumination by a cold LED light source (SL150, Karl Storz).

FIG. 1.

Experimental setup showing the acoustic levitation of a droplet illuminated by a cold LED source. A diffuser plate is used for uniform imaging of the droplet. A CCD camera fitted with the zoom lens assembly is used for illumination. The schematic is not to scale.

A set of ten images at a burst speed of 30 fps is acquired every 2 s for the entire duration of the droplet lifetime. The spatial resolution of the images was ≈1 μm/pixel. The temporal evolution of the diameter of the evaporating droplet was extracted from the images using the “Analyze Particles” plugin in ImageJ (open source platform for image processing). The final precipitate was carefully collected on carbon tape and observed in the dark-field mode under a reflecting microscope (Olympus BX-51). A range of initial droplet diameters varying from 300 µm to 1000 µm were investigated in experiments.

A REACTION MECHANISM TO MODEL THE PANDEMIC

In this section, we model the infection spread rate using the collision theory of reaction rates, well known in chemical kinetics. The connection between droplets and the outbreak will be established later. In this model, we adopt the following nomenclature: P represents a Covid-19 positive person infecting a healthy person(s) susceptible to infection. The healthy person is denoted by H (who is initially Covid-19 negative), and R represents a person who has recovered from Covid-19 infection and hence assumed to be immune from further infection, while X represents a person who dies due to Covid-19 infection. We consider one-dimensional head on collisions, and the schematic of a collision volume is shown in Fig. 2. Here, one healthy person denoted by H with the effective diameter σ is approached by a Covid-19 positive person P of the same effective diameter with an average relative velocity . σ can be considered as the diameter of the hemispherical volume of air that is drawn by H during each act of inhalation, which comes out to be approximately 0.124 m.

FIG. 2.

A schematic of the collision rate model for the infection to occur. Infected person P ejects a cloud of infectious droplets D denoted by small red dots, and the cloud approaches a healthy person H with a relative velocity to infect them. The figure also shows the collision volume swept by the droplet cloud D and H with their respective effective diameters.

It is widely believed that Covid-19 spreads by respiratory droplets resulting from breathing, coughing, sneezing, or talking. Thus, we assume that a volume in front of P is surrounded by a cloud of infectious droplets exhaled by P. The droplet cloud is denoted by D, and the maximum cloud diameter is given by σ. Clearly, σ should be determined by the smaller of evaporation or settling time of the droplets ejected by P, the horizontal component of the velocity with which the droplets traverse, as well as the dispersion characteristics. In each such cloud, we assume that there are numerous droplets containing the active Covid-19 virus. The velocity of this droplet cloud relative to H is given by . In such a scenario, we assume that in a unit volume, there are n infected persons and n healthy persons. For a collision to be possible, the maximum separation distance between the centers of D (the droplet cloud) and H is given by

The collision volume—the volume of the cylinder within which a collision between the droplet cloud of P and air collection volume of H should lie for the collision to occur in a unit time—is given by . Thus, the number of collisions between H and the droplet cloud D of P, per unit time per unit volume, that will trigger infections, is given bywhere n and n represent the number of P and H, respectively. Now, given that each collision between P (basically, its droplet cloud D) and H results in conversion of the healthy individual to the infected individual, we can write

Now, we can define [P] = n/n and [H] = n/n, whereas n is the total number of people those who are capable of transmitting the infection, as well as accepting the infection per unit volume, in that given volume. This implieswhere

Here, ω is the reaction rate. Furthermore, if we assume that the mortality rate is about 3% for the ongoing Covid-19 pandemic, we can convert the kinetics of infection spread to a complete reaction mechanism given by the following:

It is to be recognized that H does not become P immediately on contact with the droplet cloud. The virus must proliferate for a finite time after contact to render a person infectious. A person who has just come in contact with the virus and does not have the capability to infect others yet is denoted by P*. k1, k2, and k3 are the rate constants of reactions [R1], [R2], and [R3], respectively. All rate constants must have dimensions of [T]−1 (inverse of time). Clearly, k1 > k3 for the rapid outbreak to occur. It is to be recognized that this framework implies that k1, the rate constant of the second order elementary reaction [R1] resulting from collisions between the droplet cloud from an infectious individual and healthy individual, is purely controlled by physical effects. The rate constants k2 and k3 of the other two first order elementary reactions [R2] and [R3] are essentially decay rates emerging from the time by which the respective concentrations reach e−1 levels of the initial concentration for the respective reactions. Thus, k2 and k3 are purely determined by interaction between the virus and the human body. We know that the approximate recovery time from the Covid-19 disease is about 14 days. Thus, we can assume k3 = 1/14 day−1. We also assume the latency period (not incubation period) to be 1 day; hence, k2 = 1 day−1. Given the importance of k1 in determining the outbreak characteristics, we will refer to k1 as the infection rate constant. The major contribution of this work is imparting a rigorous physical interpretation to k1 and calculating it.

Using Eq. (4), we can write the system of ODEs for d[P]/dt and d[P*]/dt as

In this paper, we are interested in modeling the initial phases of the outbreaks where [H] ≫ [P]. Hence, we can safely assume [H] ≈ [H]0, i.e., the concentration of healthy people remains approximately constant during the early phase of the outbreak and is equal to the initial concentration, which is very close to unity at t = 0, i.e., at the onset of the outbreak. The time of the beginning of the outbreak denoted by t = 0 for a particular location can be assumed to be the day when the number of Covid-19 positive persons equaled 10. [P]0 is [P] at t = 0. Then, [P] can be solved as an eigenvalue problem and is given by

C1 and C2 are constants to be determined from the eigenvectors and the initial conditions [P]0 and . λ1,2 are the eigenvalues. These can be termed growth parameters and are given by

By Eq. (5), . As mentioned before, k2 = 1 day−1 and k3 = 1/14 day−1, which yields . If k2 → ∞, i.e., a healthy person becomes infectious immediately on contact with an infectious person, λ1 → k1 − k3.

Clearly, this model does not yet account for the preventive measures such as “social distancing,” “quarantining” after contact tracing, and population wide usage of masks. We will call this “social enforcement.” However, it can be included by accounting for the time variation in [H]. Social enforcement measures reduce the concentration of healthy, susceptible individuals from [H0] to [H] where the concentration of healthy population susceptible to infection after implementing strict social distancing (at time t = t) [H] < [H0]. In the case of social enforcement, [P] will be given by

Here, [P] = [P] at t = t and λ1,, λ2, are the eigenvalues from Eq. (6) with [H] = [H]. k1, the infection rate constant, remains to be completely determined. It is to be recognized that two of the key inputs of k1 are σ and V since by Eq. (5). As already mentioned, σ is the diameter of the hemisphere from which breathable air is inhaled. σ is the diameter of the droplet cloud. The aerodynamics of the respiratory droplets needs to be analyzed to evaluate these quantities.

MODELING AERODYNAMICS OF RESPIRATORY DROPLETS

The droplets ejected during respiratory events, such as sneezing and coughing, co-follow the volume of air exhaled during the event. Studies have confirmed that due to entrainment, the exhaled air volume grows in diameter, while its kinetic energy decays with time. Specifically, Bourouiba et al. showed that initially, for a short duration, the droplets evolve inside a turbulent jet, while in later stages, the jet transitions to a puff. Recognizing that the ejected droplets during the respiratory event is surrounded by this dynamically evolving air volume and that the motion of the droplets will be strongly coupled due to the aerodynamic drag, we first model the surrounding air in two parts using the analytical results of the turbulent jet and puff, respectively. The axial location, axial velocity, and radial spread of a transient turbulent jet and puff can be expressed, respectively, asandwhere subscripts j and pf denote the jet and puff, respectively. R0 and U0 are the radius and axial velocities at a distance x0. K is a characteristic constant for the turbulent jet and is reported to be 0.457. At the inception of the respiratory event (t = 0), the jet is assumed to have a velocity U = 10 m/s and a radius R = 14 mm—the average radius of human mouth. The characteristic constants for a puff are a ≈ 2.25 and m = (xa)/(3R). Since the continuous ejection of air from the mouth lasts only for the duration of a single respiratory event, the jet behavior persists only for this period and beyond which the puff behavior is observed. The average duration of such events is roughly 1 s. Hence, the velocity and the radial spread of the air surrounding the exhaled droplets will beThe horizontal displacement (X) of the exhaled droplet and its instantaneous velocity (U) due to the drag can be solved withwhere R is instantaneous radius of the droplet, ρ and ρ are gas phase and liquid phase densities, μ is gas phase dynamic viscosity, and C is the drag coefficient, which can be taken as 24/Re for the gas phase Reynolds number, Re = (2ρ|U − U|R)/μ < 30. As it will be stated later, Re for the respiratory droplets were found to be mostly less than 0.1.

By solving Eqs. (10)–(13) over the droplet lifetime, τ, the axial distance traveled by the droplets, X, can be evaluated. The average velocity of the droplet cloud relative to P is V = X/τ. The diameter of the droplet cloud ejected by P can be approximated as twice the radial spread of the exhaled air, σ = 2R(t = τ). It is to be recognized that while the above equations are analytically tractable, given the complexities of the associated turbulent jet/puff, a detailed description of the motion of the droplets necessitates time resolved Computational Fluid Dynamics (CFD) simulations in three dimensions. This has been recently reported in Ref. 28, which simulated dispersion of water droplets using a fully coupled Eulerian–Lagrangian technique including the wind effects. In this paper, we worked with salt solution droplets, accounting for salt crystallization, but did not include wind effects to retain analytical tractability. Nevertheless, the results presented in Subsection VI B are qualitatively consistent with the CFD results.

Due to evaporation or settling, the droplet is present only for a short time τ after it has been ejected. Therefore, the steady state k1 can be defined asJust like in collision theory, not all molecules are energetic enough to effect reactions; in our case, the droplet cloud is not always present. The last fraction (τ/t) is the probability that the droplet cloud with the average diameter σ is present. t is the average time period between two vigorous expiratory events. V = (V + V) + V. We can assume V = V. It is thus apparent that τ appears in σ, V, and in the last fraction in Eq. (14), thereby emerging as a critical parameter of the entire pandemic dynamics. Hence, τ merits a detailed physical understanding. Given the composition of the respiratory droplets, modeling τ is highly non-trivial and is taken up in Sec. V.

MODELING RESPIRATORY DROPLET EVAPORATION

It is well documented in the literature that an average human exhales droplets (consisting of water, salt, proteins, and virus/bacteria) in the range of 1 µm–2000 µm. In this section, we offer a detailed exposition of the evaporation dynamics of such droplets as ejected during the course of breathing, talking, sneezing, or coughing.

The small droplets (<2 µm–3 µm) have a very short evaporation timescale. This implies that these droplets evaporate quickly (<1 s) after being ejected. However, the same conclusion does not hold for slightly larger droplets ejected in the form of cloud (>5 µm). These droplets exhibit longer evaporation time, leading to increased chances of transmission of the droplet laden viruses. In particular, when inhaled, these droplets enable quick and effective transport of the virus directly to the lungs airways causing a higher probability of infection. In general, the smaller droplets (<30 µm) have low Stokes number, thereby allowing them to float in ambient air without the propensity to settle down. For larger droplets (>100 µm), the settling timescale is very small (∼0.5 s). In effect, based on the diameter of the exhaled droplets, there are three distinct possibilities:

Small droplets (<5 µm) evaporate within a fraction of second.

Large droplets (>100 µm) settle within a small time frame (<0.5 s), limiting the radius of infection.

Intermediate droplets (∼30 µm) show the highest probability of infection due to a slightly longer evaporation lifetime and low Stokes number.

In this work, we particularly focus our attention to the modeling of droplets over a large range of diameters from 1 µm to 100 µm. Based on the available literature, we assume that the droplets exhaled during breathing are at an initial temperature of 30 °C. The ambient condition, however, vary strongly with geographical and seasonal changes, etc. Hence, in the following, we conduct a parametric study to determine the droplet lifetime across a large variation of temperature and relative humidity conditions. The droplet evaporation physics is complicated by the presence of non-volatile salts (predominantly NaCl) as present in our saliva. We would also look into simultaneous desiccation of the solvent and crystallization of such salts Subsections V A and V B. Once exhaled and encountering ambience, the droplet will evaporate as it undergoes simultaneous heat and mass transfer.

Evaporation

For the modeling purpose, the exhaled droplets are assumed to evaporate in a quiescent environment at a fixed ambient temperature and relative humidity. In reality, during coughing, talking, or sneezing, the droplets are exhaled in a turbulent jet/puff. However, as shown in Eqs. (10) and (11), the puff rapidly decelerates due to entrainment and lack of sustained momentum source, rendering the average V to be less than 1% of the initial velocity. Furthermore, since the Prandtl number, defined as ratio of kinematic viscosity and thermal diffusivity, is approximately unity (Pr = ν/α ≈ 0.71) for air, we can safely assume that the temperature and relative humidity that the droplets in the puff experience are on average very close to that of the ambient. At the initial stages, the puff will indeed be slightly affected by buoyancy, which will influence droplet cooling and evaporation dynamics. Quantifying these effects accurately, merit separate studies, see for e.g., Ref. 32 for buoyant clouds. In a higher dimensional model, these could be incorporated. Nonetheless, the evaporation rate of the droplet is driven by the transport of water vapor from the droplet surface to the ambient far field. Assuming the quasi-steady state condition, the evaporation mass flux can be written asHere, ṁ1 is the rate of change of the droplet water mass due to evaporation, R is the instantaneous droplet radius, ρ is the density of water vapor, D is the binary diffusivity of water vapor in air, and α is the thermal diffusivity of surrounding air. B = (Y1, − Y1,∞)/(1 − Y1,) and B = C(T − T∞)/h are the Spalding mass transfer and heat transfer numbers, respectively. Here, Y1 is the mass fraction of water vapor, while subscripts s and ∞ denote the location at the droplet surface and at the far field, respectively. The numerical subscripts 1, 2, and 3 will denote water, air, and salt, respectively. C and h are the specific heat and specific latent heat of vaporization of the droplet liquid. For the pure water droplet, the vapor at the droplet surface can be assumed to be at the saturated state. However, as indicated earlier, the exhaled droplets during talking, coughing, or sneezing are not necessarily pure water; rather, they contain plethora of dissolved substances. The existence of these dissolved non-volatile substances, henceforth denoted as solute, significantly affects the evaporation of these droplets by suppressing the vapor pressure at the droplet surface. The modified vapor pressure at the droplet surface for binary solution can be expressed by Raoult’s Law, P(T, χ1,) = χ1,P(T), where χ1, is the mole-fraction of the evaporating solvent (here water) at the droplet surface in the liquid phase and χ1, = 1 − χ3,. The far field vapor concentration, on the other hand, is related to the relative humidity of the ambient. Considering the effects of Raoult’s law and relative humidity, the vapor concentrations at the droplet surface and at the far field can be expressed as

and denote the molecular weights of water and air, respectively. For evaporation, the droplet requires latent heat, which is provided by the droplet’s internal energy and surrounding ambient. It has been verified that the thermal gradient in the liquid phase is rather small. Therefore, neglecting the internal thermal gradients, the energy balance is given bywhere T is instantaneous droplet temperature, and are the instantaneous mass and surface area of the droplet, ρ and e are the density and specific internal energy of the binary mixture of salt (if present) and water, and k is the conductivity of gas surrounding the droplet. is the thermal gradient at the droplet surface and can be approximated as (T − T∞)/R, which is identical to convective heat transfer for a sphere with a Nusselt number of 2. As such, including aerodynamic effects, the Nusselt number is given by . The droplet Reynolds number, Re, was observed to be mostly less than 0.1, and as such, the aerodynamic enhancement of the Nusselt number, i.e., the second term in the right-hand side, is ignored.

Crystallization

Evaporative loss of water leads to an increase in the salt concentration in the droplet with time. As shown before, P(T, χ1,) is a function of the salt concentration in the droplet, which thus must be modeled using the species balance equation, as shown in the following equation:Here, Y3 is the dissolved solute (salt) mass fraction in the droplet. ṁ3,, which represents the rate at which solute (salt) mass leaves the solution due to crystallization, is modeled below. Clearly, Eq. (18) shows that as water leaves the droplet, Y3 increases. When Y3 is sufficiently large such that the supersaturation ratio S = Y3/Y3, exceeds unity, crystallization begins. Here, we use Y3, = 0.393 based on the efflorescent concentration of 648 g/l reported for NaCl–water droplets in Ref. 33. The growth rate of the crystal could be modeled using a simplified rate equation from

Here, l is the crystal length. Following Ref. 35, for NaCl, we find the constant C = 1.14 × 104 m/s, the activation energy E = 58 180 J/mol, and constant g = 1. Using this, the rate of change of the crystal mass, which equals ṁ3,, is given byWe note that while crystallization process could involve complex kinetics of solute, solvent, and ions; a well-studied single-step crystallization kinetics has been used here for tractability. It will be shown that this model is able to predict the experimentally studied droplet lifetime reasonably well.

The governing equations [Eqs. (15)–(20)] manifest that several physical mechanisms are coupled during the evaporation process. For T > T∞, the droplet undergoes rapid cooling from its initial value. The droplet temperature, however, should eventually reach a steady state limit (wet bulb). This limit is such that the droplet surface temperature will be lower than the ambient, implying a positive temperature gradient or heat input. The heat subsequently transferred from the ambient to the droplet surface after attaining the wet-bulb limit is used completely for evaporating the drop without any change in sensible enthalpy. For a droplet with pure water, i.e., no dissolved non-volatile content, the mole-fraction of the solvent at the surface remains constant at 100%, and at the limit of steady state, the droplet evaporation can be written in terms of the well-known D2 law,where

However, for a droplet with the binary solution, the evaporation becomes strongly dependent on the solvent (or solute) mole-fraction, which reduces (or increases) with evaporative mass loss. The transient analysis, thus, becomes critically important in determining the evolution of the droplet surface temperature and instantaneous droplet size. During evaporation, the mole-fraction of the solute increases and attains a critical super-saturation limit, which triggers precipitation. The precipitation and accompanied crystallization dynamics, essentially, reduce the solute mass dissolved the in liquid phase, leading to a momentary decrease in its mole-fraction. This, in turn, increases the evaporation rate as mandated by Raoult’s law, which subsequently increases the solute concentration. These competing mechanisms control evaporation at the latter stages of the droplet lifetime. At a certain point, due to continuous evaporation, the liquid mass completely depletes and evaporation stops. The droplet after complete desiccation consists only of salt crystals, probably encapsulating the viruses and rendering them inactive. If the SARS-CoV-2 virus would remain active within the salt crystal, also known as droplet nucleus or aerosol, Covid-19 could spread by aerosol transmission in addition to that by droplets. In this paper, we focus on infection spread exclusive by respiratory droplets since the role of aerosols is not clear for transmission of Covid-19.

RESULTS AND DISCUSSIONS

Experimental validation

To validate the model, few targeted experiments were conducted to observe isolated levitated droplets evaporating in a fixed ambient condition. Particularly, the droplets with (1% w/w) NaCl solution vaporized to shrink to 30% of its initial diameter during the first stage of evaporation, as shown in Fig. 3. Hereafter, a plateau-like stage is approached due to increased solute accumulation near the droplet’s surface, which inhibits the diameter from shrinking rapidly. However, as shown in Fig. 3, shrinkage does occur (until D/D ≈ 0.2) as the droplet undergoes a sol–gel transformation. The final shape of the precipitate is better observed from the micrographs presented in Fig. 3.

FIG. 3.

Instantaneous droplet images taken by a CCD camera (top left panel) and dark field micrograph of the final salt precipitate (top right panel). Comparison of experiments and simulations in the bottom left and right panels. Evolution of the normalized droplet diameter as a function of time for pure water (left panel) and the salt water solution droplet with 1% NaCl (right panel).

Figure 3 shows the final precipitate morphology for the desiccated droplets. The precipitates display a cuboid shaped crystalline formation, which is consistent with the structure of the NaCl crystal. The size and crystallite structure does show some variation, which could be linked with the initial size of the droplet. Only precipitates from larger droplets could be collected since smaller sized precipitates tend to de-stabilize and fly-off the levitator post-desiccation. While the precipitate from larger sized droplets tend to yield larger and less number of crystals, smaller droplets seem to degenerate into even smaller crystallites. However, this work does not investigate the dynamics of morphological changes of crystallization in levitated droplets.

Figure 3 also displays the comparison between results obtained from experiments and modeling. Experiments were performed with both pure water droplets as well as with droplets of 1% salt solutions. Experiments have been described in Sec. II. For the pure water cases shown in the left panel, simulation results follow the experiments rather closely. In the pure water case, classical D2 law behavior could be observed. For the salt water droplets, a deviation from the D2 law behavior occurs, and droplet evaporation is slowed. This is governed by Raoult’s law where the reduced vapor pressure P(T, χ1,) on the droplet surface results from the increasing salt concentration with time. The evaporation rate approaches zero at about D/D for 0.3 (experiments) and 0.25 (simulations), respectively. However, the salt concentration attained at this stage exceeds the supersaturation S ≥ 1 required for onset of crystallization. Thus, the salt crystallizes, reducing its concentration and increasing P(T, χ1,) such that evaporation and water mass loss can proceed until nearly all the water has evaporated and only a piece of solid crystal as shown in Fig. 3 is left. It can be observed from Fig. 3 that in all cases, the final evaporation time is predicted within 15% of the experimental values. This suggests that despite the model being devoid of complexities (associated with inhomogeneities of temperature and solute mass fraction within the droplet and simple one step reaction to model the crystallization kinetics), it demonstrates reasonably good predictive capability. It is prudent to mention again that although we have done the analysis for the single isolated droplet, in reality, coughing or sneezing involves a whole gamut of droplet sizes in the form of a cloud.

Humans expel respiratory droplets while sneezing, coughing, or talking loudly. Such droplets have a size range from 5 μm to 2000 µm, while the dispersion could depend on the severity of the action. For example, while talking, an average human being will expel ∼600 droplets in a size range of 25 µm–50 µm, but this number goes upto ∼800 in the case of sneezing. For any given act (sneezing, coughing, or talking), the highest number of droplets fall in the range between 25 µm and 50 µm, while the smaller or larger droplets (than the above) are comparatively fewer in number. Nevertheless, the expelled volume of air contains a very small fraction of liquid droplets. To illustrate this point, the droplets are assumed to be in a uniform dispersion. The total volume of air expelled by a human being is estimated to be V ∼ 0.0005 m3. The total volume of liquid for a given mean droplet size is simply (NV), where N is the total number of droplets of size D and V is the volume of such a droplet. The total volume occupied by droplets of different sizes in a given act of coughing or sneezing is . Thus, the volume fraction (ϕ) of liquid droplets during a given act is . Using size distribution, reported in Ref. 36, one can show that ϕ for sneezing, coughing, coughing with covered mouth, and talking loudly is 59 ppm, 549 ppm, 361 ppm, and 263 ppm, respectively. These numbers indicate that respiratory spray is rather sparse, implying that collective evaporation of droplet clusters may not be significant. This also justifies the modeling based on an isolated, contact free droplet.

Ambience specific droplet lifetime

Next, we set out to use this model to predict the droplet lifetime characteristics over a wide range of ambient conditions. The droplet evaporation time t(D) is calculated from the analysis presented in Secs. V A and V B. Indeed, the droplet evaporation competes with gravitational settling. The settling time t(D) is calculated by accounting for the decreasing diameter using the equation for the Stokes settling velocity,The settling time is estimated as that time by which the droplet gets out of the radius from which breathable air is collected—already defined as σ/2 in Sec. III. Mathematically, t is obtained by the following equation:

Clearly, for any condition, while t monotonically increases with D, t monotonically decreases with D. In view of this, it is necessary to estimate the maximum time an exhaled droplet can remain within the collection volume without being evaporated or settled. Such a time can be estimated by defining a characteristic droplet lifetime τ, whereτ is essentially the time where the two curves t, t as a function of D intersect and represents the maximum time a liquid droplet of any size can exist before it is removed either by evaporation or gravity. For a given ambience specified by the ordered pair (T∞, RH∞), D corresponding to τ can be defined as D. Droplets with D > D settle due to gravity, while D ≤ D evaporate, earlier than the droplets with D = D. While droplets with D ≠ D can certainly transmit the disease, those with D establishes the boundaries in terms of the lifetime, cloud diameter, and maximum distance traversed. D is dependent on ambient conditions, i.e., temperature and relative humidity. The distribution of D over a wide range of relevant ambient conditions is shown in Fig. 4(a). Interestingly, at high T∞ and low RH∞, where the evaporation rate is very fast, even a large droplet rapidly shrinks before it can settle. By the same argument at low T∞ and high RH∞, a relatively smaller droplet cannot evaporate quickly; therefore, t = t is attained for smaller droplet sizes. This explains why we observe large D at high T∞, low RH∞, and small D at low T∞ and high RH∞. The evaporation time of the droplet of diameter D, which has been established as the characteristic lifetime τ of the droplet set, is shown in Fig. 4(b). Despite the different initial sizes, we find that τ is minimum at high T∞ and low RH∞ conditions, whereas it is maximum at low T∞ and high RH∞. At the same time, longer lifetime, i.e., large τ, allows the droplet cloud to travel a longer distance axially (X) and disperse radially (σ). Thus, X and σ are observed to be smaller (larger) for high (low) T∞ and low (high) RH∞ conditions, as shown in Figs. 4(c) and 4(d). This also shows how the minimum required “social distance” represented by X is not constant but depends on ambient conditions. Combining these results, it can be concluded that the size, lifetime, distance traveled, and the radial dispersion of the longest surviving droplet is not constant and is a strong function of ambient conditions. In particular, low temperature and high RH enhance the droplet lifetime significantly. Relative humidity strongly affects the droplet lifetime compared to temperature. An increase in droplet lifetime also implies that such droplets stay in the ambient for longer periods and hence travel longer distances as reflected by X. This implies that such droplets can lead to higher infection propensities. The critical size, droplet lifetime, distance traveled, or size of the droplet cloud for any practical condition are readily obtainable from Figs. 4(a)–4(d), respectively.

FIG. 4.

(a) D, (b) τ, (c) σ, and (d) X as a function of T∞ and RH∞.The black dots A and B denote two typical conditions. Case A represents (T∞, RH∞) = (8, 55), while case B represents (T∞, RH∞) = (28, 77). Color bars have been clipped at reasonable values to show the respective variations over the wider region of interest.

In Fig. 5, we look into the evolution of the normalized mass and temperature of the droplets for two cases named as case A and case B, also identified with black dots in Fig. 4. For case A, (T∞, RH∞) = (8, 55), while for case B, (T∞, RH∞) = (28, 77). The unit of T∞ is °C and that of RH is %. These conditions have been specifically chosen to loosely represent the spring weather in North America and South East Asia, respectively, for the early days of the Covid-19 pandemic at both these locations. Figure 5 clearly explains why τ > τ. Indeed, higher RH at case B implies that the temporary hiatus in evaporation due to reduced vapor pressure is reached at a higher water mass load of the droplet than at the case A condition. In both cases, this occurs at about 3 s. However, the higher temperature in B results in faster crystallization kinetics due to the Arrhenius nature of the equation given by Eq. (19), which causes an eventual faster crystallization rate than at A. Figure 5 clearly shows that although the knee in the solvent depletion profiles are attained at the same time, it is the crystallization and simultaneous desiccation dynamics that governs the eventual difference in lifetime of the droplets at two conditions. It remains to be seen whether this result holds for a detailed crystallization reaction mechanism. The temperature evolution plots shown in the right panel of Fig. 5 also reveals how the droplet initially exhaled at 30 °C rapidly cools to the corresponding wet-bulb temperature to subsequently allow heat transfer into the droplet leading to evaporation. However, as the salt concentration reduces due to crystallization, the temperature rises subsequently above the corresponding wet-bulb limits.

FIG. 5.

Evolution of the normalized mass of water in the droplet (left panel) and droplet temperature (right panel) as a function of time for cases A and B.

Calculated growth parameters and growth rates

With the droplet lifetime available over a wide range of conditions, the corresponding infection rate constant and eigenvalues given by Eqs. (5) and (6) could be evaluated. Just to recapitulate, τ determines the infection rate constant k1 by Eq. (14). In turn, k1 affects the exponents of the time dependent infection equation [Eq. (7)], the growth parameters—eigenvalues λ1, λ2 through Eq. (8). The contours of λ1, λ2, and k1 as a function of T∞ and RH∞ are shown in Figs. 6(a)–6(c), respectively. The direct correspondence between τ and k1 is immediately apparent upon comparing the respective Figs. 4(b) and 6(c). The infection rate constant is highest at low T∞ and high RH∞ where the droplet evaporation is slowed due to the slow mass loss rate and enhanced crystallization time. On the other hand, faster droplet evaporation leads to small infection rate constant values at high T∞ and low RH∞. The temperature and relative humidity dependency of the eigenvalues λ1 and λ2 are shown in Figs. 6(a) and 6(b), respectively. The direct correspondence of Figs. 6(a) and 6(b) with Fig. 6(c) and Fig. 4 are established through given by Eq. (8). It should, however, be noted that due to the inherent negative sign of λ2, its influence on determining the growth rate of the infected population is rather limited. It is λ1 that primarily drives the growth of the infected population as apparent from Eq. (7). From Fig. 6(a), we observe that for a fixed T∞, λ1 increases with RH∞, while for a fixed RH∞, λ1 decreases with T∞. Furthermore, we observe that the iso-λ1 contour lines bend and converge at RH∞ > 75%. This means that for RH∞ > 75%, large λ1 > 0.4 is expected over a wider range of temperatures between 5 °C < T∞ < 20 °C. This is a manifestation of the greatly reduced evaporation potential—the difference between water vapor concentration on the droplet surface and in the ambient at high RH∞ conditions. This is further reflected in the dramatic difference in the rate ratio—the ratio of the cumulative number of positive cases on a particular day to the cumulative number of positive cases seven days before: N/N in Fig. 6(d). This figure has been arrived at by assuming the local population density to be 10 000 km−2. Furthermore, we calculate t in Eq. (14) as t = 3600 × 24/Nexp. Nexp is the number of infecting expiratory events per person per day and is assumed to be 3 based on the coughing frequency of 0–16 in normal subjects. We find that purely based on ambient conditions, implying all other factors have been held constant, the rate ratio can be different by an order of magnitude between (T∞, RH∞) = (5, 75) vs (35, 20). Practically, such a contrast might be less apparent in real data in which other important factors such as population density, social enforcement, travel patterns, and susceptible supply exert significant influence.

FIG. 6.

Contours of calculated eigenvalues (a) λ1 and (b) λ2 as a function of T∞ and RH∞. (c) Infection rate constant, k1, and (d) rate ratio over seven days as a function of T∞ and RH∞. Case A represents (T∞, RH∞) = (8, 55), while case B represents (T∞, RH∞) = (28, 77). The rate ratio for A and B are 16.60 and 10.33, respectively. Color bars have been clipped at reasonable values to show the respective variations over the wider region of interest.

SUMMARY

Respiratory flow ejected by human beings consists of a polydisperse collection of droplets. In this paper, we have presented a model for the early phases of a Covid-19 like pandemic based on the aerodynamics and evaporation characteristics of respiratory droplets. The model and its inter-dependencies on the different physical principles/sub-models are summarized in Fig. 7. To our knowledge, this is the first model that utilizes the structure of a chemical reaction mechanism to connect the pandemic evolution equations with respiratory droplet lifetime by first principles modeling of the reaction rate constant. However, it must be recognized that the model assumes conditions where transmission occurs solely due to inhalation of infected respiratory droplets alongside many other simplifying assumptions. The evolution of the droplets is characterized by a complex interaction of aerodynamics, evaporation thermodynamics, and crystallization kinetics. As such, after being ejected, smaller droplets attain the wet-bulb temperature corresponding to the local ambience and begin to evaporate. However, due to the presence of dissolved salt, the evaporation stops when the size of the droplet reaches about 20%–30% of the initial diameter, but now, the droplet salt concentration has increased to levels that trigger onset of crystallization. Of course, these processes compete with settling—the process by which larger droplets fall away before they can evaporate. The smaller of the two, complete evaporation time and settling time, thus dictates the droplet lifetime τ. The infection rate constant derived using collision theory of reaction rates is shown to be a function of the respiratory droplet lifetime (τ), where τ is sensitive to ambient conditions. While the infection rate constant in reality is dependent on numerous parameters, the present approach allows us to compute its exclusive dependence on ambient conditions through respiratory droplet modeling. We find that the respiratory droplets exclusively contribute to the infection growth parameters and infection growth rate, which decrease with ambient temperature and increase with relative humidity. As such, the model could be used for providing fundamental insights into the role of respiratory droplets in Covid-19 type viral disease spread. Furthermore, the model could be used, with extreme caution and in cognizance of its limitations, toward estimating the risk potential of infection spread by droplet transmission for specific ambient conditions of interest from purely physics based calculations.

FIG. 7.

Flow-diagram outlining the interconnections of the model developed.

DATA AVAILABILITY

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