What’s in an age? Calculation and interpretation of ages and durations from U-Pb zircon geochronology of igneous rocks

Meaningful interpretation of geochronologic data critically depends on the contrast between the analytical precision of individual measurements and the duration of the geologic event of interest. Since the development of modern geochronology 70 years ago, the ratio between these values was typically assumed to be large, indicating that the geologic event of interest (pluton emplacement, volcanic eruption, etc.) occurred on time scales much shorter than the available analytical precision. This assumption allowed geochronologists to treat the process as instantaneous relative to the date’s reported precision and led to the widespread use of weighted mean ages and adoption of the mean square of weighted deviates (MSWD; Wendt and Carl, 1991) to assess the validity of this assumption. While this assumption may remain valid for studies employing relatively low-precision geochronology techniques or for geologic events of extremely short duration, it is increasingly common for high-precision geochronologic data sets to be “overdispersed,” that is, to show interanalysis variance greater than that expected from analytical errors alone (e.g., Rioux et al., 2012; Wotzlaw et al., 2013; Samperton et al., 2015; Andersen et al., 2017; Ellis et al., 2017). These data indicate one or more of the following: (1) analytical uncertainties have been underestimated; (2) the chronometer has experienced open-system behavior; (3) the analyzed sample includes material that crystallized prior to the event of interest (i.e., xenocrystic grains or inherited cores); and/or (4) the measurement uncertainties represent a duration of time shorter than the duration of the geologic event of interest.

Ongoing community-driven advances have helped to reduce or eliminate the first two causes of overdispersion in U-Pb zircon geochronology. Underestimated analytical uncertainties can be ruled out as a source of excess scatter by properly assessing the long-term reproducibility of isotopic standards and natural reference materials within and between individual laboratories (e.g., Košler et al., 2013; Schaltegger et al., 2021). Likewise, pre-analytical treatments such as chemical abrasion can effectively isolate domains less impacted by open-system behavior (Mattinson, 2005; Widmann et al., 2019). Application of these approaches yields increasingly precise, and presumably accurate, U-Pb zircon dates that are frequently overdispersed and interpreted to be the result of geologic processes that occurred over time scales much longer than the reported analytical uncertainties. These resolvable age distributions can provide additional constraints on the durations and rates of geologic processes, and the community has begun to apply increasingly complex methods to leverage this information (Caricchi et al., 2014, 2016; Glazner and Sadler, 2016; Keller et al., 2018; Ratschbacher et al., 2018; Vermeesch, 2018).

As analytical advances have minimized some causes of overdispersion, interest in understanding the remaining (geologic) overdispersion has increased, and several factors have come into focus in the interpretation of U-Pb zircon geochronology. These include factors that affect how an individual zircon date should be interpreted, the conditions necessary for zircon dates to meaningfully record duration within a single igneous system, and the effect of sample size on estimated durations (e.g., Frazer et al., 2014; Samperton et al., 2015; Glazner and Sadler, 2016; Rivera et al., 2016; Kent and Cooper, 2018; Curry et al., 2021; Gaynor et al., 2022). In this article, we examine each of these variables as they relate to interpretation of the record of zircon crystallization in igneous rocks. Through this discussion, we identify a set of best practices for estimating the duration of zircon crystallization. Although the examples we use are specific to zircon geochronology in magmatic rocks, many of the insights from our analyses are also relevant to other applications of U-Pb zircon geochronology, to U-Pb geochronology of other minerals, and to other geochronologic systems.

The age distribution of zircon produced in a crystallizing magma is a function of bulk composition and cooling rate, which, respectively, determine the temperature at which zircon saturates and the time spent between zircon saturation and the solidus. The conditions needed for zircon saturation are well known from experimental studies (e.g., Boehnke et al., 2013), while the rates at which intrusions cool are a function of a number of factors, including the magma temperature and emplacement rate, the background geotherm, the total size of the intrusion, and, in shallow systems, the efficiency of associated hydrothermal systems (e.g., Gelman et al., 2013; Annen et al., 2015; Karakas et al., 2017; Kelly et al., 2021). In most exhumed plutonic systems, the properties that governed magma residence are difficult to know a priori, and attempts to constrain them from zircon age distributions have become more common. However, it is necessary to rigorously interpret the significance of measured dates from individual grains as well as the significance of age distributions within populations of grains to properly utilize these data.

Significant recent progress has been made toward modeling the distribution of zircon ages as a function of the thermal evolution of a magmatic system (Caricchi et al., 2014, 2016) or using models of zircon crystallization and growth limited by Zr diffusion in the crystallizing melt (Sorokin et al., 2022). However, these distinct modeling approaches currently do not produce consistent predictions for zircon nucleation and growth (cf. Keller et al., 2018; Sorokin et al., 2022), and some models make predictions only for the total mass of zircon crystallized, which does not directly relate to the zircon age distribution without additional constraints on zircon growth and nucleation rates. Further, existing models do not account for inclusion of zircon grains by other crystallizing phases (e.g., Bacon, 1989; Clarke et al., 2022), which can halt individual grain growth significantly prior to the system cooling through the solidus. Therefore, to better understand the relationships between these factors and their bearing on the interpretation of zircon ages, we simplify this problem by constructing idealized synthetic zircon populations generated using a set of prescribed parameters.

In these synthetic zircon populations, we assume that individual zircon grains grow concentrically as right square prisms at a volumetric growth rate [G(t)] defined over the nondimensional time interval [0,1] (definitions of all variables listed in Table 1). For a zircon that nucleates at time τ₀ (where 0 ≤ τ₀ ≤ 1) and stops growing at time τ₁ (where 0 ≤ τ₁ ≤ 1 and τ₀ ≤ τ₁), the zircon volume is defined as:

and, given an aspect ratio h (defined as l/x, i.e., the length of the major axis relative to the minor axis, which corresponds to the length of a side of the square base), the length of the zircon major axis as a function of time is defined as

Finally, if zircon U concentration is assumed to be constant, the volume-integrated mean age of a zircon can be calculated as

[Note that this definition of is equivalent to the expected value of G(t), because the growth rate acts as the probability density function for zircon age.] This formulation is also readily extended to account for variable U incorporation during zircon crystallization. The effect of varying U concentrations is considered briefly below, but for simplicity, we otherwise assume that zircon grains crystallize with constant U concentrations. Additionally, while varying the zircon aspect ratio will influence microbeam analyses of zircon, it does not affect volume-integrated mean ages. Therefore, for simplicity, we model all zircon with an aspect ratio of 5, as appropriate for the elongate, needle-shaped morphology often observed for zircon that crystallized in relatively shallow igneous systems (Corfu et al., 2003), but we also evaluate the effect of varying aspect ratios in Figure S1¹.

Using these model zircon grains, we first examine how differences in G(t) and composition of a single zircon grain affect age measurements of either dissolved whole zircon grains or grain fragments (i.e., modern chemical abrasion–isotope dilution–thermal ionization mass spectrometry [CA-ID-TIMS] measurements), or age measurements from microbeam spot analyses of zircon polished to half thickness parallel to the major axis, analogous to typical measurements made using laser ablation–inductively coupled plasma–mass spectrometry (LA-ICP-MS) or secondary ion mass spectrometry (SIMS). We consider multiple forms for G(t), including uniform with time, and linearly increasing or decreasing growth rates (Table 2). Spot analyses are assumed to be cylindrical and to intersect the polished zircon face at a right angle and can have variable radii relative to the zircon size. Finally, we assume that the depth of the cylinder formed by spot analyses is either comparable to the radius for LA-ICP-MS analyses (appropriate for typical spot sizes of 30–35 μm with ~30 μm pits) or one hundredth of the spot radius for synthetic SIMS analyses. All volume-weighted ages from both spot analyses and zircon fragments are calculated numerically by discretizing zircon growth over 1000 time steps.

Following this investigation of single zircon, we extend our analysis by generating synthetic zircon populations as a function of the above parameters and a similarly prescribed nucleation rate, J(t), as well as an armoring probability, P_armor, and we use these data to explore the measurement of crystallization duration from zircon populations. Several different nucleation rates and armoring probabilities are explored in separate simulations (Table 2). In these simulations, J(t) defines the quantity of zircon crystallizing at a given τ₀, and each grain stops growing either upon armoring, or when the system reaches the solidus, or quenches, at t = 1. When generating zircon populations, we assume that a single G(t) is applicable to all zircon grains present in the system. This is undoubtedly a simplification but will be appropriate for systems where zircon growth is controlled primarily by the bulk zircon saturation in the system.

For the end member where zircon growth is instantaneous (i.e., τ₀ = τ₁ or P_armor = 1), the age of any grain corresponds to its time of nucleation, and normalizing J(t) by the total number of zircon grains crystallized gives the probability density function (PDF) of zircon ages (t):

For noninstantaneous zircon growth, this function defines the distribution of zircon nucleation events. If P_armor = 0, then the PDF of volumetric weighted mean zircon ages, , can be calculated as:

where E⁻¹ is the inverse function that uniquely maps to τ₀. Finally, the normalized instantaneous volume (or equivalently mass) of zircon crystallizing in the system (V_sys) as a function of time can be calculated as:

For some G(t) cases that are always positive (i.e., for systems that are always zircon saturated), E⁻¹ can be found analytically, and both and Vsys can be calculated explicitly. However, in many instances, and for any situation with 0 < P_armor < 1, must be found from numerical simulations of zircon populations.

The simulated zircon populations created here are functions of three parameters (growth rate, nucleation rate, and armoring probability) that we assume vary independently (Table 2). We examine a range of shapes for each of these parameters to explore the impact of each variable on the resulting zircon age distributions. While all three of these variables are relatively poorly constrained in natural systems, when possible, we have attempted to include scenarios suggested by previous studies. Explored nucleation rates include uniform nucleation or linearly increasing nucleation, as well as a nucleation rate defined by a beta distribution that approximates previously suggested zircon mass crystallization distributions (e.g., Samperton et al., 2017; Keller et al., 2018) and an exponentially decreasing rate that approximates kinetic nucleation rate laws (e.g., Marsh, 1988; Cashman, 1993). Growth rates considered included uniform, linearly increasing, and linearly decreasing rates, and an exponentially decreasing growth rate suggested by diffusion-limited growth (e.g., Watson, 1996; Sorokin et al., 2022). Finally, armoring rates are likely a function of both the (volume) melt fraction and the chemistry and kinetics of crystallizing phases (e.g., Clarke et al., 2022). To account for some of these factors, we explore three scenarios for P_armor: constant armoring probability, and either linearly increasing with time or more steeply increasing proportional to the cube of time. In the latter two scenarios, the increasing armoring probability as the system approaches the solidus tests how zircon populations evolve as the sensitivity to melt fraction increases.

For simplicity, we assume in all simulations that all zircon grains are autocrystic grains that do not contain antecrystic or xenocrystic components (cf. Miller et al., 2007; Gaynor et al., 2022), and that zircon grains grow continuously while in contact with melt and that resorption never occurs [i.e., G(t) > 0]. Further, we assume that zircon grains have not experienced Pb loss, which can produce additional age dispersion that, particularly in Mesozoic- and Cenozoic-aged zircon, can be difficult or impossible to identify based on age discordance. However, we note that in the relatively young systems most frequently used to study the time scales of igneous processes, resolvable Pb loss can be mitigated, or eliminated, using chemical abrasion (Mattinson, 2005; von Quadt et al., 2014; Watts et al., 2016; Widmann et al., 2019). Finally, we do not consider issues of systematic age biases that may occur when comparing ages between different methods or isotopic systems, or from comparisons between laboratories. Similarly, we assume that analytical uncertainties are appropriately estimated, as can be empirically demonstrated through repeated analyses of isotopic standards and comparably aged natural reference materials. Since each of these assumptions represents an additional source of uncertainty that can complicate age interpretations, our conclusions represent idealized scenarios that can provide intuition for how to interpret age data. We caution the reader that care should be taken to thoroughly evaluate each of these potential complications when interpreting real data.

All U-Pb zircon geochronology techniques integrate the U and Pb from an analyzed volume to measure ²⁰⁶Pb/²³⁸U and ²⁰⁷Pb/²³⁵U ratios. Thus, it is important to evaluate the significance of the volume-averaged dates produced by various techniques. Several variables will affect volume-averaged zircon dates: the volume sampled, the rate at which the grain grew [G(t)], and, potentially, variations in U content between different growth zones.

In many CA-ID-TIMS studies, following chemical abrasion, single zircon grains are completely dissolved and dated. In these studies, each measured date reflects r̄, the volumetrically integrated average age. If P_armor = 0 (implying that all zircon continued to grow until the system cooled to the solidus), then the possible mean ages for many different zircon growth histories are strongly biased toward the second half of the zircon growth interval (Fig. 1; also see Curry et al., 2021). Zircon grains with mean ages that record the early portion of zircon crystallization are only expected for early crystallizing zircon in systems with rapid early zircon growth, or when P_amor ≫ 0 (see discussion in subsequent section).

While typical CA-ID-TIMS analyses represent integrated whole grain ages, microbeam methods such as LA-ICP-MS and SIMS can analyze isolated regions within a single grain. Isolation of regions of zircon through mechanical breaking of zircon is also increasingly common in CA-ID-TIMS studies (e.g., Samperton et al., 2015). To model the ages derived from subsampling individual grains, we generated zircon with constant, linearly increasing, and linearly decreasing G(t) (Figs. 2G–2I) and calculated the integrated ages of analytical volumes typical of microbeam analyses conducted by LA-ICP-MS and SIMS, as well as volumes comparable to zircon fragments that could be analyzed by CA-ID-TIMS.

As the analytical volume increases relative to the size of the zircon grain, the measured age inevitably converges to the integrated whole grain age (Fig. 2). Specifically, we find that as the radii of LA-ICP-MS analysis spots increase above ~25% of the grain’s minor axis, the maximum possible measured age dispersion decreases to ≤40% of the total zircon growth history, and this result is consistent independent of G(t) (Figs. 2A–2C). Thus, in many studies, typical laser spot sizes represent a significant limitation on the ability to measure intragrain zircon growth duration in small- to medium-sized zircon using LA-ICP-MS, as any zircon with a minor axis ≤ 80 μm would require spots with radii ≤ 20 μm to resolve the majority of intragrain growth. More typical 30–35-μm-diameter spots, which represent up to ~50% of the minor axis for many zircon grains, will produce ages similar to the integrated whole grain age, regardless of whether the core or rim of a grain is analyzed. Laser-ablation line scans, rasters, or depth profiling may better isolate individual growth domains but would require carefully designed and implemented analytical and standardization routines (e.g., Woodhead et al., 2004; Steely et al., 2014; Marsh and Stockli, 2015; Wall et al., 2021).

The shallow depths of typical analytical volumes sampled by SIMS yield improved age resolution, particularly for analyses of zircon cores. Integrated dates from shallow pits in zircon cores diverge significantly from the mean whole grain age, even as the radius of the spot approaches the grain’s minor axis radius (Figs. 2A–2C). However, little additional age resolution is achieved using shallow pits near zircon rims, as spots intersect a range of thin, concentrically grown age domains regardless of pit depth. Analyses of zircon exterior surfaces may be one way to overcome this limitation (e.g., Matthews et al., 2015). These conclusions are generally applicable to a reasonable range of zircon aspect ratios (2–10; Fig. S1). As the aspect ratio increases, the ability to resolve near-rim ages increases for a given spot size relative to the zircon short axis. However, for typical spot sizes, the simulated grains with high aspect ratios will also be quite large.

The integrated ages produced by analyses of zircon fragments also capture a limited age range (Figs. 2D–2F). Regardless of G(t), the integrated age of a zircon rim fragment will be closer to the mean grain age than the final crystallization age for any fragment that is broken more than ~10% from the grain edge. Again, this represents a significant limitation on real-world analyses, as the necessary fragment size for typical zircon will approach the limits of material sizes that can be reliably manipulated in a laboratory (10–20 μm). If a zircon is broken into only two fragments (i.e., a “rim” piece and the remainder of the grain are both analyzed), then the larger fragment inevitably differs little from the mean grain age (Figs. 2D–2F). However, even efforts to isolate the core of zircon by breaking off both grain tips result in integrated ages that are much closer to the mean zircon age than the age of the onset of zircon crystallization (Figs. 2D–2F), as breaking off only the grain tips does not isolate the zircon core from the concentric overgrowths. Recent efforts to mill small flakes of zircon and baddeleyite using a focused ion beam and date them by thermal ionization mass spectrometry (TIMS) (Kovacs et al., 2020; White et al., 2020) may improve the ability to isolate zircon cores. Finally, in contrast to microbeam techniques, the ability to resolve intragrain dispersion from dissolution of zircon fragments does not vary with zircon aspect ratio (Fig. S1). Based on these findings, well-resolved age differences between analyses of asymmetrically fragmented zircons (e.g., Samperton et al., 2015) may indicate even longer single-crystal growth intervals, but they may also highlight the importance of more complex nonconcentric zircon growth, which is not considered here. Detailed studies of zircon megacrysts using a combination of LA-ICP-MS analyses and CA-ID-TIMS analyses of grain fragments can also be used to test these models (e.g., Schaltegger et al., 2015). However, as these megacrysts are exclusively found in unusual magmatic systems, including kimberlites and strongly alkaline pegmatites, it is unclear how closely their formation mirrors zircon crystallized from more typical magmas.

The mass of uranium incorporated in zircon is controlled by the U concentration in the melt and the melt-zircon U distribution coefficient, which is sensitive to temperature, melt composition, and fO₂ (Claiborne et al., 2018). Therefore, as zircon crystallizes from a cooling and differentiating melt, it may retain uranium compositional gradients that can influence integrated zircon ages. To evaluate this possibility, we repeated the calculations shown in Figure 2 for a zircon that grows at a constant rate but incorporates either linearly increasing or linearly decreasing U over its growth history (Fig. 3). These calculations show that for reasonable total core-rim variations in U content (total variation considered in Figure 3 is two orders of magnitude), varying U concentrations during zircon growth produces changes in whole grain mean age but has only a minor effect on the distribution of core and rim ages around the mean age in comparison to the constant U case.

The calculations presented in Figures 2 and 3 and discussed in the preceding paragraphs represent an idealized scenario that treats analytical uncertainties as negligible. However, real analyses have nonzero uncertainties, and the likelihood of measuring resolvable age differences within a single grain will be limited as the magnitude of the uncertainty approaches or exceeds the duration of grain growth. We find that in all modeled scenarios, measuring resolvable zircon growth duration from a single pair of core and rim analyses is only likely when the analytical uncertainty is <30% of the zircon growth duration (Fig. 4A), and that this constraint grows more severe as the radius of the analytical volume increases.

Figure 4 also highlights the need for caution when interpreting zircon growth durations based on paired core and rim analyses (e.g., Wu et al., 2022; Zimmermann et al., 2018). If analytical uncertainties are assumed to be normally distributed, then even when the analytical uncertainty is equal to the true duration, there is an ~10% chance of producing core and rim ages with resolvable dates (Fig. 4A). However, in these 10% of cases, the measured duration results from measurements that are outliers at the two standard deviation level and has no relation to the actual growth duration. The “resolvable” durations calculated in these 10% of measurements are also typically more than three times greater than the total true zircon duration (Fig. 4B). In general, as measurement uncertainties increase to greater than ~20% of the true duration, the measured duration calculated from any pair of analyses with resolvable dispersion will be strongly biased to longer time scales than the true duration.

Our analysis shows that microbeam methods can reliably resolve intragrain dispersion only if (1) the zircon is significantly larger than the radius of the analytical volume, and (2) the duration of zircon growth is significantly greater than the scale of analytical uncertainty. Further, measurements of intragrain duration based on a single pair of core and rim analyses should be treated with great caution because it is difficult to parse the contributions from analytical uncertainty relative to true age dispersion. Instead, dispersion calculations based on either multiple core and rim analyses of large single zircon grains or paired core-rim analyses of multiple grains that show consistent durations are necessary to convincingly show that the measured intragrain duration is robust. In contrast, more precise zircon U-Pb analyses by CA-ID-TIMS have the potential to characterize intragrain growth. However, without more sophisticated microsampling techniques (e.g., Kovacs et al., 2020; White et al., 2020), these dates are inevitably biased toward the zircon integrated mean age.

The preceding section demonstrates that, in many instances, it is difficult to accurately measure the duration over which an individual zircon grew, except for microbeam measurements where analytical volumes are significantly smaller than the total zircon volume and the analytical uncertainties are at least a factor of ~5 shorter than the “true” duration. Therefore, in most instances, characterizing the duration of zircon crystallization within an individual sample will require making inferences from the distribution of volumetrically averaged zircon ages measured from multiple grains. Systems amenable to accurately determining crystallization durations with this approach need to have: (1) protracted zircon crystallization, where zircon saturation is reached early in the crystallization history; and (2) punctuated zircon growth, where the bulk of growth for individual grains occurs over a restricted duration relative to the total duration that the magma spent between zircon saturation and the end of zircon crystallization (e.g., solidus or quench). However, magmas that were emplaced below zircon saturation temperature may result in a greater abundance of antecrystic and/or xenocrystic zircon (Miller et al., 2007), which will obfuscate the autocrystic crystallization history. Zircon saturation is controlled by magma temperature and composition (Watson and Harrison, 1983; Boehnke et al., 2013), and we do not explicitly consider these controls here and, as mentioned above, assume that only autocrystic zircon grains are present. Therefore, we focus only on the fidelity with which a zircon age distribution reflects the total duration of crystallization from zircon saturation to the solidus as a function of the three variables G(t), J(t), and P_armor.

The mean zircon age distributions produced when P_armor = 1 and P_armor = 0 represent two end-member scenarios (Fig. 5). When P_armor = 1, zircon growth is instantaneous, and the zircon age distribution reflects the distribution of zircon nucleation [J(t)]. In contrast, when there is no armoring, all grains continue to grow until the system cools to the solidus, and the distribution of zircon mean ages can be determined according to Equation 5. The distribution of ages produced in this scenario is a function of both the zircon growth and nucleation curves as described with Equation 5. This age distribution is inevitably biased toward the time when the system cools through the solidus, with the oldest possible date limited by the mean age of grains that nucleate at onset of zircon saturation (Figs. 1 and 5).

The nucleation and growth curves considered in Figure 5 are relatively simple and may not be representative of real magmatic systems. However, even these simple curves highlight the complexity in age distributions that can be generated through combinations of nucleation and growth processes. Critically, these calculations show that many systems characterized by similar distributions of total crystallized zircon mass can have significantly different zircon age distributions. Similarly, many simple growth rates only have the effect of scaling the zircon age distribution relative to the total growth duration (see, e.g., the decreasing, and uniform growth rates for each nucleation rate) and would therefore be impossible to distinguish in natural data without additional context. Thus, great care is needed when employing modeling approaches that interrogate zircon age distributions without explicitly considering nucleation and growth (e.g., Keller et al., 2018; Caricchi et al., 2014).

Observations from many natural systems show that zircon armoring is an incomplete process: Zircon grains are commonly found armored by other silicate phases but are also observed at grain boundaries/interstices and/or in volcanic glass (Pamukçu et al., 2013). To investigate the impacts of incomplete armoring on zircon mean age, we simulate the zircon mean age distributions of 1,000,000 synthetic zircon grains generated by each of the G(t) and J(t) combinations shown in Figures 5A and 5D, in combination with three different P_armor functions: no armoring, constant armoring probability, and a cubic increasing armoring probability (Fig. 6). Each of these armoring probabilities is defined so that a zircon crystallized at t = 0 in any simulation has equal cumulative armoring probability.

These simulations show that changes to armoring rate can significantly modify the distribution of zircon mean ages. The exact behavior depends on the specific parameters, but there are two broadly consistent effects: (1) an increase in the number of zircon grains with mean ages from early in the crystallization interval, and (2) a shift in the mean zircon age and the shape of the age distribution. The first effect is well illustrated by many of the simulations with increasing growth rates. Zircon mean ages dominantly record the second half of the crystallization interval without armoring but record the entire crystallization interval when P_armor is greater than 0. However, zircon grains recording the earliest period of crystallization are typically still rare. Similarly, the shapes of the age distributions change significantly with armoring: Flat uniform distributions without armoring shift to positively sloped distributions with armoring, while the characteristic shape of beta distributions evolves to more Gaussian-like distributions with increasing armoring rates. These shifts highlight the potential of armoring to further complicate the interpretation of magma dynamics from zircon age distributions in terms of only the distribution of zircon mass crystallized.

In reality, zircon armoring is likely a complex process that depends on numerous variables not considered here, including nucleation processes, the evolution of the major crystallizing phases along a liquid line of descent, the nonlinear relationship between temperature and melt fraction in a crystallizing magma, and magma dynamics. The importance of processes including synneusis (Vance, 1969), heterogeneous nucleation of either zircon or the entraining phase (Hammer et al., 2010), or zircon saturation in boundary layers formed during growth of major phases (Bacon, 1989) in relation to zircon armoring efficiency remains to be determined and may have nonintuitive impacts on zircon age distributions. Regardless, our simulations show that better understanding of zircon armoring and nucleation will be critical to using zircon age distributions to understand magmatic processes. High-precision geochronologic studies of contextually constrained zircon (e.g., Barboni and Schoene, 2014; Chambers et al., 2020), perhaps in combination with crystal size distributions and trace-element compositions (e.g., TIMS with trace-element analysis [TIMS-TEA]; Schoene et al., 2010), are an important next step in evaluating these variables. These new data will provide important constraints for future model developments that aim to directly predict the physical and chemical processes controlling zircon crystallization (cf. Caricchi et al., 2014; Ratschbacher et al., 2018; Sorokin et al., 2022).

Most high-precision CA-ID-TIMS studies analyze ≤10 grains per sample. If we assume that the selection of zircon was not biased based on zircon appearance or other criteria, the analyzed grains comprise a random sample of the true zircon distribution, and the measured age distribution will be a combination of the zircon mean age distribution and the analytical uncertainty. Therefore, the ability to make meaningful inferences about the true zircon distribution (and its underlying controls) depends on both the number of zircon grains analyzed and the relative scales of the analytical uncertainties and the zircon age distributions.

To illustrate these dynamics, we focus on one of the first questions an investigator might ask: “Do the zircon grains in this sample record dispersion that cannot be attributed to analytical uncertainty?” Answering this question is a necessary first step, as a distribution of ages that can be explained by analytical uncertainty alone provides no additional information about zircon crystallization duration, and further investigations into the zircon age distribution will be fruitless. Geochronologists have a widely used tool for answering this question: the mean square of weighted deviates (MSWD; Wendt and Carl, 1991), which, for single-dimensional data like ²⁰⁶Pb/²³⁸U age populations (i.e., as opposed to two-dimensional isochrons), can be thought of as a test for the normality of the sample population. If the MSWD value for a sample falls between calculated thresholds [with a minimum defined as and a maximum threshold of , but often assumed to be simply MSWD near ~1], then the observed variance between measurements can be fully explained by the analytical uncertainties (i.e., the sample distribution is consistent with a normal distribution with variance defined by measurement uncertainties). By contrast, samples significantly below the minimum MSWD threshold (often assumed as MSWD ≪ 1) have significantly less dispersion than is expected from the calculated uncertainties. This result implies that either the analytical uncertainties are overestimated, or the sample analyses are inappropriately filtered, introducing a bias (e.g., removing more scattered analyses). In contrast, when the MSWD is above the maximum threshold (typically assumed as MSWD > 1), the dispersion in the data cannot be fully explained by analytical uncertainty, indicating that either the analytical uncertainty is underestimated, or there is additional nonanalytical scatter in the data. While this additional scatter can be generated by several factors, including zircon inheritance and Pb loss, we focus on the situation where both of those possibilities have been eliminated, and any additional scatter is attributed only to an extended duration of crystallization.

To simulate this scenario, we tested a series of models where we randomly simulated 100,000 samples of zircon measurements with a sample size n. The age distributions were generated using one of three underlying zircon age distributions: a standard uniform distribution, a beta distribution, and a truncated exponential distribution. An additional Gaussian error was added to each measurement, with the standard deviation of the error scaled relative to the true duration (Figs. 7A–7C). We calculated the MSWD of each population and found the percent of samples that exhibited overdispersion. As expected, regardless of the sample size and the underlying age distribution, when the analytical uncertainty is comparable to the true zircon duration, the dispersion observed in nearly all populations is well explained by analytical uncertainty (Figs. 7D–7F). This effect is well illustrated in Figures 7A to 7C: When the analytical uncertainties are equal to the true zircon duration, they overwhelm the underlying zircon distribution, and the distribution of measured ages is indistinguishable from a normal distribution.

Alternatively, if the analytical uncertainty is smaller than the true duration, the underlying distribution is more evident in the measurement distribution (Figs. 7A–7C). For a uniform underlying distribution, which is the most easily resolved of the three distributions considered, when the true duration is twice the analytical uncertainty, even a 50 zircon sample has a <40% probability of displaying overdispersion (Fig. 7D), and a 95% probability of identifying overdispersion is only reached when the true duration is roughly four times longer than the uncertainty. Using a more typical 10 zircon sample size, 95% of overdispersed samples are identified only once the true duration is roughly six times longer than measurement uncertainties. The probability of resolving overdispersed samples from the beta distribution is worse (Fig. 7E): For a 10 zircon sample, there is a 95% probability of identifying overdispersion only when the analytical uncertainties are ≤7% of the true duration. Probabilities for samples drawn from the truncated exponential distribution are intermediate between these two distributions. Thus, identification of overdispersion, regardless of underlying age distribution, requires either an extremely large n or a system in which the duration of zircon crystallization was much longer than the analytical uncertainty. The systematics of the underlying distribution also influence how easily overdispersion is identified: Overdispersion becomes more difficult to identify when the underlying distribution approaches a normal distribution or is fat-tailed, while it becomes more readily identified when the underlying distribution is less similar to a Gaussian distribution, with a broad mode, or it is asymmetric.

The uncertainties in Figures 7C to 7F are scaled relative to the true duration of zircon ages. However, when analyzing zircon from an unknown sample, a researcher can only know the measured duration of the analyzed population. When uncertainties are significantly smaller than the true duration, this measured duration necessarily underestimates the true duration (Glazner and Sadler, 2016). However, as the uncertainties increase, the measured duration is increasingly controlled by the analytical uncertainties. Importantly, when uncertainties are small relative to true duration, the measured duration will approach the true duration as sample size increases, but when uncertainties approach the true duration, the measured duration increases without bound with increasing sample size, as the likelihood of multiple-standard-deviation outlier measurements also increases. We show how this dynamic modifies the likelihood of correctly identifying dispersed ages in Figures 7G to 7I. Here, we use the same approach as in Figures 7C to 7F, but the x axis shows the ratio of the analytical uncertainties to the average duration measured for the simulated samples. This presentation has two initially unexpected results. First, the effect on large sample sizes and small sample sizes is inverted: The probability of identifying overdispersed data decreases for large n samples but increases for small n samples. This sensitivity can be understood as a result of the described dynamics: Because it is much less likely to observe a multiple-standard-deviation outlier measurement in a small n sample, significant overdispersion in small n samples is more likely to be reflective of the underlying distribution. Second, the differences between the different underlying distributions observed in Figures 7D to 7F diminish greatly, and most curves for n > 3 converge for all underlying distributions and steeply increase when uncertainty/measured duration is <~0.20. This behavior points to a relatively simple rule of thumb: When the difference between the youngest and oldest measured grains is offset by roughly five times the standard deviation, it is highly probable that this spread reflects overdispersed data. However, as stated previously, this approach has been developed with the assumption that the measured ages are produced only from autocrystic grains, an assumption that needs to be evaluated before applying this rule of thumb to real data.

This exercise highlights the need for high-precision geochronology to study short-duration magmatic events. To reliably identify overdispersion in a 10 zircon population, the duration needs to be at least five times the analytical uncertainty. Thus, even in a relatively young system (ca. 20 Ma) and with 0.05% precision at the one standard deviation level, only magmatic events with zircon crystallization intervals lasting at least ~50 k.y. will be accessible. This conclusion will vary somewhat depending on the true zircon age distribution, but many tailed distributions will impose even stricter requirements.

For samples where the distribution of zircon ages is overdispersed and can be attributed to the real duration of zircon crystallization, a simple way to report both the age of the sample and the crystallization duration is desirable. This reporting should be as generalizable as possible, because, as we have shown, it is typically difficult to determine the zircon age distribution within a sample. Therefore, we advocate for an approach that does not make any assumptions about the underlying distribution.

Given these constraints, we believe that the best approach is to simply report the standard mean age and to include the age dispersion greater than that attributed to analytical uncertainties alone. This dispersion places a lower bound on the true crystallization duration, regardless of the underlying distribution. Assuming the analytical uncertainties are normally distributed, the contribution to the measured dispersion from uncertainties can be estimated as the expected value (i.e., mean) of the range of a sample set drawn from a normal distribution, and the “true” measured duration can be calculated by subtracting this value from the measured duration. Although there is no closed-form solution to this statistic, values for a given n can be found through numerical integration or Monte Carlo simulation, and values have been tabulated up to large n (Tippett, 1925). For simplicity, previous workers have presented relatively simple equations that approximate this value for a range of sample sizes (e.g., Tukey, 1955; Walter et al., 2022). Here, we simulate the mean of the range for 2–50 samples and fit a power-law equation that is accurate to within 0.1% of the calculated values for n between 4 and 40, spanning the relevant range of sample sizes for nearly all high-precision geochronologic studies (Fig. 8A):

Equation 7 estimates the size of the expected range due to measurement uncertainties alone as a multiple of the standard deviation of the error for a sample size of n zircon grains. As such, it is closely related to the studentized range statistic, but it differs in that it is calculated based on the population standard deviation (as estimated by the measurement errors) and not the sample standard deviation. Equation 7 and Figure 8A show that the sample range from a normal distribution increases rapidly as the sample size increases from ~2 to 10 analyses and continues to increase monotonically with increasing sample size. This increase reflects the increasing likelihood of multiple-standard-deviation outlier measurements as the number of analyses increases.

When the MSWD of a sample is less than or equal to the maximum threshold, it is not possible to estimate the minimum duration, because the scatter in the data is consistent with an instantaneous event overprinted by analytical uncertainty. However, these data still provide constraints on the duration of the event, because it is possible to place an upper bound on the total duration: If the zircon crystallization was not instantaneous, then the probability of a measured sample having an MSWD less than the maximum threshold will decrease as the sample size increases and as the crystallization duration increases relative to the analytical uncertainty (Fig. 7A). Based on this rationale, we again can use Monte Carlo simulations to create samples of varying size from a uniform distribution with an added Gaussian uncertainty. Further, in each simulation, the standard deviation of the uncertainty is varied relative to the total duration. Using these simulations, we calculate for a given n the ratio of the standard deviation of the analytical uncertainty to the total duration at which 95% of the simulations return MSWD values above the maximum threshold (Fig. 8B), which provides a robust constraint on the maximum duration. We then fit a power-law equation to these data as a function of the sample size:

This equation calculates the 95th percentile of the maximum duration relative to the measured analytical uncertainty. However, unlike estimating the minimum duration using Equations 7 and 8, this equation is dependent on the zircon age distribution, and therefore some care is required. For most nonuniform distributions, Equation 9 will underestimate the maximum duration, and for long-tailed distributions, the maximum duration can be an underestimated by more than a factor of 2. Thus, this approach will likely need to be revisited as better constraints on zircon age distributions are obtained.

Equations 7 and 9 were both calculated from synthetic data where the standard deviation of the added error is constant within each simulation. However, in real data sets, uncertainties vary between individual analyses due to a range of factors, including changing proportions of common and radiogenic Pb, signal intensity, and instrument performance. Thus, it is necessary to calculate the average measurement analytical uncertainty to properly scale the results obtained using Equations 7 and 9. This value is distinct from the commonly reported standard deviation of the mean age and is calculated as the root mean square (i.e., quadratic mean) of the analytical uncertainties of the individual measurements. This adjustment compared to the synthetic data with constant analytical uncertainty has a small effect for typical well-behaved zircon data sets, where analytical uncertainties show limited variation (e.g., a factor of 2 or 3). However, we advise caution when attempting to interrogate crystallization duration from data sets with highly variable analytical uncertainties, as single grains with unusually large analytical uncertainties can significantly distort the conclusions obtained using Equation 7 or 9.

In Figure 9, we use three examples of recently published CA-ID-TIMS zircon data from single samples of Mesozoic and Cenozoic plutons to illustrate these constraints on zircon crystallization duration. The data set shown in Figure 9A is for the reference material GHR1 (Eddy et al., 2019). Zircon ages from this sample are tightly clustered, yielding an MSWD of 1.52. As this is below the calculated maximum MSWD threshold (1.58 for a sample size of 25 analyses), these data are consistent with instantaneous emplacement and crystallization, with the observed dispersion consistent with the dispersion expected from analytical uncertainty alone. This is further evidenced by the expected value of the duration from analytical errors alone, which is calculated using Equation 7 and indicated with the gray band in Figure 9A. This duration is slightly greater than the observed scatter in the measured data. Although a minimum duration cannot be constrained, we can estimate an upper bound on the duration using Equation 9. Based on this equation and the reported measurement uncertainties, the maximum crystallization duration is estimated to be 108 k.y.

In contrast to Figure 9A, Figures 9B and 9C show two samples with calculated MSWD values above the maximum threshold. Thus, the dispersion in both samples is greater than the estimated range from analytical uncertainties calculated using Equations 7 and 8, and a minimum duration can be estimated. In Figure 9B, the data are only slightly overdispersed relative to the expected analytical range, and the calculated minimum duration is relatively short, at 33 k.y. In contrast, the data in Figure 9C are more significantly overdispersed, and a longer minimum duration of 102 k.y. is estimated. We emphasize here that the exact location of the crystallization duration calculated using Equation 8 relative to the weighted mean age is not constrained without knowledge of the underlying zircon age distribution. (If we consider the zircon age distributions shown in Figure 5 and Figure 6, the location of the mean or expected value zircon age varies within the zircon age range for each distribution.) In Figure 9, we show our estimated range centered on the weighted mean ages only for illustration purposes. However, this presentation should not be understood to imply that the analytical dispersion is centered on the weighted mean and that the excess duration is found in the tails of the data. Instead, it is most likely that the excess crystallization duration is located within the bulk of the data, and that the data near the tails relate to the analytical dispersion.

These results present a simple addition to existing standard data reporting to address the case where data are overdispersed. Even for overdispersed zircon populations, it is often useful for the community to have a single “date” to use, for example, in tectonic reconstructions or when relating pluton emplacement to significantly longer-scale processes. In these cases, particularly when a high-precision date is not necessary, we argue that it is still best to use the weighted mean age, as this is an age that has an intuitive meaning and is familiar to most researchers. Further, we suggest that the standard data description of weighted mean age ± uncertainties (using the x/y/z uncertainties format; see, e.g., McLean et al., 2011), sample size, and MSWD should be supplemented by the minimum or maximum duration calculated using Equation 8 or 9. The duration should not be reported as an additional error but noted separately so as not to imply that the duration is centered on the mean age, but simply that the crystallization duration includes the mean age.

The main advantages of our proposed approach to constraining crystallization duration are its simple implementation and the minimal number of underlying assumptions. Recently, two additional approaches have also been proposed for estimating crystallization duration: one described by Vermeesch (2018) and implemented in the widely used open-source data reduction software IsoplotR, and a second approach described by Ratschbacher et al. (2018). The IsoplotR approach shares some similarities with our method, with one critical difference. Our approach estimates the duration expected for a given sample size and analytical uncertainty and attributes additional observed duration to “real” overdispersion. IsoplotR aims to determine how much additional dispersion is required to produce normally distributed data that are consistent with the combined uncertainties and observed dispersion. This approach implicitly assumes that the observed age distribution can be described by the sum of two Gaussian distributions, one constrained from analytical uncertainties and a second that is attributed to “geologic scatter” (i.e., age dispersion). The assumption that the true zircon age distribution follows a Gaussian distribution is the key distinction between our approaches, and one that our population modeling (Fig. 5) does not support. By contrast, the method proposed by Ratschbacher et al. (2018) takes a very different approach: It models a magma’s liquid line of descent using AlphaMELTS, and from this calculation, it uses the predicted relationships among temperature, melt composition, and zircon mass crystallized as priors for Bayesian statistical approaches that find the time of zircon saturation and solidus most consistent with the priors and the measured zircon age distributions.

In Figure 10, we show the predictions of each of these three approaches for the data originally analyzed in Ratschbacher et al. (2018). All five of these samples have MSWD values above the maximum threshold, and thus it is appropriate to apply our minimum duration estimate using Equations 7 and 8. In each of the five samples, the duration estimated using the Ratschbacher et al. (2018) approach is significantly larger than the estimates produced using either our technique or the IsoplotR method. Further, for all five samples, the Ratschbacher et al. (2018) crystallization duration exceeds the measured duration. In contrast, our approach and the IsoplotR method produce more conservative estimates for the minimum duration and average duration, respectively. However, our approach consistently produces longer durations than IsoplotR due to the different assumptions regarding the underlying distribution. Figure 10 highlights the different duration estimates that can be produced based on the approach. Our new method and that of IsoplotR both take similar, conservative and descriptive approaches and seek to understand how much of the observed dispersion is not explained by the analytical uncertainties. In contrast, the approach presented by Ratschbacher et al. (2018) aims to estimate the true duration, but it is built upon several currently unsupported assumptions. First, it requires that the magma composition from which the sample crystallized is known, which is challenging in plutonic rocks where whole-rock analyses frequently do not provide a liquid composition (e.g., Gelman et al., 2014; Barnes et al., 2016; Cornet et al., 2022). Second, it directly equates zircon age distributions to zircon mass crystallization distributions, a relationship that we have shown is not justified for integrated, whole-grain zircon ages unless P_armor = 1, representing an extreme end-member scenario unlikely to occur in natural systems (Figs. 1, 5, and 6). Ultimately, approaches similar to that presented by Ratschbacher et al. (2018) will be necessary to constrain zircon crystallization duration because more conservative duration estimates that make no assumption about the underlying age distribution (like the one presented in this study) will inevitably only constrain the minimum duration. However, the assumptions that go into these approaches need to be carefully tested and validated before they can be confidently applied to real systems.

In this contribution, we developed a theoretical framework for relating zircon growth, nucleation, and armoring rates to zircon ages and to age distributions within zircon populations. Using this theoretical framework, we showed that there are limitations to all typical analytical approaches. Whole-grain ages produced by CA-ID-TIMS studies are inevitably biased toward the later portion of the zircon crystallization interval. In contrast, subsampling methods, either through the analysis of zircon fragments or microbeam methods, can only accurately measure intragrain zircon dispersion when the analytical volume is significantly smaller than the zircon grain size and when the analytical precision is significantly smaller than the duration time scale. Given these limitations to microbeam methods, we focused on the interpretation of zircon mean age populations. We showed that the rate of zircon mass crystallization does not uniquely define the distributions of zircon ages and that, instead, careful consideration of zircon nucleation, growth, and armoring rates is required. Finally, we showed that despite these challenges, zircon age populations can provide robust constraints on either the minimum or maximum crystallization duration.

Through this work, we have been occasionally disheartened by the conclusions we reached. Our modeling shows that interrogating magmatic processes from zircon ages is difficult even in ideal circumstances and that zircon growth and/or nucleation rates are very poorly resolved by relatively small n geochronologic data sets. Robust constraints from typical zircon data sets are relatively weak and have limited power in constraining magmatic processes. However, rather than discouraging researchers, we hope that this study will motivate more careful and creative approaches to constraining the time scales of magmatic processes. While laboratory advances resulting in smaller analytical errors and increasing analytical data set sizes will always improve the quality and power of zircon data sets, we argue that there is a wealth of nongeochronologic data that can be leveraged to further constrain magmatic processes. Models that couple zircon geochronologic data with additional petrologic data, including trace-element data and magmatic thermometers, can potentially better constrain underlying zircon age distributions. Similarly, interpretation of high-precision zircon geochronology within a textural context and interrogation of data from multiple samples within a well-constrained microstructural framework can improve constraints on the thermal and magmatic evolution of a system.

We thank J.-A. Olive for helpful discussion regarding the model. Discussion with O. Müntener and constructive feedback from two anonymous reviewers and editor B. Singer significantly improved this manuscript. B.Z. Klein acknowledges funding support from the Université de Lausanne. M.P. Eddy was supported through internal funding provided by Purdue University.