AgeSpectraAnalyst: A MATLAB based package to model zircon age distributions in silicic magmatic systems

In the last decade, improvements in the analytical precision achievable by zircon U-Pb geochronological techniques have allowed to resolve complexities of zircon crystallization histories in magmatic rocks to an unprecedented level. A number of studies have strived to link resolvable dispersion in zircon age spectra of samples from fossil magmatic systems to the physical parameters of their parent magma bodies. However, the methodologies developed have so far been limited to reproduce the effect of simple thermal histories on the final distribution of zircon ages. In this work we take a more nuanced approach, fine-tuning a thermodynamics-based zircon saturation model to predict the relative distribution of zircon ages in samples from silicic magma reservoirs experiencing open-system processes (e.g. heat/mass addition, mechanical mixing). Employing the MATLAB package (AgeSpectraAnalyst) presented in this contribution:• Users can forward model the effect that diverse thermal histories and mechanical mixing processes characteristic of silicic magma bodies have on zircon age distributions as measured by high-precision, chemical abrasion thermal ionization mass spectrometry (CA-ID-TIMS) U-Pb geochronology.• Zircon CA-ID-TIMS datasets from silicic magmatic systems can be easily compared with model output to gain semi-quantitative information on thermo-mechanical history of the system of interest.• We demonstrated (Tavazzani et al., in press) that distribution of high-precision zircon ages in crystallized remnants of shallow (∼ 250 MPa), silicic magma reservoirs can discriminate between systems that experienced catastrophic, caldera-forming eruptions and systems that underwent monotonic cooling histories.


a b s t r a c t
In the last decade, improvements in the analytical precision achievable by zircon U-Pb geochronological techniques have allowed to resolve complexities of zircon crystallization histories in magmatic rocks to an unprecedented level.A number of studies have strived to link resolvable dispersion in zircon age spectra of samples from fossil magmatic systems to the physical parameters of their parent magma bodies.However, the methodologies developed have so far been limited to reproduce the effect of simple thermal histories on the final distribution of zircon ages.In this work we take a more nuanced approach, fine-tuning a thermodynamics-based zircon saturation model to predict the relative distribution of zircon ages in samples from silicic magma reservoirs experiencing open-system processes (e.g.heat/mass addition, mechanical mixing).Employing the MATLAB package (AgeSpectraAnalyst) presented in this contribution: • Users can forward model the effect that diverse thermal histories and mechanical mixing processes characteristic of silicic magma bodies have on zircon age distributions as measured by high-precision, chemical abrasion thermal ionization mass spectrometry (CA-ID-TIMS) U-Pb geochronology.• Zircon CA-ID-TIMS datasets from silicic magmatic systems can be easily compared with model output to gain semi-quantitative information on thermo-mechanical history of the system of interest.• We demonstrated (Tavazzani et al., in press) that distribution of high-precision zircon ages in crystallized remnants of shallow ( ∼ 250 MPa), silicic magma reservoirs can discriminate between systems that experienced catastrophic, caldera-forming eruptions and systems that underwent monotonic cooling histories.

Background
The temporal evolution of magmatic systems can be revealed by analyses of refractory minerals incorporating radiogenic isotopes that decay with time.U-Pb in zircon has long been recognized as the most robust system for this purpose [3] and the improvement of analytical protocols and instrumentation design (e.g.[4][5][6][7] ) now allows determination of absolute ages with best-case internal analytical uncertainties better than 0.02 % for bulk-grain analyses of zircon through chemical abrasion-isotope dilution-thermal ionization mass spectrometry (CA-ID-TIMS; [8] ).Increase in precision means that for individual samples dispersion of single analyses outside of analytical uncertainty is commonly observed.The analysis of zircon age distributions to reconstruct the magmatic flux in a particular magmatic system is now a well-established branch of zircon petrochronology (e.g.[9][10][11][12] ).However, some of the major limitations of this approach are: (1) the modeling of zircon mass distribution at the scale of an entire magmatic system, rather than at the sample scale, inevitably blurring the heterogeneous evolution of incrementally assembled reservoirs; (2) the assumption of a constant zircon crystallization rate during the lifetime of a magmatic system or that (3) this crystallization rate can be parametrized based on theoretical assumptions [13] .In order to circumvent these limitations, some authors have employed thermodynamicsbased zircon saturation models, which integrate temperature, crystallinity, major and trace element evolution of a magma batch to assign petrologic meaning to dispersed zircon age populations observed in rock samples (e.g.[ 1 , 2 , 14-16 ]).So far, these modeling efforts have been successful but have seen their application limited to relatively simple cases of monotonic cooled magma batches.A new generation of models, which include the effect of open-system process on zircon saturation, is needed in order to interpret the relative distribution of zircon mass (i.e.zircon age distributions) in samples collected from the crystalized remnants of long-lived, thermochemically complex magmatic systems.

Rhyolite-MELTS zircon crystallization distribution
Model zircon crystallization distributions are calculated from rhyolite-MELTS [ 25 ] major-element equilibrium crystallization (EQ), fractional crystallization (FC) and recharge fractional crystallization (RFC) simulations (run on the software Magma Chamber Simulator, https://mcs.geol.ucsb.edu; [ 26 , 27 ]) using the saturation model of Boehnke et al. [ 28 ].A rhyolitic composition (from [ 29 ], given in Table 1 ; Fig. 1 ) is used as starting melt composition at a pressure of 250 MPa.This pressure corresponds to the optimal depth range for magma chamber nucleation and growth as suggested by physical models [ 30 ] and a compilation of natural datasets [ 31 ].Initial magmatic water content is set to 2.5 wt% H 2 O and ƒO 2 to the QFM (quartz-fayalite-magnetite) buffer.Increase of H 2 O content above 2.5 wt% does not substantially affect the distribution of saturated zircon mass ( Fig. 1 ).The effect on zircon saturation generated by the presence of a mixed volatile magmatic phase (e.g. with additional CO 2 , H 2 S, SO 2 , etc.) is not evaluated in this work.In RFC simulations, rhyolitic and basaltic andesite compositions (also from [ 29 ], given in Table 1 ) are used as liquid compositions of silicic and mafic recharges, respectively.A complete list of crystallization simulations explored in this work is provided in Table 2 .Each crystallization simulation is allowed to reach solidus temperature or alternatively is interrupted at 750 °C ( c. 50 vol.%melt; Fig. 1 ), to simulate the abrupt termination of crystallization caused by an eruptive event.
At each crystallization simulation step the concentration of zirconium in the residual melt fraction ( [ Zr ] melt calc ) is calculated through the Rayleigh equation for fractional crystallization:
Where F is the melt fraction value output from rhyolite-MELTS, [  ]   the zirconium concentration in the initial liquid ( Table 1 ), and D the bulk Zr partition coefficient determined using rhyolite-MELTS modal mineral abundances and Zr partition coefficients compiled from the Geochemical Earth Reference Model (see Supplementary Table 1).The direct relationship between zircon saturation [ 28 ], temperature, and melt composition is expressed as: where the T is temperature (in K, output from rhyolite-MELTS) and M the parameter describing degree of silicate melt polymerization (i.e. the balance between network modifying cations: Na, K, Ca, Mg, Fe, and silicate melt network formers: Si, Al, which affects zirconium solubility in a melt) calculated after Watson and Harrison [ 18 ] as: where Na, K, Ca, Al, and Si are expressed as molar, normalized element components output from rhyolite-MELTS.Zircon growth initiates when [ Zr ] melt calc > [ Zr ] melt sat and the mass of zircon crystallized is determined at each simulation interval as: where m liq is the mass of liquid phase output from rhyolite-MELTS (in grams).Trivial differences in crystallized zircon mass distributions are observed when tying dependence of zircon saturation on melt temperature and composition ( Eq. ( 2)) with the most recent Crisp and Berry [ 20 ] experimental calibration.This workflow produces a distribution of zircon growth as a function of temperature and melt fraction in batch crystallization and fractional crystallization simulations (i.e.zircon crystallization spectra, Fig. 1 ; [ 1 , 14 ]).However, in complex crystallization simulations (RFC), recharge-induced fluctuations of intensive parameters cause a non-linear relationship between mass of zircon crystallized and temperature (or melt fraction), with individual temperature values (or melt fraction values) corresponding to multiple values of crystalized zircon mass.To address this problem, we express the zircon crystallization distribution as a function of the number of steps (n) from crystallization initiation ( Fig. 2 ).In rhyolite-MELTS simulations, if an equal amount of enthalpy (e.g.1Jg − 1 ) is extracted at each step from the magma during solidification [ 34 , 35 ] the number of modelling steps becomes proportional to the time spent by the magma within a given interval of temperature.Under this assumption we thus have a relative crystallization chronology, which provides a time-integrated record of zircon growth/dissolution between each scenario of interest.

Stochastic sampling of rhyolite-MELTS zircon crystallization distributions
The limited number of analyses per specimen obtained with high-precision geochronological techniques (e.g.CA-ID-TIMS) may skew natural zircon age distributions towards intervals of high zircon crystallization at the expense of relatively low zircon growth periods.To provide a realistic comparison between natural (i.e.CA-ID-TIMS) and synthetic zircon age populations, we developed a workflow that allows the random sampling of an arbitrary number of zircons from a model zircon crystallization probability density function (generated from rhyolite-MELTS simulations following the numerical approach described in the previous section).To this end, the MATLAB R ○ script StochasticZirconCrystallization.m implements the stochastic sampling of rhyolite-MELTS derived zircon crystallization distribution (listed in Table 2 ).  2 ) and (d), (e), (f) recharge fractional crystallization (simulation ID: 14-3d from Table 2 ) simulations, respectively.
The probability of generating a zircon at a particular step ( n ) of a crystallization simulation is defined as: where dm   is the relative change in available zircon mass m   ( Eq. ( 4)) between two successive simulation steps: which is normalized to the overall variation in crystallized zircon mass ( ∑   =0    , where N is the total number of steps in a crystallization simulation).As our model does not account for negative zircon growth (i.e.zircon dissolution), any simulation step defined by negative change in zircon mass is set to dm   = 0 .
A zircon crystallization probability density function    (  ) is thus defined for each simulation, with non-uniform distribution: A number of zircons is randomly sampled from    (  ) to simulate the stochastic process of zircon selection during geochronological analysis.The cumulative distribution function    (  ) of the target probability density function    (  ) is obtained as: and the variable of interest (in this case n , a specific crystallization step) is randomly selected from the cumulative distribution function    (  ) as:  2 ) and (e), (f), (g), (h) recharge fractional crystallization (simulation ID: 11-5f from Table 2 ) simulations.From left to right the number of stochastically sampled zircon crystallization steps (n) per simulation is increased from N Zr = 30 to N Zr = 10000.The random selection algorithm undergoes 10 3 iterations for each simulation.
The n variables randomly selected from a non-uniform zircon crystallization probability density function    (  ) correspond to specific simulation steps at which zircons are observed ( Fig. 2 ).This set of equations is implemented using the MATLAB R ○ function randsrc().The number of randomly selected n variables (N Zr ) can be defined by the user ( Figs. 2 , 3 ).For each simulation, as to ensure robustness of the random selection process, the random selection algorithm undergoes 10 3 iterations.

Visualization of model output -zircon crystallization age density functions (natural and synthetic)
Zircon age distributions are expressed as relative zircon crystallization density functions f xtal (t r ) (after [2] ), where t r is relative time, scaled from zircon saturation ( t sat ) to full plutonic solidification ( t sol ) or eruptive truncation ( t er ): While in natural samples t sat and t sol or t er correspond to absolute radiometric dates, in the case of model zircon distributions t sat corresponds to the first simulation step in which [ Zr ] melt calc > [ Zr ] melt sat and t sol or t er correspond to the simulation step at which the simulation is interrupted.
The form of relative crystallization distributions f xtal (t r ) estimated from natural or synthetic data are then visualized through a truncated kernel density estimate of the scaled crystallization times    where: To produce the kernel density estimate of f xtal (t r ) from    we use the MATLAB R ○ function ksdensity() with a Gaussian kernel and bandwidth determined by Gaussian approximation.The resulting kernel density estimate is truncated at t r = -0.05for natural volcanic samples and model f xtal (t r ) produced by simulations interrupted at 750°C (i.e.modeling pre-eruptive zircon crystallization).Kernel density estimates are then scaled between x = 1, which denotes the initiation of zircon crystallization at t sat and x = 0 the termination of crystallization at t sol or t er ( Figs. 2 , 3 ).This methodology (first introduced by [2] ) allows a reliable estimate of f xtal (t r ) , which can reproduce broad fluctuations in the zircon relative crystallization rate, while maintaining the abrupt cutoff in the crystallization distribution that should be caused by an eruption.The kernel density estimates of natural, CA-ID-TIMS datasets represent averaged zircon crystallization distributions of all the samples collected from a single magmatic system.In synthetic datasets, kernel density  2 ), and for (b) a recharge fractional crystallization simulation (ID: 14-3d from Table 2 ).(c), (d) Box plots of temperature values from the same simulations.
In panel (e) and (f), the randomly sampled zircon crystallization distributions are displayed as kernel density estimates of the scaled simulation steps (n * ).Synthetic zircon age and temperature distributions are scaled from initiation (x = 1) to termination (x = 0) of zircon crystallization.Box plots show median, interquartile ranges (IR) and extreme values for bins corresponding to 1/10 of the total zircon crystallization history.
estimates are obtained for both zircon age distributions from individual stochastic outputs and zircon age distributions averaged over 10 3 stochastic outputs.This is done to better display the range of possible zircon crystallization distributions and associated kernel density estimates that are obtainable from the same underlying    (  ) through a stochastic sampling process.

Zircon crystallization temperature distribution
In the StochasticZirconCrystallization.m routine, each simulation step (n) that is randomly selected from a non-uniform zircon crystallization probability density function    (  ) has a corresponding temperature value.This value represents the temperature of zircon crystallization at a specific simulation step, allowing us to track absolute changes in zircon crystallization temperature in a rhyolitic magma batch, in different crystallization pathways.To take advantage of this perk of the rhyolite-MELTS engine, we implemented an additional routine [CrystallizationTemperatureEvolution.m] that returns, in addition to the kernel density estimate of a synthetic zircon population, the distribution of zircon crystallization temperatures for the same simulation ( Fig. 4 ).To account for the analytical uncertainties associated with the application of Ti in zircon thermometry to natural samples (c.f.[ 36 ]), the user can add randomized noise (e.g.analytical scatter, uncertainty on titanium activity, on the order of ± 0 to 25 °C; [ 11,37 ]) to temperature values obtained for each simulation step.

Uncertainty parametrization in the generation of CA-ID-TIMS synthetic zircon populations
To parameterize variation of U-Pb zircon age uncertainties with absolute age of zircon crystallization [ 38 , 39 ] and thus simulating an actual CA-ID-TIMS output, we: (1) Compiled ∼1700 single-crystal, Th-corrected CA-ID-TIMS zircon 206 Pb/ 238 U dates from 1 to 1000 Ma published between 2011-2021 (Supplementary Table 2; see reference list under Fig. 5 ) together with their reported analytical uncertainties in Ma (2 ).(2) Applied a linear regression (least square, LSQ) to the compiled dataset in the interval 1 to 1000 Ma to obtain an empirical relation between 206 Pb/ 238 U age and analytical 2  uncertainty.As analytical uncertainties in U-Pb dates are symmetric, we only plotted one data point per analysis.Given the spread of the published uncertainties within individual time bins (up to one order of magnitude), the goodness of fit is lacking (r 2 = 0.40).Fitting the dataset with higher order polynomial functions does not seem to improve the regression procedure ( Fig. 5 a).A least square regression of this large dataset provides a first order indication of the evolution of typical analytical uncertainties in CA-ID-TIMS 206 Pb/ 238 U dates with zircon age.[ 21,22,37,.Each line represents a different polynomial fit to the dataset, and points on a line corresponding to interpolated 2  uncertainties for three magmatic systems of increasing crystallization ages ( ∼ 7 Ma, Capanne pluton, [ 37 ]; ∼ 95 Ma, Mt.Stuart pluton, Matzel et al. [ 66 ]; ∼285 Ma, Laghi granites, [ 41 ]) (b), (c), (d) Visualization of the output of recharge fractional crystallization simulations (simulation ID: 14-3d from Table 2 ) with added analytical uncertainty corresponding to the three age brackets at 7, 95 and 285 Ma (thin grey lines) compared to an average of 10 2 iterations for the same crystallization scenario without added uncertainty (bold red line).
(3) To add the effect of analytical uncertainty corresponding to zircons crystallized at a specific time, the model calculates the median age of the input dataset and uses that value to calculate the characteristic uncertainty from the linear regression ( Fig. 5 a).A random value picked from the Gaussian uncertainty distribution is then added to each randomly selected zircon date from the cumulative distribution function    (  ) .
The increasing absolute uncertainty with increasing zircon crystallization age (here only 206 Pb/ 238 U ages are discussed) is reflected in the scatter of resampled populations.In simulations of zircon crystallization in younger systems ( < 10 Ma) the effect of accounting for analytical uncertainty during probability distribution resampling is nearly insignificant ( Fig. 5 b).This effect is still very limited in Late Mesozoic magmatic systems ( < 100 Ma; Fig. 5 c), however it becomes important in simulating zircon age distributions in Early Mesozoic to Paleozoic contexts ( Fig. 5 d).

Modeling mechanical mixing of zircon populations
In the run-up to or during a large eruption, mechanical mixing and homogenization of distinct zircon populations (i.e.generated from the cooling of isolated magma batches, or through incorporation of older magmatic products during eruption or transport within a pyroclastic density current; Fig. 6 a) can create zircon age populations that do not follow a simple distribution expected from crystallization of a single magma batch.For example, pre-eruptive recharge of an upper crustal magma reservoir may lead to entrainment of various crystallized portions of the system (e.g.[ 24 ]).Alternatively, during an eruptive event, the zircons hosted in rocks in proximity of eruptive vents could be taken as crystal-cargo ( "chimney sweeping " effect; [ 67 ]).The zircon populations that can be created through such mixing will have age distributions that cannot be reproduced with a purely thermodynamic model.Here the stochastic nature of the mechanical mixing process is simulated by juxtaposing and randomly sampling a number of zircon populations that have a parametrized temporal distribution, with an additional MATLAB R ○ routine [Resampling_for_mixing.m] that allows the simulation of complex mixing of zircon populations: (1) The model is based on zircon growth curves for a number of temporal intervals (n = 7) that correspond to CA-ID-TIMS dated samples of the Bergell pluton ( Fig. 6 a, b; [ 68 ]).The temporal distribution of samples in this intrusion represents a benchmark example of incrementally added magma batches, and nearly every sample shows a zircon age distribution that is expected for monotonic cooling of a magma batch [ 14 ].(2) From a number of zircon populations (user selected) the model creates a cumulative crystallization probability density function    _  (  ) .(3) The    _  (  ) is randomly sampled to extract a synthetic zircon population in a similar fashion to that with FC and RFC simulations (previous section), with the user being able to select the number of randomly selected n variables (N Zr ).
The model results show that the number of independent zircon populations involved in the mixing process is key to the final shape of the zircon age distribution.In the case of only two populations sampled in the mechanical mixing process, the output will Fig.6.Simulation of zircon age spectra created by mechanical mixing of individual zircon populations.(a) Natural distribution of zircon dates from different samples of the Bergell pluton [ 68 ] and (b) bootstrapped zircon age distributions for each one of these samples.The three side panels show the effects of mixing (c) 2, (d) 3 or (e) 6 individual zircon age populations with age distributions identical to that observed in the Bergell samples (thin orange lines, 50 simulations for each panel).be very close to the starting (i.e.Bergell samples) zircon distributions ( Fig. 6 c) with a shape of zircon age spectra similar to that of monotonically cooled magma batches, while a mixing that involves nearly all zircon populations will produce a completely different zircon distribution ( Fig. 6 e), with an apparent shift in zircon relative abundance.

Method validation
The chief goal of this method is to gain a better understanding of the genetic processes beyond the diverse zircon age distributions that are observed in granitic and rhyolitic lithologies.To this end, we compare (a) mass spectrometry-derived zircon age distributions from a range of natural samples and (b) stochastically derived zircon age distributions, sampled from synthetic zircon crystallization distributions.As a proof of concept, we selected three long-lived ( > 1 Myr) silicic magmatic systems (Bergell, [ 68 ]; Bingham Canyon, [ 40 ]; Sesia, [ 41 , 69 ]) that have been recently characterized both by CA-ID-TIMS geochronology and Ti-in-zircon thermometry.
As a first example, we ran a Kolmogorov-Smirnov (KS; [ 70 ]) statistical test to quantify the match between synthetic zircon crystallization density functions ( f xtal (t r ) ) obtained from model zircon crystallization distributions and natural f xtal (t r ) from the Bergell [ 68 ] and Valle Mosso (Tavazzani et al.,in press [ 69 ]) intrusive complexes.The KS test determines the probability that two data sets, in this case a synthetic and a natural zircon age distribution, are generated from the same data distribution by comparing their density functions.This test yields an acceptably high KS probability (usually p > 0.05) for both data sets being drawn from the same underlying distribution, whereby the null hypothesis that they are drawn from an identical density function is rejected when p < 0.05 [ 70 , 71 ].The performance of a modeled crystallization scenario is quantified in terms of percent of stochastic outputs that pass a KS test when compared to empirical zircon age distributions ( Fig. 7 ).In this example, the averaged zircon age distribution of Bergell pluton samples [ 68 ] is compared against outputs of a fractional crystallization simulation ( Fig. 7 a-d) and the average zircon age distribution of the Valle Mosso pluton [ 69 ] is evaluated against the synthetic zircon age distributions produced by a recharge fractional crystallization model ( Fig. 7 e-h).In the case of the Bergell pluton, the good match between simulated fractional crystallization zircon age distributions and measured zircon age distribution ( Fig. 7 c, d) agrees with overwhelming evidence for a simple, monotonic cooling of the different units of the intrusion [ 14 ].Conversely, the good match between outputs of open-system (RFC) zircon crystallization simulations and the averaged zircon age spectra from the Valle Mosso pluton ( Fig. 7 g, h) points to a more complex thermal history for this upper crustal intrusion, which likely experienced fluctuations in mass of zircon available for precipitation during its lifetime.Fig. 7. Statistically evaluated performances of two zircon crystallization simulations compared to empirical zircon age distributions from two plutonic systems.Panel (a), (e) display visual comparisons and (b), (f) results of Kolmogorov-Smirnov (KS) test application to cumulative zircon crystallization distributions of averaged stochastic outputs and plutons age spectra.The degree of similarity between model and empirical age spectra is quantified in terms of percent of stochastic outputs that pass (c), (g) or fail (d), (h) a KS test with significance level at 80 % (i.e.p < 0.2).Natural zircon age distributions are from (a), (b), (c), (d) the Bergell pluton (data from [ 68 ]) and (e), (f), (g), (h) the Valle Mosso pluton [ 69 ].
Secondly, we compared synthetic zircon crystallization temperature evolution obtained for open and closed system crystallization scenarios to the time-resolved record of temperature-modulated trace element concentration (Ti-in-zircon; [ 72 ]; Fig. 8 ) measured in natural zircon populations (e.g.[ 24 , 40 ]).Monotonically decreasing zircon crystallization temperatures characteristic of FC simulations ( Fig. 8 a) agree with observed zircon temperature distributions in natural systems characterized by emplacement and monotonic cooling of an evolved magma body at upper crustal conditions (e.g.Bingham Canyon intrusion, [ 40 ]; Fig. 8 c).Conversely, the convergence between systematic fluctuations in crystallization temperatures obtained from RFC simulations ( Fig. 8 b) and Ti-inzircon temperature evolution observed across the Valle Mosso pluton (upon application of appropriate activity parameters: a TiO2 = 0.6, a SiO2 = 1.0; from [ 29 ]; Fig. 8 d), suggest a complex history with open system processes (i.e.recharge) dominant during the lifetime of this reservoir.

Method strengths and limitations
In contrast to other studies (e.g.[ 10 , 12 ]), the model presented here is effectively adimensional (i.e.0D).It does not provide absolute temporal constraints on the crystallization process and does not contribute information about the geometry of the system or its spatial heterogeneities.The adimensional aspect of our model with respect to absolute crystallization time is exploited to compare saturation modeling results to zircon age spectra from systems of various sizes and operating over different timescales.In MELTSbased simulations, if an equal amount of enthalpy (e.g., 1 Jg − 1 ) is extracted at each step from magma during crystallization, then the number of modeling steps becomes proportional to the time spent by the magma within a given temperature interval [ 34 , 35 ].Under this assumption, a relative crystallization chronology of zircon mass distributions is obtained in our saturation simulations.
In the saturation model, crystallizing zircons are assumed to be in equilibrium with the surrounding melt at every step of the simulation.We have not included a quantification of the effect that kinetic, disequilibrium processes (i.e.diffusion) have on zircon growth and stability [ 73 ].Although the role of dissolution [ 74 ] is neglected, we can qualitatively discuss the effect that transient zircon undersaturation has on the temporal distribution of zircon mass in a periodically recharged magmatic system.Overall, the net effect of periodic zircon undersaturation in RFC scenarios is to bolster the relative zircon crystallization peak close to the end of the magmatic life cycle, through partial loss of early crystallized zircon mass.An exception presents itself when undersaturation can effectively dissolve all early crystallized zircons.In this case, the shape of zircon age spectra will reflect the dominant crystallization process operating in the system after the last magmatic recharge, thus corresponding to an incomplete record of zircon crystallization.However, this occurrence should be rarely observed in evolved subalkaline melts ( sensu [ 75 ]).In melts of this chemical affinity, zircon saturation happens at relatively low Zr concentrations and remelting of zircon-bearing, early cumulate phases during recharge events tends to perpetuate the saturation behavior [ 76 ].Only substantial additions of mafic to intermediate melts, which are rare or subtle in erupted and intrusive products of large silicic systems (e.g.[ 77 ]), can produce significant undersaturation.On account of this, results of the saturation model presented in this work should only be applied to evolved, subalkaline systems (e.g.metaluminous  2 ) and (b) recharge fraction crystallization (simulation ID: 3-1b from Table 2 ) scenarios.Two contrasting natural zircon crystallization temperature distributions are also shown: (c) Bingham Canyon intrusion [ 40 ] and (d) Sesia Magmatic System [ 69 ].To allow direct comparison, natural and synthetic data are scaled from initiation (x = 1) to termination (x = 0) of zircon crystallization.Box plots show median, interquartile ranges (IR) and extreme values for bins corresponding to 1/10 of the total zircon crystallization history.and peraluminous rhyolites) that ensure long-term preservation of zircon crystals throughout the recharge and remelting cycles experienced by an upper-crustal magma reservoir [ 76 ].
Finally, the effect of competition between growth of existing crystals and nucleation during intervals of zirconium oversaturation [ 74 ] is not included in our model.All the zircon mass available for precipitation contributes to crystallization of new zircon grains (i.e.varying nucleation rates but instantaneous growth; [ 78 ]).Moreover, in natural datasets, the inherent mixing of age domains within whole grains produced by CA-ID-TIMS analyses may result in slightly shorter crystallization ranges than expected (i.e.averaging effect; [ 22 , 68 , 78 ]), adding a degree of uncertainty in comparisons between time-resolved natural and synthetic zircon crystallization temperature histories ( Fig. 8 ).

Concluding remark
The method developed in this study allows users to simulate zircon crystallization trajectories in upper crustal magma reservoirs for a variety of thermal histories.The outputs of this model are tailored to allow a direct comparison with high-precision U-Pb zircon geochronology (CA-ID-TIMS) datasets.Specifically, the stochastic sampling of synthetic zircon age populations mimics the random process of selection of a limited number of zircons from a rock sample, while the uncertainty on the model output is parametrized based on a comprehensive dataset of high-precision zircon ages.The model is inclusive of statistical tools that allow the user to quantitatively compare simulation results and data from natural magmatic systems, thus informing on the dominant thermomechanical processes taking place in a crystallizing magma batch.The provided set of MATLAB R ○ codes contains the set of instructions necessary to run any set of simulations, as well as inputting new high-precision zircon datasets to be compared with the various model outputs.

Fig. 2 .
Fig. 2. Visualization of the stochastic zircon sampling workflow.Synthetic zircon age and temperature distributions are scaled from initiation (x = 1) to termination (x = 0) of zircon crystallization.Panels (a) and (d) represent cumulative zircon crystallization probability distribution functions    ( n * ) calculated from crystallization probability density function    (n) through Eq. (8) .Panels (b) and (e) represent the result of random sampling from    ( n * ) , where the 30 randomly selected crystallization steps (n) are shown as frequency histograms (upper boxes) and rank-order plots (lower boxes).In panel (c) and (f), the randomly sampled zircon crystallization distributions are displayed as kernel density estimates of the scaled simulation steps (n * ).Plots in (a), (b), (c) correspond to fractional crystallization (simulation ID: 1-fc from Table2) and (d), (e), (f) recharge fractional crystallization (simulation ID: 14-3d from Table2) simulations, respectively.

Fig. 3 .
Fig. 3. Effect of over-and under-sampling on the shape of synthetic zircon age spectra.Synthetic zircon age and temperature distributions are scaled from initiation (x = 1) to termination (x = 0) of zircon crystallization.Kernel density estimates (KDEs) of individual and averaged (bold lines) stochastic outputs from (a), (b), (c), (d) fractional crystallization (simulation ID: 1-fc from Table2) and (e), (f), (g), (h) recharge fractional crystallization (simulation ID: 11-5f from Table2) simulations.From left to right the number of stochastically sampled zircon crystallization steps (n) per simulation is increased from N Zr = 30 to N Zr = 10000.The random selection algorithm undergoes 10 3 iterations for each simulation.

Fig. 4 .
Fig. 4. Visualization of zircon crystallization temperature output from stochastic zircon sampling workflow.Individual temperature values that correspond to randomly selected zircon crystallization steps are shown in (a) for a fractional crystallization simulation (ID: 1-fc from Table2), and for (b) a recharge fractional crystallization simulation (ID: 14-3d from Table2).(c), (d) Box plots of temperature values from the same simulations.In panel (e) and (f), the randomly sampled zircon crystallization distributions are displayed as kernel density estimates of the scaled simulation steps (n * ).Synthetic zircon age and temperature distributions are scaled from initiation (x = 1) to termination (x = 0) of zircon crystallization.Box plots show median, interquartile ranges (IR) and extreme values for bins corresponding to 1/10 of the total zircon crystallization history.

Fig. 8 .
Fig. 8.Comparison of synthetic and natural zircon crystallization temperature distributions.Synthetic zircon crystallization temperature distribution patterns obtained from (a) fractional crystallization (simulation ID: 1-fc from Table2) and (b) recharge fraction crystallization (simulation ID: 3-1b from Table2) scenarios.Two contrasting natural zircon crystallization temperature distributions are also shown: (c) Bingham Canyon intrusion[ 40 ] and (d) Sesia Magmatic System[ 69 ].To allow direct comparison, natural and synthetic data are scaled from initiation (x = 1) to termination (x = 0) of zircon crystallization.Box plots show median, interquartile ranges (IR) and extreme values for bins corresponding to 1/10 of the total zircon crystallization history.

Table 2
Crystallization simulation scenarios (for compositions refer to Table1).