Bénédicte Coulm
email: benedicte.coulm@aphp.fr

Cécile Henquell

Praticien

Sous-dochead La périnatalité à l'épreuve du Covid-19 Conséquences de l'infection à Sars-CoV-2 chez les femmes enceintes et les nouveau-nés

Keywords: 19, foetus, grossesse, infection, risque foetal, risque maternel

L'épidémie de Covid-19 survenue au premier semestre 2020 a fait naître de nombreuses craintes vis-à-vis des conséquences de l'infection virale chez les femmes enceintes et leurs nouveau-nés. Au début de l'été, les données disponibles ne mettaient pas en évidence de surrisque d'infection et de complications chez les femmes enceintes par rapport à la population générale. En l'absence d'atteinte maternelle sévère, celles concernant le foetus et le nouveau-né étaient également rassurantes.

responsable du Covid-19, a été à l'origine d'une crise sanitaire sans précédent. La susceptibilité des femmes enceintes à certains types d'infections et l'expérience d'autres épidémies, tant liées aux coronavirus qu'à d'autres virus respiratoires, comme ceux de la grippe, ont engendré de nombreux questionnements quant aux risques encourus par la mère et l'enfant à naître. Les publications sur ce sujet évoluent de manière permanente et il n'est ainsi pas aisé, au jour de rédaction de cet article, d'évaluer de manière détaillée l'impact de cette infection virale lorsqu'elle survient durant une grossesse. Les données publiées concernant les femmes enceintes infectées par le Sars-CoV-2 au premier ou au deuxième trimestre de grossesse sont rares, et les publications portent principalement sur des infections maternelles au troisième trimestre. Nous avons toutefois souhaité présenter les principaux éléments pouvant être dégagés des données disponibles 1 .

T1 Une émergence d'un nouveau virus et d'une nouvelle maladie

TEG1 Le 31 décembre 2019, la Chine informait l'Organisation mondiale de la santé (OMS) de la survenue de cas groupés de pneumonie d'origine inconnue chez des personnes ayant été exposées à des animaux vivants vendus sur un marché de la ville de Wuhan [START_REF]Organisation mondiale de la santé. Flambée de maladie à coronavirus[END_REF]. TEG1 Le 11 mars 2020, l'OMS a qualifié l'épidémie de Covid-19 de pandémie, du fait de sa diffusion mondiale [START_REF]Organisation mondiale de la santé. Flambée de maladie à coronavirus[END_REF]. En date du 3 juin 2020, 6 348 900 cas avaient été confirmés dans le monde depuis le 31 décembre 2019 et 380 810 décès avaient été comptabilisés. En France, depuis les premières identifications, le 24 janvier 2020, 151 677 cas confirmés (figure 1) et 29 021 décès étaient recensés [START_REF]COVID-19[END_REF].

TEG1 Comme les autres coronavirus, le Sars-CoV-2 provoque des infections des voies respiratoires, du tube digestif et du système nerveux. Dans sa forme classique, la plus courante (80 % des personnes atteintes), le Covid-19 est une infection respiratoire aiguë de type syndrome grippal. Toutefois, des tableaux beaucoup plus sévères de pneumopathie avec syndrome de détresse respiratoire aiguë ou défaillance multiviscérale sont observés, majoritairement chez les sujets âgés et les personnes présentant des comorbidités (diabète, obésité, hypertension, maladies chroniques respiratoires). TEG1 L'infection est principalement propagée par les gouttelettes de salive (toux, éternuements, etc.) et par contact du visage avec des mains souillées, d'où l'importance des gestes barrières [START_REF] Aumeran | Prévenir le risque infectieux en milieu de soins grâce aux mesures barrières[END_REF], particulièrement le port d'un masque et le lavage régulier des mains. La période d'incubation est de deux à quatorze jours.

T1 Un surrisque chez les femmes enceintes ? TEG1 Les modifications physiologiques liées à la grossesse, d'ordres immunologique et cardiorespiratoire, rendent les femmes enceintes plus vulnérables aux infections respiratoires et à leur évolution vers des formes sévères infectieuses. Dans le cadre de la grippe, ce risque de complications, à l'origine d'hospitalisations, augmente à mesure de l'avancement de la grossesse, passant d'un facteur 1,7 à 7,9 [START_REF] Dodds | Impact of influenza exposure on rates of hospital admissions and physician visits because of respiratory illness among pregnant women[END_REF].

TEG1 Cette vulnérabilité avait été mise en évidence lors de précédentes épidémies liées à des coronavirus (Sars-CoV, Mers-CoV), mais les données publiées concernant le Covid-19 sont plutôt rassurantes : il n'existerait pas de surrisque d'infection ni de complications chez les femmes enceintes par rapport à la population générale [START_REF] Yu | Clinical features and obstetric and neonatal outcomes of pregnant patients with COVID-19 in Wuhan, China: A retrospective, single-centre, descriptive study[END_REF][START_REF] Chen | Clinical characteristics and intrauterine vertical transmission potential of Covid-19 infection in nine pregnant women: A retrospective review of medical records[END_REF]. Quelques cas de décès maternels ont bien été recensés, mais chez des patientes présentant des comorbidités ou d'autres pathologies déjà sévères 1 .

T1 Quels risques pour les foetus ? TEG1 Les publications ne mettent pas en évidence de cas de transmission verticale maternofoetale du Covid-19 1 . De même, à quelques exceptions près, il n'a pas été observé de signes cliniques entrant dans le tableau du Covid-19 chez les enfants nés vivants de mères infectées.

TEG1 Chez les quelques nouveau-nés ayant présenté un tableau évocateur, confirmé par un test pharyngé positif au Sars-Cov-2, aucun élément ne permet d'établir si l'infection a eu lieu in utero ou après la naissance. En outre, tous les symptômes (lymphopénie, fièvre, pneumopathie, troubles digestifs) ont spontanément régressé en quelques jours sous antibiothérapie [START_REF] Zeng | Neonatal early-onset infection with SARS-CoV-2 in 33 neonates born to mothers with COVID-19 infection in Wuhan, China[END_REF][START_REF] Dong | Possible vertical transmission of SARS-CoV-2 from an infected mother to her newborn[END_REF][START_REF] Khan | Association of COVID-19 with pregnancy outcomes in health-care workers and general women[END_REF][START_REF] Alzamora | Severe COVID-19 during pregnancy and possible vertical transmission[END_REF][START_REF] Elshafeey | A systematic scoping review of COVID-19 during pregnancy and childbirth[END_REF].

TEG1 Chez quatre enfants bien portants négatifs pour le Sars-CoV-2, des immunoglobulines M ont été détectées à la naissance. L'interprétation de ces résultats est complexe du fait de la faible sensibilité des tests ; l'éventualité de faux positifs ne peut pas être écartée [START_REF] Zeng | Neonatal early-onset infection with SARS-CoV-2 in 33 neonates born to mothers with COVID-19 infection in Wuhan, China[END_REF][START_REF] Dong | Possible vertical transmission of SARS-CoV-2 from an infected mother to her newborn[END_REF].

TEG1 Des prélèvements effectués sur liquide amniotique, vaginaux et sang de cordon à la recherche du virus (une trentaine) sont tous revenus négatifs 1 , ce qui va dans le sens d'une absence de transmission maternofoetale du virus. [START_REF]Groupe de pathologie infectieuse pédiatrique[END_REF] concernant la prise en charge des nouveau-nés de mères présentant une infection confirmée ou en attente du résultat de dépistage du virus par réaction en chaîne par polymérase.

TEG1 Des recommandations ont été émises par la Société française de néonatalogie [12] et le Groupe de pathologie infectieuse pédiatrique

T1 Quid de la morbidité obstétricale ? TEG1 Par analogie avec d'autres infections virales respiratoires (grippe, mais aussi Sars-CoV et Mers-CoV), on peut émettre l'hypothèse qu'une infection maternelle sévère due au Sars-Cov-2 pourrait être à l'origine d'un plus grand nombre de naissances prématurées ou de retards de croissance intra-utérins.

TEG1 Sur une petite série de patientes infectées par le Covid-19 hospitalisées, des auteurs ont retrouvé un taux de prématurité de 41 %, mais ce résultat est difficile à interpréter car nous ne disposons pas d'informations sur le contexte de la naissance ; il est possible qu'il s'agisse d'une prématurité induite par un état maternel aggravé [START_REF] Mascio | Outcome of coronavirus spectrum infections (SARS, MERS, COVID-19) during pregnancy: A systematic review and meta-analysis[END_REF].

TEG1 Quelques rares morts foetales in utero ont également été rapportées, en lien avec un état maternel très péjoratif 1 . Les patientes présentaient toutes une infection sévère avec hypoxémie majeure.

T1 Que penser de l'allaitement ? TEG1 Les données disponibles sur l'allaitement et le Covid-19 sont très peu nombreuses, mais aucun cas de transmission de l'infection via l'allaitement maternel n'a été décrit 1 . Le risque de contamination de l'enfant semble essentiellement lié à la promiscuité avec sa mère infectée. Sur une vingtaine de prélèvements, le virus n'a été retrouvé dans le lait que dans un cas, positif vingt-quatre heures après la naissance, puis négatif au troisième jour. TEG1 Les différentes institutions internationales et sociétés savantes ont précisé que l'allaitement pouvait être poursuivi, même si la mère était infectée, à la condition de prendre soin de respecter quelques précautions standard : port d'un masque lors des contacts avec l'enfant ; lavage fréquent et soigneux des mains [START_REF]Propositions de la Société française de néonatalogie & de la Société française de pédiatrie concernant les nouveau-nés dans le contexte d'épidémie à[END_REF][START_REF]Questions-réponses sur la Covid-19 et l'allaitement[END_REF].

T1 Conclusion

Le risque d'une infection par le Covid-19 en cours de grossesse semble principalement reposer sur la gravité de l'atteinte maternelle 1 , ce qui est cohérent avec les observations décrites chez les femmes enceintes atteintes par d'autres virus, comme celui de la grippe. Il paraît important de prévenir autant que possible la contamination d'une femme enceinte, particulièrement en présence de comorbidités susceptibles d'entraîner une évolution vers une forme grave (obésité, pathologie sous-jacente notamment cardiovasculaire, etc.). Le respect des gestes barrières au quotidien est primordial, ainsi que la vaccination, lorsqu'un vaccin est disponible et que son administration n'est pas contre-indiquée en cours de grossesse. À quelques semaines du début de la campagne vaccinale contre la grippe en France, nous pouvons espérer que des leçons seront tirées de cette épidémie de Covid-19 et que les professionnels de santé et les femmes enceintes seront bien plus nombreux que les années passées à se faire vacciner.

Déclaration de liens d'intérêts

Les auteurs déclarent ne pas avoir de liens d'intérêts.

Note

1 Les données retranscrites sont celles qui étaient disponibles au moment de la rédaction de cet article.

Encadré 1 .

 1 Le Sars-CoV-2, un nouveau coronavirus Le Sars-CoV-2 est un virus à acide ribonucléique qui fait partie de la famille des coronavirus, comme le Sars-CoV (Severe Acute Respiratory Syndrome Coronavirus) et le Mers-CoV (Middle-East Respiratory Syndrome Coronavirus). Les personnes infectées présentent des symptômes caractéristiques : hyperthermie, fatigue, toux, difficultés respiratoires, dont la gravité est variable, ou encore anosmie. L'infection est principalement propagée par les gouttelettes de salive et par contact du visage avec des mains souillées. © B. Coulm

Figure 1 .

 1 Figure 1. Nombre cumulé de cas confirmés de Covid-19 (système d'information pour le suivi des victimes [SI-VIC] + trois laboratoires) et cas confirmés supplémentaires du Système d'information de dépistage [SI-DEP], par date de signalement à Santé publiqueFrance du 27 mars au 28 mai 2020 (données arrêtées le 28 mai 2020)[START_REF]COVID-19[END_REF].

TEG1 Dès le 7 janvier 2020, l'agent responsable de cette nouvelle maladie, dénommée

  Covid-19 pour Coronavirus Infectious Disease 2019, était identifié : il s'agissait d'un coronavirus, d'abord appelé 2019-nCoV. Le séquençage de son génome complet a rapidement montré son lien de parenté avec un autre coronavirus apparu en Chine en 2002, le Sars-CoV pour Severe Acute Respiratory Syndrome Coronavirus,

ce qui a conduit à sa dénomination officielle de Sars-CoV-2, pour Severe Acute Respiratory Syndrom Coronavirus- 2 [1] (encadré 1). TEG1 L'hypothèse la plus probable aujourd'hui permettant d'expliquer l'émergence de ce nouveau virus est celle

	d'une transmission à l'homme à partir d'un réservoir animal naturel,
	la chauve-souris, après un passage par un réservoir intermédiaire encore non formellement
	identifié. La comparaison, au niveau moléculaire, de coronavirus animaux et du Sars-CoV-2
	semble désigner, en tant qu'hôte intermédiaire possible, le pangolin, petit mammifère des
	régions tropicales.

© IMAGE POINT FR -LPN/BSIP


Merlin Méheut

Mauricio Ibañez-Mejia

François L H Tissot

Drivers of zirconium isotope fractionation in Zr-bearing phases and melts: The roles of vibrational, nuclear field shift and diffusive effects

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Zirconium isotopes have long been studied in extraterrestrial materials as they hold clues to the chronology of, and stellar contributions to, the early Solar System (e.g., [START_REF] Yin | Supernova Sources and the 92Nb-92Zr p-Process Chronometer[END_REF][START_REF] Akram | Zirconium isotope constraints on the composition of Theia and current Moon-forming theories[END_REF]. Recently, investigations of mass-dependent isotopic effects in igneous samples have revealed significant δ 94/90 Zr (i.e., [( 94 Zr/ /90 ZrSample/ 94 Zr/ /90 ZrStandard)-1]*10 3 ) variability at both the mineral [START_REF] Ibañez-Mejia | Extreme Zr stable isotope fractionation during magmatic fractional crystallization[END_REF][START_REF] Zhang | Determination of Zr isotopic ratios in zircons using laser-ablation multiple-collector inductively coupled-plasma mass-spectrometry[END_REF][START_REF] Guo | Significant Zr isotope variations in single zircon grains recording magma evolution history[END_REF] and bulk-rock scales [START_REF] Inglis | Isotopic fractionation of zirconium during magmatic differentiation and the stable isotope composition of the silicate Earth[END_REF][START_REF] Feng | High-precision stable zirconium isotope ratio measurements by double spike thermal ionization mass spectrometry[END_REF][START_REF] Tian | The zirconium stable isotope compositions of 22 geological reference materials, 4 zircons and 3 standard solutions[END_REF], suggesting Zr isotopes may be useful tracers of magmatic differentiation. Yet, the mechanism(s) driving the observed variability and their petrogenetic significance remain poorly understood. In particular, the direction and magnitude of equilibrium Zr stable isotope fractionation between Zr-bearing phases and silicate melts is unknown.

At thermodynamic equilibrium, heavy isotopes typically populate shorter and stiffer bonds with the lowest coordination number [START_REF] Schauble | Applying Stable Isotope Fractionation Theory to New Systems[END_REF][START_REF] Young | High-temperature equilibrium isotope fractionation of non-traditional stable isotopes: Experiments, theory, and applications[END_REF][START_REF] Blanchard | Equilibrium Fractionation of Non-traditional Isotopes: a Molecular Modeling Perspective[END_REF]. Since Zr coordination in the melt (i.e., coordination number VI) is expected to be lower than in baddeleyite (VII) or zircon (VIII) [START_REF] Farges | Structural environments of incompatible elements in silicate glass/melt systems: I. Zirconium at trace levels[END_REF][START_REF] Louvel | Zr complexation in high pressure fluids and silicate melts and implications for the mobilization of HFSE in subduction zones[END_REF], the first-order expectation is that these common Zr-rich phases would be isotopically light relative to the melt from which they crystallize. Such a scenario has recently been invoked by [START_REF] Inglis | Isotopic fractionation of zirconium during magmatic differentiation and the stable isotope composition of the silicate Earth[END_REF] to explain a positive correlation between δ 94/90 Zr and SiO2 wt. % in bulk volcanic rock samples from the Hekla Volcano, Iceland. In stark contrast, a detailed investigation [START_REF] Ibañez-Mejia | Extreme Zr stable isotope fractionation during magmatic fractional crystallization[END_REF] of single zircon and baddeleyite crystals from a well-characterized gabbroic igneous cumulate (FC-1 anorthositic gabbro from the Duluth Complex, USA) revealed that these Zr-rich phases are isotopically heavy (by up to ~1 ‰) relative to the starting melt composition (as represented by the bulk rock value). Adding to this conundrum, microbeam [START_REF] Kirkpatrick | Zr Isotope Variations in Zircon as an Indicator of Magmatic Differentiation[END_REF][START_REF] Zhang | Determination of Zr isotopic ratios in zircons using laser-ablation multiple-collector inductively coupled-plasma mass-spectrometry[END_REF] and solution [START_REF] Tompkins | Zirconium stable isotope analysis of zircon by MC-ICP-MS: methods and application to evaluating intra-crystalline zonation in a zircon megacryst[END_REF] δ 94/90 Zr data in zircon from several reference localities (e.g., Mud Tank, 91500) display homogeneous compositions, which is at odds with isotope fractionation as being driven by zircon crystallization. Recent LA-ICPMS data by [START_REF] Guo | Significant Zr isotope variations in single zircon grains recording magma evolution history[END_REF] further complicate the picture, reporting lighter composition in the core (compared to the rims) of zircon from the Gangdese Arc in Tibet, which the authors interpreted as indicative of negative and variable equilibrium fractionation factors between zircon and melt (from -0.12 to -0.45 ‰).

Here, we approach these conflicting observations from a theoretical standpoint: (i) We performed first-principles, ab initio calculations to quantify the direction and magnitude of Zr isotope fractionation between silicate melts and Zr-rich phases at thermodynamic equilibrium.

(ii) We explored the potential of Nuclear Field Shift effects for driving significant δ 94/90 Zr variations, and established a maximum allowable contribution of this effect to the published data. (iii) We evaluated the potential for diffusion-driven (i.e., kinetic) effects as drivers of Zr isotope fractionation. Altogether, these investigations place important constraints on the drivers of Zr isotope variability in Zr-bearing phases and liquids at magmatic temperatures, and highlight that diffusion in silicate melts may be the major mechanism of Zr isotopic fractionation.

Zirconium bonding environments in silicate melts and Zr-rich phases

Differences in the local bonding environment of Zr in silicate liquids and co-existing minerals are the main drivers of vibrational equilibrium isotope fractionation. The speciation of Zr in fluids at high-P-T and in variably hydrated silicate melts of different compositions is fairly well known from studies using molecular dynamics, X-ray absorption spectroscopy near-edge structure (XANES), extended X-ray absorption fine structure (EXAFS) and ab initio calculations [START_REF] Farges | Structural environments of incompatible elements in silicate glass/melt systems: I. Zirconium at trace levels[END_REF][START_REF] Galoisy | Evidence for 6-Coordinated Zirconium in Inactive Nuclear Waste Glasses[END_REF][START_REF] Montorsi | Molecular dynamics study of zirconia containing glasses[END_REF][START_REF] Wilke | Zircon solubility and zirconium complexation in H2O+Na2O+SiO2±Al2O3 fluids at high pressure and temperature[END_REF][START_REF] Louvel | Zr complexation in high pressure fluids and silicate melts and implications for the mobilization of HFSE in subduction zones[END_REF][START_REF] Wang | Coordination of Zr4+/Hf4+/Nb5+/Ta5+ in silicate melts: insight from first principles molecular dynamics simulations[END_REF]. These studies have found Zr in silicate melts to exhibit short-range order and to be mainly in octahedral (6-fold) coordination, with six O nearest neighbors and mean Zr-O bond distances (RZr-O) from 2.05 to 2.13 Å for the first coordination shell. A second shell contribution arises from interactions with tetrahedralnetwork Si(Na/Al) atoms located at around 3.7 Å [START_REF] Farges | Structural environments of incompatible elements in silicate glass/melt systems: I. Zirconium at trace levels[END_REF]Farges, 1996[START_REF] Louvel | Zr complexation in high pressure fluids and silicate melts and implications for the mobilization of HFSE in subduction zones[END_REF]. Although Zr 4+ solubility in silicate melts is a strong function of temperature, pressure, and melt composition [START_REF] Watson | Zircon saturation in felsic liquids: Experimental results and applications to trace element geochemistry[END_REF][START_REF] Harrison | Kinetics of zircon dissolution and zirconium diffusion in granitic melts of variable water content[END_REF][START_REF] Farges | Structural environments of incompatible elements in silicate glass/melt systems: I. Zirconium at trace levels[END_REF][START_REF] Boehnke | Zircon saturation re-revisited[END_REF][START_REF] Louvel | Zr complexation in high pressure fluids and silicate melts and implications for the mobilization of HFSE in subduction zones[END_REF], experiments at conditions relevant to magmatic environments have found Zr speciation to be independent of temperature and pressure [START_REF] Louvel | Zr complexation in high pressure fluids and silicate melts and implications for the mobilization of HFSE in subduction zones[END_REF], as well as water or halogen (F-and Cl-) contents [START_REF] Farges | Structural environments of incompatible elements in silicate glass/melt systems: I. Zirconium at trace levels[END_REF]Farges, 1996;[START_REF] Louvel | Zr complexation in high pressure fluids and silicate melts and implications for the mobilization of HFSE in subduction zones[END_REF][START_REF] Wang | Coordination of Zr4+/Hf4+/Nb5+/Ta5+ in silicate melts: insight from first principles molecular dynamics simulations[END_REF]. Only at high-P-T in fluids with very low Si-Na contents (i.e., <35 wt.%) has Zr 4+ been observed to occur in 8-fold coordination in a liquid phase [START_REF] Wilke | Zircon solubility and zirconium complexation in H2O+Na2O+SiO2±Al2O3 fluids at high pressure and temperature[END_REF][START_REF] Louvel | Zr complexation in high pressure fluids and silicate melts and implications for the mobilization of HFSE in subduction zones[END_REF]. On the other hand, Zr 4+ coordination in common Zr-bearing minerals, varies. It is 8-fold in the zircon structure (tetragonal ZrSiO4; Robinson et al., 1971), 7-fold in baddeleyite (monoclinic ZrO2;Smith and Newkirk, 1965) and can also replace 6-coordinated Ti in Ti-oxide minerals such as ilmenite (FeTiO3) and rutile (TiO2) [START_REF] Abrahams | Rutile: Normal Probability Plot Analysis and Accurate Measurement of Crystal Structure[END_REF]Farges et al., 1996).

Mean RZr-O in the first coordination shell of Zr 4+ are 2.14-2.28 Å for zircon and 2.04-2.18 Å for baddeleyite [START_REF] Louvel | Zr complexation in high pressure fluids and silicate melts and implications for the mobilization of HFSE in subduction zones[END_REF].

In this work, two mineral structures are considered for Zr-bearing phases: zircon (tetragonal ZrSiO4), in which Zr is found in 8-fold coordination, and baddeleyite (monoclinic ZrO2), where Zr is found in 7-fold coordination (Fig. 1a,b;[START_REF] Robinson | The Structure of Zircon: A Comparison With Garnet[END_REF][START_REF] Smith | Crystal structure of baddeleyite (monoclinic ZrO2) and its relation to the polymorphism of ZrO2[END_REF]. To model Zr in the silicate melt, several mineral structures were computed: gittinsite (Fig. 1c; ZrCaSi2O7), sabinaite (Fig. 1d; Na4Zr2TiC4O16), and vlasovite (Fig. 1e; Na2ZrSi4O11), where Zr is found in, respectively, 6-fold, 8-fold and 6-fold coordination.

Vlasovite was chosen here as the best proxy for understanding Zr bonding properties in silicate melts, because Zr 4+ ions in this structure occupy 6-fold coordinated sites in alkalizirconosilicate clusters just as those observed in silicate liquids (e.g., [START_REF] Louvel | Zr complexation in high pressure fluids and silicate melts and implications for the mobilization of HFSE in subduction zones[END_REF], and Zr-O bond distances in this phase (ca. 2.09 Å) are similar to those found in silicate liquids (i.e., 2.05 to 2.13 Å; [START_REF] Farges | Structural environments of incompatible elements in silicate glass/melt systems: I. Zirconium at trace levels[END_REF][START_REF] Louvel | Zr complexation in high pressure fluids and silicate melts and implications for the mobilization of HFSE in subduction zones[END_REF]. To more precisely assess the domain of variation of Zr fractionation properties, Zr-bearing minerals of contrasting compositions and Zr bonding characteristics and coordination number, yet representative of magmas, were further considered. Only experimentally resolved and ordered structures, with unit cells sufficiently small to avoid prohibitively difficult computation, were retained.

Gittinsite and sabinaite were found to be the best choices for this, among other existing minerals that were considered (e.g., catapleiite Na2ZrSi3O9.2H2O, armstrongite CaZrSi6O17H4, tumchaite Na2(Zr,Sn)Si4O11.2H2O, painite CaZrBAl9O18, elpidite Na2ZrSi6O18H6, georgechaoite NaKZrSi3O9.2H2O, or calzirtite Ca2Zr5Ti2O16).

Methods

The equilibrium fractionation of Zr isotopes between two phases is obtained by combining the β-factors of both phases. The β-factor corresponds to the isotopic fractionation factor between the phase, and a perfect gas of Zr atoms. This quantity can be computed within the harmonic approximation from the vibrational frequencies of the phase of interest [START_REF] Urey | The thermodynamic properties of isotopic substances[END_REF]. If a is a crystalline solid, and Y the element of interest, the β-factor of phase a for the Y/Y * isotopic substitution can be written as:

Eq. 1 where Q(AY * N) is the partition function of the system having all the Y atoms substituted by Y * .

The partition function of a crystalline solid can be estimated from its harmonic phonon frequencies:

Eq. 2 where νq,i are the frequencies of the phonon with wavevector q and branch index i = 1,3N at .

N at is the number of atoms in the unit cell, N is the number of Y atoms in the unit cell, T is the absolute temperature, k is the Boltzman constant and h is Planck's constant. The product is performed on a sufficiently large grid of Nq q-vectors in the Brillouin zone. The * symbol above the product indicates that the three translational modes with ν0,i = 0 are not considered.

The phonon frequencies are computed from first principles within density functional theory (DFT) [START_REF] Hohenberg | Inhomogeneous Electron Gas[END_REF][START_REF] Kohn | Self-Consistent Equations Including Exchange and Correlation Effects[END_REF], using the local density approximation to the exchange-correlation functional of Perdew and Zünger (PZ) [START_REF] Perdew | Self-interaction correction to density-functional approximations for many-electron systems[END_REF]. This approach is necessary because the gradient-corrected functionals typically used in previous studies, such as BLYP or PBE (e.g., [START_REF] Dupuis | Fractionation of silicon isotopes in liquids: The importance of configurational disorder[END_REF], have been shown inappropriate to describe Zr-containing materials such as the ZrO2 polymorphs (Fadda et al., 2010a, b). A plane-wave basis set, and atomic pseudopotentials were used. The pseudopotentials for Si and O were generated for the PZ functional using the same characteristics as those described in the electronic annexes of [START_REF] Méheut | Equilibrium isotopic fractionation in the kaolinite, quartz, water system: Prediction from first-principles density-functional theory[END_REF], which were generated for the PBE functional. The pseudopotentials for Na, Ti, C and Ca are taken from the PSlibrary [START_REF] Dal Corso | Pseudopotentials periodic table: From H to Pu[END_REF]. For Zr, we considered two pseudopotentials. The first one is taken from the PSLibrary and is therefore ultrasoft. This permits a numerically efficient calculation, but the computation of Raman intensities is not implemented for this type of pseudopotentials in the Quantum Espresso software package. In order to compute Raman intensities for baddeleyite and zircon, and compare precisely their calculated gamma-point frequencies with experimental spectra, we generated a Zr norm-conserving pseudopotential [START_REF] Troullier | Efficient pseudopotentials for plane-wave calculations[END_REF] in the Kleinman-Bylander form [START_REF] Kleinman | Efficacious Form for Model Pseudopotentials[END_REF], similar to the ones used for O and Si. For this second Zr pseudopotential, hereafter referred to as "Zrtm", we used the 4p 6 4d 4 4s 2 configuration with core radii of 1.20 a.u. for each orbital and the local potential is s. For the computations of all the minerals not using the Zrtm pseudopotential, electronic wave-functions are expanded in plane-waves up to an energy cutoff εcut = 80 Ry and the charge density cut-off is set to 4 εcut. For the calculations of zircon and baddeleyite using Zrtm, electronic wave-functions are expanded in plane-waves up to an energy cutoff εcut = 200 Ry and the charge density cut-off is set to 4 εcut. Electronic integration is done by sampling the Brillouin zone with a 3×3×3 k-point grid for zircon, baddeleyite and gittinsite, a 3×1×3 grid for sabinaite and a 2×2×2 grid for vlasovite [START_REF] Monkhorst | Special points for Brillouin-zone integrations[END_REF].

Phonon frequencies were computed using linear response theory [START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF], with the Quantum-Espresso package [START_REF] Giannozzi | QUANTUM ESPRESSO: a modular and opensource software project for quantum simulations of materials[END_REF]. Interatomic force-constants are obtained from the dynamical matrices computed exactly (within DFT) on a n×n×n grid of qvectors (n=3 for baddeleyite, zircon, sabinaite, gittinsite, n=2 for vlasovite). Long-range effects are taken into account by computing Born effective-charges and static dielectric constants [START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF]. Dynamical matrices and thus phonon frequencies, can then be obtained in any point of the reciprocal space by Fourier-interpolation of the force constants.

For all materials, the vibrational partition function is converged with a 6×6×6 interpolation grid. In order to help correlate calculated frequencies at gamma-point with experimental spectroscopic (Raman or infrared (IR)) data, we computed the intensities (IR or Raman crosssections). IR cross-sections were computed for all the modes whereas raman cross-sections can only be computed for the Zrtm calculation. For Raman-active modes, non-resonant Raman cross-sections were computed using second-order response as in [START_REF] Lazzeri | First-Principles Calculation of Vibrational Raman Spectra in Large Systems: Signature of Small Rings in Crystalline ${\mathrm{S}\mathrm{i}\mathrm{O}}_{2}[END_REF]. For IR-active modes, cross-sections were computed from the Born effective charges and static dielectric constant (see [START_REF] Balan | First-principles modeling of the infrared spectrum of kaolinite[END_REF] for a more detailed description). Note that we used these intensities in a rather qualitative way, in order to distinguish between "important" and weak modes. A proper modeling of IR or Raman spectra would require treatment of polarization, medium and particle shape effects (for IR: [START_REF] Balan | First-principles modeling of the infrared spectrum of kaolinite[END_REF]for Raman: Lazzeri and Mauri, 2003), which is beyond the scope of the present study.

Results

Structure and vibrational properties

The equilibrium computed structures of the new materials considered in this study are reported in Table 1 and compared with experimental data. In all cases, atomic positions and cell parameters are obtained after relaxation at zero pressure, until the residual forces are less than 10 -3 Ry/Å (10 -12 N). Except for the c-axis of sabinaite, all cell parameters of the relaxed structures underestimate experimental values by 0.5-1.5%. These calculated structural parameters are further compared with selected previous theoretical studies, representative of the range of values obtained for calculations made with the same functional (Local Density Approximation (LDA), Table 1). In both cases, our calculations fall within the range of previous results. For baddeleyite, we refer to the detailed discussion of [START_REF] Fadda | First-principles study of the structural and elastic properties of zirconia[END_REF] for comparison with other studies.

In the case of sabinaite, the c-axis length is underestimated by 2.3%. Sabinaite structure consists of (Zr,Ti)O layers alternating with NaCO3 layers. While Zr-O distances are normally underestimated by the calculation (within 0.5 -1.5%), we note that Na-Na distances are more notably underestimated (by 3.9%). Similarly, calculation of calcite structure within LDA results in c-axis lengths (and therefore Ca-Ca distances) are underestimated by 3.7% (e.g., [START_REF] Valenzano | Ab Initio Study of the Vibrational Spectrum and Related Properties of Crystalline Compounds; the Case of CaCO3 Calcite[END_REF]. The unusual discrepancy with respect to the experimental data available for sabinaite is therefore interpreted as a failure of LDA to treat cation-CO3 2-bonds consistently with other bonds. Although this observation implies that the sabinaite calculations should be taken with more caution than those from other phases, the obtained result is useful because it informs on the effect that a change in the chemical environment of Zr might have on its fractionation properties.

The computed phonon frequencies of baddeleyite and zircon at the zone center (wave vector = 0) are compared with spectroscopic measurements in Table 2 (Raman data) and Table 3 (IR data), and illustrated in Figure 2. All calculated mode frequencies, along with IR and Raman intensities, are reported and most are assigned to a corresponding mode.

Theoretical frequencies show very good agreement with measured ones.

Fractionation properties

The β-factors for zircon, baddeleyite and the mineral structures representative of Zr in silicate melts are shown in Figure 3. Predicted equilibrium isotope fractionations between all calculated phases relative to zircon (i.e., calculated from differences in β-factors as ΔA-B ≃ 10 3 lnαA-B = 10 3 lnβA-10 3 lnβB) are shown in Figure 4. The coefficients of the polynomial expansion of 10 3 lnβ in terms of even powers of the inverse of temperature (in K): 10 3 lnβ = A1/T 2 + A2/T 4 + A3/T 6 , are reported in Table 4, from which the mean force constant, noted <F> (in N/m), of the chemical bond between Zr and its coordination atoms in each phase, are calculated (see Table 4 footnote for details). For a given phase, at magmatic temperatures, the equilibrium isotope fractionation factor can be easily calculated using the force constant (e.g., [START_REF] Dauphas | A general moment NRIXS approach to the determination of equilibrium Fe isotopic fractionation factors: application to goethite and jarosite[END_REF]Dauphas et al., 2018 and references therein), as:

1000 lnβ = <F> 1000 (1/mlight -1/mheavy) ħ 2 / (8 kB 2 T 2 ),
where mlight and mheavy are the masses of the light and heavy isotopes (in kg), respectively, ħ is the reduced Planck's constant (J⋅s), kB is Boltzmann's constant (J/K), and T is the absolute temperature. For 94 Zr/ 90 Zr ratios, this expression reduces to 1000 lnβ( 94/90 Zr ≈ 2081.6 <F> / T 2 .

At thermodynamic equilibrium, heavy isotopes typically populate shorter and stiffer bonds with the lowest coordination number [START_REF] Schauble | Applying Stable Isotope Fractionation Theory to New Systems[END_REF][START_REF] Young | High-temperature equilibrium isotope fractionation of non-traditional stable isotopes: Experiments, theory, and applications[END_REF][START_REF] Blanchard | Equilibrium Fractionation of Non-traditional Isotopes: a Molecular Modeling Perspective[END_REF]. Consistent with these first-order expectations, vlasovite (VI) and gittinsite (VI) are isotopically indistinguishable from each other (<F> = 407.8 N/m for both phases), and heavier than zircon (VIII) (<F> = 392.7 N/m). On the other hand, baddeleyite (VII) (<F> = 372.4 N/m) and sabinaite (VIII) (<F> = 380.7 N/m) are lighter than zircon (VIII).

In order to better understand why materials with higher coordination can be isotopically heavier (e.g., zircon or sabinaite w/r to baddeleyite), we analyzed the vibrational contributions to fractionation properties for these various materials (methodology described in [START_REF] Méheut | Structural control over equilibrium silicon and oxygen isotopic fractionation: A first-principles density-functional theory study[END_REF]. The results, shown in Figure 5, reveal a complex situation. For zircon and baddeleyite (Fig 5a), the main contribution to Zr isotopic fractionation comes from low (typically 100-400 cm -1 ) frequencies. However, zircon also presents some significant contributions at frequencies higher than baddeleyite. In particular, the contribution of the highest frequencies of zircon (~900 cm -1 ) is non-negligible (0.08‰ in total). These high frequencies correspond to Si-O stretching modes. In addition to coordination, this observation suggests the importance of phase chemical composition in influencing Zr isotope fractionation properties, and more specifically of the presence of Si in the structure. The presence of Si in the mineral leads to high frequencies modes (Si-O stretching modes), and, assuming that these would present a non-negligible coupling with Zr displacement, high frequency Zr displacement would occur, leading to high frequency contribution to the βfactor. The two other materials containing Si (gittinsite and vlasovite, Fig. 5c) present a very similar contribution of SiO stretching modes (~0.08‰ in total). Now, looking at sabinaite (Fig. 5b), which is a carbonate, we can also distinguish a contribution of C-O stretching frequencies (~1350-1650 cm -1 ). This contribution is smaller than for Si-O (~ 0.04 ‰), which means that the coupling of Zr movement with these displacements might be smaller.

Nevertheless, this feature is once again able to explain the heavier isotopic composition of sabinaite with respect to baddeleyite (despite of their respective coordination), as the integrated contribution to the β-factor of sabinaite becomes greater than that of baddeleyite above this spectral domain.

Perhaps more important than the relative isotope fractionation factors is the fact that, at magmatic temperatures, the magnitude of these isotope effects are extremely small. Even at 800°C, the maximum equilibrium zircon-melt or baddeleyite-melt fractionation is ≤ 0.05 ‰.

At 1200°C these effects are ≤ 0.02-0.03‰.

Discussion

Negligible equilibrium Zr isotope fractionation vs. observed magmatic trends

The results of our ab initio calculations indicate that vibrational equilibrium Zr isotope effects can only yield minute zircon-melt Zr isotope fractionation; i.e., Δ 94/90 Zrzircon-melt = (δ 94/90 Zrzircon-δ 94/90 Zrmelt) ~0.05 ‰ at 800°C. This very limited degree of fractionation is consistent with the inter-and intra-crystal Zr isotope homogeneity observed in some reference zircon localities: e.g., 91500, Mud Tank and Plešovice (LA-ICPMS: [START_REF] Zhang | Determination of Zr isotopic ratios in zircons using laser-ablation multiple-collector inductively coupled-plasma mass-spectrometry[END_REF]SIMS: Kirkpatrick et al., 2019;double-spiked solutions: Tompkins et al., 2020).

These small fractionation factors are, however, at least one order of magnitude smaller than those required to explain the data obtained for multiple igneous systems (Fig. 6). Ibañez-Mejia and Tissot (2019) inferred a Δ 94/90 Zrzircon-melt ≈ +1 ‰ through direct comparison of δ 94/90 Zr values from 70 individual crystals of Zr-rich accessory phases (zircon and baddeleyite) and the bulk rock of a well-characterized gabbroic igneous cumulate (FC-1, Duluth gabbro). [START_REF] Inglis | Isotopic fractionation of zirconium during magmatic differentiation and the stable isotope composition of the silicate Earth[END_REF] hypothesized that the positive correlation between δ 94/90 Zr and SiO2 wt. % in bulk volcanic-rock samples from the Hekla Volcano (Iceland) was the result of zircon removal with a Δ 94/90 Zrzircon-melt ≈ -0.50 ‰. The variable (~-0.12 to -0.45 ‰) fractionation factors reported by [START_REF] Guo | Significant Zr isotope variations in single zircon grains recording magma evolution history[END_REF] for zircons in the juvenile Gangdese Arc (southern Tibet), come from laser ablation determination of zircon crystal zoning, which revealed lighter compositions towards the core and heavier compositions towards the rims of the crystals. Interpretations regarding the driving mechanisms behind these effects also varied between studies: [START_REF] Inglis | Isotopic fractionation of zirconium during magmatic differentiation and the stable isotope composition of the silicate Earth[END_REF] and [START_REF] Guo | Significant Zr isotope variations in single zircon grains recording magma evolution history[END_REF] attributed these large effects to mass-dependent equilibrium isotope fractionation during Zr crystallization, whereas Ibañez-Mejia and Tissot (2019) highlighted the fact that equilibrium effects may not be the sole driver of these isotope fractionations. The results of the ab initio calculations presented here confirm the latter hypothesis, indicating that none of these data can be explained by massdependent, vibrational, equilibrium isotope effects. In the next sections, we consider alternative drivers of Zr isotope fractionations.

Nuclear Field Shifts as drivers of Zr isotope fractionation?

The first-principle calculations presented above quantify mass-dependent isotope fractionations stemming from differences in the vibrational properties of the chemical bonds in Zr isotopologues [START_REF] Bigeleisen | Calculation of Equilibrium Constants for Isotopic Exchange Reactions[END_REF][START_REF] Urey | The thermodynamic properties of isotopic substances[END_REF], but do not account for the so-called Nuclear Field Shift (NFS) effects, which arise from differences in the volume and shape of atomic nuclei [START_REF] Bigeleisen | Nuclear Size and Shape Effects in Chemical Reactions. Isotope Chemistry of the Heavy Elements[END_REF][START_REF] Knyazev | Specific Effects of Heavy Nuclei in Chemical Equilibrium[END_REF][START_REF] Schauble | Role of nuclear volume in driving equilibrium stable isotope fractionation of mercury, thallium, and other very heavy elements[END_REF][START_REF] Yang | Nuclear field shift effects on stable isotope fractionation: a review[END_REF]. Below, we examine the possibility that NFS could result in significant Zr isotope fractionation at magmatic temperatures.

NFS effects manifest as, and are named after, shifts in the emission/absorption frequencies of isotopic atoms that cannot be explained by mass differences and hyperfine splitting [START_REF] King | Isotope Shifts in Atomic Spectra[END_REF]. To first order, they scale proportionally to | (0)| 2 δ <r 2 Nuc>, where | (0)| 2 is the total electron density at the nucleus, and δ<r 2 Nuc> is the isotopic change in the mean-square nuclear charge radius (e.g., [START_REF] King | Isotope Shifts in Atomic Spectra[END_REF][START_REF] Schauble | Role of nuclear volume in driving equilibrium stable isotope fractionation of mercury, thallium, and other very heavy elements[END_REF]. Because nuclei with odd (or magic) numbers of neutrons, tend to have smaller nuclear charge radii than expected based on simple interpolation between their neighboring isotopes, NFS effects tend to display diagnostic, mostly odd-even staggering, patterns (e.g., Hg isotopes, [START_REF] Bergquist | Mass-dependent and -independent fractionation of hg isotopes by photoreduction in aquatic systems[END_REF]. As mass-dependent fractionation decreases with increasing mass, and NFS increase with electron density at the nucleus, NFS effects are most pronounced in heavy elements (i.e., with large nuclei), where they can become the dominant fractionation mechanism (e.g., NFS during U redox reactions are 3x larger than mass-dependent effects, [START_REF] Bigeleisen | Nuclear Size and Shape Effects in Chemical Reactions. Isotope Chemistry of the Heavy Elements[END_REF][START_REF] Schauble | Role of nuclear volume in driving equilibrium stable isotope fractionation of mercury, thallium, and other very heavy elements[END_REF]Abe et al., 2008). Furthermore, while mass-dependent fractionation scales proportionally to 1/T 2 , NFS effects scale with 1/T, which implies that their relative contribution to the total isotope fractionation increases in high-T (e.g., igneous) environments.

To determine the magnitude of equilibrium NFS effects, one can either (i) perform ab initio calculations that directly model the electronic structures of the nuclei while explicitly varying nuclear charge radii (e.g., [START_REF] Knyazev | Specific Effects of Heavy Nuclei in Chemical Equilibrium[END_REF][START_REF] Schauble | Role of nuclear volume in driving equilibrium stable isotope fractionation of mercury, thallium, and other very heavy elements[END_REF]Abe et al., 2008), or (ii) perform equilibrium exchange reactions and measure the resulting isotopic effects, which can be broken down into mass-dependent and NFS contributions (e.g., [START_REF] Fujii | Nuclear field shift effect in the isotope exchange reaction of cadmium using a crown ether[END_REF], provided both odd and even isotopes are analyzed. At this time, neither type of data are available for Zr. Nevertheless, some general considerations can be made. Because electrons in s-orbitals can cause larger electron density at the nucleus than p-, d-, or f-orbitals, heavier isotopes will preferentially be enriched in the species with fewer s electrons, or more p, d and f electrons [START_REF] Schauble | Role of nuclear volume in driving equilibrium stable isotope fractionation of mercury, thallium, and other very heavy elements[END_REF]. This, for instance, results in 238 U enrichments in reduced, U 4+ -[Rn]5f 2 , species relative to oxidized, U 6+ -[Rn]5f 0 , species, which has led to the development of U isotopes as tracers of redox conditions in the ocean (e.g., [START_REF] Weyer | Natural fractionation of 238U/235U. Geochim[END_REF][START_REF] Tissot | Uranium isotopic compositions of the crust and ocean: Age corrections, U budget and global extent of modern anoxia[END_REF], and possibly in magmatic environments (e.g., [START_REF] Tissot | 238U/235U measurement in single-zircon crystals: implications for the Hadean environment, magmatic differentiation and geochronology[END_REF]. In all phases considered here, zirconium is present as Zr 4+ , removing the possibility that NFS effects arise from redox changes and associated differences in the population of electron orbitals. Instead, differences in the electronic structure of the nuclei must therefore mainly be driven by differences in size and shape of the Zr isotope nuclei, and thus changes in the mean-squared nuclear charge radius, δ<r 2 >. Figure 7a shows the δ<r 2 > of natural Zr isotopes, and Fig 7b shows δ<r 2 > variations normalized to 90 Zr and 94 Zr. The isotopic patterns for Zr shown in Fig. 7 do not conform to the typical odd-even staggering observed for many other elements (e.g., Mo) for two reasons: (i) 90 Zr has a magic number of neutrons (N= 50), which means it consists of entirely closed shells of neutrons [START_REF] King | Isotope Shifts in Atomic Spectra[END_REF], resulting in a particularly compact nucleus which removes most of the expected 90 Zr-91 Zr-92 Zr odd-even staggering. (ii) 96 Zr has the peculiar property of being doubly-magic [START_REF] Angeli | N and Z dependence of nuclear charge radii[END_REF], with filled subshells of both protons (Z=40) and neutrons (N=56), resulting in a more compact nucleus, making 96 Zr deviate from the expected mass-dependent relationship of even isotopes. Because NFS effects scale proportionally to δ<r 2 > values, Zr NFS-induced effects, normalized to 90 Zr-94 Zr, are expected to show significant and opposite anomalies on, respectively, 96 Zr and 92 Zr (Fig. 7b).

Currently, detection of NFS effects in natural samples is limited to the heaviest elements (e.g., Hg, Tl, U), and calculations of NFS effects for elements in the ~90-100 amu mass range predict only small effects [START_REF] Knyazev | Specific Effects of Heavy Nuclei in Chemical Equilibrium[END_REF][START_REF] Schauble | Role of nuclear volume in driving equilibrium stable isotope fractionation of mercury, thallium, and other very heavy elements[END_REF]. With no available data on the magnitude of NFS effects for Zr isotopes, we use molybdenum (Mo) as a proxy to draw first-order estimates on Zr effects. Molybdenum and Zr are both transition metals, with three isobaric isotopes (92-94-96), and the variations in the root mean-square nuclear charge radii of their isotopes are similar, if only slightly larger for Mo (Fig. 7b).

Liquid-liquid isotope exchange experiments performed at 295 K (Fujii et al., 2006b) reported isotope fractionation patterns consistent with NFS equilibrium effects of ~ 0.08 ‰/amu for Mo. Assuming a simple 1/T temperature dependency, this translates to effects of ~ 0.02 ‰/amu at magmatic temperatures of 800-1000°C. Based on the similar extent of variability of the mean-square nuclear charge radii of Mo and Zr, it is reasonable to expect similar magnitudes of NFS effects for Zr, which would result in a δ 94/90 Zr of ~0.08‰ at magmatic temperatures. While 2 to 3x larger than the -vibrational-equilibrium mass-dependent effects we predict (Fig. 4), this value is still nowhere near the ~+1‰ fractionation factor needed to explain the isotope effects in the individual zircons and baddeleyite of the FC-1 gabbro (Fig. 6, Ibañez-Mejia and Tissot 2019).

An upper limit can be independently placed on the magnitude of NFS effects potentially affecting the published zircon and baddeleyite data from Ibañez-Mejia and Tissot (2019) and [START_REF] Tompkins | Zirconium stable isotope analysis of zircon by MC-ICP-MS: methods and application to evaluating intra-crystalline zonation in a zircon megacryst[END_REF]. In both studies, a double-spike (DS) method was used to quantify the mass-dependent Zr isotope fractionation of the samples (δ 94/90 Zr). The presence of mass-independent anomalies, if unaccounted for during DS data reduction, would result in inaccurate estimates for the extent of mass-dependent fractionation [START_REF] Eugster | Isotopic analyses of barium in meteorites and in terrestrial samples[END_REF][START_REF] Niederer | Absolute isotopic abundances of Ti in meteorites[END_REF][START_REF] Burkhardt | Evidence for Mo isotope fractionation in the solar nebula and during planetary differentiation[END_REF][START_REF] Greber | Titanium stable isotopic variations in chondrites, achondrites and lunar rocks[END_REF][START_REF] Hu | Double-spike data reduction in the presence of isotopic anomalies[END_REF].

Typically for double-spike data reduction, only 4 isotopes (3 isotope ratios) are considered, providing an exact solution to the three unknowns being solved for (i.e., the proportion of spike, instrumental mass-bias and natural mass bias between the sample and the standard).

For such scenarii, [START_REF] Hu | Double-spike data reduction in the presence of isotopic anomalies[END_REF] derived analytical equations that rigorously estimate the corrections on delta values required by the presence of isotopic anomalies. In the Ibañez-Mejia and Tissot (2019) and [START_REF] Tompkins | Zirconium stable isotope analysis of zircon by MC-ICP-MS: methods and application to evaluating intra-crystalline zonation in a zircon megacryst[END_REF] work, however, all 5 isotopes of Zr are used in a minimization scheme to perform the DS data reduction, implying that an exact solution not always necessarily exists and that the impact of mass-independent effects on delta values can be assessed numerically. As detailed in [START_REF] Tompkins | Zirconium stable isotope analysis of zircon by MC-ICP-MS: methods and application to evaluating intra-crystalline zonation in a zircon megacryst[END_REF], a 4-ratio minimization scheme provides a way to check for the presence of mass-independent effects (relative to the normalizing standard), which would manifest as isotopic anomalies even after data reduction. For all samples measured in these two studies, all isotope ratios were found to be mass-dependent within uncertainty (~±0.01 ‰/amu), thus allowing calculation of a maximum permissible NFS effect to these data, as detailed below.

Assuming that the pattern of NFS isotopic effects scales as the 90 Zr-94 Zr normalized δ<r 2 > shown in Fig. 7b, isotopic anomalies of increasing magnitude were (numerically) added to the isotopic composition of the Zr standard to simulate the presence of NFS effects, and these anomalous standards were mixed in ideal proportion with the DS. For each mixture in this hypothetical dataset, the DS data reduction was performed using the minimization approach described in [START_REF] Tompkins | Zirconium stable isotope analysis of zircon by MC-ICP-MS: methods and application to evaluating intra-crystalline zonation in a zircon megacryst[END_REF], and the difference between the algorithm solution and standard composition was calculated. The results of these calculations are plotted in Figures 8 and9. The observed mass-dependency of all Zr isotope ratios within ±0.01‰/amu in the Ibañez-Mejia and Tissot (2019) and [START_REF] Tompkins | Zirconium stable isotope analysis of zircon by MC-ICP-MS: methods and application to evaluating intra-crystalline zonation in a zircon megacryst[END_REF] datasets, places tight constraints on the apparent mass-dependent fractionation that can be attributed to NFS effects. Considering the calculations shown in Fig. 8, at most 0.001-0.002 ‰/amu effects are permissible, above which the 92 Zr/ 90 Zr ratio would no longer be mass-dependently related to the other Zr isotope ratios. This indicates that (i) NFS effects were not a significant driver of isotopic fractionation in these samples, which display much larger mass-dependent effects (from 0.014 to 0.250‰/amu), and (ii) that NFS-induced mass-independent effects in the samples, expressed in epsilon notation as 91 Zr, 92 Zr and 96 Zr anomalies ( 9x Zr= [( 9x Zr/ /90 ZrStandard+NFS/ 9x Zr/ /90 ZrStandard)-1]*10 4 ), were smaller than, respectively, ~-4.6, +0.4 and -5.2 -units (Figs. 8 and9).

Kinetic effects as drivers of Zr isotope fractionation

In addition to the thermodynamic (equilibrium) effects described above, stable isotope fractionation amongst phases in high-T environments can also be driven by the massdependency of reaction rate constants or self-diffusion coefficients in solids and liquids (i.e., kinetic isotope effects; see reviews by Richter et al., 2009a;[START_REF] Watkins | Kinetic Fractionation of Non-Traditional Stable Isotopes by Diffusion and Crystal Growth Reactions[END_REF]. Mineral growth from a crystallizing magma is a first-order process known to result in compositional (e.g., [START_REF] Albarede | Kinetic disequilibrium in trace element partitioning between phenocrysts and host lava[END_REF] and isotopic (e.g., [START_REF] Watson | Non-equilibrium isotopic and elemental fractionation during diffusion-controlled crystal growth under static and dynamic conditions[END_REF] disequilibrium effects for major and trace elements. These fractionations arise from a kinetic imbalance between crystal growth, chemical diffusion in the melt, and the isotopic mass dependency of diffusion. As shown by [START_REF] Albarede | Kinetic disequilibrium in trace element partitioning between phenocrysts and host lava[END_REF], for elements with equilibrium fractionation coefficients (Keq) between solid and melt smaller than unity, a 'pileup' of the species of interest occurs near the moving solid-liquid interface when diffusion is too slow to maintain a uniform concentration throughout the liquid. Conversely, for elements with Keq > 1 in the crystallizing solid, growth of the mineral will deplete that species in the surrounding melt, inducing a 'draw down' near the moving solid-liquid interface. The concentration gradients that develop as the solid grows lead to diffusive fluxes within the liquid [START_REF] Tiller | The redistribution of solute atoms during the solidification of metals[END_REF][START_REF] Smith | A mathematical analysis of solute redistribution during solidification[END_REF][START_REF] Albarede | Kinetic disequilibrium in trace element partitioning between phenocrysts and host lava[END_REF], and such regions are known as diffusive boundary layers (DBLs) [START_REF] Watson | Non-equilibrium isotopic and elemental fractionation during diffusion-controlled crystal growth under static and dynamic conditions[END_REF][START_REF] Watkins | Kinetic Fractionation of Non-Traditional Stable Isotopes by Diffusion and Crystal Growth Reactions[END_REF].

Zirconium is a moderately to highly incompatible lithophile element with Keq values smaller than unity for most rock-forming phases growing from a silicate melt (see data compilation of [START_REF] Keller | Temporal variation in relative zircon abundance throughout Earth history[END_REF]. Thus, as a magma cools and fractionally crystallizes, this incompatibility results in a progressive Zr enrichment in differentiated liquids. For instance, in the case of the Duluth Complex FC-1 anorthositic gabbro, magma crystallization is dominated by the initial growth of plagioclase, followed by intra-cumulus growth of clinopyroxene and Fe-Ti oxides within trapped melt pockets (Ibañez-Mejia and Tissot, 2019).

Zr is extremely incompatible in plagioclase (Keq < 0.01; e.g., [START_REF] Bindeman | Trace element partitioning between plagioclase and melt: investigation of dopant influence on partition behavior[END_REF][START_REF] Nielsen | Trace element partitioning between plagioclase and melt: An investigation of the impact of experimental and analytical procedures[END_REF], strongly to moderately incompatible in clinopyroxene (Keq = 0.12 -0.27; e.g., [START_REF] Watson | Partitioning of zirconium between clinopyroxene and magmatic liquids of intermediate composition[END_REF][START_REF] Hart | Experimental cpx/melt partitioning of 24 trace elements[END_REF][START_REF] Huang | Effect of melt structure on traceelement partitioning between clinopyroxene and silicic, alkaline, aluminous melts[END_REF], and moderately incompatible in ilmenite (Keq= 0.33 -0.80; e.g., [START_REF] Nakamura | Hf, Zr, and REE partition coefficients between ilmenite and liquid: Implications for lunar petrogenesis[END_REF][START_REF] Klemme | The partitioning of trace elements between ilmenite, ulvospinel, armalcolite and silicate melts with implications for the early differentiation of the moon[END_REF]. Therefore, growth of all these phases from a magma will lead to an overall Zr enrichment in the residual liquid, as well as the development of markedly enriched DBLs around the crystallizing phases. To illustrate this, Figure 10 shows an example calculation of the [Zr] profiles that would develop in a solid (left) and co-existing liquid (right) during clinopyroxene crystal growth at isothermal conditions, as a function of variable boundary migration velocities (v) expressed in cm/Myr. These calculations clearly illustrate that the magnitude of the concentration gradient that develops, and therefore the resulting diffusive flux, are proportional to v (or the v/D ratio). Furthermore, DBLs such as those shown in Figure 10 will be the most likely regions of zircon nucleation during magmatic crystallization, owing to their much higher [Zr] and the concentration dependence of zircon saturation in silicate liquids [START_REF] Boehnke | Zircon saturation re-revisited[END_REF]. Kinetic theory predicts that the ratio of the diffusivities of two isotopes in an ideal gas is proportional to the square root of the ratio of their masses (i.e., D1/D2 = [ 1/ 2]), but for condensed systems this formulation overstates the mass-dependency of diffusion due to non-negligible interactions between the diffusing species and their nearest neighbors [START_REF] Richter | Isotope fractionation by diffusion in molten oxides[END_REF]. Thus, for aqueous solutions and silicate melts it has been customary to express the ratio of the self-diffusivities using an inverse power-law function:

2 1 = 1 2 Eq. 3
where Dx are the diffusivities of the two isotopes in question, mx are their respective masses, and is a dimensionless empirical parameter that expresses the efficiency of diffusive isotope fractionation [START_REF] Richter | Isotope fractionation by diffusion in molten oxides[END_REF]. This last term should not be confused with the -factors presented in Figs. 3 and5, which characterize the equilibrium isotopic fractionation between two phases. Although no experimental data yet exists for constraining (in Eq. 3) in silicate melts for elements heavier than Fe [START_REF] Richter | Isotope fractionation by chemical diffusion between molten basalt and rhyolite[END_REF][START_REF] Richter | Magnesium isotope fractionation in silicate melts by chemical and thermal diffusion[END_REF](Richter et al., , 2009b;;[START_REF] Watkins | Diffusion of multiisotopic chemical species in molten silicates[END_REF], or for trace elements other than Li [START_REF] Richter | Isotope fractionation by chemical diffusion between molten basalt and rhyolite[END_REF][START_REF] Richter | Isotope fractionation of Li and K in silicate liquids by Soret diffusion[END_REF], some useful approximations can be made from the existing data: (i) the overall range of factors for various elements in silicate melts is between 0 and 0.22 (see recent review in [START_REF] Watkins | Kinetic Fractionation of Non-Traditional Stable Isotopes by Diffusion and Crystal Growth Reactions[END_REF]; and (ii) the efficiency of diffusive fractionation scales proportionally to the diffusivity of the species of interest normalized to the diffusivity of the solvent, or Si in the case of silicate liquids [START_REF] Watkins | Influence of liquid structure on diffusive isotope separation in molten silicates and aqueous solutions[END_REF]. [START_REF] Baker | Chemical interdiffusion of dacite and rhyolite: anhydrous measurements at 1 atm and 10 kbar, application of transition state theory, and diffusion in zoned magma chambers[END_REF] performed chemical interdiffusion experiments between dacitic and rhyolitic melts at 1 and 10 kbar, and simultaneously determined the diffusivities of Zr and Si at varying SiO2 concentrations. He found the Zr diffusivity to be invariably faster than that of Si, with average DZr/DSi of 1.43, 1.62 and 2.01 for SiO2 contents of 65, 70 and 75 wt.%, respectively. When compared with the solvent-normalized diffusivities and coefficients experimentally determined for Ca, Mg and Fe, these experimental DZr/DSi ratios fall within a range where the coefficients for major elements vary between 0.03 and 0.1 (Fig. 11). As low abundance elements exhibit, in general, higher mass discrimination (i.e., higher coefficients) than major elements [START_REF] Richter | Isotope fractionation by chemical diffusion between molten basalt and rhyolite[END_REF][START_REF] Richter | Magnesium isotope fractionation in silicate melts by chemical and thermal diffusion[END_REF](Richter et al., , 2009b;;[START_REF] Watkins | Liquid composition-dependence of calcium isotope fractionation during diffusion in molten silicates[END_REF][START_REF] Watkins | Influence of liquid structure on diffusive isotope separation in molten silicates and aqueous solutions[END_REF], the experimental DZr/DSi values shown in Fig. 11 suggest that the values for Zr in silicate melts may in fact be non-negligible.

The FC-1 (Duluth Gabbro) case: can diffusion explain its extreme δ 94/90 Zr variability?

Given the faster diffusion of lighter isotopes, the intuitive expectation for a phase growing from a melt under simple diffusion-controlled conditions (e.g., [START_REF] Bindeman | Zircon Survival, Rebirth and Recycling during Crustal Melting, Magma Crystallization, and Mixing Based on Numerical Modelling[END_REF] is that its bulk composition should be lighter than that of the initial liquid from which it precipitated. In stark contrast, the zircon and baddeleyite crystals from the Duluth complex (Ibañez-Mejia and Tissot, 2019) display extreme δ 94/90 Zr variability, and are mainly heavier, compared to their host rock (assumed to represent the starting magma composition).

To explore the potential role of diffusive fractionation as a driver of large massdependent Zr isotope variations, such as those observed in the FC-1 gabbro (Ibañez-Mejia and Tissot, 2019), we performed a series of exploratory calculations using parameters relevant to FC-1 crystallization. The diffusivity of Zr at 900˚C, which is the zircon crystallization maxima for this sample as determined from Ti-in-zircon thermometry [START_REF] Hiess | Ti-in-zircon thermometry applied to contrasting Archean metamorphic and igneous systems[END_REF] and thermodynamic modelling (Ibañez-Mejia and Tissot, 2019), was calculated using the 'dry' experimental parameters of [START_REF] Mungall | Chemical diffusivities of 18 trace elements in granitoid melts[END_REF]. To constrain a reasonable range of moving boundary velocities, we observe that the maximum timescale of zircon crystallization determined from high-precision U-Pb dating of FC-1 is ca. 220 kyr and that intra-cumulus phase assemblages (i.e., clinopyroxene, ilmenite) in this rock reach up to ca. 0.25 cm in diameter [START_REF] Ibañez-Mejia | Extreme Zr stable isotope fractionation during magmatic fractional crystallization[END_REF]. Thus, a maximum boundary migration velocity during the interval recorded by zircon crystallization should be on the order of ~1 cm/Myr for the largest melt pockets, or possibly slower for smaller pockets. To account for this variability, calculations were performed for v values from 1.0 down to 0.2 cm/Myr. Similarly, and based on the observations summarized in Fig. 11, Zr values ranging from 0.01 to 0.05 were used to account for the uncertainty of this entirely unconstrained parameter. Two cases were considered: a moving plagioclase-melt boundary (Keq = 0.004; [START_REF] Bindeman | Trace element partitioning between plagioclase and melt: investigation of dopant influence on partition behavior[END_REF] and a clinopyroxene-melt boundary (Keq = 0.18; [START_REF] Huang | Effect of melt structure on traceelement partitioning between clinopyroxene and silicic, alkaline, aluminous melts[END_REF].

The results of these calculations are summarized in Fig. 12, which show the isotopic profiles that would develop in a melt as a function of the identity of the growing solid, and the variable values of v and Zr. Intuitively, one might expect that, as the DBL becomes enriched in Zr (e.g., Fig. 10), faster diffusion of lighter isotopes away from the boundary would render the DBL isotopically heavy towards the mineral-melt interface, and isotopically light away from the interface (i.e., towards the tail of the DBL). The results shown in Fig. 12 confirm that this is indeed the case, and that reasonable combinations of v, Zr, and Keq have the potential to produce per-mil level δ 94/90 Zr variations. For similar v/D values, plagioclase-melt boundaries (i.e., lower Keq interfaces) would, in general, result in larger magnitudes of stable isotope fractionation within the DBL relative to phases with higher Keq, such as clinopyroxene. Likewise, faster boundary migration velocities (higher v/D) would result in larger magnitudes of isotope fractionation, as would also be the case for larger magnitudes of Zr (i.e., lower D H /D L ). Because zircon saturation in the melt is directly dependent on [Zr] [START_REF] Boehnke | Zircon saturation re-revisited[END_REF], this mineral would be more likely to nucleate within the DBLs rather than away from them, as 'background ' [Zr] in the liquid would remain systematically lower. This observation implies that, throughout the fractional crystallization history of FC-1, zircon would be more likely to grow in regions of the melt that are enriched in heavy Zr isotopes rather than in isotopically 'normal' regions, as the latter would be accompanied by significantly lower [Zr]. Therefore, even if the equilibrium fractionation factor Δ 94/90 Zrzirconmelt is essentially negligible (Sections 4.1 and 4.2), kinetic isotopic fractionation effects within the DBLs that develop during fractional crystallization can explain the preferential formation of isotopically heavy zircon, and thus the distribution of δ 94/90 Zr values in FC-1 zircon and baddeleyite (Ibañez-Mejia and Tissot, 2019).

The case of negligible (Mud Tank) and negative (Gangdese) fractionation values

Kinetic controls could also explain the δ 94/90 Zr results obtained in zircon from the Mud Tank carbonatite [START_REF] Tompkins | Zirconium stable isotope analysis of zircon by MC-ICP-MS: methods and application to evaluating intra-crystalline zonation in a zircon megacryst[END_REF] and the Gangdese arc, Tibet [START_REF] Guo | Significant Zr isotope variations in single zircon grains recording magma evolution history[END_REF], which suggest negligible to negative Δ 94/90 Zrzircon-melt values, respectively (Fig. 6). Unlike FC-1 zircon, which likely grew within isotopically heavy diffusive boundary layers developed during extensive magmatic fractional crystallization (section 5.3.1.), zircon crystals from the Gangdese arc, Tibet [START_REF] Guo | Significant Zr isotope variations in single zircon grains recording magma evolution history[END_REF], exhibit apparent light-to-heavy δ 94/90 Zr profiles from core to rim. Such profiles would be generally consistent with simple diffusion of Zr from the melt towards the crystallizing solid (e.g., [START_REF] Bindeman | Zircon Survival, Rebirth and Recycling during Crustal Melting, Magma Crystallization, and Mixing Based on Numerical Modelling[END_REF], which would result in lighter isotopes to be concentrated near the core of a growing crystal and progressively heavy compositions towards the rim (e.g., [START_REF] Watson | Non-equilibrium isotopic and elemental fractionation during diffusion-controlled crystal growth under static and dynamic conditions[END_REF]. [START_REF] Guo | Significant Zr isotope variations in single zircon grains recording magma evolution history[END_REF] discarded the possibility that their measured profiles could result from a kinetic control to Zr isotope fractionation arguing that: (i) Zr diffusivities in magmas (D) and typical radial zircon growth rates (R) usually yield R/D <0.01 cm -1 , which would be expected to produce only minor diffusive isotope fractionation effects [START_REF] Watson | Non-equilibrium isotopic and elemental fractionation during diffusion-controlled crystal growth under static and dynamic conditions[END_REF], and (ii) the absence of 'W' shaped profiles from their zircon transects indicates the lack of a kinetic control. There are multiple issues with these interpretations, but most importantly: (i) radial growth rates for zircon, which have been suggested to be on the order of 10 -13 -10 -17 cm/s [START_REF] Watson | Dissolution, growth and survival of zircons during crustal fusion: Kinetic principals, geological models and implications for isotopic inheritance[END_REF][START_REF] Schmitt | Acigöl rhyolite field, Central Anatolia (part 1): high-resolution dating of eruption episodes and zircon growth rates[END_REF], when divided by Zr diffusivities in a liquid at relevant magmatic temperatures for the Gangdese samples (e.g., ~3.5×10 -18 cm 2 /s at 740 ˚C, calculated using the diffusion parameters of [START_REF] Mungall | Chemical diffusivities of 18 trace elements in granitoid melts[END_REF], result in R/D values on the order of 2.9×10 5 to 2.8×10 1 cm -1 , which are well within the range where large diffusive isotope fractionations can be expected based on the calculations of [START_REF] Watson | Non-equilibrium isotopic and elemental fractionation during diffusion-controlled crystal growth under static and dynamic conditions[END_REF]; (ii) the initial drop in isotopic composition away from the geochemical core of a growing zircon crystal -which would produce the W-shaped pattern -would be extremely steep as a function of radial position for the large R/D and Keq values involved during zircon growth. This means that a clear 'W' shape profile would be difficult -if not impossible -to resolve at the spatial resolution of the [START_REF] Guo | Significant Zr isotope variations in single zircon grains recording magma evolution history[END_REF] LA-ICP-MS analyses, potentially leading to convoluted signals and poorly-resolved W shapes such as those observed in several of their profiles.

Given these considerations, the observations made by [START_REF] Guo | Significant Zr isotope variations in single zircon grains recording magma evolution history[END_REF] cannot in fact rule out a kinetic control on the measured δ 94/90 Zr and calculated αzircon-melt values. This, in addition to the much smaller Δ 94/90 Zrzircon-melt found from our ab initio calculations, clearly point to the [START_REF] Guo | Significant Zr isotope variations in single zircon grains recording magma evolution history[END_REF] profiles as being kinetically controlled rather than produced via equilibrium isotope fractionation.

Unlike the FC-1 and Gangdese arc zircon, which formed from silicate melts, the Mud Tank zircon studied by [START_REF] Tompkins | Zirconium stable isotope analysis of zircon by MC-ICP-MS: methods and application to evaluating intra-crystalline zonation in a zircon megacryst[END_REF] is a ~732 Ma megacryst crystallized from a carbonatite melt [START_REF] Crohn | The Mud Tank Carbonatite, Strangways Range, central Australia[END_REF]. This zircon displays remarkably homogeneous δ 94/90 Zr values, suggesting a lack of equilibrium and/or kinetic isotopic fractionation between zircon and melt. The lack of equilibrium fractionation is well explained by our ab initio calculations (section 4.2), but the lack of kinetic effects deserves further discussion.

Differences in kinetic fractionation effects might stem from the peculiar composition of carbonatite magmas, the diffusive properties of Zr, and/or the timescales of megacrystic zircon growth. For instance, if the diffusivity of Zr in carbonatite melts is much larger than in silicate melts, such that any diffusive boundary layers in a crystallizing carbonatite are rapidly homogenized relative to the zircon growth velocity (i.e., negligible v/D), this would lead to insignificant kinetic fractionation effects to be recorded in the solid [START_REF] Watkins | Kinetic Fractionation of Non-Traditional Stable Isotopes by Diffusion and Crystal Growth Reactions[END_REF].

Alternatively, the contrasting abundance of SiO2 between carbonatite and silicate melts could influence zircon formation. In silicate melts, Si is in the majority of cases the most abundant cation, and so zircon growth should be limited by the supply of Zr to a growing seed, rather than by Si availability and diffusivity. On the other hand, in carbonatite melts, SiO2 can be as low as 2-5 wt.%, and thus its availability may also play a role in the kinetics of zircon crystallization. If Zr diffuses faster than Si in carbonatite magmas as it does in silicate melts (Fig. 11), then zircon growth in carbonatites could be controlled by Si diffusion rather than Zr diffusion, thus allowing more time for Zr to maintain isotopic equilibrium throughout the melt and with the crystallizing phases. To our knowledge, however, Zr and Si diffusion coefficients have so far not been constrained in carbonatite melts, so the plausibility of this hypothesis cannot yet be verified. These remarks underline the determining importance of the status of an element (trace vs. major) on its fractionation properties. This importance was emphasized and discussed in a very different situation, for the fractionation of trace (Ba, Sr, B, Li) and major (Mg, Ca) elements during the precipitation of carbonates (calcite, aragonite) at low temperatures [START_REF] Mavromatis | Kinetics of Mg partition and Mg stable isotope fractionation during its incorporation in calcite[END_REF].

In summary, although no experimental data yet exist for determining Zr in silicate liquids, the exploratory calculations made here suggest that reasonable combinations of v, Zr, and Keq have the potential to produce per-mil δ 94/90 Zr values in crystallizing magmas. As such, kinetic effects are at present the most likely explanation for the large (and variable) massdependent Zr isotope variations observed amongst zircon and baddeleyite from silicate melts (see Fig. 6), and might also explain the lack of isotopic zoning in megacrystic zircon from the Mud Tank carbonatite.

Future directions and concluding remarks

The limited magnitude of Zr isotope fractionation predicted by our ab initio calculations and inferences on NFS effects are consistent with the homogeneous δ 94/90 Zr values of some zircon localities, such as Mud Tank, 91500 and Plesovice. However, these fail to explain the large isotope variations observed in several other igneous systems such as the Duluth Complex, Hekla volcano, and the Gangdese arc batholith (Fig. 6). On the other hand, Zr diffusion in response to concentration gradients developed within silicate liquids by rockforming phase crystallization, and/or zircon nucleation and growth, seems to provide a more plausible mechanism to explain large isotopic effects in igneous minerals.

Building from the observations made here, future studies should target the assumptions and testable predictions of the diffusive fractionation hypothesis. In particular: 1) Experimental determination of Zr values in silicate liquids, as well as Zr values and the relative diffusivities of Zr and Si in carbonatite liquids, will be needed to evaluate the hypotheses outlined here regarding the Duluth gabbro and Mud Tank systems.

2) Major phases with low Zr partition coefficient (e.g., plagioclase, clinopyroxene) should display diagnostic concentration and isotopic profiles, with both [Zr] and δ 94/90 Zr values increasing in the direction of solid-melt boundary propagation. Coupled with adequate modeling, spatially-resolved (e.g., LA-ICPMS) investigations could provide a robust test of the diffusion hypotheses put forth here.

3) Zircon and baddeleyite formed near the interface between major phases and melt should display elevated δ 94/90 Zr values, while crystals formed further away from the boundary may record either unfractionated or lower δ 94/90 Zr values. In-situ isotope measurements of Zr-bearing phases using LA-ICPMS or SIMS techniques can be used to test such spatial distributions of δ 94/90 Zr values in a petrologic context.

Another aspect needing further work is the possibility that some rock forming minerals which can incorporate a significant amount of Zr (e.g., clinopyroxene, ilmenite), could fractionate Zr isotopes. Ab initio calculations could not be performed for these minerals, owing to two main limitations: (i) First is the general lack of knowledge on the host site of Zr in these mineral structures. Taking the example of pyroxene, Zr 4+ could substitute either in the M1 or M2 site [START_REF] Morimoto | Nomenclature of Pyroxenes[END_REF] for elements like Mg 2+ , Ca 2+ , Na + or Li + . This would require, at least Zr 4+ -Mg 2+ substitution coupled to, for instance, two Al 3+ -Si 4+ substitutions (or one Al 3+ -Si 4+ plus one Na + -Ca 2+ ) to preserve charge balance. The large number of possible substitutions, challenges associated with proper computation of the isotope properties of such non-isoelectronic cation substitution, as well as the extremely large unit-cell to consider to obtain realistic results (to account for the dilution of Zr) make such calculations prohibitively difficult. (ii) Assuming sufficient computational power, the lack of constraints on the vibrational (and other) properties of Zr in these minerals would make it virtually impossible to validate the model results. Based on the large variety of Zr environments tested here, and on the very low fractionations obtained in all instances, we do not expect these additional calculations to yield fractionation coefficients that are sufficiently large to explain the natural observations by equilibrium mechanisms alone. Yet, to test these assumptions and evaluate whether Zr isotope fractionation is associated with Zr uptake into major rock forming minerals, future studies could perform high-precision Zr isotope analyses on clean mineral separates. Such investigations in rock samples from well-characterized igneous systems could help establish whether such fractionation factors are dependent on the composition, temperature, water content, crystallization timescales and/or redox conditions of the melt.

The nascent field of mass-dependent Zr stable isotope variations is rapidly expanding.

As a growing number of studies document the natural variations of δ 94/90 Zr in igneous, metamorphic and sedimentary systems, elucidating the mechanisms responsible for the partitioning of Zr isotopes will become increasingly critical. Our work suggests that kinetic processes, particularly diffusive fractionation, are likely to play a dominant role in driving variability in Zr-bearing phases in high-T environments. This is a critical first step which, together with future studies, will help establish a robust framework to interpret the petrogenetic significance of fractionated Zr isotopic compositions. [START_REF] Ibañez-Mejia | Extreme Zr stable isotope fractionation during magmatic fractional crystallization[END_REF], Hekla volcano bulk rocks [START_REF] Inglis | Isotopic fractionation of zirconium during magmatic differentiation and the stable isotope composition of the silicate Earth[END_REF], Mud Tank carbonatite zircon megacryst [START_REF] Tompkins | Zirconium stable isotope analysis of zircon by MC-ICP-MS: methods and application to evaluating intra-crystalline zonation in a zircon megacryst[END_REF] and Gangdese arc zircons [START_REF] Guo | Significant Zr isotope variations in single zircon grains recording magma evolution history[END_REF]. For many of the systems studied, equilibrium massdependent effects do not seem to be the main driver of isotopic fractionation.

Figures

Figure 8: (Top) Offset in mass-dependent Zr isotope composition after DS data reduction (expressed as δ 9x/90 Zrδ 94/90 Zr, in ‰/amu, with x =1, 2, or 6) and (bottom) NFS induced mass-independent effects (in epsilon notation, 9x Zr= ( 9x/90 Zr)standard+NFS/( 9x/90 Zr)standard-1)*10 4 ) as a function of the δ 94/90 Zr obtained using a 4 ratio minimization DS reduction. The mass-dependency of all ratios, within ±0.01‰/amu in zircons measured to date [START_REF] Ibañez-Mejia | Extreme Zr stable isotope fractionation during magmatic fractional crystallization[END_REF][START_REF] Tompkins | Zirconium stable isotope analysis of zircon by MC-ICP-MS: methods and application to evaluating intra-crystalline zonation in a zircon megacryst[END_REF] indicates that NFS effects were not the dominant driver of isotopic fractionation in these samples. The maximum permissible offset of 0.02‰/amu between δ 9x/90 Zr and δ 94/90 Zr values, provides (red arrows) upper limits of -4.6 91 Zr, +0.4 92 Zr and -5.2 96 Zr on the NFS induced anomalies in these samples. [START_REF] Huang | Effect of melt structure on traceelement partitioning between clinopyroxene and silicic, alkaline, aluminous melts[END_REF], the Zr diffusion parameters in dry melts of [START_REF] Mungall | Chemical diffusivities of 18 trace elements in granitoid melts[END_REF], and the analytical solution to the linear 1-D diffusion equation with a moving boundary of [START_REF] Smith | A mathematical analysis of solute redistribution during solidification[END_REF]. Note the break in vertical and horizontal scales between the two plots. (Top) Illustration of how growth of a major rock-forming phase with a low (<1) Keq for Zr (e.g., clinopyroxene) leads to preferential zircon nucleation within high [Zr]diffusive boundary layers in the melt.

53 Watkins et al. 2017 and references therein). Gray shaded region shows the range of experimental DZr/DSi in dacitic to rhyolitic melts (i.e., SiO2 wt.% from 65 to 75) at 10 kbar pressure (after [START_REF] Baker | Chemical interdiffusion of dacite and rhyolite: anhydrous measurements at 1 atm and 10 kbar, application of transition state theory, and diffusion in zoned magma chambers[END_REF], suggesting that βZr might be non-negligible (between 0.01 and 0.10).

54 a LDA refers to the present calculation. b Int.: calculated Raman cross-section in atomic units. c For zircon, measurements (exp) come from [START_REF] Kolesov | The dynamic properties of zircon studied by single-crystal X-ray diffraction and Raman spectroscopy[END_REF]. For baddeleyite, measurements (exp) correspond to the study of [START_REF] Carlone | Raman spectrum of zirconia-hafnia mixed crystals[END_REF] at 15K. Other experimental data at 300 K, when different by more than 5cm -1 , are listed below. §

The peak at 310.5cm -1 measured at 15K by [START_REF] Carlone | Raman spectrum of zirconia-hafnia mixed crystals[END_REF] is found at 305 cm -1 at 300K, by [START_REF] Carlone | Raman spectrum of zirconia-hafnia mixed crystals[END_REF], [START_REF] Ishigame | Temperature Dependence of the Raman Spectra of ZrO2[END_REF] and [START_REF] Quintard | Comparative Lattice-Dynamical Study of the Raman Spectra of Monoclinic and Tetragonal Phases of Zirconia and Hafnia[END_REF]. § §

The peak at 745cm -1 measured at 15K by Carlone, 1992 corresponds to a broad feature difficult to discriminate, and measured at 757cm -1 by [START_REF] Quintard | Comparative Lattice-Dynamical Study of the Raman Spectra of Monoclinic and Tetragonal Phases of Zirconia and Hafnia[END_REF], and 780cm -1 by [START_REF] Carlone | Raman spectrum of zirconia-hafnia mixed crystals[END_REF], at 300K. The mean force constants are calculated using the relationship A1=1000 (1/mlight -1/mheavy) ħ 2 / (8 kB 2 ) ✕ <F> (see [START_REF] Dauphas | A general moment NRIXS approach to the determination of equilibrium Fe isotopic fractionation factors: application to goethite and jarosite[END_REF].

57

Figure 1 :

 1 Figure 1: Mineral structures modeled in this study: (a) zircon VIII ZrSiO4; (b) baddeleyite VII ZrO2 ; (c) gittinsite VI ZrCaSi2O7 ; (d) sabinaite Na4 VIII Zr2TiC4O16 ; (e) vlasovite Na2 VI ZrSi4O11 . Green: Zr atoms; red: oxygen atoms; deep blue: Si atoms; yellow: Na atoms; light blue: Ti atoms; grey blue: Ca atoms. Roman numerals denote coordination number.

Figure 3 :

 3 Figure 3: Temperature dependence of ln β values (in ‰) for 94 Zr/ 90 Zr in various minerals investigated here. See text for details on calculation methods.

Figure 4 :

 4 Figure 4: Temperature dependence of mineral-zircon isotope fractionation factors (1000 ln ⍺mineralzircon ≈ Δ 94/90 Zrmineral-zircon) for 94 Zr/ 90 Zr for various minerals investigated here.

Figure 5 :

 5 Figure 5: Comparison of the spectral contributions to the zirconium β-factors of (a) zircon and baddeleyite; (b) sabinaite and baddeleyite; (c) vlasovite and gittinsite.

Figure 6 :

 6 Figure 6: Comparison of the zircon-melt fractionation factors (Δ 94/90 Zrzircon-melt= δ 94/90 Zrzirconδ 94/90 Zrmelt, in ‰) predicted by ab initio calculations at 700°C (this work) with those derived from the analyses of igneous samples: Duluth gabbro zircons[START_REF] Ibañez-Mejia | Extreme Zr stable isotope fractionation during magmatic fractional crystallization[END_REF], Hekla volcano bulk rocks[START_REF] Inglis | Isotopic fractionation of zirconium during magmatic differentiation and the stable isotope composition of the silicate Earth[END_REF], Mud Tank carbonatite zircon megacryst[START_REF] Tompkins | Zirconium stable isotope analysis of zircon by MC-ICP-MS: methods and application to evaluating intra-crystalline zonation in a zircon megacryst[END_REF] and Gangdese arc zircons[START_REF] Guo | Significant Zr isotope variations in single zircon grains recording magma evolution history[END_REF]. For many of the systems studied, equilibrium massdependent effects do not seem to be the main driver of isotopic fractionation.

Figure 9 :

 9 Figure9: NFS induced mass-independent effects (expressed using the epsilon notation) as a function of the offset in mass-dependent Zr isotope composition after DS data reduction (expressed as the absolute value of δ 9x/90 Zrδ 94/90 Zr, in ‰/amu, with x =1, 2, or 6).

Figure 10 :

 10 Figure10:[Zr] profiles during crystal growth in a melt, calculated using a clinopyroxene-melt partition coefficient (Keq) of 0.18 and crystal growth velocities ranging from 0.2 to 1.0 cm/Myr. (Left) Concentration profiles developed within the growing solid. (Right) Concentration profiles developed within the liquid, as driven by concurrent solid growth and diffusion. Calculations used the clinopyroxene/melt Keq of[START_REF] Huang | Effect of melt structure on traceelement partitioning between clinopyroxene and silicic, alkaline, aluminous melts[END_REF], the Zr diffusion parameters in dry melts of[START_REF] Mungall | Chemical diffusivities of 18 trace elements in granitoid melts[END_REF], and the analytical solution to the linear 1-D diffusion equation with a moving boundary of[START_REF] Smith | A mathematical analysis of solute redistribution during solidification[END_REF]. Note the break in vertical and horizontal scales between the two plots. (Top) Illustration of how growth of a major rock-forming phase with a low (<1) Keq for Zr (e.g., clinopyroxene) leads to preferential zircon nucleation within high[Zr]diffusive boundary layers in the melt.

Figure 11 :

 11 Figure 11: Efficiency of diffusive isotope fractionation (β) in silicate melts vs. solvent-normalized diffusivity (modified fromWatkins et al. 2017 and references therein). Gray shaded region shows the range of experimental DZr/DSi in dacitic to rhyolitic melts (i.e., SiO2 wt.% from 65 to 75) at 10 kbar pressure (after[START_REF] Baker | Chemical interdiffusion of dacite and rhyolite: anhydrous measurements at 1 atm and 10 kbar, application of transition state theory, and diffusion in zoned magma chambers[END_REF], suggesting that βZr might be non-negligible (between 0.01 and 0.10).

Figure 12 :

 12 Figure 12: ẟ 94/90 Zr profiles developed in a silicate melt due to kinetic (diffusive) isotopic fractionation in the liquid driven by crystallization of plagioclase (top) and clinopyroxene (bottom). (Left) Effects of variable solid growth velocity (v) at a constant Zr=0.05. (Right) Effects of variable Zr at a constant solid growth velocity (v= 1 cm/Myr). Preferential zircon saturation in high [Zr] regions of the diffusive boundary layers (i.e., nearest to the mineral-melt interface; see Fig. 10) will lead to the formation of isotopically heavy zircon.

Table 2 :

 2 Comparison of model and measured Raman phonon frequencies in Baddeleyite and Zircon. LDA a (cm -1 ) Int bExp c (cm -1 ) LDA a (cm -1 ) Int b Exp c (cm -1 )

	Baddeleyite
	Ag	103.43	6	101	Bg	177.29	6	na
	Ag	183.29	13	177	Bg	227.49	6	222.5
	Ag	190.79	17	189.5	Bg	315.96	2	na
	Ag	319.92	37	310.5 §	Bg	335.84	30	331.6
	Ag	349.61	23	348.1	Bg	388.08	46	378
	Ag	386.31	17	384.9	Bg	495.18	24	501.1
	Ag	473.96	180	474.1	Bg	537.47	32	537.5
	Ag	552.48	53	556.3	Bg	607.39	102 614.4
	Ag	638.66	145	635.0	Bg	754.83	22	745 § §
	Zircon
	A1g	437	59	439	B2g	257	0.5	266
	A1g	960	25	974	Eg	195	2	202
	B1g	224	3	215	Eg	225	3	225
	B1g	390	3	393	Eg	362	24	356
	B1g	626	1	641	Eg	535	0.1	546
	B1g	1005	251	1008	Eg	914	0.02 923

Table 3 :

 3 Comparison of model and measured IR phonon frequencies (TO: transversal optic; and LO: longitudinal optic modes a ) in baddeleyite and zircon.

		LDA b (cm -1 ) Int. c	Exp d (cm -1 )		LDA b (cm -1 )	Int. c	Exp d (cm -1 )
	Baddeleyite TO-LO		TO		TO-LO		TO
	Au	183	0.1	na	Bu	229-241	16-7	224 §
	Au	245-246	1	na	Bu	307	1	na

Table 4 :

 4 Fits of 10 3 ln β( 94/90 Zr)mineral for 135-1200°C, based on A1x+A2x 2 +A3x 3 with x=10 4 /T 2

	mineral	<F> (N/m) §	A1	A2	A3
	gittinsite	407.8	84.88	-61.74	99.92
	vlasovite	407.8	84.88	-63.52	101.78
	zircon	392.7	81.75	-50.39	82.53
	sabinaite	380.4	79.18	-53.31	106.62
	baddeleyite	372.4	77.53	-41.00	49.53
	§

Acknowledgement

This research was supported by NSF-EAR grants 1823748 (to MIM) and 1824002 (to FT) and start-up funds to MIM provided by University of Rochester and to FT provided by Caltech.

This work was performed using HPC resources from CALMIP (Calcul en Midi-Pyrénées; Grant 2019-P1037). We are grateful for critical evaluations from Editor Mark Rehkämper, Jörn-Frederik Wotzlaw, and two anonymous reviewers, which helped to improve the clarity of this manuscript.

Tables