1 Introduction

Detrital geochronology is often the only way to estimate the depositional age of siliclastic rocks in the absence of fossils or volcanic ash layers. Detrital zircon U–Pb geochronology in particular has become a popular technique to obtain maximum depositional ages (MDAs). Numerous MDA estimation algorithms have been proposed over the years (Nelson, 2001; Barbeau et al., 2009; Dickinson and Gehrels, 2009; Tucker et al., 2013; Chen et al., 2016; Zhang et al., 2016; Ross et al., 2017; Herriott et al., 2019; Copeland, 2020). In a detailed review paper that was recently published in this journal, Coutts et al. (2019) compared and contrasted the most popular among these methods, namely:

1. the youngest single grain (YSG);
2. the mode of the youngest graphical peak on a probability density plot (YPP);
3. the youngest grain cluster at 1σ (YGC1σ) or 2σ (YGC2σ);
4. the youngest detrital zircon estimated by Ludwig (2003)’s Monte Carlo resampling algorithm (YDZ);
5. the outcome of Ludwig and Mundil (2002)’s TuffZirc algorithm;
6. the weighted mean of the youngest three (Y3Z) or four (Y4Z) zircons;
7. the weighted mean of the grains in the youngest peak of a probability density plot (τ); and
8. the minimum age of grains selected by Gehrels (2003)’s AgePick algorithm.

Additionally, they also introduced a new estimator:

1. the ‘Youngest Statistical Population’ (YSP), which groups the youngest sub-sample of more than two grains that pass a Chi-square test for homogeneity.

Using numerical simulations and synthetic age distributions, Coutts et al. (2019) found that all but three of these methods gradually drift to younger ages with increasing sample size. The only exceptions were the YPP, τ and YSP methods, which converged to ages that were systematically too old. The failure of YPP and τ to retrieve the correct depositional age reflects the underlying flaws of the probability density plots (PDPs) on which they are based. The shortcomings of PDPs are explained in detail by Vermeesch (2012) and Vermeesch (2018a) and won’t be discussed further in this paper. The positive bias of the YSP method will be briefly discussed in Section 7. The remainder of the paper will focus on the undesirable drift of the remaining six MDA estimation algorithms towards geologically unreasonable young values.

Section 2 will prove that it is theoretically impossible for the YSG, YGC1σ/2σ, YDZ and Y3Z/Y4Z methods to converge to the correct solution even in a best case scenario. Section 3 shows that the ‘minimum age model’ of Galbraith and Laslett (1993) does not suffer from this problem. This maximum likelihood based method was developed for detrital fission track thermochronology but can be equally useful for zircon U–Pb studies in a modified form developed by Galbraith (2005, p.107). We will refer to the minimum age model as the ‘maximum likelihood age’ (MLA) method to avoid confusion with minimum depositional age estimation, which is an entirely different subject.

Section 4 applies the MLA algorithm to a detrital zircon U–Pb dataset from the Colorado Plateau Coring Project. This case study validates the different MDA methods by comparing measurements obtained by Laser Ablation Inductively Coupled Plasma Mass Spectrometry (LA-ICPMS) with measurements obtained on the same samples by Chemical Abrasion Thermal Ionisation Mass Spectrometry (CA-TIMS). The new MLA algorithm is the only method that passes this stringent test.

Despite the advantages of the MLA method over all existing MDA estimation algorithms, it is not a silver bullet for samples that have a fat tail at the young end of the age spectrum, or that have been affected by post-depositional Pb-loss. Section 6 discusses these limitations and provides suggestions to mitigate their effects.

2 On the failure of existing MDA algorithms to converge

The YSG, YGC1σ/2σ, YDZ and Y3Z/Y4Z methods essentially amount to averaging the n youngest grains obtained by some selection criterion. The selection criterion can be simple (e.g., YSG) or complex (YDZ). But all these algorithms boil down to analysing the young tails of detrital age spectra. This is a difficult thing to do with ad hoc methods.

Detrital age spectra are the convolution of two probability distributions: (1) the distribution of the true (but unknown) zircon U–Pb ages; and (2) the distribution of their analytical errors. Analytical uncertainty obscures the true ages and must be accounted for during MDA estimation. Geochronologists generally make the reasonable assumption that the analytical errors follow a normal distribution with mean μ = 0 and standard deviation σ. Under this condition it is well known that the actual U–Pb age is approximately 95% likely to fall in a ±2σ interval around the measured value. However this is only true when considering a single analysis.

When multiple measurements are considered together, then the likelihood that the measured value is more than ±2σ away from the true value increases. In statistics this is known as a ‘Type-1 error’. As the number of analyses increases, so does the probability that at least one value falls outside the ±2σ interval. Given a sufficiently large sample size, there inevitably comes a point when the data contain values that differ from the true value by 3σ or more. This is true for both tails of the normal distribution. In the context of MDA estimation, it means that the youngest date in a huge sample may be less than the actual depositional age.

To further explore this important point, let us consider the simplest case of an MDA estimation exercise. Suppose that 100% of the grains in a sample are concordant and syn-depositional at 10 Ma, with normally distributed uncertainties of 1 Ma at 1σ. The solution to this toy example is trivial. The depositional age is simply given by the mean of all the dates. Yet six of the nine MDA estimation algorithms reviewed by Coutts et al. (2019) fail to retrieve it (YPP, τ and YSP don’t, but will fail in nearly all non-trivial examples).

Figure 1 shows why existing MDA estimation methods fail the toy example. It plots the expected age¹ of the n^th youngest grain among a sample of N grains drawn from our normal distribution with mean μ = 10 and standard deviation σ = 1. For example, consider a sample containing N = 5 grains, then the expected age of the youngest grain (n = 1) is the 1/5^th quantile of the normal distribution, which is 9.16 Ma. When the sample size is increased to 20 grains, then the expected value for the youngest age is given by the 1/20^th quantile, which is 8.36 Ma.

Thus, the youngest single grain age monotonically decreases with increasing sample size, and so does the second youngest grain (n = 2), the third youngest grain (n = 3) and so forth. Also the average of the youngest n grains decreases with increasing sample size, and all the other grain selection criteria do so as well. The failure of existing MDA estimation algorithms to solve the simplest toy example undermines their credibility in more complex scenarios.

3 The maximum likelihood model of Galbraith and Laslett (1993)

The failure of existing MDA estimation algorithms to find the correct solution for the toy example is diagnostic of the fact that they are ad hoc algorithms that lack a formal statistical basis. It is desirable for statistical estimation techniques to converge to the true solution with increasing sample size. Unfortunately none of the existing method manage to do so. The conventional way to solve this issue is to parameterise the problem and determine the parameters by Maximum Likelihood Estimation (MLE). Geological examples of this approach include isochron regression (Titterington and Halliday, 1979; York et al., 2004), concordia age estimation (Ludwig, 1998) and finite mixture modelling (Galbraith and Green, 1990; Sambridge and Compston, 1994). The MLE approach can also be used for MDA estimation. In fact, such an algorithm already exists.

Galbraith and Laslett (1993) introduced a minimum age estimator for fission track thermochronology that is equally applicable to U–Pb geochronology. The model asumes that the data were drawn from a two component mixture, in which a fraction π_γ of the population is derived from a discrete minimum age peak at t_m = exp[γ] and the remaining fraction (1 -π_γ) of the grains were drawn from a (log)normal distribution with location parameter μ and dispersion parameter σ, truncated at t_m (Figure 2). The four model parameters (γ, π_γ, μ and σ) can be estimated by maximising the likelihood function using an appropriate error model for the analytical uncertainties.

The original fission track implementation of Galbraith and Laslett (1993)’s maximum likelihood model assumed binomial counting errors (van der Touw et al., 1997). An alternative formulation using normal errors was formulated by Galbraith (2005, p.107). This version of the model is most appropriate for MDA estimation by detrital U–Pb geochronology. The minimum age model will be referred to as the MLA (Maximum Likelihood Age) method in the remainder of this paper, so as to avoid any confusion with minimum depositional age estimation.

Standard MLE theory also provides a way to estimate the uncertainties of the model parameters. This is done by inverting the negative matrix of second derivatives (i.e. the Hessian) of the log-likelihood function with respect to the parameters. Experience shows that the uncertainties of γ and π_γ are generally smaller than those of μ and σ. Because μ and σ are not directly relevant to the MDA estimation problem anyway, there is little harm in reducing the number of model parameters from four to three, by requiring that γ = μ. The resulting gain in numerical stability benefits small datasets whilst having only a minor effect on the accuracy of the MDA estimates.

It can be shown that among all statistical estimators, the MLE approach is the most efficient, which means that it yields the most precise estimates for any given sample size. MLE algorithms are also consistent, which means that they are guaranteed to converge to the true solution with increasing data size (assuming that the model assumptions are met). Thus the MLA method fixes the key issue with existing MDA estimation methods. Applying it to the toy example asymptotically yields parameter values of π_γ = 1, σ = 0 and γ = 2.3, which corresponds to the correct answer of t_m = 10 Ma.

4 Application to U–Pb data from the Colorado Plateau

Leaving the toy example behind and moving on to a geologically more realistic example, we will now apply the MLA model to a recently published dataset of detrital zircon U–Pb ages obtained from the Colorado Plateau Coring Project (CPCP). The dataset contains 13 detrital zircon samples that were extracted from a ~520 m drill core in Petrified Forest National Park (Arizona, USA).

The 13 samples belong to the Permo-Triassic Coconino, Moenkopi and Chinle Formations. Between 221 and 308 randomly selected zircon grains from each sample were dated by LA-ICPMS, yielding 13 U–Pb age spectra ranging from 189 to 3428 Ma (Gehrels et al., 2020). Subsequently, the youngest 2 – 19 grains from each sample were extracted from the LA-ICPMS grain mount and re-analysed by high precision CA-TIMS (Rasmussen et al., 2020). This paired LA-ICPMS + CA-TIMS dataset allows us to compare the performance of the different MDA estimators across a range of sample sizes and analytical precision.

Figure 3 shows the LA-ICPMS and CA-TIMS age estimates for sample 297 on so-called radial plots. These are a graphical device that was invented by Rex Galbraith, the creator of the MLA model (Galbraith, 1988, 1990). The radial plot is designed to visualise heteroscedastic datasets (i.e., datasets with unequal uncertainties). It is the most elegant way to visualise the results of the MLA model, which explicitly takes into account this heterescedasticity. See the Appendix for further details.

The LA-ICPMS data of sample 297 follow a bimodal age distribution. This is evident in Figure 3 as two linear arrays on the radial plot radiating from the origin towards the radial scale at ~220 Ma and ~1800 Ma, respectively. Unsurprisingly, the radial plot of the CA-TIMS data looks very different. It exhibits a unimodal distribution with ages ranging from 220 to 226Ma.

From the radial plots it is evident that age dispersion exceeds analytical uncertainty for both the LA-ICPMS and CA-TIMS datasets. In fact none of the 13 samples pass a Chi-square test for age homogeneity, indicating the presence of geological dispersion. The overdispersion of the LA-ICPMS data is expected because these were meant to capture multiple provenance components. However the excess scatter of the CA-TIMS data merits further discussion.

The overdispersion of the CA-TIMS data could either reflect the presence of multiple detrital components of similar age, or it could be caused by protracted pre-eruptive residence of zircon in a volcanic magma chamber. In either case, the CA-TIMS data present a similar statistical challenge as the LA-ICPMS data in that it is not immediately obvious how to estimate the MDA.

To meet this challenge, all the MDA estimate algorithms were applied to both the LA-ICPMS and CA-TIMS datasets. Thus we can objectively compare the effects of the different algorithms and of the different analytical techniques (Section 5).

5 Empirical comparison of the different MDA algorithms

Both radial plots and MLA calculations were made in Vermeesch (IsoplotR, 2018b). See Appendix B for further details.

Applying the MLA method to sample 297 yields MDA estimates of 219.68±0.46 Ma and 221.09±0.99 Ma for the LA-ICPMS and CA-TIMS data, respectively. The small (-0.6%) age difference between the two estimates contrasts starkly with existing methods such as YSG (-17.4%), YGC1σ (-5.4%), YGC2σ (-4.5%), Y3Z (-7.8%), Y4Z (-5.2%), YDZ (-17.5%) and YSP (-2.1%). Repeating the same exercise for the remaining 12 samples confirms this picture.

Figure 4 compares the MDA estimates of the LA-ICPMS and CA-TIMS data for all the methods discussed in this paper. The results are shown as bivariate scatter plots with 95% confidence intervals shown as error bars. These plots can be divided into three zones.

1. The 1:1 line marks samples whose LA-ICPMS and CA-TIMS based MDA estimates agree within error. This is the expected scenario for compositionally homogenous zircon crystals without common Pb.
2. The area below the 1:1 line groups LA-ICPMS based MDA estimates that are older than their CA-TIMS based counterparts. This may indicate complications such as common Pb or the presence of small scale growth zones.
3. The area above the 1:1 line is a ‘forbidden zone’. In general, one would not expect CA-TIMS ages to exceed the corresponding LA-ICPMS ages, unless the zircons have experienced post-depositional Pb-loss. It is unlikely that such Pb-loss is a common occurence in nature, for reasons that will be discussed in Section 6.

All but two of the ten ad hoc MDA estimation methods shown on Figure 4 spill over into the ‘forbidden zone’. The YSG, YGC1σ, YGC2σ, Y3Z, Y4Z and YDZ methods are the worst offenders because all their LA-ICPMS based MDA estimates are younger than their CA-TIMS based MDA estimates. This troubling result is likely caused by the sample-size effect discussed in Section 2: because the LA-ICPMS datasets are two orders of magnitude larger than the CA-TIMS datasets, they are much more likely to drift towards unreasonably young values.

The MLA method does not suffer from this problem. All its MDA estimates either fall on or below the 1:1 line. The only two ad hoc methods that fare reasonably well are YPP and TuffZirc. However the YPP method does not provide a measure of uncertainty, while the TuffZircon dates include one sample in the forbidden zone and one sample whose LA-ICPMS age is three times older than the CA-TIMS age. All things considered, the MLA method provides the most sensible results.

The first two panels of Figure 4 show the results of the 3-parameter and 4-parameter versions of the MLA algorithm, respectively. The difference between the two models is negligible (≪1%), reflecting the insensitivity of the method to the parametric model assumptions.

6 Limitations

Any statistical model is only as good as the assumptions that it makes. For the minimum age model, these assumptions are that

1. the true age distribution approximately follows the functional form shown in Figure 2; and
2. the analytical uncertainties are well characterised.

The first assumption is generally easy to verify. If the young end of the age spectrum is marked by a cluster of nearly concordant ages, then the minimum age model will appropriately average these. And it will do so even if the old end of the spectrum does not look like a (log)normal distribution. It is not uncommon for detrital zircon U–Pb age distributions to be fat tailed at both ends of the spectrum. However, the similar results obtained by the 3- and 4-parameter formulations of the MLA algorithm indicates that this is generally not a problem.

If the age difference between the youngest and the second youngest grains in a sample is significantly greater than their respective analytical uncertainties, then the minimum age model will simply return the youngest age as a result. In other words, for fat tailed age spectra, the minimum age model reduces to the YSG model. This is the most sensible solution from a statistical point of view. In the absence of a discrete youngest age peak, the objections to the YSG model raised in Section 2 and Figure 1 do not apply. However whether age of the youngest grain is also the most sensible solution from a geological point of view is a different matter.

Post-depositional Pb-loss is one way to violate assumption 2. Zircon grains that have experienced Pb-loss can yield precise yet inaccurate U–Pb age estimates that plot in the ‘forbidden zone’ of the LA-ICPMS vs. CA-TIMS plot (Section 5). The minimum age model is unable to detect this problem. But neither, of course, are any of the existing MDA estimation algorithms. In the absence of a statistically sound solution, it is useful to inspect any textural or compositional data for clues to identify grains that have been affected by Pb-loss. For example, if the young outliers are all characterised by extremely high U- and/or Th-concentrations, then this could be taken as evidence for radiation damage induced Pb-loss. Such additional information can be visualised on the radial plot as optional fill colours.

However, given the low diffusivity of Pb in zircon at the temperatures found in most sedimentary basins, post-depositional Pb-loss is probably unlikely to be a common problem in detrital geochronology (Cherniak and Watson, 2001; Copeland, 2020). It is possible that many anomalously young ages that have been attributed to common Pb-loss to are instead just ordinary outliers violating assumption 2. The probability that a dataset contains such anomalously young values increases with sample size. In principle these values could be removed by outlier detection algorithms. But all such algorithms are heuristic by nature and are not guaranteed to produce sensible results.

A better albeit more onerous solution for fat tailed age distributions is to re-analyse the youngest grains either by internal isochron analysis (Nemchin and Cawood, 2005), or by CA-TIMS geochronology as shown in this paper. After such validation, even a single isolated grain of zircon can provide a robust MDA constraint.

7 Discussion

This paper has shown that the MLA model of Galbraith and Laslett (1993) and Galbraith (2005) is superior to all the MDA estimation algorithms reviewed by Coutts et al. (2019). It is the only method that is built on solid statistical foundations, and the only method that can converge to the correct solution with increasing sample size. In the absence of complications with recent Pb-loss or young outliers, the MLA method is entirely objective and completely hands off. It does not require any arbitrary decision such as the number of grains to average (Y3Z, Y4Z), or the probability cutoff to use (YGC1σ, YGC2σ).

Although this objectivity is appealing in several ways, it lacks any geological context. Copeland (2020) argues that detrital zircon grains cannot be grouped into discrete (maximum depositional) age components without additional geochemical justification. However whilst it is true that MLA estimates may not correspond to a particular geological event, it should be emphasised that they are not meant to do so. The principal purpose of any MDA algorithm is to obtain an upper limit for the depositional age. Whether this upper limit corresponds to a particular geological event, or to a mixture of multiple overlapping events, is a different issue.

At first glance, the MLA model may seem similar to Coutts et al. (2019)’s YSP method, because both algorithms assume that the minimum age belongs to a discrete age component. However there is a crucial difference between the two techniques. The YSP method assumes that it is possible to make a clean separation between grains whose true ages belong to the youngest age peak, and grains whose true ages belong to older age components. In reality this distinction is blurred by the analytical uncertainties. The rank order of single grain age estimates may not necessarily be the same as that of the true U–Pb ages. This makes the YSP method biased relative to the minimum age model. The latter does not arbitrarily assign the grains to two groups, but jointly considers all the age estimates when fitting the model parameters.

Besides its applications in fission track thermochronology and detrital geochronology, the MLA model can be used in other Earth Science applications as well, such as the determination of eruption ages from collections of volcanic zircon ages. Volcanic rocks often exhibit positively skewed age distributions with a short tail of syn-eruptive U–Pb ages and a long tail of pre-eruptive or xenocrystic ages. In this context the MLA model would be a good alternative to ad hoc approaches such as YSG and YSP that have previously been used, or to the more sophisticated Bayesian approaches that have recently been proposed (Keller et al., 2018). It would also be useful to visualise such volcanic datasets on radial plots instead of the interval plots that are currently used.

8 Acknowledgments

P.V. would like to thank George Gehrels, Peter Copeland and Daniel Coutts for constructive reviews that significantly improved the paper. This research was supported by NERC standard grant #NE/T001518/1 (‘Beyond Isoplot’).

Appendix A: about radial plots

Like the MLA model itself, radial plots were originally developed by Rex Galbraith (1988, 2005) for the purpose fission track geochronology. And like the MLA model, radial plots can be equally useful for U–Pb and other geochronometers (Vermeesch, 2009, 2018b). Let z_j and s_j be the log-transformed values of N single grain ages t_j and their standard errors σ_j (for 1 ≤ j ≤ N):

(1)

then a radial plot is a bivariate (x_j,y_j) scatterplot with:

(2)

where z_∘ is a reference value such as the weighted mean of the z_j values. Radial plots allow the observer to simultaneously assess both the magnitude and the precision of quantitative data. Each point on the diagram represents a single grain analysis. Old grains plot at high (positive) angles to the origin, whereas young grains plot at low (negative) angles. The analytical uncertainty can be obtained by extrapolating lines from the origin to the radial scale through the top and the bottom of an imaginary 2σ-error bar added to each sample point. Thus, precise measurements plot towards the right hand side of the diagram, whereas imprecise measurements plot towards the left. Drawing two parallel lines at 2σ distances from either side of the origin allow the analyst to visually assess whether all the single grain ages within a sample agree within the analytical uncertainties.

With its ability to visualise heteroscedastic datasets (i.e., datasets with unequal uncertainties), the radial plot realises the unfulfilled promise of the probability density plot (Vermeesch, 2012). It is the most elegant way to visualise the results of the MLA model, which explicitly takes into account this heterescedasticity. Discrete components (including despositional age components) in a mixed age distribution can be easily spotted as linear arrays radiating from the origin towards the radial scale (left panel of Figure 3). Radial plots and the MLA model have been implemented in Java (RadialPlotter, Vermeesch, 2009) and R (IsoplotR, Vermeesch, 2018b). Detailed instructions for the latter application are provided in Appendix B.

Appendix B: Implementation in IsoplotR

IsoplotR is an R package for radiometric geochronology (Vermeesch, 2018b). Figure 5 shows its online user interface. The software implements two ways to plot radial plots and calculate minimum ages. The first of these is to (1) select other as a data type, (2) paste the U–Pb ages and their uncertainties into the input window, (3) choose radial plot as an output device, (4) select Options → Finite mixtures → minimum, and (5) click PLOT. A second way to use the minimum age model is to (1) select U-Pb as the data type, (2) paste the U–Pb data into the input window, (3) proceed as before. Both methods produce identical results.

Alternatively, the same calculations can also be performed from the command line with the R programming language. Using the ages and their uncertainties:

or using the U–Pb compositions:

where 297ages.csv and 297isotopes.csv are two comma separated text files with the age data and isotopic compositions, respectively. These files are provided in the Supplementary Information. 297ages.csv specifies the uncertainties as absolute values (standard errors), whereas 297isotopes.csv contains relative uncertainties (coefficients of variation). This difference is handled by the ierr argument to the read.data function.

References