First‐Principles Diffusivity Ratios for Kinetic Isotope Fractionation of Water in Air

Kinetic isotope fractionation between water vapor and liquid water or ice depends on the ratio of the diffusivities of the isotopic species in air, but there is disagreement as to the values of these ratios and limited information about their temperature dependence. We use state‐of‐the‐art intermolecular potential‐energy surfaces for the water‐nitrogen and water‐oxygen pairs, along with the kinetic theory of molecular gases, to calculate from first principles the diffusivities of water isotopologues in air. The method has sufficient precision to produce accurate diffusivity ratios. For the HDO/H2O ratio, we find that the often used hard‐sphere kinetic theory is significantly in error and confirm the 1978 experimental result of Merlivat. For the ratios involving 17O and 18O, the simple kinetic theory is relatively close to our more rigorous results. We provide diffusivity ratios from 190 K to 500 K, greatly expanding the range of temperatures for which these ratios are available.


Introduction
Stable water isotopes, in particular the molecules HDO, H 2 17 O, and H 2 18 O, are widely used to model processes involving the atmosphere, ocean and fresh water, and ice (Gat, 1996). In many situations, isotopic fractionation between the atmosphere and a condensed phase is determined not only by equilibrium thermodynamics but also by a kinetic effect that depends on the relative diffusivities of the isotopic species in air. Evaporation and precipitation in environments where diffusion affects fractionation are described by models that incorporate both equilibrium and kinetic effects depending on the details of the process (Casado et al., 2016;Craig & Gordon, 1965;Gonfiantini et al., 2018;Horita et al., 2008;Jouzel & Merlivat, 1984;Lamb et al., 2017;Nelson, 2011). While the equilibrium fractionation is fairly well understood, at least for vapor-liquid equilibria (Horita et al., 2008;Japas et al., 1995), there is significant disagreement, especially for D/H fractionation, regarding the correct diffusivity ratio for the kinetic effect. It is the purpose of this paper to resolve these disagreements.
Since there seems to be no standard notation for these diffusivity ratios, for this work we define the relative diffusivities D r; HDO ≡ D HDO =D H2O , D r; 17 ≡ D H2 17 O =D H2O , and D r; 18 ≡ D H2 18 O =D H2O , where D i is the diffusivity of species i in air (or in a different gas; nitrogen is sometimes used as a proxy for air). We note that, in some of the literature, the reciprocals of these ratios are used instead.
D r,HDO and D r,18 were reported in air at 20°C by Ehhalt and Knott (1965), but no information about the experimental method was given. The first well-described experiments were those of Merlivat (1978), who reported D r,HDO and D r,18 at 21°C in nitrogen. Cappa et al. (2003) is often cited for diffusivity ratios, but the values in that paper were not obtained from experiment but rather were calculated from the simple ©2020. The Authors.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

RESEARCH LETTER
10.1029/2020GL089999 kinetic theory described below (and then shown to be consistent with their experiments in the context of other modeling assumptions). Barkan and Luz (2007) reported D r,17 and D r,18 in air at 25°C and 40°C; the same group (Luz et al., 2009) subsequently reported D r,HDO and D r,18 at four temperatures from 10°C to 69.5°C. To illustrate the lack of consistency among reported results, values for D r,HDO near 20°C are given as 0.9852 (Ehhalt & Knott, 1965), 0.9757 (Merlivat, 1978), 0.9839 (Cappa et al., 2003), and 0.9775 (Luz et al., 2009). Since it is the difference between D r and unity that affects the fractionation, these differences are significant.
The uncertain temperature dependence is also a problem, particularly for D r,HDO . The only temperature-dependent experimental study is that of Luz et al. (2009), who found values of D r,HDO that increase strongly with temperature. Other work, such as that of Cappa et al. (2003), has recommended a single value (in their case obtained from simple kinetic theory) that is constant with temperature. This situation creates uncertainty when performing calculations for ice and snow, which require diffusivity ratios extrapolated to temperatures far below the lowest measured temperature of 10°C; the difference between extrapolating the temperature dependence of Luz et al. (2009) and assuming a constant value can greatly affect the calculated kinetic fractionation.
The relationship between the fractionation of different isotopes is also of interest. The ratio of the D/H diffusive fractionation to that of 18 O is described by the quantity φ, which in our notation can be written as There is a wide variation of reported φ in the literature. The diffusive fractionation of 17 O is traditionally defined relative to that of 18 O, using a logarithmic ratio θ diff ¼ lnD r;17 lnD r;18 : Attempts to interpret data have relied on a simplified kinetic theory of gases, derived for mixtures of hard spheres at low density. In the first-order approximation, the ratio of the diffusivity of an isotopic species (subscript i) to that of the reference species (subscript 0) in gas G is (Chapman & Cowling, 1970;Merlivat, 1978) where Γ is the diameter of a molecule and M is its molar mass. Under the reasonable assumption that different isotopologues have the same diameter, the first factor is unity and Equation 3 reduces to a simple function of the molar masses. With M air ¼ 28:96546 g mol −1 (Picard et al., 2008), Equation 3 yields 0.9836 for D r,HDO and D r,17 and 0.9687 for D r,18 . The value of φ given by the simple kinetic theory is 0.525, and the value of θ diff is 0.5183.
Deviations of experimental results from these kinetic-theory values have led authors to discuss whether the "diameter" of a water molecule varies with isotopic substitution (Barkan & Luz, 2007;Cappa et al., 2003;Horita et al., 2008;Merlivat, 1978). The simplifying assumptions of Equation 3 have largely gone unquestioned (an exception is Luz et al. (2009), who noted the possible inapplicability of simple kinetic theory for polar gases).
Physically, the water molecule is very far from being a hard sphere, so one would not expect Equation 3 to work well. The D/H substitution might be particularly poorly described, because the mass asymmetry of the HDO molecule significantly changes the rotational dynamics and those dynamics are completely absent from the hard-sphere theory.
Modern kinetic theory can do much better. For molecular gases that can be modeled as rigid, the relevant collision integrals (which are often referred to as generalized cross sections) can, with sufficient computer time, be calculated essentially exactly from the full intermolecular potential-energy surface. In this work, we use state-of-the-art pair potentials for H 2 O-N 2 and H 2 O-O 2 interactions to calculate the diffusivity ratios D r,HDO , D r,17 , and D r,18 . We perform these calculations as a function of temperature, providing diffusivity ratios at conditions where they have not been measured.

Methods and Results
The present kinetic theory calculations are a direct extension of those performed recently by one of us for H 2 O in N 2 (Hellmann, 2019a) and H 2 O in O 2 and in air (Hellmann, 2020), which are based on new and highly accurate H 2 O-N 2 and H 2 O-O 2 pair potentials developed using state-of-the-art quantum-chemical ab initio approaches. We used these pair potentials without modification, thus assuming that isotopic substitution in the water molecule affects the collision dynamics mainly through the changes to the total molecular mass and to the moment of inertia tensor. In the supporting information (Text S1 and the associated Table S1), we provide an analysis, based on calculations with the CFOUR (Stanton et al., 2019) and ORCA (Neese, 2012) quantum-chemistry packages, that shows that the error introduced by using the unmodified pair potentials likely does not exceed 0.1% for D r,HDO and should be completely negligible for D r,17 and D r,18 . This approach has the important advantage that any inaccuracy in the intermolecular potential that would cause D H2O to be in error would have a similar effect on D for the substituted isotopologues, making the diffusivity ratios insensitive to such errors.
Here, we provide a very brief summary of the methodology of the kinetic theory calculations and refer the interested reader to Hellmann (2019aHellmann ( , 2020 for a more detailed description as well as for details on the intermolecular potentials. The relevant generalized cross sections (or collision integrals) for calculating the diffusivity ratios were extracted from classical trajectories describing binary collisions of HDO, H 2 17 O, and H 2 18 O with N 2 and O 2 . The trajectories were calculated assuming rigid molecules by solving Hamilton's equations numerically from precollisional to postcollisional asymptotic conditions. Generalized cross sections at a constant collision energy can be expressed for these molecules as 11-dimensional integrals over the initial states of the trajectories, which necessitated a Monte Carlo integration approach involving the calculation of typically a few million trajectories for each collision energy. The generalized cross sections as a function of temperature, from which the diffusivities in N 2 and O 2 can be directly computed, were obtained from those at constant collision energies by an appropriate thermal averaging procedure. The range of investigated collision energies allowed us to obtain the generalized cross sections, and thus the diffusivities, at temperatures from 190 K to 2000 K. The calculations of the energy-and temperature-dependent generalized cross sections were carried out using an in-house version of the TRAJECT code (Heck & Dickinson, 1996), which, unlike the original code, is not restricted to linear molecules.
The diffusivities were computed for water mole fractions x w in the limit x w → 0, which is the most sensible choice for atmospheric applications. In this limit, the diffusivities depend only on the unlike-species interaction potentials. Therefore, we did not need any models for H 2 O-H 2 O, N 2 -N 2 , and O 2 -O 2 interactions in this study. We note that the variation of the diffusivities with x w does not exceed a few tenths of a percent at any given temperature, an effect that should almost completely vanish when taking the diffusivity ratios.
While the kinetic-theory calculations provide the product ρ m D in the low-density limit (ρ m is the molar density), pressures in the atmosphere are low enough that this product will not differ significantly from its low-density value. Also, any small finite-pressure effects that might exist will cancel to first order when diffusivity ratios are taken.
To obtain the diffusivities in air, we weighted the diffusivities in N 2 and O 2 using the appropriate first-order kinetic theory relation (Marrero & Mason, 1972), where x N2 and x O2 are the respective mole fractions in dry air, with the value for N 2 also accounting for Ar (and all other minor components). This is justified because the diffusivities of water in argon and nitrogen are very similar (O'Connell et al., 1969). The mole fractions used in Equation 4 are x N2 ¼ 0:790603 and x O2 ¼ 0:209397 (Hellmann, 2020).
The calculated diffusivity ratios D r,HDO , D r,17 , and D r,18 in air are listed for selected temperatures up to 500 K in Table 1. They have expanded statistical uncertainties (related to the Monte Carlo integration over the initial conditions of the trajectories) of less than 0.05% (k ¼ 2, roughly corresponding to a 95% confidence interval). The expanded uncertainties listed in the table are the total expanded uncertainties, which also take into account another potential error source, namely, the neglect of quantum effects on the generalized cross sections. Quantum effects depend on the masses and moments of inertia of the molecules in addition to the pair potential and temperature, which is why they will not fully cancel out in the ratios. Note that our estimate for the influence of quantum effects is an educated guess based on experience, which should be quite conservative. It is supported, for example, by the fact that the viscosity of dilute water vapor calculated from classical generalized cross sections for H 2 O-H 2 O collisions, for which the neglect of quantum effects is expected to be more severe than for H 2 O-N 2 and H 2 O-O 2 collisions, differs from the best experimental data at and above ambient temperature by less than 1% (Hellmann & Vogel, 2015).

Comparison With Literature Data
Figure 1 compares our calculated results with those from the literature for D r,HDO (a) and D r,18 (b). The shading depicts our expanded uncertainty at the 95% confidence level. For D r,HDO , the difference from the simple kinetic theory, and from the datum of Ehhalt and Knott (1965), is quite large. On the other hand, we are in excellent agreement with the datum of Merlivat (1978). This excellent agreement remains if we adjust for the fact that Merlivat's experiments were in nitrogen instead of air; the value we calculate for D r,HDO in N 2 (see supporting information) differs from that in air by only about 0.0005. The data of Luz et al. (2009) are also in agreement near room temperature, but their temperature dependence, while having the correct sign, is much too strong. We note that the error bars plotted for experimental sources are those reported in the original papers, which probably are not complete uncertainty estimates in the modern sense. For example, those of Barkan and Luz (2007) and Luz et al. (2009) are described as the "precision" of their experiments, suggesting that possible systematic uncertainties are not included.
Regarding the work of Merlivat (1978), we note that her measured value for the absolute diffusivity of H 2 O in N 2 (which has a stated uncertainty of 1.6%) differs by only −0.9% from the calculated value (Hellmann, 2019a). The same approach used for H 2 O in N 2 and O 2 yields similar levels of agreement between theory and experiment for other gas pairs; for example, the best experimental data for the diffusivity of CO 2 in N 2 (with a stated uncertainty of less than 0.3%) agree within 0.2% with the first-principles results of Crusius et al. (2018).
For D r,18 , the simple kinetic theory lies only slightly below our more rigorous calculations. The same experimental sources reported data as for D r,HDO ; in addition Barkan and Luz (2007) reported an averaged value from experiments at 25°C and 40°C, which we plot at 32.5°C. In this case, our results are reasonably close to all the experimental data except for the highest temperature point of Luz et al. (2009).
For D r,17 , one value was reported by Barkan and Luz (2007) of 0.9856, averaged from their experiments at 25°C and 40°C. This is in fair agreement with our result of 0.9847; the simple kinetic theory (0.9836) is somewhat further from experiment. Barkan and Luz (2007) also reported a value for the logarithmic ratio θ diff of 0.5185, which agrees well with our value of 0.5188. While it is often assumed that θ diff is constant with temperature, we find a slight temperature dependence, with θ diff (which we recommend computing from the equations in section 4) decreasing from 0.5207 at 190 K to 0.5167 at 500 K.
We can also examine the relative diffusive fractionation between D and 18 O, defined as φ in Equation 1. As shown in Figure 2, our results for this quantity differ greatly from those of the simple kinetic theory. At room temperature, we agree well with the two most recent experimental studies (Luz et al., 2009;Merlivat, 1978), while again there may be a problem with the temperature dependence of Luz et al. (2009).

Discussion
It is clear from Figure 1a that the simple kinetic theory is significantly in error for D r,HDO compared to more rigorous calculations. Physically, this is not surprising, since the replacement of an H atom by D greatly changes the three principal moments of inertia and the orientations of the two principal axes in the molecular plane, an effect that is missing in the hard-sphere model. Rotational dynamics are important in molecular  diffusion, and Figure 1a suggests that the rotational effect of the D/H substitution on the diffusivity of HDO in air is almost as large as the effect of the mass difference. The much smaller deviation from simple kinetic theory for 17 O and 18 O substitution also makes physical sense, because these substitutions are much closer to the center of mass of the molecule and therefore have little impact on the rotational dynamics.
The small temperature dependence for D r,HDO , and the even smaller dependence for D r,18 , arise solely from the temperature dependences of the collision integrals, which are due to the relative contributions of collisions of various energies and which would cancel when taking the ratios of the diffusivities of different isotopologues only if rotational dynamics were absent (e.g., for noble gases and the hard-sphere model system). At low temperatures, there is a greater contribution from low-energy collisions in which attractive intermolecular forces (dispersion, dipole-quadrupole interactions, etc.) play a larger role. At higher temperatures, high-energy collisions become more important; these are mostly determined by repulsive forces. One would therefore expect the behavior to become more similar to that of the hard-sphere model at high temperatures, which is indeed the case. Since D r,18 differs less than D r,HDO from the temperature-independent hard-sphere model, the temperature dependence of D r,18 is also weaker than that of D r,HDO . We are more confident in our temperature dependence than the uncertainty shading in our figures might suggest; errors in our diffusivity ratios (e.g., due to missing quantum effects) would mainly be systematic in nature, so that any displacement of the true values within the shaded uncertainty would probably lie entirely on one side or the other of our curves. Clearly, the temperature trend of D r,HDO from Luz et al. (2009) is incompatible with our results; we have no hypothesis for why their experiments show such a large apparent temperature dependence. However, we note that the temperature trend of our calculated D r,HDO values is physically more reasonable, since the deviations from the hard-sphere result decrease monotonically with increasing temperature, whereas the data of Luz et al. (2009) cross the constant hard-sphere value sharply at a quite moderate temperature.
To illustrate the importance of the temperature dependence, we consider the well-known model of Jouzel and Merlivat (1984) for kinetic fractionation in snow formation. The overall kinetic factor α k for D/H fractionation is given by where S is the relative saturation (the amount by which S exceeds unity is the fraction by which the vapor is supersaturated) and α eq is the equilibrium fractionation ratio. For definiteness, we consider a temperature of 230 K, where α eq is roughly 1.23 (Merlivat & Nief, 1967), and a relative saturation S of 1.2, which is typical for polar snow formation (Casado et al., 2016;Jouzel & Merlivat, 1984). From Table 1, we obtain D r; HDO ¼ 0:9745, and Equation 5 yields α k ¼ 0:958. However, if we attempt to extrapolate the values of Luz et al. (2009) to 230 K (see Figure 1a), a value of roughly D r; HDO ¼ 0:96 would be obtained, yielding α k ¼ 0:955. While this difference does not seem large, Luz et al. (2009) showed (working with the related quantity φ) that differences of about this magnitude in temperature extrapolation can significantly alter the interpretation of variations in D and 18 O in Antarctic ice cores. A similar calculation with Equation 5 for D r,18 yields a negligible difference, because our results for that quantity are (except for one high-temperature point) in fairly good agreement with those of Luz et al. (2009).
For convenience in practical applications, we fitted simple correlation functions to our calculated diffusivity ratios in air for the temperature range from 190 K to 500 K using the symbolic regression software Eureqa (Schmidt & Lipson, 2009 where T * ¼ T=ð100 KÞ. The correlations reproduce the calculated ratios within ±2 × 10 −5 and thus well within their uncertainties. These correlations are also recommended to be used for the calculation of φ as defined by Equation 1 and θ diff as defined by Equation 2. Similar calculations could be performed for other atmospheric compositions. One interesting possible application is the atmosphere of Mars, where scientists are beginning to use isotopic information to study the planet's water cycle (Krasnopolsky, 2015;Montmessin et al., 2005;Vos et al., 2019) but so far have not included diffusion fractionation in their models. The required calculations for the diffusion of water isotopes in CO 2 could be performed with the recent H 2 O-CO 2 pair potential of Hellmann (2019b).

Conclusion
We performed first-principles molecular kinetic-theory calculations of the diffusivities of water isotopologues in air and used the results to calculate diffusivity ratios for kinetic isotope fractionation. Our results demonstrate that the frequently used hard-sphere kinetic approximation is significantly in error for D/H fractionation, while the experimental result of Merlivat (1978) is accurate. Our calculations provide diffusivity ratios over a wide range of temperature; the temperature dependence is much smaller than that obtained in one study that measured at multiple temperatures (Luz et al., 2009). Because of this discrepancy with the only temperature-dependent experiments, and the importance of kinetic fractionation for ice and snow, it would be desirable for an independent experiment to validate the behavior of D r,HDO at a temperature well below the 10°C limit of existing experimental data.
Our results are described by Equations 6-8, which we recommend to replace the existing experimental data and simple kinetic-theory estimates in all relevant applications.

Data Availability Statement
The intermolecular H 2 O-N 2 and H 2 O-O 2 potentials used in this study are available as Fortran 90 codes in the supporting information of Hellmann (2019aHellmann ( , 2020.  Table S2 are also provided in the NIST Public Data Repository (https://doi.org/10.18434/M32271).