Uncertainty in helium diffusion in zircon limits thermochronologic resolution: application to the Great Unconformity

Thermochronology provides a unique perspective on the timing and magnitude of erosion during the generation of unconformities. Recently, thermochronology has been used to reinvigorate a long-standing debate about the origin of the Great Unconformity, a geological feature representing a period of missing time of almost a billion years from the stratigraphic record at the end of the Precambrian. Our ability to resolve thermal histories is fundamentally limited by how well we know parameters describing the temperature sensitivity of a thermochronometric system. The (U-Th)/He in zircon system is particularly well suited to examine the origin of the Great Unconformity because it is very sensitive to long durations of time at relatively low temperatures (< 200–250 °C). Here we determine uncertainties in the Zircon Radiation Damage and Annealing Model (ZRDAAM, Guenthner et al. 2013) that describes changes in ⁴He diffusion kinetics as a function of radiation damage accumulation and annealing. We show that the dispersion in predicted zircon (U-Th)/He ages for a given thermal history can be 100's of Ma for a specific amount of radiation damage highlighting that thermal histories are less well resolved than previously appreciated. Additional diffusion experiments and calibration with natural laboratories would provide better constraints on diffusion kinetic parameters.

1 Introduction

Thermochronometry is widely used to constrain the evolution of Earth's surface and upper crust, transforming our understanding of the magnitudes of sedimentation and erosion linked to climate and tectonic change (Reiners and Brandon, 2006; McDannell and Flowers, 2020; Gautheron et al., 2022). The principle behind thermochronometry is that rocks experience temperature changes over geological timescales: rocks closer to Earth's surface are cooler than deeper ones and so exhumation leads to cooling. The time scales associated with these changes in temperature are determined using the concepts of geochronology and diffusion. The daughter products of radioactive decay used for thermochronometry have a temperature dependent rate of loss from the target mineral. At high temperatures the daughter products are effectively lost instantaneously. By contrast, as the rock cools to lower temperatures, the rate of diffusive loss decreases and daughter products are progressively retained until there is effectively no diffusive loss. Therefore a measured age reflects the duration of residence at low temperature (Dodson, 1973; Zeitler et al., 1987; Reiners, 2005; Fox and Carter, 2020), but many thermal histories, that may incorporate reheating events, will produce the same age.

One of the most widely used methods for deep-time thermochronometry is (U-Th)/He in zircon (ZHe) (Reiners, 2005). The basis of this method is that the radioactive elements uranium and thorium are incorporated into crystals of zircon and the decay product, helium, is trapped in the crystal lattice. Helium (⁴He) diffuses from zircon at high temperatures, but is retained at lower temperatures. The exact temperature range at which the transition from open-system to closed-system behaviour occurs is dependent on damage to the crystal lattice produced primarily from recoil during decay processes (Guenthner et al., 2013; Ketcham et al., 2013). Radiation damage accumulates at a rate that depends on the amount of uranium (U) and thorium (Th) in the crystal and may anneal if temperatures increase. This means that two zircon crystals from the same rock, experiencing the same thermal history, could have very different thermal sensitivity. In turn, models accounting for this variable temperature sensitivity as a function of radiation damage can be used to leverage more complex thermal histories than constant kinetic parameter models.

Thermochronometry has been used recently to understand the development of the Great Unconformity (e.g., DeLucia et al., 2018; Flowers et al., 2020; McDannell et al., 2022), a global feature that marks the boundary between the Precambrian and Phanerozoic. The time-period that the Great Unconformity spans varies depending on location, but on the North American craton, in the Grand Canyon, erosion removed stratigraphy between 1300 to 250 million years (Thurston et al., 2022). The origin of the Great Unconformity has been debated for over 125 years, with recent contributions using thermochronometry to gain new insight (Flowers et al., 2020; McDannell et al., 2022; Thurston et al., 2022). This approach has resulted in papers that support the hypothesis that the Great Unconformity is the result of glacial erosion (Keller et al., 2019), and also that it is the result of diachronous tectonic events (Flowers et al., 2020).

Both those in favour of a globally, reasonably, synchronous glacial origin (McDannell et al., 2022) and regionally diachronous tectonic origins (Flowers et al., 2020) of the Great Unconformity have interpreted ZHe data using the Zircon Radiation Damage And Annealing Model (ZRDAAM) of Guenthner et al. (2013). Importantly, it is clear that it is always better to combine multiple thermochronometric systems as this helps resolve thermal histories and reduces the negative effects of systematic uncertainties for specific systems (McDannell et al., 2022). In ZRDAAM model, a crystal is composed of undamaged and damaged parts that combine to give bulk diffusion kinetics as a function of the amount of radiation damage. Accumulated radiation damage is calculated based on the concentrations of the parent elements and the thermal history of a sample. The damage accumulates at low temperatures, but can be annealed at higher temperatures, calculated using fission track annealing kinetics (Yamada et al., 1995; Rahn et al., 2004; Tagami, 2005; Yamada et al., 2007). Therefore, the diffusivity of helium at a specific temperature is a function of the past thermal history. This makes the overall problem very non-linear so that changing the temperature at some time in the thermal history can have unexpected effects on the resulting age, as also shown for the (U-Th)/He in apatite system by Fox and Shuster (2014).

Using the zircon radiation damage and annealing model (ZRDAAM) with inverse models, researchers have resolved tight temperature constraints on thermal histories over billion-year timescales.

Recently, contrasting interpretations of Great Unconformity evolution in North China have highlighted the sensitivity of inferred thermal histories to modelling assumptions. Zhan et al. (2026) interpreted multichronometer datasets as evidence for predominantly Paleoproterozoic exhumation, whereas McDannell et al. (2026) argued that thermal overprinting limits the ability of ZHe data to discriminate between competing scenarios. These contrasting conclusions raise broader questions concerning the propagation of uncertainty through diffusion kinetic models used in thermal-history inversion.

In the Eastern Grand Canyon, Thurston et al. (2022) inferred a 1700 Ma thermal history from ZHe ages. Parts of this history were reported to within less than 10 °C between 700 and 250 Ma and then again from 15–7 Ma. It is unclear whether the data really provide such tight constraints on temperatures in the past or whether these are at least partly the consequence of not sufficiently sampling parameter space, model assumptions and/or, potentially, overconfidence in the adopted diffusion kinetic parameters. More generally, recent discussions have highlighted that inferred thermal histories may depend strongly on modelling assumptions and the way geological constraints are incorporated into inverse models (Flowers and Peak, 2025). Here, we focus specifically on uncertainty in the diffusion kinetic parameters underpinning ZHe thermal-history reconstruction. Specifically, we want to know how uncertainty in the estimates of diffusion kinetic parameters may impact on the inference of thermal histories and the associated uncertainties.

Calibration of ZRDAAM has been carried out using measured diffusion kinetics of crystals with well constrained amounts of radiation damage (Guenthner et al., 2013). However, the accuracy and precision of this model has not been assessed. In particular, it is unclear how the propagation of uncertainties to model parameters affects the dispersion or sensitivity of predicted thermochronometric ages. Here we show that the uncertainties in the radiation damage model make it challenging to accurately infer the timing and magnitude of unconformities in the deep past. We begin by highlighting why ZRDAAM needs to be calibrated by accounting for uncertainties and present a new calibration based on the same underlying data. We then propagate uncertainties from this model calibration into time temperature path uncertainty.

2 The existing calibration of the radiation damage model

The rate of diffusion is controlled by the diffusivity and the curvature of the concentration of the diffusant (Fick's Law). Although the production distribution of the diffusant (helium in our case) can be important in some scenarios, diffusion tends to smooth the distribution. More significant is the fact that the diffusivity can vary by orders of magnitude with variations in temperature. The diffusivity at any temperature is given by the Arrhenius equation of the diffusivity D (cm² s⁻¹):

$$\begin{array}{}\text{(1)}& D\left(T\right)=D_{\mathrm{0}}{e}^{\frac{-E_{\mathrm{a}}}{RT}}\end{array}$$

where D₀ is the frequency factor (cm² s⁻¹), E_a is the activation energy (kJ mol⁻¹), R is the gas constant ($\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{mol}}^{-\mathrm{1}}$) and T is the temperature in Kelvin. Taking the logarithm of Eq. (1), gives:

$$\begin{array}{}\text{(2)}& \ln \left(D\right(T\left)\right)={\displaystyle \frac{-E_{\mathrm{a}}}{R}}{\displaystyle \frac{\mathrm{1}}{T}}+\ln \left(D_{\mathrm{0}}\right)\end{array}$$

so that the slope of the line between ln (D(T)) and $\mathrm{1}/T$ gives $E_{\mathrm{a}}/R$ and the intercept of the line provides ln (D₀). Diffusion experiments in which a crystal is step-wise degassed in vacuo are used to calculate D(T) for different combinations of specific temperatures and time (Fechtig and Kalbitzer, 1966). The resulting plot can be used to determine the Arrhenius parameters (D₀,E_a). Importantly, estimates of the two model parameters extracted from this linear inversion covary with one another, i.e., ln (D₀) is strongly correlated with E_a.

Analyzing different zircon crystals with known radiation damage values allows us to assess how diffusion kinetics vary with damage. As it is challenging to visualize both model parameters (D₀ and E_a) for each crystal as a function of radiation damage it is common to combine the parameters and calculate a diffusivity at a specific temperature or a closure temperature at a specific cooling rate and grain size. By combining the parameters, however, information on how the two parameters (D₀ and E_a) are correlated is lost. In turn, our ability to accurately infer thermal histories is reduced.

The results of Guenthner et al. (2013)'s diffusion experiments highlight two general trends. At low damage values the closure temperature increases with increasing damage. At higher damage values, the closure temperature decreases with increasing damage. This general behaviour has been reproduced in numerical models conducted at a range of scales (Ketcham et al., 2013; Gautheron et al., 2020). To interpret these trends in ZRDAAM, a model is used in which the diffusion kinetics for a specific radiation damage value are a combination of a theoretical minimally damaged crystal and an extremely damaged crystal. The diffusion kinetics of these end-member crystals need to be estimated. The frequency factor of the minimally damaged, theoretical, crystal, (^zD₀) was estimated by extrapolating the frequency factors of measured crystals down two orders of magnitudes using a power-law relationship. The activation energy for the theoretical crystal (^zE_a) is set as the average of the activation energies of minimally damaged crystals (see Guenthner et al., 2013 for details). Extrapolating values to a minimally damaged crystal, however, will add uncertainties and there is no obvious way to account for these in the power-law relationship. Crucially, this approach does not account for the correlations between the model parameters. This is important because the correlations provide additional information that can yield more precise estimates of model parameters and allow propagation of uncertainties into model predictions. The diffusion kinetics for the extremely damaged crystal are estimated using sample N17 (Guenthner et al., 2013), and also involve correlated model parameters (^N17D₀ and ^N17E_a). However, the accuracy of this model has only been assessed by looking at general trends in model predictions (Guenthner et al., 2013). Here we attempt to formally quantify the uncertainty in these damage model parameters and how these uncertainties translate to uncertainties in temperature sensitivity of the ZHe system, and in particular predicted ZHe ages. We note that additional work, informed by our exploration, could be carried on quantifying uncertainties in other aspects of the ZRDAAM.

3 A new calibration of the zircon radiation damage and annealing model

In order to account for the correlation between the frequency factor and the activation energy, we model the measured helium diffusivities directly. We use the same diffusion dataset, alpha damage values and parameterisation as Guenthner et al. (2013). However, in contrast with Guenthner et al. (2013), we determine the diffusion of the end member crystals using the radiation damage model directly, rather than by non-linear extrapolation from high to low radiation damage levels. We focus on the Guenthner et al. (2013) model and data, as opposed to the revised Guenthner (2021) model, as the original model was used in the papers focused on Great Unconformity. We also highlight that our approach can be adopted to constrain any damage model with different datasets. Our goal is not to simply increase the accuracy of the model parameters but to determine their precision. By tracking correlations in model parameters, that are the result of correlations within individual diffusion experiments, we can simulate (U-Th)/He in zircon ages accounting for uncertainties in the original diffusion experiments.

The model we fit through each diffusion data set is given by equation 8 of Guenthner et al. (2013) and describes the diffusivity D as a function of the amount of damage:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{1}}{D_e/a^{\mathrm{2}}}}={\displaystyle \frac{f_{\mathrm{c}}^{\prime }}{\frac{\mathrm{1}}{{\left(\frac{l_{\mathrm{int},\mathrm{0}}}{l_{\mathrm{int}}}\right)}^{\mathrm{2}}}\left(\frac{D_z}{(a{f}_{\mathrm{c}}^{\prime }{)}^{\mathrm{2}}}\right)}}+{\displaystyle \frac{f_{\mathrm{a}}^{\prime }}{\frac{D_{\mathrm{N}\mathrm{17}}}{(a{f}_{\mathrm{a}}^{\prime }{)}^{\mathrm{2}}}}}\end{array}$$

where $f_{\mathrm{c}}^{\prime }$ and $f_{\mathrm{a}}^{\prime }$ are the crystalline and amorphous fractions, respectively, l_int and l_int,0 are parameters describing how far a helium atom can travel within a crystal lattice without encountering damage in a minimally damaged and extremely damaged crystal, respectively, and a is the grain size, or more accurately, the equivalent spherical radius. D_z and D_N17 can be calculated using the diffusion kinetics of the undamaged and damaged theoretical crystals, using Eq. (4). Therefore, for every sample with an estimated amount of damage, we can calculate different diffusivity values for degassing steps using four model parameters ^N17D₀ and ^N17E_a for D_N17 and ^zD₀ and ^zE_a for D_z. In this way, the model predicts the diffusivity at a specific temperature of the data using the model parameters and we do not rely on linear regressions through individual Arrhenius relationships. The log likelihood function is defined as the negative least-squares misfit when fitting all the diffusion data for a given set of model parameters. However, the degassing experiments of Guenthner et al. (2013) each have different numbers of steps. The likelihood function involves a summation, related to the number of steps, and so number of data points in each experiment, experiments with more steps would tend to dominate our results. Similarly, key experiments which might have fewer steps would have far less influence on the model parameter estimation. To account for this problem, the least squares misfit for each experiment is weighted accordingly, so that the log-likelihood (LL) function is:

$$\begin{array}{}\text{(4)}& \text{LL}=-\mathrm{0.5}\cdot \sum _{i=\mathrm{1}}^N{\displaystyle \frac{\mathrm{1}}{M_i}}\sum _{j=\mathrm{1}}^{M_i}{\left({\displaystyle \frac{D_{i,j}-P_{i,j}}{\mathit{\sigma }}}\right)}^{\mathrm{2}}\end{array}$$

where N is the number of crystals analyzed, M_i is the number of degassing steps used in the inversion for a specific crystal, D_i,j is the observed diffusivity for a specific crystal at a specific degassing step, P_i,j is corresponding predicted diffusivity calculated with Eq. (4), and σ is the estimated uncertainty set to  $\mathrm{1}\ln (\mathrm{1}/s)$ here, based on reported uncertainties.

We use the Bayesian Markov Chain Monte Carlo (MCMC) method incorporating the Metropolis-Hastings algorithm to sample the full posterior distribution of the model parameters (^N17D₀ and ^N17E_a for D_N17 and ^zD₀ and ^zE_a). We tune the proposal distributions (which determine how far to perturb model parameters between iterations) to ensure that approximately 20 % of the proposed models are accepted as this represents an efficient balance between exploring parameter space and sampling the parameter values (Gelman et al., 1997; Gallagher, 2012). The Markov Chain is initialised with the model parameters of Guenthner et al. (2013) and the algorithm runs until 1 million sets of model parameters have been accepted. In this way, at each iteration we: (1) propose four model parameters by perturbing the current model parameters; (2) calculate an Arrhenius relationship for each diffusion experiment using the four proposed model parameters and calculate the likelihood; (3) accept the proposed model parameters if they improve the likelihood or accept with a probability determined by the ratio of the likelihoods; (4) update the current model parameters to the proposed model parameters if they are accepted. The advantage of this approach over simply extracting diffusion parameters from individual experiments is that that uncertainty on diffusion data and model correlations are propagated through to uncertainty on model parameters. Preliminary experiments highlighted that the sampling was relatively insensitive to the diffusion kinetics of the extremely damaged N17 crystal. To ensure that the diffusion kinetics of this unique and crucial end member crystal was accurately captured, we reduced the uncertainty of the diffusion data for N17 to $\mathrm{0.1}\ln (\mathrm{1}/s)$. If this uncertainty is not decreased, the frequency factor and activation energy of this crystal are not predicted within error.

Figure 1 Key parameters controlling how radiation damage controls diffusivity have been inferred from fitting Arrhenius relationships from step-degassing experiments. The datapoints are shown by circles and the lines represent model fit. Both are coloured by the alpha dose. (A) The fit to the step-degassing experiments for the model parameters inferred by Guenthner et al. (2013). (B) The fit to the step-degassing experiments using model parameters extracted from the data using a Markov Chain Monte Carlo analysis. The model parameters are the minimum misfit model parameters and represent a single realization of the parameter set we infer. A key point of the figure is to show the fits to the data and the data that are used to predict the model parameters, and to highlight that both models do not predict all the data, highlighting the source of uncertainty.

Figure 2 The sampled marginal posterior distributions for the four diffusion parameters representing the two hypothetical crystals. (A) and (B) are the frequency factor (^N17D₀) and activation energy (^N17E_a), respectively, for an extremely damaged crystal. (C) and (D) are the frequency factor (^zD₀) and activation energy (^zE_a), respectively, for a minimally damaged crystal. The red lines show the values of these parameters used by Guenthner et al. (2013).

Table 1 Covariance matrix of the four diffusion-kinetic parameters. Rows and columns are ordered identically. Diagonal terms are variances and off-diagonal terms are covariances.

The comparison of the predicted and observed diffusivities for the best fitting model is shown in Fig. 1B. We also show the model fit with the original Guenthner et al. (2013) parameters in Fig. 1A. From these results, it is clear that neither set of model parameters fit the data perfectly and this is expected with noisy data. However, our model parameters do provide a better fit to the data based on our misfit criterion in Eq. (5). Results of our analysis are plotted as 4 histograms showing the original model parameters and our inferred model parameters (Fig. 2). These are histograms of model parameters sampled by our MCMC algorithm. These are known as the marginal posterior distributions for each model parameter. The probability of the model parameters is proportional to the height of the model histograms. For the extremely damaged crystal, our maximum a posteriori model parameter values are close to the original values, shown in red. However, for the low-damaged crystal, the model estimates are quite different, although we note that the original values fall along the strong correlation trend seen in the 2D probability distribution. The original values are a little displaced from the peak (Fig. 3) showing that these parameters do not fit the data quite as well, but the combination of values are consistent with the linear relationship between ^zE_a and ^zD₀. We note also that while there is a strong correlation between ^N17D₀ and ^N17E_a and between ^zD₀ and ^zE_a the 2 sets of parameters are not strongly correlated (i.e. ^N17D₀ and ^N17E_a are effectively independent of ^zD₀ and ^zE_a). Approximating the posterior distribution as a four-dimensional Gaussian distribution provides the mean values of the parameters, the covariance matrix and the correlation matrix. Below, the parameters are ordered as ^N17E_a, log 10(^N17D₀), ^zE_a and log 10(^zD₀). The mean values for the four parameters are: 5.61 × 10⁴ J mol⁻¹, −3.75 log (cm² s⁻¹), 1.42 × 10⁺⁵ J mol⁻¹, 3.54 log (cm² s⁻¹) with a model covariance matrix shown in Table 1.

Figure 3 Inferred model parameters from diffusion data are strongly correlated. Our approach to infer the diffusion kinetics of the hypothetical crystals using the radiation damage model maintains this correlation. This is clearly illustrated with the diffusion parameters for the undamaged crystal. The pink spot shows the diffusion kinetics inferred by Guenthner et al. (2013) and while it is reasonably far from the center of our inferred distribution, it still falls on the clear correlation trend defined by the sampling. The colours are proportional to posterior probability with oranges/yellows reflecting highest probabilities.

Table 2 Correlation matrix corresponding to the covariance matrix in Table 1.

In addition to the covariance matrix, it is useful to inspect to the correlation matrix as it easier to see how parameter values correlate without the variable units. The correlation matrix shows strong correlations across some parameters but very low correlations across other parameters, Table 2.

The correlation matrix reveals an almost perfect positive correlation between activation energy and frequency factor within each endmember crystal. Correlation coefficients are 0.999 for the N17 endmember and 0.991 for the low-damage z endmember. In contrast, the diffusion parameters of the two endmembers are essentially independent, with correlation coefficients close to zero. This indicates that uncertainty in the Arrhenius parameters of one endmember has little influence on uncertainty in the parameters of the other endmember.

4 Propagating model uncertainties

To assess the importance of the uncertainties in the radiation damage and annealing model parameters, we predict ages using a simple thermal model. The time-temperature path is chosen to resemble that of the Minnesota samples from Miltich (2005) and McDannell et al. (2022), and represents a typical inferred time temperature path of a Deep Time target locality. Here the rocks have been below 600 °C since 1.5 Ga. From 700 to 650 Ma the rocks cooled from 200 to 150 °C, and then gradually to 0 °C by 200 Ma before experiencing reheating at 50 Ma to 100 °C. Between 50 Ma and the present the rocks cooled linearly to 0 °C. Radiation damage accumulates throughout this history such that some zircon crystals that transitioned from open to closed system behaviour during the cooling event at 700 Ma, transitioned back to open behaviour simply due to the accumulation of radiation damage. Some crystals however, with intermediate temperature sensitivity only record the final cooling event. Other crystals, with low-temperature sensitivity, also record cooling associated with the 100 Ma burial event. In terms of constraining thermal history models, this potential to have a wide range of temperature sensitivities within a single sample makes the ZHe method very powerful.

To calculate thermochronometric ages we use the radiation damage and annealing model of Guenthner et al. (2013) with our updated model parameters. Note, the numerical methods used to implement the diffusion equations and damage models have been used previously to interpret zircon ${}^{\mathrm{4}}\mathrm{He}{/}^{\mathrm{3}}\mathrm{He}$ thermochronometric data (Tripathy-Lang et al., 2015). Twenty different crystals are simulated spanning an effective U concentration ([eU] = [U]+0.24[Th]; Cooperdock et al., 2019; Gastil et al. (1967)) interval from 31 to 2828 ppm.

Initially, the Guenthner et al. (2013) model parameters were used to simulate ages. Grain sizes are all set to 70 µm. A continuous age-[eU] relationship is produced with no variability in age at a specific [eU] value reflecting the lack of uncertainty in the model parameters (Fig. 4).

Figure 4 Propagating uncertainties in the radiation damage model produces a wide range of ZHe ages for a specific amount of damage. The colours highlight the relative probability of obtaining a specific age for a given [eU] value. The red line shows the predicted age-[eU] relationship using the canonical values of the radiation damage and annealing model (Guenthner et al., 2013). The continuous thermal history used to produce the result is shown in the inset.

Next, uncertainties from our radiation damage model calibration were propagated through to predict age distributions. Twenty different crystal ages with different [eU] values were calculated 1000 times with different model parameters for ^N17D₀, ^N17E_a, ^zD₀ and ^zE_a . To do this, we extracted every 10th model from the posterior ensemble of models generated during the MCMC algorithm. This ensures that we are sampling model parameter space in proportion to probability but also that the model correlations are reliably captured. As before the grain sizes were all set to 70 µm.

The results highlight the large spread in predicted ages for a representative thermal history, typical for exploiting the ZRDAAM (Fig. 4). The overall spread in age between the youngest age and the oldest age is expected given the different temperature sensitivity of the crystals due to radiation damage. However, even for a specific amount of radiation damage there is still a large dispersion in the predicted ages. For example, at [eU] values of about 1600 ppm, ages are expected to vary between 50 and 550 Ma. If a range of grain sizes were also modelled for a specific [eU], the spread would be even larger (Whipp et al., 2022). Furthermore, parent isotope zonation, bad neighbours, broken grains, inclusions or variable pre-deposition thermal hisotries will also contribute to this dispersion (see Fox et al., 2019 for a discussion).

5 Incorporating model parameter uncertainty in inverse models

To assess the effect of imperfectly defined kinetic parameters in the ZHe radiation damage model of Guenthner et al. (2013), we modified QTQt to allow sampling of the 4 parameters (^N17D₀, ^N17E_a, ^zD₀ and ^zE_a) from the joint posterior distribution. The resampling draws values from the 4-dimensional Gaussian with the empirical mean vector and covariance matrix calculated from the joint posterior obtained, as described above. The assumption of Gaussian distributions is not strictly true given the form of the sampled distributions, particularly log (^N17D₀) and ^N17E_a, butthe goal is to sample over the range of values appropriately, and then this approximation is good enough. We sample from a zero mean, unit variance normal distribution, use a Cholesky decomposition of the covariance matrix to rescale those values and then add the means, using the triangular factorization method described in Barr and Slezak (1972). Note the covariance matrix incorporates the correlation between the parameters and thus so do the sampled values. We used age uncertainties of 25 % for the ZHe age values reported in McDannell et al. (2022) as this uncertainty allowed us to reproduce the characteristics of the thermal history reported in McDannell et al. (2022). We ran QTQt for 2 million iterations with the first million being discarded as burn-in. We ran the model with the original kinetic parameters (no resampling) and with resampling of the 4 kinetic parameters from the joint posterior every 500 iterations (and over the 2 million iterations this resampling effectively reproduces the original posterior distribution of these mode parameters). The resampling reproduces the distributions and correlations for the 4 kinetic parameters. However, we do not attempt to infer the optimal value of these parameters, instead we simply incorporate this model parameter uncertainty into the inverse thermal history models and their predictions. Importantly for a given thermal history we use the same radiation damage parameter values for all the grains but due to the different degrees of damage, the resulting diffusion parameters are different for each crystal. The results are given in Fig. 5a and b, which has a summary of the inferred thermal histories. The modelled ages predict the observed ages for the expected models as shown in McDannell et al. (2022): the expected model without resampling has a log posterior value of −376.4 and the expected model with resampling has log posterior value of −351.5. We see that the expected and maximum likelihood thermal histories for the two different kinetic parameter models have a similar form in terms of the timing of cooling events, some variation in the magnitude of cooling and a different range for the credible intervals around the first period of cooling (800–700 Ma). The most significant differences are in the form of the maximum posterior models and the distribution of accepted thermal histories. In particular, the resampling model suggests a second mode or group of solutions that do not have a second period of cooling (around 400–350 Ma). This does not indicate that the solution has not converged, instead the solution has converged to reveal these two groups. We note that, for the “no resampling” model, the maximum likelihood and maximum posterior models are similar, while these two models tend to fall on two separate modes for the resampling model. In this latter case, the maximum posterior model is a little simpler (fewer time-temperature points) and does not reproduce the observations (log likelihood $=-\mathrm{192.1}$) quite as well as the no resampling model (log likelihood = −164.8). A of the model fits to the data is shown in Table 3. The effect of resampling the kinetic parameters from their posterior distribution is to add uncertainty to the inferred thermal history models. This is perhaps not as significant as we might think as the pairs (D and E_a) of kinetic parameters are very strongly correlated (Mialhe et al., 1988). Then, resampling tends to produce combinations of the kinetic parameters that provide similar predictions for a given thermal history. However, the differences are still enough to add dispersion to the distribution of the thermal histories and increase uncertainty. The degree of additional uncertainty will depend on how the data are distributed in age and accumulated radiation damage.

Figure 5 Thermal histories inferred using QTQt with (A) the Guenthner et al. (2013) diffusion parameters and (B) resampling of the 4 diffusion kenetics parameters. The dataset used for this analysis is reported in McDannell et al. (2022). The resampling leads to a wider distribution of possible thermal histories, wider credible intervals, and also identifies 2 different forms, or modes in the posterior distribution (as seen by the Max Posterior and Max. Likelihood models). Note in this case, the expected, or average, model falls between the two modes and so makes less good predictions relative to the other two models (see log-likelihoods in Table 3).

Table 3 Model comparisons for no resampling and resampling damage parameters. In both inversions, the expected (weighted average) model does not fit the data quite as well as the max likelihood and max posterior models due to the averaging process tending to underestimate maximum temperatures during reheating when the timing of reheating is not well resolved. The log posterior values depend on the prior and therefore, there is no value for the expected models. The likelihoods for the max likelihood and max posterior models with no resampling are similar to those of the resampling models. However, they are slightly better because, for the same number of iterations and no resampling of damage parameters, the inversion can more easily find the best fitting models.

6 Implications

Our method to propagate ZRDAAM uncertainty highlights how variable the age-[eU] relationship might be for a given thermal history. In particular, our results suggest that the uncertainty of ZRDAAM-based thermal history inversions may be significantly underestimated. This may have major implications for our ability to differentiate between subtle differences in temperature at specific times. In turn, it may be challenging to resolve cooling histories sufficiently to attribute the Great Unconformity to Cryogenic Glaciations or geodynamic processes using thermochronology alone.

The potential to underestimate age uncertainty for thermal modelling has been discussed by McDannell et al. (2022) and to some extent this can be accounted for in the inverse modelling software QTQt (Gallagher, 2012). For example, if two dates have the same [eU] but their measured uncertainties do not overlap, QTQt can sample additional uncertainty for the measurement age to allow for this excess dispersion. However, if the two ages do not have the same [eU] concentration, the situation is more difficult. There are two options. Either additional uncertainty can be assigned to the measurements by resampling a scaling factor (> 1) that multiplies the input errors and allows the predicted age-[eU] relationship to pass through the observed data+resampled uncertainty. Alternatively, the thermal history can be adjusted to change the predicted age-[eU] relationship to try and ensure that the predictions fit the data, at least to within the error. The first option tends to produce simpler thermal histories than the second option, as the data fitting criterion is less strict. For example, McDannell et al. (2022)'s results for Pikes Peak highlight how models that ignore overdispersion appear to resolve a 700 Ma cooling signature, which is smoothed out when the overdisperion is effectively reduced by adding excess uncertainty on some of the data. Importantly, Flowers et al. (2022) argued that the data do not have the ability to measure the timing of cooling irrespective of how overdispersion is treated.

We have shown the continuous spread of ages as a function of [eU] as a probability heat map for a specific history (Fig. 4). In reality, most thermochronometric studies typically analyse 5–30 crystals for each sample. We can illustrate the effect of model uncertainty by comparing two simulated datasets of 15 ages generated from the same thermal history. We produce these datasets by sampling our age probability distribution randomly (Fig. 6c). For this example, each crystal in the dataset might have a different combination of model parameters drawn from the parameters we have estimated. Although these two datasets display overall similarities, there are subtle differences between their age-[eU] relationship over specific [eU] values. To accurately capture the spread in age for a single radiation damage value, many more thermochronometric samples would need to be collected. To illustrate this point, we draw random samples from the probability distribution in Fig. 4, for [eU] values ranging from 1500 to 2000 ppm, and investigate the spread in age (Fig. 6). The dispersion of the ages varies greatly with increasing sample size, converging to the predicted frequency distribution of Fig. 4. In our specific example, the distribution stabilises for sample sizes of 40 or more crystals. The need to accurately capture spread are especially important if ages need to be averaged within [eU] bins to find acceptable paths as the uncertainty for the mean age is determined by the standard deviation (Flowers et al., 2020; Peak et al., 2021; Thurston et al., 2022). Ault et al. (2018) showed that simple visual identification under the microscope of the degree of metamictization is useful for obtaining good [eU] coverage and this approach could be adopted to ensure that multiple ages for the same [eU] are measured to get an idea of the spread in age. Zonation of eU concentration can also complicate the interpretation of thermochronometric data (Anderson et al., 2017; Weisberg et al., 2018; Anderson et al., 2018), potentially because different zones can evolve to have completely different diffusion kinetics (Fox et a., 2014). One promising approach to estimate the amount of zonation and anomalous diffusion behaviour could be to screen crystals using ramped degassing experiments (Idleman et al., 2018; Guo et al., 2024; McDannell et al., 2018). Alternatively, zircon ${}^{\mathrm{4}}\mathrm{He}{/}^{\mathrm{3}}\mathrm{He}$ thermochronology (Tripathy-Lang et al., 2015) could be exploited to identify spatial variability in ⁴He and diffusivity.

Figure 6 Model realizations and expected ranges of ages. (A, B) random samples are drawn from the probability distribution in Fig. 4 to highlight the sorts of datasets that are expected given the typical number of ages measured on a single sample. The simulated ages are different between the two realizations of a typical dataset. The gray regions show the [eU] range inspected in panel (C). (C) Many ages need to be sampled in order to accurately capture the spread in ages over a specific [eU] bin.

The large uncertainties on the parameters controlling helium diffusion in zircon and the potentially dramatic impact this has on temperature sensitivity highlights that this aspect is important to consider. Currently, it may not be practical to incorporate diffusion kinetic uncertainties in inverse models directly because this dramatically increases the volume of the parameter space that needs to be searched and could lead to significantly longer run times with current software. However, with the development of faster computers and parallelized inverse methods, this will not be a major obstacle. More importantly, to ensure that we sample crystals with the same [eU] values to resolve diffusion kinetic parameters, many more grain ages per sample need to be analysed. For example, to accurately capture the spread in age for a relatively narrow [eU] range of 1500–2000 ppm, 40 crystals from this interval were required (Fig. 6), and therefore for a sample with a wide range of eU values many multiples of this number would be required. A practical solution to avoid measuring so many crystals per sample and running millions of simulations in an inversion is to use forward modelling. To do this, a single thermal history that is close to what might be expected for a specific area given prior geological knowledge could be used to assess expected age spread. This expected age spread could then be added to the age uncertainties used for inverse modelling. This procedure could be iterated to produce realistic uncertainties.

ZRDAAM has been calibrated using a limited number of diffusion experiments. Additional work is required to develop this dataset to capture diffusion kinetics at different radiation damage values (Ginster, 2018). These experiments could also aim to replicate diffusion kinetics at previously measured radiation damage values to quantify the degree of dispersion. In addition, natural laboratories could be utilised to resolve diffusion parameters: areas with known thermal histories can be exploited to predict ZHe ages by varying diffusion parameters; complementary thermochronometers can be leveraged to find thermal histories and diffusion parameters that match the observed data. Ultimately, by reducing the uncertainties in helium diffusion kinetics using the constraints from man-made and natural laboratories, the timings of cooling events in the past can be resolved with improved accuracy and precision.