Two-dimensional mathematical framework for evaporation dynamics of respiratory droplets

In the context of awareness and understanding the spread of the deadly SARS-CoV-2 virus, which is responsible for the COVID-19 pandemic, many research studies^1–4 have shown that airborne droplet is the main route of disease transmission. Other viruses and bacteria-like influenza virus,⁵ Mycobacterium tuberculosis,⁶ and Pseudomonas aeruginosa⁷ also transmit through respiratory aerosol causing the spread of disease. Studies have also proven that aerosols are involved in the spread of SARS, MERS, H1N1, and many other diseases.^8–10 Wells^11,12 was the pioneer to provide one of the first mechanistic description of the disease transmission route through respiratory droplets. The pathogens are transmitted from infected to healthy person through droplet cloud, produced by the respiratory actions during cough, sneeze, speech, etc. This cloud consists of droplets of varying size¹³ where larger droplets settle on ground due to gravity,⁴ while smaller droplets or aerosols can travel larger distances causing elongated exposure. Recent works on COVID-19 have shown that respiratory events such as cough and sneeze are predominantly made of aerosol particles contained inside a multiphase turbulent gas (a puff) cloud that entrains ambient air.^4,14,15

Aerosols are unstable accumulation of liquid or solid particles of different sizes that constantly varies with external velocity field, gravitational settling, temperature, and mass transfer with the surrounding gas medium. The particles in a bio-aerosol are generally 0.3–100 $μm$ in diameter, and its ability to remain in the air and travel larger distances is of primary concern.¹⁶ Moreover, Gorbunov¹⁷ suggested that aerosols can ensemble in regularly used public spaces mainly hospitals, public transports, and supermarkets, facilitating direct human-to-human respiratory bio-aerosol transmission. Enhanced transmission rate in confined spaces is also discussed by several researchers.^18–20 To minimize the transmission, the world scientific and medical community suggests the use of face masks and social distancing of 2 m.^21–23 However, this social distancing norm of 2 m does not hold for aerosol transmission,¹⁷ depending on various atmospheric conditions.

One of the crucial factors in the transport phenomena of these respiratory droplets in the air is the evaporation of the respiratory droplets, which is dependent on the ambient temperature and relative humidity. As presented in review by Pica and Bouvier,²⁴ pathogen infection spread and transmission can highly depend on factors such as relative humidity, ambient temperature, and surrounding airflow. However, the composition of respiratory fluid is highly complex, containing organic compounds, salts, proteins, surfactants, and pathogens (virus or bacteria).²⁵ Therefore, the evaporation dynamics of these complex fluid droplets in the air are complicated. These complex fluid droplets behave significantly different from pure water droplets⁴ as the nonvolatile components present in the fluid affect evaporation dynamics.²⁶ As shown by Vejerano et al.,²⁶ respiratory droplet primarily contains ∼1% (w/w) NaCl dissolved in water (solvent). Moreover, it is also noted that the pathogens (virion) do not affect the evaporation or precipitation dynamics as shown in experimental data¹⁴ and thus neglected in modeling the global transport processes.

Furthermore, during any respiratory event, the droplets are exhaled in the ambient, consisting of a large volume of air, which initially act as a turbulent jet but gradually transforms into a buoyant puff.²⁷ As the droplets within the cloud transmit through jet and puff, the current model incorporates the evolution of these respiratory jet and puff in the mathematical description for realistic approach.

Because of the complexities described above, in this paper, we have numerically simulated a two-dimensional model based on an isolated droplet from the cloud, assuming one-way coupling between the droplet and surrounding air. In this case, the droplet–droplet interaction or effects of surrounding droplets are neglected as the liquid phase volume (droplets) is too small to affect the gas phase (air) for respiratory events. More popular approach with the homogeneous 0D (zero) model⁴ assumes that the droplet is homogenous, and as such, no liquid phase transport is considered. As an important modification to the homogeneous model, we have solved a detailed transport model for liquid phase inside the droplet and investigated a spatiotemporal solute–solvent level heat and mass transfer, which affects the evaporation rate and further nucleation of the droplet. The current model provides a reasonably agreeable results validated with experimental and the homogeneous model data.⁴ Similar to the homogeneous model, here, the temporal or spatial evolution of flow in gas phase is not solved in detail, as it requires a large scale either direct numerical transport (DNS) or large-eddy simulations (LES) to capture the turbulent characteristics.^28–33 Although such perspective could provide a greater insight, they need more computational resources and simulations for varying ambient conditions, which are difficult to achieve. These DNS and LES studies, on the other hand, do not solve the liquid phase transport due to prohibitive computational costs. Thus, findings of our study with detailed liquid phase transport complement the previous approaches. We will show that varying atmospheric temperature and relative humidity based on the geographical location, the two-dimensional model successfully predicts the corresponding evaporation time until nucleation and distance traveled by the aerosol in nominal-windy conditions globally. However, the axial distance increases by multiple folds having a constant evaporation time when average wind velocity is considered.

As mentioned in the introduction, during any respiratory episode such as breathing, sneezing, or coughing, bio-aerosols are ejected in clouds which initially for a short duration evolve as turbulent jet and then transform to a puff as discussed in.⁴ These respiratory droplets can either remain in the air or settle on a surface as fomites. Evaporation of these bio-aerosols dispersed with the cough-cloud is complicated due to the heat and mass transfer of the solvent and nonvolatile salts (predominantly NaCl) present in our saliva (Fig. 1). These conditions strongly influence the evaporation dynamics of a droplet. The two-dimensional model gives an absolute analysis encompassing the gas phase coupled with the liquid phase responsible for complex evaporation phenomenon of the airborne droplet kinetics in the ambient environment. The modeling framework as advocated by Abramzon and Sirignano³⁴ has been adopted here. They introduced a new approximate model of vaporization of a moving fuel droplet extending the classical droplet vaporization model including important effects as variable physical properties and non-unitary Lewis number in the gas phase, influence of the Stefan flow (blowing) on heat and mass transfer, and the effect of the transient liquid heating inside the internally circulating droplet. The schematic of the lifetime of an isolated exhaled droplet in the ambient atmosphere is given in Fig. 1.

The model involves three-part analysis considering the motion of the droplet in the ambient, gas phase around the droplet and liquid phase of the droplet. The motion of the droplet is solved to determine the velocity of the droplet and its trajectory. The droplet trail in the ambient air is a function of ambient temperatures and relative humidity that the droplets experience. Here, we assumed the cough jet velocities to be dominant in the axial direction; therefore, the radial velocity component of the cough jet is neglected.

The mean axial location and velocity together with mean radial spread of the droplet gas phase divided into transient turbulent jet and puff can be expressed,^{4,15,27,35,36} respectively, as

and

where subscripts j and pf mean jet and puff, respectively. R₀ and U₀ are the initial radius and axial velocities at a distance x₀. For a fleeting turbulent jet, K is a characteristic constant based on,³⁷ which is estimated to be 0.457. Since jet cloud initiates from the mouth itself so at t = 0, the radius of the jet event is approximated to be similar to human mouth radius, 14 mm. However, the ejection velocity of the jet cloud U_j,0 depends on the type of respiratory event. Taking examples of various respiratory event, U_j,0 for sneezing is on the order of 40 m/s; for coughing, it is approximately 10 m/s, while for speaking and singing, it is about 2–3 m/s. Thereafter, to scale the turbulent puff, the characteristic constants considered are a ≈ 2.25 and m = (x_p,0a)/(3R_p,0) ≈ 4.³⁵ 

However, the current paper is restricted to only cough events (U_j,0 = 10 m/s), but can easily be adapted for other respiratory events. A significant factor for the jet event must be noted that for a single cough event, the ejection of jet lasts only for approximately 1 s.³⁶ Beyond this time, the turbulent jet undergoes transition to the turbulent puff together constituting the gas phase described as

The momentum equations responsible for the droplet motion and the droplet evaporation rate are defined as

where $Xp$ is the horizontal displacement of the exhaled droplet and $Up$ is its instantaneous velocity. $Vp$ and $Zp$ are the droplet velocity and displacement in the traverse direction, respectively. Here, $Ug$ and $ρg$ are the droplet gas phase velocity and density of the gas phase; $rs$ is instantaneous radius of the droplet; $ρL$ is liquid phase densities, $μg$ is gas phase dynamic viscosity, and $CD$ is the drag co-efficient, which can be taken as 24/$Rep$ for the gas phase Reynolds number, $Rep=2ρvUg−Uprs/μg<30,$ which is mostly less than 0.1 for respiratory droplets. The model assumes that the respiratory jet and the droplets within that are exhaled in equivalent velocity. Therefore, to refrain from the singularity in calculating CD, the initial droplet velocity is considered as 99% of the jet velocity. The variation in humidity and temperature within the droplet cloud is not regarded in the current approach, while it can be found in Ref. 15.

The average temperature is chosen as T_mean = (2 $Ts$ + $T∞$⁠)/3 as suggested by Hubbard et al.,³⁸ where $Ts$ and $T∞$ are the droplet surface temperature and temperature of the gas phase at far-field, respectively. The solution of Eqs. (4) and (5) is determined by the vaporization rate. This vaporization rate is obtained from the coupled liquid and vapor phase analyses, where the solution for vapor phase surrounding the liquid droplet is assumed to be in a quasi-steady-state condition. Thus, the expressions for vaporization mass flux are given as (heat and mass transfer limits)

where $ṁ$ is the rate of change of droplet water mass due to evaporation, $ρv$ is density of water vapor, and $Dυ$ is the binary diffusivity of water vapor in the air, and $αg$ is the thermal diffusivity of surrounding air. Here, $Nu*$ and $Sh*$ are the corrected Nusselt and Sherwood numbers, respectively, according to the relative changes of film thickness due to surface blowing effect. The Spalding heat and mass transfer numbers are defined as $BM$ = $(YW,s−YW,∞) /(1−YW,s)$ and $BT$ = $Cp,l(T∞−Ts)/(hfg+Ql/ṁ)$⁠. Here, $YW$ is mass fraction of water vapor, while subscript s and $∞$ denote location at droplet surface and at far field, respectively. The subscripts W, A, and N will denote water, air, and salt, respectively. $Cp,l$⁠, and $hfg$ are the specific heat and specific latent heat of vaporization of the droplet liquid, and $Ql$ is the heat transferred into the liquid. In the case for respiratory droplets, the evaporation kinetics gets affected as the nonvolatile solutes (1% w/w NaCl) suppress the vapor pressure at the droplet surface. This is taken into account by considering an ideal solution obeying the Raoult's Law,³⁹ $PvapTs, χW,s=χW,sPsat(Ts)$⁠, where $χW,s$ is the mole fraction of evaporating solvent (here water) at droplet surface in the liquid phase⁴⁰ and $χW,s=1−χN,s$⁠. Following the above considerations, the vapor concentrations at droplet surface and at far field can be expressed as

After determining the mass vaporization rate and the heat transferred into the droplet interior $Ql$⁠, the temperature and concentration distributions inside the vaporizing droplet are solved with the help of governing equations for conservation of mass and energy.

Now, if we look inside the liquid phase of the drop, moving droplets produces convective motion due to the intensive liquid circulation caused by the surface friction. Two-dimensional flow is assumed by considering internal flow structures described by the well-known Hill's spherical vortex.⁴¹ 

The flow around the droplet is assumed to be axisymmetric such as the flow characteristics are independent of θ but varies with time. Using the Hill spherical vortex solution, the radial and angular velocity field in the spherical coordinate system $r,θ$ inside the moving evaporating droplet is given as

where $Us=(1/32)Ug−Up(μg/μl)RepCF$ is the velocity of the liquid at the liquid–vapor interface and is calculated by the continuity of the shear stress across the interface; $μl$ is liquid phase dynamic viscosity. Using the correlation given by Renksizbulut and Yuen⁴² that the Stefan flow reduces the friction drag coefficient by a factor (⁠$1+BM$⁠), the skin friction coefficient for an evaporating sphere is given by

The thermal energy equation inside the circulating and evaporating droplet is obtained by non-dimensionalizing the transient convection–diffusion equation in polar coordinates:⁴³ 

Here, $rs¯=rs/r0$ is the dimensionless radius of the droplet; $η=r/rs$ is the dimensionless radial coordinate; $Vr¯=Vr/Us$ and $Vθ¯=Vθ/Us$ are the dimensionless radial and angular velocities, respectively; $T¯=(T−T0)/T0$ is the dimensionless temperature; $τ=αlt/r02$ is the dimensionless time; $β=0.5∂rs¯/∂τ$ is the dimensionless parameter proportional to the surface regression rate of the droplet; and $αl$ is thermal diffusivity of the liquid phase. The initial and boundary conditions are given as

Similarly, the conservation of species equation can be written as⁴³ 

Here, $Y¯N=(YN−YN,0)/YN,0$ is the normalized mass fraction of the solute. The corresponding initial and boundary conditions are

$Kl$ is the thermal conductivity of the liquid phase and $Dv,za$ is the mass diffusivity of NaCl in water. $PeL$ and $LeL$ are the Peclet number calculated based on liquid properties and Lewis number of the liquid phase, respectively. It is to be noted we assumed that temperature and concentration along with droplet surface are uniform; thus, the Marangoni stress terms are not considered in our equation. At the onset of the respiratory event (t = 0), the gas phase velocity of the bio-aerosol cloud is assumed to be 10 m/s and a radius of the cloud as 14 mm, which is the mean human mouth opening.

The numerical simulation has been carried out by solving Eqs. (4)–(16), both for external vapor and internal liquid regions. However, the internal liquid phase solution requires the boundary conditions (14) and (16), obtained by solving the vapor phase. Accordingly, Eq. (4) is solved by numerical integration method using forward marching scheme. The mass flux $ṁ$ calculated from Eq. (6) generates the solution of Eq. (5); hence, the instantaneous droplet radius can be computed by a forward marching scheme. This mass flux is subsequently used as the boundary condition for the conservation of species equation. Following the completion of the mass transfer part of the vapor phase, heat transfer in the vapor phase is formulated implicitly.⁴³ The calculated surface heat flux $QL$ is employed in the boundary condition for the energy equation of the liquid phase, similar to the mass transfer solution. A detailed stepwise algorithm of the numerical procedure ³⁴ is given in the supplementary material with diagram representing the dependence of each equation as Fig. S1.

After obtaining the boundary conditions from the vapor phase, the thermal energy Eq. (13) and species Eq. (15) along with the boundary conditions (14) and (16) in the liquid phase are computed numerically by fully implicit iterative finite difference scheme called standard second-order Peaceman–Rachford ADI method,⁴⁴ which is a double time step algorithm. The resulting spatiotemporal approximations considered in this method are of second order, and the finite-difference grid spacings (⁠$Δη$⁠, $Δθ$⁠) are uniform. Due to the consideration of an implicit scheme, this method is unconditionally stable, where $Δτ$ = 0.0057, $Δη$ = 0.05, and $Δθ$ = 0.157 were taken for the entire simulation process. Moreover, it is investigated that, if we further decrease $Δτ$⁠, no significant change in the numerical results is exhibited.

Due to the jet and consecutive puff profile of the droplet cloud, it affects both droplet velocity and vapor phase velocity. As discussed earlier, the injection of droplet velocity is very close to the injection of respiratory jet, the droplet velocity U_p = 0.99U_g, and subsequently, the difference becomes smaller due to drag. This provides a similar profile of the droplet velocity and vapor phase velocity with time given in Figs. 2(a) and 2(b). The kink in both the figures at t = 1 s arises from the transition of gas phase solution from the jet to the puff behavior as shown in Eq. (3). Both the velocities in the puff phase slowly approach to zero until the final evaporation time before nucleation of the droplets is reached at 28 °C ambient temperature and 48% relative humidity. Figure 2(c) depicts the droplet velocity (⁠$D0$ = 100 $μm$⁠) along with the axial direction $xp$ of the droplet moving in the air. This gives us the idea that a droplet with initial diameter 100 $μm$ can travel upto 3.5 m in non-windy scenario at the above-mentioned ambient temperature and humidity conditions. The trajectory of a single droplet with three different initial diameters moving in the air is generated in Fig. 2(d). The starting point of the droplet is made at a radial location equivalent to the height of the mouth in an average human at 1.65 m. The bigger droplet (⁠$D0$ = 100 $μm$⁠) axially travel upto certain distance due to the jet flow and gradually fall in the ground due to gravity. However, smaller droplets (⁠$D0$ = 50 and 20 $μm$⁠) evaporate in the air and never settles on the ground. These droplets after nucleation may travel longer distances enhancing the spread of disease. Furthermore, under windy conditions, such as in locations of average wind speed of 2 m/s, the axial spread of the droplets (⁠$D0$ = 100 $μm$⁠) increases 20 times than in a non-windy condition, maintaining identical evaporation time. This refutes the standard social distancing norm of 2 m or 6 feet as bio-aerosol can spread up to several meters depending on the size of aerosol and wind speed of the location. Moreover, in this study, the Weber number $≪$ $1$⁠, considering the surface tension of water (very low salt concentration). Thus, breakup or deformation of the droplet is neglected. Apart from the locations close to the mouth, the number density of the droplets in exhaled volume is rather small. Hence, the entire system can be considered to be dilute. As the surrounding jet-puff of the air expands in diameter and the droplets disperse, the probability of collision and subsequent coalescence should be rare.

For validation of the numerical model, few oriented experiments were conducted using an isolated levitated droplet evaporation of 1% w/w NaCl solutions to replicate the actual respiratory droplets in 28 ± 0.2 °C ambient temperature and 48 ± 0.002% relative humidity.⁴ Figure 3(a) exhibits a comparison between experimental outcome, the 0D model,⁴ and the current 2D model for initial droplet diameter as 593 $μm$⁠. The model shows reasonably good agreement with experimental results and 0D model with inconsequential deviations. The current model is terminated immediately with the trigger of precipitation as solute concentration at the surface reaches the critical value $YN,S=0.393.$⁴⁵ Moreover, the validation of the results is also done using varying grid sizes as shown in supplementary material Fig. S2.

In Figs. 3(b) and 3(e), comparisons between the 0D and 2D model are presented for two initial diameters (⁠$D0$ = 50 and 100 $μm$⁠). The ambient T_∞ = 28 °C and RH_∞ = 0.48 were maintained, and initial droplet temperature was T_s,0 = 30 °C. Due to the minimal contrast in the initial phase of injection velocity of the droplets and ambient gas-phase velocity (U_p,0 = 0.99U_j,0), followed by similar profile [Figs. 2(a) and 2(b)], induces a momentary Hills vortex and hence, we expect less inhomogeneity inside the droplet. Thus, the $D/D0$ plot with time for both initial diameters in the 2D model agrees with the 0D homogeneous model. The 2D model predicts a slightly higher wet bulb temperature compared to 0D model. Furthermore, transport of heat and mass from the surface of the droplet to the center indeed takes finite time (an aspect which is ignored in homogeneous model). Nevertheless, the comparison illustrates the fact that the 0D model⁴ demonstrates reasonably good predictive capability. However, under a different circumstance (with large thermal and velocity difference) such as high-speed ambience or droplet injection in plasma jet, the relative velocity is enhanced. Such high relative velocity and temperature differences will lead to significant errors in the 0D homogeneous model in which the current 2D model can handle efficiently.

As expected, for T_s,0> T_∞, the droplet undergoes rapid cooling and eventually reach a steady state limit (wet bulb) as the evaporation progresses. As shown in Figs. 3(d) and 3(e) in both droplet sizes, both models predict a wet bulb temperature around 21 °C (2D model predicts slightly higher temperature), a limit which ensures T_s,0 < T_∞. Consequently, after attaining the wet-bulb limit, we have positive heat transfer from the ambient to the droplet surface, which is used completely for evaporating the droplet. For a respiratory droplet with the dissolved nonvolatile salts, the evaporation is dependent on the solvent (water here) mole fraction at the droplet surface. Due to evaporative mass loss during evaporation, the solvent concentration is reduced at the droplet surface, while the mole fraction of the solute increases. Once the mole fraction of the solute reaches the critical supersaturation limit, precipitation begins.

Tables I and II, in Ref. 46, represent the ambient temperature and relative humidity of different cities for winter and summer/monsoon, respectively. A–F represents outdoor conditions of regions like Berlin, Delhi, Kolkata, London, Nairobi, and New York, respectively.

The evolution of normalized diameter $D/D0$ and surface temperature with time t for different initial diameter $D0$ and seasons are represented in Fig. 4 for the cities given in Tables I and II. It is observed that in cities having higher relative humidity, the droplet takes larger time to evaporate. Due to the larger relative humidity, the vapor pressure reduces and reaches a point of equilibrium between the evaporation of water mass from the droplet surface and condensation of water mass from the vapor layer back to the droplet. This delays the evaporation process for which high humidity cities could not evaporate beyond $DD0∼0.4$⁠. The right-side panel of Fig. 4 shows the evolution of surface temperature of the droplet exhaled initially at 30 °C. The surface temperature rapidly drops to the analogous wet-bulb temperature allowing the heat from the ambient to enter the droplet leading to evaporation. The left and right columns in Fig. 5 show the evolution of temperature and concentration contours, respectively, for an evaporating surrogate respiratory droplet at an ambient condition of T_∞ = 28 °C and RH_∞ = 0.48. The initial droplet diameter and temperature were $D0$ = 100 $μm$ and T_s,0 = 30 °C. The initial solute concentration $YN,0=0.01$⁠. The primary goal of the 2D model is to understand the temperature and solute concentration distribution within the droplet and ascertain the degree of inhomogeneity of the scalars inside the liquid phase. As the solvent vaporizes from the droplet surface, the solute mass fraction near the surface of the droplet is maximum. Messing et al.⁴⁷ proposed that precipitation could occur in various conditions depending on the saturation concentration and temperature. However, in this study, we have considered a homogeneous nucleation phenomenon that assumes a uniform and instantaneous precipitation of salt in the zones that achieved critical supersaturation concentration $YN,S=0.393.$⁴⁴ Figure 3(e) illustrates the droplet temperature plot with time that depicts the initial transient phase up to 100 ms where rapid cooling of a droplet takes place due to evaporation, while heat is supplied from the center of the droplet (⁠$tI−tII$⁠). In this region, the temperature contours follow the flow streamlines line to some extent, thereby exhibiting some convection effects due to relative velocity between droplet and jet/puff. However, as the temperature reaches the wet bulb (⁠$tIII−tVI$⁠), both temperature and concentration contours become concentric exhibiting minimal convection. This is because, at later stages, as the droplet velocity reaches the velocity of the puff cloud, diffusive transport is most dominant. The homogeneous temperature distribution inside the liquid phase is shown in Fig. 5. The temperature distribution contours are plotted for each time step (⁠$tI−tVI$⁠) denoted in Fig. 3(e) to illustrate this effect. However, the concentration distributions exhibit concentric lines at slightly later stages due to low mass diffusivity compared to thermal diffusivity. However, in the case of spray combustions in liquid-fueled engines and in plasma jets, the convective effects become significant. Such phenomena cause distortion of the spherical symmetry of the temperature distribution inside the droplet.^48–55 It is also observed that the salt concentration enhances near the droplet surface due to evaporation. As the droplet moves in the air, its size reduces and the surface NaCl concentration reaches the critical supersaturation limit $%YN,S=39.3%$⁠. Here, % concentration of salt is presented taking the initial concentration as 1%. Salt accumulation slows down the rate of vaporization significantly and crystallization starts from the droplet surface. We observe similar trends for different ambient conditions.

However, the model can capture the inhomogeneity in temperature and concentration gradient in bigger droplets as shown in Fig. 6. For droplet of initial diameter $D0$ = 593 $μm$⁠, the left and right panel shows the temperature and concentration contours, respectively, at ambient temperature of 28 °C and relative humidity 48%. The initial droplet temperature were T_s,0 = 30 °C with initial solute concentration of $%YN,0=1%$⁠. Unlike smaller droplets as shown in Fig. 5, temperature gradient in larger droplets takes longer time of about 6.39 s to attain homogeneity. During this period, the droplet surface temperature reduces due to evaporation and finally attains a wet bulb temperature. This causes rapid cooling of the droplet, while heat is supplied from the center of the droplet. Moreover, both the temperature and concentration contours exhibit stronger convection effects due to relative velocity between the droplet and the ambient (1 m/s). This elucidates the significance of the 2D model to capture such convection effects in larger droplets, which could not be possible in a homogeneous model.

Figure 7 reveals the region in which solute concentration $YN,S>0.393$ in a droplet of initial diameter $D0$ = 100 $μm$⁠. In this region, the solute concentration becomes higher than the supersaturation concentration, which triggers the precipitation formation. The salt precipitate is uniform along the surface of the droplet as homogeneous nucleation is considered. In a nutshell, low relative velocity and hence low convective effects ensure homogeneous temperature and concentration distribution within the droplet in the 2D model. This model is terminated as the solute concentration reaches the supersaturation concentration. However, to upgrade the pseudokinetic model of crystallization,⁴ the current 2D model can be expanded to understand the nucleation phase inside a droplet.

Figure 8 gives the distribution of axial distance $xp$ traveled by the droplet in the air, over a wide range of ambient temperature and relative humidity for two different initial droplet diameter (⁠$D0$ = 20 and 100 $μm$⁠). For both diameters, interestingly we observe that at high ambient temperature and low relative humidity, the distance traveled by the droplet is minimum. This interprets that the time of evaporation $τ$ must be smaller in cases such as this.⁴ However, for low ambient temperature and high relative humidity, the axial distance $xp$ increases as the time $τ$⁠, the droplet takes to evaporate is much larger. The contour plots [Figs. 4(a) and 4(c)] are produced for winter season and [Figs. 4(b) and 4(d)] for summer/monsoon. Different cities across the globe are denoted on the contour plots (Tables I and II) for both $D0$⁠. The axial spread of the droplets just prior to nuclei formation in outdoor conditions is observed to be maximum in Berlin and London during the winters, followed by New York, Delhi, and finally Nairobi and Kolkata. However, during summer/monsoon, Kolkata will have the highest axial spread, followed by London, Berlin, Nairobi, and finally New York and Delhi. Once droplets form nuclei, they can remain suspended in the air for hours and can travel much further distance. In cases of such aerosols, the distance traversed by the droplets prior to crystallization/nuclei formation becomes insignificant and largely redundant. However, for diseases where droplets are the principle mode of transmission rather than aerosols, the contour plots and associated distances as shown in Fig. 8 are significant. In addition, for short-range spread of diseases through droplets, the above data may be significantly useful.

In this paper, we have presented a two-dimensional modeling approach that aids to a detailed scientific-based explanation on the transport and translation process in respiratory bio-aerosol (smaller droplets). Respiratory pathogens such as coronavirus, influenza virus, Mycobacterium tuberculosis, and Pseudomonas aeruginosa exhaled during any respiratory event from an infectious person are carried within these droplets and get transmitted. Therefore, it becomes necessary to understand the evolution of the exhaled droplets characterized by complex interaction of spatiotemporal aerodynamics and evaporation kinetics. Once a droplet cloud is exhaled in any respiratory event, the droplets experience aerodynamical drag, increasing fluid shear rate and causing internal circulations inside the spherical droplet. The main motivation of this approach is to understand detailed flow dynamics in the liquid phase of the droplet while moving in the ambient air. In cases such as turbulent windy conditions causing high relative velocity and greater temperature differences, this model is significantly useful. The transient heat and mass transport equations in liquid phase provide significant knowledge on the internal flow dynamics and concentration gradient. The presence of nonvolatile solutes in the respiratory fluid generates this concentration gradient that enhances at the droplet surface due to evaporation process. The model is simulated until the surface concentration reaches a critical value after which droplet nucleation starts. Moreover, it is observed that the axial spread for bio-aerosols with initial diameter less than 50 μm is higher as they start nucleating in the air and does not settle on the ground. These aerosol particles can float in the air for longer time (hours) enhancing the spread of pathogens. The current model is capable of successfully predicting the evaporation time and axial spread of the droplets in any geographical location on the globe until the point of nuclei formation.

See the supplementary material showing grid independence study and numerical schemes adopted to solve the equations.

The authors thank Dr. Prasenjit Kabi for the experimental data and crystal image. S.B. gratefully acknowledges the funding received from the DRDO Chair Professorship Award. D.C. gratefully acknowledges Tata Innovation fellowship and Astra Chair fellowship. We also acknowledge the esteemed reviewers for their valuable suggestion.

The data that support the findings of this study are available within the article.