Rapid heterogeneous assembly of multiple magma reservoirs prior to Yellowstone supereruptions

Large-volume caldera-forming eruptions of silicic magmas are an important feature of continental volcanism. The timescales and mechanisms of assembly of the magma reservoirs that feed such eruptions as well as the durations and physical conditions of upper-crustal storage remain highly debated topics in volcanology. Here we explore a comprehensive data set of isotopic (O, Hf) and chemical proxies in precisely U-Pb dated zircon crystals from all caldera-forming eruptions of Yellowstone supervolcano. Analysed zircons record rapid assembly of multiple magma reservoirs by repeated injections of isotopically heterogeneous magma batches and short pre-eruption storage times of 103 to 104 years. Decoupled oxygen-hafnium isotope systematics suggest a complex source for these magmas involving variable amounts of differentiated mantle-derived melt, Archean crust and hydrothermally altered shallow-crustal rocks. These data demonstrate that complex magma reservoirs with multiple sub-chambers are a common feature of rift- and hotspot related supervolcanoes. The short duration of reservoir assembly documents rapid crustal remelting and two to three orders of magnitude higher magma production rates beneath Yellowstone compared to continental arc volcanoes. The short pre-eruption storage times further suggest that the detection of voluminous reservoirs of eruptible magma beneath active supervolcanoes may only be possible prior to an impending eruption.

[2] Analytical techniques For the majority of zircons analysed in this study we use analytical protocols that allow the analysis of oxygen isotopes, trace elements, U-Pb and Hf isotopes on the same crystals employing four different mass spectrometers. The entire sequence on analyses is illustrated in Fig. S1 and the protocols of individual analytical techniques is given in detail below.
[2.1] Oxygen isotope analysis Sample preparation and secondary ion mass spectrometry (SIMS) were carried out at the Canadian Centre for Isotopic Microanalysis (CCIM), University of Alberta. Polished zircon mid-sections of unknowns and zircon reference materials were exposed within a 25 mm diameter epoxy mount (M1174) using diamond grits. The mount was cleaned with a lab soap solution, and de-ionized H 2 O. Prior to scanning electron microscopy (SEM), the mounts was coated with 5 nm of high-purity Au. SEM characterization was carried out with a Zeiss EVO MA15 instrument equipped with a high-sensitivity, broadband cathodoluminescence (CL) detector and a backscattered electron detector. Typical beam conditions were 15 kV and 3-4 nA. A further 25 nm of Au was subsequently deposited on the mount prior to SIMS analysis. CL images of all zircons are given in Fig. S2. Oxygen isotopes ( 18 O, 16 O) in zircon were analyzed using a Cameca IMS 1280 multicollector ion microprobe. A 133 Cs+ primary beam was operated with impact energy of 20 keV and beam current of ~ 3.0 nA. The ~12 µm diameter probe was rastered (15 x 15 µm) for 30 s prior to acquisition, and then 5 x 5 µm during acquisition, forming rectangular analyzed areas ~15 x 18 µm across and ~2 µm deep. The normal incidence electron gun was utilized for charge compensation. Negative secondary ions were extracted through 10 kV into the secondary (Transfer) column. Transfer conditions included a 122 µm entrance slit, a 5 x 5 mm pre-ESA (field) aperture, and 100x sample magnification at the field aperture, transmitting all regions of the sputtered area. No energy filtering was employed. The mass/charge separated oxygen ions were detected simultaneously in Faraday cups L'2 ( 16 O-) and H1 ( 18 O-) at mass resolutions (m/∆m at 10%) of 1900 and 2250, respectively. Secondary ion count rates for 16 O-and 18 O-were typically ~3 x 10 9 and 6 x 10 6 counts/s utilizing 10 10 Ω and 10 11 Ω amplifier circuits, respectively. Faraday cup baselines were measured at the start of the analytical session. A single analysis took 4 minutes, including pre-analysis rastering, automated secondary ion tuning, and 90 s of continuous peak counting. Instrumental mass fractionation (IMF) was monitored by repeated analysis of a zircon primary reference material (RM), S0081 (UAMT1) with δ 18 O VSMOW = +4.87 (R. Stern, unpublished laser fluorination data from Ilya Bindeman, University of Oregon) and a secondary zircon RM, S0022 (Temora-2) zircon with δ 18 O VSMOW = +8.20 ‰ (ref. 49). One analysis of the primary and secondary RM was taken after every 4 and 8 unknowns, respectively. The data set of 18 O-/ 16 O-for S0081 zircon was processed collectively for each of four analytical sessions, yielding a standard deviation of 0.07 -0.10 ‰ after corrrection for systematic within-session drift (<0.3 ‰). Overall correction for IMF was <0.3 ‰. The individual spot uncertainties at 95% confidence for δ 18 O VSMOW reported include errors relating to within-spot counting statistics, between-spot (geometric) effects, and correction for instrumental mass fractionation, and average ±0.25 ‰ (Tab. S1). Accuracy and reproducibility of these protocols was assessed by repeat analyses of secondary reference zircon Temora-2 (S0022); twenty-six analyses, processed as unknowns, yielded a mean δ 18 O VSMOW = +8.260 ± 0.061 (MSWD = 1.6, n = 26; Fig. S3).

[2.2] Trace element analyses of zircons and melt inclusions
Trace elements in zircons were measured by laser ablation inductively coupled plasma mass spectrometry (LA-ICPMS) at the University of Lausanne employing a NewWave UP-193FX ArF excimer laser ablation system coupled to a Thermo Element XR sector-field ICPMS. The laser ablation system was operated with a repetition rate of 5 Hz, a spot diameter of 25 µm and an on-sample energy density of 3 J/cm 2 . The glass standard reference material (SRM) 612 of the National Institute of Science and Technology (NIST) was used for external standardization and all intensities were internally normalized to 29 Si. During data reduction, care was taken to avoid parts of the time-resolved signals that were affected by ablation of mineral or melt inclusions. In some cases not all of these signal parts could be avoided leading to elevated non-stoichiometric bivalent elements as well as elevated LREE (mainly related to feldspar and/or apatite inclusions). Repeat analyses of secondary zircon reference material Mud

[2.3] Uranium-lead geochronology
Following in situ oxygen isotope and trace element analyses, zircons were recovered from the grain mounts and analysed for U-Pb isotopes employing chemical abrasion isotope dilution thermal ionization mass spectrometry (CA-ID-TIMS) techniques at the University of Geneva following established procedures (Refs. 21,22,57). All crystals were annealed by heating to 900°C for 48 h in a muffle furnace prior to in-situ analyses. Individual crystals were subsequently extracted from the grain mounts under a binocular microscope, transferred into 3 ml Savillex beakers, rinsed with 3 N HNO 3 and loaded into 200 µl Savillex microcapsules for partial dissolution (i.e., "chemically abraded"; ref. 51) in HF + trace HNO 3 at 180°C for 15 h in Parr bombs. After partial dissolution, residual crystals/crystal fragments were transferred back into 3 ml Savillex beakers, rinsed with water, fluxed for several hours in 6 N HCl and ultrasonically cleaned in 3 N HNO 3 . Zircons were then reloaded into their respective 200 µl microcapsules, spiked with 3-5 mg of the EARTHTIME 202 Pb-205 Pb-233 U-235 U tracer solution (hereafter referred to as ET2535; tracer calibration v. 3.0; http://www.earth-time.org/; ref. 52) and dissolved in ~70 µl HF at 210°C for >72 h in Parr bombs. After dissolution, samples were dried down and redissolved in 6 N HCl at 180°C overnight, dried down again and redissolved in 3 N HCl. U and Pb were then separated using a modified HCl-based singlecolumn anion exchange chemistry 58 . The U-Pb fractions were dried down with a microdrop of 0.02 M H 3 PO 4 and loaded on outgassed single Re-filaments with a Si-Gel emitter 59 . U and Pb isotopic measurements were performed on a Thermo TRITON thermal ionization mass spectrometer. Pb was measured in dynamic mode on a MasCom secondary electron multiplier. Analyses of ET2535-spiked samples were corrected for instrumental mass fractionation using the fractionation factor derived from the measured 202 Pb/ 205 Pb ratio relative to a true value of 0.99924. U was measured as U-oxide in static mode on Faraday cups equipped with 10 12 Ω resistors. Measured isotopic ratios were corrected for isobaric interferences of 233 U 18 O 16 O on 235 U 16 O 2 using an 18 O/ 16 O of 0.00205 and for mass fractionation using the measured 233 U/ 235 U ratio relative to a value of 0.99506 for both tracers and a sample 238 U/ 235 U of 137.818 ± 0.045 (2σ; ref. 60). Total procedural blanks measured during the same period as the samples averaged 0.17 ± 0.09 pg (n = 17) and the average blank isotopic composition (see Tab. S2) was used for blank correction of zircon analyses. Common Pb in excess of blank was corrected using an estimate computed using the two-stage Pb evolution model of Stacey and Kramers 61 . Note that the computed common Pb composition is well within uncertainty of the measured blank composition. Data reduction was performed using Tripoli and U-Pb_Redux software 53 that employs algorithms of McLean et al. 54 . U-Pb ratios and dates were calculated relative to a 235 U/ 205 Pb ratio of 100.23 ± 0.046 % (2σ) and using the decay constants of ref. 62. All 206 Pb/ 238 U dates were corrected for initial 238 U-230 Th disequilibrium in the 238 U-206 Pb decay chain that arises from preferential exclusion of Th relative to U during zircon crystallization (e.g., ref. 63). We here use a model correction assuming that variations in Th/U of analyzed zircons are due to variations in the Th/U of the magma in equilibrium with the respective zircon during crystallization and not due to variations in relative partitioning between Th and U (i.e., D Th /D U remains constant). We use a partition coefficient ratio D Th /D U derived from published empirical and experminetal data 64-68 and our analyses of melt inclusion-zircon pairs (Fig. S5). The average D Th /D U derived from our melt inlusion-zircon pairs is 0.214 ± 0.094 (Fig. S5). Including the published data yields an average of 0.214 ± 0.074 (Fig. S5), that we use for all analyses of zircons in this study and for recalculation of published data 24,27 . This correction results in a constant increase of 206 Pb/ 238 U dates by 86.5 ± 8.6 ka. Uncertainties associated with this correction were propagated into the uncertainty of 206 Pb/ 238 U dates of individual zircons.
Tank yielded external reproducibilities (2 R.S.D.) between 2.9 % for Hf and 18.1 % for Ti (Tab. S1; Fig. S3). Temora-2 is heterogeneous with respect to most trace elements but the analyses are in good agreement with previously reported data (Ref. 49). Mud Tank is homogeneous but has significantly lower trace element concentrations compared to our samples, suggesting that the reproducibility can be considered a maximum uncertainty for the samples. Th/U ( 238 U and 232 Th) was measured in zircon-hosted melt inclusions of some of the analysed crystals (n=5; Tab. S1; Fig. S5) using the same instruments as for zircons. Due to the small size of the analysed melt inclusions, the ablation system was operated with a pit diameter of 10 µm at a repetition rate of 7 Hz. We employed a fast scanning routine with static magnet to reduce magnet settling times and depth penetration. Mass bias was quantified using repeat analyses of NIST SRM 612. Time-resolved signals were carefully evaluated to avoid signal parts affected by ablation of host-zircon material. Th/U was calculated from 232 Th/ 238 U assuming a 238 U/ 235 U of 137.8.

[2.4] Hafnium isotope analysis
Hafnium isotope analyses were conducted at the University of Geneva employing a THERMO Neptune Plus MC-ICP-MS equipped with nickel cones (Thermo X-series). Analyses were performed under dry plasma conditions employing a Cetac Aridus-2 desolvation unit and a Teflon nebulizer (~50µl /min uptake rate). The 9 Faraday cups were configured to measure the following masses at low mass resolution (m/∆m ~ 450): 172 (L4), 173 (L3), 175 (L2), 176 (L1), 177 (C), 178 (H1), 179 (H2), 180 (H3), 181 (H4). All analytes were diluted in 2 % HNO3 with traces of HF in order to minimize memory effects during analyses and stabilize Hf ions in solution. The certified JMC475 standard was used as a reference throughout the sessions, and was run at a concentration of ~30 ppb. A series of JMC475 solutions doped with variable amounts of Yb with concentrations ranging from 0.5 to 5 ppb (i.e., 0.001< 173 Yb/ 177 Hf<0.5) was prepared in order to optimize the 176 Yb interference correction during postmeasurement data processing. Natural zircon standards Temora-2 50 and Plešovice 69 were measured as secondary reference materials. Two aliquots of Plešovice were also doped with Yb ( 173 Yb/ 177 Hf ~ 0.23). Samples were prepared from the Zr-Hf-trace element fractions of U-Pb column chemistry that were dried down and brought into solution using ~500 µl of the same batch of 2% HNO 3 used for analytes dilution and wash cycles (~400 s between analytes). Measurements included an 80-100 s uptake, baseline control over 30 s and 330 s (80 cycles of 4.2 s) of analysis, for a total consumption of ~450 µl. The βYb and βHf mass bias coefficients were calculated using an exponential law 70 from the measured 172 Yb/ 173 Yb and 179 Hf/ 177 Hf respectively and using natural abundances reference values: 172 Yb/ 173 Yb = 1.3534 and 179 Hf/ 177 Hf = 0.7325 (ref. 71). Mass bias correction on Lu was calculated using βYb and the isobaric interference of 176 Lu was removed using the natural abundances of 175 Lu (0.97416) and 176 Lu (0.02584). The isobaric interference of 176 Yb is the most significant and was evaluated from Yb doped JMC475 solutions by calculating 176 Yb/ 173 Yb using βYb and using a reference value that is empirically determined for each run by minimizing the dependence of corrected 176 Hf/ 177 Hf ratios on the measured 173 Yb/ 177 Yb (Fig. S3;ref. 72). Correction for in situ 176 Hf in-growth due to 176 Lu λ-decay has been calculated using λ 176 Lu of Scherer et al. 73 and the 206 Pb/ 238 U date of the respective zircon. The 176 Hf/ 177 Hf ratio for all JMC475 over the whole campaign of measurements is 0.282152 ± 22 (2 S.D., n=35), corresponding to a reproducibility of ± 0.78 ε-units (Fig. S4). All analyses were normalized to the recommended value of 0.282160 (ref. 74). Reported uncertainties of samples and secondary reference zircons include the within-run precision and the reproducibility of the JMC475 measurements propagated by quadratic addition. Accuracy and reproducibility of these protocols were assessed by repeat analyses of secondary zircon standards Plešovice and Temora-2 (Fig. S4). Analyses of Plešovice zircon solutions yielded a weighted mean εHf t of -3.33 ± 0.26 (2σ; MSWD = 1.3; n = 10; t = 338 Ma) and repeat analyses of Temora-2 yielded a weighted mean εHf t of 5.90 ± 0.20 (2σ; MSWD = 1.2; n = 18; t = 417 Ma) that are both in good agreement with previously reported data 69,75 (Fig. S4). All compositional and isotopic data for samples and reference materials are given in Tab S2.

[2.5] Zircon trace element modelling
We modelled the effect of fractional crystallization on zircon trace element compositions employing a simple mass balance model. The model uses the average compositon of zircon cores as the starting composition and models changing zircon trace element compositions as a consenquence of co-crystallization of a typical Yellowstone mineral assemblage. The model was run as a Markov Chain Monte Carlo simulation that randomly varies the input parameters within the specific limits with 10 4 simulations per step (10% steps; displayed are only median values of modelled populations). The fractionating assemblage of the model includes: Quartz (30-40%), sanidine (20-30%), plagioclase (20-30%), clinopyroxene (5-8%) and accessory zircon (0.01-0.04%), allanite and chevkinite (treated together; 0.02-0.2%). Partition coefficients for all modelled mineral phases were varied widely covering the full range of values of refs. 76-82.
All uncertainties are reported at the 95% confidence level. Single crystal dates are reported and plotted only with their analytical uncertainties excluding systematic contributions. However, for accurate comparison with 40 Ar/ 39 Ar dates, systematic uncertainties associated with the tracer isotopic composition and the 238 U decay constant have to be taken into account. Figure S1. Illustration of our analytical protocol that permits analyses of oxygen isotopes and trace elements in-situ prior to zircon dissolution followed by U-Pb and Hf isotope analyses on the same volume of dissolved zircons after ion exchange chemistry (modified from ref. 83). Oxygen isotope and trace element analyses were performed on the exact same location within a given crystal or within the same textural domains of zircons.      Single-crystal hafnium isotopic data for Yellowstone zircons with reduced chi-squared statistics (MSWD) to assess hafnium isotopic heterogeneities within samples. Also shown are the data for secondary reference zircons for direct comparison. Note the significantly larger scatter in the samples compared to the secondary standards, suggesting real isotopic diversity in the samples. Figure S8. Comparison of temperature estimates for Yellowstone tuffs. Shown are liquidus temperatures and zircon saturation temperatures 55 recalculated using the calibration of Boehnke et al. 84 in comparison to median Ti-in-zircon temperatures of zircon cores and rims. Ti-in-zircon temperatures were calculated using the calibration of Ferry and Watson 85 with a SiO2 =1 and a TiO2 =0.55. Table S1. In-situ oxygen isotopic (SIMS) and trace element (LA-ICPMS) data for analysed Yellowstone zircons, secondary zircon reference materials, and Th/U of zircon-hosted melt inclusions with calculated D Th/U .         Liq. T Zrc Ti T, core   b Corrected for initial Th/U disequilibrium using the average Th/U partition coefficient ratio derived from melt inclusion analyses and published data ( Figure S5). d Th contents calculated from radiogenic 208 Pb/ 206 Pb assuming concordance between U-Th and Pb systems. e Total mass of radiogenic Pb. f Total mass of common Pb. g Ratio of radiogenic Pb (including 208 Pb) to common Pb. h Th/U ratio of magma from which mineral crystallized. i Measured ratio corrected for fractionation and spike contribution only. j Measured ratios corrected for fractionation, tracer, blank and initial Pb. Blank was corrected using the average composition of total procedural blank (TPB; average=0.17 ± 0.09 pg; n=17) measurements peformed throughout the course of this study: 206Pb/204Pb=18.102 ± 0.378, 207Pb/204Pb=15.352 ± 0.469; 208Pb/204Pb=38.01 ± 1.14 (all uncertainties are 1 S.D.).