# Isotope Fractionation by Thermal Diffusion in Silicate Melts

Recent experiments on silicate melts have shown that isotopes fractionate significantly along a temperature gradient, such that heavy isotopes accumulate in the cold end and light isotopes accumulate in the hot end [1–4]. The magnitude of this fractionation can be quite large; for example, in molten basalt an isotopic enhancement of several tenths of a percent is observed for ${}^{26}\mathrm{Mg}$ with as little as $50\text{\hspace{0.17em}}\xb0\mathrm{C}$ difference in temperature [2]. Given the precision with which isotopic ratios of many elements can now be measured, it may be possible to detect thermally driven isotope fractionation in magmatic systems [5]. Indeed, the observation of increasingly heavier isotope ratios with magmatic differentiation may provide evidence for such an effect [6,7]. However, in order to understand the origin of observed isotopic variations in magmatic systems, and use them as tools to infer the magmatic history, it is necessary to understand the physical basis for thermally driven isotope fractionation, and the controls on its magnitude and direction.

The physical origins of this isotope fractionation are not well understood. A phenomenological model for this fractionation based on a quantum zero-point effect has recently been put forth [8], but it is not clear that this model can reproduce experimental results with physically reasonable parameters [9]. We show here that the fractionation arises from classical mechanical effects, and that a simple scaling relation can quantitatively predict behavior in silicate melt systems.

The steady-state variation in the ratio of a heavy (mass ${m}_{h}$) and light (mass ${m}_{l}$) isotope along a temperature gradient in silicate melts can be characterized by $\widehat{\delta}$,

where ${C}_{h}(x)$ and ${C}_{l}(x)$ are the concentrations of the heavy and light isotope at position $x$, and ${x}_{0}$ is a reference position (our definition of $\widehat{\delta}$ is not on a per mil basis, as is common). The strength of the isotope fractionation at steady state is often expressed in terms of the parameter $\widehat{\Omega}$:

where ${T}_{0}$ is the temperature for which $\widehat{\delta}=0$. While $\widehat{\Omega}$ may vary with temperature [10], this dependence is found to be small for silicate melts [11]. It has been recognized that the values of $\widehat{\Omega}$ may depend on mass and the atomic interactions, and an empirical scaling was proposed in which $\widehat{\Omega}$ is correlated with the isotope masses, ionic charge, and ionic radius [4].

Here we propose a *nonempirical* scaling that is based on Chapman-Enskog theory [12–15]. This is a rigorous theory with the only assumptions being that the interactions between atoms are binary elastic collisions described by classical mechanics, and that the interatomic forces are spherically symmetric. In Chapman-Enskog theory, the relation for steady-state isotope fractionation along a temperature gradient to leading order is [12–15]

where the value of ${\alpha}_{0}$ is predicted to be of order 1. For example, for hard-sphere systems, ${\alpha}_{0}=105/118$ in the dilute gas limit [14] and it increases with increasing density to a value of $\sim 5$ for a liquid near its maximum density (at the glass transition) [15].

We hypothesize that this scaling relation [Eq. (3)] describes thermal isotope fractionation for complex fluids such as silicate melts. To stringently test the scaling behavior of isotope fractionation, new experiments and molecular dynamics (MD) simulations were carried out. The experiments generated results for isotope fractionation of Sr, Hf, and U in silicate melts (see below for details)—the new data provide a strong test for theory as it expands the range of isotope fractionation data to elements covering most of the periodic table. Our MD simulations on liquid ${\mathrm{MgSiO}}_{3}$ are the first to address isotope fractionation by thermal diffusion in a silicate melt. The MD simulations provide a strong test for the scaling behavior as the simulations must be carried out under “extreme” conditions (in comparison to laboratory experiments) in order to obtain adequate signal-to-noise for the local isotope concentrations with the relatively small number of atoms and short times that are accessible computationally. In particular, the simulations (a) use heavy isotopes that have masses that are 400% of the normal atom, in contrast to natural differences of isotopic masses that are typically 5%–10% or less; (b) use temperature gradients of $\sim {10}^{11}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{K}/\mathrm{m}$, in contrast to gradients of $\sim {10}^{4}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{K}/\mathrm{m}$ in experiments; and (c) are carried out at mean temperature ${T}_{0}=4000\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{K}$, in contrast to experiments carried out at ${T}_{0}\sim 1500\u20132000\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{K}$ (see below for details).

The new experimental data presented here come from thermal migration experiments performed in the University of Illinois, Urbana-Champaign experimental petrology laboratory. The Sr, Hf, and U data represent analyses of solutions derived from dissolving spatially located subsamples from the hot melt-rich end of a thermal migration experiment involving a basaltic starting material. The experiment was performed at 0.5 GPa by placing a graphite capsule containing nominally anhydrous BCR-2 (USGS basalt standard powder) doped with several trace elements at 1000 to 2000 ppm levels into the temperature gradient of a ${\frac{3}{4}}^{\u201d}$ piston cylinder for 34 days. Using the spinel diffusion profile method to determine the temperatures [16], we estimate temperature over the entire capsule ranged from $1260\text{\hspace{0.17em}}\xb0\mathrm{C}$ to $800\text{\hspace{0.17em}}\xb0\mathrm{C}$. Like previous thermal migration experiments, the basaltic starting material evolves to an all melt region at the hot end, a middle region of crystals (orthopyroxene, garnet, clinopyroxene, and ilmenite) plus melt and fine grained crystalline material at the cold end which has not changed greatly in composition relative to the starting material. Probably because the lower half of this experiment has undergone little chemical transport (it appears to be entirely crystalline), no isotope fractionation relative to the starting material is observed such that the calculated fractionations used here only apply to subsamples from the upper half of the experiment (analogous to the Mg and Fe isotopic variations previously observed [17]). The experimental charge was cut into 7 equal $\sim 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ thick sections and each section dissolved in $\mathrm{HF}\mathrm{\text{-}}{\mathrm{HNO}}_{3}$ to make a master solution for that section. Aliquots were removed for Sr, U, and Hf isotope analysis. Sr was purified using standard Sr spec methods with SRM987 also processed simultaneously. The purified Sr solution was run in dry plasma mode using a DSN-100 desolvating nebulizer coupled to a Nu Pasma HR MC-ICPMS (at UIUC) with the chemically processed SRM987 used as a bracketing standard between each solution analysis. ${}^{84}\mathrm{Sr}$, ${}^{86}\mathrm{Sr}$, ${}^{87}\mathrm{Sr}$, and ${}^{88}\mathrm{Sr}$ were collected and offsets for ${}^{88}\mathrm{Sr}/{}^{86}\mathrm{Sr}$ between standard and solutions were determined. Estimated precision is 0.015% ( $2s$), much smaller than the total offset observed which was 0.237%. Because of interference issues on ${}^{84}\mathrm{Sr}$ from Kr, we were unable to obtain satisfactory analyses of ${}^{84}\mathrm{Sr}/{}^{88}\mathrm{Sr}$. Hf was purified using established procedures in the chemistry laboratory at the University of Iowa. An Alfa ICPMS Hf solution was used as a bracketing standard for dry plasma analyses on the Nu Plasma HR (UIUC). Masses 177, 178, 179, 180 were analyzed and mass dependence of the fractionation assessed. This is much larger than the estimated precision of $<0.005\%$. U was analyzed by a double spike technique after purification by standard anion exchange methods. The total offset in ${}^{238}\mathrm{U}/{}^{235}\mathrm{U}$ is 0.075%, significantly larger than the estimated precision of 0.01% ( $2s$). More details of this experiment and analyses are provided elsewhere [18].

Molecular dynamics (MD) simulations use the classical equations of motion to follow the trajectories of atoms as they move under the influence of interatomic forces. Our simulations are carried out for liquid ${\mathrm{MgSiO}}_{3}$, with interatomic forces parameterized by potentials that have previously been shown to accurately model real liquids [19]. Our simulations are carried out for a system composed of 2160 atoms in a cubic simulation box with length $L=3.2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ (corresponding to the density 21 830 ${\mathrm{MgSiO}}_{3}\text{\hspace{0.17em}}\mathrm{moles}/{\mathrm{m}}^{3}$); periodic boundary conditions are used to remove edge effects. The equations of motion are integrated numerically with a time step of 1 fs; simulations are run for durations of 68 ns (this long simulation time is needed to obtain sufficient signal-to-noise to resolve the small differences in the local concentrations for the isotope pairs). To address isotope effects, 23% each of Si and Mg atoms and 8% of O atoms are replaced by “heavier isotopes” with masses 4 times as large. To generate a temperature gradient, the nonequilibrium MD method is used whereby a “cold slice” and “hot slice” of the simulation box are defined, of thickness $d=0.13\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ and oriented perpendicular to the $x$ axis at $x/L=0.25$ and $x/L=0.75$, respectively; at each time step an energy $\Delta \epsilon =0.2{k}_{B}T$ is transferred from the atoms in the cold slice to the atoms in the hot slice, by rescaling the velocities of these atoms [20]. The simulations were carried out with the GROMACS software package [21], which we modified to include the thermal gradient code.

First, we address the accuracy of our MD simulations and relevance to experiment. Our MD results for $\widehat{\delta}$ are compared with experiments in Fig. 1, using the scaling relation of Eq. (3)—the agreement is very good, considering that the MD results follow from fundamental physics with no fitting of any kind to the experimental results. This agreement suggests that the properly scaled results are not strongly dependent on the values of the atomic mass differences and temperature gradients. The MD and experimental results for ${\alpha}_{0}$ are compared in Table I. The MD results correctly reproduce the relative magnitudes of the steady-state isotope fractionation, with $\mathrm{Si}<\mathrm{O}<\mathrm{Mg}$. The MD result for O is within experimental error, and the MD result for Si is within 30% of the experimental value. The MD result for Mg differs by a factor of 2.5 from the experimental result (perhaps this larger difference is due to the high temperature of the simulations, which would act to make all elements behave more similarly because the influence of energetic interactions scales as $1/\mathrm{kT}$).

The scaling of Eq. (3) is tested in Fig. 2 using all available data for silicate melts—these results include experimental fractionations of Si, O, Fe, Mg, Ca in supraliquidus basaltic to andesitic liquids and subliquidus andesite from Refs. [1–4], our new data for Sr, Hf, and U in subliquidus basalt, and our MD results for Si, O, and Mg in ${\mathrm{MgSiO}}_{3}$. This wide range of data largely collapses when the scaling of Eq. (3) is applied. It is noteworthy that the same scaled fractionation occurs whether the system is at super- or subliquidus conditions and regardless of the direction of chemical fluxes. Table I gives the values for ${\alpha}_{0}$ and $\widehat{\Omega}$ (by comparing Eqs. (2) and (3), ${\alpha}_{0}=\widehat{\Omega}({m}_{h}+{m}_{l}){T}_{0}$). The range of ${\alpha}_{0}$ values in silicate melts, $0.5<{\alpha}_{0}<3.5$, is very similar to that from Chapman-Enskog theory for hard-sphere systems, i.e., $1<{\alpha}_{0}<5$ [15]. Furthermore, there is a distinct difference in the values of ${\alpha}_{0}$ for network formers (Si and O) and network modifiers (Mg, Fe, Ca, Sr, U, Hf), with smaller values for the network formers and little variation within each of these two groups; this difference is not evident in terms of the values of $\widehat{\Omega}$.

Since both Chapman-Enskog theory (even when evaluated for hard-sphere systems) and MD simulations give quantitative agreement with experiment with no adjustable parameters, it can be concluded that the origins of isotope fractionation lie in the classical mechanical collisions between pairs of particles (and are not solely a product of quantum effects, as has been suggested [8]). The mathematical complexity of the Chapman-Enskog analysis obscures the underlying physical picture.

Therefore, we use a simpler system to demonstrate the origin of the effect. Consider a head-on collision between a heavier atom and lighter atom with short-ranged interactions. We can define a persistence fraction, $p$, as the fraction of the collisions in which the motion of an atom persists in the same direction after the collision. From simple considerations following from the conservation of energy and momentum, it can be shown that the motion of the heavy atom may persist in the same direction if its velocity is sufficiently high, but the light atom can never do so. We consider the collision of a heavier ( $h$) and lighter ( $l$) particle in one dimension, where the interparticle interaction is short-ranged. The collision satisfies the conservation of momentum,

and the conservation of energy,

where ${m}_{k}$, ${v}_{k,0,}$, and ${v}_{k,f}$ are the mass, initial velocity, and final velocity of particle $k$ (the initial and final velocities correspond to positions at which the particles are farther apart than the range of their interaction). The final velocities can be determined from these two conservation equations,

(and an analogous equation holds for the light particle). The motion of the heavier particle persists in the same direction after the collision when the initial velocities of the two particles satisfy the criterion

(note the lighter particle can never persist in the same direction). In a fluid at temperature $T$, particles have a distribution of velocities—thus there will be a probability ${p}_{h}$ that the motion of the heavy particle persists in the same direction after a collision. The velocity distribution is a Gaussian function with mean of zero and standard deviation ${\sigma}_{i}=\sqrt{2k{T}_{i}/{m}_{i}}$ (the Maxwell distribution in one dimension). Using the fact that the probability distribution for the ratio of two Gaussian distributed variables is described by a Cauchy distribution [22], the probability that the criterion in Eq. (7) is satisfied is given by the cumulative Cauchy function,

which can be simplified in the limit that $\frac{{m}_{h}-{m}_{l}}{({m}_{h}{m}_{l}{)}^{1/2}}$ is small, to yield

Equation (9) shows that the motion of a heavier particle is more likely to persist in the same direction after a collision with a lighter particle when it comes from a region with higher temperature. While this analysis was for a one-dimensional collision involving only two particles, evidence for the higher persistence fraction of heavier atoms has been found in MD simulations of bulk systems, via differences in the velocity autocorrelation function [23]. For lighter isotopes, the velocity autocorrelation function became negative at intermediate times, due to the recoil after a collision. In contrast, for the heavier isotopes the velocity autocorrelation function never became negative, showing that recoil after a collision occurred less frequently for the heavy isotope.

This mass-dependent persistence effect can be visualized by an analogy to American football—in a collision between a (heavier) lineman and a (lighter) cornerback, the lineman can push his way through the cornerback *if he has enough speed*. In the same way, Eq. (9) shows that it is more probable that a heavy particle will move from the hot side to the cold side than vice versa. For this reason, the cold side will become enriched in heavy isotopes, leaving the hot side enriched in light isotopes, as has been found in all experimental studies.

## References

- T. K. Kyser, C. E. Lesher, and D. Walker, Contrib. Mineral. Petrol.
**133**, 373 (1998). - R. M. Richter, E. B. Watson, R. Mendybaev, F.-Z. Teng, and P. Janney, Geochim. Cosmochim. Acta
**72**, 206 (2008). - F. M. Richter, E. B. Watson, R. Mendybaev, N. Dauphas, B. Georg, J. Watkins, and J. Valley, Geochim. Cosmochim. Acta
**73**, 4250 (2009). - F. Huang, P. Chakraborty, C. C. Lundstrom, C. Holmden, J. J. G. Glessner, S. Kieffer, and C. E. Lesher, Nature (London)
**464**, 396 (2010). - C. C. Lundstrom, Geochim. Cosmochim. Acta
**73**, 5709 (2009). - J. A. Schuessler, R. Schoenberg, and O. Sigmarsson, Chem. Geol.
**258**, 78 (2009). - P. S. Savage, R. B. Georg, H. M. Williams, K. W. Burton, and A. N. Halliday, Geochim. Cosmochim. Acta
**75**, 6124 (2011). - G. Dominguez, G. Wilkins, and M. H. Thiemens, Nature (London)
**473**, 70 (2011). - D. J. Lacks, J. A. Van Orman, and C. E. Lesher, Nature (London) (to be published).
- F. M. Richter, Nature (London)
**472**, E1 (2011). - F. Huang, P. Chakraborty, C. C. Lundstrom, C. Holmden, J. J. G. Glessner, S. Kieffer, and C. E. Lesher, Nature (London)
**472**, E2 (2011). - D. Enskog, Phys. Z.
**12**, 533 (1911). - S. Chapman, Phil. Trans. R. Soc. A
**217**, 115 (1918). - W. Furry, R. Jones, and L. Onsager, Phys. Rev.
**55**, 1083 (1939). - J. M. Kincaid, E. G. D. Cohen, and M. Lopez de Haro, J. Chem. Phys.
**86**, 963 (1987). - E. B. Watson, D. A. Wark, J. D. Price, and J. A. Van Orman, Contrib. Mineral. Petrol.
**142**, 640 (2002). - F. Huang, C. C. Lundstrom, J. Glessner, A. Ianno, A. Boudreau, J. Li, E. C. Ferré, S. Marshak, and J. DeFrates, Geochim. Cosmochim. Acta
**73**, 729 (2009). - C. J. Bopp IV, Ph.D. thesis, University of Illinois, Urbana 2010.
- D. J. Lacks, and J. A. Van Orman, Geochim. Cosmochim. Acta
**71**, 1312 (2007). - P. Jund and R. Jullien, Phys. Rev. B
**59**, 13 707 (1999). - B. Hess, C. Kutzner, D. van der Spoel, and E. Lindhal, J. Chem. Theory Comput.
**4**, 435 (2008). - E. W. Weisstein, “Normal Ratio Distribution.” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/NormalRatioDistribution.html;; “Cauchy Distribution.” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/CauchyDistribution.html.
- N. Kiriushcheva and S. V. Kuzmin, Physica (Amsterdam)
**352A**, 509 (2005).