The distribution of stable and radiogenic isotopes within and among phases provides a critical means of quantifying the origin, residence and cycling of materials through terrestrial reservoirs (Wahl and Urey 1935; Epstein and Mayeda 1953; Johnson et al. 2004; Eiler 2007; Porcelli and Baskaran 2011; Wiederhold 2015). While isotopic variability is globally observable, the mechanisms that govern both their range and distribution occur largely at atomic (e.g., radioactive decay), molecular (e.g., the influence of mass on the free energy of atomic bonds) and pore (e.g., diffusive transport to reactive surface) scales. In contrast, the vast majority of isotope ratio measurements are based on sample sizes that aggregate multiple pathways, species and compositions. Inferring process from such macro-scale observations therefore requires unraveling the relative contribution of a variety of potential mechanisms. In effect, the use of isotopes as proxies to infer a specific parameter, such as temperature (Urey 1947) or residence time (Kaufman and Libby 1954), carries the implicit requirement that one mechanism is the primary influence on the measured isotopic composition of the composite sample.

In the present chapter, we consider a wide variety of macro-scale observations of isotope partitioning across fluid–solid phase boundaries. For this purpose we define the continuum scale as a representation in which interfaces are averaged over elementary volumes, as opposed to the pore scale in which these interfaces are explicitly resolved. Throughout this review it will be demonstrated that observations of isotope partitioning across fluid–solid boundaries require some representation of the isotopic composition of the solid surface and surrounding fluid distinct from ‘bulk’ or ‘well mixed’ reservoirs. For example, this distinction is necessary in order to (1) quantify the partitioning of radioactive and radiogenic species, (2) describe transport limitations that may impact the macroscopic partitioning of isotope ratios, (3) explain observations of transient fractionation due to dissolution through preferential release at the solid surface, and (4) parameterize apparently variable fractionation factors during precipitation. The ability to describe isotopic partitioning specific to the phase boundary then influences the accuracy of simulations from highly variable, field scale systems (e.g., Druhan et al. 2013) to highly controlled laboratory experiments (e.g., Tang et al. (2008)). These observations imply that quantifying the composition of solids and fluids as averages across a given representative volume carries some inherent loss of information that isotopes appear to be sensitive to. This conclusion leads us to consider mechanistic descriptions of isotope partitioning that could be significantly improved by modeling approaches that are capable of resolving spatial zoning within individual solids (e.g., Li et al. 2006; Tartakovsky et al. 2008; Molins et al. 2012; Molins 2015, this volume; Yoon et al. 2015, this volume).

The suite of experimental data and quantitative approaches described herein support a conceptual framework in which macroscopic observables, such as an apparent fractionation factor, are the emergent result of multiple interacting processes that are strongly influenced by the physical and chemical characteristics of the fluid–solid boundary. In this sense we consider the pore scale as a characteristic length over which the unique physical and chemical properties of the phase interface may be described as distinct from bulk or aggregate values. This definition of the pore scale provides a critical reference frame over which molecular scale mechanisms combine to yield the macroscopic observables (Steefel et al. 2013).

A conceptual model of isotope partitioning at the pore scale

The flux of solutes in porous media is influenced by the development of spatial and temporal gradients resulting from the combined effects of transport and reactivity. How the pore-scale nature of these processes emerges into continuum scale observations of fluid–solid interaction is yet unresolved and requires multiscale (e.g., fractal) upscaling methods that remain a challenge. For example, Darcy scale fluid flow may occur along a fixed gradient, but the size, shape and distribution of individual solid grains is such that at the pore scale velocity vectors vary in both magnitude and direction. As a result, the principle mechanism of solute transport in some areas is advective, while in others it is diffusive. Mixing between these distinct regimes is approximated at the continuum scale by either a heterogeneous conductivity field (Li et al. 2010; Sudicky et al. 2010), a hydrodynamic dispersion coefficient (Gelhar and Axness 1983; Steefel and Maher 2009) or a non-uniform fluid travel time distribution (Maloszewski and Zuber 1982; Bellin and Tonina 2007).

Reactions that occur between fluid and solid phases are influenced by these local pore-scale transport regimes. For example, where transport of solutes between a solid surface and a surrounding fluid is accomplished primarily by diffusion, the rate of reaction across that interface may be limited by either the delivery of solute to the reactive surface or the approach to equilibrium. When advection is dominant, the same process may be governed by the relative rates of reactivity and flow (Rolle et al. 2009; Maher 2010; Hochstetler et al. 2013). The factors influencing isotopic partitioning are then equally variable across the pore scale (Fig. 1). In areas of the domain governed by diffusion, the partitioning of isotopes may reflect a transport limitation or diffusive fractionation, whereas in areas where flow is relatively fast, the isotopic composition is likely governed by the partitioning associated with a reaction (Lemarchand et al. 2004; Fantle and DePaolo 2007; DePaolo 2011). The accumulation of newly formed solids is also anticipated to vary with the local transport regime, which in turn exerts an influence on the isotopic composition of the fluid–solid interface.

This relationship between the dominant mechanisms of transport and fractionation suggests that over a representative volume of porous media a variety of fractionation factors may be observed. For example a fractionation factor characteristic of diffusive transport between the solid–fluid surface and the well-mixed fluid, or a fractionation factor associated with the difference in isotope composition between the fluid and solid at each reactive surface. The continuum-scale or observable fractionation factor over the timescale of the reaction (e.g., seconds to years) is then subject to the location and volume of material sampled. Over much longer periods of time, processes like recrystallization and solid-state diffusion homogenize spatial zoning within the solid. This implies that transient isotopic partitioning should be observable over significant periods of time, and that the observed macroscopic fractionation factor is potentially (1) a combination of multiple, distinct mechanisms and (2) variable as a result of parameters such as saturation state and flow rate.

From this perspective, variability in the magnitude of an observed fractionation factor is in some ways analogous to the discrepancy in rate constants observed across natural systems (Malmstrom et al. 2000; White and Brantley 2003; Maher et al. 2006b). This wide range in what should in principle be a fixed value has been attributed to the relative influence of a variety of processes, such as changes in reactive surface area (VanCappellen 1996), transport limitation (Steefel and Lichtner 1998; Maher 2011; Li et al. 2014); and even climate variations (Kump et al. 2000; Maher and Chamberlain 2014). Similarly, the apparent isotopic partitioning observed in flux-weighted or homogenized natural samples is often distinct from that obtained under controlled experimental conditions. These effects have been noted in a variety of contexts, including oxygen and nitrogen isotope fractionation in marine sediments (Brandes and Devol 1997), selenium isotope fractionation in wetlands (Clark and Johnson 2008), compound-specific stable isotope analysis (CSIA) of organics (Van Breukelen 2007) and chromium isotope fractionation in a contaminated aquifer (Berna et al. 2010).

To the extent that pore structure influences the distribution of isotope ratios in reactive systems, using volume-averaged sample measurements to quantify related parameters, such as reaction progress or mixing, can result in significant uncertainty. For example, both analytical and numerical solutions have demonstrated that neglecting the effects of hydrodynamic dispersion on observed stable isotope ratios can lead to an underestimation of reactivity in through-flowing systems (Abe and Hunkeler 2006; Van Breukelen and Prommer 2008). As noted above, dispersion is one method of parameterizing mixing at the continuum scale, and thus approximating the contribution of distinct fluid isotopic compositions at the pore scale (Fig. 1). Improved accuracy should then be achievable by describing distributed pore-scale isotopic compositions as the result of a variety of mechanisms, e.g., multiple fractionating reactions, a difference in the diffusion coefficient of the isotopologues of a compound or the dampening of a kinetic fractionation as a result of transport limitation. In the current article, we consider observations of isotopic partitioning across fluid–solid boundaries associated with a variety of fractionating mechanisms. In reviewing these examples we emphasize the aspects of each process that lead to pore-scale heterogeneity and thus a disconnect in behavior between the scale of mechanism and the scale of observation. These categories are by no means an exhaustive list of all processes that result in isotopic partitioning across reactive interfaces, but serve as examples in which interpreting the response of continuum-scale isotopic values to variations in external parameters is improved by consideration of pore-scale isotopic distributions.

Organization of article

The structure of the article is broadly divided into three sections. First, we provide a brief explanation of the notation used to describe isotope partitioning, with an accompanying discussion of primary mechanisms and models for fractionation. Second, we describe experimental observations for four examples of isotopic partitioning across fluid–mineral interfaces: α-recoil, diffusion, dissolution, and precipitation. Across this wide range of processes, a common observation is that the interface between fluid and solid phases (1) governs material transfer between reservoirs and (2) displays isotopic compositions distinct from bulk values. Associated models for the mass balance at phase boundaries are described with particular emphasis on the mechanisms necessary in order to quantify macroscale observations. This leads to the second section in which current modeling techniques for the description of transient isotopic final and the development of zoned mineral grains are discussed, concluding with a novel application of pore-scale modeling techniques to describe the distribution of stable isotopes across a fluid–mineral boundary as an initially oversaturated system establishes equilibrium. Throughout this review the principal intent is to describe experimental and modeling studies of isotopic partitioning with reference to the ways in which the composition and gradients across fluid–mineral boundaries is distinct from bulk-averaged values and uniquely influences continuum-scale observations.

The exchange of material across fluid–mineral surfaces influences both radiogenic and stable isotope distributions. Unlike (most) stable isotope fractionation, variations in radiogenic isotopes are mass-independent and arise due to the radioactive decay of a parent nuclide to an intermediate radioactive nuclide or a stable, radiogenic nuclide. The uranium decay series (Fig. 2) illustrates the relationship amongst parent and daughter isotopes, and the terminology is applicable to other radioactive-radiogenic systems.

The 238U isotope is radioactive and for N nuclei at time t in a closed system the number of decaying nuclei dN in the time interval dt follows a homogeneous Poisson statistical law and is proportional to N:


where lambda (λ) is the decay constant. The expression on the right side is a probability per unit time of a radioactive decay event. Integration of Equation (1) yields an exponential equation for radioactive decay:


where the integration constant (N0) is equal to the initial number of nuclei at t = 0. The decay constant is related to the nuclide half-life


For intermediate daughter products such as 234U and 226Ra, isotopic abundance is commonly discussed in terms of activity (A) where:


and thus the initial activity at t = 0 is A0 = λN0. In a closed system that starts with no daughter isotopes, for example the 238U- 234U system, the initial activity of 234U is zero and it will increase until A238U = A234U. At this point the supply of 234U from 238U decay is balanced by the decay of 234U to 230Th, representing a steady state where the ratio of λ limits the maximum possible intermediate daughter activity. When A1 = A2 this is a special condition known as secular equilibrium. In the remainder of the chapter we will use the notation 234U/238UAR = A234U/A238U where the subscript AR denotes the activity ratio.

Finally, for stable radiogenic daughter isotopes (e.g., 206Pb, in the 238U series), accumulation of isotope N2 is related to the initial abundance of the parent isotope (N0) and time:


such that as t goes to infinity N2 goes to N0. Krane (1988) presents an exhaustive discussion of radioactive decay including considerations for branched decay, and very short-lived intermediate daughter isotopes. The subsequent discussion will use the 238U–234U system to consider the influence of pore-scale effects in quantifying α-recoil damage and associated alterations to solid surfaces and solution chemistry.

As noted above, there can be ambiguity in the notation used to describe the partitioning of stable isotopes. In the subsequent text we primarily adhere to the guidelines put forth by Coplen (2011). The isotope ratio (R) of a particular reservoir is defined as


where N is the number of atoms of iE and jE, the isotopes i and j of the element E in substance P. This value is commonly reported relative that of a standard ratio (std), referred to as the delta value (δ), and defined as


This delta value is nondimensional, but is commonly multiplied by 1000 to report values as parts per thousand or per mil (‰). The difference between the isotope ratios (Δ) of two compounds or phases (P and Q) is then


The apparent or net isotopic fractionation factor does not make use of delta notation, but is defined as


Conversion between ΔiEP/Q and αiEP/Q may be approximated as ΔiEP/Q ≈ ln αiEP/Q. From these relationships it is noted that both the difference between the isotope ratios (Eqn. 8) and the apparent fractionation factor (Eqn. 9) between two reservoirs are obtained from direct measurement regardless of the variety of reaction pathways necessary to produce them, and are therefore a function of the representative volume of the sample in all but the most simplified systems. In contrast, a fractionation factor associated with a specific reaction pathway may be derived in reference to the rate law for that reaction. These types of fractionating processes are discussed in the subsequent section. Finally, an isotopic mole fraction (X) may be defined as the ratio of the amount of a particular isotope in a given species, compound, or reservoir divided by the total amount of that element in the same group


This value is used in the derivation of isotope-specific reversible rate expressions that involve a solid phase.

A note on fractionation

In this chapter, stable isotope fractionation will be discussed in terms of equilibrium and kinetic processes. For this purpose chemical equilibrium is described as a dynamic state that occurs when two elementary reactions, the forward reaction from reactants to products and the backward reaction from products to reactants, are in balance. For example, the hydroxylation of dissolved CO2:


which is at equilibrium when in a closed system the distribution of species is invariant in time. The ratio of the elementary forward and backward rates is then termed the equilibrium constant (Keq). Isotopic equilibrium may be described similarly, for example considering carbon isotope equilibrium between the carbon-bearing species in Equation (11):


leading to a separate equilibrium constant to account for this isotopic partitioning. A comparable relationship could be written for the partitioning of oxygen isotopes between CO2 and HCO3, or between HCO3 and OH. While all of these equilibria are described using the same law of mass action, it is important to note that they are not predicated on one another. In other words, these descriptions allow for the establishment of chemical equilibrium without necessarily requiring isotopic equilibrium. For reactions involving simple stoichiometric relationships, where one atom of the element of interest occurs in all reactant and product species, the isotopic equilibrium constant is equivalent to the fractionation factor (Eqn. 9); however, the relationship is often more complex (Schauble 2004). In general, the equilibrium fractionation factor is temperature dependent and larger for low mass elements and for isotopes of the same element that have large differences in mass. Typically, the partitioning of stable isotopes between two phases at equilibrium preferentially incorporates the heavier isotope in the phase with lower bond energy.

An imbalance between the forward and backward rates leads to a net accumulation of either the product or reactant, and the rate of this mass transfer is described by kinetics. Many reactions take place under conditions in which the reverse reaction is in some way prohibited, or the system is very far from equilibrium, such that isotopic partitioning is entirely kinetic. A kinetic fractionation factor is expressed as the ratio of the isotopic composition of the instantaneously generated product species (Pinst) and the residual reactant (Q) through a single reaction pathway:


This αkin can be related to the kinetic rate constant of the reaction (k), depending on the order of the reaction. For example, a first order dependence on concentration (e.g., dP/dt = kQ), leads to an expression for αkin = ik/jk (Mariotti et al. 1981). From an observational perspective it is often difficult to categorically identify equilibrium vs. kinetic effects on isotope partitioning. In low-temperature systems equilibration can be extremely slow, and, in addition, open system conditions may support the establishment of a steady state that appears balanced but is not necessarily equilibrated. In this sense the distinction between a specific state that is dynamic equilibrium, and an observable net rate of precipitation or dissolution does not imply an exclusive influence of equilibrium vs. kinetic fractionation. The approach towards an equilibrated system implies instead a continuum between pure kinetic fractionation and equilibrium fractionation. Later sections of this chapter explore the consequences of such a model and the extent to which pore scale treatments are capable of improving upon current approaches.

Under certain conditions, it is possible to develop theoretical models to predict the evolution of the isotopic composition of products and reactants. These models offer the advantage that they are simple to use, but are generally limited by particular assumptions and do not reflect the complexity of the reaction pathway or the relationship between transport and reaction in the isotopic mass balance description. Rayleigh fractionation, for instance, assumes an open system distillation process, where the reactant is progressively consumed, such that the isotope ratio of the reactant follows (Rayleigh 1902):


where the isotope ratio (Eqn. 6) is equal to the product of the initial isotope ratio (R0) and the fraction of reactant remaining relative to the initial concentration (f) raised to the power (α−1). This relationship requires a constant fractionation factor α, and produces an exponential relationship between reaction progress and isotopic partitioning. This model is only strictly intended for systems in which the reactant is continually supplied (f cannot go to zero) and the product of the reaction is instantaneously removed or segregated from the reactive system (Criss 1999). In practice the Rayleigh model is used in a wide variety of systems because it offers a simple relationship between reaction progress and fractionation without requiring knowledge of the reaction pathway or transport mechanisms. As a result, several studies have pointed out limitations to the Rayleigh model in application to hydrogeochemical systems (Brandes and Devol 1997; Abe and Hunkeler 2006; Van Breukelen and Prommer 2008). A goal of the current chapter is to describe a variety of common fractionating mechanisms that often undermine the assumptions of such simplified models and promote new quantitative methods for explicit treatment of reactivity and transport in the description of isotope partitioning.

Alpha recoil

Prior to the development of modern radioactive decay counting and mass spectrometry techniques it was widely assumed that the 234U/238UAR could not deviate from secular equilibrium. This was assumed to be the case because the small mass difference (~1.7%) of the two isotopes would not produce stable isotope effects like those observed for hydrogen and oxygen. Careful study of natural rocks and minerals in the 1950’s, however, revealed variations from secular equilibrium (Chalov 1959). Since the chemical behavior of 234U and 238U should be nearly identical, and yet variations in 234U/238UAR can exceed 500%, researchers suggested that the energy associated with 238U decay could directly eject the 234Th daughter isotope from the mineral surface (10’s of nm) or damage the mineral lattice allowing preferential leaching of the daughter isotope (Rosholt et al. 1963; Kigoshi 1971).

Kigoshi (1971) carried out a pioneering study where fine-grained zircon sand was suspended in dilute aqueous solutions and the addition of 234Th and 234U to the solution was quantified. The 234Th activities of the solutions were consistent with predicted α-recoil injection to the solution based upon a spherical grain model and an α-recoil distance of 55 nm. Earlier inferences (Turkowsky 1969) and later laboratory studies (Fleischer 1988) suggest that the recoil distance is closer to 30–40 nm. Both the 234Th and 234U in the fluid increased with time at a rate greater than the increase in 238U, demonstrating that 234Th has a recoil distance of tens of nm and can be ejected from the mineral structure or preferentially leached from lattice defects (Kigoshi 1971; Fleischer and Raabe 1978).

Kigoshi (1971) calculated the expected rate of addition of 234Th to a fluid (Q) due to α-recoil ejection from a mineral grain based on:


where L is the α-recoil range, S is the surface area, ρ is the solid density, and N238 and λ238 are the number of 238U atoms per gram of solid and the decay constant, respectively. The activity of 234Th in solution is a production–decay equation (combining Eqns. 4 and 5) of the form:


All of the terms on the right side of the equation are known or can be measured independently, however, in natural systems the grain surface area (S) and its influence in directing α-recoil to the fluid phase is arguably the most important variable. For example, Kigoshi (1971) assumed spherical grain geometry for their 1–10 μm diameter zircon sand and then fit the measured 234Th activities to calculate a characteristic recoil distance of 55 nm. Once L is known from similar experiments or determined by other methods, the probability of an α-particle (fα) being ejected from a mineral grain assuming a spherical geometry can be approximated:


Subsequent papers have explored the evolution of the solid phase and the effects of non-ideal grain geometry on the production of daughter products to the fluid phase (Kigoshi 1971; Fleischer and Raabe 1978; Fleischer 1980; 1988; Maher et al. 2004; DePaolo et al. 2006, 2012; Lee et al. 2010; Handley et al. 2013). A general conclusion from these studies, particularly those of Lee et al. (2010) and Handley et al. (2013), is that the chemical treatment of sediments and the model of grain surface structure need to be carefully considered and normalized across multiple laboratories in order for results to be broadly useful. Similar observations are noted, with important caveats for the differing chemical behavior, in the elements Th, Ra and Rn (Torgersen 1980; Semkow 1990; Sun and Semkow 1998; Porcelli and Swarzenski 2003). For example Ra is strongly adsorbed to mineral surfaces at low ionic strength but soluble at high ionic strength, thus the solution chemistry is critical for interpreting Ra activities (Moore 1976).

The general equations describing the evolution of uranium series isotopic ratios in pore fluids and minerals are presented by Ku et al. 1992; Porcelli et al. 1997; Henderson et al. 2001; Tricca et al. 2001; Porcelli and Swarzenski 2003; Maher et al. 2004. The evolution of pore fluid composition is related to primary mineral dissolution, secondary mineral precipitation reactions, α-recoil and preferential leaching of daughter isotopes to the pore fluid, sorption–desorption reactions and diffusion–advection of fluid in and out of a pore. Alpha-recoil loss, preferential leaching near the mineral surface, and solid-state diffusion, will all affect the solid composition with respect to time. However, these processes operate at different timescales, allowing mineral grains to evolve distinct 234U/238UAR domains. There is considerable discussion in the uranium series literature about the relative roles of preferential leaching and direct α-recoil leading to daughter isotope accumulation in the pore fluid (Rosholt et al. 1963; Vigier et al. 2005; DePaolo et al. 2006; Dosseto et al. 2006). Preferential leaching may occur due to mineral lattice damage associated with 234Th recoil (Rosholt et al. 1963) or due to preferential oxidation of the 234U in damaged mineral lattice by aqueous fluids (Kolodny and Kaplan 1970). Regil et al. (1989), Roessler (1983, 1989) and Adloff and Roessler (1991) present detailed models of 234U oxidation due to α-recoil.

It is difficult to ascertain the exact mechanism that transfers 234U preferentially from the solid phase to the pore fluid, but this process has important implications for the interpretation of uranium series isotopes in mineral–fluid systems. DePaolo et al. (2006) suggest that based on fine-grained alluvial sediments, the leaching depth into grains is not appreciably greater than the recoil distance. For the purposes of this discussion we proceed under the assumption that direct α-recoil is the primary mechanism of 234U transfer, but acknowledge this may be unwarranted, particularly in coarser-grained or uranium rich minerals.

The mass conservation equations presented below are applied to the 234U and 238U isotopes and illustrated schematically for distinct pore fluid, solid surface, and solid interior compositions in Figure 3. This formulation can be applied to other intermediate daughter products with additional consideration for differing chemical behavior (e.g., strong sorption/ secondary mineral partitioning of thorium and radium under certain conditions). At the scale of a single pore surrounded by mineral grain surfaces, the 234U/238UAR of the fluid will evolve based on the dissolution–precipitation reactions and α-recoil ejection to the fluid.

The 234U/238UAR evolution of mineral grains with time is:


It is apparent from Equation (18) that probability of α-particle ejection (fα) and time are the critical variables that affect the activity ratio of the solid grains (DePaolo et al. 2006). Estimating fα requires knowledge of the mineral volume surface area and recoil distance. DePaolo et al. (2006, 2012) calculated fα as a function of grain diameter and shape (surface area), demonstrating, for example, that mineral grains of 10-μm diameter could have greater than a factor of 10 variability in fα(Fig. 4).

A subsequent study by Bourdon et al. 2009 suggests employing a combined BET sorption surface area measurements and a fractal dimension model to overcome the uncertainty of mineral surface areas in natural sediments.


where SBET is the measured surface area, a is the molecule diameter for the adsorbate and D is the fractal dimension of the grain surface, which is measured independently. Non-ideal (spherical) surface area is particularly important at the pore scale, where roughness (Lee et al. 2010) can increase the surface area of a mineral grain by a factor of 17. The sorption data of Bourdon et al. (2009) have been successfully applied to studies of soils (Oster et al. 2012) and to date ice using trapped fine-grained sediments (Aciego et al. 2011).

The evolution of 234U/238UAR of pore fluids is described by Equation (20):


where C is the concentration and A is the 234U/238UAR of the solid (s) and fluid (f), Rd is the rate of dissolution in inverse time, Rf is the retardation factor accounting for sorbed uranium and Ms is the solid/fluid volume ratio (Tricca et al. 2001; Maher et al. 2004).

The conservation equation can be rearranged to yield the following approximation for the pore fluid at steady state (Af(ss)):


This relationship can then be manipulated to use the activity ratios to calculate the solid dissolution rate:


Maher et al. (2004) used this relationship to calculate the dissolution rates of deep-sea sediments over timescales longer than those accessible by laboratory experiments.

The models described above are a small sample of the work on uranium series isotopes in hydrogeologic systems. At the pore scale additional complexities may be considered, particularly when the steady-state assumption is not valid. For example, rapid changes in the dissolution rate could dissolve the 234U depleted layer of mineral grains, lowering the 234U/238UAR in the pore fluid. Additionally, where mineral grains are in direct contact with one another the ejected daughter isotope may be implanted into an adjacent mineral grain and not accumulated in the pore fluid. The net effect of significant implantation would be to underestimate the sediment dissolution rate (Eqn. 22) or the age of fine-grained sediments (DePaolo et al. 2006). In general the aforementioned studies demonstrate that uranium series disequilibria is largely controlled by the surface structure of individual grains that affect the isotopic composition of fluids considered at the continuum scale. Uncertainty in the relative contributions of α-recoil and preferential leaching to isotopic fractionation in the uranium series and precise measurements of mineral surface area and recoil tracks currently limit the fidelity of existing models. Future directions in this area may include incorporating daughter isotopes with different recoil distances such as 230Th and 226Ra to better constrain the fα parameter and more sophisticated imaging techniques to quantify mineral surface structure. Additionally, precise measurements of the α-recoil distance (L) and the relationship between surface area and fα will be necessary to improve the application of pore scale 234U/238UAR to diverse problems such as sediment dating, mineral dissolution rates and vadose zone transport.

Diffusive fractionation

The change in concentration (C) of a solute due to a net divergence of the flux (J) in an arbitrary volume


is often approximated with the linear transport theory that relates the diffusive flux to the gradient in chemical potential or concentration, i.e., Fick’s first law:


Here the transport coefficient D is the molecular diffusivity of the solute and has units of length squared per time.

Under ideal conditions, kinetic theory provides an expression for the diffusion coefficient of molecules in a gas


as the product of the mean free path (l) between collisions, and the square root of the product of the gas constant (R) and temperature (T) divided by the mass (m) of the molecule or atom. This relationship suggests that the mass of the particle influences the transport coefficient and by extension the flux to the reactive surface. Therefore the lighter the particle the larger the diffusion coefficient, which leads to the expectation that if two particles differ in their masses (m1 and m2), then for the same conditions the difference in their diffusion coefficients would be


However, this relationship is specific to an ideal gas, and in liquids the factors that contribute to the value of an individual species diffusion coefficient become much more complex (Moller et al. 2005; Yamaguchi et al. 2005). For instance, in electrolyte solutions, the influence of electrochemical forces on the interaction of diffusing species cannot be neglected, and thus requires accounting for cross-diagonal terms representing the influence of one species’ or isotopic activity gradient or electrochemical potential gradient on another. This effect may be quantified in the absence of charged surfaces by including treatment of electrochemical migration using the Nernst–Plank equation (Steefel and Maher 2009; Steefel et al. 2015; Rasouli et al. 2015), while treatment of surface charge requires further inclusion of a complete electrical double layer model (Tournassat and Steefel 2015, this volume).

In the context of stable isotope fractionation in aqueous solutions, this complexity has only recently been considered. Richter et al. (2006) conducted a series of experiments to determine the fractionation of isotopes associated with diffusion in comparison with previously estimated values. Their experimental design allowed a small volume of aqueous solution containing a dissolved salt (V1) to be suspended inside of a much larger volume of dilute fluid (V2). The length and aperture of the tube connecting the two reservoirs was selected such that the solution inside of V1 would remain well mixed over the timescale necessary for dissolved solutes to diffuse into V2. Richter et al. (2006) considered departure from Equation. (26) by recasting the power as a parameter (β), which they then fit to their experimental data.


The results of their study determined β values less than 0.5, but greater than 0.0, thus providing evidence that kinetic fractionation due to diffusion could contribute to observed variations in stable isotope ratios in aqueous solutions. However, their method required the implicit assumption that an inverse power law (Eqn. 27) is an accurate description of the dependence of a diffusion coefficient on the mass of the solute. Bourg and Sposito (2007) sought to directly test this assumption by utilizing molecular dynamics (MD) simulations. The use of a numerical method allowed them to obtain the diffusion coefficient for more than two isotopes of the same species, and thus test the generality of the β values reported by Richter et al. (2006). Their results showed good agreement with the Richter et al. (2006) values, demonstrating that the value of β for a given mass pair in aqueous solutions is typically less the value of 0.5 predicted by kinetic theory for an ideal gas (Eqn. 26). Following on the agreement of these results, Bourg et al. (2010) extended both the experimental and numerical results to include additional ions. The results of these studies are summarized in Table 1.

Two key aspects of these observations are that the monovalent solutes show a larger contrast in the diffusivity of their isotopes than divalent solutes, and this difference does not follow a direct correlation with the ratio of the masses. Bourg et al. (2010) demonstrated that this relationship occurs as a result of (1) the size of the solute radius, and (2) the strength of attractive interactions between the solute and solvating water (effective mass of the diffusing species), both of which influence coupling of motion between the solute and solvent. Using the residence time of water (τ) in the first solvation shell of the ion obtained from the MD simulations as a proxy for this coupling, Bourg et al. (2010) showed a clear correlation between the inverse of τ and β (Fig. 5).

The aforementioned experimental observations and MD simulations indicate that pore scale gradients in both concentration and isotopic abundance between a reactive surface and the surrounding fluid (Fig. 1) are influenced by diffusion coefficients that differ as a result of the size and charge of ions and the difference in masses of their isotopes. The effects of ion size and charge are explicitly treated at the resolution of MD simulations, but must be approximated at larger scales. Therefore accurate representation of these diffusive effects on isotopic partitioning at the pore and continuum scale requires independent knowledge of the magnitude of fractionation associated with the isotopes of each ion of interest (Table 1). Furthermore, the extent to which observed differences in diffusivity relate to fractionation associated with dispersion at larger scales is unknown.

From the results of Richter et al. (2006), Bourg and Sposito (2007) and Bourg et al. (2010), the value of β relating the ratio of the masses to the ratio of the diffusion coefficients is highly variable among ions and solution compositions. In principle, the ratio of the diffusion coefficients could be directly obtained for two isotopologues, specific to a given solution composition, through a simple diffusion experiment. For example, implementing a generic tracer in a discretized domain using the CrunchFlow (Steefel et al. 2015) reactive transport code such that at the start of the simulation an elevated concentration exists in the center of a closed system, one may obtain the spatial partitioning of the tracer as it diffuses through the system (Fig. 6A). If this total tracer is broken into two components, where one ‘isotope’ composes the majority of the tracer concentration and the other is a trace amount, then a slight difference in the diffusion coefficients of the two species results in a spatial profile of their ratio through time (Fig. 6B). We emphasize that this simple modeling example demonstrates a proof of concept wherein an experimental system could provide a time series of the tracer concentration (Fig. 6C) and isotope ratio (Fig. 6D) at a fixed observation point. For any system in which variation in the isotope ratios are large enough to be detected beyond measurement error, a time series dataset of this nature could be used to obtain unique values for the diffusion coefficients. If the difference in the masses of the isotopologues is known, these values could then be related to obtain direct estimates of β for a given system.

The sensitivity of isotope ratio measurements are such that diffusive fractionation may influence the interpretation of observed fractionation, but these effects have only recently begun to be considered in hydrologic studies. Analytical and numerical models have demonstrated that even in the absence of true diffusive fractionation, neglecting the effects of diffusion and dispersion leads to a lower overall estimate of reaction progress based on stable isotope ratios than when these factors are considered (Abe and Hunkeler 2006; Van Breukelen and Prommer 2008). These results have lead to a recent focus on the effects of dispersive mixing on stable isotope ratios, with particular emphasis on compound-specific stable isotope labels used to track contaminant degradation. Hydrodynamic dispersion is a fundamentally distinct process from molecular diffusion, arising from the fluctuations in velocity within and among connected pores during flow. Under certain assumptions, dispersion has been shown to display a diffusive behavior. In practice, most efforts to directly simulate fractionation due to dispersion in through-flowing systems have made use of estimated relationships between the ratio of the masses and the ratio of the dispersion coefficients. Dispersive isotope fractionation has been simulated at the continuum scale for both hydrogen and carbon isotopologues of solutes during flow through heterogeneous porous media using a variety of approximations for the difference in isotope-specific dispersivities. LaBolle et al. (2008) calculated isotope-specific values of dispersivity by substituting the reduced masses of the solutes into Equation (26). Rolle et al. (2009), Eckert et al. (2012) and Van Breukelen and Rolle (2012) used an empirical correlation from Worch (1993) close to Equation (26) but with a power of 0.53. Thus far it is not evident what, if any relationship may exist between diffusive and dispersive fractionation, or even if the power law relationship between the ratio of the masses and ratio of the diffusion coefficients validated by Bourg and Sposito (2007) holds true for the ratio of two isotopic dispersion coefficients.


Fractionation directly associated with the dissolution of material from the solid phase is often considered negligible over the reactive timescales represented by many natural samples. This assumption is justified in that many of the mechanisms thought to contribute to isotopic partitioning, such as transport (discussed above), and changes in solvation (Hofmann et al. 2012), are restricted within a solid phase. Furthermore, where apparent shifts in isotopic ratio are observed during net dissolution of natural samples, it is difficult to distinguish preferential mobilization of individual, isotopically distinct mineral phases from true fractionation (Ryu et al. 2011). However, a wide variety of studies now appear to demonstrate observable fractionation during dissolution, usually during the initial stages. These observations are significant in that they implicitly require that the solid become isotopically differentiated as material is removed from the reactive surface. This implies the development of a surface that is compositionally distinct from the interior of the grain, a process that is not readily compatible with bulk or volume fraction representations of solids.

The earliest observations of fractionation during dissolution involved dissimilatory reduction of ferric iron by anaerobes (Beard et al. 1999, 2003; Brantley et al. 2001, 2004), which naturally led to the issue of whether observed iron isotope partitioning was a vital effect or if it could be reproduced in abiotic systems. From this starting point, a wide range of experimental conditions, including both pure mineral phases and whole rocks, as well as a variety of stable isotope ratios, have shown evidence for partitioning of isotopes during net dissolution (Table 2). The compilation in Table 2 is limited to low temperature and ambient pressure experimental conditions in which efforts have been made to minimize fractionation of fluid phase isotopes due to additional effects, such as secondary mineral precipitation.

Despite wide variability in experimental conditions, a common observation is that during dissolution the isotopic composition of the fluid phase tends to be lighter than that of the initial solid, i.e. the light isotope tends to be preferentially dissolved. Furthermore, in studies where the isotopic value of the fluid is monitored as a function of time, the magnitude of this fractionation appears to vary. When a solid is introduced to an aqueous solution such that the system is undersaturated, the fractionation between the fluid and solid phase tends to exhibit a maximum value relatively soon after the experiment is initiated (Fig. 7). For example, the dissolution of goethite at pH 3 (Wiederhold et al. 2006) promoted by ligand complexation leads to a maximum fractionation between dissolved Fe and the bulk solid of −1.83‰, but this difference lessens as dissolution continues (Fig. 7A). A similar temporal trend is noted in iron isotopes during proton-promoted dissolution of granite (Kiczka et al. 2010)( Fig. 7b), in zinc isotopes during both proton and ligand promoted dissolution of biotite (Weiss et al. 2014) and in silicon isotopes during proton promoted dissolution of Hawaiian basalt (Ziegler et al. 2005). At higher temperatures and pressures, the same trend is noted in lithium isotopes during the dissolution of basalt (Verney-Carron et al. 2011) and in magnesium isotopes during the dissolution of basalt and forsterite (Wimpenny et al. 2010), but magnesium isotopes show the opposite behavior during magnesite dissolution (Pearce et al. 2012). A low-temperature exception to this general trend may be the fractionation of copper isotopes associated with oxidative dissolution of sulfide minerals. Mathur et al. (2005) reported Δ65Cu values of +2.74‰ and +1.32‰ for the abiotic dissolution of chalcocite and chalcopyrite, respectively. They suggested this isotopic enrichment of the solute was associated with the accumulation of secondary phases, however, a subsequent study demonstrating +2.0‰ values for dissolution of whole rock containing chalcopyrite argued that secondary mineral formation was negligible (Fernandez and Borrok 2009).

The general observation of an initially negative Δfluid–bulksolid value that trends towards zero with further reaction progress is common to a variety of reactive pathways, minerals and isotope systems, and is observable even during dissolution of natural samples. For example, the biotite and chlorite enriched component of a granite sample dissolved in the presence of hydrochloric and oxalic acids demonstrated negative Δfluid–solid values that were most pronounced early in the reaction (Kiczka et al. 2010) (Fig. 7B). This study explored a range of dissolution rates by varying the potassium concentration, which served to slow down the rate of iron release in the HCl experiments, and observed a greater enrichment in the light isotopes of the fluid phase with decreasing dissolution rate. In contrast, addition of potassium to the oxalic acid experiments only made a small difference in the Fe release rates, and little difference in the Δfluid–bulksolid values were observed. Reaction of a biotite granite with the same acids yielded similar temporal trends, though the maximum Δfluid–bulksolid values showed lighter isotope composition for the fluid phase with more acidic solutions (Chapman et al. 2009). These observations suggest that the partitioning of stable isotopes may potentially be utilized as an indicator of the mechanisms of mineral weathering, leading to development of a variety of conceptual and quantitative models for fractionation during congruent dissolution.

The preferential removal of light isotopes during the initial stages of a dissolution reaction is commonly considered a kinetic fractionation associated with enrichment of the solid surface in the heavy isotope. Some studies have considered whether or not the initial stage in which this surface layer is established could be described as a distillation or Rayleigh process. Wiederhold et al. (2006) discussed this conceptual model in the context of iron isotopes and Weiss et al. (2014) in the context of zinc isotopes. Both studies point out that a Rayleigh model would carry the implicit assumption that the isotopes of the solid phase are well mixed, and thus homogeneously released. However, the development of either a depth into the solid that is isotopically zoned or a non-uniform distribution of reactive surface sites along a monoatomic surface layer both contradict this requirement. The consistent observation of a maximum Δfluid–bulksolid early in the reaction progress that becomes less depleted with time (Fig. 7) then requires consideration of an additional mechanism.

The development of an isotopically distinct solid surface has been modeled in terms of the propagation of a reacted surface layer into the unreacted solid. This is based on the idea of a preferentially reacted or ‘leached’ layer as commonly used to describe non-stoichiometric dissolution of silicates, e.g., Casey et al. (1993); Casey and Ludwig (1996); Hellmann et al. (1990); Oelkers (2001). Brantley et al. (2004) derived a model for isotopic fractionation during dissolution based on a transfer of material from the bulk solid through a layer altered by dissolution to the fluid. The depth of the altered surface layer (x), the flux of material from the bulk solid to this layer (F) and the rate of dissolution from the surface layer to the fluid (R = kC(t)) all contribute to the transient value of Δfluid–bulksolid (Fig. 8A). The model is given as a ratio of the rates of 54Fe and 56Fe entering the fluid phase:


where 54k and 56k are the rate constants for the dissolution of material from the surface layer to the fluid for 54Fe and 56Fe, respectively. Similarly 54Co and 56Co are the respective concentrations of the isotopes of iron in the bulk solid and t is the time since dissolution started. As discussed in Brantley et al. (2004), this equation effectively describes the period of time from the start of dissolution to the establishment of a steady state value for a given surface layer depth (x) (Fig. 8B and C). At very early time, relative rates of dissolution of the two isotopes are governed by the ratio of the rate constants, but as time progresses, the flux from the unreacted solid to the surface layer becomes more influential and eventually controls the composition of material transferred to the fluid. Conceptually, this transition occurs because at t = 0, the surface layer has not been formed yet and the fluid is in contact with the bulk solid, allowing for preferential release of one isotope vs. another as a function of the ratio of the rate constants. As the surface layer develops and eventually reaches steady state, the transfer of material to this layer from the underlying, unreacted solid becomes the governing rate. In the Brantley et al. (2004) model, this steady state value is achieved quickly if the thickness of the surface layer is large, if the migration of this reacted surface layer into the unreacted solid is slow, or if rate of dissolution from the surface into the fluid is fast.

This model approach clearly yields a temporal trend comparable to those observed during many of the dissolution experiments discussed in the observations section (Table 2). Furthermore, it suggests that the establishment of a steady state value is independent of the fraction of solid that is mobilized, but is rather a function of the establishment of a steady state leach layer. This implies that the trend towards a Δfluid–bulksolid of approximately 0.0‰ as dissolution progresses is not the result of a depletion in the abundance of the preferentially released isotope (i.e., a mass balance restriction) but rather a shift from the influence of dissolution (into the fluid) to the influence of material transfer (into the leach layer). There are some aspects of the assumptions built into this model that may be further developed. For example, the thickness of the leach layer is fixed and allowed to propagate inward throughout the simulation, but one may conceptualize a period at the very beginning of the model in which the distinct surface layer has yet to establish and thus x = 0. This would then require consideration of a period over which the surface layer grows into the solid and establishes a fixed depth. Furthermore, in Equation (28) an eventual Δfluid–bulksolid value of 0.0‰ results from a flux of Fe from the bulk solid to the surface layer (F) that is non-fractionating. This distinction between the mechanism through which iron is supplied to the leach layer and the mechanism by which it is dissolved into solution are potentially difficult to reconcile, and may benefit from consideration of solid state diffusion through the reacted zone (Verney-Carron et al. 2011).

Recent analytical developments have supported increasingly high-resolution observations of the structure and chemical composition of solid surfaces and a corresponding refinement in the conceptual model of reacted zones. Nano-scale resolution of reactive fronts are becoming more commonly described as a dissolution–reprecipitation process (Hellmann et al. 2012, 2015) leading to formation of a distinct, often amorphous phase at the solid surface. The extent to which the Brantley et al. (2004) model might be reworked to consider fractionating dissolution in systems where a reacted surface layer does not form has only recently begun to be addressed. Wiederhold et al. (2006) quantified stable isotope partitioning during dissolution of a solid phase that does not form an observable reacted surface layer using a model they describe as conceptually similar to the Brantley et al. (2004) approach. In their derivation for iron isotope partitioning during goethite dissolution, the reacted surface depth was replaced by a percentage of the total iron contained in the solid phase that was considered to exist within a reactive, monoatomic surface layer. The isotopes of iron could be preferentially removed from the reactive surface sites at the step-edge of this receding monoatomic layer, but are replenished by the isotopic composition of the bulk solid. Their derivation results in two free parameters, the fractionation associated with dissolution of Fe from the reactive surface layer and the percentage of the total solid phase Fe composing the reactive surface site pool. Importantly, they fit these two parameters and obtained a best estimate for the reactive surface site pool (2.4%) that matched quite closely with the value they independently quantified based on the surface site density, the measured surface area, and the molar mass of goethite (2.7%). This model was successfully extended by Kiczka et al. (2010) to consider the contribution of two independent types of reactive surface site in a more complex granite consisting of enriched biotite and chlorite fractions.

Macroscopic or Darcy-scale reactive transport models are now becoming a common means of quantifying stable isotope partitioning by separating the isotopes of interest into individual ‘species’ that follow parallel reaction networks. The fractionating reaction in these simulations is commonly a homogeneous aqueous process (Dale et al. 2009; Gibson et al. 2011; Druhan et al. 2012, 2014; Jamieson-Hanes et al. 2012; Wanner and Sonnenthal 2013). To date, models designed to simulate the influence of dissolution on the isotope ratio of fluids rely on the assumption that the solid-phase isotopic composition remains fixed over time. In doing so, dissolution of the solid acts as a source to the fluid phase that is isotopically invariant with time (Steefel et al. 2014a; Wanner et al. 2014). This approach is necessary given that at the continuum scale mineralogy is described as a volume fraction of the total solid in each grid cell, and thus does not readily support the development of an isotopically zoned mineral surface. Some approximation to this zoning might be achievable through the use of a nested-grid or multiple interacting continua approach (e.g., Xu and Pruess 2001), but this would severely restrict the size of the spatial domain and would require independent, a-priori knowledge of the length scales of zoning. Specifying a fixed isotopic composition of a solid phase also restricts the range of saturation states that can be accurately simulated. If a system were to become supersaturated with respect to a solid phase where the isotopic composition of that solid was pre-specified, accumulation of new precipitate would induce fractionation of the surrounding fluid. Stated differently, specifying a fixed solid phase isotopic composition effectively means that the rate law describing the dissolution of that solid is no longer inherently reversible.

The observed temporal trend in fractionation associated with dissolution indicates that these effects are likely negligible to long term processes, but potentially influential over shorter time periods. This implies that models for isotope fractionation applied to contaminant remobilization, microbially catalyzed oxidative dissolution and system response to short term environmental perturbations may benefit from an improved capability to simulate the isotopic partitioning characteristic of the initial stages of dissolution. Based on current modeling approaches, it appears that accounting for fractionation due to dissolution of a solid within a reactive transport framework requires a means of describing the development of a distinct solid surface that influences isotopic partitioning. This implies the need to resolve the geometry and surface of individual grains, i.e., a pore-scale model. Precipitation is less problematic, and has been more extensively described by Darcy-scale simulations, but as will be discussed in the next section this process is also influenced by the unique composition of the reactive interface.


In the previous section concerning studies of net dissolution, experimental conditions largely prohibited the accumulation of secondary minerals. Where secondary solids composed of the fractionating element were able to form, precipitation likely influenced the observed stable isotope ratio. For example Mathur et al. (2005) observed partitioning of copper stable isotopes during abiotic oxidative dissolution of chalcocite and chalcopyrite where the Δ56Cufluid–bulksolid was positive, in conjunction with mass balance calculations that suggested removal of copper from solution. Similarly, Wimpenny et al. (2010) report positive Δfluid–bulksolid values for lithium when the pH range allowed for concurrent secondary mineral precipitation.

The formation of a precipitate from solution is broadly recognized as a fractionating mechanism in a variety of stable isotope systems. For example, partitioning of the lighter 6Li isotope into the solid phase was predicted by ab initio calculations (Yamaji et al. 2001) and observed during carbonate precipitation (Marriott et al. 2004), clay formation (Williams and Hervig 2005; Vigier et al. 2008) and accumulation of secondary silicates (Wimpenny et al. 2010). These observations suggest that precipitation of secondary minerals contributes to enrichment in heavy lithium isotope in rivers relative to primary mineralogy in a wide variety of weathering environments (Vigier et al. 2009; Lemarchand et al. 2010; Millot et al. 2010; Dellinger et al. 2014). Similar partitioning of silicon isotopes associated with precipitation was predicted (Meheut et al. 2007, 2009) and inferred from observations of progressive basalt weathering (Ziegler et al. 2005; Georg et al. 2007) and riverine composition (Georg et al. 2006; Opfergelt et al. 2013). Recently, Geilert et al. (2014) demonstrated a temperature dependent preferential incorporation of the lighter isotope of silicon into amorphous silica, while Oelze et al. (2015) demonstrated a rate dependence of silicon fractionation during sequential precipitation of amorphous silica in the presence of dissolved aluminum. Johnson (2011) noted that in general the isotopic partitioning of redox sensitive elements are governed by pathways that alter the redox state of reacting species, and thus precipitation in the absence of electron transfer tends to impart comparatively minor fractionation in these systems. However, Skulan et al. (2002) demonstrated that under extremely rapid precipitation rates kinetic effects were observable in iron isotopes.

The Skulan et al. (2002) study offers a closed-system kinetic end-member example of the isotopic partitioning associated with precipitation. They utilized experimental conditions in which dissolved ferric iron was rapidly precipitated from solution to form hematite, resulting in approximately 80% removal of Fe(III)aq in one hour. These conditions lead to a large kinetic fractionation (Fig. 9A), which was observable in both the residual aqueous phase and the bulk hematite stable isotope ratios of iron. This observation was considered a kinetic effect due to the rapid rate of accumulation prohibiting isotopic equilibration between phases and because long-term exposure of hematite to solution resulted in small Δ56Fe values. Furthermore, enrichment in the heavy isotope composition of the fluid as the reaction progressed followed a Rayleigh model (Eqn. 14), adhering to an apparently constant fractionation factor of α = 0.99868 ± 12 and suggesting the influence of a single fractionating process. This distillation effect was pronounced enough that the shift in isotopic ratio as a function of reaction progress was also observable in the Δ56Fe of the bulk hematite. Although this cumulative value only increased by approximately 0.8 ‰ compared with a fluid enrichment of roughly 4.4 ‰, mass balance dictates that the isotope ratio of the hematite surface forming from the surrounding fluid varied as much as the fluid phase, leading to formation of a solid that was isotopically heterogeneous, or zoned. This zoning was evident when hematite crystals synthesized by this rapid process were subsequently dissolved in HCl (Fig. 9B). Despite a bulk hematite δ56Fe composition of −0.29 ±0.05 ‰, dissolution of the first 6.8 % of the solid yielded a fluid value of +1.16 ± 0.11 ‰, which subsequently decreased with continued dissolution to reach the value of the bulk solid.

Two aspects of the Skulan et al. (2002) results are relevant to the current discussion. First, a solid phase can preserve an isotopic signature associated with the kinetic effects of formation, and second in a closed system the kinetic signal imparted during formation of the solid phase varies as a function of the progress of the reaction, which in turn produces an isotopically heterogeneous or zoned precipitate. Dissolution of such a solid could then lead to the appearance of fractionation as an artifact of zoning. This concern was noted in several of the dissolution studies discussed in the previous section, for example Crosby et al. (2007) conducted separate partial dissolution measurements to ensure the homogeneity of their hematite and goethite samples. Based on the Skulan et al. (2002) results, if the light isotope is preferentially incorporated into the solid during precipitation, and this solid is subsequently dissolved, zoning would seem to counteract the general trend of light isotopes preferentially mobilizing in the fluid phase during dissolution (Table 2). The extent to which preferential mobilization of the light isotope influenced the values shown in Fig. 9B cannot be disentangled from the effect of the zoned hematite composition.

If the Skulan et al. (2002) results represent one extreme, in which the isotopic composition of a precipitate is entirely governed by unidirectional, kinetic fractionation, then the other end member is the formation of a solid at slow enough rates that isotopic equilibrium is maintained between phases. Many natural systems exist in a range of saturations close to chemical equilibrium (Helgeson et al. 1984) and solids are thus often anticipated to form over timescales sufficient to reflect isotopic equilibrium. For example, the use of stable isotope ratios as paleoproxies commonly requires that the timescales associated with formation of solids are sufficient to maintain isotopic equilibrium. This has been a particular issue in the interpretation of carbon and oxygen isotope compositions in carbonates, where multiple studies have demonstrated a combination of equilibrium and kinetic effects (Romanek et al. 1992; Kim and Oneil 1997; Chacko et al. 2001; Mickler et al. 2006; Dietzel et al. 2009; Gabitov et al. 2012; Gabitov 2013; Watkins et al. 2013). More recently, the same issues have been considered in the mid-mass stable isotope ratios of cations in carbonates (Fantle and Tipper 2014). For example both magnesium and calcium isotopes have been shown to partition during carbonate mineral formation, typically such that the fluid is isotopically heavy relative to the solid (Galy et al. 2002; Lemarchand et al. 2004; Gussone et al. 2005; Tang et al. 2008; Kisakuerek et al. 2009; Immenhauser et al. 2010; Rustad et al. 2010; Reynard et al. 2011; Wombacher et al. 2011; Mavromatis et al. 2013).

Observations of both light- and mid-mass stable isotope partitioning during carbonate precipitation seem to indicate a combination of kinetic and equilibrium effects. These observations imply that the mechanisms of fractionation should not be viewed as a binary system, but rather as a continuum between the two processes. As a result, a range of reaction rates should exist near equilibrium where the apparent or observed fractionation factor associated with formation of a solid phase is rate dependent (Fig. 10). Within a closed system, this effect leads to the accumulation of an isotopically heterogeneous solid as a function of both changes in reaction rate and distillation of the reactant reservoir. Within open systems, in the absence of any supply or transport limitations, the isotopic composition of the solid surface may still vary as a result of changes in the saturation state of the system. Thus experiments designed to quantify the magnitude of isotopic partitioning between a fluid and precipitating solid often utilize open systems in which the chemical composition of the fluid is continually titrated to maintain a fixed saturation state (Tang et al. 2008; Watkins et al. 2013), thereby ensuring that the accumulated solid phase is isotopically homogeneous. In natural systems, lack of such controls on saturation state or reaction rate lead to the formation of solids such that the observed fractionation factor appears to vary both spatially (e.g., Tipper et al. 2012) and temporally (e.g., Druhan et al. 2013). Therefore models that describe the exchange of material between a coevolving fluid and solid surface should be capable of describing this variability.

Several models have been proposed for the partitioning of stable isotopes during net precipitation under controlled experimental conditions. Tang et al. (2008) utilized the surface entrapment model (SEMO) developed by Watson and Liang (1995) and Watson (2004) to describe their variable fractionation factor as a function of precipitation rate (Fig. 10). The SEMO model was developed to describe zoning of trace elements and isotopes in solids that were thought to form under such slow growth rates that the fluid and solid surface must reflect equilibrium partitioning and an absence of aqueous transport limitations. In order to explain such observations, the SEMO model assumes that although the isotope ratios of the interior of the solid and the bulk fluid reflect the equilibrium fractionation factor, there exists a surface layer of the solid in which the rare isotope is depleted, with a maximum partitioning defined by a surface entrapment factor (F). Under growth conditions, this depleted surface layer may become ‘entrapped’ in the crystal beneath newly formed solid, leading to a disequilibrium fractionation. This effect is offset by re-equilibration due to diffusion within the crystal lattice. A key aspect of the model is the requirement that ion mobility, in this case cast as the solid-state diffusion coefficient governing re-equilibration (Watson 2004), increases from the crystal interior to a higher value at the solid surface. In application this approach was able to accurately reproduce the trend in fractionation as a function of precipitation rate for the Tang et al. (2008) dataset. However, the diffusion coefficient at the solid surface required to achieve this result was roughly 16 orders of magnitude higher than that estimated for trace element diffusion in calcite.

An alternate approach to the SEMO model was proposed by DePaolo (2011), which considers a kinetic effect in the fluid boundary layer surrounding the precipitating solid (Fantle and DePaolo 2007). This model represents a fundamental shift from both the SEMO approach and the models used to describe dissolution discussed in the previous section, in that there is no description of a depth or volume percentage considered to represent a reactive surface. Instead, the only aspect of the solid considered is the surface in contact with the surrounding fluid (Fig. 11A), which dictates the isotopic composition of the backward reaction. As discussed in the introduction it is the cumulative effect of an irreversible forward addition of material to the mineral surface and an irreversible backward removal of material from the mineral surface that influences the observable fractionation factor, despite the fact that the overall system is supersaturated and thus precipitate is accumulating. In principle, this creates a coupled parameter set, in that the isotope ratio of the fluid and the solid surface are both evolving and dependent on one another. DePaolo (2011) simplifies this relationship by assuming that the isotopic composition of the solid surface is at steady state, as in the case of the buffered titration experiments in which a fixed supersaturation is maintained (Tang et al. 2008; Watkins et al. 2013). In doing so, the formulation leads to a closed form solution for the observed fractionation factor (αobs) at steady state


where Rf and Rb are the forward and backward reaction rates (note Rf may be recast as Rnet + Rb to utilize the net observed reaction rate), αf is the fractionation factor associated with the forward addition of material to the precipitating mineral and αeq is the fractionation factor associated with dynamic equilibrium between the solid surface and surrounding fluid. In order to apply this model, the value of the backward rate of material removed from the solid surface, as well as the two fractionation factors αf and αeq must be supplied. DePaolo (2011) demonstrated that for an Rb of 2.16×103 μmol/m2/hr roughly corresponding to the values reported by Chou et al. (1989) and fractionation factors αf = 0.9984 and αeq = 0.9995 the model captures the trend observed by the Tang et al. (2008) study (Fig. 11B).

This model is based on the same concept of a dynamic system discussed in the introduction section. An implicit assumption of this affinity-based approach is that reactivity and fractionation are describable based on the composition of the fluid phase relative to equilibrium with the solid. In practice, the parameters required to apply this type of model often vary between systems, for example requiring a higher order dependence on the departure from equilibrium, or a distinct kinetic rate constant (Stack 2014). DePaolo (2011) points out that the value of αeq used to reproduce the Tang et al. (2008) dataset (0.9995) is inconsistent with independent estimates (Fantle and DePaolo 2007; Jacobson and Holmden 2008). An alternate parameter set is also able to fit the Tang et al. (2008) data using an αeq of 1.000 and a variable backward reaction rate (Rb), which DePaolo (2011) calculates as a function of the fluid saturation state.

Recent efforts have focused on replacement of the TST derivation with crystal growth models that include a description of the type and prevalence of reactive sites at the mineral surface in an effort to provide a more mechanistic and predictive capability (Wolthers et al. 2012; Bracco et al. 2013). This type of crystal growth model has been applied in combination with the model of DePaolo (2011) to derive solutions for calcium isotope partitioning during calcite precipitation as a function of solution stoichiometry (Nielsen et al. 2012) and for oxygen isotope partitioning during calcite precipitation as a function of pH (Watkins et al. 2013). The key development of the recent modeling approaches is to allow the influence of factors in addition to reaction affinity to influence the reaction rate and thus the isotopic compositions of the fluid and solid. However, these derivations still largely incorporate a steady state solution as demonstrated in Figure 11. Recently, Steefel et al. (2014a) discuss the need to explicitly describe the time-variant evolution of both the fluid and solid surface in order to model isotopic re-equilibration of a kinetically fractionated system. In the following section, we describe current modeling approaches to resolve the time dependence of isotopic exchange between fluids and solid surfaces in reactive systems during net precipitation, and show that aspects of this problem are uniquely resolvable at the pore scale.

Extending the concept of a variable apparent fractionation factor to through-flowing, field scale systems has thus far been accomplished primarily through the use of continuum-scale reactive transport models. As discussed in the dissolution section, these models have commonly been applied to simulate fractionation during homogeneous (aqueous) reactions. Accumulation of a solid phase that is composed of two isotopes of the same element is problematic in that any change in the isotope ratio of the fluid through time needs to be incorporated into the solid phase without imposing an erroneous fractionation factor. Druhan et al. (2013) accommodated a variable solid phase isotopic composition in the context of calcium isotope partitioning during carbonate mineral precipitation using a solid solution. In common practice, the overall rate of carbonate mineral precipitation:


where k is the kinetic rate constant, Keq is the equilibrium constant and a represents the activity of the individual species, is simplified by the fact that the activity of a pure solid phase is equal to unity (aCaCO3(s) = 1). In the Druhan et al. (2013) model, the isotopes of calcium are simulated as separate ‘species’, which are coupled through their precipitation rates as:



The key to this approach is the recognition that if aCaCO3(s) = 1, then the activities of the individual isotopologues of calcite must sum up to unity. Assuming an ideal solid solution, these activities can be rewritten as mole fractions of the total mineral:



where n is the number of moles of each isotopologue of calcite. Expressing the mole fractions of the solid phase in the calculation of 40Rnet and 44Rnet effectively couples the two isotopes to one another in an overall expression for carbonate mineral accumulation as a solid solution model. The problem then becomes the definition of the mole fractions X. At one extreme, the back reaction between the solid phase and fluid may be suppressed by setting X equal to the instantaneous mole fraction of the isotopes in solution. This results in a system where the solid phase is isotopically segregated from the fluid, and the only mechanism of fractionation is through an irreversible αk, that is the ratio of the kinetic rate constants. Such an approach would reproduce the isotopic dataset of the Skulan et al. (2002) hematite precipitation experiment (Fig. 9A). In systems where the back reaction is significant, the DePaolo (2011) conceptual framework would dictate that values of X represent the surface of the solid phase in contact with the surrounding fluid. In practice this value is not explicitly tracked in continuum scale models, which represent the composition of solid phases as volume fractions within a given discretized domain. Druhan et al. (2013) applied their model to a dataset of fluid calcium isotope ratios measured through time during unbuffered carbonate mineral precipitation. They noted that these data did not conform to a fixed fractionation factor through time (Fig. 12A), but an accurate simulation of the observed isotopic partitioning was obtained using Equations (31a,b) by approximating X as the isotopic composition of the bulk calcite as it accumulated through time (Fig. 12B).

Versions of this derivation in which isotopic fractionation was considered irreversible have been applied to chromium isotopes (Wanner et al. 2014) and sulfur isotopes (Druhan et al. 2014; Hubbard et al. 2014). Steefel et al. (2014a) described a variable fractionation factor for calcium isotopes in marine carbonates as a result of reversible reactivity using the same model, and provided a detailed derivation emphasizing the capacity to resolve isotopic re-equilibration unique to a model for the transient evolution of the coupled fluid and solid surface. Still, the accuracy of bulk-solid values of X presents a significant challenge because the composition of an isotopically distinct solid surface is homogenized with the preexisting solid over the volume of each gridblock. Druhan et al. (2013) and Steefel et al. (2014a) both provided approximations to the solid surface mole fractions by tracking a running average of the isotopic composition of material precipitating from solution as a function of time, resulting in an equally good representation of the measured fluid isotope ratios for a different value of αk. Similarly, one may track the isotope composition of a distinct solid surface layer and allow for exchange/equilibration with that zone of the solid based on the mineral volume fraction, surface area and a representative length scale of each zone. However, as discussed in the dissolution section, this approach requires some independent, a priori knowledge of the length scales of zoning as well as tracking of the isotopic composition of multiple zones within the total volume of a solid.

A common theme emerging from the previous sections is that partitioning of isotopes across a phase boundary requires some parameterization of the isotopic composition of the solid surface and surrounding fluid distinct from ‘bulk’ or ‘well mixed’ reservoirs. The ability to describe isotopic partitioning unique to the phase boundary appears to influence the accuracy of simulations from highly variable, field scale systems (e.g., Druhan et al. 2013) to highly controlled laboratory experiments (e.g., Tang et al. 2008). This observation implies that .mechanistic descriptions of isotope partitioning could be significantly improved by modeling approaches that are capable of resolving spatial zoning both within a single phase and within individual crystals or grains (e.g., Li et al. 2006; Tartakovsky et al. 2008; Molins et al. 2012; Molins 2015, this volume; Yoon et al. 2015, this volume).

One such method is the use of a lattice Boltzmann approach (Yoon et al. 2015, this volume), which is capable of explicitly describing the location and shape of a boundary between independently discretized fluid and solid phases. This approach has previously been applied to describe heterogeneous reactivity (Kang et al. 2002, 2003, 2004) including multicomponent systems (Kang et al. 2006, 2007, 2010; Yoon et al. 2012). Recently, Huber et al. (2014) developed an approach to describe the transfer of mass across a phase boundary that is independent of the shape of the solid surface. This approach is well suited to the description of fractionating reactions that are coupled through isotope partitioning across a phase boundary. As a proof of concept, Figure 13 demonstrates application of the Huber et al. (2014) model to a closed system in which the fluid phase surrounding a grain of calcite is initially supersaturated with respect to the solid. Two isotopes of calcium are explicitly tracked using the formulations given in Equation (31a,b) over time as the system reaches chemical and isotopic equilibrium. Figure 13 depicts a particular snapshot in time, during which the fluid phase has not yet reached equilibrium with the solid. This disequilibrium is evident in the gradients of reactive species as well as the fluid isotope ratio with distance from the solid phase boundary. Importantly, the surface of the solid phase, which is composed of pure CaCO3(s) is also isotopically variable, and supplies a more accurate estimate of the values of the mole fraction for each isotope X (Eqn. 32) than previous bulk averaged estimates. This is not readily observed in the spatial distribution of isotope ratios because only a very small portion of the solid phase surface is considered in contact with the surrounding fluid. However, a time series plot of the isotope composition of the calcite surface clearly demonstrates this variability, which feeds back into the rate at which the fluid phase established isotopic equilibrium with the reactive solid surface.

This proof of concept exercise suggests that significant improvements could be gained from pore scale descriptions of isotopic partitioning during precipitation, dissolution, transport and decay. For example, application of the Druhan et al. (2013) derivation in a system where X (Eqn. 32) is only describable as the bulk solid phase value would fail to reproduce observed isotope partitioning during dissolution. The same is true for the closed form solutions of DePaolo (2011) and Brantley et al. (2004). The result is a situation in which reactions that generate isotope partitioning, which are formulated as a reversible departure from equilibrium, cannot be described across the range of saturation states from net dissolution to net precipitation with a single model. The use of pore scale discretization may provide a means of circumventing this discontinuity.

This material is based in part upon work supported as part of the Subsurface Science Scientific Focus Area at Lawrence Berkeley National Laboratory funded by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research and Office of Basic Earth Sciences, under Award Number DE-AC02-05CH11231. The authors wish to thank Kate Maher and Carl Steefel for their thoughtful comments, which greatly improved the chapter.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.