Novel insights from Fe-isotopes into the 1 lithological heterogeneity of Ocean Island Basalts 2 and plume-influenced MORBs

8 The extent of lithological heterogeneity in the Earth’s convecting mantle is highly debated. Whilst the 9 presence of pyroxenite in the mantle source regions of Ocean Island Basalts (OIBs) has traditionally 10 been constrained using the minor-element chemistry of olivine phenocrysts, recent studies have 11 shown that the Ni and Mn contents of primitive olivines are influenced by the conditions of mantle 12 melting, as well as magma chamber processes. Nevertheless, constraining the lithological properties 13 of the mantle is important due to it’s influence on the P-T path followed by solid mantle material 14 during adiabatic ascent, as well as the density of upwelling mantle plumes. We have therefore 15 explored the use of Fe-isotopes as a novel method of tracing lithological heterogeneity in the mantle 16 source regions beneath plume-influenced segments of the global Mid-Ocean Ridge system as well as 17 OIBs. 18 We present new Fe-isotope (δFe) and trace-element data for 26 basaltic glasses from the 19 plume-influenced Galápagos Spreading Centre to investigate the relative roles of pyroxenite and 20 peridotite in the mantle source region of oceanic basalts. Our data reveals significant heterogeneity 21 in the Fe-isotope composition of the Galápagos Spreading Centre basalts (+0.05 +0.25‰ δFe), 22 which correlates with key majorand trace-element parameters (e.g. CaO(8)/Al2O3(8), [La/Sm]n). 23 Application of new models developed to calculate Fe-isotope fractionation during mantle melting, 24 alongside Monte Carlo simulations for melting of a 2-component peridotite mantle, show that this 25 This manuscript has been accepted for publication in Earth and Planetary Science Letters. Please cite this article as: Gleeson et al. (2020). Novel insights from Fe-isotopes into the lithological heterogeneity of Ocean Island Basalts and plume-influenced MORBS. Earth and Planetary Science Letters. https://doi.org/10.1016/j.epsl.2020.116114 variation cannot be caused by changes in melting processes and/or oxygen fugacity of a peridotitic 26 mantle. Instead, our new δFe data is best explained by variations in the proportion of isotopically27 heavy pyroxenite-derived melt that contributes to the GSC basalts, and conclusively shows that 28 lithological heterogeneity exists in the Galápagos mantle plume. Our findings have implications for the 29 moderately-heavy δFe compositions measured in plume-influenced basalts from the Society Islands, 30 Rochambeau Ridges of the Lau back-arc basin, and the FAMOUS segment of the Mid-Atlantic Ridge, 31 which we suggest may also represent contribution from pyroxenite-derived melts. 32


INTRODUCTION
Earth's mantle, and its lithological properties (i.e. relative modal proportions of olivine and pyroxene), 39 have important implications for mantle dynamics and the density of upwelling mantle plumes 40 (Shorttle et al., 2014), the nature and abundance of these components is poorly constrained. 41 Incorporation of recycled oceanic lithosphere into the convecting mantle is widely believed to result 42 in the presence of lithologically-distinct components via high-pressure melting of eclogite 43 (metamorphosed remnants of recycled slabs), and subsequent reaction of these melts with 44 surrounding peridotite (Sobolev et al., 2007). This results in the formation of pyroxene-rich (and highly 45 fusible) components in the mantle (pyroxenite; Sobolev et al., 2007;Yaxley and Green, 1998). Fe-isotope, major-and trace-element composition of our samples to 8 wt% MgO using published 262 mineral-melt trace-element partition coefficients and isotope fractionation factors (Sossi et al., 2016;263 full details can be found in Appendix A). 264 Due to the considerable uncertainties regarding the influence of crystallisation on the Fe-isotope 265 composition of basaltic lavas, other fractionation factors between olivine/pyroxene and the melt 266 phase were tested to ensure that our results are independent of our fractional crystallisation 267 correction (e.g. from Dauphas et al., 2014; Appendix A). In addition, by applying a correction to 8 wt% 268 MgO, within the range displayed by the GSC basalts, rather than Mg#~70 (i.e. in equilibrium with 269 mantle peridotite), which would require significant extrapolation, we minimise the propagated error 270 that results from this correction. 271 Our fractionation-corrected dataset displays considerable δ 56 Fe heterogeneity and confirms that 272 magmatic differentiation does not have a major influence on the variation observed in the δ 56 Fe 273 composition of the GSC basalts. This demonstrates that significant heterogeneity in the composition 274 of primary mantle melts beneath the GSC must exist. Below we consider whether this variability is 275 related to melting processes (e.g. melt fraction, presence of garnet), or due to heterogeneity in the 276 mantle source (oxidised/lithologically distinct components). 277 6.2 MELT FRACTION AND fO2 278 As outlined above, heavy Fe isotopes (i.e. 56  predicted to possess a higher proportion of Fe 3+ than high-fraction melts (Dauphas et al., 2014(Dauphas et al., , 2009. 283 In addition, an increase in the fO2 of the source (and therefore an increase in the Fe 3+ available to 284 enter the melt phase) is predicted to result in a larger fractionation of Fe-isotopes during mantle 285 melting (Dauphas et al., 2014(Dauphas et al., , 2009 The proportion of enriched melt generated in each model, as well as the proportion of melt that is 317 derived from the garnet stability-field, is recorded alongside the deviation between the trace-element 318 composition predicted by each model and the REE composition of each sample. These results are used 319 to generate a probability distribution for the proportion of melt from the enriched mantle component 320 that contributes to the formation of each sample along the GSC ( Fig. 5; Appendix A) and reveal a strong 321 correlation between δ 56 Fe and the mass fraction contribution of enriched melt (Fig. 6a) Sossi and O'Neill, 2017). This is an attractive explanation for the heterogeneity observed in the GSC 340 basalts, due to the very strong correlation observed between δ 56 Fe and [Sm/Yb]n (Fig. 3), as well as 341 the strong correlation between δ 56 Fe and the proportion of melt derived from the garnet-stability field 342 that is estimated by our MCMC models (Fig. 6b). compositions for mantle pyroxenites were tested and the results indicate that the δ 56 Fe composition 381 of the pyroxenite-derived melt that is estimated below represents a maximum value (i.e. by using 382 other pyroxenitic compositions less heterogeneity in the Fe-isotope composition of the mantle source 383 is required to explain our results; Fig S.4). 384 Our Melt-PX based mantle melting model was coupled to a MCMC algorithm (as described above) in 385 order to estimate the proportion of pyroxenite-derived melt that contributes to each GSC basalt ( are assumed to be slightly more oxidised than those of a depleted peridotite (Fig. 4). Nevertheless, to 400 generate the highest δ 56 Fe compositions displayed by the GSC basalts we require the pyroxenitic 401 source beneath the GSC to be isotopically-heavy (δ 56 Fe=+0.18 -+0.20‰; Fig As a result, we suggest that the heterogeneity displayed by the GSC basalts is primarily driven by 406 changes in the amount of pyroxenite-derived melt that is delivered to the sub-ridge magmatic 407 sills/mush region. However, the presence of highly-enriched basalts, found only ~20km away from 408   The influence that these factors have on melt chemistry was investigated using a Markov Chain Monte Carlo algorithm (with 5000 iterations). In each model the maximum upwelling velocity, proportion of enriched material in the mantle source, width from which melts are pooled, and depth to the top of the melt column were randomly sampled within pre-set bounds. For each model the differences between the REE concentrations predicted by the mantle melting model (see below) and the observed REE concentrations in the GSC basalts were calculated. These differences were used to calculate an 'Acceptance Ratio' where 1 represents a perfect match between the model and data and 0 indicates no match between the model and data. The acceptance ratio is calculated using the following method: Where [C] i is the concentration of trace-element i in the model; [C] i is the concentration of the same trace-element in the sample under consideration; [Ĉ] i is the analytical uncertainty for that trace-element; and n is the number of trace-elements under consideration. The results of the 5000 model iterations are then used to construct density distributions for the different parameters (i.e. determine which combination of parameters gives the best match to the trace-element composition of that sample). This method also allows us to construct an estimate (including uncertainties) for parameters such as the proportion of enriched melt that is required to explain the trace-element signature of each sample.

Constraining melt compositions (calculating aggregated melts)
In this section we describe the method by which aggregate melt compositions are calculated following near-fractional melting of a two-component mantle (we assume that the instantaneous composition of the melts formed at each depth are known and/or previously calculated).
In our model we simulate melting of a two-component mantle in a triangular melting region beneath a mid-ocean ridge. As recent studies have shown that melts produced at large horizontal distance from the ridge axis may not be transported to the axial magma chamber and effectively pooled (Behn and Grove, 2015) we split our melting region into an upper triangular region and a lower rectangular region, with the transition at depth d c . This effectively simulates a scenario where melts produced at horizontal distances greater than d c * tan(θ) from the ridge axis are not pooled, and therefore do not contribute to the composition of the erupted melts. θ represents the angle between the mid-plane of the melting region and the base of the lithosphere. This scenario is shown schematically in Fig 8. The aggregated melt compositions are then calculated using the following method: Firstly, we define a height δh, which represents the pressure interval at which the melt-fraction is calculated in pHMELTS and Melt-PX (10M P a). This equates to approximately 309m if a density of 3300kg/m 3 is used for the mantle. We then denote the volume of mantle between height h and h + δh as v i . We define v i in two ways. For depths greater than d c : v i = 2 * δh * d c * tan(θ) and for depths shallower than d c : v i = 2 * δh * (d c − (n + 0.5) * δh) * tan(θ) where n = 0 at z = d c , and represents the number of height increments (δh) above the transition from a triangular to a rectangular melting regime.
Additionally, we also define the voume of the triangular melting region (V top ) and rectangular melting region (V bottom ) as: where Z o is the maximum depth at which melting of either mantle component occurs.
We then ratio the volume at each incremental height to the volume of the triangular melting regime.
For depths shallower than d c : For depths greater than d c : By calculating the ratio between the volume of the triangular melting regime and the volume of the entire melting region, we are then able to raito the volume at each height increment to the total volume of the melting region: Therefore, for depths shallower than d c : and for depths greater than d c : From this we can then calculate the fraction of melt supplied from each depth interval for the depleted (F rac d ) and enriched (F rac e ) components: where U r i represents the relative rate of mantle upwelling at the depth corresponding to point i, and is assumed to be constant at each depth (i.e. no lateral change in the rate of mantle upwelling). If U r i = 1 the no active upwelling is present and the rate of mantle upwelling is equal to that generated in response to plate spreading. F d i and F e i represent the fraction of melt produces at that pressure interval.
The final compositino of the melt is derived by multiplying the instantaneous compositions of the melt from each depth interval (C d i ) with the fraction of melt derived from that depth: where P d represents the total fraction of melt derived from the depleted mantel component and is calculated as: where X d is the mass fraction of depleted material in the mantle source. The δ 56 F e composition of the melt formed at each step was calculated using a mass balance approach where the δ 56 F e of each phase (including the melt) can be calculated if the proportions of each phase, concentration of FeO in each phase, and the isotopic fractionation factors between each phase are known. For simplicity we treat Fe as a trace-element so that is can be assigned partition coefficients during mantle melting. The partition coefficients chosen were those used in the calculations of Williams and Bizimis (2014) so that our results could be directly compared. At each pressure interval in our model theorectical Fe-O bond force constants for the 5 mantle mineral phases were calculated using the method of Sossi and O'Neill (2017). These force constants were sued to define mineral fractionation factors at each step (termed ∆ 56 F e ol−min = δ 56 F e ol − δ 56 F e min ) and melt-olivine fractionation factors (∆ 56 F e melt−ol = δ 56 F e melt − δ 56 F e ol ) were calculated using the Fe-O force constants for basaltic melt that are estimated from NRIXS measurements by Dauphas et al. (2014).
Once these fractionation factors have been calculated the δ 56 F e composition of the melts formed at each pressure increment can be calculated using the following method: For each component the δ 56 F e composition of the melt formed at each pressure interval can be calculated as:

Optimisation of LA-ICP-MS analysis of basaltic glass
Systematic testing of basaltic glass analysis by the new 193 laser-ablation ICP-MS system at the University of Cambridge was required in order to optimise analytical conditions. We carried out systematic tests on spot size, fluence, and repetition rate. Matrix effects were also investigated by varying the calibration material between NIST SRM 612 and BCR-2G standard glasses, but similar results were observed in both cases. Results are presented for a range of trace-elements including rare earth elements (REE), first row transition elements (FRTE), large ion lithophile elements (LILE), and high field strength elements (HFSE).
A similar study on the accuracy and precision of LA-ICP-MS analysis was carried out by Jenner and O'Neill (2012). Their study determined that a precision of ¡4% can commonly be achieved for a suite of 20 elements, when the 'optimal ablation diameter' is used. The precision and accuracy of analysis depends on a number of variables, including the ablation diameter, laser repetition rate and fluence, and other factors kept constant in this analysis (e.g. ablation time, counting times per element; choice of internal and external standard materials Jackson et al., 1992). Counting statistics, dependent on how much material is bought into the ICP-MS, and matrix effects during ablation both control the accuracy and precision of analysis (Jenner and O'Neill, 2012).

Spot size
In order to test the effects of spot size we varied the ablation diameter from 30 to 140 µm at a constant fluence (8J/cm 2 ) and repetition rate (10 Hz). Results for the full suite of REE in BHVO-2G are shown in Fig. S.13. These results reveal that (for most REE) a RSD of < 5% can commonly be achieved when the ablation diameter is > 60µm. Below this value the precision of analysis rapidly decreases (precision of LREE < 10% and HREE < 20% at an ablation diameter of 40µm). A similar effect occurs when the FRTE, LILE, and HFSE, are considered although the decrease in precision is not nearly as dramatic for the FRTE as it is for the REE. The accuracy of analysis is also dependent on the ablation diameter, and recovery is shown to fall closest to 1 in the range 60 − 100µm regardless of calibration material used. Overall, we conclude that the best precision and accuracy of analysis is achieved when the ablation diameter is between 80 and 100µm regardless of calibration material. The drop off in precision occurs at ∼ 60µm for enriched material (i.e. BHVO-2G) but at slightly larger ablation diameters for more depleted material (i.e. BIR-1G).

Fluence
Next we tested the effect of fluence (J/cm 2 ) on the precision and accuracy of our results. Fluence of 3 − 13.19J/cm 2 were investigated. Similarly to ablation diameter, precision is shown to decrease at low fluence. At > 6J/cm 2 the precision of all REE and FRTEs are commonly < 6% and usually < 4%, but this increases slightly at lower fluence (< 14% at 4J/cm 2 ; Fig. S.14). Accuracy is near 100% for analysis at 8 − 10J/cm 2 regardless of calibration material used. These results relate to REE, HFSE, LILE, and FRTE in both depleted and enriched samples (BIR-1G and BHVO-2G) and we therefore conclude that a fluence of 8 − 10J/cm 2 represents the optimal analytical conditions. Whilst lower residuals are sometimes seen at higher fluence (12 − 13.19J/cm 2 ), variations in the ablation characteristics between the calibration material NIST SRM 612 and the unknown glass materials cause a decrease in the accuracy of the results.

Repetition rate
We varied the repetition rate between 5 and 20 Hz whilst using a fluence of 8J/cm 2 and ablation diameter of 100µm. The precision of analysis is observed to decrease with repetition rate, and the lowest RSD values are generally achieved at 20 Hz, although a precision of < 6% is seen for most elements at a repletion rate of 10 Hz (Fig. S.15). An accuracy of near 100% is achieved at a repetition rate of 10-20 Hz, but the concentration Fig S.14: Precision and accuracy determined for analysis of REE in a BHVO-2G glass standard using different fluences. The best precision is seen at > 6J/cm 2 , whereas the best accuracy is seen between 8 and 10J/cm 2 . Grey shaded area shows the conditions that give the highest accuracy and precision.
of REE in the unknown material is over predicted at lower repetition rates. We therefore suggest that a repetition rate of 10-20 Hz should be used for analysis of basaltic glass. Results suggest that below 10 Hz, the accuracy and precision of analysis is severely affected. Grey shaded area represents the conditions which give the best accuracy and precision.