Isotopic fractionation of zirconium during magmatic differentiation and the stable isotope composition of the silicate Earth

Abstract High-precision double-spike Zr stable isotope measurements (expressed as δ94/90ZrIPGP-Zr, the permil deviation of the 94Zr/90Zr ratio from the IPGP-Zr standard) are presented for a range of ocean island basalts (OIB) and mid-ocean ridge basalts (MORB) to examine mass-dependent isotopic variations of zirconium in Earth. Ocean island basalt samples, spanning a range of radiogenic isotopic flavours (HIMU, EM) show a limited range in δ94/90ZrIPGP-Zr (0.046 ± 0.037‰; 2sd, n = 13). Similarly, MORB samples with chondrite-normalized La/Sm of >0.7 show a limited range in δ94/90ZrIPGP-Zr (0.053 ± 0.040‰; 2sd, n = 8). In contrast, basaltic lavas from mantle sources that have undergone significant melt depletion, such as depleted normal MORB (N-MORB) show resolvable variations in δ94/90ZrIPGP-Zr, from −0.045 ± 0.018 to 0.074 ± 0.023‰. Highly evolved igneous differentiates (>65 wt% SiO2) from Hekla volcano in Iceland are isotopically heavier than less evolved igneous rocks, up to 0.53‰. These results suggest that both mantle melt depletion and extreme magmatic differentiation leads to resolvable mass-dependent Zr isotope fractionation. We find that this isotopic fractionation is most likely driven by incorporation of light isotopes of Zr within the 8-fold coordinated sites of zircons, driving residual melts, with a lower coordination chemistry, towards heavier values. Using a Rayleigh fractionation model, we suggest a αzircon-melt of 0.9995 based on the whole rock δ94/90ZrIPGP-Zr values of the samples from Hekla volcano (Iceland). Zirconium isotopic fractionation during melt-depletion of the mantle is less well-constrained, but may result from incongruent melting and incorporation of isotopically light Zr in the 8-fold coordinated M2 site of orthopyroxene. Based on these observations lavas originating from the effect of melt extraction from a depleted mantle source (N-MORB) or that underwent zircon saturation (SiO2 > 65 wt%) are removed from the dataset to give an estimate of the primitive mantle Zr isotope composition of 0.048 ± 0.032‰; 2sd, n = 48. These data show that major controls on Zr fractionation in the Earth result from partial melt extraction in the mantle and by zircon fractionation in differentiated melts. Conversely, fertile mantle is homogenous with respect to Zr isotopes. Zirconium mass-dependent fractionation effects can therefore be used to trace large-scale mantle melt depletion events and the effects of felsic crust formation.


INTRODUCTION
Over the past twenty years the field of non-traditional stable isotope geochemistry has proliferated different isotope systems developed and applied to various geo-and cosmo-chemical problems (see reviews by Johnson et al., 2004;Teng et al., 2017). The stable isotopes of the highfield strength elements (HFSE; Zr, Hf, Ti) have received relatively little attention to date, however, with only Ti being explored in terms of natural mass-dependent stable isotope variations (e.g. Millet et al., 2016;Greber et al., 2017;Deng et al., 2018Deng et al., , 2019. The HFSE are of particular interest for understanding processes such as melt genesis because they behave incompatibly during partial mantle melting events (Woodhead et al., 1993;Johnson, 1998), and, as a consequence, are highly enriched in crustal rocks (for Zr $ 200 ppm; Rudnick and Gao, 2003) relative to the mantle (for Zr $ 4 ppm; McDonough and Sun, 1995). Furthermore the HFSE, including Zr, are relatively insoluble in sub-critical aqueous fluids and are highly refractory, meaning that they are resistant to later modification or resetting by metamorphism or alteration (Brenan et al., 1994). As such, Zr can serve as a high-fidelity tracer of mantle depletion events and the genesis of ancient igneous rock and crust (Condie, 2005).
Zirconium isotope geochemistry has seen limited exploration in terms of mass-dependent variations. Almost all studies have been limited to exploring nucleosynthetic anomalies (e.g. Akram et al., 2015;Akram and Schö nbächler, 2016). The only mass-dependent study of Zr within natural samples was limited to a series of geological standard reference materials (SRM), spanning bulk major element compositions between basalt and granite ([SiO 2 ] = 49.9 wt% to 69.9 wt%) (Inglis et al., 2018). It was found that d 94/90 Zr IPGP-Zr , which represents the deviation of the 94 Zr /90 Zr ratio from the IPGP-Zr standard in parts per thousand, increases with increasing SiO 2 content from 0.044 ± 0.044‰ to 0.186 ± 0.035‰, demonstrating not only that resolvable Zr stable isotope variations exist in nature, but also that Zr isotopes can potentially serve as sensitive tracers of magmatic differentiation. If this is the case then it is anticipated that Zr stable isotopes can be used to tracer ancient magmatic differentiation and provide insight into the emergence and evolution of differentiated crustal lithologies throughout Earth history.
To examine the suitability of Zr stable isotopes as a tracer of magmatic processes and to estimate the Zr isotope composition of the terrestrial mantle, new high-precision mass-dependent Zr isotope data are presented for a range of mid-ocean ridge (MORB) and ocean island (OIB) basalts, along with two well-characterised magmatic differ-entiation suites from Kilauea Iki lava lake, Hawaii, and Hekla volcano, Iceland.

Ocean island basalts (OIB)
To examine possible isotopic heterogeneities that may be preserved within the mantle sampled by OIB, fourteen whole rock OIB lavas were analysed. These OIB were selected from different intraplate volcanic islands from the Atlantic and Pacific oceans and have been previously studied for major & trace element and isotope systematics. Samples were selected to give a range of different radiogenic isotope ''mantle flavours" such as EM-1, EM-2 and HIMU. Specifically, the OIB samples analysed here comprise: three samples from São Nicolou island of the Cape Verde archipelago   (Millet et al., 2008;Mourão et al., 2012) (Hart and Jackson, 2014), two samples from El Hierro of the Canary Islands (EH15 and EH17) (Day et al., 2009(Day et al., , 2010, one sample from Tristan da Cunha (TDC -BM1962, 128[114]) (Baker et al., 1964) and one sample from Gough Island (ALR 40G) (le Roex, 1985).
Where available the literature major and trace element data for these samples are presented in Table 1 and Supporting Information Table 1. For a number of samples trace element data was not available. In this instance trace element abundance measurements were performed on the same sample aliquots used for Zr isotope measurements and are presented alongside the literature elemental data in Supporting Information Table 1. 2.1.3. Magmatic differentiation suites 2.1.3.1. Hekla volcano (Iceland). The Hekla volcano is an active volcanic fissure situated in the South Iceland Volcanic Zone. The volcano has been active since 1104 A.D. and has seen 18 historic eruptions, with the most recent occurring in 2000 A.D (Hö skuldsson et al., 2007). The volcano is one of the most active on Iceland and, unlike most others on the island, produces an array of lavas spanning a bulk compositional range between primitive basalt to evolved rhyolite ($45-72 wt% SiO 2 ) (Sigmarsson et al., 1992;Chekol et al., 2011). Differing models of magmatic evolution for Hekla have been proposed. Sigmarsson et al. (1992) suggested a multi-stage model whereby basaltic melt pools towards the bottom of a shallow crustal magma chamber undergoes fractional crystallisation and evolve along a line of liquid descent to basaltic andesite. This melt triggers melting of crustal lithologies to form a magma of dacitic composition, which in turn mixes with the basaltic andesite magma to form a melt of andesitic composition. This hybrid melt progresses along a line of liquid descent to generate the rhyolites. This model has been recently challenged based on coupled radiogenic isotope data and geophysical observations, that the differentiation of the basaltic parent melt to intermediate compositions of basaltic-andesites could be accounted for by simple closed system fractional crystallisation (Chekol et al., 2011). This model differs from that of Sigmarsson et al. (1992) in that it suggests that these early basic to intermediate melts were generated within a deep seated magma chamber, this in turn supplied magma to a series shallower, crustal melt lenses which underwent a process of assimilation fractional crystallisation (e.g. DePaolo, 1981). Despite these two differing models no fundamental difference exists with respect to the application of the Hekla suite to examine the magmatic behaviour of Zr isotopes owing to the nature of the potential crustal assimilant.
In recent years Hekla has been used as a natural laboratory to examine the mass-dependent stable isotope behaviour of several elements (Schuessler et al., 2009;Savage et al., 2011;Chen et al., 2013;Prytulak et al., 2017). The nine Hekla samples analysed here are taken from the study of Savage et al. (2011) and represent magmatic differentiates with SiO 2 contents spanning from 46.4 to 72.1 wt% and MgO contents of between 0.08-5.6 wt%.
2.1.3.2. Kīlauea Iki (Hawaii). The Kīlauea Iki lava lake formed during the 1959 eruption of the Kīlauea volcano, where large volumes of basaltic lava was erupted and ponded within the Kīlauea Iki pit crater and underwent relatively rapid cooling and differentiation as part of a closed system (Helz and Thornber, 1987). The resulting differentiation sequence comprises lower olivine cumulates to evolved andesitic melt segregation veins (Helz, 1987). As part of an extended campaign, the United States Geological Survey (USGS) conducted an extensive drilling program of the solidifying lava lake. This campaign resulted in the recovery of >1500 m of drill core, which sampled a range of lithologies between the most primitive to most evolved magmatic differentiates. These samples have been wellcharacterised in terms of petrography, major and trace element geochemistry (Helz, 1987;Helz et al., 1994;Helz and Taggart, 2012).
Much like Hekla, Kīlauea Iki has been used as a case study for examining the effect of magmatic differentiation on various stable isotope systems (Tomascak et al., 1999;Teng et al., 2007;Teng et al., 2008;Chen et al., 2013;Badullovich et al., 2017;Kato et al., 2017;Amsellem et al., 2018). The nine samples analysed as part of this study were originally characterised by Helz et al. (1994) and Helz and Taggart (2012), but have been analysed as part of the aforementioned stable isotope studies. They represent a range of different SiO 2 contents of between 46.7 and 57.1 wt% and MgO contents of 2.4 and 13.5 wt%.

Zirconium isotope measurements
Zirconium isotope measurements were performed using the method described in Inglis et al. (2018) that is summarised here. Depending on the Zr content of the samples, between 20 and 150 mg of whole-rock sample powder was weighed directly into clean Savillex PFA Teflon squared bodied 7 mL beakers, herein referred to as Teflon bombs, to obtain $1-2 lg of sample Zr. The sample powders were spiked prior to dissolution with a 91 Zr-96 Zr double spike at a Zr proportion of 43:57 spike to sample. Initial sample decomposition was achieved by the addition of concentrated HF (29 M) and HNO 3 (16 M) acids at a ratio of 1:2 and heating in high walled hotplate at 165°C for 5 days. This dissolution technique has been demonstrated to fully dissolve zircon grains and high silica rocks efficiently for Zr isotope measurements (Inglis et al., 2018). Following this initial decomposition, the samples were evaporated to dryness and residues were re-dissolved in 6 M HCl and 16 M HNO 3 to ensure complete dissolution of any fluoride phases that may have formed during the initial HF step. To ensure that the data obtained for the samples dissolved by this method were not biased by incomplete dissolution of refractory phases (i.e. zircon) repeat dissolution tests were made using high pressure-temperature steel digestion vessels (Parr bombs). For the samples digested in this manner, the sample powder was weighed into 3 ml Savillex hex beakers, spiked accordingly before concentrated HF and HNO 3 acids were added at a ratio of 1:2. These sample beakers were then placed inside the Teflon Ò interior chamber of the Parr bomb and $3 ml of 1:2 HF:HNO 3 were added to enable vapour exchange between the Teflon chamber and sample beakers when subjected to high P-T. The sealed Parr bombs were then placed inside an oven and heated  Hart and Jackson (2014) at 220°C for 7 days, after which samples were refluxed in 6 M HCl and 16 M HNO 3 prior to column chemistry. The isolation of Zr from the sample matrix was achieved using a two-stage ion exchange chromatography. The first stage of this chemistry was performed on BioRad TM Poly-Prep columns filled with AG1-X8 (200-400 mesh) anion exchange resin. Samples were loaded onto the column and the matrix eluted in 4 M HF, while the fraction containing the sample Zr was recovered in 6 M HCl + 0.01 M HF. This sample fraction was then evaporated to dryness and re-dissolved in 12 M HNO 3 . For the second stage of the chemistry procedure the sample solutions were loaded onto BioRad TM PolyPrep columns packed with 1 mL of Eichrom DGA resin. The matrix was then eluted from the column, first in 12 M HNO 3 and then in 3 M HNO 3 . Finally, the Zr fraction was recovered from the column in 3 M HNO 3 + 0.2 M HF. The Zr fractions were evaporated to dryness before being brought back into solution in 0.5 M HNO 3 prior to mass spectrometry.
All Zr isotope measurements were performed on a Thermo Scientific Ò Neptune Plus multiple-collector inductively coupled plasma mass spectrometer (MC-ICPMS) housed at the Institut de Physique du Globe de Paris, France (IPGP). Samples were introduced to the instrument, which was run in low mass-resolution mode, via a quartz SIS spray-chamber and PFA 50 lL min À1 nebuliser at a running concentration of 200 ng mL À1 total Zr solution, which gave a total Zr beam intensity of $20 V. The intensity of all the isotopes of Zr were collected ( 90 Zr + , 91 Zr + , 92 Zr + , 94 Zr + and 96 Zr + ) alongside the ones of 88 Sr + and 95 Mo + -the latter in order to monitor for interferences of 94 Mo + and 96 Mo + on 94 Zr + and 96 Zr + , respectively. Between all samples, baseline corrections (on peak zeros) were performed in clean 0.5 M HNO 3 , the effects of which were also negligible. After each measurement a 15 min wash-out was performed with clean 0.5 M HNO 3 in order to return to background beam intensities. Because it is known that Zr is not particular stable in dilute HNO 3 tests were made to examine the effect of running samples in 0.5 M HNO 3 + 0.1 M HF. For these tests, all samples, standards, on peak zeros and wash solutions were prepared in 0.5 HNO 3 + 0.1 M HF and introduced into the instrument using a Savillex PFA cyclonic spray chamber and an inert sapphire injector. Instrument operating parameters are the same as when running in 0.5 M HNO 3 except washout times were reduced to 7 min. Several total procedural blanks were measured as part of this study. These were found to range between 100 and 200 pg total Zr, which is negligible when compared to the $<2 lg of sample Zr processed.
All of the Zr isotope data presented (given in Table 1) represent the average of four individual sample measurements, with the reported error being the twice standard deviation (2 sd) of these four measurements, except for the case of the standard reference material data discussed in 3.1, where individual measurements are used. Data was reduced via a series of double-spike inversion calculations using the double-spike data reduction tool IsoSpike (Creech and Paul, 2015). All data is reported in standard delta notation as d 94/90 Zr IPGP-Zr , which represents the 2.2.2. Trace element abundance measurements Trace element abundance measurements were performed for all twelve of the MORB samples and eight of the OIB samples. Roughly 50 mg of sample powder was weighed directly into clean Teflon bombs or round bottomed 7 mL Savillex Ò beakers, added to which was 1 mL 29 M HF and 2 mL 16 M HNO 3 . These were then heated in a high walled hotplate for four days at 130°C, before being evaporated to dryness, re-dissolved in 2 mL of 6 M HCl and refluxed at 130°C for a further three days; this final HCl reflux was to ensure that all fluoride phases formed by the addition of HF were fully brought into solution. Prior to analysis, samples were once again evaporated to dryness and brought back into solution in 0.5 M HNO 3 and diluted by a factor of 2000.
Trace element concentrations were measured using an Agilent 7900 ICP-QMS housed at the IPGP in low resolution mode. Sample introduction was achieved with a micro-nebulizer (MicroMist, 0.2 mL/min) through a Scott spray chamber. Masses between 23 (Na) and 75 (As) were measured using a collision-reaction interface with helium gas (5 mL/min) to remove polyatomic interferences. Scandium and indium internal standards were injected after inline mixing with the samples to correct for signal drift and matrix effects. A mix of certified standards was measured at concentrations spanning those of the samples to convert count measurements to concentrations in the solution. Uncertainties on sample concentrations are calculated using algebraic propagation of blank and sample counts uncertainties. Additionally, BHVO-2 was analyzed as an unknown. Comparison of the values obtained for this reference material with certified concentration data demonstrates an external reproducibility of between 5-10% for all elements given in Supporting Information Table 1. Fig. 1. Data for the standard reference materials basalt BHVO-2 and granite GA. Data points represent individual measurements of five separate repeat dissolutions for BHVO-2 and two repeat dissolutions for GA. The hollow symbols for both BHVO-2 and GA represents data presented in (Inglis et al., 2018), while the filled symbols are data generated as part of this study. For both samples, the grey solid line represents the mean of n and the shaded box represents the 2sd of n. Fig. 2. The Zr stable isotope composition (d 94/90 Zr IPGP-Zr ) of products of igneous differentiation from basaltic magmas (Kīlauea Iki and Hekla, OIB and MORB), and standard reference materials, analysed as part of this study. Each data point represents four measurements (n) of the same sample aliquot, with the error bar being the 2sd of n. The grey solid line surrounded by the grey box represents the primitive mantle value and the 2sd uncertainty, as discussed in the text.

Standard reference materials
As part of this study a series of external standard reference materials (SRM) were analysed in order to validate the accuracy of data collected in the campaign against previously published values for these SRM (i.e. Inglis et al., 2018). Two separate repeat dissolutions and ion exchange separations of basalt BHVO-2 were carried out, as well as a single dissolution of the granite GA. Additionally the rhyolite, RGM-1 was also analysed for the first time for Zr stable isotopes. The data for these standard reference materials is provided in Table 1.
The individual d 94/90 Zr IPGP-Zr values for repeat measurements of BHVO-2 and GA that have been analysed here and elsewhere (Inglis et al., 2018) are presented in Fig. 1

Zirconium isotope composition of MORB, OIB and igneous differentiates
The Zr isotope data for the MORB and OIB samples are given in Table 1 and presented alongside all other Zr isotope data from this study in Fig. 2

The Zr stable isotope composition of MORB and OIB
The MORB and OIB samples have been analysed to determine the Zr isotopic composition of mantle derived basaltic melts. Of the two sample types, the OIB show a restricted range of d 94/90 Zr IPGP-Zr values (mean d 94/90 Zr -IPGP-Zr = 0.046 ± 0.037‰; 2sd, n = 13) and all plot within analytical uncertainty of one another (Fig. 2), and of previous Zr isotope measurements of the Hawaiian basalt SRM BHVO-2 (0.044 ± 0.044‰; Inglis et al., 2018). It is apparent that these samples are isotopically homogenous, with no variation between geographical area or relationship with isotopic mantle flavour (e.g., HIMU, EM). Based on this observation, we suggest that the fertile mantle source regions sampled by OIB is uniform with respect to Zr stable isotope composition.
The MORB samples studied here have been grouped accordingly to their chondrite normalized La/Sm ratio (La/Sm N ) based on the classification of (Schilling, 1975;Schilling et al., 1983). By this scheme samples with La/ Sm N < 0.7 are termed normal or N-MORBs and are suggested to represent products of melting of the depleted MORB mantle (DMM). Samples with La/Sm N > 1.8 are termed enriched or E-MORBs, whilst samples with La/ Sm N ratios between the two end members are classed as transitional or T-MORB. This grouping assumes that both T-and E-MORBs represent increased melting of a more fertile, primitive mantle source.
The variation of d 94/90 Zr IPGP-Zr values seen within the MORB samples (-0.045 ± 0.018 to 0.074 ± 0.023‰ for d 94/90 Zr IPGP-Zr ) is resolvable outside of analytical uncertainty for many of the samples. The MORB samples ana-  (Schilling et al., 1983). The error bars on individual isotope measurements represents the 2sd of n. The solid horizontal black line surrounded by the grey box represents the estimate of the primitive mantle (0.048 ± 0.032‰; 2sd, n = 48). lyzed as part of this work were selected to give a broad geographical distribution and thus representative of different classifications of MORB type (e.g. N-, T-and E-MORB; Schilling et al., 1983). The Zr isotopic composition of T-and E-MORB types are within analytical uncertainty of one another (see Fig. 3) and also plot within the range of values seen in OIB. In fact, the mean d 94/90 Zr -IPGP-Zr for T-and E-MORB types (0.053 ± 0.040‰; 2sd, n = 8) is indistinguishable within error of the OIB mean value (0.046 ± 0.037‰; 2sd; n = 13), consistent with the concept that these melts are likely tapping a mantle source region, which, much like the ones sampled by OIB, is isotopically homogenous with respect to Zr stable isotopes. The origin of incompatible trace element enrichments in MORB is debated; some argue that it is melting of a deeper, more fertile mantle region (i.e. plume source) that can account for such trace element enrichments within E-MORB (Schilling et al., 1983), while others have evoked recycling of oceanic crust and minor amounts of continental material (Hofmann and White, 1982). More recently it has been suggested that incorporation of metasomatically enriched sub-crustal lithosphere could be a viable mechanism for the formation of E-MORB (Workman et al., 2004). The new Zr isotope data for the E-and T-MORB samples presented here does not provide definite insight into the mechanism of E-MORB formation, but does imply melting of a more fertile mantle source given the uniform, unfractionated Zr isotope data displayed for these samples.
The Zr isotope values obtained for the N-MORB samples show a resolvable variation from the T-and E-MORB samples (Fig. 3). With the exception of one sample from the Mid Atlantic Ridge (RD87 DR18-102), the N-MORB samples are offset towards lighter values. These represent the lightest values for any basalts analyzed to date. To a first order, it would thus appear that the melting regime experienced during melting of depleted MORB mantle is responsible for producing an isotopically light melt, relative to T-, E-MORB, or OIB, or that the source of the N-MORB had been previously depleted in the heavier isotopes by continued melt extraction. It is well known that mantle melting and melt extraction processes can exert a strong effect on other stable isotope systems (e.g. Weyer and Ionov, 2007;Huang et al., 2017). However, the lack of Zr isotope measurements of depleted MORB mantle (DMM) samples (e.g. abyssal peridotites), precludes us from commenting exactly as to the origin of this isotopic variation within these samples. One possible scenario is that partial melting and extraction of the melt could preferential incorporate heavy Zr isotopes, progressively driving the source region of these N-MORBs towards lighter values. It is possible that this may result from incongruent mantle melting and citing of isotopically light Zr in an octahedral coordination within mantle orthopyroxene (i.e. Niu, 1997). There is a negative correlation between Zr isotopic composition and Zr/Sm in the N-MORB samples Fig. 4) indicating isotopically light Zr is coupled with high Zr/Sm in N-MORB samples. Orthopyroxene has similar partition coefficients (K D ) for both Zr (0.014-0.036) and Sm (0.016-0.039), with an average Zr/Sm fractionation factor of $0.9, and significant possible extremes (0.3-2.2) (Salters and Longhi, 1999). The correlation with Zr isotopes and Zr/Sm is plausibly consistent with increased influence of orthopyroxene melting in some of the isotopically light N-MORB samples. On the other hand, by investigating komatiites (25-40% partial melts of the mantle), Deng et al. (2018) has highlighted a progressive depletion of both heavy Ti isotopes and incompatible trace elements in the mantle during the late Archean and suggest that this signature is present within the mantle sources of modern MORB samples. The authors considered that such a depletion of heavy Ti isotopes in the depleted MORB mantle could not be solely attributed to partial melting of the mantle (i.e. mantle depletion by the extraction of mafic crust), and that a reworking of mafic crust and recycling of the melt residues into the mantle would be needed to account for the formation of the depleted MORB mantle. This is also a plausible mechanism for the correlation between Zr isotopes and Zr/Sm observed here (Fig. 4), since the melt residues from the crust would preferentially incorporate Sm over Zr (Foley et al., 2002) and light Zr isotopes.

Fractionation of Zr stable isotopes during igneous differentiation
A major objective of this study is to examine the effect of magmatic differentiation on Zr isotopes. This has been tested by studying two separate, cogenetic differentiation suites -Kīlauea Iki and Hekla. The Kīlauea Iki suite displays a restricted range in d 94/90 Zr IPGP-Zr values, with a mean of 0.044 ± 0.023‰; 2sd, n = 9, and displays no apparent covariation with indices of differentiation (Fig. 5). Indeed, the Zr isotope data for this suite are indistinguishable within uncertainty of the mean values for T-and E-MORBs, and also for OIB samples. We conclude from this that the melts produced at the Kīlauea Iki volcano are sampling a mantle source region that, like other OIB, is homogenous with respect to Zr isotopes and that basaltic magmatic differentiation in this setting leads to little or no Zr isotope fractionation. None of the phases formed during continued differentiation of the Kīlauea Iki magmas are responsible for fractionating Zr stable isotopes.
Different systematics are observed for the Hekla sample suite. These samples span a much greater compositional range, and represent a continuous differentiation suite. Consequently, they are ideally suited to examine the effect of more extreme degrees of magmatic differentiation. d 94/90 Zr IPGP-Zr values increases together with the SiO 2 content (Fig. 5) with a steep inflection of the d 94/90 Zr IPGP-Zr values at $65 wt% SiO 2 , whereby values diverge from the range seen in OIB, T-and E-MORB to significantly heavier values, up to 0.52‰, suggesting that there is a control on Zr isotopes during the dacitic stages of differentiation. To understand the mechanism for this isotopic fractionation it is important to consider the elemental behavior of Zr during differentiation. Typically, Zr behaves as an incompatible element during melting, and consequently is concentrated into the melt during continued fractional crystallization (i.e. Woodhead et al., 1993). This elemental behavior of Zr is evident within the Hekla sample suite and is shown on Fig. 5. It is apparent that during evolution between the primitive basaltic melt at >$45 wt% SiO 2 and more evolved dacitic melts at $65 wt% SiO 2 , elemental Zr is continually enriched in the melt. At SiO 2 contents >$65 wt% a steep decrease in [Zr] is observed, which is assumed to represent the precipitation of a Zr rich phase, most likely zircon, which is removed from the residual melt. Indeed, previous studies have suggested the zircon saturation was reached in the Hekla melts (Sigmarsson et al., 1992;Bindeman et al., 2012;Weber and Castro, 2017). Within the Hekla suite it is likely that Zr concentration ([Zr]) would increase linearly until the point of saturation is reached. Owing to the lack of samples spanning the compositional range of SiO 2 contents between 59.64 -68.31 wt %, within which the point of Zr saturation is likely to occur (Watson and Harrison, 1983), this has been estimated by plotting trendlines pre and post zircon saturation and using the intercept of these two lines to infer the point of Zr saturation (Fig. 4). Based on this, the point of Zr saturation within the Hekla melt is suggested to occur at $66.75 wt% SiO 2 . Given the density (q) of zircon ($4.5 g cm À3 ) relative to that of rhyolitic melts ($2.2 g cm À3 ; (Murase and McBirney, 1973)), it is plausible that after the appearance of zircon on the solidus a degree of crystal settling likely occurred, and the zircon was partitioned into the cumulate phase, thus removing it from the residual melt. Strikingly, there is a strong relationship between the behavior of elemental [Zr] and d 94/90 Zr IPGP-Zr , whereby the point of suggested zircon crystallisation is marked by an enrichment of the heavier isotopes of Zr within melts. Precipitation of zircon within the Hekla suite is the likely driver for the large isotopic fractionation seen within these samples.
Theory predicts that equilibrium stable isotope fractionation is largely controlled by the bonding environment of a given element, with the magnitude of the resulting isotopic fractionation being roughly proportional to the bond stiffness between the two equilibrated phases, with the heavier isotopes being concentrated where bond stiffness is greatest (e.g. Urey, 1947;Schauble, 2004). The bond stiffness in such a system can be controlled by a number of factors including oxidation state, the type of elements involved and coordination numbers. Indeed, the shortest and strongest (and thus stiffest) bonds are associated with lower coordination numbers, thus it is predicted that low coordination number lattices will be enriched in the heavier isotopes. The coordination of Zr in silicate melts has been studied by Farges et al. (1991), who suggested that generally Zr is found with 6-fold coordination in silicate melts, but an increase in 8-fold coordinated Zr accompanying Zr saturation is noted. Furthermore, Zr K edge X-ray absorption near edge structure measurements have confirmed the 8fold coordination of Zr ions within the lattice of zircons, relative to melts with lower coordinated Zr (Tobase et al., 2015). In this scenario the crystallisation of 8-fold coordinated Zr ions within zircon would incorporate the lighter Zr isotopes, whilst the residual melt is enriched in the heavier isotopes.
Building on this interpretation that zircon crystallisation is responsible for driving the isotopic fractionation seen within the most magmatically evolved rocks presented here, we have modelled the d 94/90 Zr IPGP-Zr evolution of a residual melt during the progressive fractional crystallisation of zircon. A Rayleigh distillation equation is used (Eq. (2)) to calculate the d 94/90 Zr IPGP-Zr of the residual melt decreasing concentration of Zr, which is assumed to reflect continued zircon crystallisation.
Where f represents the fraction of Zr remaining in the melt, a is the isotopic fractionation factor between the precipitating zircon and the melt, and d 94 Zr (initial) represents the starting d 94 Zr of the model. The model presented here is based on the data from the four Hekla samples that are suggested to represent post-zircon crystallization melts (see Fig. 5). These samples were selected because they represent a relatively simple melting and crystallisation regime from one parent melt. From these data the a zircon-melt has been derived empirically, given a distillation curve that best fits the measured data (Fig. 6). Based on this we suggest a a zircon-melt of 0.99950 best describes the results observed here.

An estimate of the Zr isotope budget of the primitive mantle
Despite the evidence for Zr isotope fractionation during igneous differentiation, isotopic fractionation is only apparent at the onset of Zr saturation in melts (>$65 wt% SiO 2 ; Fig. 5). In addition, is has been demonstrated that melt extraction processes occurring within the DMM could possibly impart a control on Zr isotope compositions, but the absence of peridotite data prohibit us from elucidating this process fully. Furthermore, owing to its incompatible nature, it is likely that during mantle melting Zr is strongly partitioned into the silicate melt, suggesting that partial melting should exert little effect on Zr isotope systematics. As such it is appropriate to use basaltic samples to give an estimate of the primitive mantle value for Zr isotopes. Considering the effect of igneous differentiation and MORB melt extraction we can use the mafic samples presented here and by Inglis et al. (2018) to place a constraint on the primitive mantle Zr isotope composition.
Despite covering a wide range of temporal and spatial settings, the data for the basaltic materials (OIB, T-and E-MORB, Kīlauea Iki and the Hekla samples containing <65 wt% SiO 2 ) are all isotopically undistinguishable from one another within analytical uncertainty, ultimately suggesting that the mantle sampled by these melts is isotopically homogenous with respect to Zr isotopes. Following the assumption that the d 94/90 Zr IPGP-Zr of the primitive mantle is estimated based on the average and 2sd of mafic samples, the primitive mantle composition can be estimated at 0.048 ± 0.032‰; 2sd, n = 43. The data used to calculate the primitive mantle value are given in Supporting Information Table 2 and the average and 2sd of these measurements is presented as the black solid line and grey box on Fig. 2. Since Zr is a purely lithophile element with an estimated concentration in the Earth's mantle that is chondritic (e.g. Fig. 6. A Rayleigh distillation model for the evolution of d 94/ 90 Zr IPGP-Zr during the progressive removal of Zr from the residual melt, which is taken to represent zircon crystallisation. Three different distillation curves are given for 3 different a zircon-melt . An empirically derived a zircon-melt of 0.99950 can best account for the Hekla data presented here (orange circles). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Palme and O'Neill, 2014) the isotopic composition of the Earth's mantle represents the isotopic composition of the bulk Earth and we can conclude that the best estimate of the d 94/90 Zr IPGP-Zr of the bulk Earth is equal to 0.048 ± 0.032‰.

CONCLUSIONS
A range of basaltic igneous rocks and differentiates been analysed in order to understand the Zr stable isotope composition of primitive mafic melts and to examine the effect of igneous differentiation on Zr isotopes. The OIB samples display a restricted range of d 94/90 Zr IPGP-Zr values, which are unresolved from previous measurements of basalts (Inglis et al., 2018), or from T-and E-MORB samples. Of the MORB samples, those with a (La/Sm) N ratio < 0.7 display a resolvable variation in Zr isotopes. Because these MORB samples are from melting of mantle source regions that have seen significant prior melt extraction, it is probable that these petrogenetic process exert a strong control on Zr isotopes. To further test this hypothesis the Zr isotopic composition of DMM samples (peridotites) would have to be analysed.
Two igneous differentiation suites have also been examined. It was found that samples from the Kīlauea Iki volcano, which span SiO 2 contents of between 46.7 and 57.1 wt%, showed no resolvable isotopic variation, and that the mean value obtained for these samples is indistinguishable within error from the mean value obtained from the OIB, T-and E-MORB samples. Unlike Kīlauea Iki the samples from Hekla volcano span a much greater range in SiO 2 contents (46.4-72.1 wt%), which correlates with d 94/90 Zr IPGP-Zr values for these samples. Linking the d 94/90 Zr IPGP-Zr values with Zr content of the samples, we conclude that the isotopic fraction seen within the Hekla samples occurs as a result of zircon crystallisation in the melt, whereby isotopically light Zr is incorporated into zircon and removed from the melt, driving this residue towards a heavier isotopic composition.
Based on these observations lavas originating from the effect of melt extraction from a depleted mantle source ((N-MORB) or that underwent zircon saturation (SiO 2 > 65 wt%) are removed from the dataset to give an estimate of the primitive mantle Zr isotope composition of 0.048 ± 0.032‰; 2sd, n = 48. These data show that major controls on Zr fractionation in the Earth result from incongruent melting of the depleted MORB mantle and by zircon fractionation in differentiated melts. Conversely, fertile mantle is homogenous with respect to Zr isotopes. Zirconium mass-dependent fractionation effects can therefore be used as tracers to examine large-scale mantle melt depletion events and the effects of felsic crust formation.

ACKNOWLEDGMENTS
We thank the ERC under the European Community's H2020 framework program/ERC grant agreement # 637503 (Pristine)) and for the UnivEarthS Labex program (no. ANR-10-LABX-0023 and ANR-11-IDEX-0005-02). Parts of this work were supported by IPGP multidisciplinary program PARI, and by Region île-de-France SESAME Grant (no. 12015908). Thibault Sontag and Pascale Louvat are thanked for assistance in the chemistry labs and with the MC-ICPMS facility at IPGP, while Pierre Burckel is thanked for his help with the ICP-MS measurements. Manuel Moreira is acknowledged for providing the MORB samples. This paper was greatly improved by insightful reviews from Fang Huang and two anonymous reviewers alongside careful editorial handling from associate editor Qing-Zhu Yin.