Figures

top

Letter

There is an unprecedented increase in use of non-traditional stable isotopes fuelled by analytical advances in geochemistry and cosmochemistry (Watkins et al., 2017

Watkins, J.M., DePaolo, D.J., Watson, E.B. (2017) Kinetic Fractionation of Non-Traditional Stable Isotopes by Diffusion and Crystal Growth Reactions. Reviews in Mineralogy and Geochemistry 82, 85–125. https://doi.org/10.2138/rmg.2017.82.4

). After a certain system becomes analytically possible, development proceeds in the following sequence: 1) to find measurable ranges of variations, 2) to interpret their origin, 3) to find exciting applications that solve urgent geological problems (Aarons et al., 2021

Aarons, S.M., Johnson, A.C., Rader, S.T. (2021) Forming Earth’s Continental Crust: A Nontraditional Stable Isotope Perspective. Elements 17, 413–418. https://doi.org/10.2138/gselements.17.6.413

), 4) to critically appraise the new results, and finally 5) “settling” the new isotopic system as “conventional” with well-established measurement protocols and range of uses. Zircon is a unique mineral that is being used extensively for U-Th-Pb geochronology, trace element variations, O, Li and radiogenic Hf variations, and as such it has seen many of these “rises” and “falls”. As we demonstrate below, this may characterise their Zr isotope variations, a subject of much recent research.

Zirconium has five stable isotopes and the 94/90 ratio expressed as δ^94/90Zr relative to a standard is commonly used to understand its variations (Inglis et al., 2018

Inglis, E.C., Creech, J.B., Deng, Z., Moynier, F. (2018) High-precision zirconium stable isotope measurements of geological reference materials as measured by double-spike MC-ICPMS. Chemical Geology 493, 544–552. https://doi.org/10.1016/j.chemgeo.2018.07.007

, 2019

Inglis, E.C., Moynier, F., Creech, J., Deng, Z., Day, J.M.D., Teng, F.-Z., Bizzarro, M., Jackson, M., Savage, P. (2019) Isotopic fractionation of zirconium during magmatic differentiation and the stable isotope composition of the silicate Earth. Geochimica et Cosmochimica Acta 250, 311–323. https://doi.org/10.1016/j.gca.2019.02.010

). This ratio can be measured with good sub-per mille precision in rocks, major minerals, and in situ inside zircons by either SIMS or LA-ICP-MS methods, permitting core to rim, or face to core profiles (Ibañez-Mejia and Tissot, 2019

Ibañez-Mejia, M., Tissot, F.L.H. (2019) Extreme Zr stable isotope fractionation during magmatic fractional crystallisation. Science Advances 5, eaax8648. https://doi.org/10.1126/sciadv.aax8648

; Zhang et al., 2019

Zhang, W., Wang, Z., Moynier, F., Inglis, E., Tian, S., Li, M., Liu, Y., Hu, Z. (2019) Determination of Zr isotopic ratios in zircons using laser-ablation multiple-collector inductively coupled-plasma mass-spectrometry. Journal of Analytical Atomic Spectrometry 34, 1800–1809. https://doi.org/10.1039/C9JA00192A

; Guo et al., 2020

Guo, J.-L., Wang, Z., Zhang, W., Moynier, F., Cui, D., Hu, Z., Ducea, M.N. (2020) Significant Zr isotope variations in single zircon grains recording magma evolution history. Proceedings of the National Academy of Sciences 117, 21125–21131. https://doi.org/10.1073/pnas.2002053117

; Tompkins et al., 2020

Tompkins, H.G.D., Zieman, L.J., Ibañez-Mejia, M., Tissot, F.L.H. (2020) Zirconium stable isotope analysis of zircon by MC-ICP-MS: methods and application to evaluating intra-crystalline zonation in a zircon megacryst. Journal of Analytical Atomic Spectrometry 37, 1167–1186. https://doi.org/10.1039/C9JA00315K

). Inglis et al. (2019)

Inglis, E.C., Moynier, F., Creech, J., Deng, Z., Day, J.M.D., Teng, F.-Z., Bizzarro, M., Jackson, M., Savage, P. (2019) Isotopic fractionation of zirconium during magmatic differentiation and the stable isotope composition of the silicate Earth. Geochimica et Cosmochimica Acta 250, 311–323. https://doi.org/10.1016/j.gca.2019.02.010

demonstrated a 0.5 ‰ positive shift in a single silicic magmatic suite of Hekla volcano (Iceland), which the authors interpreted as the result of igneous fractionations of isotopically light zircons. At the same time, Ibañez-Mejia and Tissot (2019)

Ibañez-Mejia, M., Tissot, F.L.H. (2019) Extreme Zr stable isotope fractionation during magmatic fractional crystallisation. Science Advances 5, eaax8648. https://doi.org/10.1126/sciadv.aax8648

studied zircons in the Duluth gabbro and found most of them isotopically heavier than the melt from which they crystallised (thus discovering an opposite trend of Zr isotopic fractionation), yet they also found a few extremely low δ^94/90Zr grains, with the range covering almost 5 ‰. Guo et al. (2020)

Guo, J.-L., Wang, Z., Zhang, W., Moynier, F., Cui, D., Hu, Z., Ducea, M.N. (2020) Significant Zr isotope variations in single zircon grains recording magma evolution history. Proceedings of the National Academy of Sciences 117, 21125–21131. https://doi.org/10.1073/pnas.2002053117

measured δ^94/90Zr in individual zircon grains and found lighter cores and heavier rims, while Tompkins et al. (2020)

Tompkins, H.G.D., Zieman, L.J., Ibañez-Mejia, M., Tissot, F.L.H. (2020) Zirconium stable isotope analysis of zircon by MC-ICP-MS: methods and application to evaluating intra-crystalline zonation in a zircon megacryst. Journal of Analytical Atomic Spectrometry 37, 1167–1186. https://doi.org/10.1039/C9JA00315K

found no isotopic differences in a large zircon megacryst that crystallised from carbonatitic magma. Given the large δ^94/90Zr ranges of several per mille discovered, Tian et al. (2021)

Tian, S., Moynier, F., Inglis, E.C., Rudnick, R.L., Huang, F., Chauvel, C., Creech, J.B., Gaschnig, R.M., Wang, Z., Guo, J.-L. (2021) Zirconium isotopic composition of the upper continental crust through time. Earth and Planetary Science Letters 572, 117086. https://doi.org/10.1016/j.epsl.2021.117086

suggested that Zr isotopes can be used to monitor continental crust evolution since the Archean, and also to provide a new tool to investigate the origin of the earliest crust. The field is rapidly developing and thus requires better understanding of what controls Zr isotopic variations in nature. This question has been recently addressed by Chen et al. (2020)

Chen, X., Wang, W., Zhang, Z., Nie, N.X., Dauphas, N. (2020) Evidence from Ab Initio and Transport Modeling for Diffusion-Driven Zirconium Isotopic Fractionation in Igneous Rocks. ACS Earth and Space Chemistry 4, 1572–1595. https://doi.org/10.1021/acsearthspacechem.0c00146

, Méheut et al. (2021)

Méheut, M., Ibañez-Mejia, M., Tissot, F.L.H. (2021) Drivers of zirconium isotope fractionation in Zr-bearing phases and melts: The roles of vibrational, nuclear field shift and diffusive effects. Geochimica et Cosmochimica Acta 292, 217–234. https://doi.org/10.1016/j.gca.2020.09.028

, and Tissot and Ibañez-Mejia (2021)

Tissot, F.L.H., Ibañez-Mejia, M. (2021) Unlocking the Single-Crystal Record of Heavy Stable Isotopes. Elements 17, 389–394. https://doi.org/10.2138/gselements.17.6.389

, among others. Theoretical ab initio modelling (Méheut et al., 2021

Méheut, M., Ibañez-Mejia, M., Tissot, F.L.H. (2021) Drivers of zirconium isotope fractionation in Zr-bearing phases and melts: The roles of vibrational, nuclear field shift and diffusive effects. Geochimica et Cosmochimica Acta 292, 217–234. https://doi.org/10.1016/j.gca.2020.09.028

) returns very small <0.1 ‰ equilibrium melt-zircon fractionations at realistic high temperature (>600 °C) conditions, when zircon is able to crystallise from magma, down to the low temperature of water saturated pegmatites. These authors, Tissot and Ibañez-Mejia (2021)

Tissot, F.L.H., Ibañez-Mejia, M. (2021) Unlocking the Single-Crystal Record of Heavy Stable Isotopes. Elements 17, 389–394. https://doi.org/10.2138/gselements.17.6.389

and Aarons et al. (2021)

Aarons, S.M., Johnson, A.C., Rader, S.T. (2021) Forming Earth’s Continental Crust: A Nontraditional Stable Isotope Perspective. Elements 17, 413–418. https://doi.org/10.2138/gselements.17.6.413

suggested that kinetic effects likely dominate the Zr isotope variations in nature, via faster diffusion of light isotopes of zirconium during zircon crystallisation from magmas, and in front of crystallising major minerals.

We here present the result of numerical modelling of Zr isotopes and Zr/Hf ratios during zircon crystallisation, based on the further development of the Bindeman and Melnik (2016)

Bindeman, I.N., Melnik, O.E. (2016) Zircon Survival, Rebirth and Recycling during Crustal Melting, Magma Crystallization, and Mixing Based on Numerical Modelling. Journal of Petrology 57, 437–460. https://doi.org/10.1093/petrology/egw013

and Melnik and Bindeman (2018)

Melnik, O.E., Bindeman, I.N. (2018) Modeling of trace elemental zoning patterns in accessory minerals with emphasis on the origin of micron-scale oscillatory zoning in zircon. American Mineralogist 103, 355–368. https://doi.org/10.2138/am-2018-6182

models. In short, the model is realised in spherical coordinates, in which zircon is growing from the centre and other minerals are forming on the outer cell boundary (Fig. 1) by following a T-X (% crystals) phase diagram. Upon cooling and solidification, zircon grows outward and other minerals advance inward, decreasing the proportion of melt in the cell. The melt compositional parameter, M factor (Harrison and Watson, 1983

Harrison, T.M., Watson, E.B. (1983) Kinetics of zircon dissolution and zirconium diffusion in granitic melts of variable water content. Contributions to Mineralogy and Petrology 84, 66–72. https://doi.org/10.1007/BF01132331

), that determines zircon solubility is taken into account via its temperature-compositional dependence as in Bindeman and Melnik (2016)

Bindeman, I.N., Melnik, O.E. (2016) Zircon Survival, Rebirth and Recycling during Crustal Melting, Magma Crystallization, and Mixing Based on Numerical Modelling. Journal of Petrology 57, 437–460. https://doi.org/10.1093/petrology/egw013

. Zirconium, its isotopes, and other trace element diffusion coefficients are parameterised as a function of temperature and water content based on Zhang et al. (2010)

Zhang, Y., Ni, H., Chen, Y. (2010) Diffusion Data in Silicate Melts. Reviews in Mineralogy and Geochemistry 72, 311–408. https://doi.org/10.2138/rmg.2010.72.8

. Since each zirconium isotope has a large abundance, for this work we have rewritten saturation boundary conditions at the zircon/melt boundary, so they add together to total [Zr] at C_sat at each T and M factor (Fig. 1). Mass conservation in the entire melt cell is obeyed (Fig. S-1).

Full size image

We assume that different isotopes of zirconium diffuse with different rates in proportion to their mass, following a mass dependent relationship of diffusion coefficients via a kinetic law with a beta exponent parameter β = 0.05 (Watkins et al., 2017

Watkins, J.M., DePaolo, D.J., Watson, E.B. (2017) Kinetic Fractionation of Non-Traditional Stable Isotopes by Diffusion and Crystal Growth Reactions. Reviews in Mineralogy and Geochemistry 82, 85–125. https://doi.org/10.2138/rmg.2017.82.4

). Variations in this parameter within a reasonable ±0.02 range are less important than isotopic effects unravelled due to, for example, cooling rates and other crystallisation conditions described below. The outer boundary on which major minerals crystallise has a condition of Zr (and Hf) partitioning in accordance with their Zr partition coefficient, K_d, set by the user. As Zr and Hf are highly incompatible in most minerals, K_d is <<1; however, we also experimented by making this K_d around 1 for the case of amphibole, ±micro-zircon or sphene inclusion co-precipitation. We also implemented and monitored the behaviour of U, Th, Y, P, and Ti (Fig. S-4).

We performed a series of computations (Figs. 2, 3 and S-3 to S-5) of zircon crystallisation from magma during linear cooling from 900 to 650 °C by varying 1) zircon crystallisation duration increasing from 150 years to 0.5 Myr, 2) zircon crystallisation in low and high temperature ranges, achievable by increasing and decreasing Zr concentrations in the initial melt, 3) co-crystallisation of minerals with different Zr partition coefficients, and 4) crystallisation of magmas with different water contents. These parameters capture most of the conditions in nature. Figure 2 presents diffusion profiles through the melt cell initially containing 100 ppm Zr, showing Zr isotopic and Zr/Hf distribution in a typical run at realistic geologic parameters. Notice evolving sigmoidal shapes explained by the crystallisation of zircon and major minerals. We additionally consider Hf concentration and monitor Zr/Hf ratio. Hf is an element twin for zirconium with an identical ionic radius, 4+ charge, and thus extremely close chemical properties (partitioning with K_d ≈ 10³–10⁴). However, because of its greater mass (174–180), Hf is expected to diffuse more slowly in melts, similar to the heaviest Zr isotope.

Full size image

Full size image

The main result of this study is that we observe positive, negative, and inverse δ^94/90Zr values and Zr/Hf ratios (Figs. 3, 4). Notice that this diversity of isotopic results can explain the range, as well as both positive and negative δ^94/90Zr values observed in nature. Furthermore, a simultaneous change of these parameters and more complicated cooling and heating paths can lead to further complexity of Zr isotope behaviour, which can enhance or diminish gradual trends in Figure 3. Some general statements can also be made. Given the faster diffusion of ⁹⁰Zr, zircon is predicted to decrease in δ^94/90Zr from core to rim upon its sole growth in a cell. These effects are very small (<0.1 ‰) to almost non-existent at high T, high H₂O and large melt cells when diffusion coefficients for both ⁹⁰Zr and ⁹⁴Zr are high, leading to lack of any δ^94/90Zr zoning. Melt cells around zircons evolve toward higher δ^94/90Zr values (Fig. 2). Crystallisation of the other minerals on the outer boundary (set here not to discriminate Zr isotopes) leads to piling up (bulldozing) effects and higher δ^94/90Zr values and Zr/Hf ratios due to effects of zircon crystallisation on the inside. It is thus expected that major minerals too, will have diversity of δ^94/90Zr and Zr/Hf ratios due to zircon co-crystallisation despite having overall low Zr and Hf concentrations. Therefore, simultaneous crystallisation of zircon and major minerals leads to a complex δ^94/90Zr distance profile (Fig. 2) and diverse δ^94/90Zr isotopic trends (Fig. 2).

Full size image

In nature, after zircon saturation, both Zr and Hf decrease in melts, but Hf is decreasing less, and thus the Zr/Hf ratio in rocks decreases by a factor of 2–3 from basalt to rhyolites with the progress of igneous differentiation (Claiborne et al., 2006

Claiborne, L.L., Miller, C.F., Walker, B.A., Wooden, J.L., Mazdab, F.K., Bea, F. (2006) Tracking magmatic processes through Zr/Hf ratios in rocks and Hf and Ti zoning in zircons: An example from the Spirit Mountain batholith, Nevada. Mineralogical Magazine 70, 517–543. https://doi.org/10.1180/0026461067050348

). Figure 4 plots the Zr/Hf ratio in zircons at different conditions and shows them to be a strong function of crystallisation parameters. During a long duration of crystallisation, in melts with greater water concentrations, Zr and Hf diffusion toward zircon’s growing boundary are near equilibrium, generating nearly identical Zr/Hf profiles within a zircon. If the crystallisation of zircon is rapid, from water-poor melt, with slow overall diffusion, and fast growth rates, then differences in Zr and Hf diffusion coefficients create kinetic disequilibrium conditions, and the zircon develops positive, negative, and inverse Zr/Hf profiles, with secondary kinks at the end of crystallisation and toward the zircon rim explained by co-crystallisation of Zr-and Hf-poor phases on the outside of the cell (Fig. 1).

Overall, the Zr/Hf chemical ratio follows the δ^94/90Zr trends (compare Fig. 4 to Fig. 3), positive at fast crystallisation and negative at slow. Both δ^94/90Zr and Zr/Hf ratios in growing zircons become lower for greater zircon growth rates, where the diffusive boundary layer around zircon favours light isotopes and Zr over Hf. Again, this is realised with greater rates of cooling, lower water, larger cells, and lower Zr concentrations. In our previous model targeting element partitioning into zircon, Hf concentration decreases toward the rim as a strong function of the rate of cooling (see Fig. 7 in Melnik and Bindeman, 2018

Melnik, O.E., Bindeman, I.N. (2018) Modeling of trace elemental zoning patterns in accessory minerals with emphasis on the origin of micron-scale oscillatory zoning in zircon. American Mineralogist 103, 355–368. https://doi.org/10.2138/am-2018-6182

), matching natural observations (Claiborne et al., 2006

Claiborne, L.L., Miller, C.F., Walker, B.A., Wooden, J.L., Mazdab, F.K., Bea, F. (2006) Tracking magmatic processes through Zr/Hf ratios in rocks and Hf and Ti zoning in zircons: An example from the Spirit Mountain batholith, Nevada. Mineralogical Magazine 70, 517–543. https://doi.org/10.1180/0026461067050348

, 2018

Claiborne, L.L., Miller, C.F., Gualda, G.A.R., Carley, T.L., Covey, A.K., Wooden, J.L., Fleming, M.A. (2018) Zircon as Magma Monitor: Robust, Temperature-Dependent Partition Coefficients from Glass and Zircon Surface and Rim Measurements from Natural Systems. In: Moser, D.E., Corfu, F., Darling, J.R, Reddy, S.M., Tait, K. (Eds.) Microstructural Geochronology: Planetary Records Down to Atom Scale. Geophysical Monograph 232, American Geophysical Union, Washington, D.C., John Wiley and Sons, Inc., Hoboken, NJ, 1–33. https://doi.org/10.1002/9781119227250.ch1

). It is expected that Hf stable isotopes (¹⁷⁶Hf/¹⁷⁷Hf ratio) will behave in a similar manner to ⁹⁴Zr/⁹⁰Zr isotopes during zircon crystallisation; Hf is a highly compatible element in zircon (K_d > 1000) and thus isotopic preferences for isotopically light ¹⁷⁶Hf are set at the zircon-melt interface in the same way as for Zr isotopes. Computed affects (Fig. S-5) are smaller than 1 ‰, but because Hf isotopes are measured in ɛ(¹⁷⁶Hf/¹⁷⁷Hf) units (which are 10× the δ scale), these are not negligible in affecting primary radiogenic Hf ratios. Supplementary Figure S-4 plots other trace element ratios that are easy to compute within the present model. The behaviour of highly compatible elements in zircon (U, Th, Y, P) is comparable to that of Hf (with piling up and inversions at fast crystallisation) with some nuances related to their diffusion coefficients and K_d.

Incompatible trace elements, such as Ti behave much more simply without inversions. Considering isotopes of other incompatible trace elements that partition in zircons, e.g., ¹⁴³Nd/¹⁴⁴Nd and ⁷Li/⁶Li, requires further development of the present model, but since these are incompatible, no kinetic isotope effects are expected due to zircon growth. Crystallisation of major minerals with different K_d on the outer cell boundary may further influence zircon growth and its element behaviour, however, due to the small size of zircon and the cubic relation of melt cell volume vs. its remaining radius, only very advanced degrees of crystallisation (>90 %) would produce piling up effects reflected in zircons.

The main result of the present work is a strong dependency of δ^94/90Zr and Zr/Hf ratios on crystallisation parameters. The software that we offer as a part of this paper (https://github.com/crystalworkshop/isotop_diffusion) can be used by practitioners to interrogate their specific situations of magmatic cooling and zircon growth, including other trace elements (Fig. S-4). As δ^94/90Zr will soon become an easy to obtain parameter in zircon research for a wide number of users (similar to how Ti in zircon, for example, is now used), the present approach will help to understand variations in both δ^94/90Zr and Zr/Hf ratios. Specific conditions of magmatic crystallisation, as well as the presence or absence of Rayleigh fractionation, in which an incrementally different δ^94/90Zr zircon or melt are removed or added to the system (such as in the classic example of Hekla) will explore the wealth of natural situations. A few general statements as a result of the modelling can be mentioned: higher temperature of crystallisation, slower rates and a longer duration of zircon growth, higher water contents of magma, and higher K_d(Zr) of co-crystallising phase all lead to minimal isotopic differences; the opposite processes lead to initially more negative δ^94/90Zr, which then invert with values spanning several per mille.

The present model provides an easy test for two prevailing end member models of silicic magma petrogenesis: 1) slow crystallisation with melt segregation in large mushy reservoirs vs. 2) rapid crystallisation in individual dike and sills, accreting during batholith/magma chamber construction, followed by storage with limited zircon core recrystallisation. The first scenario should result in no core to rim δ^94/90Zr zoning, while the second case should result in significant δ^94/90Zr core to rim variations (and Zr/Hf, plus other trace elements zoning). Trace element concentrations and ratios (e.g., Figs. 4, S-4) should also reflect these processes, and our model provides a platform for investigating their behaviours in zircons from different magmatic regimes and throughout geologic history.

top

Data and Materials Availability

The MATLAB Code and User Guide are available from GitHub: https://github.com/crystalworkshop/isotop_diffusion.

top

Acknowledgements

We thank Francois Tissot and Mauricio Ibañez-Mejia for constructive reviews and Helen Williams for editorial comments. Supported by NSF Grant EAR1833420.

Editor: Helen Williams

top

References

top

Supplementary Information

The Supplementary Information includes:

Mathematical Model

MATLAB Code User Guide

Figures S-1 to S-4

Supplementary Information References

Download the Supplementary Information (PDF)