Zircon age spectra to quantify magma evolution

Zircon geochronology is the gold standard for defining the geologic timescale (e.g., Schaltegger et al., 2015; Schoene, 2014). The U-Th-Pb system in zircon, including the chronometer based on intermediate decay chain isotopes ²³⁰Th and its indirect parent, ²³⁸U, ranks nominally at the top of closure temperatures for geochronology, whereby the dissolution of zircon typically outpaces the diffusive mobilization of parent and daughter isotopes within crystalline zircon at reasonable temperatures and timescales for crustal magmatic systems (Cherniak and Watson, 2001). Hence, U-Th-Pb zircon dates can be straight forwardly interpreted as recording its crystallization in zircon-saturated melts (e.g., Watson and Harrison, 1983), although zircon can also form in metamorphic and hydrothermal environments (e.g., Kohn et al., 2015; Schaltegger, 2007). Three main zircon geochronology methods are used, and each has its individual strengths and limitations (Table 1). Isotopic analysis of U-Th-Pb in zircon by isotope dilution–thermal ionization mass spectrometry (ID-TIMS) is now achieved at <0.1% relative age uncertainty for single crystals or µgsized fragments (e.g., Schaltegger et al., 2015) because of the routine pretreatment of zircon by chemical abrasion (CA), which successfully mitigates the effects of Pb-loss due to metamictization (e.g., Kryza et al., 2012; Mattinson, 2005; von Quadt et al., 2014). These improvements were paralleled by progress in high spatial resolution analysis at the scale of subgrain domains (with <1 ng material consumed) using laser ablation–inductively coupled plasma–mass spectrometry (LA-ICP-MS), and at even higher spatial resolution with secondary ionization mass spectrometric (SIMS) techniques (Table 1). U-Pb zircon geochronology has a lower limit for precise age determination due to a combination of low radiogenic signals and required corrections for initial disequilibrium in the U decay chains (e.g., Schaltegger et al., 2015). Extending accessory mineral geochronology into the Holocene requires application of U-series methods, in which disequilibrium between ²³⁰Th and its indirect parent, ²³⁸U, is measured (U-Th geochronology). The absolute age resolution for the U-Th method, however, diminishes when approaching secular equilibrium, which effectively limits its applicability to zircon of <300 ka (e.g., Schmitt, 2011). Applying these methods, sometimes in combination, has increasingly helped to resolve complex crystallization histories and inheritance in zircon, whereby the highest absolute age precision is achievable in samples of young age and/or by the analytically most precise method (Table 1).

The general expectation of most geochronology users is that applying any of these different and often complementary methods to zircon will lead to a welldefined age—the equivalent of a point on the time axis (Fig. 1)—where assigned uncertainties only reflect the analytical and systematic uncertainties of the method (e.g., uncertainties in isotope ratio measurements, and decay constant errors, respectively). This age may then be interpreted as corresponding to a specific geologic event, such as the eruption of a volcanic rock (Fig. 1) or magma emplacement of a plutonic rock, with these events being thought to have happened faster than the timespan indicated by stated uncertainties. Uncertainties can be further reduced by averaging results from replicate analyses of individual crystals and/or spots from the same sample or closely related samples. Normally distributed ages from replicate analyses are generally vetted using the reduced chisquare statistics, with the mean square of weighted deviates (MSWD) expressing the degree of dispersion of analyses around an average relative to the magnitude of the analytical errors. Data sets in which the MSWD is within acceptable limits of unity are reasonably interpreted to be derived from the same normal distribution. In this case, the uncertainty for the age is the standard error, i.e., the standard deviation divided by the square root of the number of analyses. Averaging replicate analyses results in much improved precision compared to the individual ages (e.g., Crowley et al., 2007). Data sets with MSWD values >>1 indicate data dispersion exceeding what is explainable by analytical uncertainties alone. These are often interpreted as either affected by Pb-loss or zircon inheritance, which leads to younger and older ages in a population, respectively.

With increased analytical precision and improved spatial selectivity of available geochronological methods, however, zircon dates even from hand sample-sized volcanic or plutonic rocks often display a spectrum of ages rather than defining a point in time. Zircon age spectra thus result from dates being distributed over a time interval in which at least the extremes are distinctly resolved beyond analytical uncertainties. Age overdispersal relative to analytical uncertainties ranging between 103 years and >10⁶ years has been documented for both U-Th disequilibrium and U-Pb geochronometers at different timescales, and it is found regardless of the analytical method applied (TIMS, LA-ICP-MS, and SIMS) or the geological environment of the sample (e.g., for volcanic and plutonic zircon). Examples are detailed in section 2.6 below. The veracity of zircon age spectra is also upheld if stringent criteria for identifying potential Pb-loss or unrecognized inheritance (i.e., the incorporation of recycled zircon typically from wall rocks) are upheld, although both effects are certainly relevant for ancient zircon, or for some TIMS dates due to the larger crystal volume sampled using TIMS as compared to LA-ICP-MS or SIMS (e.g., Schaltegger et al., 2015). Furthermore, smoking-gun evidence against zircon age spectra being simply caused by unrecognized Pb-loss or inheritance comes from the U-Th zircon geochronology of late Pleistocene–Holocene volcanic rocks, which are too young to be affected by metamictization, and where U decay products, including ²³⁰Th, are essentially immobile in the crystal lattice (Cherniak et al., 1997). Thus, post-crystallization disturbance can be ruled out for young zircon crystals, and spatially selective analysis (e.g., through SIMS depth profiling) can distinguish growth domains at the micrometer scale that show significant rimward younging. Zircon growth rates (Table 1) have been determined by this method (e.g., Storm et al., 2012; Friedrichs et al., 2021).

Thus, it has become evident over the past 25 years or so that zircon crystallization in magmas must be conceptualized as prolonged rather than instantaneous, and this may significantly predate a volcanic eruption or emplacement of a shallow intrusion (Fig. 1; Reid et al., 1997). This has severe implications for interpreting the timing of a geological event such as a volcanic eruption from zircon age data, especially when compared to low closure temperature chronometers such as ⁴⁰Ar/³⁹Ar (Fig. 1; e.g., Oberli et al., 1990; Keller et al., 2018). In this regard, the resulting zircon age spectra in magmatic rocks are analogous to those found in sedimentary rocks, and for interpreting the geologic significance of these heterogeneous ages, igneous petrologists face problems similar to those of sedimentologists attempting to extract depositional ages and information on provenance from detrital zircon data (e.g., Coutts et al., 2019). Notable exceptions are (per-)alkaline magmas, where zircon solubility is extremely high (e.g., Gervasoni et al., 2016) and zircon can crystallize rapidly (and sometimes even as millimetersized megacrysts) within the resolution of U-Th and U-Pb chronometers (Table 1). Such rocks display comparatively uniform zircon age distributions in individual samples, but not necessarily for different samples from the same volcano (e.g., Campi Flegrei: Fedele et al., 2006; Gebauer et al., 2014; Laacher See: Schmitt et al., 2010, 2017). Other cases where zircon crystallizes rapidly involve hightemperature, mafic- to- intermediate magma compositions, where eruption temperatures are close to zircon saturation conditions (e.g., in mid-oceanic-ridge settings; Hayman et al., 2019; Schmitt et al., 2011). Nonetheless, even such zircon age data sets, when viewed system-wide instead of at the scale of an individual hand-sample, can be complex and indicative of protracted magmatic evolution (e.g., Schmitt et al., 2010).

Here, we review challenges in zircon geochronology to reliably interpret zircon age spectra in the context of the different methods by which they were generated. As our focus is primarily on magmatic systems, we largely disregard the possibility of subsolidus (“hydrothermal”) zircon crystallization in this review (e.g., Troch et al., 2018) and refrain from any discussion of zircon crystallization during regional metamorphism; instead, we refer readers to excellent recent reviews (Kohn, 2016; Rubatto, 2017). Because zircon crystallization in magmas is directly linked to the thermochemical conditions of the melt that are faithfully recorded by zircon (e.g., Claiborne et al., 2006, 2010), zircon age spectra systematics, when properly understood, can reveal processes of magma migration, accumulation, and cooling timescales for natural systems at timescales unattainable by other methods (Fig. 1). Calvin Miller, to whom this themed issue is dedicated, pioneered the investigation of accessory minerals to gain insights into magmatic processes that operate at scales ranging from individual volcanoes or plutons to planetary crusts. In this tradition, zircon age spectra are increasingly harnessed to constrain the thermal history of magmatic systems and the longterm magma fluxes that sustain them (e.g., Caricchi et al., 2014, 2016; Ratschbacher et al., 2018; Friedrichs et al., 2021; Liu et al., 2021; Lukács et al., 2021; Tierney et al., 2016; Weber et al., 2020; Melnik et al., 2021; Angeles-De La Torre et al., 2023). Used in this sense, they are complementary to crystalscale diffusion chronometers that provide magmatic timescale information if diffusion parameters are known, and which—unlike zircon—also can be applied to mafic rocks where zircon is commonly absent, but with absolute time information (Fig. 1; see Costa et al., 2020, for a recent review).

Zircon crystallizes primarily upon cooling and differentiation of magma as first experimentally quantified by Watson and Harrison (1983). Because magmatic systems evolve dynamically, with magma reservoirs undergoing episodes of recharge and eruptive venting (e.g., Annen et al., 2006), thermal variability in these reservoirs over time and space is expected. At low temperature and in subalkaline melts, the kinetics of zircon dissolution are sluggish; shortlived prograde pulses and the associated chemical changes resulting from magma recharge only cause incomplete resorption (Watson, 1996). The propensity of zircon to survive such heating events is underscored by the preservation of rounded crystals in arc-volcanic rocks with bulk compositions and eruptive temperatures at which zircon is undersaturated (e.g., Weber et al., 2020). If not erupted, preexisting zircon, albeit potentially partially resorbed, can become the nucleus for renewed zircon crystallization during retrograde episodes (Szymanowski et al., 2020; Bindeman and Melnik, 2016). The shape of the natural distribution of zircon ages in individual samples thus can be highly variable, and it also will depend on the method applied to determine zircon ages (see section 3). Different approaches to visualizing zircon ages, especially those that are over-dispersed relative to analytical uncertainty, have been developed, and these are described in the next paragraph and sections 2.2–2.4.

Most intuitively, ages can be displayed in a histogram with individual analytical results assigned to discrete bins. The appeal of a histogram is its simplicity (the shape is only affected by the bin width) and the faithful reproduction of the data with a low risk of over-interpretation. However, the use of a single bin width (e.g., equivalent to the analytical uncertainty) for an entire distribution is precisely what makes data visualization difficult for large data sets that cover a comparatively large age range relative to analytical uncertainties. Moreover, analytical uncertainties for young zircon are primarily controlled by low signal intensities, unlike older zircon crystals, where other sources of uncertainty affecting the results dominate, such as correction factors for isotopic mass fractionation in ID-TIMS or relative sensitivities in LA-ICP-MS and SIMS (e.g., Schaltegger et al., 2015; Schoene, 2014). Hence, uncertainties for young high-U zircon are often much smaller than those for low-U zircon. One should note that as the number of analyses, n, approaches infinity in a data set, the shape of a histogram converges with that of a continuous curve. However, time and cost are limiting factors for generating such high-n data sets; so, for most samples, an approach is required where an acceptable approximation can be extrapolated from a limited number of dates (e.g., n = 20–100). The generation of such an approximation can be achieved in multiple ways, and the resulting representation of zircon age spectra is dependent on several key assumptions (see section 2.5).

PDPs are another straightforward way to generate age spectra in addition to histograms. The age and uncertainty of each measurement provides the mean and standard deviation of a Gaussian distribution, and the sum of these overlapping Gaussian distributions generates the age spectrum; alternatively, the probabilities can be summed into a cumulative PDP. The utility of this approach, however, has been called into question if analytical uncertainties are invariable and small compared to the range of values in the data set. Vermeesch (2012) noted that PDPs can “punish” less precise measurements by under-representing such measurements in detrital spectra. Namely, the sampling of a zircon population with higher inherent error (e.g., due to lower U concentration) will generate a broader Gaussian distribution. Then, the PDP of this distribution will be correspondingly smoothed due to the higher uncertainties of the individual data points. For these reasons, PDPs are sometimes foregone in favor of kernel density estimates (KDEs).

KDEs are similar to PDPs in basic principle: Each measurement is arranged into a distribution (referred to here as a “kernel,” for which different shapes can be prescribed, including Gaussian), where the width of the distribution is not determined by the analytical uncertainty (as in PDPs), but by a “bandwidth” that is algorithmically adjusted based on point density and therefore variable across the age range. As noted by Vermeesch (2012), KDEs have advantages for provenance analysis using detrital zircon age spectra where analytical uncertainties and ages vary considerably.

Drawbacks of KDEs are the complexity of the formulas used and the potential for the arbitrary selection of bandwidths. Whereas histograms and PDPs are easily produced with minimal input parameters and subsequently easily reproduced by other researchers, KDE formulas are typically hidden from view and can vary in how users of the program (e.g., DensityPlotter; Vermeesch, 2012) adjust the bandwidth. For example, consider a synthetic bimodal data set with two prescribed peaks at 5 Ma and 7 Ma (7000 ages; Fig. 2A). A 250-point random sample of the population is selected, representing the collection of data from a much larger population. The algorithm of Botev et al. (2010) yields multiple peaks that actually represent noise in the original data set (Fig. 2B). The estimate of DensityPlotter (Vermeesch, 2012), with an adaptive kernel bandwidth, avoids such overinterpretation and faithfully reproduces the original population shape, which indicates that this is the algorithm of choice for KDEs (Fig. 2C). However, the most notable results of this demonstration are that (1) KDEs can be strongly model-dependent; (2) it is essential to report metadata with KDEs (e.g., constant versus adaptive bandwidths); and (3) the algorithm of Vermeesch (2012) appears to be suitable for modeling zircon age spectra with realistic uncertainties, provided that no major changes are introduced to the source code.

Another method of data visualization that is currently underutilized in U-Th-Pb zircon geochronology is the radial plot, which was introduced by Galbraith (1990) to represent fission-track data. Here, an additional dimension is added, allowing the user to visualize the amount of relative uncertainty on each data point on the x-axis, as well as the distribution of ages on the y-axis (Figs. 3B and 3C). Individual age peaks are calculated, as with KDEs (Fig. 3B), with the bonus that xenocrystic populations that may be difficult to spot in a histogram appear alongside the rest of the data at a lower point density (see Miller et al., 2007, for definitions of antecrystic and xenocrystic as used here). A modified version of this plot, the abanico plot (named after the Spanish term for a folding fan), includes a PDP on the right side of the graph (Vermeesch, 2012), combining properties of Figs. 3A and 3B.

The advantages of the radial plot being applied to magmatic zircon populations are, firstly, that the consideration of analytical uncertainty immediately informs the user about the relative value of each data point, decreasing the risk of overinterpretation even at a glance. Secondly, the existence of many zircon analyses of similar age but varying uncertainty works to create vectors in the radial plots that are easily differentiated from one another, and the slope of these vectors anchored at the origin indicates the age (marked along the radial part of the plot). Lastly, the addition of another data visualization method increases the confidence in the user’s interpretation, particularly because the interpretation of age spectra is only a semiquantitative method that requires a careful qualitative assessment as opposed to strict quantitation.

The attractiveness of the radial/abanico plot is that it easily and clearly captures the entire age spectrum, making it applicable to volcanic and plutonic zircon suites that may have a large xenocrystic component. This type of plot additionally excels at finding peaks of high sample density, which is particularly applicable to LA-ICP-MS/SIMS data sets. Application to TIMS data sets, on the other hand, is complicated by their relatively small size (typically <20 zircons) and high precision, which leads to the accentuation of individual zircon ages in a radial plot as opposed to a representative age of abundant zircon crystallization.

The central supposition behind zircon age spectra analysis is that practical sampling of the zircon population (which is typically limited to tens to hundreds of zircons) still allows generation of an accurate representation of the true distribution. This is, however, far from trivial, because selection bias due to the method of zircon extraction (which may hinder, e.g., crystals included in phenocrysts from being liberated) or picking (where typically large, euhedral crystals are preferred) can hardly be avoided. With this caveat, we address the more tangible statistical problem of how many of the crystals that are available should be analyzed after a rock is subjected to a certain zircon identification and extraction procedure (Table 1). The larger the data set, the more accurate the representation of the true spectrum, but in practice, there are upper limits on the number that can be analyzed given time and cost considerations. The optimal size of n depends on (1) how many modes there are in a spectrum, as well as (2) the separation of the peaks and the analytical uncertainty, and (3) the model used to generate the spectral representation.

The number of modes in a zircon age distribution is effectively unconstrained, and the shape of a natural population will theoretically always appear spiky, provided there is sufficient temporal resolution (i.e., there could be as many age peaks as zircon crystals or crystal subdomains analyzed). One must therefore make a rough judgment about the number of crystals required to accurately capture the variability of their data set. For provenance analysis in sedimentary records (where age ranges may span billions of years), it is necessary to look at data sets numbering 100–1000+ zircons to resolve low-abundance age populations (Dodson et al., 1988; Pullen et al., 2014; Vermeesch, 2004). For magmatic studies, where the age range is much smaller, a data set of 100 zircons may encompass the natural variability of the population. A simple subsampling experiment demonstrates this concept (Fig. 4): Using subsamples of 20, 50, and 100 zircons from a data set of Miocene volcanic zircon crystals from Pannonia (Lukács et al., 2018), one can see that the overall variability of the entire data set is adequately captured by subsampling 50–100 zircons, which is in line with estimates from Caricchi et al. (2016). Here, the larger subsample also begins to resolve input from xenocrystic zircons at ca. 25 Ma as well as a shoulder in the data at ca. 18–19 Ma, which may reflect another population too close in age to be resolved by LA-ICP-MS. Importantly, small data sets (10–20 zircons) demonstrate no bimodality, which suggests that the KDE adequately protects against over-assigning individual peaks. For such small data sets, sampling only suggests an age of roughly 15–17 Ma (Figs. 4A and 4B). In larger subsamples, however, one can more confidently argue for the existence of two age peaks and perhaps identify different proportions of xenocrysts and the existence of subsidiary age peaks (Figs. 4C and 4D). However, these latter two interpretations should be considered strictly qualitative, and they require corroborating evidence. A statistically robust approach for identifying distinct populations in such mixtures was developed by Sambridge and Compston (1994) that is based on deconvolving multiple Gaussian peaks from a polymodal PDP.

The distance between age peaks (Δt) and analytical uncertainties plays a large role in determining whether age peaks can be resolved. Using the convention of a 95% confidence interval (abbreviated as 2σ, where σ refers to the standard error of an individual analysis; Table 1), it follows that for two age peaks to overlap at less than 5%, their separation must be at least 4σ (2σ for each peak). This is consistent with the Bayesian modeling of Keller et al. (2018), who determined that the ideal separation of age peaks should fall within the 3–5 σ range (i.e., Δt/σ = 3–5). In the previous example, with large age peaks at 14.8 Ma and 17.2 Ma and an analytical uncertainty amounting to ca. 0.7 Ma at the 4σ level (the mean of the 3–5σ range identified by Keller et al., 2018), the corresponding Δt of 2.4 Ma translates into a Δt/σ of ~14. The apparent shoulder at ca. 18.5 Ma presents a conundrum, however. While on the one hand, the Δt between this “peak” and the one at ca. 17.2 Ma is 1.3 Ma (well outside the 4σ criterion), the distribution appears to be continuous toward older ages, which suggests that this period experienced prolonged zircon crystallization and may in fact be composed of many smaller, and ultimately unresolvable, peaks. The test applied here suggests that the minimal sample size for resolving bimodality in the Pannonian zircon example is ~50 zircon analyses, but reliable estimates of the percentage of each mode may require a larger number of analyses (e.g., 100; Fig. 4D). In the data set used here, the range between the modes is sufficiently large, and analytical uncertainties are sufficiently low. Under other circumstances, these criteria might not be met, and it is difficult to define a set of rules that applies universally, apart from two basic principles: (1) Δt/σ must be at least 3–5, and (2) a consistent algorithm that is commonly available must be used (e.g., Vermeesch, 2012).

Slightly different considerations hold for U-Th zircon ages, which cover an age range between ca. 0 ka and 300–350 ka. Because of the exponential relationship between isochron slope m (which for single zircon ages is usually defined by two points: the measured zircon composition and a constant model melt composition that is typically derived from U-Th wholerock or glass analyses) and time, the absolute age uncertainty increases as secular equilibrium is approached. Kent and Cooper (2018) concluded that >100 zircon analyses are required to characterize U-Th zircon age spectra and to correctly interpret them in the context of a magma system’s thermal history, but as ages approach secular equilibrium, the absolute age uncertainty, even for the same analytical precision, becomes too large to meaningfully distinguish different age populations. In favorable cases (e.g., high-radiogenic Pb yields), U-Pb geochronology can seamlessly complement U-Th ages.

Obtaining large zircon data sets is not always practical, especially for ID-TIMS geochronology, and we therefore showcase welldocumented examples where zircon age spectra have been obtained at high-n sampling. To illustrate, we have summarized in Figs 5–8 PDPs, KDEs, and radial plots from four data sets of U-Pb and U-Th zircon data covering a wide age range. Data have been collected over the past few years at the University of California at Los Angeles/Heidelberg University (UCLA-HD) and ETH Zürich using SIMS and LA-ICP-MS methods, respectively, and include: (1) The Belfond lava dome from the late Pleistocene Soufrière Volcanic Center (Saint Lucia, Lesser Antilles Arc); (2) the Miocene plutonic complex of La Gloria (Central Andes); (3) the Oligocene volcanic province of the Southern Rocky Mountains volcanic field (also previously called the San Juan volcanic field), western USA; and (4) the Ordovician plutonic complex of the Sierra de Famatina–Valle Fértil in NW Argentina.

The Soufrière Volcanic Center data (Fig. 5) comprise zircon from samples of the dacitic Belfond lava dome and cogenetic granodioritic enclaves that erupted at ca. 13.6 ka (Barboni et al., 2016; Schmitt et al., 2010). The data set includes zircon rim analyses of pristine growth surfaces, as well as interior analyses of sectioned crystals, but we do not differentiate between those in this comparison (cf. Schmitt et al., 2010). Lateral overlap of the ~30-µm-diameter ion beam used in SIMS analysis for crystal interior analyses could lead to mixed ages, whereas this is largely mitigated for the rim analyses in which the SIMS sputtering is parallel to the growth direction, integrating to a depth of <5 µm. Most of the Belfond zircon crystallization ages significantly predate the eruption age (dashed lines in Fig. 5) which was independently constrained by (U-Th)/He zircon dating. Interestingly, the dominant zircon age mode at ca. 28–29 ka is similar in the zircon populations of both the host lava and enclosed plutonic enclaves, which suggests recycling of zircon from crystalrich domains of an intrusive reservoir during a younger magmatic episode (Fig. 5). These age modes predate the dome eruption by ~14 k.y. and the last major explosive eruption in the region at ca. 20 ka by ~8 k.y. Several minor peaks in decreasing abundance are present between ca. 40 ka and 114 ka, as well as some older crystals, in particular in the plutonic enclaves. The central age of the population (ca. 62 ka and 49 ka for lava and plutonic enclaves, respectively) has no geologic significance, but the significance of the older zircon crystals in terms of continuous versus punctuated crystallization is more difficult to distinguish, especially when the ion beam might have overlapped different growth domains (Fig. 5). The predominance of ages between ca. 40 ka and 114 ka, however, corresponds well with major eruption pulses in the Soufrière Volcanic Center and a possible hiatus between ca. 250 ka and 110 ka, which again supports crystal recycling in a longlived magma system (Schmitt et al., 2010).

La Gloria U-Pb zircon LA-ICP-MS data show a nearly Gaussian distribution, with a single peak at ca. 10.4 Ma (Fig. 5), which agrees within error with a published U-Pb zircon LA-ICP-MS age of 10.3 ± 0.2 Ma (Deckart et al., 2010). However, higherprecision ID-TIMS data on the same samples establish that the pluton crystallized zircons for almost 2 m.y., between 10.2 Ma and older than 12.0 Ma (Gutiérrez et al., 2018). Hence, it would be erroneous to interpret the relatively simple “unimodal” population seen by the LA-ICP-MS data as all zircon crystallizing during a single event corresponding to the weighted mean age of ca. 10.4 Ma. That such an interpretation is inadequate is also underscored by the elevated MSWD for the entire data set.

At even larger spatial and temporal scales than for individual volcanoes or plutons, zircon age spectra are also found for entire volcanic fields or arc segments. Geological background and additional age information is presented in the references provided, and here, we only summarize zircon ages from these cited studies. Large zircon age data sets for the Southern Rocky Mountains volcanic field (Sliwinski et al., 2022) and the Famatina Province show distinct peaks (Figs. 7 and 8). Again, both data sets display elevated MSWD values and central ages that lack any significance for dating a specific geological event. For the Southern Rocky Mountains volcanic field, two main peaks stand out at ca. 33–34 Ma and ca. 27–29 Ma that correspond to the peaks of activity in the northern and southern parts of the volcanic province, respectively. There is also a younger peak at ca. 23 Ma, which corresponds to the beginning of the Rio Grande Rift in the region, and the emplacement of some large silicic units (in particular, the Sunshine Peak Tuff that erupted at 22.93 ± 0.02 Ma; Lubbers et al., 2020; Mehnert et al., 1973).

For the Famatina region, there is a broad main peak, centered on an age of ca. 463 Ma, and a younger tail (<350 Ma) that may represent an episode of Carboniferous magmatism. However, as both the Southern Rocky Mountains volcanic field and the Famatina regions are relatively old, the LA-ICP-MS uncertainties likely mask more complexity within the zircon age population. Based on highprecision ID-TIMS data, we know that zircon from single ignimbrite units of the Southern Rocky Mountains volcanic field has crystallization intervals on the order of 0.5–1 m.y., with potential zircon populations hidden within those peaks (e.g., Curry et al., 2021). Moreover, zircon ages systematically correlate with geochemical indicators such as Yb/Dy, which reflects accessory mineral saturation and hence the temperature evolution of the magma (Bachmann et al., 2007; Wotzlaw et al., 2013). Similarly, preliminary unpublished ID-TIMS data on some of the Famatinian units suggest peaks of zircon growth at ca. 475–480 Ma and ca. 465–470 Ma, hidden in the broad Ordovician peak of the comparatively low-precision LA-ICP-MS data. Hence, as pointed out in section 2.5 above, distinguishable events must be separated by >3–5σ to be isolated. With typical errors of ~1–2% (at 2σ) for U-Pb zircon LA-ICP-MS data, this means that igneous pulses for San Juan and Famatina need to be separated by more than 1–1.5 m.y. and 10–15 m.y., respectively, to be statistically distinguishable.

However, not every zircon age data set where excess scatter beyond analytical uncertainties exists (commonly inferred from an elevated MSWD value) should be considered a true zircon age spectrum. Sampling bias and analytical problems can produce excess scatter, and common issues are often specific to the analytical method applied. For example, even with chemical abrasion pretreatment, zircon domains that experienced Pb-loss may still contribute to an ID-TIMS analysis (e.g., Widmann et al., 2019). Unrecognized inheritance, especially when xenocrystic contaminants are only marginally older than the magmatic episode, can also shift ID-TIMS dates to older ages (e.g., Widmann et al., 2019). Similar considerations hold for data from spatially selective methods, where CA pretreatment also shows promise but may not always be the silver bullet to eliminate zircon affected by Pb-loss (e.g., Kryza et al., 2012; Sliwinski et al., 2017; von Quadt et al., 2014).

Bias can also affect spatially selective methods: LA-ICP-MS dates are sensitive to fractionation corrections that depend on the depth of the ablation crater, and because ablation characteristics vary depending on the structural state of zircon, this may result in apparent age variability (Marillo-Sialer et al., 2016; Sliwinski et al., 2017). SIMS analysis of high-U zircon is known to yield apparent older ages when the relative sensitivity factor (RSF) is calibrated on normal-U zircon references (e.g., Schaltegger et al., 2015; White and Ireland, 2012). Lastly, assigning errors is modeldependent in all analytical methods, and not all sources of error are always identifiable. Thus, it is essential to monitor analytical performance by analyzing matching secondary zircon references to identify potential under- or over-estimation of analytical error. Thus, we urge caution in invoking “protracted zircon crystallization” as a carte blanche for explaining away excess analytical scatter, particularly for geologically old igneous rocks, because the magmatic lifetimes typically observed for young volcanic systems are but a minor fraction of the overall analytical uncertainty. Analytical reliability, e.g., from analyzing homogeneous and well-characterized zircon references, is a prerequisite for identifying true zircon age spectra.

In the second part of this review, we focus on existing and upcoming forward models to generate distributions of zircon crystallization ages with the goal of quantifying rates and durations of magmatic input into volcanic plumbing systems (e.g., Caricchi et al., 2014, 2016; Friedrichs et al., 2021; Liu et al., 2021; Lukács et al., 2021; Tierney et al., 2016; Ratschbacher et al., 2018; Weber et al., 2020; Melnik et al., 2021). The main idea behind this modeling approach is to identify for which set of input parameters a computed zircon distribution—including age and possibly crystal age zonation as well as chemistry—is the closest to measured zircon from volcanic or plutonic rocks. These models are based on five tenets: (1) the thermal history of a magmatic intrusion can be adequately predicted, with reasonable assumptions regarding the shape of an igneous body (e.g., Annen et al., 2015); (2) zircon ages record magmatic crystallization; (3) zircon crystallizes between the zircon saturation temperature (T_zr) and the solidus temperature (T_s), and the zircon fraction per unit mass of magma decreases from T_zr to T_s according to an empirical zircon saturation model (e.g., Baker et al., 2002; Boehnke et al., 2013; Borisov and Aranovich, 2019; Gervasoni et al., 2016; Shao et al., 2020; Watson and Harrison, 1983); (4) the number of zircon crystals crystallizing in a magmatic plumbing system at each time interval is proportional to the change in the temperature distribution between T_zr and T_s; and (5) the fraction of measured zircon of different ages is representative of the proportions of zircon that crystallized within each age interval.

This last assumption deserves further comment, as the age distribution ultimately measured depends on the type of sample investigated (i.e., plutonic or volcanic) and the analytical method used (Table 1). Firstly, volcanic and plutonic zircons will sample different volumes and evolutionary stages of a magmatic system as plutonic zircons can grow in areas of very high crystallinity (and possibly hydrothermally at subsolidus conditions), whereas rheological thresholds for magma being able to erupt require that volcanic zircon is largely derived from the liquid-dominated part of a system at the time of eruption. As a telling example, basalt is typically devoid of zircons, but many gabbros contain zircon that crystallized interstitially from residual melt at near-solidus conditions (e.g., Hayman et al., 2019). Secondly, the analytical technique used for dating can introduce bias (i.e., high spatial resolution on sectioned crystals in SIMS and LA-ICP-MS versus whole zircon by TIMS; Table 1). We first introduce these general concepts behind the modeling and then discuss the robustness of the critical assumptions in the following section (3.2).

To explain the fundamentals of the modeling approach, we can consider the simple case of a spherical magmatic intrusion cooling at depth (Fig. 9) before addressing more complicated but realistic scenarios in which a magma system undergoes recharge. In the simplest case, magma is injected instantaneously at its liquidus temperature, (T₀), and the intrusion cools starting from its outer portions. Zircon crystallizes between T_zr down to T_s, and the rate of zircon crystallization either decreases linearly or (more realistically) exponentially with temperature (Fig. 9A). The exponentially decaying mass fraction of zircon is consistent with zircon saturation calibrations, where ~40% of the total amount of zircon that can form in a batch of magma crystallizes within 25 °C below the onset of zircon saturation, and only <10% crystallizes over the last 25 °C before reaching the solidus temperature (Harrison et al., 2007). Several observations can be made from this simple scenario. The distributions of ages calculated assuming linear or nonlinear (exponential) relationships between T and the amount of zircon crystallized differ, with linear models slightly overrepresenting younger zircon crystals compared to the nonlinear (exponential) model (Fig. 9B). In the absence of dynamic processes such as convection, the age distributions in the different portions of a solidified intrusion would be different, and the age of the youngest measurable zircon would decrease from the outer to the inner portions of the intrusion (Fig. 9B). If, hypothetically, the different shells were dynamically mixed prior to eruption (e.g., due to recharge and remobilization; Bachmann and Bergantz, 2006; Huber et al., 2009), a volcanic rock could contain a mix of all ages and different compositions (e.g., Pamukçu et al., 2022).

Considering a geologically more plausible system assembled by the injection of discrete pulses (e.g., Coleman et al., 2004), a zircon age distribution that is noticeably different from the previous scenario must be expected (Fig. 10). Importantly, the calculated distribution of measurable zircon ages would also look different if instead of considering the time-temperature paths of different portions of the intrusion, we were to only track the evolution of the average temperature. This implies that zircon age distributions calculated using an “average” time-temperature path cannot be used to retrieve reliable information about the thermal evolution of magmatic systems assembled by episodic magma injection (Figs. 10A and 10B).

The type of sample under investigation is also relevant, as zircon age distributions from plutonic or volcanic rocks will not have the same characteristics. For instance, if an eruption occurs during the assembly of a crustal magma reservoir, zircon sampled by the eruption will record a shorter age range with respect to zircon collected from the solidified magma reservoir (Fig. 10B versus Fig. 10C). An exception could result from the remelting of a completely solidified portion of the reservoir, thus releasing its zircon cargo to the melt. The zircon age distributions in plutonic and volcanic rocks are different because an eruption must occur before the full solidification of a magmatic reservoir, and therefore the youngest measurable age in volcanic zircons must be older than the youngest measurable age in plutonic zircon. Additionally, some of the injected magma continuously cools below the solidus during assembly of a magmatic reservoir (e.g., Annen, 2009; Annen et al., 2006, 2015; Karakas et al., 2017), and zircon locked within the solidified portion of the magmatic system is less likely to be sampled by an eruption (e.g., Shane et al., 2012). Finally, zircon ages determined by bulk analysis techniques such as TIMS are average ages that record integrated down-temperature zircon crystallization from a zircon-saturated melt that may have lasted longer than analytical uncertainties indicate (Table 1). High spatial resolution ages from LA-ICP-MS and SIMS analyses, by contrast, permit the resolution of different stages of the magma’s thermal history at 10³–10⁴ year intervals (Table 1). The example introduced earlier in this section of a single batch of magma cooling in a thermally uniform crustal environment is of course unrealistically simplified, as longlived magmatic systems are sustained by frequent magmatic recharge. Thermal models for such systems that are open for recharge exist, but different approaches are rarely vetted against one another. Two of these models that are in use at the University of Geneva and UCLA-HD are therefore described here in detail to assess their similarities and differences regarding the inferred magmatic thermal history and to compare their predictions for zircon age spectra, primarily in volcanic settings.

There has been a long tradition of modeling the thermal evolution of magmatic systems following the repeated intrusion of dike- or sill-shaped magmatic batches. In these models, dikes or sills are kinematically inserted into a model space while moving host rocks to accommodate magmatic input, and the subsequent thermal evolution of the crust is monitored while taking partial melting and effects such as latent heat into account (e.g., Dufek and Bergantz, 2005; Annen et al., 2006; Caricchi et al., 2014; Tierney et al., 2016; Ratschbacher et al., 2018; Weber et al., 2020; Melnik et al., 2021). All models show that the first injections cool rapidly and a certain time is required to build up a partially molten system. Yet, the details of the different models vary, and therefore results of three numerical codes for two end-member cases are compared.

One model that has been successfully applied for the forward modeling of zircon age spectra is the finite elements (FE) approach of Caricchi et al. (2014, 2016), in which thermal diffusion was solved using an FE approach with periodic sill insertion and passive tracers recording the temperature evolution over time for different portions of the magmatic system. A subsequent version of this model used the finite difference (FD) method to discretize the diffusion equation, whereby the temperature evolution is recorded in each numerical cell that contains magma above the solidus (Weber et al., 2020). These Geneva models track the thermal evolution of a magmatic system assembled by the injection of sills (Fig. 11A), which was previously shown to be chiefly controlled by the thermal state of the crust and the vertical rate of magma accumulation (Annen et al., 2015; Caricchi et al., 2014). Variations of the mode of magmatic injection and geometry of the magmatic system (either randomly in a prescribed volume, or by underaccretion, in which sills are stacked in a downward growth direction) can be modified. These parameters only play a minor role in the longterm thermal evolution of magmatic systems and therefore in the calculated zircon age distributions, but they can produce strongly divergent thermal histories during the early stages of a system (Annen et al., 2006, 2015; Carrigan, 1988; Karakas et al., 2017).

In the Geneva family of models (Fig. 11A), an axisymmetric geometry is assumed, magma is injected at a set temperature (e.g., 1000 °C), and the magma crystallinity (x_f) increases with decreasing temperature (T in °C), as:

where

Alternatively, the relation between T and x_f was estimated from a fourth and fifth order polynomial fit to the experimental melt fraction temperature curves of Blatter and Carmichael (2001) and Marxer and Ulmer (2019). The time-temperature (t-T) trajectories of the system are then derived from a large number (e.g., 1000) of passive tracers sampling temperature at constant time intervals between the onset and the end of magmatic injections. Each tracer can produce an arbitrary number of measurable ages (e.g., 100), and the rate of zircon crystallization decreases with decreasing temperature following the derivative of Equation 3 for T (Tierney et al., 2016):

where the cumulative fraction of zircon (F_zr) for each unit of magma increases with decreasing T, and zr_m is the maximum fraction of zircon crystallizing between zircon saturation (T_zr) and the solidus temperature (T_sol). At the time when the zircon population is sampled (e.g., in a volcanic eruption), only tracers with a temperature higher than a minimum temperature (T_min) are considered. T_min can represent either the minimum temperature at which magma (and zircon therein) can erupt, or the temperature of the interstitial melt extracted (together with zircon) to feed a volcanic eruption. The tracer that spends the longest period within the temperature interval of zircon saturation provides the maximum possible duration of continuous zircon crystallization. Once this is determined, all measurable ages between the oldest age and the time of eruption are summed as schematically shown in Figure 10.

The second type of model (UCLA-HD; Friedrichs et al., 2021; Tierney et al., 2016) starts with a two-dimensional grid of 0.1 × 0.1 km cells over a 20 × 60 km (width × depth) rectangle that represents a mid–upper crustal block. Parameters for modeling the thermal and compositional evolution of magma and its surroundings within this grid are based on equations for recharge-assimilation–fractional crystallization (RAFC; Spera and Bohrson, 2001). Geometrically, the UCLA-HD model differs from the Geneva models in that instead of vertically stacking sills from a prescribed depth and typically displacing the crustal rocks downward (e.g., Weber et al., 2020), magmatic injection always occurs in the center of an ellipsoidal volume of magma, i.e., the hottest and therefore likely most longlived standing volume of magma in the intrusive body (Fig. 11B). Upon growth due to periodic magmatic recharge, the body maintains its originally defined aspect ratio. The thermal evolution of individual cells is monitored throughout the simulation as they migrate away from the center of magmatic recharge and are translated into a pseudo-3–D distribution assuming rotational symmetry around the vertical central axis of the magmatic body in the center of the grid.

The crystal fraction, x_f, is calculated somewhat differently from that of the Geneva model (Caricchi et al., 2014, 2016) by using the following relationship:

where T_liq and T_sol are the liquidus and solidus temperatures of the recharge magma, respectively. The temperature of magmatic recharge is selected based on reasonable geological constraints (e.g., based on thermometry of coerupted mafic rocks; Friedrichs et al., 2021). At each timestep during the model run, the temperature change of each cell is monitored: For each cell undergoing cooling within the saturation range (e.g., 800–700 °C), the fraction of zircon growth is then calculated using Equation 3. This results in a time-resolved zircon growth curve for each cell, and each curve is then summed at the time of eruption. Only cells that are above 700 °C (supersolidus) are considered and weighted according to their relative volumes when translating the 2-D grid into 3-D. These zircon growth curves distinguish between core (e.g., at 20% of each curve) or rim ages (e.g., at 80%; Friedrichs et al., 2021). An average total zircon age can also be obtained from the integration of the entire zircon growth curve.

A limitation of existing numerical modeling studies that simulate the thermal evolution of magmatic systems is that the source code is typically not readily accessible in publications, and the assumptions behind the models often vary according to a particular geological setting, which makes comparison of the different studies difficult. Moreover, numerical implementation details, such as how to deal with nonlinearities, are beyond the limits of typical publications, which limits reproducing results. Therefore, we developed a new Julia package to simulate the thermal evolution of magmatic systems (MagmaThermoKinematics.jl, or briefly MTK). The new package is open source and works in 2-D, 2-D axisymmetric, and 3-D geometries. It utilizes modern hardware such as graphics processing units (GPUs) and multicore processors and provides a simple yet extendable way to change material parameters (see the Appendix for a detailed description). Several dike injection algorithms are implemented, including one that employs an analytical solution to insert coinshaped cracks into an elastic halfspace, as well as central injection and sill under-accretion. Tracers are injected within newly added dikes and record the subsequent thermal evolution from which zircon age distribution can be calculated. We have benchmarked this third modeling approach versus the two end-member scenarios shown in Fig. 11, which demonstrate that differences in the thermal structure are generally <5 °C at the end of the simulation (see Fig. S1 in the Supplemental Material¹).

Overall, both model geometries define plausible end-members of multicyclic intrusions, where the Geneva sill geometry leads to initially more rapid heat dissipation (especially when the depth of sill emplacement is random, rather than vertically stacked) compared to the UCLA-HD model, where the ellipsoid has a smaller surface- to- volume ratio than a thin cylindrical sill, and recharge always occurs in the hottest part of the system. Crustal thermal conditions prior to intrusion are defined in both models as constant mantle and surface heat flow. Country-rock partial melting and assimilation (Spera and Bohrson, 2001), as well as radioactive heating, are irrelevant for comparatively brief model-run durations and the shallow depth of sub-volcanic igneous bodies and their low country-rock temperatures.

To understand how differences in the choice of material parameters and boundary conditions affect the thermal structure, we systematically varied many of the parameters for the two scenarios. Results show that the central intrusion setup is very robust with respect to variations in the material parameters, and that the sill under-accretion case has larger variations (Fig. S1; see ^{footnote 1}). This difference is mostly likely caused by the sill under-accretion case not having reached steady-state yet, even over 1.5 m.y. run durations, whereas this is the case for the central intrusion case.

In principle, knowing the thermal conditions in an intrusion at any time allows calculation of the number of possible zircon crystallization ages potentially available for measurement. The calculated number of measurable ages can then be converted into a PDP (normal or cumulative) for comparison with the distribution of measured ages (Fig. 12). Alternatively, a simple parameter, such as the age difference (Δt) between stated percentiles of the modeled zircon age distribution, can be used for comparison (e.g., using the 16–84 percentile difference, corresponding to ± 1 s.d., although the resulting distributions in both types of models are strongly skewed toward younger ages; Fig. 12). In the model-scenario examples (Fig. 12), Δt systematically increases with the duration of continuous melt presence and recharge rate. The input parameters of the model that best match the measured distribution of zircon ages can then provide a quantitative estimate of processes operating in a longlived crustal intrusive body, such as the rate of magma recharge. However, this would be true only if dating could be performed in each zone for a large number of zircon crystals. Instead, actual data sets are necessarily limited in number and resolution, and it is thus important to consider how to compare the modeled and natural distribution of zircon ages measured with ID-TIMS or high spatial resolution techniques.

Before discussing applications of these models, we first focus on the analysis of zircon collected in plutons and discuss the modifications to apply to the modeling when dealing with zircon from volcanic products. The properties of a natural zircon age distribution are captured by common analytical techniques such as ID-TIMS, LA-ICP-MS, and SIMS, which afford age resolution to different degrees depending on their required sampling volumes (Table 1). To demonstrate the potential differences between these types of data in an idealized case, we simulate a pool of 200 zircon crystals (assumed to be spherical for simplicity) spanning an age range of 3 m.y., with radii between 50 µm and 100 µm, and with the thickness of the different zones selected randomly but arranged to decrease from the crystal center toward the rim (a selection is shown in Fig. 13A). The U content of the different zones and the variation of the zircon zone thickness from core to rim affect the distributions of ages obtained from whole zircon versus high spatial resolution techniques (Curry et al., 2021). In this simulation, each zircon comprises 2–10 zones of decreasing age from core to rim and only zircons with 10 zones have a core with the oldest age (i.e., 3 Ma). During the prescribed 3 m.y. of magmatic history that zircon could potentially record, the mass fraction of zircon crystallized increases as the amount of magma capable of crystallizing zircon grows, but then drops toward the end of the interval to mimic the waning stages of the magma system, consistent with the expected zircon age distributions introduced in section 3.3 (Fig. 13B). We then use this simulated zircon production curve to calculate the distribution of zircon ages that would be obtained by performing ID-TIMS analyses by summing the product of the mass fraction of each sector and the respective age of the sector (Fig. 13B). Not surprisingly, this calculation shows that the range of zircon ages that would be obtained using ID-TIMS is significantly shorter than the “real” duration of the magmatic event during which zircon was crystallizing (Fig. 13; Glazner and Sadler, 2016). However, it is interesting that the median and mode for the calculated ID-TIMS analysis are relatively close to the real values, whereas the standard deviation is generally underestimated by the subsamples (Fig. 13D). Thus, for ID-TIMS age determinations, mode and median provide relatively reliable metrics for comparison of modeling results, whereas the true standard deviation of ages is underpredicted. This test implies that using ID-TIMS ages in combination with the model proposed by Caricchi et al. (2014, 2016) is appropriate and produces relatively good estimates for the rate of magmatic input into the system but underestimates the volume of the system.

Using the same simulated population of zircon crystals, we calculated the distribution of ages that would be measured by placing 1100 analysis spots of a high spatial resolution technique (LA-ICP-MS or SIMS) with a beam diameter of 30 µm (Fig. 13A). We assume that crystals are analyzed in cross-sectioned mid planes, disregarding mixing of crystal domains at different depths for simplicity. High spatial resolution analyses better capture the total duration of zircon crystallization, but because of lateral mixing of growth domains, they provide biased estimates for the mode and median of the actual distribution of zircon ages (Fig. 13D). Therefore, the comparison between the calculated distribution of zircon ages and measurements performed using high spatial resolution techniques should focus on the total spread of zircon ages (e.g., as described by the standard deviation) rather than on the mean or median. The exceptions are analyses of unsectioned crystals acquired in depth-profiling mode (e.g., Schmitt and Vazquez, 2017), in which the ion beam ideally samples individual growth domains with minimal mixing as the depth resolution of the ion beam (<1 µm) is small relative to the dimensions of the growth zones. Surface analysis of individual and unsectioned zircon crystals also allows age determination of the last growth episodes largely without mixing different growth zones, and thus not only the total spread can be modeled, but also the shape of the probability density function (i.e., the relative proportions of newly crystallized zircon relative to antecrystic grains).

Although direct comparability between the currently applied models as summarized in section 3.2 is complicated by the fundamental differences in how they are conceptualized and coded, we can demonstrate general similarities between the predicted thermal histories for magmatic systems open for recharge over typical lifetimes of volcanic systems (e.g., 1–1.5 m.y.; Weber et al., 2020). To mitigate model-dependent bias, we calculated zircon crystallization ages for both models using identical boundary conditions. For this, conditions typical of a longlived arc volcano were selected, which include a normal initial geothermal gradient of 30 °C and a magma recharge temperature of 1000 °C (Weber et al., 2020). Furthermore, both types of models were set up to produce the same total magmatic volumes with an approximately identical depth range of the accumulated intrusion after the end of each run (Fig. 12). A corollary of the cylindrical sill geometry with downward displacement is that the magma always comes into contact with comparatively cold crust, and the magma will thus be more affected by rapid initial cooling than the outward growing ellipsoid, as earlier intruded and hot magma/rock shields the interior of the ellipsoid from heat loss to cold country rocks. To compensate for a protracted initial incubation time in the Geneva models during which zircon productivity is negligible, the run duration of the UCLA-HD models was shortened to only match the zircon-producing duration of the Geneva models (as defined by the presence of ~1% of final magmatic volume). The volume of the initial intrusion that becomes inflated by continuous recharge (in pulses every 5000 years) was prescribed as ~10% of the final intrusive volume in the UCLA-HD model, and the aspect ratio of the ellipsoid was held constant (horizontal:vertical axes = 4:1). As a result, the average magmatic temperatures are always lower for the Geneva model than for the UCLA-HD model, with temperatures gradually increasing and decreasing with the run duration, respectively (Fig. 12). This is expectedly matched by the longer duration of zircon crystallization (measured as the difference between the 16^th and 84^th percentile in the cumulative zircon age distribution) in the overall cooler Geneva models relative to the hotter UCLA-HD models at the modeled magmatic recharge rates of 7.5 × 10⁻⁶ to 1.1 × 10⁻⁵ km³/yr/km² (Fig. 12). Despite these differences, the overall shape of the cumulative probability distribution of expected zircon ages is similar, which implies that under the selected conditions, sampling of younger zircon is more likely than sampling of older crystals due to progressive growth of the zircon-saturated thermal layer.

Despite conceptual and computational differences in the setup of thermochemical models for zircon crystallization in evolved magmatic systems fed and maintained by a repetitive influx of more primitive magma, they provide broadly similar results. This can yield unique quantitative constraints on longterm magmatic fluxes and magmatic crustal accretion in a specific tectonic environment. In contrast to geophysical observations that monitor presentday conditions, modeling the thermal evolution of magmatic chambers as tracked by their zircon inventory can reveal longterm variations in magmatic input in long-lived systems, and it can also be applied to fossil magmatic systems in the form of plutons.

One important insight gained from the modeling of zircon age spectra is that substantial intrusive volumes accumulate under arc volcanoes. Hence, intrusive to extrusive ratios determined by this method for volcanic arcs range from 25:1 to 50:1 for Nevado de Toluca, a stratocone in the Trans-Mexican Volcanic Belt (Weber et al., 2020), to 75:1 for lava domes in the Altiplano-Puna region, which erupted after an ignimbrite flareup episode (Tierney et al., 2016). These exemplary estimates significantly exceed the commonly inferred intrusive to extrusive ratio of 10:1 (Crisp, 1984). The intrusive to extrusive ratio also changes as a function of time and tends to decrease in thermally mature volcanic systems (Giordano and Caricchi, 2022, and references therein).

Zircon growth models also quantitatively predict volumes of eruptible magma in the subsurface that can be directly compared with those estimated from other petrological constraints (e.g., Weber et al., 2020; Friedrichs et al., 2021) or geophysical observables (e.g., Lukács et al., 2021). To illustrate how zircon age spectra and their quantitative modeling can reveal fundamental differences between the underlying magma systems, we compare data from four volcanoes recently studied (Fig. 14). Here, this can only be accomplished in an abridged form, and we refer to the cited publications for more detail (Weber et al., 2020; Friedrichs et al., 2021; Lukács et al., 2021).

One key outcome of these studies is that potentially eruptible volumes of magma can be significant in systems that are currently dormant, or even considered extinct by conventional standards, but also that volumes of potentially eruptible magma vary significantly in different systems. This is important in assessing mid- to longterm hazards from volcanoes (e.g., Lukács et al., 2021; Weber et al., 2020). The estimates for the total amount of eruptible magma (at <50% crystal content) derived from modeling zircon age spectra are comparatively small (albeit non-negligible) for the Ciomadul dome field of Romania at ~5–10 km³ (Lukács et al., 2021), whereas for Nevado de Toluca, the total amount of eruptible magma is about one order of magnitude larger (Figs. 14A and 14B; Weber et al., 2020). This agrees with the overall differences in cumulative eruptive volumes (8–10 km³ for Ciomadul versus ~60 km³ for Nevado de Toluca), and eruptive frequency (e.g., last major eruption at ca. 30 ka for Ciomadul versus ca. 14 ka for Nevado de Toluca).

Zircon age spectra can also indicate variations in magmatic recharge fluxes over time (Friedrichs et al., 2021; Angeles-De La Torre et al., 2023). This is the case for two adjacent stratovolcanoes in Central Anatolia (Mt. Hasan and Mt. Erciyes) that show different styles of continuous versus pulsed recharge throughout the late Pleistocene (Figs. 14C and 14D; Friedrichs et al., 2021). One important aspect of the study by Friedrichs et al. (2021) is that zircon crystals were separated and analyzed following identical protocols for samples from both volcanoes so that the observed differences in the age spectra are robust. Moreover, different patterns in longterm recharge rates, steady and frequent for Mt. Hasan versus fluctuating and infrequent for Mt. Erciyes, correlate with the eruptive frequency and type for each volcano, which are characterized by frequently recurring lava flows of intermediate (andesitic) composition for Mt. Hasan and initially highly explosive eruption bursts of evolved (rhyolitic) magmas from multiple vents for Mt. Erciyes (Friedrichs et al., 2021). Zircon age spectra modeling suggests that the active magmatic system underneath Mt. Hasan is larger and more frequently rejuvenated by mafic magmatic input, whereas that of Mt. Erciyes has shrunk to contain mostly residual evolved melts in what may represent the waning stages of its magmatic system, despite the overall larger volcanic edifice volume compared to that of Mt. Hasan (Friedrichs et al., 2021).

As these examples illustrate, interpretation of zircon age spectra should ideally integrate as much geological information as possible, including the eruptive pre-history of a volcanic system. The onset of magmatic input into the system is especially important for refining estimates of the rates of magmatic input at depth. This is because the dominant factor controlling the thermal evolution of the magmatic system is the total amount of sensible heat added into the crust, which is the product of the rate of magmatic input and time. The onset of magmatic input can either be obtained from the oldest (non-xenocrystic) zircon or by dating the first eruption of a volcanic system, e.g., with ⁴⁰Ar/³⁹Ar or (U-Th)/He chronometers. On one hand, both approaches will likely underestimate the duration of magmatic input into the plumbing system, as zircon only starts crystallizing once sufficiently evolved melts are generated, and products of older eruptions may be covered or destroyed. On the other hand, it is sometimes difficult to distinguish between antecrysts and xenocrysts (Miller et al., 2007), which may lead to an overestimation of the onset of magmatism.

The thermal models presented in section 3.2 illustrate how for similar zircon age distributions, the average temperature can change significantly when different geometries of magmatic input into the plumbing system are considered (Figs. 12D–12F). Such differences should be reflected in different abundances and temporal trends in the chemistry of volcanic products, as well as the trace-element characteristics of zircon. Hence, knowing the temporal and chemical evolution of volcanic products and zircon from eruptions spanning the widest possible time interval provides additional constraints for identifying the most appropriate thermal model. This, in turn, allows the retrieval of more accurate and precise estimates of the average rates of magmatic input into the volcanic plumbing system. The same approach could be applied to exposed crustal sections by focusing on plutons spanning the widest possible range of emplacement ages. Studies that complemented zircon age spectra with thermochemical indicators, such as trace element data, have demonstrated this potential (e.g., Ti- in- zircon; Samperton et al., 2017; Kent and Cooper, 2018; Weber et al., 2020; Friedrichs et al., 2021).

Models such as those shown in Figures 11 and 12 indicate that protracted zircon crystallization in a longlived system open for recharge results in a distribution of zircon ages in volcanic rocks that differs significantly from that expected from zircon saturation in a magmatic volume that (unrealistically) cools without any internal thermal gradients. Zircon age spectra, successfully modeled in terms of opensystem magmatic recharge, thus may have little in common with the simple exponential decrease in zircon production predicted from zircon saturation models for a single batch of melt undergoing monotonic cooling (e.g., Harrison et al., 2007). The models described in section 3.2 suggest zircon distributions for longlived silicic magma systems with younger zircon predominating over older zircon, the opposite of the distribution expected for cooling of a single-magma batch (Figs. 9 and 10). This implies that it is not straightforward to prescribe an a priori age distribution to zircon age spectra to extract eruption ages from zircon data. Thus, extrapolating from zircon ages to a precise eruption age (e.g., Keller et al., 2018) requires firm knowledge of the pre-eruptive zircon crystallization history of a particular magmatic system. Because zircon ages track the pre-eruptive accumulation of magma, they are an important complement to eruption ages for unraveling potential leads and lags regarding climatic variability (e.g., sea-level changes and glaciation–deglaciation) and associated lithospheric stresses that may influence the productivity of magma, recharge fluxes, and eruptive frequency (e.g., Satow et al., 2021).

To further advance zircon as a magmatic thermochronometer, improvements in analytical methodology and more realistic models for zircon crystallization, residence, and mobilization in magmas are required. Because zircon age spectra based on a few analyses may not reliably represent the underlying distribution, large-n data sets (n > 100) will potentially provide additional information, in line with the approach advocated for state- of- the- art detrital zircon studies (Pullen et al., 2014). Moreover, there are analytical limitations in resolving the true duration of zircon crystallization (Table 1). To further utilize the potential to reconstruct thermal histories from individual zircon crystals, improved spatial resolution and precision are essential. Zircon trace-element geochemistry also furnishes complementary information such as temperature (Ti- in- zircon) and melt evolution (e.g., indices for differentiation such as Zr/Hf and Eu/Eu*, as well as for crustal assimilation such as δ¹⁸O and εHf) that can be used alongside and in support of the age data. Lastly, the understanding of the significance of zircon ages is limited by the fact that, in most studies, textural relationships of zircon to other mineral phases are unknown because zircon is extracted from its matrix by physical and chemical methods. Therefore, it remains ambiguous whether an individual zircon crystal stopped crystallizing because it became entombed in major minerals and therefore was shielded from the melt, or because of subsolidus cooling and complete solidification of the magma.

Although zircon saturation is well constrained, zircon can crystallize before wholesale saturation is reached when melt composition and Zr abundance become locally conducive for zircon crystallization, e.g., when zircon crystallizes next to fastcrystallizing oxides (Bacon, 1989). A related critical aspect is that evolved (per-)alkaline system compositions currently defy modeling zircon crystallization because experimental constraints for zircon saturation in such melts are scarce (cf. Gervasoni et al., 2016). Empirically, zircon crystallizes late in such systems due to the high solubility of Zr and Hf in alkaline melts (Linnen and Keppler, 2002), but then crystallization from supersaturated, volatile-rich, and low-viscosity melts can produce large amounts of zircon, often much larger in grain size than those that crystallized from metaluminous or peraluminous melts. The processes of zircon saturation and growth in such melts, along with the formation of hydrothermal zircon (e.g., Troch et al., 2018), should receive more attention in the future.

The thermal models presented in section 3.2 are severely restricted in realistically describing mixtures between crystals and melt where relative proportions of solids and liquids are expected to fluctuate spatially and temporally over the lifetime of a magma system. Remobilization of crystal-rich marginal domains or convective overturn in a magma prior to volcanic eruptions are disregarded in these models, although both processes are likely common in nature. Rounded and partially melted plutonic enclaves are frequent ingredients of volcanic deposits, which suggests that the reheating of largely solidified magmatic crusts can liberate zircon from otherwise non-eruptible domains of a magmatic system (e.g., Schmitt et al., 2010). Furthermore, massive recharge events are expected to cause remobilization, stirring, and homogenization over large distances in a longlived crystal mush (e.g., Bachmann and Bergantz, 2006; Huber et al., 2009). This seems to be a prerequisite for indistinguishable zircon age spectra at the handspecimen scale for the same eruption despite these crystals originating in different parts of a magmatic reservoir and thus recording different thermal histories (e.g., Pamukçu et al., 2022). The frequency and intensity of these remobilization and overturning events in longlived magmatic systems are empirically poorly constrained.

Creating a petrologically grounded thermomechanical model applicable to magmatic systems is one of the biggest challenges of geodynamic modeling. Indeed, such a model has to integrate: (1) complex elastoviscoplastic rheologies for the deforming host rock; (2) several modes of magmatic transport/segregation, including mass convection, pervasive flow within an interconnected solid network (low melt/rock ratio), settling of crystals within a magmatic body (high melt/rock ratio), and diking through the hostrock (high melt/rock ratio); (3) endmember cases where either fully solidified rock or fully liquid magma is present; (4) extremely large viscosity difference between magma and solid host; (5) thermodynamically consistent phase stabilities and transitions, involving reaction kinetics, volume change, and thermal stress; and (6) the coupling of tectono-magmatic processes occurring over a wide range of temporal and spatial scales in a self-consistent and numerically stable way.

State-of-the-art approaches are currently divided into three main categories that tackle these problems from different angles. The first category aims to better couple thermomechanical deformation with twophase flow while disregarding or simplifying phase reactions (e.g., Keller et al., 2013, 2017; Keller and Katz, 2016; Turner et al., 2017; Rummel et al., 2020; Keller and Suckale, 2019; Katz et al., 2022). Although significant breakthroughs have been achieved during the last decade, to date there are no elastoplastic thermomechanical models that are able to account for volume change and phase reactions simultaneously. Instead, the second category has focused on coupling twophase flow models with phase reaction while disregarding thermomechanical deformation (e.g., Jackson et al., 2003, 2005; Annen et al., 2006; Bouilhol et al., 2011; Solano et al., 2012; Jackson et al., 2018; Riel et al., 2019). Here, several 1-D models of reactive twophase flow have been successfully conducted, but only a limited number of 2-D applications have investigated the formation of a reactive meltfilled network (e.g., Bouilhol et al., 2011), and all ignore the role of elastoplastic thermomechanical deformation of the hostrock, including diking. Moreover, the coupling with phase reactions is mainly achieved using simplified parameterizations (e.g., Jackson et al., 2003, 2005, 2018; Solano et al., 2012), which provides useful insights into the overall model behavior but makes it difficult to directly compare the simulations with geochemical data. The third category aims to better understand the relationship between magmatic chamber recharge and eruption. Using a 0-D modeling approach, these models consider up to three phases (solid, liquid, and gas), visco-elasto-plastic rheology of the surrounding crust, parameterized thermodynamically consistent phase changes within the magma chamber, as well as pressure variations related to recharge and eruptions events (e.g., Degruyter and Huber, 2014; Townsend et al., 2019; Townsend and Huber, 2020; Townsend, 2022; Liao et al., 2021). Although the latter category allows scaling relationships to be drawn based on eruption frequency and size (e.g., Degruyter and Huber, 2014; Townsend et al., 2019; Townsend and Huber, 2020; Townsend, 2022), it falls short of describing the internal dynamics of a magmatic chamber and has yet to be applied to higher spatial dimensions.

Melt segregation within magmatic chambers and cycling between magmatic recharge and extraction in and out of a magmatic reservoir imply that compositions of crystals and melts are constantly evolving. Hence, the use of phase equilibrium parameterization or look- up- tables of fixed bulk-rock composition cannot predict the chemical evolution of the magmatic chamber. Alternative methods, such as a pre-computed petrological database covering the expected compositional space (Rummel et al., 2020), or direct coupling with a Gibbs free energy minimizer (Riel et al., 2019, 2022) is a step forward. The coupling of a physical model of melt extraction and crystal accumulation in a magma chamber to a petrological database (Rummel et al., 2020) or a Gibbs energy minimizer (e.g., Riel et al., 2019), while challenging on its own, can allow the mechanical and chemical evolution of a magmatic reservoir to be modeled. With adequate zircon saturation models or a thermodynamic database that includes Zr (Kelsey and Powell, 2011), the growth and resorption of zircon can be readily predicted. However, the sluggish dissolution kinetics of zircon in silicate melts (e.g., Watson, 1996; Bindeman and Melnik, 2016) also play an essential role in controlling the available Zr pool for zircon crystallization. Therefore, the coupling of thermomechanical models for magmatic chamber dynamics and petrological reactions also needs to account for non-equilibrium thermodynamics. A possible solution could be to use stable phase equilibrium prediction as a directional step (e.g., Lasaga, 1986), while the length of that step, which corresponds to the reaction kinetics, could be constrained by experimental melting studies (e.g., Hou et al., 2020).

Zircon, although a tiny and seemingly negligible component in rocks, is nonetheless an essential probe into magmatic systems concealed at depth. It is not only a time capsule, but also a faithful recorder of the thermal and chemical conditions in magmas. For these reasons, zircon has become indispensable for studying differentiation, longevity, and pre-eruptive priming of magmatic systems, where it adds a retrospective capability that geophysical monitoring restricted to brief human timescales lacks. Although zircon age dispersion is seemingly a burden from the traditional perspective of an event-oriented geochronology, the increasing power to unravel the complexity of accessory mineral crystallization in time and space also bears tremendous potential. In this review, we emphasized the need for improved statistical tools to adequately portray and interpret geochronological data and outlined potential strategies to advance conceptual models for dynamic magma systems. Progress in this regard requires the integration of empirical studies of volcanic and plutonic rocks that combine accessory mineral geochronology with other geological and petrological information such as (chrono-)stratigraphy, as well as correlative trace-element and isotope data, and lastly refinement of numerical simulations of the thermomechanical evolution of multiphase magmatic systems in dynamic settings. Further study of the thermal and compositional evolution of magmatic systems has direct implications for assessing volcanic hazards, constraining the heat sources of magmatic geothermal systems, and understanding the formation of ore deposits associated with magmatic intrusions.

This contribution is the outcome of a workshop, “ZASSy” (Zircon Age Spectra Systematics), which was generously supported by the Terrestrial Magmatic Systems (TeMaS) Research Platform of Heidelberg and Mainz Universities, as well as by a MAGMA ERC grant to Boris Kaus. Luca Caricchi received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 677493 – FEVER) and the Swiss National Science Foundation (grant no. 200021_184632). Olivier Bachmann acknowledges funding from the Swiss National Science Foundation (grant no. 200021_178928). We thank journal reviewers Adam Kent and Blair Schoene, as well as associate editor Jonathan Miller for insightful comments. Marco Pfennig of Gau-Bickelheim, Germany, is thanked for a stimulating mid-workshop excursion.

Modeling the thermal evolution of magmatic systems has a long tradition, and a range of commercial and non-commercial codes have been employed that intrude dikes and/or sills, and simulate the resulting thermal evolution of the intrusion and surrounding rocks. As the modeling assumptions, boundary conditions, and material parameters employed often differ between publications, it is not always straightforward to compare results from different studies. In addition, none of the software currently in use is available as an open source software, which makes it difficult to reproduce results or understand the finer implementation details.

To amend this situation, we have developed a new, open source, software package—MagmaThermoKinematics.jl (https://github.com/boriskaus/MagmaThermoKinematics.jl) (MTK)— using Julia (Bezanson et al., 2017; https://julialang.org), which can be employed to simulate the thermal evolution of magmatic systems. Julia is a highlevel programming language that has some similarities to MATLAB or Python (and allows the creation of short, easy- to- understand code), but is compiled and therefore significantly faster than the above-mentioned interpreted languages. Moreover, it has outstanding support for multi-threaded parallel central processing units and GPUs, an excellent testing framework, and a fairly rich ecosystem of readily available packages that are written in a highly composable manner, which makes it easy to combine different packages.

The governing advection-diffusion equation relevant for the thermal evolution of magmatic systems with intruding magma is:

Here, T denotes temperature [K], v hostrock velocity (induced by injecting a new sill), t time, ϕ = (1−ϕ) the volume fraction of solid rock, ϕ the melt fraction (both dimensionless parameters between [0–1]), ρ density [kg m^-3], c heat capacity [J kg^-1 K^-1], k conductivity [W m^-1 K^-1], Q_L latent heat [J kg^-1], and H radioactive heat sources [W m^-3]. Many of the material properties can also be a function of temperature, which makes this a nonlinear equation.

In general, we solve Equation (A1) into two steps.

An advection step is applied only when a new sill is intruded, which happens every few 1000 years in most cases discussed here.

We employ a Semi-Lagrangian advection scheme for this (e.g., Staniforth and Côté, 1991; Becker and Kaus, 2016), which takes into account that the host rock moves to create space for the intruding magma. The hostrock velocity is either set to be a constant below a horizontal sill (in the Geneva setup), can be elliptical in mass-conserving manner (in the UCLA-HD setup), or can be set from an elastic dike intrusion algorithm.

A diffusion step is applied more frequently:

The last term in the equation involves the time derivative of the solid fraction. It turns out that this is numerically rather unstable and requires a very small time-step. Therefore, we apply the chain rule to this term and write this as:

Here, we made use of the fact that ∂ϕ_s/∂T = −∂ϕ/∂T. This formulation shows that the latent heat effect can be implemented by locally changing the heat capacity (e.g., Simpson, 2017; Gerya, 2019), which is numerically much more stable than if using Equation (A2b) directly. It requires knowledge of the derivative of the melt fraction versus temperature, ∂ϕ ∕ ∂T, which can be computed analytically for most melt parameterizations and is obviously zero below the solidus and above the liquidus. Yet, an important aspect for stable numerical implementations is that this derivative must be continuous throughout the full temperature range. Some melting parameterizations result in a significant jump in ∂ϕ ∕ ∂T around the solidus (or liquidus), which puts severe timestep restrictions on the numerical simulation. One solution is to apply a small amount of smoothening to existing parameterizations around the solidus and liquidus to enforce continuity. This is an option implemented in GeoParams. jl and allows larger time-steps to be taken while giving essentially the same results.

In MTK, the time-derivative in Equation (A3) is discretized in an explicit manner:

The spatial derivative, , is implemented using a staggered grid finite difference scheme, which defines thermal conductivity inbetween temperature grid points (e.g., Gerya, 2019). In 1-D, this would be:

where Δx is the (constant) spacing between temperature gridnodes that are defined at points indicated with a superscript i. The thermal conductivity is physically located at a point that is exactly in between two temperature grid points. In case of temperature-dependent conductivity, we first interpolate temperature to that same grid point, i.e., .

The advantage of employing the explicit approach is that it results in an algorithm that can be quite efficiently implemented on the GPU, for example by using the ParallelStencil.jl Julia package (e.g., Räss et al., 2022). A disadvantage is that there are time-step restrictions due to the Courant-Friedrichs-Lewy criteria, which state that:

where is the thermal diffusivity and Δx, Δy, Δz the grid spacing in three directions (for 3-D). In case material parameters are spatially variable or depend on temperature, one can use an upper bound to estimate κ.

Equation A4 works for the case that parameters are independent of temperature, which is generally not the case (ϕ and ∂ϕ ∕ ∂T obviously depends on T). Therefore, one typically has to employ nonlinear iterations. In the general case that ρ, c_p, k, and ∂ϕ ∕ ∂T depend on T, this gives:

We employ dampened Picard iterations, which are computed as:

where ω = [0 − 1] is a dampening parameter and the superscript ^it denotes the iteration step. Equations (A7) and (A8) are solved until the relative difference between subsequent iteration steps,

becomes smaller than a predefined tolerance (e.g., 10^-5).

In MTK, the governing Equations (A7) and (A8) are implemented in 2-D, 2-D axisymmetric, and 3-D geometries. Finite difference stencils are implemented using ParallelStencil. jl (https://github.com/omlins/ParallelStencil.jl), which is a Julia package for efficient stencil computations that allows the simulations to be run on either parallel multicore processors or on GPU cards.

Parameterizations for material properties (such as thermal conductivity, melt fraction, or heat capacity) are implemented in a separate package (GeoParams.jl, https://github.com/JuliaGeodynamics/GeoParams.jl), which allows reuse of those in independent software packages (while having to implement and test them only once). We provide a range of melt-intrusion algorithms, which include elliptical or cylindrical geometries, or follow the analytical solution for intrusion of a dill/dike in an elastic halfspace. Tracers can be injected to track the thermal evolution of the intruded magma and host rock, and there is a simple mechanism to change the type of material properties employed (e.g., using a temperature-dependent conductivity instead of a constant value). Both flux and isothermal bottom boundary conductions can be implemented.

The results of each calculation comprise, besides the thermal field, time-temperature paths on every tracer, which can be used in a postprocessing step to compute zircon age data (using a Julia implementation of R-scripts provided by Weber et al., 2020, which is provided in GeoParams.jl). The file example 2D_ZASSy.jl reproduces the results shown in Figures 11 and 12, which include both endmember models.

Petrological-thermal modeling studies of magmatic systems often use slightly different assumptions, such as constant versus temperature-dependent conductivity, 2-D versus axisymmetric geometries, or employ different melt-curve parameterizations. It is often unclear what the impact of those assumptions will be on the final temperature distributions of the magmatic system. For that reason, we performed a sensitivity analysis using the two endmember scenarios described in the manuscript (Figs. 11 and 12). In each case, we kept all parameters the same as in the reference simulation but changed one parameter. We report the final temperature profiles in the center and at the side of the domain. Results are shown in Fig. S1 (see text ^{footnote 1}) and demonstrate that the numerical parameters (such as resolution or timestep) only have a minor effect, but that latent heating and temperature-dependent conductivity can have an impact that is equal to or greater than reducing the magma flux for the sill under-accretion setup, whereas they have little impact on the central intrusion setup. This is likely because the central intrusion setup has nearly reached a steady state, where magma intrusion balances thermal diffusion and the partially molten zone remains more or less constant in size. The sill under-accretion case, on the other hand, has a partially molten zone that still grows in space after 1.5 m.y. (Fig. S1).