Growth of respiratory droplets in cold and humid air

The ambient conditions surrounding liquid droplets determine their growth or shrinkage. However, the precise fate of a liquid droplet expelled from a respiratory puff as dictated by its surroundings and the puff itself has not yet been fully quantified. From the view of airborne disease transmission, such as SARS-CoV-2, knowledge of such dependencies are critical. Here we employ direct numerical simulations (DNS) of a turbulent respiratory vapour puff and account for the mass and temperature exchange with respiratory droplets and aerosols. In particular, we investigate how droplets respond to different ambient temperatures and relative humidity (RH) by tracking their Lagrangian statistics. We reveal and quantify that in cold and humid environments, as there the respiratory puff is supersaturated, expelled droplets can first experience significant growth, and only later followed by shrinkage, in contrast to the monotonic shrinkage of droplets as expected from the classical view by William F. Wells (1934). Indeed, cold and humid environments diminish the ability of air to hold water vapour, thus causing the respiratory vapour puff to super-saturate. Consequently, the super-saturated vapour field drives the growth of droplets that are caught and transported within the humid puff. To analytically predict the likelihood for droplet growth, we propose a model for the axial RH based on the assumption of a quasi-stationary jet. Our model correctly predicts super-saturated RH conditions and is in good quantitative agreement with our DNS. Our results culminate in a temperature-RH map that can be employed as an indicator for droplet growth or shrinkage.

The flow physics of respiratory droplet is crucial to understand the airborne transmission of respiratory diseases (1,2), and recently due to the pandemic of COVID-19, the interest on this subject has been renewed (3)(4)(5)(6)(7)(8)(9)(10)(11)(12). Tiny saliva and mucus droplets, mostly of the size of tens to hundreds of micrometers (13,14), are expelled from the mouth at speaking, screaming, shouting, singing, coughing, sneezing, or even breathing. The expelled micro-sized droplets have long been recognized to be the carriers for viruses which are responsible for the viral transmission from one host to others (15).
Airborne disease transmission is a difficult subject because many intricate factors simultaneously affect the transmission and a complete understanding requires multi-disciplinary collaboration (16,17). In addition to the understanding from epidemiology (18)(19)(20) and virology (21), fluid mechanics also offers essential insight into the mechanisms of airborne disease transmission (2). However, surprisingly little is known on the transport and fate of the respiratory droplets, once they have been expelled from the mouth. For example, how the air temperature and the relative humidity influences the actual fate of the droplets remains an open question. It is a complicated process involving interactions between turbulent eddies, vapour field, and temperature field with the respiratory droplets. Therefore, accurate quantification of these effects is necessary.
On the fate of respiratory droplets in fluid flow, classical work was done by William F. Wells in 1930s (22,23). The picture that Wells developed, at that time in connection with the transmission of tuberculosis, considers only the settling and evaporation of the droplets. After droplets are expelled from the mouth, they settle to the ground by gravity, and, in this simple model, they evaporate with their surface area decreasing linearly with time, the so-called d 2 -law (24). Then for a given ambient temperature and relative humidity, one can estimate the lifetime of the respiratory droplet based on the d 2 -law.
However, in an actual respiratory event, the assumptions of settling and evaporation are insufficient to describe the full physical process. In reality, more complicated flows are possible: Abkarian et al. (25) demonstrated that speech can produce jet-like flow transport, and Bhagat et al. (11) illustrated that body or breathing plumes from a person can also spread droplet nuclei in enclosed spaces. Recently, a new paradigm has been suggested by Bourouiba et al. (1,3,26), which states that instead of the isolated droplet being considered, one should also take into account the importance of turbulent vapour puff. Inspired by this paradigm, using direct

Significance Statement
Influence of environmental conditions on airborne diseases transmission is an important issue, especially during the pandemic of COVID-19. Human-to-human transmission is mediated by the transport of virus-laden respiratory droplets. Here we investigate the problem from a fluid mechanics perspective by conducting numerical simulations to quantify the fate of respiratory droplets in a warm humid coughing puff under different ambient conditions. We reveal a non-intuitive regime with considerable growth of respiratory droplets, dominated by a super-saturated vapour field, preferentially occurring in cold and humid environments. We further propose a theoretical model that accurately predicts the condition for droplet growth. Our work should inform socializing policies and ventilation strategies for controlling indoor ambient conditions to mitigate dispersion of droplets from asymptomatic individuals.  Fig. 1. Flow visualisation snapshots from our direct numerical simulations of water droplets in a warm humid puff in ambient air at (A) θamb = 30°C and (B) θamb = 10°C, both at RHamb = 90%. Corresponding movies can be seen in Movies S1 and S2. The snapshots show (i) vertical 2D planes of the local RH fields, (ii) the instantaneous droplets spatial distribution, (iii) the heavy large droplets which already fell on the ground, and (iv) the instantaneous droplet size histograms versus distance. The local RH planes are taken from the vertical mid-plane of the puff and are plotted on the background for clarity. Droplets are colour coded by their instantaneous sizes. Initial droplet sizes are prescribed with a distribution similar to J. P. Duguid (13) (plotted in Fig. S2) and are injected evenly in time with the same local inflow velocity. The initial temperature of the droplets and puff is 34°C. Both snapshots are taken at 0.6s corresponding to the cut-off time of the puff. In the colder conditions of (B), the expelled humid puff over-saturates (seen as the lighter colored RH field and visible in Movie S2), which in turn dictates growth of smaller droplets caught within the puff. Correspondingly, the droplet counts are confined within a narrower range of sizes in (B) as compared to (A), see Fig. 2A.

D R A F T
numerical simulations, Chong et al. (27) showed that droplets have O(100) longer lifetimes in a turbulent vapour puff at high ambient relative humidity than that predicted by the d 2 -law. Past numerical studies that consider different ambient conditions typically adopt mean flow modelling approaches (28)(29)(30). The role of ambient conditions requires large-scale parameterizations, (31-33) but run the risk of obscuring the underlying physics.
In this study, we aim at quantifying the effect of ambient temperature (θ amb ) and relative humidity (RH amb ) on the respiratory droplets using direct numerical simulations of a turbulent puff. Normal coughing with an actual injection profile has been simulated with 5000 droplets being tracked by the Euler-Lagrangian approach. We reveal a regime in which the respiratory droplets considerably grow before their shrinkage due to evaporation. The droplets growth happens under specific combinations of ambient conditions, preferentially in cold and humid environment. To predict the criteria of having droplet growth, we further propose a theoretical model to estimate the local relative humidity experienced by the droplets, which agrees excellently to our simulation results. Whilst our model can help to formulate public health guidelines or indoor ventilation strategies in order to mitigate the spreading of respiratory droplets, the infectivity of virus-containing droplets and droplet nuclei, which remain after evaporation still requires further investigation, in particular on the virus dose required for infections and whether dried-out droplet nuclei may also be infectious. This issue remains controversial (15) and cannot be answered within a fluid dynamical study.

Results and Discussion
A. Effect of Ambient Temperature on Turbulent Vapour Puff and Droplet Counts. First, we highlight the differences in the local RH values for two different θ amb (10°C and 30°C) at RH amb = 90%. Note that θ amb and RH amb both constitute the control parameters of our study. The local RH field is 2s for θamb = 10, 20, 30°C. As a point of reference, the solid lines denote the initial droplet size distributions. The bin size of the histogram is 2µm. At t = 0.6s (panel A), the peak of the histogram for θamb = 30°C marginally shifts down, whereas the overall histogram for θamb = 10°C remains roughly similar to the initial distribution. The gap is more pronounced at t = 1.2s (panel B): For θamb = 30°C, the peak of the histogram is halved from the initial peak, whereas in contrast, for θamb = 10°C, the histogram still remains roughly unchanged. The relative invariance of the histogram for θamb = 10°C is due to the droplet growth dictated by the super-saturated vapour puff, which is in turn determined by cold and humid ambient conditions, as described in the main text.

D R A F T
shown in figure 1 together with the instantaneous droplet distributions and instantaneous count histogram. At θ amb = 30°C, the local RH field is almost indistinguishable from RH amb , as shown in Fig. 1A. In contrast, at θ amb = 10°C (Fig. 1B), the local RH field reaches super-saturation values larger than 100% (whitish colour shading). In the count histograms in Fig. 1B, most of the droplet diameters remain at least around 20 µm throughout the spread distance, whereas in Fig. 1A, droplet diameters can go below 20 µm. This means that droplets expelled into colder ambient surroundings shrink slower because of the super-saturated local RH field. On the other hand, droplets expelled into hotter ambient surroundings shrink much more rapidly because the local RH rapidly undersaturates.
To illustrate the temperature sensitivity of the droplet size distributions, in Fig. 2, the distributions across three θ amb , 10°C, 20°C, and 30°C are shown at t = 0.6s ( Fig. 2A) and 1.2s (Fig. 2B). In the figure, the initial droplet size distribution (solid line) is also included for comparison. In Fig. 2A, which corresponds to the cut-off time of the simulated puff, many smaller droplets 15um already exist at θ amb = 20°C and 30°C as compared to θ amb = 10°C. Interestingly, the histogram distribution of droplet diameters for θ amb = 10°C is slightly higher and shifted to the right from the initial droplet size histogram, reflecting that a portion of the smaller droplets in fact grew. The reason for this increase is that due to the cold surrounding air, which can take up less vapour than the warm one, the puff gets super-saturated. This effect will be quantified in detail in Section B. At the later time of t = 1.2s (Fig. 2B), most of the droplet sizes for θ amb = 20°C and 30°C are shifted to the left towards much smaller sizes. Remarkably, the size distribution for θ amb = 10°C remains roughly similar to the initially prescribed size distribution, with some lower counts at sizes 15µm.
After having shown the distinct droplet shrinkage behaviour under different ambient conditions, in the following section, we will further and in more detail quantify how the local fluid properties affect the shrinkage rate of the droplets.

B. Lagrangian Statistics of the Droplets.
For the following Lagrangian statistics, we focus on the droplets with initial diameter around 15µm, given that for those we observe the most pronounced temperature effect, see the count histogram in Fig. 2B. First, we track the droplet surface area ∝ d 2 versus time as shown in Fig. 3A,B. Note that droplets are expelled from the mouth at different time instants until the respiratory event stops (0.6s for the coughing event considered here). In order to compare the droplet statistics, we have shifted the time frame by their respective ejecting times.
For θ amb = 30°C (Fig. 3A), the normalized surface area of the droplets d 2 /d 2 init decreases monotonically with time (the ensemble average of the value is depicted by the solid curve), which indicates pure evaporation of the expelled droplets. The surface area decreases linearly with time with the averaged droplets surface area being halved after 1s. Although a linear decline of the surface area is consistent with the classical d 2law (24), one should note that the magnitude of the droplet shrinkage rate is still much smaller than that predicted by William F. Wells (22). Such a discrepancy is caused by the fact that the turbulent vapour puff engulfs small droplets and the puff contains higher RH than the ambient fluid (3,27).
Surprisingly, for the low θ amb of 10°C (Fig. 3B), the res- piratory droplets evolve in a very different and distinct way: Instead of the pure evaporation during the time evolution of the droplets, there exists an initial stage (from the expelled time to about 0.4s) during which the droplets can grow with the averaged d 2 increasing by 10%. This non-monotonic trend of the shrinkage markedly differs from the classical d 2 -law. It also contradicts the well-accepted view that the respiratory droplets should simply evaporate after expulsion. Such a phenomenon of droplets condensation during the respiratory event is the first time, to the best of our knowledge, to be quantified thanks to the Lagrangian statistics obtained from DNS, which has been neglected in the hitherto studies of droplet transmission in exhalation. Another indicator that for 10°C the droplet size first increases is that the droplet temperature θ drop increases relative to the fluid temperature at the location of the droplet θ fluid (Fig. 3D). As latent heat is released during droplet condensation for 10°C (Fig. 3D), the difference θ drop − θ fluid increases during the early time. In contrast, θ drop − θ fluid remains less than zero for θ amb = 30°C (Fig. 3C) needed for evaporation. Eventually, the droplet temperature saturates at a value lower than the fluid temperature. This is reasonable because the Jakob number, which is the ratio of the sensible heat to the latent heat due to the phase change, is much smaller than one (0.04 and 0.007 for θ amb = 10°C and 30°C, respectively). Therefore, the evaporation cooling effect plays a significant role here.
To understand why the droplets can grow during the early time, a key quantity to examine is the local RH experienced by the droplet. In Fig. 3E,F, we compare the time evolution of the local RH at the droplet location for the two different θ amb , showing quite different behaviours: For 30°C, the local RH decreases monotonically to RH amb which is 90% in this case. However, for 10°C, the droplets first experience super-saturation within 0.4s after expelling and the maximum RH reaches almost 115%. This duration of having supersaturation is consistent with the duration of droplet growth, which reflects that it is the local super-saturation of the surrounding air which makes the droplets grow.
In fact, our daily life experience gives us intuition about this intriguing droplet condensation phenomenon. In cold and humid weather, one may observe the "white mist" coming out from the mouth while breathing or speaking. The physical explanation is that: Warm air can contain more moisture than cold air. Therefore, the exhaled warm vapour puff becomes super-saturated when it enters the cold ambient fluid, and the droplet nuclei and dust in air favours vapour condensation. Such super-saturated conditions become more prominent if the temperature difference between the exhaled vapour and ambient fluid is larger. Here, we have clearly demonstrated the significance of super-saturation in the respiratory event that leads to the possibility of vastly different droplet dynamics.

C. Theoretical Model Predicting Relative Humidity Profiles
and Super-saturation Criteria. Motivated by the strong effect of the local RH field on the droplets, here, we propose a simple model to calculate the local RH. Several assumptions first need to be made to justify this model. First, the axial scalar quantities (vapour mass fraction and temperature) are the dominant fluid properties that determine droplet growth/evaporation, and secondly, the exhaled puff at the time of t = 0.2s exhibits jet-like properties such that the quasistationary mean properties of the puff admits self-similarity.
Here, we employ the Antoine equation in order to compute the saturated vapour mass fraction for a given temperature field (34). In short, given the local values of vapour mass fraction and temperature along the jet centreline, the axial RH can be directly computed. Following the theory for axisymmetric turbulent jets (35)(36)(37), the axial scalar concentration distribution (say Φ axial ) can be approximated by the relation, [1] where Φ inlet is the inlet scalar concentration value, and x/d inlet is the axial distance normalised by the inlet diameter. Since this model diverges as x/d inlet → 0, we further impose the requirement that max[Φ axial /Φ inlet ]= 1. Both the exhaled vapour mass fraction and temperature are modelled using Eq. (1) along the axis normal to the inlet.
To assess the validity of our assumptions, we plot the resulting model of the axial RH profiles for different θ amb in Fig. 4. The agreement with our DNS (lighter curves) is quite promising, which is consistent with our assumption that the puff retains self-similar jet-like characteristics at cut-off time. This agreement lends some confidence to our ability to model the axial RH of our flow at a specific time instant, while the cough will approximate to a turbulent puff or an interrupted jet at large t after the injection velocity becomes negligible (38).
Based on the robustness of this model, we further propose the use of the maximum local RH as a metric to indicate the likelihood of the growth for small droplets. To illustrate this point, in Fig. 5, we plot the axial RH map for various θ amb versus distance at the fixed RH amb of 90%. The dashed line shows the crossover between super-saturated and undersaturated axial RH values. One interesting observation from

D R A F T
this model is that with decreasing θ amb , the overall axial distance that experiences super-saturation increases. In other words, it is more likely that with decreasing θ amb (at RH amb = 90%) small droplets will encounter longer super-saturated distances and grow in size. This likelihood to grow is important in understanding how far the droplets can be advected by the puff in the initial stages of the cough, before evaporating and dispersing.
Finally, we extend our idea of modelling the axial RH to map out the maximum local RH values for a range of RH amb and θ amb . This map is shown in Fig. 6. From the figure, we find that the maximum local RH values are dependent on both RH amb and θ amb . Several interesting observations can also be made from the figure. For instance, at a given indoor ambient conditions of RH amb = 40-50% (within acceptable occupancy comfort levels), the maximum local RH is at most 100% as long as θ amb > 20°C. However, once θ amb < 20°C, the local RH tends to super-saturate (RH > 100%), implying that in cooler environments, smaller droplets within the humid puff will first undergo supersaturation-driven growth in the early stages of the cough before evaporating upon fall-out. Such conditions become particularly exacerbated in outdoor environment, for example, in stadiums or sports field during autumn or winter, when the temperatures tend to be colder. Similar conditions can also occur in crowded bus and train stations. Super-saturated conditions are more easily attainable (easily seen in exhaled mists) and therefore alters the respiratory droplet dynamics.  Fig. 5 (red dashed curve). The values can be used as a gauge to indicate the tendency of droplets to experience super-saturation, given that the droplets are caught within the turbulent puff. The color is coded by the upper bounds for super-saturation. The plot indicates that in cool and humid environments, small droplets within a coughing puff will tend to grow once expelled.
To assess the performance of our model, we also plot our DNS results for droplets with an initial size of ≈ 15µm and indicate growth or shrinkage using circle and cross symbols, respectively. (The detailed breakdown is plotted in Fig. S1.) A good agreement between our DNS and the RH model is seen.
To substantiate our proposed model, in Fig. 6, we indicate on the plot 12 different control parameter combinations (RH amb and θ amb ) for which droplets with an initial size of ≈ 15µm grow (indicated by a circle symbol) or shrink (indicated by a cross symbol) in our DNS and compare with the predictions of the model. (The profiles of d 2 versus time is provided in Fig. S1 of the Supplementary Material). By com-paring the symbols from the DNS and the map from the model, a clear correlation between growth and super-saturated local RH fields can be seen. Thus we propose that this map can be used as an indicator to determine the droplet size distributions of a similar cough.
It is also worth noting that we did not consider the influence of background wind, which is an important aspect both in ventilated indoor and outdoor environments. We also stress that we make no claims to correlate growth or shrinkage of droplets with virus transmissions or infectivity. The subject of the viability of droplet sizes for viral media is still the topic of ongoing investigations in other communities and are not the focus of this work.

Conclusions
Through the Lagrangian statistics of the respiratory droplets, we have identified that when droplets are expelled from the mouth to the cold environment with high enough relative humidity, the droplet will first grow instead of immediately shrink, because the turbulent vapour puff becomes super-saturated. More importantly, we have provided a theoretical framework to accurately predict (i) how the relative humidity varies with distance in a respiratory puff, and (ii) the threshold of ambient temperature (θ amb ) and relative humidity (RH amb ) for which the super-saturated vapour field can result. Our theoretical model can also be applied to other respiratory events whenever there is jet-like transport, such as speaking (25), which is extremely relevant for asymptomatic and presymptomatic spreading of the coronavirus. In the situation of continuous speaking, when the hot vapour is continuously injected into the cold environment, we expect that the super-saturation zone may remain in front of the mouth for a long period of time. Our finding implies that tiny droplets of the size of tens of micrometers can propagate further away when the weather is cold and humid due to the protection from the super-saturated vapour field. The result of the very distinct fate of respiratory droplets under different RH amb and θ amb may be instructive for developing timely strategies in mitigating the COVID-19 pandemic in different seasons. For example, during the Winter time, where the maritime climate for many coastline European countries have much higher RH amb , the transmission in outdoor environments can presumably not be overlooked. However, we note that obviously also non-fluid-dynamical aspects affect the transmission and infectivity of viruses, such as from virology, epidemiology, and sociology, to name a few. These are not covered within our study and presumably interfere with our fluid dynamical results.

Materials and Methods
Governing equations. We consider an incompressible flow (∂u i /∂x i = 0) of gas phase, with both temperature and vapour concentrations coupled to the velocity field by employing the Boussinesq approximation. The governing equations read: hmAn(θg,n − θn)δ( x − xn), [3] ∂c ∂t . [4] For droplets, we apply the spherical point-particle model, and consider the conservation of momentum (Maxey-Riley equation (39)), energy, and mass as follows: r 2 n f d +gβêy, [5] ρ l c p,l Vn dθn dt = ρ l AnL drn dt + hmAn(θg,n − θn), [6] drn The notations we used in equations are as follows: u i , u i,n , and u i,g,n are the velocity of gas, droplets, and gas at the location of droplets, respectively. θg, θn, and θg,n are in Kelvin and used to represent the temperature of gas, droplets, and gas at the location of droplets, respectively. c is the vapour mass fraction, rn the droplet radius, An the surface area of the droplets, Vn the volume of the droplets, mn the mass of the droplets. Also p denotes the reduced pressure. hm is the heat transfer coefficient. L is the latent heat of vaporisation of the liquid. ρ l and c p,l are the density and specific heat capacity of the droplets assumed to consist of water. ν air is the kinematic viscosity of air. kg is the thermal conductivity of gas which relates to the thermal diffusivity of gas Dg by kg ≡ Dgρgcp,g. ρg and cp,g are the density and specific heat capacity of gas, Note that ρgcp,g = (ρacp,a + ρvcp,v), with ρa and ρv being the densities of air and vapour and cp,a and cp,v being the specific heat capacities of air and vapour. c drop and c f luid denote the vapour mass fractions of the droplet and the fluid at the location of the droplet. β is a dimensionless measure of the droplet density relative to the fluid density, and is defined as β ≡ 3ρg/(ρg + 2ρ l ) − 1. f d is the prefactor for the drag corrections defined as f d = 1 + 0.169Re where we have droplet Reynolds number Re drop = |u i,g,n − u i,n |(2r) ν air . [10] To calculate the local relative humidity, the saturated vapour mass fraction csat,vap is determined by ideal gas law csat,vap = Psat/(ρgRθg), where R is the specific gas constant of water vapour and Psat is the saturated vapour pressure, determined through Antoine's relation as: P sat (θg) = 10 5 exp 11.6834 − 3816.44 226.87 + θg − 273. 15 . [11] Simulation details. The computational domain size is 0.37 m (spanwise length) ×0.37 m (height) ×1.47 m (streamwise length) with corresponding grid points chosen as 512 × 512 × 2048. We model the mouth as a circular inlet centred at mid-height of the domain and the cough temporal profile we apply is a gamma distribution as U cough (t) = U cough αt exp(−αt/4), [12] where α = 60.9 s −1 , such that the entire cough process lasts for about 0.6 s. Based on an experimental measurement (13), we seeded the respiratory event with droplets with initial diameters ranging between 10 µm and roughly 1000 µm. We employ a similar droplet size distribution as (13,26) with 5000 droplets. The droplets are randomly positioned at the inlet area and evenly injected in time with the velocity matching the local inlet velocity. Droplets are assumed not to collide or coalesce since the estimates of the droplet volume fractions give O(10 −6 ). These droplets also do not set the relative humidity of the puff, which at inflow is assumed to be 100%.