High-precision zircon age spectra record the dynamics and evolution of large open-system silicic magma reservoirs

The emplacement history and thermal evolution of subvolcanic magma reservoirs determine their longevity, size, and ability to feed volcanic eruptions. As zircon saturation is dependent on melt temperature and composition, quantitative analysis of zircon age distributions provides insight into the timing and magnitude of intensive parameter variations during the lifetime of a magma reservoir. Here we present chemical abrasion-isotope dilution-thermal ionization mass spectrometry (CA-ID-TIMS) U-Pb zircon crystallization ages and in-situ zircon trace-element geochemistry for a suite of caldera-related, coeval plutonic and volcanic units of the Permian Sesia Magmatic System (northern Italy). This dataset documents the protracted growth and evolution (~1 Myr) of a voluminous ( > 1000 km 3 ), upper crustal (3.5 – 2.0 kbar) silicic magma reservoir. Systematic changes in zircon composition with time reveal episodic intrusions of new magma into a single, progressively differentiating reservoir which was dominantly kept at high crystallinity conditions ( > 60 vol.%). The volcanic and plutonic units both show dispersed and heterogeneous CA-ID-TIMS age distributions. A stochastic, thermodynamics-based zircon saturation model accounting for fractional crystallization, recharge and thermal rejuvenation effects on zircon growth and stability in a rhyolitic magma body reproduces the observed distributions of zircon U-Pb ages. Model results suggest that zircon compositional variability and crystallization age heterogeneity, characteristic of the Sesia and other large silicic systems, distinguish open-and closed-system magmatic processes. Specifically, an increase of the mass of crystallized zircons over time is the hallmark of thermally and chemically mature magma reservoirs, capable of producing voluminous eruptions upon changes in the thermal and mechanical state of the system. Conversely, a decrease in crystallized zircon mass over time characterizes systems where crystallization dominates over magma mobility, hindering substantial mass release during eruptive events. We suggest that quantitative analyses of dispersed zircon crystallization age distributions obtained with high-precision techniques can be applied to the plutonic and volcanic record to identify mature silicic magma reservoirs whose properties and storage conditions allowed for catastrophic caldera-forming eruptions. This work provides the petrologic community with a new tool to investigate the genesis and evolution of zircon-bearing intermediate to silicic intrusive magmatism and caldera-forming volcanism through time and across tectonic settings.


Introduction
Time-resolved records of the thermal and compositional evolution of shallow (<3 kbar) crustal magma bodies provide insights into the processes that regulate tempo, volume, style of volcanic eruptions and the origin of silicic plutons.Zircon, a common accessory mineral in silicic igneous rocks, represents the prime tracer of physico-chemical processes unfolding during the lifetime of a subvolcanic magma reservoir, owing to its unparalleled ability as a timekeeper (Schoene and Baxter, 2017) and the slow diffusivity of most elements through its crystal lattice (Cherniak and Watson, 2003).Recent advances in zircon geochronology by chemical abrasion-isotope dilution-thermal ionization mass spectrometry (CA-ID-TIMS) allow for internal age precision better than the 0.02% level for single zircon 206 Pb/ 238 U dates (Wotzlaw et al., 2017;Szymanowski and Schoene, 2020).The integration of CA-ID-TIMS U-Pb geochronology, trace element and isotopic analyses in zircons (e.g.Schoene et al., 2012;Wotzlaw et al., 2013;Szymanowski et al., 2017) with advances in thermomechanical modeling (e.g.Annen et al. 2009;Gelman et al., 2013;Karakas et al., 2017;Huber et al., 2019) have provided invaluable insights into the nature and duration of upper-crustal magmatism.These works have contributed to the paradigm of long-lived (10 5 -10 6 yr), largely crystalline magma reservoirs that are periodically recharged by new magma inputs, as the dominant model of shallow, silicic magma storage in the Earth's crust (Sparks et al., 2019).
In U-Pb geochronological studies of upper-crustal magmatic reservoirs, the internal age resolution afforded by CA-ID-TIMS can resolve heterogeneities in zircon crystallization dates, i.e. dispersed age populations (Keller et al., 2018).Such age dispersion is interpreted as either indicative of protracted saturation and growth of autocrystic zircons (Wotzlaw et al., 2013;Samperton et al., 2015;Ratschbacher et al., 2018;Szymanowski et al., 2019;Klein and Eddy, 2023) or as a record of recycled and juxtaposed zircons crystallized in a polybaric storage system (Schoene et al., 2012;Barboni et al., 2015;Farina et al., 2018).In the case of protracted saturation, as growth and survival of autocrystic zircons is linked to the physical and chemical conditions of magma reservoirs (Watson, 1996;Boehnke et al., 2013;Bindeman and Melnik, 2016), their distribution in time has been proposed to provide quantifiable constraints on the compositional and thermal states of crystallizing magma bodies (Caricchi et al., 2014;Tierney et al., 2016;Samperton et al., 2017;Kent and Cooper, 2018;Kirkland et al., 2021) and their eruptive potential (Weber et al., 2020).Specifically, the simple distribution of zircon mass crystallized as a function of time and temperature observed in monotonically cooled plutonic systems (Samperton et al., 2017;Kirkland et al., 2021) can be reproduced by experimentally derived zircon solubility models (Watson and Harrison, 1983;Boehnke et al., 2013) and thermodynamic assessment of phase equilibria (Ghiorso et al., 2002;Holland and Powell, 2011).However, no study has yet linked the complex, heterogeneously dispersed CA-ID-TIMS zircon age datasets to open-system, crystal-melt dynamics evident from petrological studies of silicic systems (Szymanowski et al., 2020) and suggested by volcanic and plutonic zircon age spectra displaying crystallization gaps and apparent late peaks of zircon crystallization (e.g.Wotzlaw et al., 2013;Rioux et al., 2016).
To address this issue, we integrate CA-ID-TIMS U-Pb zircon geochronology and in-situ geochemistry on co-genetic, silicic intrusive and eruptive products of a Permian magmatic plumbing system (Sesia Magmatic System; Quick et al., 2009;Karakas et al., 2019) with a novel approach to zircon saturation modeling that combines mass-and energy-constrained thermodynamic tools with stochastic sampling (Keller et al., 2017;Bohrson et al., 2020;Tavazzani et al., 2023).We demonstrate that the compositional variability and heterogeneous age distribution of autocrystic zircons from an internally complex subvolcanic intrusion and coeval voluminous rhyolitic eruptions are faithful chroniclers of incremental growth, differentiation and episodic recharge in a large upper-crustal magma reservoir.Based on this successful application, quantitative evaluation of zircon age dispersion is proposed as a tool to push the analyses of open-system thermomechanical processes into deep time plutonic and volcanic records.

Geology of the Sesia magmatic system
The Sesia Magmatic System (SMS) is a transcrustal igneous system that developed in the pre-Alpine basement of the Southern Alps during early Permian extension.Due to Mesozoic-Cenozoic exhumation, the SMS is now exposed in tilted crustal blocks in the Italian Alps as part of the Ivrea-Verbano and Serie dei Laghi terranes (Fig. 1a;Fountain, 1976;Handy and Zingg, 1991).During its four-million-year lifetime (~285 to ~281 Ma; Karakas et al., 2019) the SMS produced synchronous magmatism across a thick crustal column, ranging from a stratified gabbroic intrusion in the lower crust (Rivalenti et al., 1981) to explosive volcanism at the surface (Quick et al., 2009).At the apex of its activity, the SMS generated the Sesia Caldera (>13 km in diameter) that erupted at least 400 km 3 of dominantly silicic materials (Quick et al., 2009).
The main eruptive unit of the Sesia Caldera (SC) is a densely welded, crystal-rich (25-30 vol.%) rhyolitic ignimbrite that includes patches of megabreccia (Fig. 1b).Textural and compositional features of quartz phenocrysts indicate that the eruption of voluminous crystal-rich material was triggered by thermal rejuvenation (i.e.reduction in crystallinity; Huber et al., 2011) of an evolved magma body, stored at near-solidus, high crystallinity conditions and shallow depths (2.0 to 2.5 kbar; Tavazzani et al., 2020).A layer of volcaniclastic sediments (~30 m thick) separates the main eruptive unit of the SC from a ~400 m thick, welded rhyolitic ignimbrite unit (Sbisà, 2009), characterized by high-silica content (75-77 wt.% SiO 2 ), diffuse sub-parallel eutaxitic texture and low crystallinity (~10 vol%).This laterally continuous but volumetrically limited (1-10 km 3 ) rhyolitic unit is interpreted as the product of a post-caldera eruption due to its stratigraphic position, and it is referred in this work as the upper ignimbrite unit.
Directly underneath the SC deposits, the upper-crustal segment of the SMS consists of three granitoid plutons emplaced in the Serie dei Laghi metamorphic sequence (Handy and Zingg, 1991;Brack et al., 2010).The southernmost of these bodies is the Valle Mosso Pluton (VMP), a composite silicic intrusion (~260 km 3 ; Fig. 1b) ranging from quartz-monzonite to high-silica leucogranite (65-77 wt% SiO 2 ).The bulk of the intrusion consists of variably porphyritic monzogranite which contains petrological evidence of recharge by less evolved melts and fluctuating thermal conditions (e.g., disequilibrium textures in low temperature phases; Tavazzani et al., 2017Tavazzani et al., , 2020)).Overall, the architecture, textural and bulk-rock compositional features of the VMP (Fig. 1b) reflect the evolution of a silicic magma reservoir undergoing protracted differentiation, driven by crystal-melt separation, and episodically recharged by granitic and, to a lesser extent, mafic magma (Sinigoi et al., 2016;Tavazzani et al., 2020).
The spatial proximity, overlapping compositional range and shared thermal history support a genetic link between the VMP intrusion and the eruptive products of the SC.Recent high-precision U-Pb zircon age determinations on VMP and SC samples (Karakas et al., 2019) further reinforce this relationship, providing an ideal opportunity to quantify timescales of assembly, storage and differentiation in an interconnected, upper crustal magmatic reservoir in relation to a voluminous caldera-forming eruption.

Materials and methods
Informed by recent, detailed mapping and petrologic investigation (Tavazzani et al., 2017(Tavazzani et al., , 2020)), we selected six samples for high-precision CA-ID-TIMS U-Pb zircon dating (Fig. 1b) that cover the full extent of textural and geochemical variability of the upper-crustal intrusive and eruptive products of the SMS.We compile 39 new dates with recently published zircon ages of Karakas et al. (2019) into a dataset of 96 single CA-ID-TIMS zircon U-Pb dates for upper-crustal lithologies of the SMS.

LA-ICP-MS U-Pb zircon geochronology and trace element analyses
Zircon grains were separated from thirteen crushed rock specimens using standard heavy mineral separation techniques.Selected zircons were hand-picked and annealed by heating to 900 • C for 48 h in a muffle furnace before mounting in epoxy mounts and polishing to expose central cross sections.Zircons were then texturally characterized employing cathodoluminescence (CL) imaging (Supplementary Figures S2-S4).Following imaging, zircon grains were analyzed by laser ablation-inductively coupled plasma mass spectrometry (LA-ICPMS), to determine trace element abundances and U-Pb ages from the same ablation crater, employing a 193 nm Resonetics Resolution S155 laser ablation system coupled to a Thermo Element XR, sector-field single collector ICP-MS (Guillong et al., 2014).Detailed procedures and analytical parameters are reported in Supplementary Text and Supplementary Data Table S2, following community-derived guidelines (Horstwood et al., 2016).

CA-ID-TIMS zircon U-Pb geochronology
Based on CL-imaging and LA-ICPMS results, fifty-six magmatic zircons from six samples, showing no evidence of inheritance or contamination by inclusions were selected for CA-ID-TIMS analyses.Individual zircons were extracted from the epoxy mounts using stainless steel tools and chemically abraded for 15 h at 190 • C (Mattinson, 2005), cleaned in 6 M HCl and 3.5 M HNO 3 , loaded into 200 μl Savillex microcapsules with concentrated HF and spiked with ~5-8 mg of the EARTHTIME 202 Pb-205 Pb-233 U-235 U solution (Condon et al., 2015;McLean et al., 2015).Samples were dissolved in Parr bombs at 220 • C for ~60 h.After dissolution, samples were dried down and redissolved in 6 M HCl at 190 • C for several hours and then dried down again before redissolving them in 3 M HCl for ion exchange chemistry.U and Pb were separated using a single-column HCl-based ion exchange chemistry (Krogh, 1973).
The U-Pb fractions were loaded on outgassed Re filaments with ~1 μl of Si-gel emitter (modified from Gerstenberger and Haase, 1997).U and Pb isotopes were measured using a Thermo TRITON Plus thermal ionization mass spectrometer at the Institute of Geochemistry and Petrology of ETH Zurich (Wotzlaw et al., 2017).Data reduction, age calculation and uncertainty propagation were performed using the Tripoli and ET_redux software package (Bowring et al., 2011;McLean et al., 2011).Detailed analytical procedures are reported in Supplementary Text.

Zircon saturation model
In parallel with the CA-ID-TIMS analytical effort, we developed a model to quantify and track changes in saturated zircon mass during protracted crystallization, differentiation and repeated intrusions of magmas in an upper crustal rhyolitic reservoir.In contrast to previous studies (e.g.Tierney et al., 2016;Kent and Cooper, 2018;Weber et al., 2020), this model does not rely on a parameterized relationship between temperature and fraction of zircon crystallizing (Harrison et al., 2007) but integrates temperature, crystallinity, melt major and trace element evolution to assess zircon saturation (Fig. 2).The trajectory of intensive variables during crystallization is derived using the software Magma Chamber Simulator (Bohrson et al., 2014;Bohrson et al., 2020), a modeling tool combining thermodynamic assessment of phase equilibria (based on the MELTS family of algorithms of Ghiorso et al., 2002;Gualda et al., 2012) with energy and mass conservation.
The mass of zirconium (Zr) available for crystallization (m Zr , Fig. 2a) is calculated at every step of a crystallization simulation applying the workflow developed by Keller et al. (2017): where m liq is the mass of liquid phase, [Zr] melt calc is the Zr concentration in the liquid phase and [Zr] melt sat is the Zr concentration required to saturate the melt in zircon using the calibration of Boehnke et al. (2013).At each step (n) of the simulation, relative change in zircon mass (dm Zrn = m Zrn − m Zr (n− 1) ) is used to derive zircon crystallization probability (p Zrn , Fig. 2b) which is normalized to the overall variation in crystallized zircon mass ( ∑ N n=0 dm Zrn , N represents the total number of steps in a crystallization  (3) In order to simulate the stochastic process of zircon selection during geochronological analysis, a variable number of zircons (N Zr , user defined) is randomly sampled from the cumulative distribution function F pZr (n) of the target probability density function f pZr (n) through a Monte Carlo based model (Fig. 2c).The Monte Carlo algorithm undergoes 10 3 iterations for each simulation to ensure robustness of the random selection process.Bootstrapped zircon age spectra, where age is expressed as relative zircon crystallization density functions f xtal (t r ) in the time interval (t r ) from crystallization of the first zircon (t sat ), are generated from outputs of the Monte Carlo simulation and visualized through a truncated kernel density estimate of the scaled crystallization ages Δt h is the difference between the oldest and youngest zircon dates in a sample; Δt a is the difference between oldest and youngest autocrystic zircon dates.Eruption ages are defined for each volcanic sample as the Bayesian estimate (Keller et al., 2018) using the bootstrapped distribution of Th-corrected 206 Pb/ 238 U dates as prior, with 2σ uncertainty given as internal only/with tracer calibration/ with tracer and 238 U decay constant.Individual zircon dates for crystal-rich rhyolite sample MBR and range of zircon dates from 4 granitoids of the Serie dei Laghi are reproduced from Karakas et al. (2019).(b, c) Temporal distribution of individual zircon ages reported as kernel estimates of relative zircon crystallization density functions f xtal (t r ), where t r is relative crystallization time, normalized between x = 1, denoting zircon saturation (t sat ) and x = 0, corresponding to full plutonic solidification (t sol ) or eruption (t er ).The two kernel density estimates represent zircon crystallization distributions bootstrapped over all dated samples in the eruptive (b) and intrusive (c) domains of the SMS.
Taking advantage of this flexible setup, we simulated contrasting thermal and physical trajectories for a single rhyolitic magma batch stored at upper crustal depth, namely: (1) equilibrium crystallization, (2) fractional crystallization, or (3) intrusion during continuous differentiation (via fractional crystallization) by one to five recharges of variable volume and composition (i.e.rhyolitic or basaltic).Simulations are run until the magma body reaches solidus temperature or alternatively interrupted at 750 • C (c. 50 vol.%melt), to simulate the abrupt termination of crystallization caused by an eruptive event.The wholerock composition of a monzogranite from the VMP is used as the initial liquid composition, initial magmatic water content is set to 2.5 wt % H 2 O and fO 2 to the QFM (quartz-fayalite-magnetite) buffer (Tavazzani et al., 2020).The complete list of input parameters, explored crystallization simulations, model limitations and applicability range are discussed in Tavazzani et al. (2023).The full set of open-source MATLAB® scripts necessary to run zircon saturation simulations is available at https://github.com/TavazzaniL/OpenSystemZrSat).

CA-ID-TIMS geochronology
Out of fifty-six analyses, thirty-nine Th-disequilibrium corrected 206 Pb/ 238 U dates are interpreted as magmatic crystallization ages whereas the remaining dates are affected either by the presence of xenocrystic components or post-crystallization loss of radiogenic lead (discussion of rejected fractions is available in the Supplementary Text).This dataset is augmented with published zircon data of a Sesia Ignimbrite sample from Karakas et al. (2019).All samples display resolvable age dispersion (Fig. 3a), aside from a small microgranitic porphyry body (Gp) in the VMP in which textural and petrologic evidence suggest rapid cooling from near-liquidus conditions (Tavazzani et al., 2020).The range of individual zircon dates overlaps in an interval between 283.984 ± 0.057 and 282.97 ± 0.17 Ma (2σ uncertainty, Fig. 3a) with the exception of an older (up to 284.8 Ma) component represented by several analyses from the crystal-rich rhyolite (Sei) and in two zircons from the VMP homogeneous monzogranite samples (Hwg).The zircon fraction younger than 284 Ma shows dispersion in single crystal 206 Pb/ 238 U dates of 1.012 ± 0.178 Myr for the VMP, comparable to observations by Karakas et al. (2019) on other granitoid bodies of the SMS (Fig. 3a).In the volcanic samples, excluding the older component (>284 Ma), U-Pb dates record 0.309 ± 0.080 Myr of zircon crystallization in the crystal-poor rhyolite unit and 0.477 ± 0.188 Myr of crystallization in the crystal-rich Sesia Ignimbrite sample of Karakas et al. (2019), with crystallization intervals for each sample calculated with the iterative approach of Eddy et al. (2016; see Supplementary Text for details).The Sesia and Upper Ignimbrite eruption ages can be obtained as either the Th-corrected 206 Pb/ 238 U ages of the single youngest zircon of each eruptive unit (Sesia Ignimbrite: 283.46 ± 0.14/0.16/0.35Ma; Upper Ignimbrite: 283.22 ± 0.07/0.11/0.33Ma) or as the Bayesian estimate based on the bootstrapped zircon distribution (Keller et al., 2018) that yields 283.49± 0.10/0.13/0.34Ma (Sei) and 283.21 ± 0.10/0.13/0.34Ma (Ui, 2σ uncertainty given as internal only/with tracer calibration/with tracer and 238 U decay constant; Fig. 3a).The estimates obtained by these two approaches are indistinguishable and resolve a ~300 kyr temporal gap between the two eruptive events, consistent with field observations (erosional surface and volcanoclastic interval; Sbisà, 2009).
The two kernel density estimates obtained by bootstrapping zircon CA-ID-TIMS age distributions of all samples from volcanic (Fig. 3b) and plutonic (Fig. 3c) units between zircon saturation (t sat ) and full plutonic solidification (t sol ) or eruption (t er ) are similar and show an increase in relative number of crystallized zircons with time.The distributions are asymmetrical, with zircon crystallization peaks shortly preceding eruption in the volcanic domain and full plutonic crystallization in the intrusive units.The peaks are more evident when including the age fraction >284 Ma (defined as antecrysts in Fig. 3b, c) into the kernel density estimate.

Zircon trace element composition
The dated zircons display a geochemical variability that coincides with the compositional range observed by Karakas et al. (2019) for zircons from upper crustal rocks of the SMS (Fig. 4).Overlapping zircon compositions from volcanic and plutonic units align along a trend consistent with feldspar fractionation (decreasing Eu anomaly, Eu/Eu* and increasing U concentration; Fig. 4a, c) and are strongly correlated with the temperature-modulated, Ti content (Fig. 4b; Ferry and Watson, 2007).The investigated zircons also display a clear negative correlation between Ti, REE contents and heavy to middle REE ratios (i.e.Yb/Dy; Fig. 4d, e).This relationship suggests that heavy to middle REE budgets in SMS melts are controlled by zircon and apatite crystallization without a significant role played by titanite (cf.Wotzlaw et al., 2013;Samperton et al., 2015;Szymanowski et al., 2017;Fig. 4d), consistent with the observed mineral assemblage (Tavazzani et al., 2020).Few zircons display resolvable intra-grain compositional variability (i.e.core to rim; Supplementary Figure S17).
All zircons dated by bulk-grain CA-ID-TIMS reveal a systematic evolution of trace element compositions over the ~1.5 Myr long range of U-Pb zircon crystallization ages (Fig. 5; for a complete discussion of the biases associated with matching in-situ compositions with bulk-grain U-Pb ages see Supplementary Text).The zircon fraction with U-Pb ages >284 Ma shows the least evolved compositions (e.g.Hf<10 4 ppm, Eu/ Eu* of 0.2-0.4,Ti between 5 and 15 ppm; Fig. 5a, d, e) observed in the VMP and SC zircon record.Starting at ~284 Ma, zircons display an increase in U, Hf and MREE concentrations paired with decreasing Ti and Eu/Eu*.During the remaining lifetime of the reservoir (~1 Myr), zircons are characterized by scatter around low Ti (<5 ppm; Fig. 5e), evolved compositions (e.g.Hf>10 4 ppm, Dy ~200 ppm and Eu/Eu* of ~0.1; Fig. 5a, c, d).

Zircon saturation modeling
Three end-member scenarios explore systematic variations in saturated zircon mass during the crystallization of a rhyolitic magma batch, as modulated by the dominant crystallization processes (i.e.equilibrium or fractional crystallization), magma recharge frequency, volume and composition.Simulation results are compared to SMS volcanic and plutonic zircon age distributions (Fig. 6).
The first scenario, shown in Fig. 6a and 6d, simulates the temporal distribution of zircons crystallized in a single magma batch during a fractional crystallization process.In this scenario, after saturation is reached at ~775 • C, zircons precipitate from progressively cooler residual melt until the residual melt fraction is erupted at ~750 • C (Fig. 6a) or until the solidus temperature is reached (Fig. 6d).In saturation calculations truncated at eruption temperatures, early uniform, moderate crystallization rates persist after saturation and decrease sharply just before the eruptive event (Fig. 6a).At plutonic conditions, outputs display high crystallization rates close to zircon saturation temperature, followed by a steady decline approaching the solidus (Fig. 6d).Generally, in volcanic simulations, zircon distributions lack the pronounced asymmetry observed in plutonic simulations.
The changes in temporal distribution of zircon mass induced by magma injection in a crystallizing magma body are explored in two distinct simulations.In the first recharge scenario, three basaltic andesite magma batches are injected in a rhyolitic magma body at nearsolidus, high-crystallinity conditions (i.e.730 • C and c. 75 vol.%crystals) during monotonic cooling and differentiation.Each recharge has the same size (i.e.1/10 of the initial magma mass) and it is added to the crystallizing rhyolite body when residual melt fraction drops below 25 vol.% (Fig. 6b, e).In this scenario, the negative inflection in rate of zircon crystallization is an expression of the transient zircon undersaturation (sensu Bindeman and Melnik, 2016) triggered by stark compositional contrast between host (i.e.rhyolitic) and recharging (i.e.basaltic andesite) magmas.Increase in zircon growth rates is observed only after the last mafic recharge enters the system (Fig. 6b, e).As a result, zircon distributions resulting from these simulations are asymmetrical, with the majority of zircons precipitated from residual melt shortly before solidus is reached or the residual melt fraction is erupted.
The last model setup (Fig. 6c, f) simulates zircon crystallization in a rhyolitic magma body, replenished at near-solidus, high-crystallinity conditions (i.e.730 • C and c. 75 vol.%crystals) by magma batches compositionally identical to the starting material.In this scenario, three large (i.e.1/2 of the initial magma mass), equally sized rhyolitic recharges are introduced in a fractionating magma batch to simulate the incremental growth of a granitic, upper-crustal body.The impact that a sequence of large, rhyolitic recharges has on the temporal distribution of precipitated zircon mass is overall comparable to the effect caused by small mafic injections.However, due to the limited thermal and compositional contrast between newly injected and crystallizing magma, zircon-saturated conditions are maintained throughout the whole simulation.Zircon crystallization rates increase steadily until they reach the highest values near eruption (Fig. 6c) or full plutonic crystallization (Fig. 6f).

Time-resolved record of inception, maturation and eruption(s) of an upper-crustal magma reservoir
The overlapping crystallization ages (Fig. 3a), matching temporal distributions (Fig. 3b, c) and coherent geochemical trends (Figs. 4 and 5) of zircons sampled across volcanic and intrusive lithologies of the SMS suggest zircon saturation and precipitation from melts stored in a common upper crustal reservoir over ~1 Myr.In this context, we rely on precise U-Pb ages (2σ uncertainties on individual dates of 40-190 kyr; Fig. 3) and paired in-situ geochemical analyses to track temporal variations in melt temperature, composition and degree of crystallization over the lifespan of the reservoir (e.g.Schoene et al., 2012;Samperton et al., 2015;Melnik and Bindeman, 2018;Large et al., 2021).
Early crystallized zircons (i.e.>284.0Ma; Fig. 3a) display low incompatible element concentrations (i.e.Hf, Dy), high Ti and Eu/Eu*.The composition of these zircons is comparable to that of zircons crystallized in lower to middle crustal dioritic and gabbroic magmas of the SMS (Karakas et al., 2019;Fig. 5).We interpret these zircons as antecrysts, crystallized during the initial stage of magmatism in the SMS at lower to middle crustal depths from melts of mafic to intermediate compositions.Ascent, progressive cooling and feldspar dominated  2007) using an a TiO2 of 0.6 and a SiO2 of 1.0 (see Tavazzani et al., 2020 for details on activity coefficients determination).In (c), (f) fractional crystallization models are superimposed on trace elements co-variation diagrams.The median composition of zircons with the highest Ti-in-zircon crystallization temperature is selected as initial composition.Percent crystallinity (20%, 50%, etc.) is indicated along the fractionation trajectories.Partition coefficients used and details of each crystallization models are provided in the Supplementary Material.The two models in (c) use either constant partition coefficients ("Constant D") or partition coefficients that vary with temperature ("Variable D") during the course of crystallization.Crystallizing mineral assemblage constitutes of sanidine, plagioclase and quartz in near-eutectic proportions (Luth et al., 1964).In (f), REE evolution scenarios are evaluated as a function of the amount of zircon (zrn), apatite (ap) and titanite (ttn) in the crystallizing mineral assemblage.
fractionation of these parental melts, together with variable incorporation of crustal melting-derived melts, is likely the source of upper crustal granitoid-and rhyolite-hosted zircons of the SMS (Fig. 4; see also Sinigoi et al., 2016;Storck et al., 2021).During a prolonged interval (ca.283.9 to 283.2 Ma) of continuous or punctuated zircon crystallization, pronounced zircon Eu anomaly and high HREE (e.g.Dy) contents indicate the existence of a highly fractionated magma body, likely stored at high crystallinity conditions (estimated >60 vol.% crystals based on simple geochemical modeling; Fig. 4c, d) and invariant, nearly eutectic temperatures (Ti: 5 ± 1 ppm corresponding to 728 ± 20 • C at a SiO2 = 1.0 and a TiO2 = 0.6; Ferry and Watson, 2007).Acknowledging the limits of fractional crystallization modeling and the possibility of higher melt fractions at lower temperatures (Jackson et al., 2018), a similar record might be produced by a less crystallized reservoir whereby the compositional homogeneity is caused by frequent overturns or stirring of the magma body (e.g.Trubač et al., 2017).Regardless of the details of physical state of the upper crustal reservoir, the homogeneous, differentiated melts that characterize the system after 284 Ma (Fig. 5) indicate the establishment of a large, unified magma reservoir or an interconnected network of reservoirs in the upper crust underneath the SC, similar to that inferred for plumbing systems at a number of young calderas (e.g.Yellowstone, Huang et al., 2015;Wotzlaw et al., 2015;Toba, Jaxybulatov et al., 2014).In the case a high crystallinity reservoir, deviations from the trend, in the form of individual zircon analyses or groups of grains showing increased Eu/Eu* (Fig. 5d), temperature-dependent Ti (Fig. 5e) and reduced HREE contents (e.g.Dy; Fig. 5c) can be interpreted as signals of new magma influxes in the reservoir, adding heat and mass to the system with direct effect on melt temperature and mineral phase stability (e.g.Szymanowski et al., 2019).Alternatively, in a scenario where higher melt fraction can be sustained at lower temperatures and the reservoir is well stirred, the variability in composition could reflect the sampling of peripheral part of the system during overturning or rejuvenation events.
Contrasting compositions and age distributions of zircon grains extracted from the two SC ignimbrite units analyzed in Karakas et al. (2019) and in this work (i.e.Sesia Ignimbrite and Upper Ignimbrite) suggest different processes involved in melt mobilization and eruption.The long crystallization interval (~500 kyr, excluding antecrystic components; Fig. 3a) and broad geochemical variability of zircons from the crystal-rich Sesia Ignimbrite (Fig. 4 and Fig. 5), which includes compositions corresponding to both highly differentiated and relatively unevolved melts (Fig. 5d, e), can be explained by thermal rejuvenation of an internally heterogeneous magmatic reservoir (Huber et al., 2011) in the build-up to the SC climactic eruption (>400 km 3 ; Fig. 1b).This process was potentially accompanied by mechanical mixing of distinct zircon populations entrained by peripheral, variably crystallized portions of the system (e.g.Szymanowski et al., 2019;Curry et al., 2021).The narrow crystallization interval and the geochemical homogeneity of zircons from the volumetrically limited, high-silica rhyolite unit (Upper Ignimbrite: 1 to 10 km 3 ; Fig. 1b) suggest a different eruption priming mechanism.In this unit, zircon crystallization ages span <300 kyr and display compositions compatible with crystallization from a highly differentiated melt (Fig. 4 and Fig. 5).Evolved melt domains, isolated in a crystalline framework, are likely buffered at the haplogranitic minimum (Schaen et al., 2018), thus explaining the limited compositional variability of zircons from this unit (at least in compatible elements; Fig. 5d, e).Moreover, the extraction and ascent of these highly differentiated, residual melts through a viscous crystal mush (Bachmann and Huber, 2019), rather than destabilization and homogenization of the entire mush column, provide a mobilization mechanism that accounts for the limited temporal variability of zircons sampled in the Upper Ignimbrite unit.
This set of observations combines with existing textural (e.g.cryptic contacts) and bulk-rock geochemical data (Tavazzani et al., 2020) to convey the image of a large, mostly crystalline and periodically recharged upper crustal reservoir, containing isolated and/or Vertical lines indicate eruption ages defined for each volcanic sample as the Bayesian estimate of the zircon Th-corrected 206 Pb/ 238 U dates bootstrapped distribution (Keller et al., 2018) with 2σ uncertainty given as internal only/with tracer calibration/with tracer and 238 U decay constant.Ti-in-zircon crystallization temperatures estimates use the model of Ferry and Watson (2007) and activity values from Tavazzani et al. (2020).
interconnected melt domains capable of zircon precipitation, which persisted in the upper-crust for ~1 Myr and produced at least one catastrophic caldera forming eruption.This picture supports models emerging from a number of subvolcanic systems tapped by voluminous rhyolitic eruptions (e.g.Wotzlaw et al., 2013;Rioux et al., 2016;Szymanowski et al., 2017) or porphyry intrusions (e.g.Chelle-Michou et al., 2014;Buret et al., 2016) and agrees with timescales derived from thermal models (Gelman et al., 2013;Karakas et al., 2017), which imply storage at dominantly high-crystallinity conditions of episodically replenished silicic magma reservoirs for 10 5 to 10 6 years in the upper crust.

Zircon age dispersion in the SMS
The theoretical trajectory of zircon growth during monotonic cooling and crystallization of an isolated magma batch has been codified in a seminal work by Watson (1996) and recently reproduced integrating thermodynamically derived phase proportions, residual melt compositional evolution and empirical zircon saturation equations (Keller et al.,  (c, f) Recharge and fractional crystallization (RFC) simulation using a rhyolitic recharge magma in three consecutive injection episodes.Zircon age distributions, expressed as relative zircon crystallization density functions f xtal (t r ), where t r is relative time, scaled from zircon saturation (t sat ) to full plutonic solidification (t sol ) or eruptive truncation (t er ), are shown as kernel density estimate truncated at ± 0.5 kernel bandwidth (Keller et al., 2018).Each synthetic zircon distribution (i.e.bold black lines) is the product of random sampling of 30 zircons from synthetic FC and RFC distributions averaged over 10 3 stochastic iterations.Temperature/crystallinity evolution and increments of magma mass over scaled crystallization time are shown in the top panels of each figure.Crystallinity and solidus temperature are based on rhyolite-MELTS (Gualda et al., 2012) derived temperature crystal-fraction curve (i.e.solidus-liquidus relation) for starting granitic composition used in the simulation (Tavazzani et al., 2023).2017; Kirkland et al., 2021;Klein and Eddy, 2023;Tavazzani et al., 2023).Results of these studies indicate that in a monotonically cooled magma batch both equilibrium and fractional crystallization pathways produce non-linear zircon growth.Specifically, a rapid increase in precipitated zircon mass upon saturation in zirconium is followed by a steady decline in zircon crystallization rates until solidus is reached (Fig. 2d; Fig. 6d), with zircon crystallizing from progressively cooler residual melt.These predictions are backed by evidence from silicic intrusions characterized by comprehensive high-precision U-Pb zircon datasets, which display decreasing mass of crystallized zircon with time and linear decrease in temperature (e.g.Samperton et al., 2017;Large et al., 2021) In intrusive and volcanic units of the SMS, the relative distribution of zircon crystallization ages (i.e. increase in relative number of crystallized zircons with time; Fig. 3b, c) is incompatible with zircon growth from a monotonically cooled, isolated magma batch (Fig. 6).Additionally, the SMS zircon geochemical record shows non-linear variations in temperature-modulated Ti contents (Fig. 5e), not consistent with crystallization in a monotonically cooled system.A good match between the shape of zircon age spectra from VMP and SC units and modeled zircon growth is obtained by simulating repeated recharges during fractionation (RFC; Fig. 6c, f), either by small basaltic andesite melt batches (Fig. 6b, e) or large rhyolite injections (Fig. 6c, f).The equilibration of newly injected magma with the resident crystallizing magma causes fluctuations in zircon saturation parameters (e.g.temperature, Zr concentration, degree of melt polymerization) and is reflected in non-linear zircon crystallization probabilities (p Zrn ; Fig. 2b).Moreover, synthetic zircon populations from RFC simulations display a broad spread and systematic changes in crystallization temperatures that match Ti-inzircon temperature evolution observed across VMP and SC zircons (see Tavazzani et al., 2023), when applying activity parameters appropriate for the SMS (i.e. a TiO2 =0.6; a SiO2 =1.0 from Tavazzani et al., 2020; Fig. 5e).Overall, the two RFC simulations produce comparable results, with the rhyolitic recharge scenario (Fig. 6c, f) favored by quantitative comparison of synthetic and empirical zircon age spectra via Kolmogorov-Smirnov statistical test (Press et al., 1988;see Tavazzani et al., 2023, for a statistical match example).

Implications for the interpretation of zircon age dispersion in large silicic systems
Heterogeneous age spectra with late peaks in zircon crystallization are a distinctive feature of rhyolitic ignimbrites generated in calderaforming eruptions (Fig. 7a; e.g.Huckleberry Ridge Tuff, Singer et al., 2014;Wotzlaw et al., 2015;Kneeling Nun Tuff, Szymanowski et al., 2017, 2019).The same pattern of age dispersion is observed in upper-crustal granitoid intrusions that are contemporary and genetically related to centers of rhyolitic volcanism (inset in Fig. 7a; e.g.Organ Mountains, Rioux et al., 2016;SMS, Karakas et al., 2019, this work).Here, based on results from zircon saturation modeling, we present two non-mutually exclusive pathways of late zircon crystallization peak generation in plutonic and volcanic products of silicic magmatic systems.
The first pathway is through zircon growth in an interconnected, continuously recharged reservoir, wherein a significant fraction of melt (>25 vol.%) remains in the system between successive magmatic pulses (Fig. 7b).In this scenario, each new magma injection contributes heat to the system and dilutes Zr concentration in the residual melt (i.e.assuming new injections are undersaturated in zircon) causing delay in zircon mass precipitation.A sharp increase in zircon growth rates is attained only at the end of the recharge cycle.This pathway agrees with the prediction that sizeable eruptions (>0.1 km 3 ) are preceded by incremental assembly and chemical evolution of shallow magma reservoirs (Forni et al., 2018) in systems characterized by intermediate magmatic fluxes (e.g.recharge rates of ~0.001 km 3 yr − 1 ), which allow survival of melt between emplacement of successive batches (Gelman et al., 2013;Tierney et al., 2016;Townsend and Huber, 2020).A second pathway, whose applicability is limited to the interpretation of zircon age spectra in rhyolitic ignimbrites, is related to mechanical mixing of zircon populations crystallized in sequentially emplaced and monotonically cooled magma bodies shortly before (i.e.homogenization; Huber et al., 2011Huber et al., , 2012) ) or during a large eruptive event (simultaneous or sequential tapping; e.g.Pearce et al., 2020).If the majority of magma bodies are emplaced close to eruption, pre-or syn-eruptive mixing can generate volcanic age spectra with an apparent late peak in zircon crystallization (Fig. 7c; Tavazzani et al., 2023).The reported increase in magma flux approaching the climactic eruption stage of caldera-forming magmatic centers (Bouvet de Maisonneuve et al., 2021) make this pathway a viable alternative for the generation of late zircon crystallization peaks in rhyolitic ignimbrites generated in caldera-forming eruptions.
In summary, late zircon crystallization peaks may reflect either thermochemically modulated zircon saturation during open-system crystallization of an interconnected magma reservoir, or the pre/syneruptive mechanical mixing of zircon populations crystallized in monotonically cooled, isolated magma batches.Overall, the increase in crystallized zircon mass with time, when observed in the plutonic or volcanic record of silicic systems (Fig. 7a), reflects long-term magma storage conditions, magma fluxes and short-term perturbations, which would allow the discharge of large volumes of mobile melt in single, catastrophic eruptions.
By contrast, age spectra of silicic intrusive suites not known to be directly associated with large ignimbritic eruptions display a relative decrease in zircon crystallization with time (Fig. 8a).This characteristic shape of plutonic zircon age spectra is matched by synthetic zircon distributions generated in equilibrium and fractional crystallization simulations of monotonically cooled, isolated magma batches (Fig. 8b, c; Samperton et al., 2017;Kirkland et al., 2021).As true closed-system crystallization is unlikely in natural magmatic systems (Cashman et al., 2017), an alternative process capable of producing these zircon age distributions is incremental growth characterized by high instantaneous flux (recharge rates >0.01 km 3 yr − 1 ; e.g.Eddy et al., 2016) in a long-term averaged low flux environment (i.e.<0.001 km 3 yr − 1 ; Caricchi et al., 2016) where no melt remains in the system between successive batches.The absence of melt hinders substantial thermal and compositional interaction between pulses (Barboni et al., 2015) and leads to the monotonic cooling of each separate magma batch, approximating a closed-system crystallization scenario.
It follows that those plutonic systems that are characterized by zircon age spectra showing relative decrease in saturated zircon mass with time (Fig. 8a) are unlikely to have experienced a developmental stage corresponding to a thermochemically interconnected magma reservoir (i.e.mature reservoir, Bouvet de Maisonneuve et al., 2021).The absence of late zircon crystallization peaks also eliminates the possibility of a homogenization event, related to increasing magma flux into the system (Fig. 7c).Overall, the absence of reservoir maturation or late-stage homogenization signals in the plutonic record of many silicic systems implies that long-term and/or transient conditions favorable for the evacuation of large volumes of magma in a single, catastrophic event were never established.If present, volcanic activity associated with these plutons would have accounted for volumetrically limited, distributed eruptions fed by individual magma batches intruded and differentiated in the upper crust.

Implications for eruption age estimation of volcanic units
High-precision geochronology of eruptive products hosted in sedimentary successions is a fundamental tool to calibrate the Earth's geologic history.However, the use of zircon U-Pb dates to define the eruption (or emplacement) age is invariably biased by pre-eruption zircon crystallization, Pb loss and methodological biases (e.g.undersampling of the zircon population, over-representation of late-stage crystallization events; Klein and Eddy, 2023;this work).Therefore, increasingly sophisticated models that rely on statistical likelihood of a group of zircons representing an eruptive event have been proposed and tested on a number of geological objects (Schoene et al., 2019;Keller et al., 2018).These models rely on assumptions about the distribution of zircon mass during crystallization of the source magma body (Keller et al., 2018).A key observation from our study, although it does not specifically treat the kinetics of zircon crystallization from a multi-phase magma (e.g.not-constant growth rate, armoring effects; Bindeman and Melnik, 2016;Klein and Eddy, 2023), is that the mass of zircon crystallized in upper crustal magmatic bodies feeding large eruptions is not normally distributed, with relatively more zircon mass generated towards the end of the lifetime of the system (i.e.eruption or complete crystallization).In light of this observation, we suggest that models such as that of Keller et al. (2018), which allows the selection of a-priori distribution of zircon mass in the melt, should be used with a maximally informed prior distribution reflecting this observation.For large-n samples (e.g.LA-ICP-MS datasets) a prior distribution generated by bootstrapping the data may be the most appropriate; however, in cases where few zircon dates are available to estimate eruption/deposition age (e.g.CA-ID-TIMS datasets), significant biases may be introduced via bootstrapping, or by choosing an inappropriate prior.Our forward modeling (Fig. 6) and data compilation (Figs. 7 and 8) suggest that in absence of other information, the best choice to estimate eruption age of large volcanic explosive events might be approximating the time evolution of mass of zircon crystallized with a truncated normal distribution (Keller et al., 2018), as this distribution reproduces the late peak of zircon crystallization characteristic of long-lived, interconnected magmatic reservoirs (Fig. 7a).

Quantification of undersampling effect on zircon age distribution
The success of quantitative analyses of zircon age spectra relies on the acquisition of a large enough number of zircons per sample, such that the entire history of zircon saturation is represented (Klein and Eddy, 2023).In this context, the stochastic aspect of our modeling template can quantify the effect that the limited number of zircons per sample routinely analyzed with the CA-ID-TIMS technique have on the representativeness of zircon age distributions in natural systems (Samperton et al., 2015).A compilation of zircon age spectra generated by CA-ID-TIMS displays the influence that undersampling has on apparent relative distribution of zircon in time (Figs.7 and 8).Assuming L. Tavazzani et al. equivalent crystallization histories, well characterized systems (i.e.>20 zircons per sample) display accentuated asymmetries in zircon crystallization distribution, whereas systems with a low number of zircons analyzed (i.e.<10 zircons per sample) show smoother age spectra.Our workflow can reproduce the "smoothing" effect by controlling the number of zircons randomly sampled from a synthetic probability distribution (Fig. 7b, c and Fig. 8b, c).Such observation confirms the stochastic nature of undersampling in high-precision U-Pb and roughly indicates a minimum number of 10 zircon dates per igneous sample in order to assign petrological significance to any age dispersion.

Conclusion
Comparison of high-precision (CA-ID-TIMS) zircon age spectra from volcanic and plutonic products of large silicic systems with results of a thermodynamic-based zircon saturation model highlights the link between magmatic flux, zircon crystallization history and eruptive potential of shallow silicic magmatic reservoirs: (1) The dominant crystallization process (i.e.equilibrium or fractional crystallization) controls the shape of zircon age spectra during monotonic cooling and crystallization of isolated magma batches.In systems characterized by extensive, supra-solidus interaction between successive magma recharges the frequency, volume and composition of new magma injections affect zircon growth curves.
(2) In evolved, subalkaline systems, an increase in relative number of crystallized zircons with time denotes 'active', thermochemically connected magma reservoirs, where increase in heat supply can produce system-wide crystallinity variations conducive to evacuation of large volumes of magma in a single, voluminous eruption.Conversely, a decrease in crystallized zircon mass with time identifies 'passive' systems, which are incrementally assembled by rapid emplacement of individual magma batches, but never experience system-wide crystallinity variations allowing for the segregation and release of large mobile magma volumes.(b, c) Synthetic zircon age spectra from fractional crystallization (FC) and equilibrium crystallization (EQ) simulations.Zircon crystallization probability (p Zr ) evolution is shown in the top panels of each figure and each roman numeral represents a new episode of recharge during the evolution of the system.The number of zircons (N Zr ) randomly sampled from synthetic FC and EQ probability density functions f(p Zr ) is variable between 5 and 30.Conceptual sketches are adapted from Szymanowski et al. (2017) and Forni et al. (2018).References: 1. Bergell pluton (Samperton et al., 2015); 2. Mt.Capanne pluton (Barboni et al., 2015); 3. La Gloria pluton (Gutiérrez et al., 2018); 4. Golden Horn Batholith (Eddy et al., 2016).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Location of the study area and sampling sites.(a) Schematic geologic map of the Sesia Magmatic System (SMS) and its location in the Ivrea-Verbano and Serie dei Laghi units of southern Alpine basement (northern Italy), redrawn after Handy and Zingg (1991).(b) Geologic map of the upper crustal segment of the SMS to the South of the Cremosina Line with localities of samples used in hugh-precision geochronology study, after Tavazzani et al. (2020).

Fig. 2 .
Fig. 2. Zircon saturation modeling workflow, exemplified for closed-system (equilibrium crystallization -EQ, fractional crystallization -FC) and open-system (recharge fractional crystallization -RFC; Bohrson et al., 2014) crystallization scenarios.(a) Normalized zircon saturated mass distribution m Zr (n * ) derived for the three crystallization scenarios.(b) Cumulative zircon crystallization probability distribution functions F pZr (n *) derived from integration of m Zr (n * ) over the zircon saturation interval.(c) Three synthetic zircon age datasets (number of zircons, N Zr = 10) extracted from F pZr (n * ) curves of EQ, FC and RFC simulations.In these datasets the ratio of true crystallization timescale (ΔT) to analytical uncertainty (σ) is equal to 20 to simulate results of CA-ID-TIMS analyses (e.g.ΔT = 1 Myr and σ = 50 kyr, corresponding to the typical precision of our analytical data).(d) Relative zircon crystallization distributions f xtal (t r ) (i.e.age spectra) scaled from initiation to termination of zircon crystallization, shown as kernel density estimates of EQ, FC and RFC simulations at N Zr = 30.The solid black line in (d) corresponds to the kinetic zircon crystallization model ofWatson (1996), based on zirconium diffusion constraints.(For interpretation of the colors in the figures, the reader is referred to the web version of the article).

L.
Tavazzani et al.  simulation).A zircon crystallization probability density function f pZr (n) is defined for each simulation, with non-uniform distribution:

Fig. 3 .
Fig. 3. Geochronology of the Valle Mosso Pluton and related silicic volcanic rocks.(a) Rank-order plot of CA-ID-TIMS zircon dates from the SMS.Colored bars represent Th-corrected CA-ID-TIMS 206 Pb/ 238 U dates of individual zircon crystals with their 2σ uncertainty.Analyses interpreted as antecrysts are semi-transparent.Δt h is the difference between the oldest and youngest zircon dates in a sample; Δt a is the difference between oldest and youngest autocrystic zircon dates.Eruption ages are defined for each volcanic sample as the Bayesian estimate(Keller et al., 2018) using the bootstrapped distribution of Th-corrected 206 Pb/ 238 U dates as prior, with 2σ uncertainty given as internal only/with tracer calibration/ with tracer and 238 U decay constant.Individual zircon dates for crystal-rich rhyolite sample MBR and range of zircon dates from 4 granitoids of the Serie dei Laghi are reproduced fromKarakas et al. (2019).(b, c) Temporal distribution of individual zircon ages reported as kernel estimates of relative zircon crystallization density functions f xtal (t r ), where t r is relative crystallization time, normalized between x = 1, denoting zircon saturation (t sat ) and x = 0, corresponding to full plutonic solidification (t sol ) or eruption (t er ).The two kernel density estimates represent zircon crystallization distributions bootstrapped over all dated samples in the eruptive (b) and intrusive (c) domains of the SMS.

Fig.
Fig. Zircon trace element co-variation diagrams.LA-ICPMS trace element data for zircons also analyzed by CA-ID-TIMS (a), (d) overlap with SMS upper-crustal zircon compositions from Karakas et al. (2019).Data are color-coded Ti-concentration in zircon in (b), (e).Zircon crystallization temperatures are calculated based onFerry and Watson (2007) using an a TiO2 of 0.6 and a SiO2 of 1.0 (seeTavazzani et al., 2020 for details on activity coefficients determination).In (c), (f) fractional crystallization models are superimposed on trace elements co-variation diagrams.The median composition of zircons with the highest Ti-in-zircon crystallization temperature is selected as initial composition.Percent crystallinity (20%, 50%, etc.) is indicated along the fractionation trajectories.Partition coefficients used and details of each crystallization models are provided in the Supplementary Material.The two models in (c) use either constant partition coefficients ("Constant D") or partition coefficients that vary with temperature ("Variable D") during the course of crystallization.Crystallizing mineral assemblage constitutes of sanidine, plagioclase and quartz in near-eutectic proportions(Luth et al., 1964).In (f), REE evolution scenarios are evaluated as a function of the amount of zircon (zrn), apatite (ap) and titanite (ttn) in the crystallizing mineral assemblage.

Fig. 5 .
Fig. 5. Temporal variations of trace elements in the upper crustal part of the SMS zircons, with mid-and lower-crustal data from Karakas et al. (2019) shown for comparison.Each point indicates one in-situ LA-ICPMS analyses (i.e.core or rim) acquired prior to bulk-crystal CA-ID-TIMS U-Pb dating, plotted against individual zircon crystallization age.Box plots show median, interquartile ranges (IR) and extreme values for bins of 100 kyr (7 youngest bins) or longer (200 kyr for dates >283.80Ma) calculated on rim and core compositions.Vertical lines indicate eruption ages defined for each volcanic sample as the Bayesian estimate of the zircon Th-corrected 206 Pb/ 238 U dates bootstrapped distribution(Keller et al., 2018) with 2σ uncertainty given as internal only/with tracer calibration/with tracer and 238 U decay constant.Ti-in-zircon crystallization temperatures estimates use the model ofFerry and Watson (2007) and activity values fromTavazzani et al. (2020).

Fig. 6 .
Fig. 6.Comparison of SMS zircon CA-ID-TIMS crystallization age distributions and synthetic zircon distributions obtained through zircon saturation modeling from three different magma crystallization scenarios.(a, d) Result of fractional crystallization (FC) simulations.(b, e) Recharge and fractional crystallization (RFC) models, employing a sequence of three basaltic andesite recharges.(c,f) Recharge and fractional crystallization (RFC) simulation using a rhyolitic recharge magma in three consecutive injection episodes.Zircon age distributions, expressed as relative zircon crystallization density functions f xtal (t r ), where t r is relative time, scaled from zircon saturation (t sat ) to full plutonic solidification (t sol ) or eruptive truncation (t er ), are shown as kernel density estimate truncated at ± 0.5 kernel bandwidth(Keller et al., 2018).Each synthetic zircon distribution (i.e.bold black lines) is the product of random sampling of 30 zircons from synthetic FC and RFC distributions averaged over 10 3 stochastic iterations.Temperature/crystallinity evolution and increments of magma mass over scaled crystallization time are shown in the top panels of each figure.Crystallinity and solidus temperature are based on rhyolite-MELTS(Gualda et al., 2012) derived temperature crystal-fraction curve (i.e.solidus-liquidus relation) for starting granitic composition used in the simulation(Tavazzani et al., 2023).

(Fig. 7 .
Fig. 7. Conceptualization of the effect of thermal and mechanical processes on the shape of CA-ID-TIMS zircon age spectra in subvolcanic magma reservoirs.(a) Empirical kernel density estimates of zircon crystallization distributions from four rhyolitic ignimbrites generated in catastrophic caldera-forming eruptions (N Zr : range of number of zircons analyzed per sample for each eruption).In the inset panel, zircon crystallization distribution form two uppercrustal intrusions contemporary and genetically related to caldera-forming silicic volcanism.(b) Synthetic zircon age spectra output from fractional crystallization (RFC) simulations employing a sequence of three rhyolitic recharges in a cooling and differentiating magma reservoir.(c) Age spectra resulting from pre-or syn-eruptive mixing of zircon populations from sequentially emplaced, distinct magma batches during a caldera-forming event.The kernel density distributions are generated through the randomized sampling of zircon populations from four individual rhyolitic magma batches emplaced, monotonically cooled and differentiated in the upper crust (FC).Zircon crystallization probability (p Zr ) evolution is shown in the top panels of each figure and each roman numeral represents a new episode of recharge during the evolution of the system.The number of zircons (N Zr ) randomly sampled from RFC and composite FC+homogenization probability density functions f(p Zr ) is variable between 5 and 30.Conceptual sketches are adapted from Szymanowski et al. (2017) and Forni et al. (2018).References: 1. Kneeling Nun Tuff (KNT, Szymanowski et al., 2019); 2. Huckleberry Ridge Tuff (HRT, Singer et al., 2014; Wotzlaw et al., 2015); 3. Sesia Caldera (SC, Karakas et al., 2019, this work); 4. Organ Needle caldera (ONC, Rioux et al., 2016); 5. Valle Mosso pluton (VMP, this work); 6. Organ Needle pluton (ONP, Rioux et al., 2016).

Fig. 8 .
Fig.8.Summary of petrologic implications of the shape of zircon age spectra in silicic intrusive suites.(a) Bootstrapped zircon crystallization distributions from four silicic intrusive suites (N Zr : range of number of zircons analyzed per sample for each intrusion).(b, c) Synthetic zircon age spectra from fractional crystallization (FC) and equilibrium crystallization (EQ) simulations.Zircon crystallization probability (p Zr ) evolution is shown in the top panels of each figure and each roman numeral represents a new episode of recharge during the evolution of the system.The number of zircons (N Zr ) randomly sampled from synthetic FC and EQ probability density functions f(p Zr ) is variable between 5 and 30.Conceptual sketches are adapted fromSzymanowski et al. (2017) andForni et al. (2018).References: 1. Bergell pluton(Samperton et al., 2015); 2. Mt.Capanne pluton(Barboni et al., 2015); 3. La Gloria pluton(Gutiérrez et al., 2018); 4. Golden Horn Batholith(Eddy et al., 2016).