Water evaporation from solute-containing aerosol droplets: Effects of internal concentration and diffusivity profiles and onset of crust formation

Evaporation from solute-containing droplets is an important process in many industrial and scientific applications ranging from pharmacology, agriculture, food, and cosmetics production to medical, biochemical, material, and soil sciences.^7–12 Depending on the rate of the solvent evaporation from the droplet surface and the diffusive transport rate of solutes inside the droplet, the drying process can lead to particles with different morphologies and thus different properties and functionalities. It is, therefore, of importance to investigate the morphological evolution during evaporation, which is affected by solute diffusion inside the droplet. These processes are controlled by various parameters, such as temperature and relative humidity of the ambient air, the physical–chemical properties of solvent and solutes, the initial volume fraction of solutes, and the initial size distribution of droplets. Theoretical investigations are, therefore, required to explain how the droplet evaporation kinetics depends on all the relevant parameters.

In addition, the droplet-evaporation process plays an important role in the transmission of infectious diseases and respiratory viruses via the airborne route,^6,13–16 which is the main motivation of the present study. It is known that saliva is primarily composed of water, but additionally includes a variety of organic and inorganic substances such as salt, proteins, peptides, mucins, virions, etc.¹⁷ The presence of such nonvolatile components produces a lower limit for the water content of a saliva droplet that can be reached by evaporation. Accordingly, water evaporation from a saliva droplet might eventually produce a droplet nucleus, which is a small light particle with a minimal moisture level that stays floating in air for a long time.¹⁸ The creation of such droplet nuclei, which is more likely in the case of small droplets due to their sufficiently long sedimentation time,^19–21 can significantly influence the infection risk from virus-containing respiratory droplets, especially in indoor environments,^14,18 by increasing the sedimentation time of droplets.⁶ It is, therefore, of utmost importance to investigate how respiratory droplets dry out and how their evaporation kinetics is affected by different factors such as the droplet composition and its initial size, the ambient relative humidity, temperature, and the surrounding air flow.

Although there are many studies^11,22–25 concerning the drying process of droplets placed on solid surfaces, also known as sessile droplets, fewer investigations have been conducted on free droplets surrounded by air. Indeed, designing and conducting experiments on free droplets faces many challenges, mostly since such droplets are ephemeral and difficult to follow.⁸ So far, a few experiments have been conducted and are used alone or in combination with theoretical analysis to measure time-dependent properties of evaporating or condensing free droplets, such as their size,^26,27 viscosity,^26,28 and hygroscopicity,^29–31 as well as the solvent diffusion coefficient inside droplets.^32,33 Also, the morphology of aerosols has been investigated using different experimental methods.^27,34 Various physical–chemical effects and mechanisms that affect the droplet evaporation process must be accounted for in the interpretation of these experiments. Therefore, analytical and numerical modeling methods must be used to better understand the details of the evaporation mechanism that is relevant for sedimenting droplets. Although modeling of free droplets does not face complexities due to substrate–droplet interactions that control the droplet shape at the solid/liquid/air interface^35,36 and Marangoni flow,^37,38 which occur in the case of sessile droplets, other factors such as evaporation-induced concentration gradients inside the droplet³⁹ and the possibility of crust formation on the droplet surface,⁴⁰ which are consequences of the increasing solute surface concentration during evaporation,⁴¹ cause difficulties even in the modeling of free droplets in the presence of solutes. In addition, physical and chemical properties of the drying droplets, such as the internal viscosity,^42,43 the diffusivity of solvent and solutes in the liquid phase,⁴³ and the activity coefficient of the solvent,⁴⁴ are dependent on the local concentration of solutes (and, consequently, on both position and time), which makes the problem rather complex. A challenge in modeling evaporating droplets is, therefore, to account for the concentration gradients created and developed within droplets during the solvent evaporation process. Such concentration gradients are not only relevant for solutions containing slowly diffusing species, e.g., polymers and proteins, but also for relatively quickly diffusing solutes such as NaCl.⁸ In fact, internal concentration gradients have been experimentally observed for drying respiratory fluids suspended on superhydrophobic substrates using optical and fluorescence microscopy.⁴⁵ 

There are a few studies that address concentration gradients inside a solute-containing drying droplet^39,46–48 and propose analytical and numerical methods for modeling the evaporation process before^47–50 and after^40,49,51 crust formation. Handscomb et al.⁵² proposed a numerical model consisting of a set of advection–diffusion equations coupled to ordinary differential equations describing the particle size and temperature to simulate the drying of droplets that contain solid particles prior to crust formation. Solute nucleation and crust formation in drying free droplets were numerically modeled by Robinson et al.⁵³ using a single-particle approach based on the diffusion equation. Diffusion through atmospheric aerosols has been numerically simulated using core–shell multilayer models such as the KM-GAP model proposed by Shiraiwa et al.⁵⁴ based on the PRA (Poschl-Rudich-Ammann) framework,⁵⁵ the model proposed by Fowler et al.⁵⁶ based on the Maxwell–Stefan diffusion framework, and the Fi-PaD model proposed by O'Meara et al.⁵⁷ based on the Fickian diffusion framework. The latter two models can be used for droplets containing both non- and semi-volatile solutes. For example, Ingram et al.²⁶ simulated evaporation kinetics of droplets that contain water and semi-volatile organic compounds using the Fi-PaD model. Also, a few analytical or semi-analytical approaches have been proposed to treat solute diffusion inside a drying droplet, most of which use a fixed droplet shrinking rate. In a series of works by Sazhin et al.⁵⁸ and Gusev et al.,⁵⁹ it was assumed that the droplet radius decreases linearly in time. This assumption, however, holds only for very small droplets in the so-called reaction-rate-limited regime.⁶ In fact, for industrial droplet-drying processes and airborne transmission of infections, the initial radius of the drying droplets is typically larger than the threshold radius below which the reaction-rate-limited scenario is valid [$∼68 nm$ for pure water droplets at 25 °C (Ref. 6)] In other droplet evaporation models,^46,60–64 the analytical calculations are based on the so-called R²-law, which assumes that the droplet temperature and thereby also the droplet evaporation rate remain constant during the evaporation process and, thus, the droplet surface decreases linearly in time.⁶⁵ Although this approximation describes the evaporation-induced shrinkage of pure solvent droplets accurately, it cannot predict the osmotic slowing down of the evaporation process due to the solute-induced water vapor-pressure reduction.⁶ We show in this paper that the osmotic effect plays a key role in determining the droplet evaporation time. Deviations from the R²-law were accounted for in a model developed by Eslamian et al.,^66,67 where the droplet evaporation rate was calculated from a combination of Fick's law⁶⁸ and kinetic theory.⁶⁹ The equations provided by this model are, however, not analytically solvable. Our model describes the evaporation process of solute-containing droplets and, in particular, accounts for the solute-induced deviation from the classical R²-law model.

The evaporation-induced concentration gradients inside drying droplets are found to strongly affect the drying process not only by decreasing the evaporation rate,⁴⁸ but also by influencing the morphological evolution of the drying droplet^40,70 and, in particular, the shape and properties of the particles produced at the end of the drying process. An important physical parameter is the ratio of the evaporation rate to the diffusive transport rate of solute particles inside the droplet,⁷¹ which is known as the Péclet number.^64,72 When the evaporation process is very slow or, alternatively, the internal diffusion is very fast, the solute particles have enough time to redistribute during the drying process. In this case, the solute particles remain evenly distributed within the droplet to form a full solid particle at the end of the drying process [Fig. 1, scenario (a)]. In the case of fast evaporation, however, the solvent evaporation from the droplet surface increases the solute concentration at the surface and creates a concentration gradient in the droplet [see Fig. 1, scenarios (b)–(f)]. If the critical supersaturation concentration at which phase separation occurs is not reached during evaporation, the concentration gradient gradually disappears as the evaporation proceeds and a full solid particle remains at the end of the drying process [Fig.1, scenario (b)]. When the critical supersaturation concentration is reached during evaporation, however, two-phase coexistence is obtained due to radial de-mixing³⁴ and the solute particles at the droplet surface form a solid crust [see Fig. 1, scenarios (c)–(e)], which is expected to mainly affect the drying mechanism and to considerably slow down the evaporation process.⁷⁰ The fate of the crust and the physical and morphological properties of the dried particle formed at the end of the evaporation process are affected by various parameters such as droplet diameter, initial solute concentration, solvent composition, and drying rate.²⁷ 

Depending on the type of solutes, the crust that forms might be either a dry crust (expected for salty droplets)^8,73 or a gel-like wet skin consisting of polymers or proteins and other suspended particles.^74,75 In the presence of a dry crust, the liquid from the internal liquid core reaches the crust surface through capillary action within the crust pores.^70,76 As the evaporation proceeds, the crust will become thicker and the crust pores will lengthen, which increases the resistance to heat and mass transfer and decreases the evaporation rate. In the case of a wet crust, water evaporation continues via diffusion through the gel crust.⁷⁴ Since the diffusion coefficient and the concentration gradient within the gel phase are both much lower than in solution, the evaporation rate is expected to fall after the formation of a wet crust as well. Regardless of the crust type, crust formation is expected to determine the final morphology of the droplet by producing a hollow structure, with a size larger than expected in the absence of a crust. Depending on the type of the solute particles within the droplet, the hollow structure produced at the end of the evaporation process may or may not contain small agglomerates of solutes⁴⁰ [Fig. 1, scenarios (c) and (d), respectively]. Also, a crust collapse might occur when the crust cannot withstand the pressure difference caused by the continued solvent evaporation⁷⁰ [Fig. 1, scenario (e)]. Another scenario accounts for non-uniform drying conditions due to amplification of concentration inhomogeneities over the droplet surface, which is probable when the solvent evaporation is very fast. In such case, an irregularly shaped dry particle, also known as a wrinkled particle, will form at the end of the drying process⁷⁰ [Fig. 1, scenario (f)].

In this paper, a droplet evaporation model is presented by solving the heat and mass diffusion equations analytically or semi-analytically to describe the various effects nonvolatile solutes have on the drying process. Starting from the simple model scenario that only accounts for the solute-induced minimal water concentration inside the droplet, we gradually add complexities to our model to address different factors associated with the presence of solutes, i.e., the solute-induced water vapor-pressure reduction, effects of solutes on evaporation cooling, evaporation-induced water concentration gradients inside the droplet, and the dependence of the internal water diffusivity on the solute concentration. We also include the possibility of crust formation in our calculations to evaluate its effect on the droplet equilibrium size. All our calculations are done in the one-phase regime, meaning from the beginning of the evaporation process up to the point where crust formation sets in at the surface, and for a single droplet at rest, which defines the so-called stagnant-flow approximation and is valid for droplet radii up to 60 μm.⁶ The time-dependent solute and solvent densities of the droplet are described in terms of the time-dependent solute volume fraction of the droplet, assuming ideal mixing of solute and solvent. Our results indicate that the presence of solutes causes a reduction of the evaporation rate and thereby also of the evaporation-cooling of the droplet not only by reducing the water vapor pressure, but also by inducing a locally reduced water concentration inside the droplet near the surface, which arises from the finite diffusivity of water and solutes in the droplet. Any factor that reduces the internal water and solute diffusivity, such as the presence of strongly hydrated solutes inside the droplet, intensifies this effect and thus further slows down the evaporation process. The time-dependent local variation of the internal water diffusivity, which occurs especially near the droplet surface due to the increased solute concentration, is, however, found to play only a minor role in determining the droplet evaporation time. Crust formation is suggested to affect the morphology of the final droplet nuclei by producing a hollow particle with a size larger than expected in the absence of a crust. The size of the hollow particle formed after crust formation follows approximately from our analytical model.

The main objective of the present study is to model water evaporation from a solute-containing aqueous droplet that is sedimenting in air. For this, the time-dependent water–vapor concentration profile outside the droplet and the liquid–water and solute concentration profiles inside the droplet are determined from the diffusion equation, which for spherical symmetry depends on the radial coordinate and reads

where the subscript $i$ denotes different molecular species (i.e., water vapor outside the droplet, liquid water, and solutes inside the droplet), $cir,t$ is the concentration profile, and $Di[cir,t]$ is the respective molecular particle diffusion coefficient that in general depends on the concentration and, consequently, on time and position. As discussed in supplementary material Sec. B.1, the coupled diffusion equations for solvent and solutes within an incompressible solution are simultaneously solvable only when the diffusivities of the two components are the same. In the presence of slowly diffusing solutes, such as polymers or proteins, the water dynamics becomes slaved to the solute dynamics and thus one can approximately describe the diffusion dynamics by a single diffusion equation with a reduced effective diffusion coefficient, which is the approach we are taking. The exact problem is much more complicated and requires the solution of two coupled diffusion equations while accounting for the finite solution compressibility.

An important factor that affects the droplet evaporation rate is the evaporation-induced cooling of the droplet, which decreases the temperature and thereby also the water vapor concentration at the droplet surface. To account for this, one needs to calculate the temperature distribution around the evaporating droplet from the heat-conduction equation in radial coordinates,

with $λ$ and $CP$ being the heat conductivity and the heat capacity of air, respectively. In the present study, we ignore temperature gradients inside the droplet and thus the heat-conduction equation is solved only outside the droplet.

In the following, the results obtained from solutions of the above equations at various approximation levels are discussed. Derivations of all equations are provided in the supplementary material.

For later comparison with our main results, we first assume that the particle diffusion inside the droplet is sufficiently rapid, so that the water concentration at the droplet surface equals the mean water concentration in the droplet. If one furthermore assumes that the water diffusion in the vapor phase is fast so that at each instant of time a stationary water–vapor concentration profile is obtained, the time-dependent variation of the droplet radius due to the water evaporation results from the stationary solution of the coupled water–vapor-diffusion and heat-conduction equations [Eqs. (1) and (2)] outside the droplet, considering total water mass conservation. As detailed in supplementary material Sec. A, the droplet radius R is determined by the differential equation

where $ψc$ is a radius-dependent factor that accounts for evaporation cooling (as described below), $Φ$ is the momentary volume fraction of solutes, $γ$ is the water activity coefficient that accounts for non-ideal effects due to water–solute and solute–solute interactions (we assume here that $γ$ does not depend on the water concentration), $Dw$ is the water diffusion constant in air, $cg$ is the saturated water vapor concentration, $vw$ is the liquid water molecular volume, $RH$ is the relative humidity, and $Rev$ is the equilibrium droplet radius reached by evaporation. Neglecting the possibility of crust formation due to phase separation at the droplet surface (which will be discussed later), the equilibrium droplet radius reads (see supplementary material Sec. A.4)

with $R0$ and $Φ0$ being the initial droplet radius and the initial volume fraction of solutes, respectively. The cooling factor $ψc$ in Eq. (3), which describes the slowing down of the evaporation process due to evaporation cooling, is given by (see supplementary material Sec. A.3)

where $ΦR=Φ0R03/R3$ is the momentary volume fraction of solutes, $εC=0.032 K−1$ is a numerical coefficient that accounts for the reduction of the water vapor concentration at the droplet surface due to the temperature depression [see Eq. (A19)] and $εT≡Dwcghevλair$⁠, with $hev$ being the molecular evaporation enthalpy of water and $λair$ the heat conductivity of air, controls the dependence of the temperature reduction at the droplet surface on the relative humidity and the solute volume fraction [see Eq. (A22)]. In the absence of solutes, Eq. (5) gives rise to a constant cooling factor, $ψcΦ0→0=1+γεcεT−1$⁠. At a room temperature of 25 °C, one obtains $εT=54$ (see the Nomenclature) and, thus, $ψcΦ0→0≈0.36$⁠, demonstrating that evaporation cooling considerably slows down the evaporation process of a pure water droplet.

To make the differential equation (3) analytically solvable, one can neglect the effect of solutes on the evaporation-cooling factor, which corresponds to replacing $ψcΦR$ by $ψcΦ0→0$ in Eq. (3). From an asymptotic analysis of the solution of the differential equation (see supplementary material Sec. A.4 for more details), the radius-dependent evaporation time $tR$⁠, i.e., the time it takes for the droplet radius to decrease from its initial value $R0$ to $R$⁠, can be written as

with $θ=2γDwvwcg1+γεcεT−1$ being a prefactor with units of a diffusion constant. We will show later that for droplets with rather low initial solute volume fractions, the solute-concentration dependence of evaporation cooling has only a minor effect on the droplet evaporation time and thus the approximation provided by Eq. (6) is rather accurate.

The logarithmic term in Eq. (6) reflects the osmotic slowing down of the droplet evaporation due to the solute-induced water vapor-pressure reduction. Figure 2(a), which shows $tR$ with and without considering the logarithmic term, reveals that this osmotic effect only becomes relevant for droplet radii close to the final equilibrium radius $Rev$⁠, where the droplet has lost most of its water content and solute effects cause a diverging evaporation time. This effect is, however, negligible for droplets with low initial solute volume fractions [see Fig. 2(a) and inset]. In this case, one can safely neglect the logarithmic term in Eq. (6), which results in the classical R²-law that is defined by $R2=R02−κt,$⁶⁵ with $κ$ being the shrinkage rate of the droplet surface, which according to Eq. (6) equals $θ(1−RHγ)$⁠. The approximate expression for the droplet evaporation time without the logarithmic term and using $R=Rev$ in Eq. (6) results in the evaporation time

where $τev$ in Eq. (7) is given by

and for ideal solutions (⁠$γ=1)$⁠, i.e., when the solute volume fraction is negligible, represents the droplet evaporation time in the absence of solutes.⁶ The bar in Eq. (7) indicates that this equation is obtained assuming that the internal concentration profile remains homogeneous throughout the evaporation process.

As discussed above, Eq. (7) neglects a few important details of the evaporation process, such as the solute-induced water vapor-pressure reduction and the effect of solutes on evaporation cooling. The presence of an internal water concentration gradient, which arises from the finite diffusivity of water and solutes inside the droplet, is also neglected. Not only can this concentration gradient affect the droplet evaporation rate, it can also lead to formation of a solid crust on the droplet surface, which affects both the evaporation rate and the equilibrium droplet radius. The importance of these factors is evaluated next.

To account for the solute-induced water–vapor pressure reduction, the logarithmic term in Eq. (6) is now included. Since near the evaporation equilibrium, the droplet radius varies very slowly with time, as follows from the logarithmic term in Eq. (6) [see Fig. 2(a)], we can no longer use the definition for the evaporation time that leads to Eq. (7) (i.e., the time at which $R=Rev$⁠). Instead, we define the evaporation time as the time at which the equilibrium droplet radius $Rev$ reaches 99% of the droplet radius $R$

From Eq. (6) and considering the definition provided by Eq. (9), the evaporation time of a solute-containing droplet follows as

Again, the bar on top of symbols indicates the assumption that the water concentration inside the droplet remains homogeneous throughout the evaporation process (evaporation-induced concentration gradients inside the droplet are discussed later). Also, the logarithmic term $τ¯ev,slog$ in Eq. (10) reflects the increase in the evaporation time due to the solute-induced water vapor-pressure reduction, which causes a distinct deviation from the R²-law,

According to Fig. 2(b), the droplet evaporation time $τ¯ev,s$ [see Eq. (10)] and its logarithmic contribution $τ¯ev,slog$ [see Eq. (11)] both exhibit non-monotonic variations as a function of the initial solute volume fraction $Φ0$⁠, while the non-logarithmic term $τ¯ev,snon−log$ [see Eq. (7)] decreases with increasing $Φ0$⁠. The reason for this can be traced to solute effects on $τ¯ev,snon−log$ and $τ¯ev,slog$⁠. As follows from Eq. (4), an increase in $Φ0$ increases the droplet equilibrium radius $Rev$⁠, i.e., the droplet radius at which the evaporation process stops, which tends to decrease both $τ¯ev,snon−log$ and $τ¯ev,slog$⁠. An increase in $Φ0$ additionally intensifies the solute-induced divergence of the logarithmic term in the function $t(R)$ given by Eq. (6) [see Fig. 2(a) and inset], which tends to increase $τ¯ev,slog$⁠. Figure 2(b) reveals that at low values of $Φ0$⁠, the second effect is dominant and thus both $τ¯ev,slog$ and $τ¯ev,s$ increase with increasing $Φ0$⁠. At higher values of $Φ0$⁠, however, the first effect becomes dominant and, consequently, $τ¯ev,s$ and $τ¯ev,slog$ begin to decrease with increasing $Φ0$⁠. By inserting $Rev$ from Eq. (4) into Eq. (10) and differentiating the resulting equation with respect to $Φ0$⁠, the value of $Φ0$ at which the evaporation time is maximal is obtained as

where $Φev=1−RH/γ$ is the volume fraction of solutes at evaporation equilibrium, i.e., when the evaporation flux [given by Eq. (A24)] vanishes. Accordingly, when the solute volume fraction inside the droplet is initially around 45% of its final equilibrium value, the droplet experiences its maximally possible evaporation time, which according to Eqs. (4) and (10) is around 1.84 times longer than the evaporation time of a pure water droplet with the same initial radius.

Figure 2(b) also indicates that at low values of $Φ0$⁠, $τev,sollog$ is negligible compared to $τev,solnon−log$⁠, meaning that the R²-law is almost valid and thus Eq. (7) is accurate. In the case of droplets with high initial solute volume fractions, however, the logarithmic terms in Eqs. (6) and (10) are important and can no longer be neglected. By equating Eqs. (7) and (11), one finds that the effect of the water vapor-pressure reduction becomes dominant when $Φ0$ becomes higher than approximately $0.13Φev$⁠.

The evaporation-induced cooling of the droplet decreases the temperature at the droplet surface compared to the ambient air by $ΔT=Ts−T0$⁠, where $Ts$ and $T0$ are the temperatures at the droplet surface and far away from the droplet, respectively. The value of $ΔT$ is obtained by solving the coupled heat-conduction and water–vapor-diffusion equations outside the droplet (see supplementary material Sec. A.3) as

where $ψcΦR$ is the cooling factor given by Eq. (5). Equation (13) reveals that not only $ΔT$ is linearly related to the ambient relative humidity $RH$⁠, as already reported in the absence of solutes,⁶ but it also depends on the momentary volume fraction of solutes $Φ$⁠, which increases with time due to the water loss by evaporation. This equation shows that at evaporation equilibrium, i.e., when the momentary volume fraction of solutes reaches its equilibrium value, $Φev=1−RH/γ$⁠, $ΔT$ goes to zero, meaning that the droplet reaches thermal equilibrium with its environment as well.

So far, we used Eq. (10) to calculate the radius-dependent evaporation time. As mentioned before, this equation neglects the effect of solutes on evaporation cooling, which corresponds to using the cooling factor in the absence of solutes, $ψcΦ0→0$⁠, in the calculations. According to Eq. (13), this approximation results in a linear decrease in $ΔT$ with an increase in the momentary solute volume fraction $Φ$⁠. The solute-concentration dependence of the evaporation-cooling factor, however, causes a non-linear relationship between $ΔT$ and $Φ$⁠, as follows from Eqs. (5) and (13). To account for this effect, one has to calculate the shrinkage rate of the droplet radius directly from Eq. (3), which has to be done numerically. Figure 3(a) shows the time-dependent droplet radii calculated from the numerical solution of Eq. (3), with and without considering the effect of solutes on evaporation cooling {corresponding to using $ψcΦR$ and $ψcΦ0→0$ in the calculations, respectively}. This figure indicates that the solute-concentration dependence of evaporation cooling slightly speeds up the evaporation process. This speed up becomes more significant as the evaporation proceeds, which can be attributed to the increased solute concentration within the droplet. This effect is found to be more pronounced for droplets with higher initial solute volume fractions [see Fig. 3(a) and inset].

Figure 3(b) shows the time-dependent variation of the evaporation-induced temperature reduction at the droplet surface $ΔT$ in the absence of solutes, and in the presence of solutes with and without considering solute effects on evaporation cooling. This figure shows that the temperature reduction at the surface of a pure water droplet remains constant throughout the evaporation process. In the presence of solutes, however, the temperature reduction decreases with time and eventually reaches zero at evaporation equilibrium. Figure 3(b) reveals that neglecting the effect of solutes on evaporation cooling slightly reduces the time-dependent decay of $ΔT$⁠, which causes a delay in the evaporation equilibrium state.

To better understand how the solute-concentration dependence of evaporation cooling affects the droplet evaporation time, the evaporation time is now calculated from the numerical solution of Eq. (3), using the definition provided by Eq. (9). The results obtained from the numerical calculations with and without considering the solute-concentration dependence of evaporation cooling {which correspond to using $ψcΦR$ and $ψcΦ0→0$ in the calculations} are, respectively, shown by solid and broken lines in Fig. 3(c). Also, dotted lines in this figure show the results obtained from Eq. (10), which is an approximate analytical solution of Eq. (3) when the effect of solutes on evaporation cooling is neglected. This figure clearly shows that the solute-concentration dependence of evaporation cooling causes a decrease in the droplet evaporation time, especially when $Φ0$ becomes close to $Φ0τ max$ [see Eq. (12)]. The approximation used in the derivation of Eq. (10) [which is described by Eq. (A32)] is, however, found to partly compensate for the difference caused by neglecting such concentration dependence, especially at high relative humidity. According to this figure, at low initial solute volume fractions $Φ0≤0.01Φev$⁠, which for $RH≤0.5$ includes the normal solute volume fraction of saliva, $Φ0≃0.05,$⁷⁷ one can safely use Eq. (10) to estimate the droplet evaporation time. At higher values of $Φ0$⁠, however, this approximation starts to become inaccurate. Figure 3(c) also indicates that the effect of the solute-concentration dependence of evaporation cooling is more significant in dry environments.

In the previous section, the particle diffusivity inside the droplet was assumed to be sufficiently rapid, so that the solute remains evenly distributed. In the fast evaporation scenario (see Fig. 1), however, the droplet exhibits an increased solute concentration at its surface, which results in a water concentration gradient inside the droplet. To account for this concentration gradient, one needs to solve the diffusion equation [Eq. (1)] not only in the vapor phase outside the droplet, but also in the liquid phase for both water and solutes. As discussed in supplementary material Sec. B.1, the diffusion equations for solvent and solutes in an incompressible fluid mixture are simultaneously solvable only if the diffusion coefficients of all components in the liquid mixture are equal. Here, we use this assumption, which approximates the behavior of fast-diffusing solutes rather accurately. Accordingly, the liquid water diffusion coefficient $Dwl$ is used to calculate the concentration profile inside the droplet. We will also show results for diffusion coefficients smaller than $Dwl$⁠, which is an approximate way of treating the effect of slowly diffusing solutes, as discussed before.

The general problem one encounters when considering diffusion in radial coordinates is that the diffusion equation [Eq. (1)] has a singularity at $r=0$⁠, at the center of the droplet. The easiest way to deal with this problem is to solve the discretized diffusion equation using numerical methods,⁷⁸ where the singularity is removed from the solution domain. To solve the diffusion equation inside the droplet analytically, we divide the droplet volume into two regions: an inner core, where the concentration profile remains uniform, and an outer shell, where the liquid phase exhibits a concentration gradient (see Fig. 4). The water concentration profile inside the droplet can therefore be calculated from the solution of the diffusion equation in the outer shell that does not include the singular point.

Additionally, we consider a two-stage evaporation model. In the first stage [Fig. 4(a)], the droplet evaporation reduces the droplet size and, at the same time, increases the thickness of the outer shell, where the water concentration has decreased due to loss of water by evaporation and the liquid phase is inhomogeneous. In this stage, both the droplet radius $R$ and the radius of the inner core $Ri$ shrink with time while the water concentration in the inner core $cwi$ remains equal to its initial value. When $Ri$ reaches a cut off value $Ri°$ (which is of molecular size), the drying process turns into the second stage [Fig. 4(b)]. In this stage, $Ri$ is fixed to $Ri°$ while $R$ and $cwi$ are considered as time-dependent parameters, which reach their equilibrium values at the end of the second drying stage.

We first assume that the water diffusion coefficient in the liquid phase $Dwl$ is independent of solute concentration (concentration dependence of $Dwl$ is discussed next). Using the two-stage model described above, the liquid water density profiles in the first drying stage $cw1st(r)$ and in the second drying stage $cw2nd(r)$ are obtained as (the calculations are detailed in supplementary material Sec. B)

$R$ and $Ri$ in the above equations are time-dependent parameters and $μ$ in Eq. (15) is a function of $R$⁠, which is given by

The numerical prefactor $α$ in Eqs. (14)–(16) is given by $α=Dwl/(Dwcgvw)$ and describes the ratio of internal and external water diffusivities. Considering the values listed in the Nomenclature, this prefactor is approximately equal to $α≈4$ at 25 °C. Using the assumption of infinitely rapid diffusion in the droplet, which is the main assumption in the previous section, $α$ goes to infinity. In this section, we aim to investigate the effect of the internal diffusion constant on the water evaporation process, and thus the calculations are done for different finite values of $α$⁠.

Neglecting the probability of crust formation due to the increased solute concentration at the droplet surface (which is discussed later), one can calculate the droplet evaporation time as

where $τev,s1st$ and $τev,s2nd$ are the times of the first and the second drying stages, respectively. As detailed in supplementary material Sec. B.4, mass conservation of water and solutes gives rise to the following differential equations that describe the time-dependent variation of the droplet radius in the two drying stages,

Equations (18) and (19) are solved using the numerical approaches explained in supplementary material Sec. B.4 to yield $t(R)$ in each drying stage and, consequently, $τev,s1st$ and $τev,s2nd$⁠. In our calculations, the radius of the inner core at which the droplet turns into the second drying stage $Ri°$ is set to $1 nm$⁠, which is small enough so that the calculated evaporation time is independent of $Ri°$ [see Fig. 5(a)]. Figure 5(b) shows the droplet evaporation time obtained from Eqs. (17)–(19) as a function of the liquid water diffusion coefficient $Dwl$⁠. This figure demonstrates that a decrease in $Dwl$ increases the droplet evaporation time. This is consistent with experimental evidence indicating that the extremely extended evaporation time of secondary organic aerosol (SOA) particles can be attributed to the reduced water diffusivity in the viscous phase.^33,79 For a high liquid water diffusion coefficient, however, $τev,s$ is found to be independent of $Dwl$⁠. Figure 5(b) also confirms that for the rather low initial solute volume fraction of $Φ0=0.01$⁠, Eq. (10) accurately predicts the evaporation time in the limit of $Dwl→∞$⁠, as discussed in Sec. II A 2.

As stated before, our diffusion model is valid only when we assume that the diffusion coefficients of water and solutes in the liquid phase are equal, $Dsl=Dwl$ (see supplementary material Sec. B.1 for more details). Using this assumption and neglecting the concentration dependence of the diffusion coefficient, which will be discussed next, the value reported in the Nomenclature for the water diffusion coefficient in pure water at 25 °C, $Dwl,25 °C=2.3×10−9 m2/s$⁠, which is indicated by the dotted vertical line in Fig. 5(b), can be used to estimate the evaporation time of a solute-containing water droplet at 25 °C. However, most relevant solutes diffuse slower than water molecules (i.e., $Dsl<Dwl$⁠), and thus the effective diffusion coefficient to be used inside the droplet is the solute diffusion coefficient, as discussed before. According to Fig. 5(a), therefore, one expects the evaporation process of droplets that contain slowly diffusing solutes to take longer than what is estimated here using the water diffusion coefficient in the calculations.

Figure 5(c) shows the evaporation times obtained from Eqs. (17)–(19) using the water diffusion coefficient at 25 °C, $Dwl=Dwl,25 °C,$ together with results obtained in the limit of $Dwl→∞$⁠, which leads to a homogeneous internal water concentration. According to this figure, although the effect of internal concentration gradients on the droplet evaporation time is negligible at low initial solute volume fractions, neglecting this effect can make a relatively large error at higher values of $Φ0$⁠, especially when $Φ0$ becomes close to $Φ0τmax$ [see Eq. (12)]. The relative error is found to be more significant in dry environments because the drier the air, the faster the water evaporates from the droplet surface and the larger the water concentration gradient becomes. This figure also shows that the presence of internal concentration gradients cannot significantly change the value of the solute volume fraction at which the droplet evaporation time is maximal, denoted as $Φ0τmax$⁠.

As mentioned before, the differential equations (18) and (19) that are used here to calculate the evaporation time are not analytically solvable and, hence, the evaporation time is calculated using numerical methods. To provide an equation for estimating $τev,s$ as a function of the relevant parameters that characterize the mechanisms discussed so far (i.e., solute-induced water vapor-pressure reduction, concentration dependence of the evaporation-cooling effect, and evaporation-induced concentration gradients inside the droplet), we construct a heuristic fit function inspired by Eq. (10),

where $aτ1$ and $aτ2$ are fit parameters and $θ′$ is a modified numerical prefactor given by

Figure 5(c) demonstrates that Eq. (20) (shown by broken lines) with fit parameters $aτ1=1.03$ and $aτ2=56$ perfectly fits the data over a wide range of $Φ0$ and $RH$⁠. It is worth mentioning that the values used for $aτ1$ and $aτ2$ are the averages of the fit parameters obtained for different values of $Φ0$⁠, $RH$⁠, and $γ$⁠.

Finally, the effect of an internal concentration gradient on the time-dependent variation of the droplet radius is investigated. The point to note here is that in the case of a droplet with a finite internal water diffusivity, the first drying stage is dominant. Indeed, when $Ri°$ is set to a sufficiently small value (on the order of a few nano-meters or smaller), the time of the second drying stage $τev,s2nd$ goes to zero [as demonstrated in Fig. 5(a)], meaning that the second drying stage becomes negligible. Therefore, the time-dependent radius of a droplet with finite $Dwl$ can be calculated from Eq. (18) only, which describes the first drying stage. Also, $R(t)$ in the limit of $Dwl→∞$ can be obtained from Eq. (3), which neglects the internal water concentration gradients. The results are shown in Fig. 5(d). This figure indicates that in the beginning of the evaporation process, when the evaporation-induced concentration gradients are not significant yet, the shrinkage rate of the droplet radius is almost independent of the internal water diffusion coefficient. As the evaporation proceeds, however, this parameter becomes important and affects the droplet evaporation time. This figure also shows that in the presence of slowly diffusing solutes (i.e., solutes with lower diffusion coefficients than water), the effect of the evaporation-induced internal water concentration gradient is more significant, as expected.

The diffusion coefficient of water is known⁴³ to decrease in the presence of strongly hydrated solutes, such as NaCl, NaI, and NaBr. In contrast, the presence of weakly hydrated solutes, such as KCl, KBr, and KI, increases the water self-diffusion coefficient, presumably by disrupting hydrogen bonding. According to Fig. 5(b), such variations in the water diffusion coefficient can slow down or speed up the droplet evaporation process. Therefore, the evaporation model introduced in Sec. II B must be modified in order to account for the solute-concentration dependence of water diffusivity inside the droplet. Here, we assume that the internal water diffusion coefficient is linearly dependent on the solute concentration, which approximates the behavior of NaCl and NaBr salts rather accurately.⁴³ We assume such a linear dependence also for other solutes such as proteins, carbohydrates, and viruses. Accordingly, we write

with $Dwsol$ and $Dwl$ being the water diffusion coefficients in an aqueous solution and in pure water, respectively, $β$ being a solute-specific constant (⁠$∼0.065 M−1$ for NaCl and $∼0.058 M−1$ for NaBr solutions⁴³), $vs$ being the molecular volume of solutes, and $csr,t=1/vs−vwcwr,t/vs$ being the solute concentration profile.

According to Eq. (22), when the solute concentration at the droplet surface reaches a critical value of $cs*=β−1$⁠, the water diffusion coefficient at the droplet surface $Dwsol(R)$ reaches zero. Therefore, in the cases where $cs*$ is lower than the equilibrium concentration of solutes $csev=(1−RH/γ)/NAvs$⁠, i.e., the solute concentration at which the droplet reaches the evaporation-equilibrium state [see Eq. (A14)], the water diffusion flux at the droplet surface $−Dwsol(R)∂cl∂rr=R$ goes to zero before the net evaporation flux $kecwR−kccv(R)$ vanishes. In such case, the reactive boundary condition given by Eq. (A3) cannot be satisfied, which prevents the numerical solution of the diffusion equation from converging. To avoid this problem, we exclusively consider the parameter range of $β$ and $RH$ so that $csev≤cs*$⁠. In this limit, the water diffusion coefficient remains non-zero during the entire drying process and minimally reaches zero in the evaporation-equilibrium state. At first, the possibility of phase separation at high solute concentrations is neglected. In reality, when the solute concentration at the droplet surface exceeds a critical supersaturation concentration $csmax,$⁸⁰ which is normally reached before $Dwsol$ goes to zero, a solid crust forms at the droplet surface. For example, $csmax$ in the case of aqueous NaCl droplets at $20 °C$ is around 11.1 M,⁸ and $β$ for a NaCl solution is $∼0.065 M−1.$⁴³ Therefore, according to Eq. (22), the water diffusion constant in a liquid droplet of NaCl solution at the onset of crust formation is $∼28%$ of the water diffusion constant in pure water. The possibility of crust formation is accounted for in Sec. II B 2.

By incorporating $Dwsol$ from Eq. (22) into the water-diffusion equation [Eq. (1)] and using the numerical methods described in supplementary material Sec. C to solve this equation along with the heat and mass conservation equations, the time-dependent radius of the evaporating droplet is determined. Panels (a)–(c) of Fig. 6 show the results in the absence of solutes and in the presence of solutes that increase or decrease the internal water diffusivity [respectively, corresponding to setting $β=0$⁠, $β<0$⁠, and $β>0$ in Eq. (22)]. This figure indicates that the solute-concentration dependence of the diffusion coefficient does not significantly change the time-dependent shrinkage of the droplet radius, even for droplet radii close to the final equilibrium radius, where the solute concentration at the droplet surface tends to its maximum value. The reason for this behavior is discussed below.

The presence of solutes inside a drying droplet has two effects on the internal water diffusivity: first, it changes the water diffusion coefficient in the droplet at the onset of evaporation, when the water concentration is homogeneous, and second, it induces a time-dependent diffusivity profile inside the droplet. Accordingly, Eq. (22) can be rewritten as

where $Dwsol,i$ is the initial water diffusion coefficient in the liquid droplet and thereby also the water diffusion coefficient within the inner core where the solute concentration equals its initial value (see Fig. 4), $csi$ is the initial solute concentration, and $ΔDwsolr,t$ indicates the time-dependent local variation of the water diffusion coefficient, which is due to the evaporation-induced concentration gradient in the outer shell (see Fig. 4). We first investigate the effect of the solute-concentration dependence of the initial water diffusivity on the droplet evaporation. According to Figs. 5(b) and 5(d), this effect can significantly change the droplet evaporation time only if the initial deviation of the internal water diffusivity from its value in pure water is sufficiently large (⁠$Dwsol,iDwl<0.5$ or $Dwsol,iDwl>10$⁠). For the parameter range studied here, however, this deviation is rather small [see Fig. 6, panels (d)–(f)], and thus the evaporation process is not affected by the solute-concentration dependence of $Dwsol,i$ [see Fig. 6, panels (a)–(c) and (j)–(l)]. This is expected to be the case for most droplets of ionic solutions, where $β$ is typically small, and for most respiratory droplets, where the initial solute concentration $csi$ is not large enough to significantly affect $Dwsol,i$⁠. As evaporation proceeds, the increased solute concentration at the droplet surface gradually changes the water diffusion coefficient at the surface from its initial value and thus the second effect, i.e., the effect of an inhomogeneous diffusivity profile, becomes important. Figure 6, however, shows that this effect does not significantly change the shrinkage rate of the drying droplet, even in the cases with $β>0$⁠, where the water diffusion coefficient at the droplet surface decreases down to zero during the evaporation process (this effect only makes a slight change in the droplet shrinkage rate at the very end of the evaporation process, when $R≃Rev$⁠). The reason for this is examined next.

As follows from Eq. (22), an increase in $β$ decreases the internal water diffusion coefficient at the droplet surface $Dwsol(R)$ [see Fig. 6, panels (d)–(f)]. On the other hand, an increase in $β$ slows down the internal mixing process and, consequently, reduces the growth rate of the outer shell, i.e., the shell within which the liquid phase exhibits a concentration gradient [see Fig. 6, panels (g)–(i)], while the solute concentration at the droplet surface remains almost independent of $β$ [see Fig. 6, insets of panels (d)–(f)]. This means that the water concentration gradient at the droplet surface $cw′R=−∂cw∂rr=R$ increases with an increase in $β$ [see Fig. 6, insets of panels (g)–(i)]. As a result, the water flux density at the liquid–vapor interface $j$⁠, which is the product of $Dwsol(R)$ and $cw′R$⁠, is found to remain almost independent of $β$ [see Fig. 6, panels (j)–(l)]. This implies that the solute-induced local variation of the internal water diffusivity has only a minor effect on the water evaporation rate and, consequently, on the time-dependent droplet radius. As a result, the term $ΔDwsolr,t$ in Eq. (23) can be safely neglected in the calculation of the droplet evaporation time and, thus, one can replace $Dwl$ in the numerical calculations described in the previous section (and in supplementary material Sec. B) by $Dwsol,i$ to estimate the droplet evaporation time in the presence of solutes and in the absence of the possibility of phase separation and crust formation. Also, Eq. (20), which neglects the concentration-dependence of the water diffusion, can be safely used to estimate the droplet evaporation time in case when $βcsi$ is small enough to not significantly change the initial water diffusion coefficient from its value in pure water, i.e., for $0.5≤Dwsol,iDwl≤10$⁠. For droplets that initially contain a relatively large solute fraction, especially solutes that dramatically decrease the water diffusion coefficient, such as saccharides,⁸¹ the significantly reduced initial water diffusivity inside the droplet is expected to slow down the evaporation process [see Fig. 5(b)] and, thus, Eq. (20) underestimates the droplet evaporation time.

So far, we examined various effects of nonvolatile solutes on the droplet evaporation process, without considering the possibility of crust formation (which is discussed next). Numerical examples for evaporation times obtained from the relevant equations are given in Table I for a droplet consisting of a NaCl solution with initial radius $R0=50 μm$ at 25 °C. According to this table, the presence of solutes on the one hand decreases the droplet evaporation time by producing a lower limit for the droplet size that can be reached by evaporation and by reducing evaporation cooling of the droplet. On the other hand, the evaporation time is found to increase due to the solute-induced water–vapor pressure reduction and due to the presence of local concentration gradients inside the droplet that arise from the finite diffusivity of solutes and water in the liquid phase. Overall, the latter two effects are dominant and thus the droplet evaporation time increases in the presence of solutes. Neglecting solute effects, therefore, causes an underestimate of the droplet evaporation time. Table I also shows that solute effects are more pronounced at high initial solute volume fractions.