## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Ion desolvation as a mechanism for kinetic isotope fractionation in aqueous systems

Edited by Mark H. Thiemens, University of California at San Diego, La Jolla, CA, and approved September 21, 2012 (received for review May 14, 2012)

## Abstract

Molecular dynamics simulations show that the desolvation rates of isotopes of Li^{+}, K^{+}, Rb^{+}, Ca^{2+}, Sr^{2+}, and Ba^{2+} may have a relatively strong dependence on the metal cation mass. This inference is based on the observation that the exchange rate constant, *k*_{wex}, for water molecules in the first hydration shell follows an inverse power-law mass dependence (*k*_{wex} ∝ *m*^{−γ}), where the coefficient γ is 0.05 ± 0.01 on average for all cations studied. Simulated water-exchange rates increase with temperature and decrease with increasing isotopic mass for each element. The magnitude of the water-exchange rate is different for simulations run using different water models [i.e., extended simple point charge (SPC/E) vs. four-site transferrable intermolecular potential (TIP4P)]; however, the value of the mass exponent γ is the same. Reaction rate theory calculations predict mass exponents consistent with those determined via molecular dynamics simulations. The simulation-derived mass dependences imply that solids precipitating from aqueous solution under kinetically controlled conditions should be enriched in the light isotopes of the metal cations relative to the solutions, consistent with measured isotopic signatures in natural materials and laboratory experiments. Desolvation effects are large enough that they may be a primary determinant of the observed isotopic fractionation during precipitation.

Nonequilibrium processes are generally recognized as important influences on isotopic fractionation during mineral precipitation, especially at temperatures below a few hundred degrees Celsius (1, 2). Recent work in “nontraditional” stable isotopes (Li, Mg, Ca, Fe, Cd, Cu, etc.) has confirmed this inference and shown that nonequilibrium fractionation must be accounted for to understand global biogeochemical cycles involving these elements (2⇓⇓–5). A key observation is that light isotopes are preferentially incorporated into precipitating solids (6⇓–8)—the opposite direction of enrichment expected for equilibrium fractionation, which depends on bond energetics (9, 10). The origin of this nonequilibrium light-isotope enrichment is a key unknown in the isotopic geochemistry of mineral growth.

Calcite grown from aqueous solutions provides an excellent illustration of light-isotope enrichment during precipitation. Synthetically precipitated and natural samples of calcite and aragonite are fractionated relative to aqueous Ca^{2+} such that δ^{44/40}Ca in calcite is lower by 0.5–2‰ (7, 11). Equilibrium Ca isotope fractionation between Ca^{2+}_{(aq)} and calcite has not been measured in the laboratory, but studies of deep-sea sedimentary pore fluids suggest that the equilibrium fractionation is negligible: ε_{eq} ∼ 0.0 ± 0.1‰ (12). Growth of calcite from seawater-like aqueous solutions is not likely to be diffusion-limited for Ca, so the isotopic effects are inferred to be due to kinetic effects occurring at the solid/fluid interface (11, 13). Diffusion in liquid water can result in kinetic isotopic fractionation (14⇓⇓–17), but in the case of Ca isotopes, the available data suggest that the largest possible light-isotope enrichment from diffusion-limited growth is ∼0.4‰ in δ^{44/40}Ca, which is too small to explain the observations (17).

The growth of crystals from solution requires that the constituent solute ions undergo three steps: diffusion of solvated ions to the crystal/liquid interface, adsorption onto the solid surface, and incorporation into the crystal lattice. The last two steps entail a progressive desolvation of the ions; therefore, these steps may be rate-limited by solute dehydration (18, 19). In support of this inference, molecular dynamics (MD) simulations of barium sulfate growth suggest that the kinetics of Ba^{2+} attachment—which govern barite nucleation and growth (20)—may be rate-limited by the partial desolvation of the metal to form an inner-sphere surface complex (21). Likewise, experimental (22, 23), theoretical (13, 24), and molecular simulation studies (25) suggest that the kinetics of metal attachment at calcite surfaces may be related to the dehydration frequency of the metal near the surface. More broadly, the rate constants of several metal ligand-exchange reactions [Mg and Ni binding to a range of organic ligands (26, 27), the dissolution of orthosilicate minerals containing a range of divalent metals (28, 29)] have been shown to correlate with the water-exchange rates of metals in liquid water, *k*_{wex} (the inverse of the residence time of water in the first solvation shell of the metal).

The correlation between metal-water exchange frequencies (*k*_{wex}) and the rates of other metal-ligand exchange reactions suggests—regardless of whether this correlation arises from desolvation being a rate-limiting step or from inherent similarities between different ligand-exchange reactions involving the same metal ion—that the isotopic mass dependence of *k*_{wex} may yield important insight into kinetic isotope fractionation during mineral growth and, more broadly, during metal biogeochemical cycling in aqueous systems. A mass dependence of *k*_{wex}, although never previously reported, is in fact quite plausible within the framework of reaction rate theory, where *k*_{wex} is the product of a transmission coefficient (κ) of unknown mass dependence and a transition state theory (TST) rate constant (*k*^{TST}):

As explained in more detail in *Results and Discussion*, *Theoretical Basis of the Mass Dependence of k*_{wex} if metal–water distance is selected as the reaction coordinate (the simplest and most widespread choice for this type of reaction), *k*^{TST} has an inverse square-root dependence on metal–water reduced mass μ_{i}; therefore, *k*_{wex} must vary with isotopic mass unless the unknown mass dependence of κ exactly cancels the strong mass dependence of *k*^{TST}.

MD simulations are routinely used to investigate the *k*_{wex} values of alkali and alkaline earth metals in liquid water (30⇓⇓⇓⇓–35) and these simulations readily allow for the determination of the mass dependence of *k*_{wex} (15, 17). Here, we use MD simulations in conjunction with reaction rate theory calculations to determine the isotopic mass dependence of *k*_{wex} for several common cations at three different temperatures, and we show that this dependence could account for the observed Ca isotopic effects attending calcite growth from aqueous solution.

## Results and Discussion

### Isotopic Mass Dependence of Metal Desolvation Rates.

For each 8-ns segment of our 16-ns simulations, we calculated the time-correlation function *C*_{desolv}(*t*), the average number of molecules that remain continuously located in the first solvation shell during a time interval of duration *t*. Water-exchange rates *k*_{wex} were obtained by applying a linear least-squares regression to the first 20 ps of the first-order decay expression

Grossfield (36) found that fitting the first 10 ps of similar autocorrelation functions was sufficient to model the lifetime of water molecules in the first solvation shell surrounding monovalent ions. Exchange rates calculated using Eq. **2** were evaluated for mass dependence by constructing log-log plots of *k*_{wex,i} vs. *m _{i}* (for which the subscript

*i*indicates the isotope) at each of the three simulation temperatures for the alkali (Fig. 1

*A*and Fig. S1

*A*) and alkaline earth metals (Fig. 1

*B*and Fig. S1

*B*). Data were fit with a linear least-square regression weighted by ±1 SE on each value of

*k*

_{wex,i}of the form

without constraints on either *B* or γ. Values of *k*_{wex,i} are compiled in Dataset S1; values of *B* and γ are reported in Dataset S2*A*. Although Eq. **3** cannot be valid for infinitely high masses, it provides an accurate mathematical description of the relation between solute mass and water-exchange rate for the mass range of most elements of interest (*m _{i}* ∼ 3–100 Da). [Within the precision of our simulations, plots of

*k*

_{wex,i}vs. cation mass and plots of

*k*

_{wex,i}vs. ion-water reduced mass μ

_{i}both yield statistically robust fits to the

*k*

_{wex}data using Eq.

**3**(values of

*B*and γ calculated using μ

_{i}instead of

*m*

_{i}in Eq.

**3**are reported as

*B*′ and γ′ in Dataset S2

*B*), but the coefficient of variation of γ is, on average, 8% larger for regressions of

*k*

_{wex,i}vs. ion-water reduced mass. We therefore discuss our results in terms of γ-values (i.e., as a function of cation mass). We note that our key findings would be almost identical if we used γ′-values except insofar as the average γ-values of the six metals studied here are all identical (as described below), whereas the average γ′-values show a significant dependence on the type of metal with no obvious trend.] Calculated γ-values are significantly greater than zero in most cases. Therefore, metal desolvation rates have a significant isotopic mass dependence that is well described with the inverse power-law relation

The parameter *B* is significantly influenced by temperature, the type of metal, and the choice of water model (Dataset S2*A*). The temperature dependence of *B* is consistent with an Arrhenius relation [*B* ∝ e^{−Ea/RT}, where *E*_{a} (kJ⋅mol^{−1}) is the phenomenological “activation energy” of *B*]. Our results yield an average *E*_{a} = 13.7 ± 2.6 kJ⋅mol^{−1} (confidence intervals are reported as twice the SEM) for all solutes and water models investigated here, in agreement with the conceptual view of metal desolvation as a thermally activated process. The solute dependence of *B* is consistent with X-ray diffraction data that show that solvating waters become more strongly oriented and more tightly bound to the central ion as surface charge density increases (30, 36⇓–38). This relationship translates into lower values of the radial distribution function, *g*(*r*), of each cation at the boundary between the first and second solvation shell (Fig. S2; further details in *SI Text* and Table S1) and lower values of *B* and *k*_{wex} for smaller, higher-valence metals. The calculated values of *k*_{wex} depend more strongly on the nature of the solute than on the solvent, as noted by Dunand et al. (39); nevertheless, they are significantly influenced by the choice of water model [an indication, as pointed out by Wang et al. (40), of the difficulty in accurately predicting water-exchange rates in molecular simulations]. The value of the ratio *B*_{TIP4P}/*B*_{SPC/E} is larger than 1 for four of the solutes studied here (*B*_{TIP4P}/*B*_{SPC/E} = 1.45 ± 0.03, 1.11 ± 0.04, 2.20 ± 0.94, 1.58 ± 0.05 for Li^{+}, K^{+}, Sr^{2+}, Ba^{2+}, respectively; average values calculated from data at 278, 298, and 323 K), equal to 1 for one solute (*B*_{TIP4P}/*B*_{SPC/E} = 1.02 ± 0.06 for Rb^{+}), and smaller than 1 for one solute (*B*_{TIP4P}/*B*_{SPC/E} = 0.66 ± 0.10 for Ca^{2+}). We hypothesize that this water-model dependence of *B* may result from a combination of dynamical and structural effects. Dynamically, water-exchange rates and solute diffusion coefficients may be positively coupled (33), causing the four-site transferrable intermolecular potential (TIP4P) water model, which predicts faster self-diffusion of liquid water, to predict higher *B* values than the extended simple point charge (SPC/E) model. Structurally, small differences in coordination may influence the lability of the solvation shell. For example, in the case of Ca^{2+}, the TIP4P water model yields an eightfold coordination that may be particularly stable (Table S1).

In contrast to the sensitivity of *B* to temperature, water model, and type of metal, γ-values calculated with Eq. **3** show no significant dependence on these factors. Therefore, despite the water model-dependent differences in the absolute values of *k*_{wex,i}, the mass dependence of *k*_{wex} is much less sensitive to simulation methodology. Ratios of γ-values calculated at different temperatures show that γ is independent of *T* [γ_{298K}/γ_{278K} = 0.83 ± 0.23 and γ_{323K}/γ_{278K} = 1.26 ± 0.40; average values calculated for all water models and all solutes except Ba^{2+}, for which Eq. **3** yielded a negative γ-value at 278 K (the largest errors on γ should occur where *k*_{wex} values are smallest, i.e., for divalent metals at 278 K)]. Ratios of γ-values obtained with different water models show that γ has little or no dependence on the choice of water model [γ_{TIP4P}/γ_{SPC/E} = 0.64 ± 0.26, 1.09 ± 0.36, 0.93 ± 0.23, 1.49 ± 1.05, 0.32 ± 0.36, 0.41 ± 4.03 for Li^{+}, K^{+}, Rb^{+}, Ca^{2+}, Sr^{2+}, and Ba^{2+}, respectively, with an average value for all solutes and all temperatures (excepting Ba^{2+} at 278 K) of γ_{TIP4P}/γ_{SPC/E} = 1.07 ± 0.37]. Averaging all gamma values obtained for the same solute, we obtain γ_{Li} = 0.057 ± 0.015, γ_{K} = 0.049 ± 0.008, γ_{Rb} = 0.051 ± 0.007, γ_{Ca} = 0.043 ± 0.015, γ_{Sr} = 0.064 ± 0.035, and γ_{Ba} = 0.032 ± 0.021. Our conclusion is that γ is independent of the type of metal cation within the precision of our simulation methodology. Furthermore, γ = 0.049 ± 0.009 adequately captures the results for all metals, temperatures, and water models investigated here, despite a range of *k*_{wex} values of over two orders of magnitude.

### Theoretical Basis of the Mass Dependence of *k*_{wex}.

Transition state theory is widely used to explain the rates of metal desolvation and similar phenomena, such as ion pair dissociation (31, 35, 41, 42). Here, recalling that the rate constant for dissociation of a metal–water complex can be written as the product of κ and *k*^{TST} (Eq. **1**), we show that this theory is consistent with our finding that *k*_{wex} decreases with metal isotopic mass, and also consistent with the magnitude of the mass dependence that comes from the MD simulations. Our object being to test whether reaction rate theory can predict a mass dependence of *k*_{wex}, we select the simplest and most widely used reaction coordinate (metal–water distance *r*) for evaluating κ and *k*^{TST} (31, 33, 41) [more conceptually complex choices of reaction coordinate are equally sound and may provide additional insights into *k*_{wex} (35, 42)]. For each isotope (*i*) within this parameterization, *k*^{TST} depends only on the metal–water reduced mass μ_{i} and the centrifugally averaged effective potential *W*_{eff,i}(*r*) {defined as *W*_{eff,i}(*r*) = *W _{i}*(

*r*) – 2β

^{−1}ln(

*r*/

*r*

^{‡}), where

*W*(

_{i}*r*) = −β

^{−1}ln[

*g*(

_{i}*r*)], β = (

*k*

_{B}

*T*)

^{−1},

*k*

_{B}is Boltzmann’s constant, and

*r*

^{‡}is the position of the top of the transition state barrier}:

The transmission coefficient κ accounts for dynamical effects that are not included in *k*^{TST}, such as solvent reorganization kinetics during the dissociation of the metal–water pair (43⇓–45). Here, we predict κ for each metal isotope using Kramers’ model (46⇓–48),

in which the reaction is modeled as a diffusive process characterized by a barrier frequency *ω _{i}* and a static friction coefficient

*ζ*. The isotope-specific barrier frequency,

_{i}*ω*, for each 8-ns simulation was calculated by fitting

_{i}*W*

_{eff,i}(

*r*) [derived as indicated above from

*g*(

_{i}*r*)] between

*r*

^{‡}± 0.5 Å with the nonlinear inverse parabolic expression

*W*

_{eff,i}(

*r*) =

*W*

_{eff,i}(

*r*

^{‡}) – 1/2μ

_{i}

*ω*

_{i}^{2}(

*r*−

*r*

^{‡})

^{2}(30, 47). The static friction coefficient was calculated with the implicit approximation that solute diffusion in water and solute-water dissociation have similar friction coefficients, and by applying the Stokes–Einstein relation

*ζ*=

_{i}*k*

_{B}

*T*/

*D*to the self-diffusion coefficient

_{i}*D*of the solute in water. The latter was obtained from its velocity autocorrelation function via the Green–Kubo relation as described by Frenkel and Smit (49) and following the procedure of Bourg and Sposito (15) for 8-ns simulations. Alternative models of κ have been proposed that account for the frequency dependence of ζ (47) as well as the anharmonicity of the reaction barrier and the dependence of ζ on reaction coordinate (50). For Na

_{i}^{+}desolvation in ambient liquid water, these more complex models predict κ values 1.7× to 5× larger than obtained with Eq.

**6**, in line with MD simulation results (31). Below, we show that Kramers’ model (Eq.

**6**) is sufficient to explain the mass dependence of

*k*

_{wex}, although it may underestimate κ by a factor of ∼3.

The *k*_{wex} values predicted using reaction rate theory (Eqs. **1**, **5**, and **6**), hereafter named *k*_{wex_RR}, are plotted in Fig. 2 and Fig. S3 and compiled in Dataset S3. The regression lines plotted in Fig. 2 and Fig. S3 (black and dashed red for SPC/E and TIP4P, respectively) are weighted linear least-squares fits to the *k*_{wex_RR} data using Eq. **3**. Errors on individual *k*_{wex_RR,i} values were determined solely from errors on κ* _{i}*, which were determined using a Monte Carlo routine to propagate errors on

*ζ*and

_{i}*ω*. The SE on

_{i}*D*during each 8-ns simulation (used to calculate the SE on

_{i}*ζ*) was calculated as σ(

_{i}*D*) = σ(

_{i}*D*)/√

_{i,b}*n*in which

_{b}*D*is the block-averaged diffusion coefficient and

_{i,b}*n*is the number of blocks. Each 8-ns simulation was divided into twenty 0.4-ns blocks to generate the statistical uncertainties. Predicted

_{b}*k*

_{wex_RR}values (Fig. 2 and Fig. S3) display the same general trends as those obtained by MD simulation (Fig. 1 and Fig. S1):

*k*

_{wex_RR}increases with temperature, decreases with the charge density of the metal, and is sensitive to the choice of water model. Water-exchange rates predicted with Eqs.

**1**,

**5**, and

**6**are consistently ∼0.5 log units lower than those obtained by MD simulation, in agreement with the observation by Rey and Hynes (31) that Kramers’ model underestimates κ by a factor of ∼3 in the case of Na

^{+}. The linear regressions in Fig. 2 and Fig. S3 are consistent with an inverse power-law mass dependence of

*k*

_{wex_RR}(Eq.

**4**) with slopes similar to those obtained by MD simulation. Thus, reaction rate theory calculations broadly support our MD simulation findings. The γ-values obtained from the linear regressions in Fig. 2 and Fig. S3 may be influenced by the approximations made in Eq.

**6**, but their average value (γ = 0.038 ± 0.021) is close to our MD simulation result (γ = 0.049 ± 0.009). The larger error on the reaction rate theory γ is due to Ba calculations having an average γ-value of −0.001 ± 0.014. Ignoring Ba in both the MD and reaction rate theory determinations results in average γ-values of 0.053 ± 0.007 and 0.046 ± 0.017, respectively.

Reaction rate theory also provides key insights into the fundamental basis of the observed mass dependence of desolvation frequency. In Fig. 3, we plot the contributions of *k*^{TST} and κ to our data on log *k*_{wex_RR,i} vs. log *m _{i}* for the representative case of Li

^{+}in SPC/E water at 298 K. The rate constant

*k*

^{TST}has an inverse square-root dependence on ion-water–reduced mass, as expected from Eq.

**5**. Hence, if κ had no isotopic mass dependence, one would expect a strong, inverse square-root dependence of

*k*

_{wex_RR}on the reduced mass μ

_{i}. Simulations by Møller et al. (33) suggest that such a strong mass dependence of

*k*

_{wex}exists for Li

^{+}in supercritical water, where κ should be close to unity. Our calculations using Eq.

**6**, however, predict that κ has a positive mass dependence in liquid water that attenuates the overall inverse mass dependence of

*k*

_{wex_RR}.

Conceptually, the transmission coefficient κ describes the fraction of trajectories—among all trajectories where a metal–water complex reaches the transition state—where the metal–water complex returns to its original state without dissociating (so-called barrier recrossings). A direct interpretation of our results is that the chance of barrier recrossing decreases as solute mass increases. Additional insight can be gained by noting that κ also describes the influence of solvent dynamics (the collective dynamics of all molecules except those directly involved in the metal–water dissociation reaction) on the reaction rate (41) [*k*^{TST} is influenced by solvent structure through *W*(*r*) in Eq. **5**, but it is independent of solvent dynamics]. Our results show that solvent dynamics have a decreasing influence on *k*_{wex} as solute mass increases, perhaps because the slower motions of heavy solutes allow more time for solvent reorganization to accommodate the dissociation of the metal–water complex.

On the broader topic of kinetic isotope effects, our results yield two important insights. First, the inverse mass dependence of *k*_{wex} (Eq. **4**) may provide a clue into the fundamental basis of the isotopic mass dependence of solute diffusion coefficients in liquid water (*D _{i}* ∝

*m*

_{i}^{−β}) (14⇓⇓–17), because solute

*D*values are positively correlated with

*k*

_{wex}, at least in the case of Li

^{+}and Na

^{+}in liquid water (33). Second, our simulations, in which isotopes have identical interatomic interaction parameters and differ only by their masses, show that—in addition to the expected quantum mechanical effects—classical mechanical effects can play an important role in causing kinetic isotope fractionations in geochemical systems, as pointed out by Young et al. (51). This finding is consistent with studies showing that Newtonian mechanics are the primary cause of carbon kinetic isotope effects associated with the dimerization of cyclopentadiene (52) and may contribute significantly to a range of oxygen kinetic isotope effects (53).

### Implications for Isotopic Fractionation in Geochemical Processes.

Our finding of the isotopic mass dependence of *k*_{wex} implies that any ligand-exchange reaction whose rate constant correlates with *k*_{wex} [mineral dissolution (28, 29), metal binding to organic ligands (26, 27), and possibly metal attachment to mineral surfaces (13, 18, 19, 21)] should exhibit significant kinetic isotope fractionation. For such reactions, the overall kinetic isotope fractionation (α_{kinetic}) has a maximum value in conditions where the forward reaction rate is much larger than the backward rate. This maximum value can be readily determined from the ratio of water desolvation rates of two different isotopes (*i* and *j*) of a given metal ion:

In cases where the backward rate is nonnegligible, calculation of α_{kinetic} is straightforward but requires knowledge of the ratio of forward to backward reaction rates (13). In cases for which these rates are not well constrained, the maximum value of α_{kinetic} calculated via Eq. 7 may be used to deduce the ratio of forward to backward rates from the measured overall kinetic isotope fractionation.

Using the global MD-derived γ-value (γ = 0.049 ± 0.009) and the isotopic masses given in Dataset S1, we can calculate the maximum value of α_{kinetic} for the metals studied here. Results show (Table 1) that the mass dependence of *k*_{wex} can account for relatively large kinetic isotope fractionations equal to or greater than the magnitude of observed light-isotopic enrichments in solids containing Li (8), Ca (7), and Ba (6). Each of these three studies (6⇓–8) documents mineral precipitation from aqueous solutions with ionic concentrations greater than infinite dilution (as considered here) and containing other ligands in addition to water. Although the magnitude and direction of the measured fractionations in (6⇓–8) coincide with those predicted by our simulations, the extent to which our results are directly transferable to these multicomponent systems—or to systems in which metal desolvation is extremely rapid relative to crystal growth rate or the metal is bound to a surface, an oligomer, or to an anion as a solvated contact ion pair—is unclear and should be the focus of further studies.

In the case of inorganic calcite precipitation, a recent study by Nielsen et al. (24) shows that Ca isotopic signatures in calcite reflect contributions from kinetic processes, which become more important at higher growth rates. The kinetic contribution is expected to be important (13, 24) because the characteristic time-scale of Ca attachment to calcite kink sites—on the order of microseconds (24, 54)—is faster than the characteristic time-scale of Ca detachment by a factor of 2 to >20 at typical supersaturation values used in laboratory experiments. Calcite growth rates and Ca isotope fractionations are consistent with a ratio of ^{44}Ca^{2+} to ^{40}Ca^{2+} binding rate constants (*k*_{44Ca}/*k*_{40Ca}) of the order 0.9920–0.9963 (24). Our molecular-scale value of α_{kinetic} = 0.9953 ± 0.0009 (^{44/40}Ca; Table 1) for a single-step desolvation is similar in magnitude to the observed effects associated with calcite growth. This similarity may be a coincidence, but our interpretation is that our simulations show that isotope effects associated with a key step in the crystal growth process (54) are of the correct sign and magnitude to explain the observed experimental results and hence are a potential source of the kinetic fractionation associated with crystal growth.

Because we calculated the desolvation rates of only a subset of alkali and alkaline earth metal cations, there is no a priori expectation that our results should apply to other solutes, or even to other metals. Nevertheless, the few available data on kinetic isotope effects associated with ligand-exchange reactions of aquated metals other than those studied here are remarkably consistent with Eq. **7** and our predicted mass dependence of *k*_{wex}. For example, Eq. **7** with γ = 0.049 ± 0.009 predicts α_{kinetic}(^{114/110}Cd) = 0.9982 and α_{kinetic}(^{65/63}Cu) = 0.9985, in agreement with the kinetic isotope fractionations measured during cadmium uptake by freshwater phytoplankton [α_{f}(^{114/110}Cd) = 0.9986 ± 0.0006] (55) and during copper binding to azurin enzyme, a kinetically controlled one-step process [α_{f}(^{65/63}Cu) = 0.99847] (56).

## Simulation Methods

Molecular dynamics simulations involving one cation and 550 water molecules in a periodically replicated, cubic simulation cell were performed using the program MOLDY 3.6 (57). Each simulation was carried out for 16 ns (with a 1-fs time step) in the NVE ensemble and was preceded by 202 ps of equilibration at the desired temperature. Due to the absence of a counter ion and the application of periodic boundary conditions, these simulations approximate conditions of infinite dilution (15⇓–17, 30⇓⇓⇓⇓–35). We chose to model all alkali and alkaline earth metal cations that have more than one naturally occurring stable isotope (Li^{+}, K^{+}, Rb^{+}, Ca^{2+}, Sr^{2+}, Ba^{2+}), with the exception of Mg^{2+}, for which preliminary calculations showed too few water-exchange events to accurately calculate *k*_{wex} via direct simulation. Each metal cation was simulated at three different temperatures (278, 298, and 323 K) and for a range of solute isotopic masses (both naturally occurring and hypothetical) from *m* = 3 to 176 Da (Dataset S1). Simulation cell volume was fixed to the corresponding liquid water densities at 1 atm (ρ = 1.000, 0.997, and 0.988 g⋅cm^{−3}, respectively). Ion-water interaction potentials (Table S2) were taken from Åqvist (58). To determine the sensitivity of simulation results to the choice of water model, we replicated each simulation with two widely used water models: the extended simple point charge (SPC/E) model (59) and the four-site transferrable intermolecular potential (TIP4P) model (60). A comparison of the models’ predictions for the physical properties of liquid water is given in *SI Text*.

Long-range interactions are treated by Ewald summation in three dimensions with a real space cutoff at *r*_{c} = 9.63 Å, a reciprocal space cutoff at *k*_{c} = 1.91 Å^{−1}, and a sensitivity/screening factor of α = 0.3151 Å^{−1}, resulting in a Ewald sum accuracy of >99.9%. Molecular trajectories are calculated using a modified form of the Beeman algorithm, the most accurate of all “Verlet-equivalent” algorithms (57) to solve the Newton–Euler equations of motion. We divided each simulation into two 8-ns blocks for data analysis. The total energy varied by less than 1% during each 8-ns simulation, indicating excellent energy conservation and justifying our decision to not couple the simulation box to an external thermostat (15, 17).

Unlike Impey et al. (30) and others (31, 32, 35) who allowed water molecules to make temporary (i.e., ≤2 ps) excursions from the hydration sphere and still be counted as present within the sphere, in our calculations of the time-correlation function *C*_{desolv}(*t*) we do not count water molecules as being present in the first hydration sphere once they move beyond a distance of *r _{min}* from the cation. [We tested the “2-ps excursion” method of Impey et al. (30) and found that it yields lower

*k*

_{wex}values (as expected), but the same mass dependence of

*k*

_{wex}in almost all cases. The exception was Ca

^{2+}, the metal with the lowest water-exchange rates of all of the solutes in this study, for which we hypothesize that the

*k*

_{wex}values calculated with the 2-ps excursion method may have been too small to accurately probe their mass dependence.] The desolvation correlation functions for all monovalent cations are qualitatively the same (i.e., e-folding times of ∼10–50 ps), regardless of temperature or water model; likewise, the desolvation correlation functions calculated for divalent cations all have e-folding times of ∼200–600 ps.

## Acknowledgments

Support for this work was provided by the Center for Isotope Geochemistry, funded by US Department of Energy Office of Basic Energy Sciences Contract DE-AC02-05CH11231, as well as the Center for Nanoscale Control of Geologic CO_{2}, an Energy Frontier Research Center.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: aehofmann{at}lbl.gov.

Author contributions: A.E.H. and I.C.B. designed research; A.E.H. performed research; A.E.H., I.C.B., and D.J.D. analyzed data; and A.E.H., I.C.B., and D.J.D. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1208184109/-/DCSupplemental.

## References

- ↵
- ↵
- Johnson CM,
- Beard BL,
- Albarède F

- ↵
- Stumm W,
- Morgan JJ

- ↵
- ↵
- Baskaran M

- ↵
- Böttcher ME,
- et al.

*Geophys Res Abstr*13:EGU2011-4982-1. - ↵
- Nielsen LC,
- Druhan JL,
- Yang W,
- Brown ST,
- DePaolo DJ

- ↵
- Wunder B,
- Meixner A,
- Romer RL,
- Jahn S

- ↵
- ↵
- Bigeleisen J

- ↵
- DePaolo DJ

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- deBoer RB

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Frenkel D,
- Smit B

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Zhu XK,
- et al.

- ↵
- ↵
- ↵
- ↵

## Citation Manager Formats

## Sign up for Article Alerts

## Jump to section

## You May Also be Interested in

### More Articles of This Classification

### Physical Sciences

### Earth, Atmospheric, and Planetary Sciences

### Related Content

- No related articles found.

### Cited by...

- Uranium isotope fractionation by abiotic reductive precipitation
- Barium isotope evidence for pervasive sediment recycling in the upper mantle
- Kinetic Fractionation of Non-Traditional Stable Isotopes by Diffusion and Crystal Growth Reactions
- Isotopic Gradients Across Fluid-Mineral Boundaries
- Molecular Simulation of CO2- and CO3-Brine-Mineral Systems