the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
The daughter–parent plot: a tool for analyzing thermochronological data
Birk Härtel
Eva Enkelmann
Data plots of daughter against parent concentration (D–P plots) are a potential tool for analyzing lowtemperature thermochronology, similar to isochron plots in radioisotopic geochronology. Their purposes are to visualize the main term of the radiometric age equation – the daughter–parent ratio – and to inspect the daughter–parent relationship for anomalies indicating influences of geological processes or analytical bias. The main advantages of the D–P plot over other data analysis tools are (1) its ability to detect systematic offsets in D and P concentrations, (2) its unambiguous representation of radiationdamagedependent daughter retention, and (3) the possibility to analyze potential age outliers.
Despite these benefits, the D–P plot is currently not used for analyzing lowtemperature thermochronology data, e.g., from fissiontrack, (U–Th) $/$ He, or zircon Raman dating. We present a simple, decisiontreebased classification for daughter–parent relationships based on the D–P plot that places a dataset into one of seven classes: linear relationship with zero intercept, cluster, linear relationship with systematic offset, nonlinear relationship, several age populations, scattered data, and inverse relationship. Assigning a class to a dataset enables choosing further data analysis steps and how to report a sample age, e.g., as a pooled, central, or isochron age or a range of ages. This classification scheme aims at facilitating thermochronological data analysis and making decisions more transparent. We demonstrate the proposed procedure by analyzing published datasets from a variety of geological settings and thermochronometers and introduce Incaplot, which is graphical user interface software that we developed to facilitate D–P plotting of thermochronology data.
 Article
(4605 KB)  Fulltext XML

Supplement
(318 KB)  BibTeX
 EndNote
The isochron plot is a universal tool for analyzing geochronological results, e.g., U–Pb, Ar–Ar, or Rb–Sr data (e.g., Nicolaysen, 1961). The main reason for its use is that the ratio of the isotope ratios (e.g., ^{87}Sr $/$ ^{86}Sr vs. ^{87}Rb $/$ ^{86}Sr) on the plot's axes is the essential term of the radiometric age equation. The slope and intercept of an isochron fitted to a dataset convey information about the age and initial isotopic composition of a sample. Furthermore, the isochron plot enables us to visualize anomalous features in the data, such as outliers or excess of radiogenic daughters.
The isochron plot's equivalent for lowtemperature thermochronology is the radiogenic daughter (D) vs. radioactive parent (P) plot (D–P plot), which several authors suggest for analyzing fissiontrack (FT), (U–Th) $/$ He (He), and zircon Raman (ZR) data (e.g., Fanale and Kulp, 1962; Green, 1981; Wernicke and Lippolt, 1993; Dunkl, 2002; Vermeesch, 2008; He et al., 2021; Härtel et al., 2022a). This plot allows us to (1) detect systematic offsets in daughter or parent concentration (e.g., Vermeesch, 2008); (2) analyze the influence of radiation damage on daughter retention, avoiding spurious associations in the age–eU (effective uranium concentration) plot (Härtel et al., 2022a); and (3) evaluate singlegrain ages in terms of a twodimensional distribution (e.g., for detecting outliers) or selecting a sample age (e.g., as a mean, pooled, central, or isochron age) without a preconceived idea about the singlegrain ages.
The D–P plot thus occupies the interface between the analytical results and more specific data analysis tools such as radial, kernel density estimate (KDE), or age–grain size plots – a tool that helps us to decide which data analysis techniques are applicable or not to a given dataset. It is therefore surprising that the D–P plot is not considered a standard tool for analyzing thermochronological data (e.g., Flowers et al., 2022; Kohn et al., 2024). Our aim is to fill this gap and provide guidance to users of lowtemperature thermochronology.
This article consists of two major parts: we first provide the theoretical background of the D–P plot, its differences from the classic isochron plot, and examples of commonly observed daughter–parent relationships (Sect. 2). Then, we illustrate how to analyze data in D–P space using a workflow for classifying daughter–parent relationships, suggest further data analysis tools for each type of relationship, and, if applicable, provide algorithms for sample age calculation (Sect. 3). We also introduce Incaplot, which is free, graphical user interface software dedicated to creating D–P plots that allows an easy implementation of our proposed analysis to FT, He, and ZR data.
2.1 Deriving the D–P plot
Using a plot of daughter (D) against parent (P) concentration rests upon the general age equation
where t is the age, λ is the decay constant, and c is a constant to balance out the units of D and P. It is evident from Eq. (1) that the age has a onetoone relationship with D $/$ P. Therefore, the position of a data point in a plot of D vs. P indicates the singlegrain age by the slope of a line connecting it to the origin of the plot (Fig. 1a). Data pairs from sameage grains plot on a line through the origin, representing a proportional relationship (Fig. 1b). The D–P plot is thus a graphical representation of the age equation.
This relationship of the D–P plot to the age equation is the same as that of the classic isochron plot, but there are two significant differences: (1) the isochron plot represents parent and daughter concentrations as isotope ratios with a nonradiogenic sister isotope as the common denominator. This creates error correlation between the two axes of the plot, which is not present in the D–P plot as it relies on independently measured daughters and parents. (2) The isochron plot assumes the initial presence of the radiogenic daughter isotope, which makes isochron fitting indispensable for age calculation. In contrast, for the D–P plot no initial daughters are assumed, enabling the analyst to examine the D–P relationship for patterns without the need for an isochron. To honor these differences, we prefer the generic term D–P plot over isochron plot for this type of diagram for FT, He, and ZR data.
The actual quantities of D and P depend on the dating method. For FT dating, the daughters represent the number or areal density of spontaneous tracks and the parents represent either that of induced tracks (externaldetector method) or U concentration (LAICPMSbased dating). The daughters for He and ZR dating are the αejectioncorrected He concentration and the radiation damage density, respectively. However, defining a parent concentration for these methods is difficult because several αemitting nuclides – ^{238}U, ^{235}U, ^{232}Th (and ^{147}Sm) – have to be considered. One solution is to express the parents as an effective uranium concentration (eU) – the sum of the parent concentrations weighted by their relative αproduction rate. This reduces the number of parents to one. Appendix A discusses the calculation of eU as a parent concentration in (U–Th) $/$ He and zircon Raman dating and the differences between existing eU equations (e.g., Cooperdock et al., 2019; Härtel et al., 2023). Appendix B provides additional discussion on the choice of daughter and parent concentration units for different dating methods.
2.2 Data patterns for multigrain samples
In practice, the analyst acquires multiple singlegrain data to extract information about a sample's thermal history. The number of these singlegrain analyses varies between methods and analytical protocol – from about 20–30 grains per sample for FT and ZR dating to only 3–5 grains per sample for wholegrain He dating. The D–P plot allows us to analyze such multigrain samples. However, real data deviate from the ideal trend in Fig. 1b. In the following, we give a short overview of the typical deviations from the ideal proportional D–P relationship and the information they contain regarding geological processes that influence rock cooling and heating or analytical biases.
2.2.1 Linear relationship with zero intercept
Figure 2a presents a synthetic example of a positive linear D–P relationship with a zero intercept, including random variation about the trend. This is similar to the proportional case, but with uncertainty in the D and P measurements. Additional variation may be the consequence of varying grain sizes or inaccurate αejection correction for He dating, intergrain chemical differences for FT or ZR dating, and parent concentration zoning for all three methods. The D–P plot in Fig. 2b shows an example of a linear relationship with a zero intercept for laserablation apatite He data from the Fish Canyon Tuff reference material (Pickering et al., 2020). Both the isochron (26 ± 7 Ma) and the pooled age (28.3 ± 0.6 Ma) overlap with the reference age at 28.4 Ma (Schmitz and Bowring, 2001).
2.2.2 Cluster
Figure 3a shows a synthetic example of clustered D–P data. This pattern is typical for data from samples with limited intergrain differences in parent (and daughter) concentrations, and usually their uncertainty intervals overlap strongly. In this case, the positive relationship between daughters and parents may be obscured by the uncertainty. Figure 3b shows an example of a D–P plot with clustered apatite FT data from sample FC1 from the Duluth Complex, Minnesota (Härtel et al., 2022a). Despite relatively large differences in track density, the uncertainties in D and P of most grains overlap. The data give a pooled age of 850 ± 30 Ma.
2.2.3 Linear relationship with systematic offset
In Fig. 4a, the synthetic data form a linear trend, which, compared to Figs. 2a and b, is offset from the origin. In He dating, such an offset may result from (1) “parentless helium” implanted by inclusions (Vermeesch et al., 2007), eUbearing grain boundary or neighboring phases (e.g., Spiegel et al., 2009; Murray et al., 2014), or (2) a consistent style of zoning across grains affecting αejection correction (e.g., Orme et al., 2015). In FT dating, it may also be due to a bias towards higher or lower track counts (see Green, 1981). In ZR dating, systematic offsets may result from damage calibration issues, asymmetric Raman bands, or compositionrelated Ramanband broadening (Kempe et al., 2018; Troch et al., 2018; Härtel et al., 2022b). Note that an overestimation or underestimation of P causes an apparent offset of opposite sign in D.
The D–P plot in Fig. 4b shows an example of negative offset in wholegrain zircon He data from a set of four closely spaced samples of Miocene leucogranite from the Greater Himalaya sequence (Orme et al., 2015). The singlegrain ages range from 9.9–14.7 Ma (weighted means: 10–12 Ma), whereas Orme et al. (2015) expected an age range of 14–17 Ma due to hostrock stratigraphy and other thermochronological data. They explained this by the zircon grains consistently showing compositional zoning with loweU cores and higheU rims: this causes more He to be lost by α ejection than accounted for by conventional Ft correction (e.g., Hourigan et al., 2005) and leads to the negative offset. They tested this assumption by adjusting Ft of some grains using zoning information from laserablation depth drilling (black circles in Fig. 4b). The ages range from 14.8 to 17.0 Ma (weighted mean: 15.6 ± 0.2 Ma). In the D–P plot, these data points fall above the main trend and show insignificant offset from the origin. The isochron ages for both unadjusted (14 ± 1 Ma) and adjusted data (16 ± 2 Ma) overlap with each other and with the expected age range.
2.2.4 Nonlinear relationship
Figure 5a showcases a synthetic example of a nonlinear D–P relationship. This may be due to the daughter retention depending on the degree of lattice damage from α decay of U, Th, and their daughters. The production of radiation damage is roughly proportional to the parent (eU) concentration. Its effect on daughter retention causes D and P to form either a concave (Fig. 5a, damageenhanced loss) or a convex (damageenhanced retention) relationship (Härtel et al., 2022a).
Figure 5b shows an example for a nonlinear D–P relationship due to radiationdamageenhanced helium loss in zircon He data from the Minnesota River Valley (Miltich, 2005). The dataset consists of several samples assumed to have shared the same thermal history since ∼ 1.8 Ga based on earlier thermochronological data (see references in Miltich, 2005). The He concentration increases approximately linearly with eU increasing up to 500 µg g^{−1} and falls at higher eU concentrations in response to radiation damage facilitating He loss from the zircon crystals. Guenthner et al. (2013) suggested a thermal history for these samples based on the zircon radiation damage accumulation and annealing model (ZRDAAM, black line), consistent with the D–P relationship.
2.2.5 Several populations
The synthetic example in Fig. 6a shows data forming two linear trends in the D–P plot, indicating different age components within the sample. Such a trend occurs if a sample contains groups of grains with a high contrast in kinetic properties. Figure 6b shows a D–P plot for an example of different age populations found in apatite FT data for a fully reset sedimentary sample from the Mackenzie Basin, Northwest Territory (sample I77; Issler et al., 2005). It displays two roughly linear trends in the data corresponding to two different ages (90 ± 12 and 220 ± 45 Ma). Colorcoding the data by the chlorine content shows a slight compositional difference between the two age populations, suggesting a chemical influence on FT annealing properties (e.g., Barbarand et al., 2003).
2.2.6 Scattered data
The D–P plot of synthetic data in Fig. 7a shows how random scatter can obscure the relationship of D and P. Such a pattern may result from different factors, e.g., heterogeneous daughter retention within the sample, e.g., a broad range of grain sizes or chemical compositions. Other reasons for scattered data might be the occurrence of microcracks, deformation, or parent zoning. In addition, scatter may arise from analytical factors, such as variably biased αejection correction, counting bias, or a combination of these factors.
Figure 7b shows an example of scattered data in the D–P plot from the multigrainaliquot apatite He data from the Bighorn Mountains, Wyoming (Shell Canyon 12; Reiners and Farley, 2001). The data show no relationship between He and eU. However, colorcoding the different aliquots by the massweighted average radius (MWAR) reveals an increasing age (i.e., D $/$ P ratio) with increasing grain size. This indicates a continuous age distribution due to different sensitivity of differently sized grains to volume diffusion of helium. Figure 7c shows another example of a scattered D–P relationship in wholegrain zircon He data from the Laxemar region on the Fennoscandian Shield (Guenthner et al., 2017). The color coding reflects the sampling depth, and the black line represents the ZRDAAM from the original publication. However, neither the depth of each sample – a proxy for their current temperature – nor the radiation damage model explains the scatter in the data. In this case, an unknown factor causes the age variation.
2.2.7 Inverse relationship
Figure 8a shows a synthetic dataset with an inverse relationship in the D–P plot. This pattern may occur due to (1) a small sample size causing a spurious relationship (Ketcham et al., 2018), (2) bias from over or undercorrecting the He concentration for α ejection, or (3) the data representing a falling segment of a nonlinear trend caused by radiation damage. Figure 8b provides an example for a negative D–P trend from wholegrain zircon He data from Proterozoic samples from Baffin Island (sample A1042, Ault et al., 2018; Armstrong et al., 2024). Ault et al. (2018) interpreted the age variation in this dataset to be due to radiationdamageenhanced He loss, as the ZRDAAM (black line) in Fig. 8b shows. Armstrong et al. (2024) provided additional Raman data on selected grains, showing that some of the zircon grains with eU ≥ 1000 µg g^{−1} were metamict, pointing to enhanced He loss compared to the lowereU grains (Guenthner et al., 2013).
2.3 Unique benefits of D–P plots
Figures 2–8 show the range of D–P patterns that occur in thermochronological data. While the mean D $/$ P ratio of the synthetic datasets shown in these figures is the same (2), each of these relationships requires different considerations for data analysis. This includes the question of whether reporting a single sample age is appropriate and, if so, which type of sample age to report. Commonly used data analysis tools such radial plots, KDE, or age–grain size plots help to trace some of the factors causing age variations, but there are unique benefits to analyzing data in the D–P plot.
First, the D–P plot is the only thermochronological data plot that enables us to detect systematic offset in daughter or parent concentrations (Fig. 4a, b). Systematically offset data pose a serious problem to many standard data analysis tools and should therefore be treated with caution: (1) singlegrain ages calculated from offset data are biased towards higher or lower ages depending on the sign of the offset. This bias propagates into calculated central tendencies (Härtel et al., 2022a) and into plots displaying the age as a variable, such as radial, KDE, age–grain size, and age–(e)U plots. (2) Offset data appear overdispersed (and fail the χ^{2} test) because the data uncertainties do not explain the spread in age. This further complicates the use of radial plots, as the spread in singlegrain ages may give way to a misinterpretation of ages as a mixture of discrete age components (see discussion in Vermeesch, 2019). (3) The overdispersion by systematic offset hampers inverse thermal history modeling, as the modeling algorithm will have to reconcile a large spread in ages without the uncertainties accounting for it (e.g., Vermeesch and Tian, 2014). The offset affects each data point differently so that expanding the uncertainties in D and P does not solve this problem (see Flowers et al., 2022). (4) Systematic offset also compromises the helioplot (Vermeesch, 2010), which determines the age from log ratios, because it disturbs all ratios derived from the D and P concentrations.
Figure 9 illustrates some of these biases in more traditional data analysis tools using the zircon He data of Orme et al. (2015, Fig. 4b). Figure 9a reproduces the D–P plot of the data in Fig. 4b showing the negative, systematic offset, which is the result of a biased αejection correction due to consistent compositional zoning (Orme et al., 2014). However, the radial plot in Fig. 9b does not show any anomaly in the data except for overdispersion (P(χ^{2})≈0; dispersion = 10 ± 3 %). The singlegrain ages (9.9–14.7 Ma) and the central age (11.4 ± 0.5 Ma) are substantially younger than the age range of 14–17 Ma established from stratigraphic and thermochronological constraints (purple area; Orme et al., 2015). In contrast, the isochron age (14±1 Ma) from the D–P plot fits this scenario well. In the age–eU plot (Fig. 9c), the data show a weak, positive association. However, this association is misleading in comparison to the usual interpretation of associations in the age–eU plot (see below) because the age variation arises from biased αejection correction and not a radiation damage effect on He retention (Guenthner et al., 2013; Gautheron et al., 2020).
Second, the D–P plot provides an unbiased indication of whether daughter retention in a sample depends on radiation damage: the D–P plot shows unambiguous nonlinear or inverse relationships for welldocumented cases of radiationdamagedependent daughter retention (e.g., Figs. 5b, 8b). The commonly used age–eU plot also shows a relationship for such cases (e.g., Miltich, 2005; Guenthner et al., 2013; Ault et al., 2018; Armstrong et al., 2024). However, not all associations observed in the age–eU plot reflect actual radiation damage effects but may be the result of spurious ratio correlation (e.g., Carter, 1990; Härtel et al., 2022a). Figure 10 illustrates this problem using the examples from Figs. 7 and 8. Figure 10a and d show the D–P and age–eU plots for the zircon He data from Fig. 8b with an inverse daughter–parent relationship due to radiationdamagedependent He retention (color coding and black line in Fig. 10a; Ault et al., 2018; Armstrong et al., 2024). The age–eU plot shows a negative association as expected for zircon with a high radiation damage density (Guenthner et al., 2013). For comparison, Fig. 10b and e show the D–P and age–eU plot for the scattered zircon He data from Fig. 7c (Guenthner et al., 2017), for which a radiation damage model does not explain the variation in daughters and parents (black line in Fig. 10b). The age–eU plot in Fig. 10e shows a similar negative association with Fig. 10d. Also, the apparent relationship between age and eU concentration masks the scatter clearly visible in Fig. 10b. Figure 10c and f show the D–P and age–eU plot for the scattered apatite He data from Fig. 7b (Reiners and Farley, 2001), whose variation is explainable by grain size differences. However, the age–eU plot in Fig. 10f also shows a negative association. In this case, the effect of grain size on He diffusion causes daughter and parent concentrations to vary, with this variation translating to a spurious association when reprojecting the data in age–eU space. The D–P plot therefore reliably detects radiation damage effects, whereas the age–eU plot often displays falsepositive age–eU associations. Härtel et al. (2022a) provide more examples and discussion on spurious age–eU associations.
Third, the D–P plot allows us to detect age outliers in twodimensional space, not only from singlegrain ages (e.g., He et al., 2021). It thus allows us to identify the relative position of outliers with respect to the rest of the data, showing whether its main deviation occurs in D or P. In addition, it is advantageous to identify outliers without directly considering the singlegrain ages, as this may bias the detection in favor of grains that fit an a priori assumption of the sample age well.
These advantages support a unique perspective of the D–P plot, which allows us to sidestep biases that other data analysis tools show towards data with certain D–P relationships, most notably systematically offset data. We therefore suggest the D–P plot as a first step for thermochronological data analysis to identify the D–P relationship. Based on this relationship, it is possible to choose unbiased, more specific data analysis tools such as radial, KDE, or age–grain size plots or thermal history modeling. The following section presents a practical approach to using the D–P plot for data analysis.
Our data analysis scheme rests on a decision tree approach to classify the daughter–parent relationship (Fig. 11). Depending on the class of the relationship, we then suggest further steps of data analysis. The following sections outline the use of the decision tree to systematically classify the data and find an appropriate description of the contained thermal history information.
3.1 Preliminary considerations
Before using the classification scheme in Fig. 11, it is essential to ensure that the analytical procedures and samples meet certain quality criteria established for each method, e.g., that suitable grains were selected for He dating, that data with asymmetric Raman bands were excluded from ZR dating, and that track counting was conducted on prismatic grain surfaces. Also, the number of analyses in the dataset is important, as fitting a regression line or splitting a dataset into age populations is not appropriate for small datasets (see Sect. 3.5). Another criterion to be considered is the geological background of the sample. For example, a crystalline bedrock sample with a simple cooling history will likely give a single age, while a metasedimentary rock may show different age populations due to chemical variation between grains, and a volcanic rock recording its eruption is expected to give a nearideal linear trend.
Radiation damage effects and accompanying nonlinear relationships are expected for old rocks with protracted or complex cooling histories, but not for young rocks that did not spend time in the temperature regime of radiation damage accumulation. The interpretation of a sample that strongly deviates from the geological expectations needs to be carried out with care.
3.2 The classification procedure
For analyzing the data, we calculate the daughter and parent concentrations according to the thermochronological method used (see Appendix B) and plot daughters against parents. The analysis proceeds by following the decision tree in Fig. 11 to classify the daughter–parent relationship. The first step separates datasets showing a positive D–P relationship from those that do not (A in Fig. 11). We expect a positive association between D and P from the radioactive production equation, but this association may be obscured by factors discussed in Sect. 2.2. In the case of data for which the D–P relationship is not clear, it is usually safe to assume that there is no positive relationship – a decision that may be revised in later steps. Data for which D and P are not positively associated are then classified as either clustered, scattered, or following an inverse relationship (B in Fig. 11).
For data with a positive D–P relationship, it is then essential to distinguish datasets containing a singleage population from those with several populations (C in Fig. 11). As in Fig. 6, multipleage populations form linear arrays with different slopes or clusters in the dataset with gaps between them. A KDE plot may reveal the presence of different populations for cases that are not clearcut.
For singlepopulation data, the next step is filtering the dataset for outliers (D in Fig. 11). Outliers stick out by a difference in singlegrain age from the other data beyond their uncertainty. However, this is not sufficient evidence to mark a data pair as anomalous: other factors such as systematic offset may also cause single grains to be significantly older or younger than the others (Figs. 4b, 9b). In the D–P plot, outliers show up as removed from the main trend or group of data points. Before considering such a measurement to be anomalous, other properties should be examined, e.g., grain size or mineral chemistry. If anomalous data are excluded from further analysis, this should be reported, e.g., by marking the excluded data point as empty symbol in the D–P plot. For ambiguous cases, it may be advantageous to carry out the further steps with and without the concerned data point. For He dating, Flowers et al. (2022) provide further strategies for treating outliers (their Sect. 3.1).
After examining the outliers, we test the data for a linear D–P association (E in Fig. 11). If it is not clear whether the data show a linear or a nonlinear trend from visual inspection alone, this can be verified by fitting a regression line to the data and examining the residuals, i.e., the deviation of the data points from the line. For a linear relationship, the residuals scatter randomly around zero, while in the case of a nonlinear relationship, there is an association between residual and parent concentration. Figure 12a shows the linear fit to a synthetic dataset in a D–P plot. Figure 12b plots the fitting residuals against P, revealing a boomerangshaped trend that points to a nonlinear D–P relationship.
If the D–P relationship is linear, it is necessary to test the data for a systematic offset (F in Fig. 11). This is achieved by fitting a regression line to the data and examining its intercept. If the intercept includes zero in its uncertainty envelope, the offset is not significant and the data may be treated as having a zero intercept. If the uncertainty envelope does not include zero, this is a sign for a potential offset. However, this uncertainty in the intercept may be an underestimate if the variation of the data strongly exceeds that expected from the uncertainties (e.g., high MSWD; Wendt and Carl, 1991; see Appendix C). Another simple test for an intercept is the comparison of the isochron age and the pooled age: if the data form a trend through the origin, the two ages should be indistinguishable because the pooled age assumes a zero intercept (see Sect. 3.4.1).
3.3 Sample age calculation
Upon arriving at a certain class of D–P relationships, the goal is to assign an age to the sample. This can either be a central tendency, such as a mean or pooled age, an isochron age for a sample with a singleage population, or a number of ages or a range of singlegrain age depending on the D–P relationship. If the given ages can be described by a single sample age, the simplest solution is to report a central tendency. Despite its simplicity, the (arithmetic) mean age does usually not provide a reliable sample age (e.g., Vermeesch, 2008; Härtel et al., 2022a; see Appendix C). A more robust alternative is the pooled age, which uses the ratio of the summed D and P concentrations (see Figs. 2b, 3b).
If the intrasample age variation is related to a certain grain property affecting radiogenic daughter retention, the ages may represent a continuous mixture of ages, with each grain recording a different age due to its individual properties. Figure 7b shows an example with He data varying with respect to grain size. Such a mixture is best described by the central age (e.g., Galbraith, 2005; Vermeesch, 2019) or by thermal history modeling taking into account the specific grain property.
Datasets that are systematically offset require a different approach, that of the isochron age, which rests on the slope of a fitted regression line through the D–P data (see Figs. 2b, 4b). If several discrete age components exist in a dataset, these can be separated by mixture modeling (e.g., Galbraith and Laslett, 1993; Vermeesch, 2019) or by treating each age component as a single sample (see Fig. 6). If the data cannot be described by a single age, multiple ages, or a continuous mixture related to grain properties, it is still possible to report the range of singlegrain ages, which does not rely on any model assumptions. Appendix C provides a more detailed discussion about mean and isochron ages, as well as discrete and continuous age mixtures.
Table 1 shows an example for reporting format using the data from Figs. 2–8. To make the process of data analysis transparent, we recommend either showing the daughter–parent plot for each sample or at least reporting the class of the D–P relationship and the type of the reported age. Table 1 shows an example for reporting format using the data from Figs. 2–8.
Note: LA – laser ablation, EDM – externaldetector method, WG – wholegrain method, MG – multigrainaliquot method, AHe – apatite He, AFT – apatite FT, ZHe – zircon He. The age uncertainties are 2 SD.
The following sections provide specific suggestions for how to treat data falling into each of the D–P classes of Fig. 11.
3.4 Classes of daughter–parent relationship
3.4.1 Linear relationship with zero intercept
If the daughter–parent relationship is linear and the intercept of its regression line is close to zero (F in Fig. 11), the pooled and the isochron age are similar (Fig. 2b). In this case, it is advantageous to report the pooled age, which is more robust and does not require the intercept as an additional parameter. As all singlegrain ages along the linear trend are roughly the same, the potential bias of the pooled age is negligible (see Appendix C). If the MSWD or spine factor of the fitted regression line (F in Fig. 11) is outside the upper confidence limit, the data are overdispersed. This points to two possible scenarios. (1) The first is analytical dispersion due to the uncertainties not reflecting the actual measurement error. This is especially a problem for He and laserablation FT dating (e.g., Fitzgerald et al., 2006; Ketcham et al., 2018; Cogné and Gallagher, 2021). In this case, the uncertainty in the pooled age may be expanded to account for the variation of the individual analyses (see Eq. C6 in Appendix C). (2) The second is geological dispersion due to heterogeneous grain properties affecting daughter retention, such as grain size and composition. This can be tested by plotting the age against these properties or by using them for colorcoding the D–P plot (Fig. 7b). If the data are dispersed due to a continuous range of grain properties, the central age describes the age distribution best (Appendix C). In this case, thermal history modeling may take into account the variation of this grain property (e.g., grain size or radiation damage).
3.4.2 Cluster
Clustered data are best summarized by the pooled age. It is advantageous to display the pooled age as a line in the D–P plot (Fig. 2b) to check if the singlegrain D and P uncertainties overlap with the pooled age. To make sure that there is no bias towards the oldest or highestD–P grains, the data should be screened for outliers (D in Fig. 11). If the data are overdispersed, e.g., failing the χ^{2} test (e.g., Galbraith, 2005), the uncertainty of the pooled age may be expanded to reflect the actual intergrain age variation (see Eq. C6 in Appendix C) or the data may be treated as scattered (Sect. 3.4.6). If a relationship exists between age and grain properties, e.g., by plotting the age against these grain properties or colorcoding the D–P plot (e.g., Fig. 7b), the age distribution may be described by a central age or, if possible, a thermal history model accounting for the effect of the specific grain property on age.
3.4.3 Linear relationship with systematic offset
Systematically offset data must be treated with caution as such data pose problems for many common data analysis tools (see Sect. 2.3). The only sample age that may appropriately describe systematically offset data is the isochron age determined from the slope of a regression line (see Fig. 4b; Sect. 2.3; Appendix C2). Another option is to verify the reason for the intercept, such as zoning, parentless helium, or a counting bias, and to develop an analytical strategy to eliminate the bias (e.g., Spiegel et al., 2009; Orme et al., 2015). The intercept of the regression line provides a firstorder estimate for the amount of offset. If the intercept is large, close to the mean daughter concentration, or if the data allow for a horizontal or vertical line fit, they could also be treated as a cluster (Sect. 3.4.2) or as scattered data (Sect. 3.4.6). If the data are overdispersed, e.g., showing an MSWD outside the confidence interval, it is possible to expand the uncertainty in the isochron age by multiplying it by $\sqrt{\left(\text{MSWD}\right)}$ (e.g., Ludwig, 2012). For a strong overdispersion (e.g., MSWD > 10), the data should be treated as scattered (see Sect. 3.4.6).
3.4.4 Nonlinear relationship
A nonlinear relationship in the D–P plot points to radiationdamagedependent daughter retention. This assumption can be tested against independent radiation damage measurements. Raman and infrared spectroscopy, or Xray diffraction, provides radiation damage estimates for zircon or titanite (e.g., Nasdala et al., 1995; Deliens et al., 1977; Holland and Gottfried, 1955; Heller et al., 2019), while optical absorption and Raman spectroscopy are potential tools to measure radiation damage in apatite (e.g., Ritter and Märk, 1984; Liu et al., 2008). Alternatively, a nonlinear D–P relationship could result from daughter retention depending on other grain properties and the different grains recording the same thermal history differently. This effect can be examined by plotting the age against these parameters or by colorcoding the D–P plot (e.g., Fig. 7b). If such a relationship exists, the dataset may be described by a central age (see Appendix C). If the decision for a nonlinear versus a linear relationship with an offset is not clear (E in Figs. 11; 12), the less complex linear model should be preferred over a nonlinear model (Sect. 3.4.3) in the absence of independent radiation damage measurements. For a nonlinear trend caused by radiationdamagedependent daughter retention, forward modeling of daughter retention and radiation damage accumulation as well as annealing provides further insights into a sample's thermal history (e.g., Flowers et al., 2009; Willett et al., 2017; Guenthner et al., 2013). In this case, the D–P plot allows us to compare the data to the D–P relationship predicted by the model, especially in the loweU region, where the model prediction connects to the origin (Härtel et al., 2022a). Figures 5b, 7c, and 8b show thermal history forward models for zircon He dating plotted as lines in the D–P plot.
3.4.5 Several populations
If the D–P plot suggests that several discrete age components are present in the sample, the KDE or radial plot is the standard tool to examine the data. The occurrence of different components should also be tested for consistency, e.g., if a mixture of populations makes sense in the geological context (Sect. 3.1) or by colorcoding according to a variable that may underlie the different populations (see Fig. 6b). The age distribution can either be described by a finitemixture model (e.g., Galbraith and Green, 1990; Galbraith and Laslett, 1993; Galbraith, 2005; Vermeesch, 2019) or by separating the data into age populations to be analyzed individually according to the procedure in Fig. 11.
3.4.6 Scattered data
Data that vary strongly in age and are scattered in the D–P plot may result from several scenarios: first, they may be a consequence of underestimating the uncertainties with respect to the variation in the singlegrain data (e.g., for He dating, Fitzgerald et al., 2006; Brown et al., 2013). Martin et al. (2023) and Zeigler et al. (2023) showed that the uncertainty related to αejection correction in wholegrain He dating is especially difficult to estimate, while the correction contributes significantly to the age error. Data with limited scatter, for which the uncertainties may be underestimated, may be treated as a cluster (Sect. 3.4.2). A second explanation for scatter is the occurrence of different age populations, which can be verified in a KDE plot (Sect. 3.4.5). Third, the scatter may also be due to each grain having slightly different daughter retention properties and recording a different age. Plotting the age against these parameters or colorcoding the D–P plot (Fig. 7b) allows us to assess this relationship; a central age may be used to describe such a continuous mixture (Sect. 3.4.2; Appendix C).
If the scatter cannot be explained by one of these scenarios (e.g., Fig. 7c), the range of the singlegrain ages should be reported. Scattered data also pose a serious problem to inverse time–temperature (t–T) modeling, as the age difference may not allow for a single t–T path to reconcile the spread in ages.
3.4.7 Inverse relationship
An inverse daughter–parent relationship runs contrary to the relationship expected from the age equation (Fig. 8). In general, two scenarios can account for this relationship without pointing to an analytical problem. If the dataset is small (e.g., n≤5), a spurious inverse trend could arise randomly (Ketcham et al., 2018) and the dataset should be treated as scattered (Sect. 3.4.6). However, the interpretation of small datasets should be carried out with caution (see Sect. 4). Alternatively, the inverse relationship may represent an inverse segment of a nonlinear trend if radiation damage controls daughter retention (Sect. 3.4.4; Fig. 8b). If there is no clear explanation for the inverse daughter–parent relationship, it is best to report the range of singlegrain ages (Appendix C).
3.5 Practical limits of D–P plotting
The data analysis workflow in Fig. 11 provides simple decision paths and criteria for assigning a dataset to a class. This has the advantage of keeping the data analysis process consistent, especially for studies involving many samples. Still, this decisionbased approach has some limits that we would like to point out.
First, our ability to evaluate the D–P relationship for a sample clearly depends on the number of data. There are several limits a small sample imposes on data analysis using the workflow in Fig. 11: (1) it is not possible to recognize different populations; (2) a single outlier may constitute a large proportion of the gathered data; (3) random variation may cause inverse D–P relationships (see Ketcham et al., 2018) or spurious associations between the age and other properties; and (4) in terms of sample ages, the small number of grains inhibits the use of isochron or central ages, which would require the fitting of several parameters (age and intercept or dispersion) to a small amount of data. While this hampers a strict classification following Fig. 11, it is still possible to use the D–P plot as a qualitative guide, e.g., to visualize the data in terms of their variation in D and P. It also enables examining the D–P direction in which a potential outlier deviates from the rest of the data. For example, this helps to decide if the pooled age is biased towards a single highD or P grain (see Appendix C1). In this case, we recommend checking this potential outlier or reporting the singlegrain age range. The number of analyzed grains is not a concern for FT and ZR dating (n>10), but it is a limiting factor for conventional wholegrain He dating (n<10). However, the recent development of laserablationbased He dating will increase the number of grains analyzed per sample and recognize D–P relationships (e.g., TripathyLang et al., 2013; Pickering et al., 2020). In addition, some cases may allow grouping together data from several small samples. This approach hinges on the condition that the different samples are comparable, e.g., that they share the same thermal history in the partial annealing–retention zone of the thermochronometer used. This strategy is often used for analyzing He data with respect to radiation damage effects (e.g., Figs. 7b, 8b; Guenthner et al., 2017; Baughman et al., 2017; Ault et al., 2018; Armstrong et al., 2024).
Second, not all datasets may be unambiguously assignable to a class. Examples may be cases of moderate variation falling between clustered or scattered data or cases in which the distinction between linear and nonlinear relationships is not clear (see Fig. 12). While Sect. 3.4 provides suggestions for alternative classifications, this problem highlights the necessity of transparent reporting on the decisions made by the analyst (Table 1).
Third, unreset or partially reset detrital samples often record a complex mixture of pre and postdepositional thermal history. They also often contain grains with different chemical composition and size. Detrital samples are therefore not expected to fit into the simple categories of Fig. 11. Extracting a sample age or an interpretation from a single sample or a single thermochronometer is usually not possible (e.g., Carter, 2019). Standard procedures for interpreting detrital thermochronological data include identifying peak ages in the singlegrain age distribution and putting them into the context of the stratigraphic age, age distributions of source areas, and catchment geometry (e.g., Malusà and Fitzgerald, 2019). While it is possible to evaluate different age populations in the D–P plot (see Sect. 3.4.5), KDE or radial plots are the more adequate tools for this task. Still, the D–P plot may hold additional information that is difficult to access with these plots. First, it may be used on a subset of the data to evaluate the daughter–parent relationship for a given age population and possibly detect a nonlinear or systematically offset relationship (Sect. 3.4.3, 3.4.4). However, this can only be done reliably if enough data (e.g., n≥10) are available in this grain population. Second, it may help to identify bias in grain selection. One of these is the problem with overlapping, uncountable fission tracks in old or Urich zircon that may skew zircon fissiontrack (ZFT) age populations towards younger ages and thus affect the interpretation in terms of sourcearea exhumation and erosion patterns (e.g., Malusà, 2019).
Figure 13a shows the D–P plot for a synthetic ZFT dataset. The dashed line marks the countability limit for the spontaneous tracks. This limit cuts off the track density distribution for an old grain population; it indicates that the sample may contain older or higherU grains not datable with the ZFT method. Another application is the visualization of different grain populations with respect to age, parent concentration, and other grain properties, e.g., grain size or composition, to highlight nuances in the composition of different age populations. Figure 13b shows the D–P plot for a synthetic apatite fissiontrack (AFT) dataset, colorcoded by the Cl $/$ F ratio and with a dashed line representing the depositional age. In this case, some of the grains in the age group slightly older than the depositional age stand out due to high induced track density (high U content) and Fdominated halogen composition. So, despite the complexity of detrital samples, there are situations in which the visualization of the data in a D–P plot can be useful.
3.6 D–P plotting in Incaplot
This section briefly describes Incaplot (Härtel, 2024), which is simple, Pythonbased graphical user interface software dedicated to producing D–P plots. Existing software (e.g., Trackkey, Isoplot Excel, IsoplotR) already provides the tools for D–P plotting, but these are often buried between other functions or are available for certain dating methods only. Incaplot is available for free at https://doi.org/10.5281/zenodo.10637446 (Härtel, 2024) as a onefile executable for Mac (MacOS 10.15 Catalina and younger) and Windows operating systems (Windows 8 and younger).
Incaplot allows creating D–P plots and calculating lowtemperature themochronometric ages. It also provides a range of visualization and customization options. Figure 14 shows Incaplot's main window (left), its data inspection tool (upper right), and its graphical output (lower right). The main window consists of three frames dedicated to (1) loading data files, (2) the input data and calculation algorithms to be used, and (3) modifying the plots and calculations.
Incaplot requires the input files to be Excel spreadsheet files in .xls or .xlsx format or commaseparated (.csv) text files. The plotting variables need to be organized as columns with the variable names in the first row. A user manual for the current Incaplot version and an example file displaying the input data format are available in Incaplot's Zenodo repository.
Incaplot provides a range of plot customization options, which include customizing markers, axes and ticks, adding line segments to plots, and colorcoding plots by discrete and continuous variables. While Incaplot was set up to handle mainly He, ZR, and FT data, it can also be used for other dating systems or generic scatterplots. The output plots are exportable in different raster (.jpg, .png, .tif) and vector formats (.svg, .pdf, .eps).
Besides D–P plotting, Incaplot contains functions for sample age calculation as pooled age, isochron fitting with different algorithms (see Appendix C2), calculation of singlegrain ages, and effective uranium concentrations (see Eq. 1 and Appendix A).
A plot of daughter vs. parent concentration (D–P plot) represents a graphical solution of the age equation in radiometric dating and is an effective tool to reveal information in lowtemperature thermochronology data. Its unique advantage is its capability to detect systematic offsets or radiation damage effects in the data, which often compromise other data analysis tools. It also enables the analyst to identify potential outliers with respect to both daughter and parent concentration rather than the singlegrain age only. These advantages make it an ideal first step for data analysis, allowing us to adapt the analysis strategy to a given data pattern. We show several published datasets exemplifying the range of possible D–P relationships and the underlying geological factors, and we propose a new workflow for using D–P plots in thermochronological data analysis. Our approach follows a stepwise examination of the daughter–parent relationship and assigns one of seven classes to it. Based on the daughter–parent relationships, it provides criteria to choose further data analysis tools and – if appropriate – calculate a sample age. The classification scheme is an attempt to make data analysis more consistent and transparent. Our classification approach has limitations, especially when applied to small or detrital datasets; however, the D–P plot itself may still provide relevant insights in these cases. We also introduce Incaplot, which is free, graphical user interface software inviting everyone to create and customize D–P plots in a straightforward way.
The effective uranium concentration (eU) is a summary of the αproducing U, Th, and Sm concentrations, rescaling them to a common decay rate of U:
with the terms in the brackets being the concentrations in units of mass, and k_{U}, k_{Th}, and k_{Sm} being coefficients for each concentration. There are currently two definitions of eU that result in slightly different coefficients. Shuster et al. (2006) and Cooperdock et al. (2019) recalculate the actinide concentrations to a concentration of total U, whereas Härtel et al. (2021) recalculate them to the decay rate of ^{238}U only. The latter approach enables us to use eU as a single parent with a welldefined decay rate for He and ZR dating. It also considers the change of the daughter production rate over geological time instead of using presenttime production rates. Härtel et al. (2023) showed that the formulation
gives accurate results for samples at $\mathrm{30}<t<\mathrm{1000}$ Ma but may be modified if the expected ages for a set of samples are consistently higher or constrained well enough to calculate them more accurately. The coefficients for eU are derived in Eqs. (A3)–(A9). The starting point is the αproduction equation:
N(α) is the number of alpha decays; N_{A} is the Avogadro constant; 8, 7, 6, and 1 are the numbers of alpha particles produced by the respective decay series; M represents the molar masses; λ presents the decay constants; and the symbols in brackets are the concentrations in units of mass. The constants used in the calculations are summarized in Table A1. Rescaling all summands to the terms of ^{238}U gives the following:
This equation can be simplified by replacing the weighted actinide concentrations in the square brackets by eU:
This results in
w_{235}, w_{238}, and w_{147} are the mass fractions of the ^{235}U, ^{238}U, and ^{147}Sm isotopes, and the terms in square brackets are element concentrations. Equations (A7)–(A9) define the coefficients in Eq. (A1) for each element:
Figure A1 shows how the normalization coefficients for each αproducing element change with respect to the age of a sample.
The time dependence in Eqs. (A7)–(A9) also allows iterative age calculation for He and ZR dating. This requires calculating eU from Eq. (A2) and then alternating between calculating the age from Eq. (1) and recalculating eU from Eqs. (A1) and (A7)–(A9) until the solutions converge.
Daughter and parent concentrations can be expressed differently in externaldetectormethod FT and wholegrain He dating. Several criteria can be considered to find the right set of units for the D–P plot.
In He dating, the pairs of daughters (He) and parents (eU from U, Th, Sm) can either be expressed in units of abundance and mass (e.g., fmol and ng) or as concentrations (e.g., nmol g^{−1} and µg g^{−1}). The difference between these units is the normalization by the mass of the analyzed grain. For nonnormalized data, the size or mass of the analyzed grains will introduce variation into D and P that is unrelated to the age of the sample. In the case that the grains differ strongly in size, this may bias the pooled age towards the largest grains and the isochron age towards the smallest or the largest ones (see Appendix C). Rescaling the units of D and P to concentrations eliminates this potential bias. Furthermore, it is advantageous to correct the He concentration for αejection correction before calculating the age: correcting for α ejection after age calculation introduces a positive bias to the age (e.g., Vermeesch, 2008). Therefore, the corrected He concentration should be used as the daughter concentration for plotting. In externaldetector FT dating, a similar question of units arises concerning the use of either the spontaneous and induced track counts or their track densities. In this case, it is advantageous to use the track densities instead of the counts to avoid bias towards big grains.
The specific units then determine the value of the constant c in Eq. (1). Rearranging it to a daughter production equation gives
For ZR dating, c results from equating Eqs. (B1) and (A8):
Given input damage densities in 10^{16}α g^{−1} and eU concentrations in µg g^{−1}, c takes a value of 0.494 [10^{−16} µg α^{−1}].
For He dating, the same relationship as for ZR dating applies, with the difference of He concentrations usually being reported in molar concentrations:
If the input He concentrations are in nmol g^{−1} and the eU concentrations in µg g^{−1}, c takes a value of 0.02976 [µg nmol^{−1}].
For FT dating, the constant c depends on measured experimental factors. This gives
for the externaldetector method, where 0.5 is the geometry factor, λ_{D} is the total decay constant for ^{238}U, ζ is the proportionality factor determined from dating an age reference material, and ρ_{D} is the dosimeter track density (see Hurford, 2019). In this case, c is dimensionless because the spontaneous and induced track count densities are expressed in the same measurement units.
Laserablation FT dating requires a slightly different value for c because no dosimeter glass is involved in parent measurement (see Vermeesch, 2019):
In this case, the dimension of c depends on the units of parent measurement, e.g., as U concentration or as element ratio, e.g., U $/$ Ca.
C1 Mean ages
For datasets showing a single age, it is tempting to report the arithmetic mean age due to its familiarity and simple calculation. However, the mean age is inadequate for summarizing most thermochronological ages. First, calculating a mean from ages determined by a logarithmic age equation as in Eq. (1) “linearizes” the age equation and causes a negative bias compared to applying the logarithmic age equation to a mean D $/$ P ratio. Second, even when directly applied to the ratio, the arithmetic mean gives a biased age estimate, as can be shown from its relationship to the pooled age (see below; Pearson, 1896; Härtel et al., 2022a):
v_{D} and v_{P} are the variation coefficients (standard deviation divided by arithmetic mean) of the daughter and parent concentrations, and r_{DP} is their correlation coefficient. Equation (C1) shows that for the ideal proportional D–P relationship (r_{DP}=1, v_{D}=v_{P}), the mean and pooled ages are the same. In a less ideal case, the measurement error in the parent concentration increases v_{P} and – as it is independent of the daughter concentration – weakens the relationship between D and P (decreasing r_{DP}). This causes the mean age to increase with respect to the pooled age. It means that the mean age is biased towards higher ages under nonideal daughter–parent relationships. This is especially problematic for the wholegrain He and laserablation FT methods, for which the analytical uncertainties are often too small to explain the observed age variation (e.g., Fitzgerald et al., 2006; Ketcham et al., 2018). Essentially, measurement error in the parent concentration creates a rightskewed age distribution, whose mean increases with increasing variance and is biased towards higher ages.
A more robust alternative for calculating a central tendency is the pooled age, i.e., treating all analyzed grains as a single grain by summing up all daughter and parent concentrations. The age is then calculated by substituting the ratio of these sums for D $/$ P in Eq. (1):
Vermeesch (2008) pointed out that in the presence of outliers with a high parent concentration or age, the pooled age is biased towards these grains. Also, Green (1981) and Galbraith and Laslett (1993) argued that the pooled age is not as appropriate as sample age if the age variation cannot be explained by the estimated uncertainties. However, in the case of clustered data (Sect. 3.4.2) or those forming a linear trend with zero intercept (Sect. 3.4.1) without outliers, the age variation is small so that the bias in the pooled age can be assumed to be negligible. The uncertainty in the pooled age can be estimated from error propagation of the singlegrain uncertainties. For He and ZR dating, this gives
with s representing the uncertainties in D, P, and t. FT dating also requires taking into account the uncertainty in c in Eq. (C2). For the EDM method, this gives (Galbraith, 2005)
where N_{s}, N_{i}, and N_{d} are the spontaneous, induced, and dosimeter track counts, respectively, and ζ and s(ζ) are the calibration factor and its uncertainty.
For laserablation FT dating, the uncertainty in the pooled age is
If the ages from a dataset are overdispersed due to the uncertainties not reflecting the variation in the data, it may be advantageous to estimate the uncertainty of the pooled age directly from the variation in D and P concentrations (e.g., Pearson, 1896):
v_{D} and v_{P} represent the variation coefficients of D and P, and r_{DP} is the correlation coefficient for the D–P relationship. Equation (C6) may give a more realistic uncertainty estimate than those in Eqs. (C3)–(C5) if the data are slightly overdispersed. For strongly scattered data, however, Eq. (C6) gives a large uncertainty, confirming that a single sample age may be meaningless.
C2 Isochron ages
For systematically offset data (Sect. 3.4.3), the singlegrain ages and the pooled age are offset in the same direction and give erroneously high or low ages (see Sect. 2.3). In this case, it is advantageous to calculate an isochron age by fitting a regression line to the D–P data and replacing D $/$ P in Eq. (1) by the slope m:
The uncertainty in the isochron age results from propagation of the slope's uncertainty. This logarithmic age equation avoids the bias of the isochron age identified by Vermeesch (2008) for a linear age equation.
Typical algorithms for fitting isochrons are uncertaintyweighted (York, 1968; Kullerud, 1991) and robust regression (Huber, 1981; Powell et al., 2020). Both of these assign weights to each data point: the former based on the measured uncertainty and the latter based on the uncertainty and the distance of each point from a linear “spine” in the data. Robust regression is therefore useful for datasets in which single grains fall off welldefined trends. However, its benefits are limited in the case of many grains deviating from the trend. These regression algorithms, together with the classic leastsquares regression, are implemented in Incaplot.
In general, data at the low and highparent ends of the distribution and data with small uncertainties have a strong influence on the isochron age, making it sensitive to outliers. Its use should therefore be limited to cases of systematic offset in the D–P relationship. Apart from the isochron age, the intercept may also contain important information for the interpretation and should be reported together with the age (Sect. 3.4.3).
The mean square weighted deviation (MSWD or the spine width for robust isochrons) of the isochron provides information on how well the isochron fits the data. An MSWD within the confidence interval (Table C1) indicates that the variation of the data about the isochron is within the range expected from the input uncertainties. A high MSWD outside the confidence interval (Table C1) denotes overdispersed data, whose variation is not explained by the input uncertainties alone – this may point to either unidentified sources of error or intergrain variation of true ages within a sample. For He and laserablation FT data, whose sources of error are not yet well understood, these metrics have to be used with caution.
A standard practice to account for overdispersed data in geochronology is to expand the uncertainty of the isochron age, multiplying it by $\sqrt{\mathrm{MSWD}}$ (e.g., Ludwig, 2012).
C3 Age mixtures
Apart from the simple cases, discrete or continuous mixtures of ages may occur. There are two strategies to deal with discrete age components in a sample (Sect. 3.4.5): mixture modeling (e.g., Galbraith and Laslett, 1993; Galbraith, 2005; Vermeesch, 2019) or splitting the data into different groups and calculating sample ages for each of them.
A continuous age mixture occurs if a sample contains grains with a wide range of kinetic properties responding differently to same thermal history (e.g., Vermeesch, 2019) – each grain then acts as single thermochronometer. An example could be the apatite FT age in a monotonously cooled plutonic rock with grains of different Cl $/$ F ratio. In this case, the intrasample age variation reflects both the measurement error and the trueage variation between grains. This distribution is best described by a “randomeffects model”, and the age to be reported is the central age (Galbraith and Laslett, 1993) – the dispersion parameter describes the variation in true ages. Note, however, that it is necessary to relate the singlegrain age to a kinetic parameter such as grain size, mineral chemistry, or measured radiation damage (Fig. 7b) to justify the use of a continuous mixture of ages. Galbraith (2005) and Vermeesch (2019) provide further discussion and calculation algorithms of the central age for FT dating and Vermeesch (2008) for He dating. For complex data that cannot be described by a discrete or continuous mixture, we suggest reporting the range of singlegrain ages, which requires no additional assumptions.
The synthetic D–P data shown in Figs. 1–8 are available as a Supplement to this article. Incaplot is available as a standalone executable for MacOS and Windows OS and as Python code at https://doi.org/10.5281/zenodo.10637446 (Härtel, 2024).
The supplement related to this article is available online at: https://doi.org/10.5194/gchron64292024supplement.
BH: conceptualization, methodology, formal analysis, visualization, writing (original draft preparation). EE: conceptualization, visualization, funding acquisition, writing (review and editing).
The contact author has declared that neither of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We thank Stephen Cox and an anonymous referee for their thoughtful reviews that greatly improved our manuscript. We also thank Tibor Dunai for editorial handling.
This research was funded by the University of Calgary Eyes High postdoctoral match funding to Birk Härtel, the Natural Sciences and Engineering Research Council of Canada (NSERC; RGPIN202403863 to Eva Enkelmann), and the American Chemical Society – PRF no. 67107ND8 (Eva Enkelmann).
This paper was edited by Tibor J. Dunai and reviewed by Stephen Cox and one anonymous referee.
Ault, A. K., Guenthner, W. R., Moser, A. C., Miller, G. H., and Refsnider, K. A.: Zircon grain selection reveals (de)coupled metamictization, radiation damage, and He diffusivity, Chem. Geol., 490, 1–12, https://doi.org/10.1016/j.chemgeo.2018.04.023, 2018.
Armstrong, E. M., Ault, A. K., Kaempfer, J. M., and Guenthner, W. R.: Connecting visual metamictization to radiation damage to expand applications of zircon (UTh) $/$ He thermochronometry, Chem. Geol., 648, 121949, https://doi.org/10.1016/j.chemgeo.2024.121949, 2024.
Barbarand, J., Carter, A., Wood, I., and Hurford, T.: Compositional and structural control of fissiontrack annealing in apatite, Chem. Geol., 203, 107–137, https://doi.org/10.1016/S00092541(02)004242, 2003.
Baughman, J., Flowers, R. M., Metcalf, J. R., and Dhansay, T.: Influence of radiation damage on titanite He diffusion kinetics, Geochim. Cosmochim. Ac., 205, 50–64, https://doi.org/10.1016/j.gca.2017.01.049, 2017.
Brown, R. W., Beucher, R., Roper, S., Persano, C., Stuard, F., and Fitzgerald, P.: Natural age dispersion arising from the analysis of broken crystals. Part I: Theoretical basis and implications for the apatite (UTh)/He thermochronometer, Geochim. Cosmochim. Ac., 122, 478–497, https://doi.org/10.1016/j.gca.2013.05.041, 2013.
Carter, A.: The thermal history and annealing effects in zircons from the Ordovician of North Wales, Nucl. Tracks Radiat. Meas., 17, 309–313, https://doi.org/10.1016/13590189(90)90051X, 1990.
Carter, A.: Thermochronology on sand and sandstones for stratigraphic and provenance studies, in: Fissiontrack thermochronology and its application to Geology, edited by: Malusà, M. G. and Fitzgerald, P. G., Springer International Publishing, New York, 259–268, https://doi.org/10.1007/9783319894218_14, 2019.
Cogné, N. and Gallagher, K.: Some comments on the effect of uranium zonation on fission track dating by LAICPMS, Chem. Geol., 573, 120226, https://doi.org/10.1016/j.chemgeo.2021.120226, 2021.
Cooperdock, E. H. G., Ketcham, R. A., and Stockli, D. F.: Resolving the effects of 2D versus 3D grain measurements on apatite (U–Th)/He age data and reproducibility, Geochronology, 1, 17–41, https://doi.org/10.5194/gchron1172019, 2019.
Deliens, M., Delhal, J., and Tarte, P.: Metamictization and UPb systematics – a study by infrared absorption spectrometry of Precambrian zircons, Earth Planet. Sc. Lett., 33, 331–344, https://doi.org/10.1016/0012821X(77)900851, 1977.
Dunkl, I.: Trackkey: a Windows program for calculation and graphical presentation of fission track data, Comput. Geosci., 28, 3–12, https://doi.org/10.1016/S00983004(01)000243, 2002.
Fanale, F. P. and Kulp, J. L.: The helium method and the age of the Cornwall, Pennsylvania magnetite ore, Econ. Geol., 57, 735–746, https://doi.org/10.2113/gsecongeo.57.5.735, 1962.
Fitzgerald, P. G., Baldwin, S. L., Webb, L. E., and O'Sullivan, P. B.: Interpretation of (UTh) $/$ He single grain ages from slowly cooled crustal terranes: a case study from the Transantarctic Mountains of southern Victoria Land, Chem. Geol., 225, 91–120, https://doi.org/10.1016/j.chemgeo.2005.09.001, 2006.
Flowers, R. M., Ketcham, R. A., Shuster, D. L., and Farley, K. A.: Apatite (UTh) $/$ He thermochronometry using a radiation damage accumulation and annealing model, Geochim. Cosmochim. Ac., 73, 2347–2365, https://doi.org/10.1016/j.gca.2009.01.015, 2009.
Flowers, R. M., Ketcham, R. A., Enkelmann, E., Gautheron, C., Reiners, P. W., Metcalf, J. R., Danišik, M., Stockli, D. F., and Brown, R. W.: (UTh) $/$ He chronology: Part 2. Considerations for evaluating, integrating, and interpreting conventional individual aliquot data, Geol. Soc. Am. Bull., 135, 137–161, https://doi.org/10.1130/B36268.1, 2022.
Galbraith R. F.: Statistics for fission track analysis, Chapman & Hall, Boca Raton, Florida, 219 pp., ISBN 9780367392796, 2005.
Galbraith, R. F. and Green, P. F.: Estimating the component ages in a finite mixture, Nucl. Tracks Radiat. Meas., 17, 197–206, 1990.
Galbraith, R. F. and Laslett, G. M.: Statistical models for mixed fission track ages, Nucl. Tracks Radiat. Meas., 21, 459–470, https://doi.org/10.1016/13590189(93)90185C, 1993.
Gautheron, C., Djimbi, C. M., Roques, J., Balout, H., Ketcham, R. A., Simoni, E., Pik, R., SeydouxGuillaume, A.M., and TassanGot, L.: A multimethod, multiscale theoretical study of He and Ne diffusio in zircon, Geochim. Cosmochim. Ac., 268, 348–367, https://doi.org/10.1016/j.gca.2019.10.007, 2020.
Green, P. F.: A new look at statistics in fissiontrack dating, Nucl. Tracks, 55, 77–86, https://doi.org/10.1016/0191278X(81)900299, 1981.
Guenthner, W. R., Reiners, P. W., Ketcham, R. A., Nasdala, L., and Giester, G.: Helium diffusion in natural zircon: Radiation damage, anisotropy, and the interpretation of (UTh) $/$ He thermochronology, Am. J. Sci., 313, 145–198, https://doi.org/10.2475/03.2013.01, 2013.
Guenthner, W. R., Reiners, P. W., Hendriks, D., and Tilberg, M.: Zircon, titanite, and apatite (UTh) $/$ He ages and ageeU correlations from the Fennoscandian Shield, southern Sweden, Tectonics, 36, 1254–1274, https://doi.org/10.1002/2017TC004525, 2017.
Härtel, B.: Incaplot (v1.37), Zenodo [code], https://doi.org/10.5281/zenodo.8233941, 2024.
Härtel, B., Jonckheere, R., Wauschkuhn, B., Hofmann, M., Frölich, S., and Ratschbacher, L.: Zircon Raman dating: Age equation and calibration, Chem. Geol., 579, 120351, https://doi.org/10.1016/j.chemgeo.2021.120351, 2021.
Härtel, B., Jonckheere, R., Krause, J., and Ratschbacher, L.: Spurious ageeU associations in thermochronology, Earth Planet. Sc. Lett., 599, 117870, https://doi.org/10.1016/j.epsl.2022.117870, 2022a.
Härtel, B., Jonckheere, R., and Ratschbacher, L.: Multiband Raman analysis of radiation damage in zircon for thermochronology: Partial annealing and mixed signals, Geochem. Geophys. Geosy., 23, e2021GC010182, https://doi.org/10.1029/2021GC010182, 2022b.
Härtel, B., Matthews, W. A., and Enkelmann, E.: Duluth Complex FC1 apatite and zircon: reference materials for (UTh)/He dating?, Geostand. Geoanal. Res., 47, 669–681, https://doi.org/10.1111/ggr.12492, 2023.
He, J., Thomson, S. N., Reiners, P. W., Hemming, S. R., and Licht, C. J.: Rapid erosion of the central Transantarctic Mountains at the EoceneOligocene transition: Evidence from skewed (UTh) $/$ He date distributions near Beardmore Glacier, Earth Planet. Sc. Lett., 567, 117009, https://doi.org/10.1016/j.epsl.2021.117009, 2021.
Heller, B. M., Lünsdorf, N. K., Dunkl, I., Molnár, F., and von Eynatten, H.: Estimation of radiation damage in titanites using Raman spectroscopy, Am. Mineral., 104, 857–868, https://doi.org/10.2138/am20196681, 2019.
Holden, N. E.: Total halflives for selected nuclides, Pure Appl. Chem., 62, 941–958, https://doi.org/10.1351/pac199062050941, 1990.
Holden, N. E., Coplen, T. B., Böhlke, J. K., Tarbox, L. V., Benefield, J., de Laeter, J. R., Mahaffy, P. G., O`Connor, G., Roth, E., Tepper, D. H., Walczyk, T., Wieser, M. E., and Yoneda, S.: IUPAC periodic table of the elements and isotopes (IPTEI) for the education community, Pure Appl. Chem., 90, 1833–2092, https://doi.org/10.1515/pac20150703, 2018.
Holland, H. D. and Gottfried, D.: The effect of nuclear radiation on the structure of zircon, Acta Crystallogr., 8, 291–300, https://doi.org/10.1107/S0365110X55000947, 1955.
Hourigan, J. K., Reiners, P. W., and Brandon, M. T.: UTh zonationdependent alphaejection in (UTh) $/$ He chronometry, Geochim. Cosmochim. Ac., 69, 3349–3365, https://doi.org/10.1016/j.gca.2005.01.024, 2005.
Huber, P. J.: Robust Statistics, John Wiley and Sons, Inc., New York, 305 pp., https://doi.org/10.1002/9780470434697, 1981.
Hurford, A. J.: An historical perspective on fissiontrack thermochronology, in: Fissiontrack thermochronology and its application to Geology, edited by: Malusà, M. G. and Fitzgerald, P. G., Springer International Publishing, New York, 3–24, https://doi.org/10.1007/9783319894218_1, 2019.
Issler, D. R., Grist, A. M., and Stasiuk, L. D.: PostEarly Devonian thermal constraints on hydrocarbon source rock maturation in the Keele Tectonic Zone, Tulita area, NWT, Canada, from multikinetic apatite fission track thermochronology, vitrinite reflectance and shale compaction, Bull. Can. Pet. Geol., 53, 405–431, https://doi.org/10.2113/53.4.405, 2005.
Jaffey, A. H., Flynn, K. F., Glendenin, L. E., Bentley, W. C., and Essling, A. M.: Precision measurement of halflives and specific activities of 235U and 238U, Phys. Rev. C, 4, 1889–1906, https://doi.org/10.1103/PhysRevC.4.1889, 1971.
Kempe, U., Trullenque, G., Thomas, R., Sergeev, S., Presnyakov, S., Rodionov, N., and Himcinschi, C.: Substitutioninduced internal strain and high disorder in weakly radiation damaged hydrothermal zircon from Mt. Malosa, Malawi, Eur. J. Mineral., 30, 659–679, https://doi.org/10.1127/ejm/2018/00302739, 2018.
Ketcham, R. A., van der Beek, P., Barbarand, J., Bernet, M., and Gautheron, C.: Reproducibility of Thermal History Reconstruction from Apatite FissionTrack and (UTh) $/$ He Data, Geochem. Geophys. Geosy., 19, 2411–2436, https://doi.org/10.1029/2018GC007555, 2018.
Kohn, B. P., Ketcham, R. A., Vermeesch, P., Boone, S. C., Hasebe, N., Chew, D., Bernet, M., Chung, L., Danišik, M., Gleadow, A. J. W., and Sobel, E. R.: Interpreting and reporting fissiontrack chronological data, Geol. Soc. Am. Bull., https://doi.org/10.1130/B37245.1, 2024.
Kullerud, L.: On the calculation of isochrons, Chem. Geol., 87, 115–124, https://doi.org/10.1016/01689622(91)90045X, 1991.
Liu, J., Glasmacher, U. A., Lang, M., Trautmann, C., Voss, K.O., Neumann, R., Wagner, G. A., and Miletich, R.: Raman spectroscopy of apatite irradiated with swift heavy ions with and without simultaneous exertion of high pressure, Appl. Phys. A, 91, 17–22, https://doi.org/10.1007/s0033900844029, 2008.
Ludwig, K. R.: Isoplot/Ex Version 3.75: A geochronological toolkit for Microsoft Excel, Special Publication, 4, Berkeley Geochronology Center, 1–75, https://www.bgc.org/isoplot/8857ffdd07944f29a547a23f863a1ffc (last access: 23 July 2024), 2012.
Malusà, M. G.: A guide for interpreting complex detrital age patterns in stratigraphic sequences, in: Fissiontrack thermochronology and its application to Geology, edited by: Malusà, M. G. and Fitzgerald, P. G., Springer International Publishing, New York, 279–293, https://doi.org/10.1007/9783319894218_16, 2019.
Malusà, M. G. and Fitzgerald, P. G.: Application of thermochronology to geologic problems: bedrock and detrital approaches, in: Fissiontrack thermochronology and its application to Geology, edited by: Malusà, M. G. and Fitzgerald, P. G., Springer International Publishing, New York, 191–209, https://doi.org/10.1007/9783319894218_10, 2019.
Martin, P. E., Metcalf, J. R., and Flowers, R. M.: Calculation of uncertainty in the (U–Th) ∕ He system, Geochronology, 5, 91–107, https://doi.org/10.5194/gchron5912023, 2023.
Miltich, L.: Low temperature cooling history of Archean gneisses and Paleoproterozic granites of southwestern Minnesota, B.A. thesis, Carleton College, Minnesota, 58 pp., 2005.
Murray, K. E., Orme, D. A., and Reiners, P. W.: Effects of UThrich grain boundary phases on apatite helium ages, Chem. Geol., 390, 135–151, https://doi.org/10.1016/j.chemgeo.2014.09.023, 2014.
Nasdala, L., Irmer, G., and Wolf, D.: The degree of metamictization in zircon: A Raman spectroscopic study, Eur. J. Mineral., 7, 471–478, https://doi.org/10.1127/ejm/7/3/0471, 1995.
Nicolaysen, L. O.: Graphic interpretation of discordant age measurements on metamorphic rocks, Ann. N. Y. Acad. Sci., 91, 198–206, https://doi.org/10.1111/j.17496632.1961.tb35452.x, 1961.
Orme, D. A., Reiners, P. W., Hourigan, J. K., and Carrapa, B.: Effects on inherited cores and magmatic overgrowths on zircon (UTh) $/$ He ages and ageeU trends from Greater Himalayan sequence rocks, Mount Everest region, Tibet, Geochem. Geophys. Geosy., 16, 2499–2507, https://doi.org/10.1002/2015GC005818, 2015.
Pearson, K.: Mathematical contribution to the theory of evolution. On a form of spurious correlation which may arise when indices are used in the measurement of organs, Proc. R. Soc. Lond., 60, 489–498, 1896.
Pickering, J., Matthews, W., Enkelmann, E., Guest, B., Sykes, C., and Koblinger, B. M.: Laser ablation (UThSm)/He dating of detrital apatite, Chem. Geol., 548, 119683, https://doi.org/10.1016/j.chemgeo.2020.119683, 2020.
Powell, R., Green, E. C. R., Marillo Sialer, E., and Woodhead, J.: Robust isochron calculation, Geochronology, 2, 325–342, https://doi.org/10.5194/gchron23252020, 2020.
Reiners, P. W. and Farley, K. A.: Influence of crystal size on apatite (UTh)/He thermochronology: an example from the Bighorn Mountains, Wyoming, Earth Planet. Sci. Lett., 188, 413–420, https://doi.org/10.1016/S0012821X(01)003417, 2001.
Ritter, W. and Märk, T. D.: Optical studies of radiation damage and its annealing in natural fluorapatite, Nucl. Instr. Meth. Phys. Res. B1, 394–397, https://doi.org/10.1016/0168583X(84)900983, 1984.
Schmitz, M. D. and Bowring, S. A.: UPb zircon and titanite systematics of the Fish Canyon Tuff: an assessment of highprecision UPb geochronology and its application to young volcanic rocks, Geochim. Cosmochim. Ac., 65, 2571–2587, https://doi.org/10.1016/S00167037(01)006160, 2001.
Shuster, D. L., Flowers, R. M., and Farley, K. A.: The influence of natural radiation damage on helium diffusion kinetics in apatite, Earth. Planet. Sc. Lett., 249, 148–161, https://doi.org/10.1016/j.epsl.2006.07.028, 2006.
Spiegel, C., Kohn, B., Belton, A., Berner, Z., and Gleadow, A.: Apatite (UThSm)/He thermochronology of rapidly cooled samples: The effect of He implantation, Earth Planet. Sc. Lett., 285, 105–114, https://doi.org/10.1016/j.epsl.2009.05.045, 2009.
Steiger, R. H. and Jäger, E.: Subcommission on geochronology: convention on the use of decay constants in geo and cosmochronology, Earth Planet. Sc. Lett., 36, 359–362, https://doi.org/10.1016/0012821X(77)900607, 1977.
TripathyLang, A., Hodges, K. V., Monteleone, B. D., and van Soest, M. C.: Laser (UTh) $/$ He thermochronology of detrital zircons as a tool for studying surface processes in modern catchments, J. Geophys. Res.Earth 118, 1333–1341, https://doi.org/10.1002/jgrf.20091, 2013.
Troch, J., Ellis, B. S., Schmitt, A. K., Bouvier, A.S., and Bachmann, O.: The dark side of zircon: textural, age, oxygen isotopics and trace element evidence of fluid saturation in the subvolcanic reservoir of the Island ParkMount Jackson Rhyolite, Yellowstone (USA), Contrib. Min. Petrol., 173, 54, https://doi.org/10.1007/s0041001814812, 2018.
Vermeesch, P.: Three new ways to calculate average (UTh) $/$ He ages, Chem. Geol., 249, 339–347, https://doi.org/10.1016/j.chemgeo.2008.01.027, 2008.
Vermeesch, P.: HelioPlot, and the treatment of overdispersed (UThSm)/He data, Chem. Geol., 271, 108–111, https://doi.org/10.1016/j.chemgeo.2010.01.002, 2010.
Vermeesch, P.: Statistics for fissiontrack thermochronology, in: Fissiontrack thermochronology and its application to Geology, edited by: Malusà, M. G. and Fitzgerald, P. G., Springer International Publishing, New York, 109–122, https://doi.org/10.1007/9783319894218_6, 2019.
Vermeesch, P. and Tian, Y.: Thermal history modelling: HeFty vs. QTQt, EarthSci. Rev., 139, 279–290, https://doi.org/10.1016/j.earscirev.2014.09.010, 2014.
Vermeesch, P., Seward, D., Latkoczy, C., Wipf, M., Günther, D., and Bauer, H.: αemitting mineral inclusions in apatite, their effect on (UTh) $/$ He ages, and how to reduce it, Geochim. Cosmochim. Ac., 71, 1737–1746, https://doi.org/10.1016/j.gca.2006.09.020, 2007.
Wendt, I. and Carl, C.: The statistical distribution of the mean squared weighted deviation, Chem. Geol., 86, 275–285, https://doi.org/10.1016/01689622(91)90010T, 1991.
Wernicke, R. S. and Lippolt, H. J.: Botryoidal hematite from the Schwarzwald (Germany): heterogeneous uranium distributions and their bearing on the helium dating method, Earth Planet. Sc. Lett., 114, 287–300, https://doi.org/10.1016/0012821X(93)900314, 1993.
Willett, C. D., Fox, M., and Shuster, D. L.: A heliumbased model for the effects of radiation damage annealing on helium diffusion kinetics in apatite, Earth Planet. Sc. Lett., 477, 195–204, https://doi.org/10.1016/j.epsl.2017.07.047, 2017.
York, D.: Least squares fitting of a straight line with correlated errors, Earth Planet. Sc. Lett., 5, 320–324, https://doi.org/10.1016/S0012821X(68)800597, 1968.
Zeigler, S. D., Metcalf, J. R., and Flowers, R. M.: A practical method for assigning uncertainty and improving the accuracy of alphaejection corrections and eU concentrations in apatite (U–Th) ∕ He chronology, Geochronology, 5, 197–228, https://doi.org/10.5194/gchron51972023, 2023.
 Abstract
 Introduction
 Background
 Classificationbased workflow for data analysis using the D–P plot
 Conclusions
 Appendix A: The effective uranium concentration
 Appendix B: Units of daughter and parent concentrations
 Appendix C: Age calculation and reporting
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Background
 Classificationbased workflow for data analysis using the D–P plot
 Conclusions
 Appendix A: The effective uranium concentration
 Appendix B: Units of daughter and parent concentrations
 Appendix C: Age calculation and reporting
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement