Modeling the dynamic behavior of a droplet evaporation device for the delivery of isotopically calibrated low-humidity water vapor

A simple model is presented that gives a quantitative description of the dynamic behavior in terms of water concentration (humidity) and isotope ratios of a low-humidity water vapor generator. The generator is based on the evaporation of a nL-droplet produced at the end of syringe needle by balancing the inlet water flow and the evaporation of water from the droplet surface into a dry air stream. The humidity level is adjusted by changing the speed of the high-precision syringe pump and, if needed, the dry air flow. The generator was developed specifically for use with laser-based water isotope analyzers in 5 Antarctica, and recently described in Leroy-Dos Santos et al. (2020). Apart from operating parameters such as temperature, pressure, water and dry air flows, the model has as "free" input parameters the water isotope fractionation factors and the evaporation rate. We show that the experimental data constrain these parameters to physically realistic values that are in reasonable to good agreement with literature values where available.

It is no overstatement to say that laser-based isotope analyzers have revolutionized the field of water isotope ratio instrumentation, until not so long ago dominated by Isotope Ratio Mass Spectrometers (e.g., Kerstel, 2004;Kerstel and Gianfrani, 2008). In particular, laser instruments have enabled continuous measurements of low-humidity atmospheric air in airborne and Antarctic field settings (see, among others, Iannone et al., 2009bIannone et al., , 2010Moyer et al., 2013;Steen-Larsen et al., 2013;Casado et al., 2016;Ritter et al., 2016;Bréant et al., 2019). In order to calibrate such instruments against international standard and reference mate-15 rials that are all in liquid form, it is necessary to bring these into the vapor phase without causing fractionation, or alternatively with well-controlled, quantitative fractionation, while also controlling the level of humidity (the volume mixing ratio). Several solutions have been proposed and developed into prototypes and commercial instruments, but few are capable of delivering a stable supply at low humidity levels (Iannone et al., 2009a;Sturm and Knohl, 2010;Gkinis et al., 2010;Tremoy et al., 2011).
One approach is that of the instrument developed in our laboratory with the specific aim of calibrating laser-based analyzers 20 deployed in Antarctica and first reported in Landsberg et al. (2014). This prototype instrument has undergone significant engineering developments in order to improve its performance and robustness (Leroy-Dos Santos et al., 2020). In this companion paper we provide a theoretical model that is capable of a quantitative description of its operation. reservoirs with volume V , water fluxes Φ, and isotope ratios R involved in the model. Since r 2 = r 2 0 + r 2 d and r d = r − h, it follows that r = r 2 0 + h 2 / (2h). (b) The isotope ratio profile over the liquid to vapor boundary (left to right, with the thin vertical lines representing the growing water surface) at three instants in time if the water flux Φ0 from the syringe (with δ0 = 0) follows a step function with Φ0(t) = 0 for t < t0, and Φ0(t) = F > 0 for t > t0. The isotope fractionation is taken to be ef f ≈ −71‰ for δ 2 H. While the droplet is growing δe < δ0.
At t = t∞ the incoming water flux Φ0 equals the evaporated water flux Φe and δe = δ0.
For completeness, we assume that only a fraction f of the droplet volume (a boundary layer) becomes enriched. This is later used to demonstrate that the best model results are obtained by assuming that the entire droplet becomes enriched (f = 1, see 45 section 3.1), an observation that is further supported by considerations involving the relative speeds of isotopic diffusion and the water flow (section 4.1). Figure 1(b) shows the radial isotope concentration profile inside the droplet and the neighboring vapor following a step function in the flow rate from zero to some fixed value at t = t 0 . The actual form of the profile is not important for our model and could just as well be approximated by a square profile. We can thus distinguish four different bodies of water: 50 1. The syringe reservoir with a constant isotope ratio R 0 and an outgoing water flux equal to Φ 0 (t) determined by the syringe pump speed, 2. The core volume of the droplet with an isotope ratio R c = R 0 and a time-dependent volume V c (t). The water flux from the core to the surface layer of the droplet is given by Φ c (t). Only in steady state Φ c = Φ 0 .
3. A fraction f (0 < f ≤ 1) of the total droplet volume V d (t) that will become enriched in the heavy isotope, V s (t), with 55 isotope ratio R s (t): 4. The evaporated water flux Φ e (t) leaving the droplet with isotope ratio R e (t) = R s (t) · α ef f . Here the relevant isotope fractionation factor is that between the vapor and liquid phase water: α ef f = (1 + ef f ) < 1.
the droplet: Figure 1 gives a schematic representation of our model, indicating the relevant water volumes and inter-volume fluxes, as well as the isotope ratios R of each volume. Solving the model ab initio is not difficult and will be shown to give a qualitatively and quantitatively satisfactory description of the dynamics under realistic conditions.

65
The free input parameters to the model are (1) the effective liquid-to-vapor fractionation factor α ef f , the estimate for which we base on studies by Cappa et al. (2003) and Luz et al. (2009), (2) the fraction f of the droplet that becomes enriched, which we will show in section 3 to be equal to unity, and (3) the evaporation coefficient k e , which in turn we base on previous work by Walton (2004) and Sefiane et al. (2009), and justify in the discussion of section 4.4.
The first task is to model the evaporated total water flux Φ e (t) as a function of a variable input water flux Φ 0 (t), driven by The evaporation flux Φ e (t) is a function of the droplet size through (2). For simplicity, we model the droplet at the tip of the needle as a partial sphere, a spherical cap. The surface area of the spherical-cap-shaped droplet is given by (see Figure 1): with the radius of curvature of the cap r, the cap height h (0 ≤ h ≤ 2r), and the inner diameter 2r 0 of the needle, all as defined in Figure 1. The volume of the droplet is equally a function of h: A s (t), and thus Φ e (t), can then be expressed in terms of V d by inversion of (5), with h(V d ) being obtained as the only real 80 root of the cubic equation, giving: where α := 3 108v + 12 12u 3 + 81v 2 and 85 u := 3r 2 0 , v := The above already permits calculating the droplet size (e.g., in terms of the droplet radius r(t) = r 2 0 + h(t) 2 / (2h)) and the evaporative water flux Φ e (t), both as a function of the time-dependent input water flux Φ 0 (t), by numerical integration of Eq. (3).
Going one step further we include a second mass balance equation in the model to account for the rare isotopologues (in 90 practice either 2 H 16 O 1 H or 1 H 18 O 1 H). We start with expressing the rare isotope fluxes (identified by ) in terms of the total fluxes and the isotope ratio of the reservoir in question. For the three relevant fluxes (see Figure 1): Here R w is the ratio of the abundance of the rare to the most abundant water isotope in the reservoir w (w = 0, c, s, e for respectively, the syringe and needle, the core of the droplet, droplet surface layer, and the evaporated water) (e.g., 2 R s = [ 2 H]/[ 1 H ]) s ). R V SM OW is the isotope ratio of the international standard material Vienna Standard Mean Ocean Water (IAEA, 2017), and δ w := (R w − R V SM OW ) /R V SM OW , the relative deviation of the abundance ratio in reservoir w with respect to that of VSMOW. Finally, the fractionation factor between the (evaporated) vapor phase water and the liquid α ef f ≈ 1 100 ( ef f 1), making the approximation made in Eq. (11) a very good one.
Similar equations as (9), (10), and (11) hold for the different water reservoir volumes, allowing us to write for the volume of the isotopically enriched evaporating surface layer: Substitution of: And using the definition: where we have used the approximation for Φ e (t) of Eq. (11).
We can now calculate the isotope ratio in the enriched fraction f of the droplet volume (using Eq. (1) and an appropriate value of f ) by integrating Eq. (15), while evaluating Eq.
(1) to (4) at each time step. The isotope ratio of the evaporated water 115 is then obtained as: with: and: 3 Results The above model has been programmed in Mathcad (PTC Mathcad, 2020) and used to simulate data that were recorded with a home-built, low-humidity water isotope spectrometer, named HiFi, described in Landsberg (2014) and Landsberg et al. (2014). Since we are specifically interested in the dynamic behavior of the water vapor source that feeds the spectrometer, 125 we need to take the response time of the spectrometer into account. This response is typically described by a double or even triple exponential. At humidity levels of several thousand ppmv (parts per million by volume) the initial (fast) response time of the bare spectrometer was determined to be in the range of 1 to 2 s for both the water concentration and the isotope ratios, with a second, slower exponential response of the order of 15 s. However, in the configuration of this study, and at the low water concentrations of a few hundred ppmv, the response time is significantly longer, especially for the δ 2 H 130 isotope ratio. These response times were measured using as a humidity source a predecessor of the isotopic humidity generator described in Leroy-Dos Santos et al. (2020), which was, however, equipped with two independent syringe pumps, enabling rapid switching between two different water sources using a 2-position, 4-port valve (Valco EUDA-4UWE) just before the spectrometer (Landsberg, 2014). The humidified air stream was sent either to the spectrometer or to a waste pump. The isotope response was determined by switching between two very different water standards, assuring a high signal-to-noise-ratio of 135 the measurements, while keeping the concentration constant at about 600 ppmv. The standard waters used were left-over  The normalized response curves of the spectrometer for switching between the GS-48 and BEW-2 isotope standards, both prepared as a mixture of ∼600 ppmv water vapor in dry air. The experimental data is fit with a double exponential yielding for the fast decay times 9.2 s and 20.7 s for δ 18 O (red curve) and δ 2 H (green curve), respectively.
The instrument isotope response curves are shown in Figure 2, while the double exponential fit parameters are summarized in

Humidity and isotope step responses
The model detailed in section 2 was first used to simulate the dynamic behavior of the combination of the LHLG and the HiFI isotope analyzer, while the LHLG was programmed to generate small humidity steps of about 200 ppmv around an absolute value of roughly 400 ppmv. The simulated water concentration response was fit to the experimental data by adjusting the 150 evaporation rate k e , the only free parameter in this case (see the top panel of Fig. 3). We will discuss the rationale for the values of k e determined in our study later in section 4.4. Having fixed the evaporation rate at an optimal value of k e = 3 µm/s, the next step is to confirm that the isotope responses are modeled correctly, taking into account that both the δ 2 H and δ 18 O simulated responses also depend on the fraction f of the droplet volume that becomes enriched, as well as the effective liquid-to-vapor fractionation factor α ef f . Since we expect the entire droplet to become enriched in the heavy isotopologues, we start with the 155 assumption that f = 1. As we will see shortly, this choice is validated by the experimental observations. It will subsequently be rationalized by theoretical considerations in section 4.1. The best fit is obtained for ke ≈ 3 µm/s, whereas a higher (lower) value results in a simulated dynamic response that is too fast (slow) compared to the measured response.
As to the fractionation factors, we consider that at the very low relative humidity of our experiment (h ≈ 0.01), the effective fractionation factors α ef f can be written as the product of a diffusion fractionation factor α dif f and an equilibrium fractionation factor α eq (Cappa et al., 2003). Moreover, the diffusion fractionation factor can be related to the ratio of the molecular  (Stewart, 1975), such that we may write: Here the label x refers to the rare isotope or isotopologue (  the case of fully turbulent flow. The equilibrium fractionation factors were accurately determined by Horita and Wesolowski (1994), and we take their values at 35°C. The diffusivities were determined by Cappa et al. (2003) and more recently by Luz et al. (2009). We will use the more recent values, but the difference is minimal for our purpose (Cappa et al. (2003) predict only slightly lower values of α ef f in the laminar limit of n = 1). Table 2 gives the values of the effective liquid-to-vapor fractionation factors for three different values of the flow parameter n and Fig. 4 shows the corresponding model simulations 170 compared to experimental data. The δ 18 O simulation shows a relatively large effect of changing n than the δ 2 H simulation. In contrast, changing the values of f has the same relative effect on both simulations (not shown in 4). With n = 0.43 and f = 1 a good fit to both isotope response curves is obtained. We thus also conclude that our data supports the theoretical finding (section 4.1) that f = 1.

175
The LHLG was modified immediately following the experiments presented in the previous section. Notably, it was deemed that the bore of the aluminum injector chamber that accepts the syringe's needle was too narrow. With an internal diameter of only 2 mm, careful guiding of the needle, and consequently a precise positioning of the syringe, was needed to avoid occasional contact of the droplet with the chamber wall. This also limited the maximum droplet size, and therewith the volume mixing ratio (humidity level) that could be attained to roughly 1000 ppmv. The injection chamber was therefore replaced by 180 a stainless steel sample cylinder with volume 75 mL and Sulfinert hydrophobic coating (Restek 304L-HDF4-75). Because the flow velocity is now significantly lower, the coating serves to minimize the memory effect due to surface adsorption of water molecules. In addition, a section of PTFE tubing was added between the syringe (Hamilton 84853) and the removable needle to make the alignment better manageable. This gave initially rise to unexpected results that were attributed to the appearance of small air bubbles in the water injection line. These problems were later resolved by re-engineering the LHLG as described 185 in Leroy-Dos Santos et al. (2020). But we report on these "useless" results here because they nicely demonstrate the ability of the model to simulate the behavior of this non-ideal instrument, and thus validate the model under a different operating regime.
During similar experiments as reported in the previous section 3.1, recording the response of the LHLG following small steps in the flow of injected water, relatively large sinusoidal oscillations were observed with a period that matched the revolution speed of the lead screw of the precision pump. We submit that these oscillations become prominently visible when small 190 imperfections of the lead screw combine with small air bubbles present in the water injection line, possibly amplified by viscous resistance of the liquid inside the water line and needle. Whatever the precise underlying mechanics, we modeled a sinusoidal variation of the water flow with a period equal to one revolution of the screw drive. The amplitude and phase of the (possibly amplified) lead screw imperfection was chosen to yield a simulation that best matched the observed amplitude of the oscillations. The only other parameter that needed adjustment was the evaporation rate. The value of k e = 1 µm/s was found 195 to produce a simulation that best matched the water concentration response when the pump was switched between different water flow rates, as seen in the upper panel of Fig. 5. The lower evaporation is due to the lower flow velocity of the air around the droplet (see section 4.4). The corresponding response of the isotope ratios is shown in the lower panel of Fig. 5. It may be clear that the correspondence between simulation and experiment is (already) satisfactory, considering that no further parameter adjustments were made. We

Droplet isotopic enrichment
Here we provide support for the observation of an enrichment in the heavy isotopologues of the entire droplet, and not just in a surface layer of limited thickness. Referring to Fig. 4 (for which n = 0.43, i.e., 18 α ef f = 0.98, and f = 1), in principle the 205 same amplitude of the modeled response can be obtained by assuming fully laminar flow (n = 1), ánd assuming that a much smaller fraction of the droplet becomes enriched in the heavy isotopes. This gives, however, a less satisfactory fit to the data, as shown in Fig. 6. Notably, the response simulated with n = 1 (i.e., 18 α ef f = 0.9650) and f = 0.5, reached the same maximum amplitude, but is clearly narrower than the experimental curve. Importantly, this is also not what is predicted based on the speed of isotopic diffusion in the droplet. To see that in our experiment f should be equal to unity we consider that the enrichment occurring at the surface of the droplet will diffuse inwards, resulting in an isotope gradient inside the droplet with a characteristic diffusion length given by Bird et al. (2006): with D the diffusion coefficient and t time. Differentiation of (20) yields the velocity of the diffusion front: The diffusion coefficients of HDO and H 18 OH in water have been measured by Horita and Cole to be 2.34 10 -9 m 2 /s and 2.66 10 -9 m 2 /s, respectively (Horita and Cole (2004)). This shows that diffusion over lengths comparable to the size of a typical droplet (0.1 mm) takes place on a time scale of the order of 1 s. It is therefore likely that the entire droplet becomes isotopically enriched, rather than just a surface layer: f = 1.

Back Diffusion
The question arises whether the diffusion is strong enough to allow the isotopic enrichment to propagate all the way to the syringe reservoir. To answer this question we compare the diffusion velocity of (21) to the flow velocity inside the syringe needle: Even at the lowest water flow rates of about 100 pL/min, the diffusion can be stopped well within the length of the needle (if necessary by reducing the needle inner diameter). We conclude that it is unlikely that the isotopic composition of the syringe 235 reservoir would change due to back-diffusion of heavier isotopologues. This was also confirmed experimentally by bringing the same liquid standard material into the vapor phase with both the LHLG and a commercial humidity generator (Picarro SDM) at time intervals of one month and not observing any difference between the measurements (within the measurement precision of 0.2‰and 1‰for δ 18 O and δ 2 H, respectively) (Leroy-Dos Santos et al., 2020).

Fractionation factors 240
The effect of the precise values of the 2 H 16 O 1 H-and H 18 2 O-isotopologue effective fractionation factors on the simulations was already discussed to some extent in section 3.1, where it was found that the best match with experiment is obtained by assuming fractionation factors that correspond to an intermediate case between laminar and turbulent flow (characterized by n = 0.43). This can be rationalized by estimating the Reynolds number for the flow around the water droplet, R e = ρvL/µ.
In the previous formula ρ ≈ 1.25 kg m −3 is the density of the air flowing around the needle and droplet; v = 1.6 m/s is the 245 velocity of the air around the droplet, inside the narrow-bore chamber (inner diameter 2 mm), given the air flow of 300 mL/min (STP); L ≈ 0.5 mm is the diameter of the droplet; and µ = 18.3 µPa·s is the dynamic viscosity of air at 35°C. We thus find R e ≈ 60. This contrasts with a value of v ≈ 0.007 m/s and R e ≈ 0.2 for the case of the about 30-mm internal diameter steel cylinder used in the modified instrument. The latter case should thus be much closer to the limit of fully laminar flow. We are shown in Fig. 7. Whereas the differences for 2 H are minor, the effect of the larger 18 O-fractionation (i.e., smaller liquid-to-vapor fractionation factor) in the laminar flow regime is clearly visible, and arguably provides a slightly better fit to the data, primarily during the water concentration changes, as can be seen in Fig. 7. It should be noted, however, that in the regions of oscillatory behavior in between the concentration steps, the fit could also have been nudged by adjusting the amplitude of the lead screw modulation.

255
Still, we conclude that the results of section 3.2 are just as well, and most likely better, described (than shown in Fig. 5) by assuming fully laminar flow.

Evaporation rate
The two experiments that we discussed here in sections 3.1 and 3.2 required rather different evaporation rates to simulate the data with our model, k e ≈ 3 µm/s and 1 µm/s, respectively. The difference is clearly related to the different Reynolds 260 numbers or, more directly, the different dry air flow velocities of 1.6 m/s and 0.007 m/s. In fact, the values are in reasonable agreement with the results reported by Walton (2004). Although his measurements were recorded at only a small number of air temperatures and flow velocities, we can estimate values applicable to our situation by linear extrapolation of the observed rates as a function of flow velocity, and fitting a (weakly) quadratic dependence on the temperature to the data collected at a fixed flow velocity of 1 m/s. In Fig. 8 we present selected data of Walton together with the estimated values for our case. We 265 thus predict a rate of 5.2 µm/s at a flow velocity v = 1.6 m/s, and of 1.3 µm/s at v = 0 m/s, higher than the values we found experimentally. So far we have assumed that the droplet is at the same temperature as the evaporation chamber, but it cannot be excluded that the actual droplet temperature is lower, especially in the high velocity case. However, the study by Sefiane et al. (2009) measured an evaporation rate at 22°C and 1 bar corresponding to ∼ 4 µm/s, very close to the value we extrapolated from the data of Walton (2004) at 25°C. It should also be mentioned that it is unlikely that the difference with the observations 270 by Walton or Sefiane are due to an under-estimation of the spectrometer humidity response time, as it is difficult to imagine a water concentration time response slower than that of the isotope ratios.

Conclusions
We have shown that the dynamic behavior of a humidity generator based on droplet evaporation can be accurately modeled.
Confrontation with experimental data of the water concentration and two isotopic ratios as a function of the injected water 275 flow, enables the determination of physically realistic values of the droplet evaporation rate and the liquid-to-vapor isotope fractionation factors. However, the signal-to-noise ratio of the analyzer at the very low humidity levels investigated is not quite sufficient to make very precise determinations of the fractionation factors. But recent developments in ultra-precise and ultra-sensitive isotope measurements (e.g., Stoltmann, 2017;Kassi et al., 2018) will enable to deliver more precise values by at least an order of magnitude. What may appear as a bit of a quixotic study of evaporating water droplets, may thus in 280 fact permit measuring not only the evaporation rate, but also the effective fractionation factors, and therewith also isotopologue dependent diffusivity ratios, in the evaporation of small sessile droplets. Apart from this potentially new application, it is highly satisfactory to be able to accurately simulate the dynamic behavior of the LHLG with few free parameters, and under rather different operating conditions.