Research ArticleGEOLOGY

Extreme Zr stable isotope fractionation during magmatic fractional crystallization

See allHide authors and affiliations

Science Advances  18 Dec 2019:
Vol. 5, no. 12, eaax8648
DOI: 10.1126/sciadv.aax8648


Zirconium is a commonly used elemental tracer of silicate differentiation, yet its stable isotope systematics remain poorly known. Accessory phases rich in Zr4+ such as zircon and baddeleyite may preserve a unique record of Zr isotope behavior in magmatic environments, acting both as potential drivers of isotopic fractionation and recorders of melt compositional evolution. To test this potential, we measured the stable Zr isotope composition of 70 single zircon and baddeleyite crystals from a well-characterized gabbroic igneous cumulate. We show that (i) closed-system magmatic crystallization can fractionate Zr stable isotopes at the >0.5% level, and (ii) zircon and baddeleyite are isotopically heavy relative to the melt from which they crystallize, thus driving chemically differentiated liquids toward isotopically light compositions. Because these effects are contrary to first-order expectations based on mineral-melt bonding environment differences, Zr stable isotope fractionation during zircon crystallization may not solely be a result of closed-system thermodynamic equilibrium.


Zirconium is a refractory lithophile and moderately incompatible transition metal whose concentration in igneous rocks has for decades been used as a tracer of differentiation in silicate magmatic systems (13). Although radiogenic and nucleosynthetic Zr isotopic anomalies in extraterrestrial samples have been extensively studied [e.g., (4, 5)], mass-dependent Zr isotope variations remain poorly explored and understood. Zirconium plays a key role in the development of accessory phases like zircon (tetragonal ZrSiO4) and baddeleyite (monoclinic ZrO2), which are fundamental to our understanding of geologic time and Earth’s crustal evolution [e.g., (6, 7)]. Therefore, unraveling the effects that high-temperature igneous processes have on the Zr stable isotope composition of bulk rocks and accessory phases can deepen our understanding of the Zr cycle in the silicate Earth, allowing us to better constrain the chemical evolution of the crust and mantle through geologic time.

Differences in the bonding environment of Zr in silicate liquids and Zr-rich minerals have the potential to result in mass-dependent isotope fractionation during fractional crystallization. Zr4+ ions in silicate melts occur mainly in octahedral (sixfold) coordination, with mean Zr–O bond distances between 2.05 and 2.13 Å (811). Although Zr solubility in silicate magmas is a strong function of temperature, pressure, and composition (12, 13), its speciation has been found to be comparatively insensitive to these variables (8, 11, 14). In contrast, Zr4+ is eight coordinated in the zircon structure (Zr–O = 2.14 to 2.28 Å), seven coordinated in baddeleyite (Zr–O = 2.04 to 2.18 Å), and can replace six-coordinated Ti in Ti oxide minerals like ilmenite (FeTiO3) and rutile (TiO2) (11, 15). Because heavy stable isotopes will typically populate shorter and stiffer bonds with the lowest coordination number at thermodynamic equilibrium (16), the expectation is that zircon and baddeleyite would be isotopically “light” relative to the melt from which they crystallize, driving the residual liquid toward complementarily “heavy” compositions as crystallization progresses.

To study the variability of Zr stable isotopes during magmatic fractional crystallization, we developed a high-precision and accuracy MC-ICP-MS (multicollector–inductively coupled plasma–mass spectrometry) method using the double-spike technique. This approach achieves analytical uncertainties of ±5 to 10 parts per million per amu. for Zr sample masses as low as 25 ng, thus enabling the study of Zr-rich phases such as zircon and baddeleyite at the single-crystal scale. A fresh rock fragment of the FC-1 anorthositic cumulate (17) from the Duluth Complex (DC) (Midcontinent Rift in North America) was selected for this study for several reasons: (i) zircons from this sample yield reproducible U-Pb dates at high precision (17, 18), indicating that these crystals are cogenetic and have behaved as closed systems since igneous crystallization; (ii) Ti-in-zircon thermometry suggests zircon crystallization took place over a >200°C temperature range (19), presumably during slow cooling; and (iii) the FC-1 locality also contains igneous baddeleyite with ages and initial 176Hf/177Hf ratios comparable to those of zircon (20), thus providing an opportunity to evaluate the distribution of Zr stable isotopes among cogenetic Zr-rich phases. Large plagioclase laths form the main orthocumulate framework of FC-1 (fig. S1), within which trapped intercumulus melt pockets differentiated internally to form abundant clinopyroxene and Fe-Ti oxides, and minor orthopyroxene, biotite, and apatite (fig. S1) (21). Zircon and baddeleyite occur exclusively within intercumulus pockets (Fig. 1, A and B, and fig. S2), indicating that Zr-rich minerals formed as late products of fractional crystallization of trapped interstitial liquids.

Fig. 1 Textural context and high-precision U-Pb geochronologic results of Zr-rich accessory phases from the Duluth Complex (FC-1 locality).

(A) Baddeleyite (bd) crystal in association with Fe-Ti oxide (ox) and clinopyroxene (cpx). (B) Zircon (zrc) crystal in association with ox and cpx. (C) 206Pb/238U dates. (D) 207Pb/206Pb dates. MSWD, mean square weighted deviation.

We determined the Zr stable isotope composition of 42 zircon and 28 baddeleyite single crystals, and two bulk-rock aliquots. Zircon and baddeleyite were digested individually in perfluoroalkoxy alkane (PFA) microcapsules at high pressure. Following complete digestion, all samples were spiked with a 91Zr-96Zr tracer at optimum levels before chemical purification and measurement by MC-ICP-MS. Nineteen zircons were treated using the “chemical abrasion” (CA) method [e.g., (18)], which involved thermal annealing at 900°C for 60 hours followed by a 12-hour partial dissolution step in 29 M HF at 215°C and high pressure, to selectively remove zircon domains with high accumulated radiation damage [(22); Supplementary Materials]. Seven of these zircon crystals were dated by U-Pb using ID-TIMS (isotope dilution–thermal ionization mass spectrometry). Seven baddeleyite crystals were also dated by U-Pb, but because no effective CA method yet exists for baddeleyite, all these crystals were dissolved “untreated” (i.e., not thermally annealed nor partially leached before complete digestion).


Zircon U-Pb analyses yield weighted mean 206Pb/238U and 207Pb/206Pb dates of 1095.97 ± 0.22 and 1099.96 ± 0.58 Ma, respectively (Fig. 1, C and D, and figs. S3 and S4). Baddeleyites, on the other hand, yield younger 206Pb/238U dates ranging from 1092.74 ± 0.53 to 1094.70 ± 0.57 Ma, likely reflecting variable Pb loss. A weighted mean 207Pb/206Pb date of 1101.41 ± 0.50 Ma for these baddeleyites is resolvably older than the zircon 207Pb/206Pb weighted mean, implying that either (i) baddeleyite crystallized before zircon or (ii) baddeleyite incorporated excess 231Pa during crystallization, making its apparent 207Pb/206Pb age older (20).

The Zr isotopic composition of all zircon and baddeleyite crystals, expressed as δ94/90ZrNIST [i.e., deviation in parts per thousand relative to a new Zr isotopic reference material being developed in collaboration with the National Institute of Standards and Technology (NIST); see Materials and Methods], covers a range of ca. 5.3‰ from +0.971 ± 0.010 to −4.278 ± 0.010‰ (Fig. 2 and Table 1). The mean bulk-rock δ94/90Zr value is −0.086 ± 0.016‰. This observed range of mass-dependent variation from a single hand specimen is one order of magnitude larger than the range recently reported for Zr in bulk-rock samples (23) and possibly the largest stable isotope variability yet observed for a transition metal in igneous systems [e.g., maximum δ49/47Ti ~1‰, δ98/95Mo ~0.6 ‰, and δ53/52Cr ~0.3 ‰ (24)]. Because zircon crystals treated using CA before Zr isotope determinations cover a similarly large range of δ94/90Zr values relative to our entire dataset (Fig. 2A), the spread we observe cannot be an artifact of radiation damage. The very low Mo/Zr ratios determined in our samples (Table 1) imply negligible interference-related δ94/90Zr offsets to our reported data (figs. S6 and S7), therefore demonstrating that the extreme δ94/90Zr variations reported here are not an artifact of unaccounted isobaric interferences.

Fig. 2 Zr stable isotope results of bulk-rock and Zr-rich phases of the FC-1 anorthositic gabbro, expressed as δ94/90ZrNIST (i.e., [[94/90ZrSample/94/90ZrNIST] − 1]*1000).

(A) Frequency histogram (bin width = 0.1‰) and optimum-bandwidth Gaussian kernel density estimates (kde) of zircon (blue; n = 42) and baddeleyite (green; n = 28) δ94/90Zr values. Gray background is the expected probability density function of δ94/90Zr values expected from an αsol-liq ≈ 1.00106, obtained using 105 Monte Carlo simulations (Supplementary Material). Arrows indicate the approximate composition of chemically abraded zircons that were U-Pb dated using ID-TIMS (dashed) and those treated using CA but not U-Pb dated (solid). (B) Interpretation of the zircon and baddeleyite δ94/90Zr data as a Rayleigh distillation process. A minimum fractionation factor (Δ94/90Zrsol-liq ≈ 1.06 ‰) was derived by assuming that the “heaviest” measured solid was the first to crystallize from a liquid with an initial composition equal to that of the bulk rock. Gray shading represents the region of the diagram that bulk-crystal compositions can occupy if individual grains represent integrated growth nuclei rather than “instantaneous” solids. See text for further discussion.

Table 1 Results of Zr stable isotope analyses of FC-1 zircon and baddeleyite crystals by double-spike MC-ICP-MS, and zircon Zr/Hf analyses by solution quadrupole–inductively coupled plasma–mass spectrometry.

At %, atomic percent; CA, chemically abraded before complete digestion; TEA, element weight ratios from trace element analyses; Untr., untreated using CA procedures before digestion.

View this table:


Origin of Zr stable isotope variations in FC-1 zircon and baddeleyite

The observed distribution of zircon and baddeleyite δ94/90Zr values is inconsistent with a closed-system equilibrium scenario because these solids are both isotopically “heavy” and “light” relative to the bulk-rock value (Fig. 2A). If Zr-rich phases formed as products of fractional crystallization become effectively isolated from isotopic exchange with the rest of the system upon their formation (e.g., due to gravitational removal or a kinetic impediment such as slow diffusion), then a Rayleigh fractionation mechanism may best describe the observed differences in δ94/90Zr (25). Isotopic isolation of Zr during zircon and baddeleyite crystallization is supported by simple arguments made from diffusion kinetics. Although no data for Zr diffusion in zircon or baddeleyite are currently available, the relationship between preexponential factor for diffusion (D0) and ionic radius of tetravalent species in zircon (26) can be used to approximate Zr4+ diffusivity to logD0 ≈ 2.67 m2 s−1. Assuming a 190 kcal/mol activation energy for diffusion, which is typical of tetravalent cations in zircon (26), the effective diffusion distance (2√Dt) for Zr in zircon at 900°C is predicted to be on the order of 10−4 μm/million years (Myr). Thus, considering the time scale of zircon crystallization in FC-1 as constrained by our U-Pb data [i.e., maximum of 220 thousand years (kyr); Fig. 1], the extremely slow diffusivity inferred for Zr implies that, once incorporated into a solid, Zr is essentially immobile within the temperature range relevant to FC-1 zircon. This kinetic impediment arising from slow diffusion should lead to negligible isotopic reequilibration between the liquid and crystallizing solids, therefore resulting in effective Zr removal from the melt and Rayleigh-type fractionation.

An interpretation of the δ94/90Zr data as a Rayleigh fractionation process is shown in Fig. 2B, according to[Z94rZ90r]=[Z94rZ90r]if(αsol-liq1)where [94Zr/90Zr]i represents the initial bulk composition of the melt, f is the Zr mass fraction removed from the liquid, and αsol-liq is an instantaneous mass-dependent isotopic fractionation coefficient (25). Although αsol-liq between zircon and silicate liquids is unknown, a minimum value can be derived from our data as the apparent fractionation coefficient needed to explain either the heaviest or the lightest solid measured relative to the bulk-rock value. Because Rayleigh-type fractionations between liquid and crystallizing solids impart extreme isotopic compositions to residual liquids as f decreases (25), using the lightest observed solid to estimate αsol-liq would imply that a significant number of zircon and/or baddeleyite with δ94/90Zr values in excess of +5‰ should be present in the sample (fig. S6), which were not observed (Fig. 2A). Likewise, adopting the more subdued but still negative fractionation factor (Δ94/90Zr) of −0.5‰ proposed by (23), based on indirect inferences from bulk-rock compositions, does not fit the data because (i) the majority of the measured zircon should be lighter than the bulk-rock value, which is contrary to our results (Fig. 2A); or (ii) solids with δ94/90Zr values heavier than 1‰ would be very likely to be present (fig. S6), but these were not observed; and (iii) this fractionation coefficient cannot explain zircon and/or baddeleyite δ94/90Zr values lighter than −0.5‰, which are present in FC-1.

On the other hand, a positive fractionation factor between mineral and melt during zircon and/or baddeleyite crystallization can readily explain the data presented here. The expected distribution of δ94/90Zr values, calculated using the heaviest observed solid to estimate a minimum value for αsol-liq, provides an excellent match to the observed distribution (Fig. 2A and fig. S6). On the basis of these observations, we argue that fractional crystallization of isotopically heavy Zr-rich phases provides the best explanation to the data and a lower limit on αsol-liq. The isotopically heaviest baddeleyite and zircon have similar δ94/90Zr values of +0.971 ± 0.010 and + 0.905 ± 0.010‰, respectively, so adoption of one value over the other to infer αsol-liq has a negligible effect on our conclusions. Using this approach, and given the whole rock δ94/90Zr value of −0.086 ± 0.016 ‰, an αsol-liq ≈ 1.00106 (or Δ94/90Zr ≈ 1.06) was used to calculate the analytical Rayleigh curves for liquid, instantaneous solid, and bulk solid shown in Fig. 2B.

Zr isotope constraints on the extent of Zr fractional removal from a crystallizing magma

Because zircon and baddeleyite crystals growing in a magma are not infinitesimally small solids but rather nucleation seeds that incorporate a fraction (if not most) of the Zr available to them in the surrounding melt (27), each crystal is unlikely to represent an instantaneous f value. Rather, each zircon or baddeleyite crystal represents a “nucleation seed,” which may integrate a significant fraction of Zr removal history from its surrounding liquid, and thus, the δ94/90Zr value of any bulk grain analysis must represent a weighted mean of its internal growth zoning. Variable growth integration histories imply that bulk crystal analyses could populate the entire compositional range between the bulk-solid and instantaneous solid curves on a Rayleigh fractionation diagram (gray region in Fig. 2B). For the purposes of illustration, all individual zircon and baddeleyite analyses from FC-1 shown in Fig. 2B are placed along a “50% solid integration” curve, which is the locus defined by nucleation seeds that integrate 50% of the Zr available to them in the surrounding melt during crystallization.

To explore the power and limitations that bulk-crystal δ94/90Zr values present for quantifying the extent of Zr removal from a magma as recorded by single zircon or baddeleyite grains, we derived analytical equations (see Materials and Methods) to quantify the effect of variable growth integration histories on the isotopic composition of bulk individual grains. Figure 3A shows the results of these calculations, where the solid line through the center represents the f value as a function of δ94/90Zr assuming 0% growth integration (instantaneous solid), while the blue and red dashed contours represent minimum and maximum permissible f values for the beginning and end of crystal growth, respectively, as a function of incremental growth integration. Crystals with δ94/90Zr values approaching the bulk-rock composition have greater uncertainties in the assignment of f values for both the initiation and termination of crystal growth (Fig. 3B), but the “tails” of the distribution yield much more powerful constraints. This is because there are less possible growth integration “paths” that could result in these bulk crystal isotopic compositions, as suggested by our isotopic results (Fig. 2A) and inferred from theoretical probability density function calculations (Fig. 2Aand fig. S4). Nevertheless, we note that even in the near-limiting case of zircon and baddeleyite bulk compositions representing the weighed result of a 0.95 fractional integration history (i.e., 95% integration uncertainty), the relative sequence of nucleation between the isotopically heaviest and lightest crystals in our FC-1 dataset can be confidently resolved (table S2).

Fig. 3 Uncertainties in the assignment of f for the initiation (f initial; blue) and end (f final; red) of zircon growth from FC-1 bulk crystal δ94/90Zr data.

(A) Permissible values of f final and f initial as a function of fractional growth integration. Note that the black solid line represents the instantaneous solid (0% integration) of a Rayleigh fractionation model (i.e., when f initial ≡ f final). Numerical labels for each dashed curve represent the fraction of Zr incorporated into each growth-zoned “seed,” relative to the amount of Zr available to that crystal in the surrounding melt at the time of crystallization. (B) Absolute uncertainty in the assignment of f final and f initial values to single zircon/baddeleyite crystals from FC-1, as a function of their bulk δ94/90Zr value.

Our analysis indicates that, despite uncertainties associated to unknown growth integration histories, the δ94/90Zr value of a zircon (or baddeleyite) can be a sensitive monitor of Zr removal from the melt, therefore allowing a relative crystallization sequence within a sample to be developed even at a resolution beyond the current resolving power of high-precision U-Pb dating (Fig. 1). Another key implication of our observations is that zircon crystals from FC-1 should be internally zoned in terms of δ94/90Zr, and therefore, in situ analyses of Zr isotopes (e.g., by secondary ionization mass spectrometry) could be a valuable complement to bulk crystal analysis to fully disentangle the complex δ94/90Zr history that may be preserved within individual minerals.

Relationships between δ94/90Zr, magma compositional evolution, and crystallization of Zr-rich phases

Although the large resolvable differences in δ94/90Zr found among zircon crystals are best explained by sequential precipitation from a fractionating magma, the fact that zircon U-Pb dates overlap within uncertainty (Fig. 1) means that the apparent sequence of zircon crystallization in FC-1 suggested by the δ94/90Zr data cannot be verified using current high-precision geochronology methods. On the other hand, the trace element composition of individual zircon crystals can provide insights into melt compositional evolution, and thus the relative extent of differentiation a magma had undergone at the time of zircon crystallization. Zr and Hf exhibit almost identical geochemical properties, and although Zr/Hf variations in igneous systems are known to exist, large fractionations are not dominantly driven by common rock-forming phases (28). During zircon crystallization, however, Zr is more effectively incorporated into the solid than Hf, thereby rapidly decreasing the Zr/Hf ratio of derivative melts (and of late formed zircon) due to the progressively increasing Hf concentration in the melt (28). Therefore, the Zr/Hf ratio of a zircon can provide a unique tracer not only of melt evolution but specifically of zircon removal from a fractionating magma. Figure 4 shows the Zr/Hf elemental concentration ratio of FC-1 zircon crystals studied here as a function of their Hf concentration and δ94/90Zr value, and compares them with the Zr/Hf ratio of crystals from well-documented localities where the relative sequence of zircon crystallization has been independently verified. In the Bergell intrusion (Central Alps), zircon Zr/Hf ratios correlate positively with their U-Pb crystallization age (29), and in the Spirit Mountain batholith (Nevada, USA), zircon Zr/Hf also exhibits a positive correlation with crystallization temperatures from Ti-in-zircon thermometry (28). The Zr/Hf ratio of FC-1 zircon decreases with decreasing δ94/90Zr values (Fig. 4), therefore supporting the hypothesis that isotopically heavy crystals formed earlier, and isotopically light crystals formed late from a highly fractionated melt.

Fig. 4 Zr/Hf concentration ratios in zircon against indices of magmatic differentiation.

Zircon crystals from the Duluth gabbro studied here overlap in Zr/Hf versus Hf concentration space with zircon from the Bergell intrusion (29) and Spirit Mountain batholith (28), which exhibit a decrease in Zr/Hf as a function of decreasing age and crystallization temperature, respectively. Zr/Hf ratios in FC-1 zircon decrease with decreasing δ94/90Zr values, supporting the interpretation that isotopically heavy zircon crystallized early, whereas isotopically light zircon crystallized from a highly fractionated melt.

To better place our δ94/90Zr results within the crystallization history of FC-1, we combined numerical simulations using rhyolite-MELTS (30) with a zircon solubility parametrization (13). Equilibrium fractional crystallization calculations were performed isobarically at 200 MPa along the quartz–fayalite–magnetite (QFM) buffer (31), using the bulk chemical analyses of our FC-1 sample as the starting composition. Figure 5 shows the forward-modeled phase assemblage as a function of decreasing temperature during crystallization, which accurately captures observed modal abundances as well as textural relationships among phases in the FC-1 sample (fig. S1). The melt escape threshold (~20% liquid) and liquid percolation threshold (~8% liquid) represent critical lower limits to melt permeability and melt pocket connectivity (32), respectively. Below ~8% liquid (>92% crystallization), melt pockets are physically isolated, preventing mass transfer between them and, thus, effectively imposing closed-system behavior for the remaining cooling history (Fig. 5A).

Fig. 5 Numerical simulations of fractional crystallization of FC-1 and zircon saturation/growth temperatures.

(A) Equilibrium fractional crystallization model at 200 MPa and fO2 = QFM using rhyolite-MELTS. Note change in vertical scale at mass fraction = 0.6. Melt escape threshold (MET) and liquid percolation threshold (LPT) limits from (32). Column on the right shows the observed mineral modes in FC-1 (Supplementary Materials). Bt, biotite; A, apatite. (B) Compositional parameters of the residual liquid and conditions for zircon saturation. The [Zr]liq curve uses the bulk-rock Zr concentration of 39 μg/g as starting point and assumes incompatibility in modal phases. The [Zr]sat-liq curve is the solubility limit of Zr in the liquid, calculated using the melt composition at each step and the solubility parametrization of (13); dashed range indicates values above the calibration range (i.e., M > 2.1). The gray shaded region represents the range of calculated temperatures for the onset of zircon saturation, which change as a function of initial H2O contents (see text for details). Red circles represent FC-1 zircon crystallization temperatures from Ti-in-zircon thermometry (19) recalculated using aTiO2 = 0.35 ± 0.05. Dotted red curve is a cumulative PDF of the Ti-in-zircon data, and dotted blue curve is an integrated PDF of zircon crystallization temperatures from our MELTS modeling within the explored range of initial H2O contents. Regardless of the assumed initial H2O, the modeled onset of zircon saturation temperature is always below the LPT, indicating intracumulus pockets behaved as closed systems during zircon crystallization.

Using the liquid composition retrieved from MELTS at each calculation step, we applied the parametrization of (13) to calculate the concentration of Zr necessary to saturate zircon in the liquid (Fig. 5B). Once the calculated Zr concentration in the residual melt (solid blue curve in Fig. 5B) reaches saturation (black curve), the amount of removed zircon is estimated at each step using mass balance and assuming that the melt is kept at Zr saturation during the remaining cooling path. Subtle variations in H2O content can noticeably influence phase stability and zircon crystallization temperatures (33). Therefore, although the initial H2O content of the FC-1 parental melt was likely below 0.5 weight % (wt %) (i.e., typical for anhydrous, tholeiitic magmas), to accommodate this uncertainty, we performed multiple calculations from 0.05 to 0.50 wt % initial H2O in 0.025–wt % increments. The predicted temperatures for the onset of zircon saturation change rapidly from ca. 975°C at 0.05 wt % H2O to 870°C at 0.25 wt % H2O, but then decrease negligibly up to 0.5 wt % H2O (Fig. 5B). A probability density function for zircon crystallization in FC-1, calculated by stacking the results of all our MELTS calculations (Fig. 5B; dotted blue curve), demonstrates that regardless of the assumptions made in our model, zircon is always predicted to saturate at melt fractions below the LPT. This implies that zircon crystallization took place in a closed system during fractional crystallization of late-stage intercumulus melt pockets, and thus, their δ94/90Zr values must reflect closed-system behavior. The temperature range of zircon crystallization obtained with our model is in good agreement with the range obtained from Ti-in-zircon thermometry in FC-1 zircon (Fig. 5B; dotted red curve), after recalculation of aTiO2 following the approach of (34).

Although a parametrization for baddeleyite saturation does not yet exist, the distribution of δ94/90Zr between zircon and baddeleyite crystals yields insight into their relative temporal relationships and the sources of the U-Pb date discrepancies between these two phases (Fig. 1). If the magma saturated in baddeleyite earlier than in zircon in its differentiation history, as suggested by the older 207Pb/206Pb mean baddeleyite date (Fig. 1D), then a systematically different δ94/90Zr should be recorded by baddeleyite relative to zircon. The δ94/90Zr range of both minerals, however, are virtually identical (see kernel density estimate –kde– curves in Fig. 2A), strongly suggesting cocrystallization. This observation aids the interpretation of the high-precision U-Pb data, as it implies that (i) the systematically younger 206Pb/238U dates of baddeleyite relative to zircon must be the result of variable Pb (or intermediate U-series nuclides) loss and (ii) the older 207Pb/206Pb weighted mean date determined for baddeleyite is most likely due to initial excess of 231Pa in baddeleyite (20) rather than truly reflecting an older crystallization age.

Implications for petrogenetic interpretations of bulk-rock δ94/90Zr values

Our Zr isotope results clearly demonstrate that FC-1 zircon and baddeleyite are isotopically heavy relative to the liquid from which they crystallized, which has important implications for the interpretations of bulk-rock δ94/90Zr signatures. Recently, an apparent positive correlation between δ94/90Zr and SiO2 wt % in lavas from the Hekla volcano, Iceland, was posited as reflecting fractional removal of isotopically light zircon during magmatic differentiation, resulting in the elevated δ94/90Zr values of evolved eruptive products (23). This interpretation is not only in stark contrast with the detailed investigation of zircon δ94/90Zr values presented here (Fig. 2) but also problematic because no zircon crystals from Hekla have been studied, making the connection between fractionations observed at the bulk-rock scale and their potential mineral driver(s) highly speculative. Furthermore, although Hekla’s more differentiated eruptive products have been shown to contain zircon, their U-Th crystallization ages are significantly older than, and trace element and isotopic compositions arguably in disequilibrium with, the lavas in which they are hosted (35, 36). Therefore, instead of representing zircon crystallized from the host eruptive products in which they are found, the vast majority of the zircons found in Hekla eruptions were inherited from older magma batches (i.e., xenocrystic). If these inherited xenocrystic zircons are isotopically heavy, assimilation could at least partially be at fault for the elevated bulk-rock δ94/90Zr of evolved Hekla lavas. Although a detailed investigation of Hekla zircons will be necessary to resolve these issues, the clearly heavy Zr isotope composition of FC-1 zircon and baddeleyite relative to the bulk-rock value indicates that the role of older crustal melting and assimilation in the evolution of Hekla’s evolved lavas [e.g., (37)] and the zircons they contain cannot be overlooked.

While our results clearly show that FC-1 zircon and baddeleyite are isotopically heavy relative to the liquid from which they crystallized, this observation is opposite to simple predictions made from shifts in Zr4+ bonding environment during crystallization. The reason for this discrepancy is currently unknown, and additional work will be necessary to explain this apparent disagreement. In particular, first principle ab initio calculations would help confirm the simple expectations based on zircon and baddeleyite crystal chemistry, and zircon growth experiments covering the compositional range of igneous systems will help in exploring the conditions and drivers of Zr stable isotope fractionation. Such fundamental studies, which are beyond the scope of the present work, will better allow bulk-rock and accessory phase δ94/90Zr values to be linked to petrogenetic processes, and the magnitude of the fractionation coefficients involved to be better quantified.

In summary, the results presented here provide unprecedented insights into the behavior of Zr stable isotopes during closed-system magmatic differentiation and allow us to begin developing a petrogenetic understanding of the significance of fractionated δ94/90Zr values in igneous environments. We found that (i) zircon and baddeleyite crystallization plays a major role in driving the isotopic evolution of Zr in silicate magmas, and (ii) these phases are isotopically heavy relative to the melt from which they crystallize, thus driving evolved liquids toward isotopically lighter δ94/90Zr values. Therefore, crystallization of Zr-rich accessory phases may not only provide a mechanism for inducing isotopic fractionation in silicate magmatic systems but also offer a robust record of magma compositional change. Investigation of the Zr stable isotope systematics in igneous rocks covering a large compositional range will allow δ94/90Zr values to become a powerful tracer of magmatic differentiation as recorded in single accessory and possibly even in major phases. As “non-traditional” stable isotope investigations continue to yield novel insights into magmatic and other high-temperature processes (3841), the ability to determine Zr stable isotope compositions at high precision from individual minerals, particularly zircon crystals, can provide a valuable new tool to probe the record of magmatic differentiation and crustal evolution from the igneous, metamorphic, and detrital zircon archive.


Duluth anorthositic gabbro sample (FC-1)

The DC encompasses a diverse association of plutonic and hypabyssal intrusions related to the 1.1-Ga Midcontinent Rift in North America, a failed intracontinental rift that was active for ca. 25 Myr (17). The best exposures are located in northern Minnesota, where gabbroic, anorthositic, and granitic intrusions of the DC underlie the North Shore Volcanic Group, which is thought to represent a comagmatic volcanic edifice to the DC (21). Within the DC, a section known as the “anorthositic series” comprises approximately 35 to 40% of the exposed complex. These anorthosites have 80 to 98% cumulus plagioclase, minor cumulus olivine (Fo66-38), and intercumulus clinopyroxene (En75-57), Fe-Ti oxides, and various late-stage and secondary minerals. For this study, fresh anorthositic gabbro samples were collected from the extensively studied FC-1 locality (17), with permission from the U.S. Forest Service. This outcrop is known for yielding abundant zircon and baddeleyite (17, 20), and zircons yield reproducible U-Pb dates. For this reason, FC-1 is widely used as a geostandard locality for zircon U-Pb dating. Hafnium (176Hf/177Hf) isotopic compositions of zircon and baddeleyite from FC-1 are also reproducible within measurement uncertainty at high precision (20, 42), which supports their cogenetic nature and makes them a valuable geostandard for Hf isotopic measurements.

Thin section observations from our studied FC-1 sample are in line with previous descriptions of this locality. Modal abundances of rock-forming mineral phases in FC-1 estimated from three separate thin sections are approximately 77% plagioclase, 13% clinopyroxene, 6% Fe-Ti oxides, 2% biotite, 2% relict olivine, <1% orthopyroxene, and accessory apatite (fig. S1). Plagioclase crystals form the main orthocumulate framework, within which most other phases crystallized. The main intercumulus minerals are clinopyroxene and Fe-Ti oxides, with minor late-stage biotite that dominantly grew around, and at the expense of, Fe-Ti oxides. A relict olivine crystal was observed in one of our thin sections [presumably a cumulus phase; e.g., (21)], but this crystal is surrounded by a symplectitic corona, which indicates disequilibrium with the surrounding melt and likely captures an incomplete peritectic reaction well known to occur in anorthosites during crystallization. Thus, olivine is omitted from the estimated modes of rock-forming phases in Fig. 3A. These petrographic observations are in good agreement with the modal abundances and crystallization sequence of phases obtained from our rhyolite-MELTS equilibrium crystallization simulations described in the text and below (see the “MELTS and zircon crystallization modeling” section in Supplementary Materials).

The thin section of FC-1 shown in fig. S1 was mapped for Zr concentration using an electron microprobe to identify the textural relationships between Zr-rich phases and major rock-forming minerals. Once zircon and baddeleyite crystals were identified, high-resolution secondary electron images and x-ray energy dispersion spectra maps of selected crystals (fig. S2) were collected using a Hitachi S-4700 II Scanning Electron Microscope at the Electron Microscopy & Microanalysis Facility, University of Nevada, Reno. Zircon and baddeleyite crystals were texturally found to occur enclosed within or around the periphery of intercumulus pockets and minerals (fig. S2), which indicate that these are late-stage phases that formed during differentiation of trapped intercumulus liquids. These textural observations are not only in excellent agreement with zircon crystallization simulations described in the text and below (see the “MELTS and zircon crystallization modeling” section in Supplementary Materials) but also in line with our broader understanding of how and when zircon crystallization takes place in mafic magmas (13).

Zirconium isotopic reference material

Because of the dearth of previous Zr stable isotopic investigations, no standards of certified isotopic composition are available. The NIST produces the Standard Reference Material SRM-3122, which is a gravimetric solution of certified Zr elemental concentration, but whose isotopic composition has not been calibrated nor certified. To remediate this, the authors initiated a collaboration with NIST scientists to develop a new Zr isotopic reference material (henceforth referred to as ZrNIST). A large batch of a concentrated gravimetric Zr solution was prepared at NIST using an ultrahigh-purity Zr metal rod as starting material, and an aliquot of that solution was used for this investigation. An exhaustive interlaboratory calibration of the original NIST solution is currently underway and will be reported in a future contribution. It should be noted, however, that regardless of the absolute isotopic composition that will be recommended from this interlaboratory calibration effort, the magnitude of mass-dependent deviations reported here (using δ notation) will remain unmodified, as they represent offsets relative to the standard composition. Therefore, the absolute isotopic composition of the new ZrNIST is immaterial for the conclusions of this study, and our results will be comparable to all future Zr isotopic investigations using the ZrNIST standard once it becomes commercially available.

Mineral separation

Minerals analyzed in this study were concentrated using standard density and magnetic separation methods. Rocks were crushed using a stainless steel mortar and sieved through a disposable plastic mesh with 375-μm openings. Approximately 2 kg of powder was then washed in a plastic gold pan to remove clay-size particles and preconcentrate the heavy minerals.

After drying, the remaining material was further separated using magnetic methods: Highly magnetic particles (e.g., magnetite) were removed using a hand magnet, and the remaining sample was passed through a Frantz LB-1 separator at progressively increasing magnetic fields. The fraction not susceptible to the field generated by 1.1 A in the Frantz was processed through methylene iodide (p = 3.32 g/cm3) purchased from GeoLiquids ( The “floats” consisted of mostly plagioclase and apatite, while the “sinks” were almost pure zircon and baddeleyite in a proportion ~20:1. Single zircon and baddeleyite crystals selected for analysis were individually handpicked under high-purity ethanol using a binocular microscope.

Chemical abrasion of zircon crystals

Roughly half (n = 19/42) of the zircon crystals selected for analysis were treated using the CA technique before full dissolution. This approach was originally developed by (43) to selectively remove zircon domains with high degrees of accumulated α-recoil radiation damage and that had therefore not behaved as isotopically closed systems. The original CA approach was in later years modified and improved for application to single-zircon crystals [e.g., (18)], and it is now the standard practice in high-precision zircon U-Pb geochronology. The key advantage of the CA routine is that by selectively removing zircon domains that may have behaved as “open systems” due to the enhanced element mobility enabled by radiation damage, solid residues representing pristine zircon domains can subsequently be dissolved and analyzed.

Zircon crystals selected for CA pretreatment were annealed at 900°C for 60 hours and then leached in 29 M HF in a Parr vessel held at 215°C for 12 hours at the University of Rochester Laboratory for Isotope Geochemistry (UR-LIG). The “chemically abraded” fragments were rinsed twice in Milli-Q H2O, fluxed on a hot plate (~90°C) in ultrapure 6.2 N HCl, rerinsed thrice in Milli-Q H2O, and loaded individually into clean Teflon microcapsules.

U-Pb geochronology

Following treatment using the CA approach, eight zircon crystals selected for CA-ID-TIMS (chemical abrasion–isotope dilution–thermal ionization mass spectrometry) were loaded individually in clean Teflon microcapsules and spiked with the EARTHTIME ET-535 (205Pb-233U-235U) isotopic tracer (44) at the Massachusetts Institute of Technology Isotope Laboratory (MIT-IL). Complete dissolution was achieved using 29 M HF and placing the microcapsules inside a Parr digestion vessel that was held at 215°C for 48 hours. After conversion to a chloride matrix, U and Pb were separated by anion exchange chromatography using AG-1X resin in microcolumns with 50 μl of stem volume. Pb and U were loaded on a single outgassed Re filament in 2 μl of a silica gel/phosphoric acid mixture, and isotopic measurements were made on an Isotopx X-62 TIMS instrument at MIT. Pb isotopes were measured by peak hopping using a Daly photomultiplier and corrected for instrumental mass fractionation using a linear law and α = 0.18 ± 0.02%/amu, which is based on repeat analyses of the NBS-981 Pb isotopic standard. Uranium was measured statically as UO2+ ions using Faraday cups with 1011 ohm amplifying resistors and corrected for instrumental mass fractionation using the known ratio of 233U/235U in the ET-535 tracer (44), assuming a sample 238U/235U ratio of 137.818 (45) and a 18O/16O of 0.00205 (44).

Baddeleyite crystals dated by U-Pb were subjected to a similar procedure with a few differences. Because no CA treatment is yet available to efficiently remove baddeleyite domains that have experienced high degrees of radiation damage, all crystals dated were dissolved untreated (i.e., not thermally annealed nor partially leached before complete digestion). Before loading in microcapsules, individual crystals were fluxed for 2 hours on a hot plate at 60°C in Teflon beakers containing 1 ml of ultrapure 20% HNO3. After this, the HNO3 solution was removed, and crystals were rinsed thrice with Milli-Q H2O. Crystals were then loaded individually into clean Teflon microcapsules and spiked with ET-535 before complete digestion. The rest of the procedure followed the same steps described above for the zircon analyses.

During U and Pb purification of zircon and baddeleyite fractions using AG-1X resin, 3 M HCl washes containing the sample’s Zr, Hf, and other trace elements were collected in precleaned 7-ml Teflon beakers using an approach similar to that of (42). These aliquots were subsequently processed for Zr stable isotope measurement using the methods described below (next section).

The total mass of non-radiogenic Pb measured in our FC-1 zircon and baddeleyite fractions is indistinguishable from the range of Pb determined in total procedural blanks, therefore indicating that the source of non-radiogenic Pb in our samples is laboratory blank rather than initial non-radiogenic Pb in the crystal structure. Therefore, all Pb isotopic ratios were corrected for laboratory blank based on the measured 204Pb abundance and using the MIT-IL laboratory blank composition (i.e., 206Pb/204Pb = 18.128949 ± 0.477810, 207Pb/204Pb = 15.283551 ± 0.298985, and 208Pb/204Pb = 37.043490 ± 0.884747; 1σ, absolute).

U-Pb dates and uncertainties were calculated from the ID-TIMS data following the algorithms of (46) and using the U decay constants of (47). 206Pb/238U ratios and dates were corrected for initial 230Th disequilibrium using a Th/U[magma] of 2.8 ± 1.0. U-Pb results are presented in table S1 and shown plotted in Wetherill concordia space in fig. S3. Weighted mean dates are shown separately in fig. S4, where quoted uncertainties are presented in the form ±X/Y/Z (18) using the following scheme: X is analytical uncertainty only, Y includes propagation of the analytical and ET-535 tracer calibration uncertainties, and Z is the combined analytical, tracer, and U decay constant uncertainty. Therefore, the first level of uncertainty X is appropriate when comparing the zircon and baddeleyite data obtained in this study, and thus, only those uncertainties are shown in Fig. 1 and discussed in the main text.

Zr stable isotope analyses

Zircon crystals not pretreated by CA and baddeleyites not dated by U-Pb were individually fluxed for 2 hours on a hot plate at 60°C in PFA beakers containing 1 ml of ultrapure 20% HNO3. After this, the HNO3 solution was removed, and crystals were rinsed twice with Milli-Q H2O before being individually loaded into clean PFA microcapsules. Full dissolution was achieved using 29 M HF and by placing the microcapsules in a large-volume Parr digestion vessel that was held at 215°C for 60 hours. Following the dissolution step, microcapsules were visually inspected using a binocular microscope to ensure that zircon and baddeleyite crystals had fully dissolved; complete digestion was achieved in all instances after the 60-hour heating period.

Solutions containing the digested samples in ca. 100 μl of 29 M HF were transferred into precleaned PFA beakers containing ~1 ml of ultraclean 3 M HNO3 + 0.5 M HF, and small aliquots ca. 50 μl (5% of sample) were taken for Zr and Hf concentration measurements using an Agilent 7700 quadrupole ICP-MS (Q-ICP-MS) (see methods below). Once the Zr concentration in each sample was known, solutions were spiked with our in-house 91Zr-96Zr isotopic tracer in optimal proportions (i.e., 0.43:0.57 spike-to-sample Zr mass ratio). Spiked aliquots were fluxed at 140°C on a hotplate overnight and dried down twice to a salt and redigested to ensure complete spike-sample equilibration, which is vital for obtaining accurate isotopic results. Zr was purified using a procedure modified from (48), which was optimized at the UR-LIG to improve Zr-Hf separation (fig. S5). In brief, we used ~300 μl of Ln-Spec resin (25- to 50-μm mesh) in 10 cm by 2 mm columns, fitted in an Eichrom vacuum box. New (precleaned) resin was loaded in the columns for each separation and cleaned using ~7 resin volumes (2 ml) of 2 M HF, ~3 resin volumes (1 ml) of MQ-H2O, and conditioned using ~7 resin volumes (2 ml) of 6 M HCl + 0.06 M HF before sample loading. Samples were loaded in two steps of 0.5 ml of 6 M HCl + 0.06 M HF, and rare earth elements (REEs) and Th were rinsed using an additional 2.5 ml of 6 M HCl + 0.06 M HF. Zr was eluted in 8.5 ml of the same acid mixture, resulting in a cut that contains >95% of Zr, undetectable REEs, and <3% of initial Hf (fig. S5). After the separation was completed, the Zr cuts were dried down, fluxed at 120°C for 2 to 3 hours in 2 ml of a concentrated HNO3-H2O2 mix to burn off any organic residues from the resin (49), dried down again, and redigested in 2 ml of 0.59 M HNO3 + 0.28 HF. A 50-μl aliquot from this solution was taken for concentration measurement on the MC-ICP-MS. Samples were then diluted to achieve a final concentration of 60 ng/g (total Zr) before MC-ICP-MS measurements.

For bulk-rock samples, which contain a much more complex matrix than pure zircon or baddeleyite, a two-step separation was used. Two replicate bulk-rock dissolutions of FC-1 were first processed using a slightly modified version of the tetraoctyl-1,5-diglycolamide (TODGA) separation method of (50). This purification approach has been shown to be very efficient at separating Ti (originally designed for this purpose) as well as Fe, Mo, and Ru, which must be removed to minimize isobaric interferences in the Zr mass range (48, 51). The only modification we applied to the method of (50) was to double the volumes of the 12 M HNO3 + 1 wt % H2O2 and 3 M HNO3 steps to ensure complete removal of Mo, Ti, and Fe from the final Zr-Hf cuts. Once high-purity Zr-Hf solutions were obtained with this method, Zr was further purified using the same Ln-Spec approach described above for our zircon and baddeleyite purifications.

Zr isotopic measurements were conducted on a Nu Plasma II MC-ICP-MS fitted with an enhanced sensitivity interface at the Massachusetts Institute of Technology. Analyses were conducted in dry plasma mode using a Cetac Aridus II desolvator nebulizer for sample introduction to the MC-ICP-MS. The average instrumental sensitivity using this setup was around 572 V/ppm of total Zr. All samples and bracketing standards were adjusted to a total Zr concentration of 60 ng/g, which resulted in total Zr beams of 30 to 35 V on average, or ~11 V for the most abundant isotope 90Zr (ca. 32% abundance in our double-spiked samples). Masses 90, 91, 92, 93, 94, 95, 96, and 98 were measured in static mode at 0.5 amu spacing in the Nu Plasma II collector block, allowing direct monitoring of all Zr isotopes and Mo interferences (masses 95 and 98). Each sample measurement consisted of an acid blank measurement followed by the sample, which was acquired over 50 cycles of 5-s integrations. On-peak-zero correction was done using the mean acid blank intensities before data processing to account for blank contribution as well as any “memory” effects from the Aridus II sample introduction system during the run. Each sample measurement was individually bracketed by measurements of the ZrNIST solution spiked at the same level as our samples and matched in concentration (60 ng/g) as well as acid matrix (0.59 M HNO3 + 0.28 M HF).

Data reduction of double-spike data was done by simultaneously solving the system of nonlinear equations of the formRmeas=[pRSpike+(1p)RStd(MxMn)α](MxMi)βwhere Rmeas, RStd, and RSpike are the isotopic ratio of the sample, reference standard, and double spike, p is the proportion of spike in the sample-spike mixture, Mx is the absolute mass of the isotope of interest, Mn is the absolute mass of the normalizing isotope, α is the natural fractionation factor between the sample and the reference standard, and 𝛽 is an instrumental fractionation factor (49, 52, 53). Because there are three unknowns in this formulation, at least four isotopes need to be measured and three equations need to be inverted simultaneously to solve the system. In the case of Zr, which has five stable isotopes, four ratios can be calculated. Data were reduced using a minimization approach implemented in Mathematica, taking into account all ratios (i.e., 91/90Zr, 92/90Zr, 94/90Zr, and 96/90Zr) with different weighs being assigned to each based on their associated uncertainty. For completeness, “exact” solutions to the system using only three isotope ratios were also performed in two ways: (i) using the 91/90Zr, 92/90Zr, and 96/90Zr ratios and (ii) using the 91/90Zr, 94/90Zr, and 96/90Zr ratios. This procedure was done as an internal consistency check, and in all instances, we found the exact solving and four-ratio minimization approaches to yield identical results within uncertainty. The data reported here (Table 1) correspond to results obtained using the four-ratio minimization approach.

After double-spike inversion, the spike-stripped isotopic ratios of samples and their bracketing standards were used to calculate delta offsets relative to ZrNIST using delta notation asδ9x/90ZrNIST=[(Z9x/90rSmp/Z9x/90rStd)1]103

Zirconium isotopic compositions in the main text are only discussed using δ94/90ZrNIST offsets, but for completeness, δ91/90ZrNIST, δ92/90ZrNIST, and δ96/90ZrNIST values are also presented in Table 1, demonstrating the mass dependency of the isotopic compositions we determined. The uncertainty assigned to each individual δ9x/90ZrNIST determination is the external reproducibility (at 2σ) of the spiked ZrNIST measurements from each run, which in all cases was similar in magnitude or slightly larger than the internal uncertainty determined from counting statistics.

To assess the effects that isobaric interferences from residual Mo in purified solutions have on high-precision Zr stable isotope measurements made on the Nu Plasma II at MIT, a doping test was performed by adding variable amounts of natural Mo to double-spiked NIST Zr standard solutions. Doped solutions in the range from 0.1 to 3.4% in Mo/Zr (expressed as atomic ratios of total metal content) were processed using the same approach described above for data collection and reduction. Results are presented graphically in fig. S6, where it is shown that at these doping ratios Mo isobaric interferences have a close-to-linear effect on δ94/90Zr values according to the relationshipOffset δ94/90ZrNIST=63.4xwhere x represents the total Mo/Zr atomic ratio measured in the sample. Using this relationship, we calculate the offsets that residual Mo isobaric interferences would cause to our reported FC-1 zircon, baddeleyite, and whole-rock data, as a function of the Mo/Zr ratios determined by MC-ICP-MS (Table 1 and fig. S7). The lack of any systematic relationships between δ94/90Zr and Mo/Zr in our measured samples (fig. S7), in addition to the very low measured Mo/Zr ratios (Table 1), demonstrates that the extreme δ94/90Zr variations reported here are not an artifact of unaccounted Mo interferences. Calculated offsets are reported in tabulated form in Table 1, but because the magnitude of such corrections is inconsequential for our conclusions (i.e., the vast majority are smaller than the reported analytical precision), Zr stable isotope results reported in Table 1 and discussed throughout the text are uncorrected for residual Mo interference.

Zircon Zr and Hf concentration measurements

Aliquots from each zircon dissolution were taken for trace element analyses before chromatographic purification and were measured by solution Q-ICP-MS using methods similar to those in (29, 42). A gravimetric standard solution was prepared at the UR-LIG using high-purity single element standards obtained from SPEX CertiPrep (, mixed in appropriate proportions to resemble those of a typical zircon (e.g., Zr/Hf ≈ 42). A series of dilutions of this concentrated standard solution were prepared gravimetrically to cover a range of Zr concentrations between 5 and 1000 ng/g, using a 3% HNO3 + 0.2% HF acid matrix and containing 2 ng/g In as an internal tracer. After drying down the sample aliquots to chloride salts, these were redigested using the same acid matrix as our calibration curve standards (i.e., 3% HNO3 + 0.2% HF with 2 ng/g In) before Q-ICP-MS measurements. Gravimetric standards were analyzed first, and the mean count rates from 60 1-s integrations were used to construct a sensitivity calibration curve for each element on the Q-ICP-MS. After rinsing in blank acid and performing several acid blank measurements, samples were analyzed using the same acquisition routine as the calibration standards, and drift corrections were applied using the mean count rate of 115In in the samples relative to those measured in the calibration standards. Total elemental concentrations of Zr and Hf in our sample solutions were determined relative to the known concentrations of the calibration standards and were ratioed to obtain the Zr/Hf composition of individual zircon crystals reported in Table 1.

Zr isotopic composition of individual crystals with integrated growth histories

In a Rayleigh-type fractionation scenario, the products of a reaction become isolated from the system as the reaction proceeds, preventing re-equilibration with the source reservoir and allowing for large isotopic effects to develop as the removal of the element of interest increases. The analytical equations describing the evolution of the source, sink and instantaneous sink reservoir are well established (25)RSource(f)=R0fα1(1)RInst. sink(f)=R0 α fα1(2)RCumul. sink(f)=R0RSource(f) f(1f)(3)

In the above equations, Rx is the isotopic ratio of interest for each particular reservoir, R0 is the initial isotopic composition of the bulk system, a is the isotopic fractionation coefficient between the source and instantaneous sink, and f is the fraction of the denominator isotope (in the ratio R) remaining in the source reservoir. When the isotopic fractionation between the source reservoir and the instantaneous sink is small (i.e., a is close to unity), f is approximately equal to the total fraction of the element remaining in the source.

In magmatic systems undergoing cooling and crystallization, single crystals are unlikely to represent an instantaneous solid composition, particularly if a kinetic impediment prevents effective re-equilibration with the melt. Instead, each individual crystal is likely to grow over a protracted period of time (i.e., integration over a finite growth history), where the addition of new growth layers records the changing chemistry of the melt. As such, the isotopic composition of a bulk crystal cannot simply be described by Eq. 2, which assumes the crystal represents an instantaneous value of f (i.e., fInst). Rather, each crystal records a growth-integrated isotopic composition, whose crystallization began and finished at two different times characterized by two different f values: fInitial and fFinal.

In the general case where a crystal begins to form sometime during the Rayleigh distillation (i.e., fInitial ≤1), fInitial and fFinal are related by the fraction of Zr incorporated into the crystal over its growth history, relative to the amount of Zr available to that crystal in the surrounding melt at the onset of crystallization. These two values can be calculated using the following equationsfInitial=(α fInstα1×p[1(1p)α])1(α1)=((δInst1000+1)×p[1(1p)α])1(α1) (4)fFinal=fInitial(1p)(5)where fInst is the value of f calculated for a given bulk crystal based on its isotopic composition (δInst) assuming instantaneous crystallization (i.e., Eq. 2), and p is the fraction of Zr incorporated into the crystal relative to the amount available in the melt at the onset of crystallization. The derivation of these equations is described in detail in the Supplementary Materials.

A special case must be considered for crystals with isotopic compositions higher than that of the bulk source reservoir (i.e., bulk rock). For these samples, the minimum value of fFinal is limited by Eq. 3, and there is a limit on the amount of Zr integration (p) beyond which fInitial values larger than 1 are calculated (which have no physical meaning). In such cases, the relationship between fInitial and fFinal shown in Eq. 5 does not hold anymore, and the minimum value of fFinal is simply the f value of the cumulative sink curve (Eq. 3) that matches the bulk-crystal isotopic composition. Althoughthe resulting equation cannot be solved analytically for fFinal, a solution utilizing a simple polynomial fit was found(see details in the Supplementary Materials).

The above equations can be used to quantify the uncertainties on the relative crystallization sequence of samples formed by a Rayleigh distillation process as a function of the growth history assumed for the grains. These equations were used to construct the uncertainty envelopes shown in Fig. 3 and the accuracy of the analytical solutions was confirmed numerically.

Correction (4 May 2020): This article has been corrected to include analytical equations that quantify the effect of variable growth integration histories on the isotopic composition of bulk individual grains. An excel file containing the implementation of the new equations has been added to the Supplementary Materials.


Supplementary material for this article is available at

Supplementary Text

Fig. S1. Textural relations among modal phases and Zr-rich accessory phases in the FC-1 anorthositic cumulate.

Fig. S2. Textural relations between Zr-rich accessory phases and modal phases in the FC-1 anorthositic cumulate.

Fig. S3. U-Pb (Wetherill) concordia diagram for FC-1 zircon (blue) and baddeleyite (green) single crystals dated by (±CA)-ID-TIMS.

Fig. S4. Apparent 206Pb/238U and 207Pb/206Pb dates of zircon and baddeleyite single-crystals from FC-1, obtained by (±CA)-ID-TIMS

Fig. S5. Elution curve of the Ln-Spec ion-exchange chemistry optimization used to purify Zr from zircon and baddeleyite samples.

Fig. S6. Doping test performed to quantify the offsets in δ94/90Zr induced by residual Mo isobaric interferences.

Fig. S7. δ94/90Zr values of zircon, baddeleyite, and bulk-rock aliquots of FC-1 gabbro as a function of residual Mo interference in the measured solutions.

Fig. S8. Results of the Monte Carlo simulations exploring the predicted distribution of bulk zircon δ94/90Zr values for FC-1 as a function of different isotopic fractionation coefficients.

Fig. S9. Plot of isotopic composition versus fraction of element remaining in the source reservoir showing a Rayleigh distillation evolution of the source, sink and instantaneous sink reservoirs (black curves).

Fig. S10. Plot showing the relationship between the f value of the instantaneous sink (fInst) and the cumulative sink (fCumulative) such that both reservoirs have the same isotopic composition.

Table S1. Results of U-Pb isotopic analyses of FC-1 zircon and baddeleyite crystals using ID-TIMS (isotope dilution–thermal ionization mass spectrometry).

Table S2. Assignment of f values for zircon and baddeleyite from FC-1, and uncertainty limits imposed by different magnitudes of fractional growth integration history according to Fig. 3.

Data file S1. General Rayleigh distillation accounting for crystal growth and quantifying the uncertainties in f.

References (54, 55)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank D. McGee for access to the MIT Nu Plasma II and J. DesOrmeau for the scanning electron microscopy imaging at U. Nevada in Reno. We thank E. Bloch and M. Méheut for comments to an earlier version of the manuscript. G. Gaetani, H. Marshall, and three anonymous reviewers provided comments that greatly helped improve the manuscript. Funding: This research was supported by NSF-EAR grants 1823748 (to M.I.-M.) and 1824002 (to F.L.H.T.), and startup funds to M.I.-M. provided by U. Rochester. Author contributions: M.I.-M. and F.L.H.T. designed the study, developed all laboratory methods, performed the measurements, interpreted the results, and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Aliquots from rock and mineral samples used in this study are stored at U. Rochester and available from M.I.-M. upon request.

Stay Connected to Science Advances

Navigate This Article