Modeling apparent Pb loss in zircon U–Pb geochronology

The loss of radiogenic Pb from zircon is known to be a major factor that can cause inaccuracy in the U–Pb geochronological system; hence, there is a need to better characterize the distribution of Pb loss in natural samples. Treatment of zircon by chemical abrasion (CA) has become standard practice in isotope dilution–thermal ionization mass spectrometry (ID-TIMS), but CA is much less commonly employed prior to in situ analysis via laser ablation–inductively coupled plasma–mass spectrometry (LA-ICP-MS) or secondary ionization mass spectrometry (SIMS). Differentiating the effects of low levels of Pb loss in Phanerozoic zircon with relatively low-precision in situ U–Pb dates, where the degree of Pb loss is insufficient to cause discernible discordance, is challenging. We show that U–Pb isotopic ratios that have been perturbed by Pb loss may be modeled by convolving a Gaussian distribution that represents random variations from the true isotopic value stemming from analytical uncertainty with a distribution that characterizes Pb loss. We apply this mathematical framework to model the distribution of apparent Pb loss in 10 igneous samples that have both non-CA LA-ICP-MS or SIMS U–Pb dates and an estimate of the crystallization age, either through CA U–Pb or ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ geochronology. All but one sample showed negative age offsets that were unlikely to have been drawn from an unperturbed U–Pb date distribution. Modeling apparent Pb loss using the logit–normal distribution produced good fits with all 10 samples and showed two contrasting patterns in apparent Pb loss; samples where most zircon U–Pb dates undergo a bulk shift and samples where most zircon U–Pb dates exhibited a low age offset but fewer dates had more significant offset. Our modeling framework allows comparison of relative degrees of apparent Pb loss between samples of different age, with the first and second Wasserstein distances providing useful estimates of the total magnitude of apparent Pb loss. Given that the large majority of in situ U–Pb dates are acquired without the CA treatment, this study highlights a pressing need for improved characterization of apparent Pb-loss distributions in natural samples to aid in interpreting non-CA in situ U–Pb data and to guide future data collection strategies.

1 Introduction

Zircon U–Pb geochronology is arguably one of the most important radiometric dating approaches used by geoscientists, with widespread application to constraining the age of Pleistocene and older geologic materials (Davis et al., 2003; Schoene, 2013; Gehrels, 2014). We rely on zircon U–Pb dates for calibrating the geological timescale (e.g., Compston, 2000a, b; Bowring and Schmitz, 2003; Gradstein et al., 2004; Kaufmann, 2006), constraining the timing of important Earth history events (Froude et al., 1983; Burgess et al., 2014), and determining the rates of Earth processes (Rioux et al., 2012; Johnstone et al., 2019). The zircon U–Pb geochronometer is particularly powerful due to the ability to assess agreement between the ²³⁸U → ²⁰⁶Pb and ²³⁵U → ²⁰⁷Pb decay chains, with ${}^{\mathrm{206}}{\mathrm{Pb}}^{\ast }{/}^{\mathrm{238}}\mathrm{U}$ and ${}^{\mathrm{207}}{\mathrm{Pb}}^{\ast }{/}^{\mathrm{235}}\mathrm{U}$ dates in agreement plotting on the concordia line, where ∗ indicates radiogenic Pb (Wetherill, 1956).

The causes and complications of open-system behavior (e.g., radiogenic Pb loss) in zircon have long been a topic of study (Tilton et al., 1955; Pidgeon et al., 1966). Although Pb-loss events may be discerned on U–Pb concordia diagrams in some circumstances and can provide useful geologic information about the thermal and/or fluid flow history of a region (Silver and Deutsch, 1963; Blackburn et al., 2011; Morris et al., 2015; Kirkland et al., 2017), recognizing that Pb loss remains a challenge when it occurs within several hundreds of millions of years of crystallization (Fig. 1; Andersen et al., 2019). For example, due to the shape of the ${}^{\mathrm{206}}{\mathrm{Pb}}^{\ast }{/}^{\mathrm{238}}\mathrm{U}$ versus ${}^{\mathrm{207}}{\mathrm{Pb}}^{\ast }{/}^{\mathrm{235}}\mathrm{U}$ concordia line, Pb loss in Phanerozoic zircon results in a “sliding along concordia” effect that can make Pb loss difficult to discern, particularly in relatively low-precision in situ (i.e., LA-ICP-MS or SIMS) datasets when the Pb loss produces concordant or only modestly discordant analyses (e.g., <10 %; Ashwal et al., 1999; Bowring and Schmitz, 2003; Ireland and Williams, 2003; Reimink et al., 2016; Spencer et al., 2016; Watts et al., 2016; Andersen et al., 2019). Such low levels of Pb loss have been termed “cryptic” and may be associated with spatial heterogeneities, including radiation-damaged U-rich zones and microstructures (Nasdala et al., 2005; Kryza et al., 2012; Watts et al., 2016). Although there are many potential causes of Pb loss in zircon, open-system behavior is often associated with elevated α-dose and associated metamictization (Silver and Deutsch, 1963; Pidgeon et al., 1966; Mezger and Krogstad, 1997; Cherniak and Watson, 2001; Marsellos and Garver, 2010). Mechanisms for Pb loss include recrystallization of metamict zircon during metamorphism (Kröner et al., 1994; Mezger and Krogstad, 1997; Orejana et al., 2015; Zeh et al., 2016) and leaching of Pb from metamict zones by hydrothermal or diagenetic fluids (Geisler et al., 2002, 2003; Willner et al., 2003; Morris et al., 2015; Kirkland et al., 2020) or during chemical weathering (Stern et al., 1966; Black, 1987; Balan et al., 2001; Pidgeon et al., 2017; Andersen and Elburg, 2022). Pb loss is thought to primarily occur at temperatures <250 ^∘C in which radiation damage in zircon is unable to be annealed over geologic timescales (Schoene, 2013).

Figure 1 Illustration of the influence of Pb loss on 250 Ma and 2.5 Ga zircon. Two Pb-loss scenarios are shown, with 25 % loss at half the age of the zircon and 50 % loss at present day (0 Ma). The approximately linear nature of the ${}^{\mathrm{206}}{\mathrm{Pb}}^{*}{/}^{\mathrm{238}}\mathrm{U}$ vs. ${}^{\mathrm{207}}{\mathrm{Pb}}^{*}{/}^{\mathrm{235}}\mathrm{U}$ concordia line near the origin results in Pb loss producing limited discordance if the Pb loss occurs within several hundreds of millions of years of crystallization. Note that a greater amount of ancient Pb loss is required to produce the same shift in ${}^{\mathrm{206}}{\mathrm{Pb}}^{*}{/}^{\mathrm{238}}\mathrm{U}$ relative to recent Pb loss. Thin, colored lines represent the path of each zircon.

Zircon domains that have lost Pb may be preferentially removed by first thermally annealing the zircon at high temperature (e.g., 800–1100 ^∘C) and then partially dissolving the zircon in a heated hydrofluoric acid (HF) solution in a technique called chemical abrasion (CA) (Mattinson, 2005). The CA treatment is now routinely applied in isotope dilution–thermal ionization mass spectrometry (ID-TIMS) analysis and has contributed to both the improved precision and accuracy of CA-ID-TIMS U–Pb data (Schoene, 2013). Although some in situ U–Pb laboratories practice thermal annealing routinely (e.g., Allen and Campbell, 2012; Solari et al., 2015), CA has been applied much less frequently (Crowley et al., 2014; von Quadt et al., 2014; Watts et al., 2016; Ver Hoeve et al., 2018; Ruiz et al., 2022). Several studies that have conducted paired analysis of non-CA and CA of the same samples via in situ U–Pb geochronology have found the non-CA U–Pb dates to skew younger than the CA U–Pb dates (Crowley et al., 2014; von Quadt et al., 2014; Watts et al., 2016). A growing number of maximum depositional age studies with tandem non-CA laser ablation–inductively coupled plasma–mass spectrometry (LA-ICP-MS) and CA-ID-TIMS dating have shown that the youngest non-CA U–Pb dates tend to be younger than expected relative to CA U–Pb dates or other geologic constraints, even when considering measurement uncertainty (e.g., Herriott et al., 2019; Schwartz et al., 2022; Howard et al., 2022; Sharman et al., 2023). However, there is a lack of quantitative constraints on the relative importance of Pb loss in influencing non-CA U–Pb date distributions acquired via in situ mass spectrometry, particularly as related to influencing depositional age constraints (Copeland, 2020).

This study builds upon past research on open-system behavior in zircon by presenting a mathematical framework for characterizing the distribution of apparent Pb loss on untreated (i.e., non-CA) U–Pb date distributions. We first suggest that U–Pb isotopic ratios that have been perturbed by Pb loss may be viewed as the convolution of two signals: a Gaussian distribution that reflects the measurement uncertainty about the true isotopic ratio and the distribution that characterizes Pb loss. We then apply this mathematical framework to model the distribution of apparent Pb loss that has affected 10 igneous samples of Miocene to Carboniferous age. Our results highlight the importance of quantifying distributions of apparent Pb-loss magnitude to better understand the potential influence on non-CA zircon U–Pb date distributions.

2 Mathematical framework

A series of n ${\mathrm{Pb}}^{*}/\mathrm{U}$ measurements that have undergone Pb loss, Z, may be modeled as the sum of the corresponding unperturbed ${\mathrm{Pb}}^{*}/\mathrm{U}$ values, X, and the amount that ${\mathrm{Pb}}^{*}/\mathrm{U}$ changed due to Pb loss for each date, Y,

$$\begin{array}{}\text{(1)}& \mathit{Z}=\mathit{X}+\mathit{Y},\end{array}$$

where Z, X, and Y are all 1-D arrays with n values and units of percentage offset from the true isotopic value (Fig. 2). Because Pb loss produces a lower ${\mathrm{Pb}}^{*}/\mathrm{U}$ ratio, the values of Y must be negative in our formulation of Eq. (1). If X is drawn from a Gaussian distribution f(t), whose mean (μ) approximates the true isotopic value and whose standard deviation (σ) reflects dispersion from the true value related to measurement uncertainty (e.g., Schoene, 2013), and if Y is drawn from a distribution that reflects Pb loss, g(t), then Z may be viewed as being drawn from the convolution of f(t) and g(t)

$$\begin{array}{}\text{(2)}& \left(f\ast g\right)\left(t\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}f\left(\mathit{\tau }\right)g\left(t-\mathit{\tau }\right)\mathrm{d}\mathit{\tau },\end{array}$$

provided that X and Y are independent (Fig. 2; Video S1 in Sharman and Malkowski, 2023). Convolution simply represents the summation of two random variables, in this case one related to analytical precision (i.e., random variation around the true isotopic value stemming from the measurement process) and the other related to the geologic process of Pb loss. We model Pb loss as percentage offset from the true ${\mathrm{Pb}}^{*}/\mathrm{U}$ value rather than deviation in absolute time (i.e., Myr) to promote a comparison of samples of different age (Fig. 2).

Figure 2 Illustration of how ${\mathrm{Pb}}^{*}/\mathrm{U}$ isotopic ratios from n zircon analyses that have been perturbed by Pb loss (Z) may be modeled as the summation of n non-perturbed ${\mathrm{Pb}}^{*}/\mathrm{U}$ ratios (X) and the amount of Pb loss encountered by each (Y). X is drawn from f(t) that reflects the Gaussian distribution of ${\mathrm{Pb}}^{*}/\mathrm{U}$ ratios that are unperturbed by Pb loss, and Y is drawn from g(t) that represents the distribution of Pb loss in the sample. The distribution that characterizes Z may be found by convolving f(t) and g(t). Although we assume that f(t) is a Gaussian distribution, the distribution type of Pb loss, g(t), shown in this example as a logit–normal distribution ($\mathit{\mu }=-\mathrm{4.5}$; σ=1.0), could take a number of discrete or continuous forms (Fig. 3). Note that in our modeling framework, values of X, Y, and Z are normalized as percentage deviation from the true isotopic ratio (i.e., the mean of f(t)), where negative values indicate that measured ${\mathrm{Pb}}^{*}/\mathrm{U}$ is lower than the true ratio. See Video S1 (Sharman and Malkowski, 2023) for an animation that illustrates the process of convolution and Video S2 (Sharman and Malkowski, 2023) for an exploration of the logit–normal distribution in μ and σ parameter space.

Equation (2) may be solved analytically for some forms of f(t) and g(t). For example, the convolution of Gaussian and exponential distributions is known as the exponentially modified Gaussian distribution (Grushka, 1972). However, (f∗g)(t) may also be solved numerically, which has the advantage of allowing both f(t) and g(t) to take any form.

3 Methods

3.1 Modeling approach

We use the mathematical framework described above to model both the distribution of apparent Pb loss, g(t), experienced by a group of cogenetic crystals and their unperturbed U–Pb date distribution, f(t). Because Pb loss is isotopically indiscriminate, Eq. (2) may be equally applied to ${}^{\mathrm{206}}{\mathrm{Pb}}^{*}{/}^{\mathrm{238}}\mathrm{U}$ and ${}^{\mathrm{207}}{\mathrm{Pb}}^{*}{/}^{\mathrm{235}}\mathrm{U}$. However, we model ${}^{\mathrm{206}}{\mathrm{Pb}}^{*}{/}^{\mathrm{238}}\mathrm{U}$ ratios, as these have much lower analytical uncertainty for the Carboniferous and younger samples analyzed in this study.

To model g(t), we allow the μ of f(t) to vary within the 95 % confidence interval associated with an independent estimate of the crystallization age. We then estimate both g(t) and σ of f(t) by iteratively solving for the combination of parameters that minimize the misfit between the measured ${\mathrm{Pb}}^{*}/\mathrm{U}$ values and the modeled distribution (f∗g)(t), using the Python scipy.optimize.minimize() function. We define the misfit as the sum of squared residuals between the empirical cumulative distribution function (ECDF) of the measured ${\mathrm{Pb}}^{*}/\mathrm{U}$ values and the cumulative density function (CDF) of the modeled ${\mathrm{Pb}}^{*}/\mathrm{U}$ distribution.

If both non-CA and CA analyses are available from the same sample, then the distribution of CA U–Pb dates may be used to constrain the parameters of f(t). For such samples, we modify the approach described above by first finding the Gaussian distribution f(t) that most closely approximates the treated ${\mathrm{Pb}}^{*}/\mathrm{U}$ distribution. We then use this best-fitting f(t) in estimating g(t), using the minimization-of-misfit technique described above. Such datasets have the advantage of providing constraints on σ of f(t), which is otherwise treated as an unknown parameter during modeling if only non-CA U–Pb dates are available.

In order to estimate g(t) as described above, we must choose one or more reasonable parametric models that are appropriate for describing distributions of Pb loss. One possibility is that all zircon crystals in the sample experienced the same amount of Pb loss, which could shift ${\mathrm{Pb}}^{*}/\mathrm{U}$ from 0 % to −100 % of its value. Such a scenario of constant Pb loss may be modeled by a discrete form of g(t) where a single parameter specifies the percentage of Pb lost. Convolution of such a discrete form of g(t) simply produces a negative shift in the ${\mathrm{Pb}}^{*}/\mathrm{U}$ values (i.e., Fig. 3b).

Figure 3 Illustration of how normally distributed zircon ${\mathrm{Pb}}^{*}/\mathrm{U}$ values may be perturbed by discrete (a–c) or continuous (d) distributions of Pb loss. The top row represents the distribution of Pb loss in the sample expressed as a percentage of the true isotopic ratio (e.g., ${}^{\mathrm{206}}{\mathrm{Pb}}^{*}{/}^{\mathrm{238}}\mathrm{U}$ or ${}^{\mathrm{207}}{\mathrm{Pb}}^{*}{/}^{\mathrm{235}}\mathrm{U}$) at the time of Pb loss, where the height of the black bar and ball indicates the relative probability of the specified ${\mathrm{Pb}}^{*}/\mathrm{U}$ offset. Three discrete scenarios are shown: (a) no Pb loss, (b) constant Pb loss, and (c) isolated Pb loss. A logit–normal distribution is shown as an example of continuous Pb loss in panel (d). Additional examples of continuous Pb-loss distributions are shown in Fig. A1 in Sharman and Malkowski (2023). The bottom row shows both the relative (above) and cumulative (below) probabilities of the unperturbed (solid black line) and Pb-loss-perturbed (dashed line) ${\mathrm{Pb}}^{*}/\mathrm{U}$ distributions.

Another possibility is that Pb loss was experienced by only a subset of crystals (i.e., isolated Pb loss). This scenario may also be modeled by assigning g(t) to a discrete distribution with two parameters, with one that indicates the fraction of Pb lost and one that specifies the proportion of crystals that underwent Pb loss (Fig. 3c). This parameterization of g(t) will produce a bimodal pattern in U–Pb values, particularly if the degree of Pb loss is significant relative to measurement uncertainty (Fig. 3c).

Instead of modeling g(t) as a discrete distribution where Pb loss is restricted to certain values, we may also consider a continuous probability distribution, where values of Pb loss can take on any value between 0 % and 100 % (Fig. 3d). Rather than assume a priori the form(s) that g(t) might take, we considered a wide range of one- or two-parameter distributions for the purposes of exploratory modeling (Appendix A in Sharman and Malkowski, 2023). Of the distribution types considered, we identified the logit–normal distribution, also known as the logistic–normal distribution, as perhaps the most reasonable for modeling Pb loss. The logit–normal distribution has the property of having a logit (i.e., the quantile function of the logistic distribution) that is normally distributed with a geometric mean of μ and standard deviation of σ (Aitchison and Shen, 1980; Mead, 1965)

$$\begin{array}{}\text{(3)}& f\left(x\mathit{\mu }\mathit{\sigma }\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}}{\displaystyle \frac{\mathrm{1}}{x(\mathrm{1}-x)}}{e}^{-\frac{(\mathrm{logit}\left(x\right)-\mathit{\mu }{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\end{array}$$

for $\mathrm{0}<x<\mathrm{1}$. The logit–normal distribution is well-suited for modeling constrained data types (e.g., compositional data; Aitchison and Bacon-Shone, 1999; Vermeesch, 2018b), in part due to it being defined over $\mathrm{0}<x<\mathrm{1}$. We invert and scale the distribution to extend from −100 % < x < 0 % to match the sign and units of ${\mathrm{Pb}}^{*}/\mathrm{U}$ offset due to Pb loss when expressed as a percentage (Fig. 3d).

Figure 4 explores the relationship of the logit–normal distribution to its two parameters (μ and σ) (see also Video S2 in Sharman and Malkowski, 2023). The distribution has a “spiky” character when σ is a very small number (e.g., 0.001; Fig. 4a), which would be a reasonable approximation for samples that underwent an approximately constant amount of Pb loss (e.g., Fig. 3a and b). Although the logit–normal distribution cannot model 0 % or 100 % Pb loss, these values may be approximated by making μ a large negative or positive number, respectively. A sample where most zircon exhibit very little Pb loss but with fewer zircon experiencing significant Pb loss could be produced by μ = −4 and σ=1.0 (Fig. 4c). Alternatively, a sample with a peak probability of ${\mathrm{Pb}}^{*}/\mathrm{U}$ offset <0 % may be modeled using moderate values of σ (e.g., 0.25–1; Fig. 4b and c). The logit–normal distribution produces bimodal distributions, where most probability is close to 0 % and −100 % when σ values are high (e.g., ≫1; Fig. 4d).

Figure 4 Exploration of the logit–normal distribution's parameter space. Note that we have rescaled the x axis of the logit–normal distribution, such that $-\mathrm{100}<x<\mathrm{0}$.

3.2 Samples

We apply the mathematical and modeling framework presented above to estimate the distribution of apparent Pb loss in 10 igneous samples that range in age from Carboniferous to Miocene, 9 of which have been published previously (Table 1). Samples CTU, RCP, and SRF are all from upper Eocene rhyolites of the Caetano caldera system of the western United States (Watts et al., 2016). These samples were split into non-CA and CA aliquots prior to analysis via secondary ionization mass spectrometry (SIMS) (Watts et al., 2016). We used the error-weighted mean age of the CA U–Pb dates as an estimate of the true crystallization age for each sample, with weighted means approximately 0.4–0.6 Myr older than the corresponding ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ sanidine ages (Watts et al., 2016). The number of analyses per aliquot (non-CA or CA) ranges from 17 to 34 for these three samples (Table 1).

Table 1 Sample summary.

^a Sanidine ${}^{\mathrm{39}}\mathrm{Ar}{/}^{\mathrm{40}}\mathrm{Ar}$ age (Snow and Lux, 1999). ^b Error-weighted mean of chemically abraded U–Pb dates. ^c Concordia age (CA-ID-TIMS). ^d Error-weighted mean five of five zircon crystals analyzed via CA-ID-TIMS. ^e U–Pb dates older than 28 Ma are excluded from analysis. ^fU–Pb dates older than 158 Ma are excluded from analysis. ^g U–Pb dates older than 360 Ma are excluded from analysis. ^h D_α is the alpha dose (events g⁻¹). N is for the number of analyses. NA is for not available. W₁ is for the first Wasserstein distance. W₂ is for the second Wasserstein distance.

We present analysis of five samples reported in von Quadt et al. (2014), including upper Oligocene andesite to trachyandesite from Macedonia (248-2, 029-5, and 059-1), upper Cretaceous dolerite from Romania (DG 026), and middle Carboniferous granite from west Bulgaria (AvQ 244). These samples were also split into non-CA and CA aliquots prior to analysis via LA-ICP-MS. For samples other than 059-1 we use concordia ages from CA-ID-TIMS analyses of between three and six crystals for the crystallization age of each sample (von Quadt et al., 2014; Table 1). For sample 059-1, we used the weighted mean of the CA U–Pb dates. The number of analyses per sample (non-CA or CA) ranged from 17–55 for this dataset (Table 1).

Sample ELM18DVTC-10 is from a Miocene ash-flow tuff from the Panuga Formation in the western United States that has 144 U–Pb dates acquired via LA-ICP-MS (Miller et al., 2022). We use a ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ weighted mean age of 15.7 ± 0.2 Ma (2σ) from the same unit as an estimate of the crystallization age of this sample (sample 592-GV1 of Snow and Lux, 1999). Sample ELM18DVTC-10 was highlighted by Schwartz et al. (2022), who noted the youngest zircon U–Pb dates to be much younger than the accepted ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age of this unit. Miller et al. (2022) also noted the presence of these young zircon and suggested that they may be a consequence of surface contamination from units higher in the section.

Sample MM20-EC-109 is a Lower Cretaceous intermediate ash interbedded within marine carbonaceous mudstone from the Rio Mayer Formation of Argentina with 68 zircon U–Pb dates acquired via LA-ICP-MS (Table A3 in Sharman and Malkowski, 2023). Laser ablation spot locations were selected on the rim and/or core of the zircon guided by cathodoluminescence (CL) images (Fig. A3 in Sharman and Malkowski, 2023), with 59 zircon crystals analyzed in total. We use a crystallization age of 144.43 ± 0.07 Ma (2σ) derived from a weighted mean of five zircon crystals analyzed via CA-ID-TIMS at the Boise State University Isotope Geology Laboratory (Table A4 in Sharman and Malkowski, 2023). This sample exhibits U–Pb dates acquired via LA-ICP-MS that are up to ∼ 60 % younger than the CA-ID-TIMS weighted mean.

3.3 Statistical analysis

To evaluate the likelihood that the measured ${\mathrm{Pb}}^{*}/\mathrm{U}$ distribution could have been drawn from the modeled (f∗g)(t), we apply the nonparametric, one-sided Kolmogorov–Smirnov (K–S) and Kuiper statistical tests that compare the ECDF with the cumulative CDF of (f∗g)(t) (Press et al., 2007). The Kuiper statistic is relatively more sensitive in differences in the tails of the distributions versus the K–S statistic (Vermeesch, 2018a). We reject the null hypothesis that the non-CA U–Pb dates were drawn from (f∗g)(t) if the K–S or Kuiper p value is <0.05 (i.e., 95 % confidence level). We thus interpret p values > 0.05 to indicate that the non-CA U–Pb dates could have been plausibly drawn from (f∗g)(t) at a 95 % confidence level (Press et al., 2007). However, it should be noted that Saylor and Sundell (2016) found that both K–S and Kuiper p values more frequently reject the null hypothesis than expected. We thus use p values as a general guideline to model the goodness of fit.

The Wasserstein distance has been recently proposed as a metric for quantifying the dissimilarity between detrital zircon U–Pb age distributions (Lipp and Vermeesch, 2023). We consider the first and second Wasserstein distances, W₁ and W₂, to be useful approximations for the total degree of negative ${\mathrm{Pb}}^{*}/\mathrm{U}$ perturbation that a set of analyses has experienced,

$$\begin{array}{}\text{(4)}& {\displaystyle }{W}_{\mathrm{1}}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left|M^{-\mathrm{1}}-N^{-\mathrm{1}}\right|\mathrm{d}t\text{(5)}& {\displaystyle }{W}_{\mathrm{2}}=\sqrt{\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\left|M^{-\mathrm{1}}-N^{-\mathrm{1}}\right|}^{\mathrm{2}}\mathrm{d}t},\end{array}$$

where M⁻¹ and N⁻¹ are the inverses of the CDFs M and N. Because values of Pb loss are restricted to between 0 % and 100 %, both W₁ and W₂ yield maximum possible values of 100 (i.e., 100 % of analyses have −100 % ${\mathrm{Pb}}^{*}/\mathrm{U}$ offset, or the U–Pb system is completely reset). The W₁ simply equates to the area beneath the cumulative probability distribution of g(t) (e.g., Fig. 3). Because the W₂ distance involves a squaring of the distance between the quantile functions, it imparts a higher cost penalty for the part of the distribution with strongly offset values. For example, the W₁ and W₂ distances are equal for a Pb-loss function characterized by constant Pb loss (e.g., 3 % Pb loss produces W₁ and W₂ values of 3; Fig. 3b). However, the W₂ distance is often much larger than W₁ for Pb-loss distributions with a heavy tail (Fig. 3d). As such, the $W_{\mathrm{2}}/W_{\mathrm{1}}$ ratio provides an approximation of Pb-loss distribution asymmetry, with values of 1 indicating constant Pb loss and values ≫1 indicating highly asymmetric Pb loss.

4 Results

Of the four primary types of Pb-loss distributions considered (Fig. 3), the logit–normal distribution yielded the lowest average misfit with a value of 3.5, followed by the isolated Pb-loss scenario (average of 4.5) and the constant Pb-loss scenario (average of 10.5) (Table A2 in Sharman and Malkowski, 2023). The scenario of no Pb loss performed the worst of any scenario that we considered, with an average misfit of 101.3 (Table A2). Correspondingly, both K–S and Kuiper p values for the no-Pb-loss scenario are ≪0.05 for all samples, except SRF, suggesting that the untreated LA-ICP-MS or SIMS U–Pb dates are unlikely to have been drawn from an unperturbed U–Pb date distribution.

Figure 5 presents a comparison of actual versus modeled U–Pb date distributions for each sample, with the best-fitting logit–normal distribution shown (Table 1; see Fig. A1 in Sharman and Malkowski, 2023, for individual plots that show the fit for each sample and distribution type). We chose to not consider discrete distributions of g(t) for the “best” fit because we consider it unlikely that Pb loss (or other processes that cause negative age offsets) would be limited to discrete values (e.g., Fig. 3). Values of μ for g(t) ranged from −2.69 to −4.73, with corresponding values of σ spanning 0.29 to 1.91. W₁ distances ranged between 1.8 (sample SRF) and 8.1 (sample AvQ 244) and W₂ distances between 2.9 and 11.1 (Table 1; Fig. 5).

Figure 5 Modeling of apparent Pb loss in zircon U–Pb dates acquired via LA-ICP-MS or SIMS. The best-fitting logit–normal distribution of apparent Pb loss is shown (Table 1; see Fig. A1 in Sharman and Malkowski, 2023, for plots of all samples and apparent Pb-loss-distribution types modeled). Empirical cumulative distribution functions (ECDFs) are shown as solid lines, while model results are shown as dashed lines. See text for further discussion of model results.

Figure 6 Distributions of apparent Pb loss when modeled as a logit–normal distribution. Samples with σ<1 are shown as a dashed line.

A plot of the best-fitting logit–normal distributions displays two distinct behaviors of g(t) (Fig. 6). (1) Four samples with $\mathit{\mu }\lesssim -\mathrm{3}$ and σ>1 and have a g(t) maximum relative probability close to 0 %, suggesting a strongly decaying rate of offset (i.e., most zircon experienced very little Pb loss, while a few have more significant ${\mathrm{Pb}}^{*}/\mathrm{U}$ offset). These samples also displayed $W_{\mathrm{2}}/W_{\mathrm{1}}\ge \mathrm{1.6}$. (2) The remaining six samples that yielded σ<1 and generally higher μ values ($>-\mathrm{4}$) displayed a tendency for the mode of g(t) to be >0 %, representing more of a bulk shift in age (e.g., most U–Pb dates have some offset, while relatively few have very little or very much age offset). These samples produced $W_{\mathrm{2}}/W_{\mathrm{1}}\le \mathrm{1.3}$.

5 Discussion

5.1 Assumptions and limitations

The mathematical and modeling framework that we present includes several underlying assumptions and limitations that should be considered.

1. Because g(t) could represent any geological or analytical process that introduces negative age offsets, we use the phrase “apparent Pb loss” when describing our modeled estimates of g(t). For instance, matrix-related systematic errors (Allen and Campbell, 2012), addition of U–Th during weathering (Pidgeon et al., 2019), and even sample contamination from younger minerals could introduce negative age shifts exclusive of loss of radiogenic Pb. Common Pb corrections, particularly the ²⁰⁷Pb correction, may also introduce a bias in ${\mathrm{Pb}}^{*}/\mathrm{U}$ values (Andersen, 2002; Andersen et al., 2019). We recommend that these additional complexities in the U–Pb system be considered when interpreting modeled estimates of g(t) as representing distributions of Pb loss.

2. Our approach of parameterizing g(t) for the purpose of exploratory modeling has the advantage of yielding results that are interpretable, while also being suitable for the relatively low-n datasets available. However, any parametric model is likely a simplification of the true g(t), and thus we consider our modeled estimates of g(t) to be first-order approximations. Analyzing a greater range of samples with a greater number of ±  CA in situ U–Pb analyses, with ideal datasets having hundreds or even thousands of analyses per sample (e.g., Pullen et al., 2014; Sundell et al., 2021), would likely improve our ability to constrain the form(s) of g(t) and evaluate whether the logit–normal distribution or other forms of g(t) are appropriate. Such datasets would also be more amenable to nonparametric solutions of estimating g(t).

3. For g(t) to represent the true distribution of Pb loss, the process of convolution must be applied to ${\mathrm{Pb}}^{*}/\mathrm{U}$ ratios at the time of Pb loss. Because Pb^* is progressively added to the crystal over time, a greater amount of ancient Pb loss is required to achieve the same reduction in ${\mathrm{Pb}}^{*}/\mathrm{U}$ relative to recent Pb loss. This point is illustrated in Fig. 1, where a 50 % reduction in Pb^* at 125 Myr after crystallization produces a similar reduction in ${}^{\mathrm{206}}{\mathrm{Pb}}^{*}{/}^{\mathrm{238}}\mathrm{U}$ when compared to zircon of the same age that lost 25 % of its Pb^* at present day. For this reason, g(t) can be viewed as a minimum estimate in the case of ancient Pb loss. If the timing of Pb loss is known or can be estimated (e.g., Morris et al., 2015), the input ${\mathrm{Pb}}^{*}/\mathrm{U}$ ratios can be adjusted prior to analysis, such that g(t) more accurately reflects the true magnitude of Pb loss.

4. The modeling framework presented above is designed for a group of cogenetic crystals with a shared crystallization age (e.g., autocrystic zircon from the same magmatic episode; Miller et al., 2007). This requirement stems from our definition of apparent Pb loss as a relative shift or percentage deviation from the true isotopic value (Fig. 2). The assumption that all zircon are coeval is a simplification, as even autocrystic zircon crystallize over a period of time, typically at 10³–10⁴-year timescales (Miller et al., 2007; Rossignol et al., 2019). Multimodal detrital samples or igneous samples with xenocrystic or inherited zircon are not easily modeled because these samples would violate our assumption of a shared crystallization age. Failure to recognize the true heterogeneity in crystallization age in such a sample could cause an incorrect interpretation of the apparent Pb-loss distribution.

5. For datasets with paired non-CA and CA measurements, our modeling approach assumes that the relative precision of the analyses is similar. This is because the Gaussian distribution that best approximates the CA U–Pb date distribution, f(t), is convolved with the apparent Pb-loss distribution g(t) to fit the non-CA U–Pb date distribution. The Watts et al. (2016) SIMS dataset shows similar relative precision, regardless of the treatment approach (non-CA versus CA). Some samples from the von Quadt et al. (2014) LA-ICP-MS dataset exhibit slightly lower relative precision for non-CA versus CA, with sample AvQ 244 yielding the largest difference with an average relative precision of 1.1 % (1σ) for non-CA dates and 0.8 % (1σ) for CA dates. We suggest that for the purposes of modeling apparent Pb loss, paired non-CA and CA U–Pb datasets should be collected on the same instrument, using similar acquisition parameters to avoid introducing large changes in measurement precision. Alternatively, the CA U–Pb dates may be used to only constrain the μ of f(t) in the model, with σ treated as an unknown parameter (e.g., for paired non-CA LA-ICP-MS and CA-ID-TIMS datasets; Fig. 5a and i).

6. For datasets with paired non-CA and CA measurements, we do not consider any imperfections of the chemical abrasion process. For example, although the CA treatment aims to completely remove all radiation damaged zones of the crystal (Mattinson, 2005), it is possible to have remaining residual zones of Pb loss following treatment (e.g., Schoene et al., 2010). Any such remaining compromised domains of the crystal will yield at least some apparent Pb loss when analyzed. For instance, Watts et al. (2016) interpreted three zircon U–Pb analyses from SRF to have some residual Pb loss that was not fully accounted for by the CA process (Fig. 5g). Incorporation of Pb-loss-perturbed U–Pb dates when modeling f(t) would likely produce an underestimate of the true magnitude of the apparent Pb loss.

5.2 Distributions of apparent Pb loss

What distribution type(s) characterize apparent Pb loss in natural samples? Our results strongly suggest that at least 9 of the 10 samples modeled have at least some systematic negative offset in ${}^{\mathrm{206}}{\mathrm{Pb}}^{*}{/}^{\mathrm{238}}\mathrm{U}$ that cannot be explained by random measurement uncertainties alone. This is because the K–S and Kuiper statistical tests are unable to reject the null hypothesis for many of the apparent Pb-loss-distribution types considered (Table A1 in Sharman and Malkowski, 2023). For example, only the no-Pb-loss scenario produced a p value < 0.05 for sample MM20-EC-109, suggesting that any of the other modeled distributions of apparent Pb loss may be statistically plausible for this sample. These results suggest that we cannot confidently distinguish between discrete (constant or isolated) or continuous distributions of apparent Pb loss in the datasets modeled. Except for ELM18DVTC-10, which has 144 non-CA LA-ICP-MS analyses, the samples we analyzed have relatively low numbers of analyses (between 17 and 68; average of 32) for a given sample and treatment category (non-CA or CA) (Table 1). We suspect that the collection of larger n datasets would allow a better resolution of which parameterizations of g(t) might be most appropriate.

Even if the specific distribution type(s) that characterizes g(t) cannot be uniquely identified, our analysis suggests two contrasting behaviors in apparent Pb loss (Fig. 6). We speculate that U–Pb dates that undergo a bulk shift (i.e., $W_{\mathrm{2}}/W_{\mathrm{1}}\cong \mathrm{1}$) may reflect a population of zircon crystals with relatively homogenous characteristics (size, U content, etc.) that have all experienced a similar post-crystallization history. Correspondingly, the population of zircon that produces U–Pb dates with a highly asymmetric distribution of age offset (i.e., $W_{\mathrm{2}}/W_{\mathrm{1}}\gtrsim \mathrm{1.5}$) may reflect heterogeneity between crystals, with variable characteristics. For example, Pb loss is thought to be promoted in small zircon crystals and in zircon with elevated U (Ashwal et al., 1999; Gehrels et al., 2020), and thus distributions of particle size and/or trace element geochemistry may influence asymmetric patterns in g(t). Collection of size measurements and trace element concentrations from zircon, in addition to measurement of the U–Pb date (e.g., Watts et al., 2016), would likely help evaluate hypotheses about the underlying factors that influence apparent Pb-loss distributions. Furthermore, given the relatively small number of samples modeled in this study, we suggest that there is a need for more samples to undergo paired non-CA and CA characterization to improve our understanding of the range of behaviors that may be typical. For example, it is presently unclear whether it is more common for samples to have their U–Pb dates bulk shifted (e.g., samples 029-5, 059-1, CTU, and DG 026) versus having relatively few U–Pb dates highly offset (e.g., samples MM20-EC-109 and ELM18DVTC-10; Fig. 5).

Why do some samples experience more overall apparent Pb loss than others? Although we anticipated that apparent Pb loss would be greater for samples with greater radiation damage due to U and Th decay, our analysis shows no clear trend by alpha dose (Table 1). However, we acknowledge that the relatively high degree of apparent Pb loss modeled in the youngest sample, ELM18DVTC-10, may be a consequence of contamination from overlying units instead of true Pb loss (Miller et al., 2022). Even the three samples from the same Eocene caldera system (CTU, RCP, and SRF) showed contrasting amounts of apparent Pb loss (W₂ ranges from 2.9 to 4.2; Table 1), as noted by Watts et al. (2016). Characterizing the overall magnitude of apparent Pb loss in a wider range of samples would likely help elucidate predictive factors (if any).

5.3 Importance of quantifying the distribution of apparent Pb loss in in situ U–Pb geochronology

The overwhelming majority of published in situ U–Pb dates from zircon, minimally > 600 000 and likely in the millions of analyses (Puetz et al., 2021), have not been treated using CA. In contrast, CA is now practiced routinely in the ID-TIMS community, which has contributed to growing precision and accuracy over the past 2 decades (Schoene, 2013). However, the strategy of mitigating Pb loss through avoidance is perhaps less easily adopted to routine in situ U–Pb geochronology. For instance, there may be practical limitations with chemically abrading large numbers of zircon crystals, including the potential loss of certain age modes that would be detrimental to provenance analysis. We thus suggest that there is a pressing need to improve quantitative characterization of apparent Pb-loss distributions in non-CA in situ U–Pb datasets to aid in interpreting these datasets and to guide strategies for future data collection.

It is somewhat concerning that 9 of the 10 samples analyzed in this study exhibited statistically significant amounts of negative age offset from the estimated true crystallization age. Even a small age offset of a few percent, or cryptic Pb loss (Kryza et al., 2012; Watts et al., 2016), has potentially important repercussions for interpreting the age and rates of geologic events and processes. For example, there is a growing awareness in the detrital geochronological community that the youngest zircon U–Pb dates often skew unexpectedly young relative to the plausible crystallization age (e.g., Herriot et al., 2019; Gehrels et al., 2020; Schwartz et al., 2022). Presently, there is no consensus on the importance of post-depositional Pb loss on influencing depositional age interpretations (e.g., Herriott et al., 2019; Copeland, 2020; Schwartz et al., 2022). Sample MM20-EC-109 illustrates the risk well; we initially interpreted the young tail on the U–Pb date distribution to suggest a depositional age of ∼ 125 Ma, based on the youngest cluster of overlapping U–Pb dates. The youngest single analysis was a 60.5 ± 2.4 Ma rim on a 135.3 ± 3.0 Ma core, with the second-youngest analysis being a 79 ± 1.2 Ma date measured from the core of a zircon crystal, with the corresponding rim yielding an older 129.8 ± 3.6 Ma date (Table A2 in Sharman and Malkowski, 2023). Interpretation of the youngest single U–Pb date or dates as the depositional age of this sample would have produced a highly erroneous estimate, off by up to −58 % of the true eruption age of 144.50 ± 0.07 (2σ) Ma, as determined by CA-ID-TIMS. Because this ash is interbedded within a sequence of organic-rich marine mudstone in the Austral Basin of Argentina, the misinterpretation in this case could have led to an erroneous depositional age model, with implications for interpreting the paleoclimatic and geodynamic context of these sediments.

Although modeling detrital samples was outside of the scope of this study, we believe that our results bear upon maximum depositional age analysis. The tendency for the youngest U–Pb dates in a sample to be affected by Pb loss (or other similar process) complicates even conservative estimates of the maximum depositional age (Dickinson and Gehrels, 2009; Coutts et al., 2019; Schwartz et al., 2022). If apparent Pb loss follows a continuous distribution (e.g., Fig. 3d), then it is ill-advised to assume that outlying U–Pb dates may be rejected, while the rest are considered unperturbed (see also discussion in Copeland, 2020). Even an interpretation based on the peak age probability of the youngest age mode is likely to be too young because the process of convolution produces a young shift in the mode of the distribution, in addition to creating a young tail (Fig. 3d and Fig. A1 in Sharman and Malkowski, 2023). Because existing methods of calculating the maximum depositional age (Dickinson and Gehrels, 2009; Coutts et al., 2019; Vermeesch, 2021) do not account for systematic negative age offsets, our analysis suggests that there is a higher probability for erroneous estimates of the maximum depositional age if (1) there are a large number of zircon crystals with crystallization ages that are close to the age of deposition, (2) the overall number of measured U–Pb analyses is very high, and/or (3) the magnitude of apparent Pb loss is high. In addition, a heavy-tailed distribution of apparent Pb loss (i.e., $W_{\mathrm{2}}/W_{\mathrm{1}}\gg \mathrm{1}$) will result in a greater probability of finding extremely offset ${\mathrm{Pb}}^{*}/\mathrm{U}$ values.

6 Conclusions

This study presents a mathematical framework for quantifying the distribution of apparent Pb loss on U–Pb date distributions, which could include true loss of radiogenic Pb or other processes that also produce a systematically negative age offset. We show that a Pb-loss-perturbed U–Pb date distribution from a set of zircon crystals with a shared crystallization age can be represented by the convolution of a Gaussian distribution that reflects measurement uncertainty in ${\mathrm{Pb}}^{*}/\mathrm{U}$, with a distribution that characterizes Pb loss, g(t). Our approach relies on analyzing differences between the untreated ${\mathrm{Pb}}^{*}/\mathrm{U}$ distribution from in situ U–Pb geochronology (i.e., LA-ICP-MS or SIMS) and an independent estimate of the true crystallization age, which could include U–Pb dates from a thermally annealed and chemically abraded aliquot of the same sample or from another geochronometer (e.g., ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$). We suggest that the first and second Wasserstein distances (W₁ and W₂) of the apparent Pb-loss distribution can be used to quantify the total degree of apparent Pb loss that a set of zircon analyses has undergone, with maximum possible W₁ and W₂ values of 100.

We apply this modeling framework to 10 igneous samples (Miocene to Carboniferous) analyzed with LA-ICP-MS or SIMS. All but one of the samples showed a high probability that the untreated U–Pb date distribution has been perturbed by Pb loss or other equivalent process. Although our analysis shows that multiple parameterizations of g(t) can achieve statistically acceptable fits (i.e., K–S or Kuiper p value > 0.05), we suggest that the logit–normal distribution may be a reasonable choice for exploratory modeling of apparent Pb-loss distributions. However, we caution that the number of analyses in the samples we analyzed was generally low (17–144; average of 39); future efforts to characterize g(t) may be promoted by collection of larger n datasets and through development of nonparametric methods of estimating g(t). Furthermore, our estimates of g(t) should be viewed as minimum estimates of the true amount of Pb lost, as we assumed present-day Pb loss in our analysis. These caveats aside, we noted two behaviors of apparent Pb loss: samples with a bulk shift in U–Pb date distributions ($W_{\mathrm{2}}/W_{\mathrm{1}}\lesssim \mathrm{1.3}$) and samples where most analyses had very little offset but fewer had much larger offsets ($W_{\mathrm{2}}/W_{\mathrm{1}}\gtrsim \mathrm{1.6}$). The overall magnitude of ${\mathrm{Pb}}^{*}/\mathrm{U}$ decrease was also found to be variable, with median values varying from −0.9 % to −6.4 %.

Given the widespread application of in situ U–Pb geochronology of untreated zircon across many disciplines of geosciences, improved characterization of both the distribution type(s) and magnitude of apparent Pb loss is warranted, particularly for Phanerozoic zircon, where cryptic Pb loss is difficult to identify. We highlight a need for increased sampling and high-n characterization of paired non-CA and CA in situ U–Pb datasets. In addition, we recommend simultaneous collection of parameters such as zircon size and trace elemental concentrations to aid in future efforts to understand the mechanisms of negative age offsets. Ultimately, we anticipate that improved characterization of the magnitude of apparent Pb loss will aid in interpreting non-CA in situ U–Pb datasets and guide strategies for future data collection.