Technical note: Disentangling pre- and post-depositional thermal histories in partially reset samples

Interpreting thermochronological data of a partially reset detrital sample is challenging because its age distribution reflects a mixture of each grain's individual pre-depositional and the shared post-depositional thermal history. A promising approach to meet this challenge is combined geo- and thermochronological dating on the same grains. We present an algorithm for disentangling the pre- and post-depositional thermal histories using data from zircon $\mathrm{U}/\mathrm{Pb}$–(U-Th) $/$ He double-dating. It differs from previous approaches in that it doesn't make any assumptions about the pre-depositional thermal histories of the grains in their source regions. First, we determine candidate post-depositional temperature-time paths by inverse thermal-history modeling of syn-depositional grains that crystallized at or near the time of deposition. We then calculate pre-depositional (U-Th) $/$ He model ages for each candidate thermal history representing cooling of the grains in their source regions. Subsequently, we evaluate which post-depositional thermal histories are plausible using the likelihood of their pre-depositional model ages falling in the expected interval between the $\mathrm{U}/\mathrm{Pb}$ crystallization ages of each grain and the time of sediment deposition. The final result is a selection of high-likelihood pairs of post-depositional thermal histories and pre-depositional ages. We illustrate this strategy by applying it to a Devonian sandstone from the Northern Canadian Cordillera. The results show the general agreement of the modeled thermal histories and pre-depositional (U-Th) $/$ He ages with existing thermochronological data. The post-depositional thermal histories are also consistent with additional zircon-Raman thermochronological data from the same grains. We discuss the limitations of our approach imposed by the need for syn-depositional grains, its application to other thermochronometers, and how targeting rather than avoiding partially reset samples enhances our capability to extract thermal-history information from the detrital record. Finally, we introduce PaRACAS, the implementation of our approach in a jupyter notebook that calculates, plots, and evaluates the pre-depositional zircon (U-Th) $/$ He model ages for post-depositional temperature-time paths determined in HeFTy.

1 Introduction

Detrital thermochronology is a powerful tool to date thermal events in sedimentary basins and in their source regions (e.g., Garver et al., 1999; Malusà and Fitzgerald, 2020). The information we obtain from detrital thermochronology depends on the sample's post-depositional history as shown in Fig. 1: an unreset sample (blue line in Fig. 1a) stayed at temperatures lower than the partial retention zone (PRZ). Its grains did not lose radioactive-decay products (daughters) during post-depositional heating. Its single-grain ages thus pre-date deposition (Fig. 1b); their variation mainly reflects different pre-depositional thermal histories and provides valuable information about cooling in the source regions (e.g., Garver et al., 1999; Rahl et al., 2003; Enkelmann et al., 2008; Bernet and Garver, 2005; Malusà and Fitzgerald, 2020). A fully reset sample was reheated well above the PRZ, erasing prior age information (yellow line in Fig. 1a). Its single-grain age distribution displays one age peak post-dating deposition (Fig. 1b). Minor within-sample age differences occur due to thermal-sensitivity differences between grains, e.g., related to properties such as grain size, composition, or radiation damage (Carlson et al., 1999; Reiners and Farley, 2001; Flowers et al., 2009; Guenthner et al., 2013). Reset samples allow us to constrain the time of post-depositional cooling in the basin (e.g., Cecil et al., 2014; Enkelmann et al., 2014; McKay et al., 2021).

Figure 1 Schematic thermal histories and single-grain age distributions of unreset (blue), partially reset (green), and fully reset (yellow) sedimentary samples in the zircon (U-Th) $/$ He (ZHe) system. (a) Example thermal histories for detrital grains with pre-depositional (grey) temperature-time (T–t) paths that represent simplified cooling paths for each grain and shared post-depositional paths (colored). The grey-shaded area marks the partial retention zone (PRZ; kinetics of Reiners et al., 2004, for non-metamict zircon grains with equivalent spherical radii of 20–90 µm and 30 Myr holding time). (b) Kernel-density estimate (KDE) curves for synthetic age distributions reflecting the thermal histories in (a).

In contrast, a partially reset sample experienced post-depositional heating into the PRZ (green line in Fig. 1a). It retained only part of its accumulated daughters after deposition and the age distribution spans across the depositional age (Fig. 1b). The single-grain ages vary due to both differences in pre-depositional thermal histories and differences between the thermal sensitivity of individual grains. With increasing post-depositional temperatures, the age spectrum approaches that of a fully reset sample. During holding at very high or low temperatures in the PRZ, the within-sample thermal-sensitivity differences even enable the co-existence of partially reset grains with unreset or fully reset grains.

In theory, partially reset samples contain information on both the pre- and post-depositional history that is inaccessible for fully reset or unreset samples. The mere notion of a sample having spent time within the well-defined temperature window of the PRZ provides a more precise estimate on the timing and temperature of maximum heating rather than just an upper (unreset) or lower (fully reset) temperature limit. Correcting the observed thermochronological ages for post-depositional heating would allow to extract the time when the source rocks cooled (e.g., Reiners et al., 2005) – information beyond the reach of fully reset samples. However, accurate interpretation of thermochronological data from partially reset detrital samples is complicated by the single-grain ages representing a mixture of the shared post-depositional and the individual pre-depositional thermal histories (Reiners et al., 2005; Flowers et al., 2022).

The most promising approach for disentangling pre- and post-depositional thermal histories is combined thermochronological and geochronological dating of the same grains (e.g., Carter and Moss, 1999; Reiners et al., 2005; Enkelmann et al., 2019b; Bootes et al., 2019; Danišik, 2019; Dunkl et al., 2024). In the popular zircon $\mathrm{U}/\mathrm{Pb}$–(U-Th) $/$ He double-dating scheme, the $\mathrm{U}/\mathrm{Pb}$ ages – commonly interpreted as crystallization ages – provide individual constraints on the start of the thermal history for each grain that can be used for forward and inverse thermal-history modeling of each grain (Reiners et al., 2005; Fosdick et al., 2015; Yonkee et al., 2019). Nevertheless, the existing interpretation strategies for partially reset samples based on double-dating rely on assumptions about the pre-depositional thermal histories. These include the time of initial source-rock cooling or the duration and magnitude of heating events in the source regions. That knowledge is often not available, especially for pre-Cenozoic sediments because the specific source lithologies have been eroded. An alternative approach is therefore necessary, particularly in the light of the in-situ ZHe–$\mathrm{U}/\mathrm{Pb}$ double-dating technique by laser ablation, which allows acquiring large numbers of detrital-zircon double-dates (e.g., Vermeesch et al., 2012; Tripathy-Lang et al., 2013; Evans et al, 2015; Forte et al., 2022; Zawacki et al., 2024; Dunkl et al., 2024).

Building on the approaches of Reiners et al. (2005) and Fosdick et al. (2015), we present a new interpretation strategy for combined ZHe and $\mathrm{U}/\mathrm{Pb}$ data from partially reset samples. Its main innovation is that it avoids assumptions regarding the pre-depositional history by relying on syn-depositional grains, i.e., volcanic or shallow magmatic zircon grains that crystallized at or near the time of deposition. Our approach proceeds in three steps: (1) finding candidate temperature-time (T–t) paths for the post-depositional history by inverse T–t modeling of the syn-depositional grains; (2) calculating pre-depositional model ages for each grain and each candidate T–t path; (3) evaluating the plausibility of post-depositional thermal histories based on the likelihood of the pre-depositional model ages to fall between the individual crystallization ages and the depositional age. The results of this procedure are selected pairs of post-depositional thermal histories and pre-depositional model ages that can be interpreted in terms of source-region cooling. We first explain the details of our new approach on a set of synthetic data and compare it to existing interpretation strategies (Fig. 2). We then apply it to a partially reset Devonian sandstone from the Northern Canadian Cordillera. The results show that our strategy allows extraction of meaningful information on the pre- and post-depositional histories from ancient sediments. We discuss (1) the limits of our method, including the requirement of syn-depositional grains, problems with nearly fully-reset samples, and simplifying assumptions in He-diffusion modeling, (2) the application of our modeling strategy to other double-dating schemes, and (3) the question if and when it is necessary to model partially reset samples. Lastly, we introduce PaRACAS (Partial Reset Analysis for Cooling Ages of Sources; Härtel, 2026), a jupyter notebook that provides the functions for determining, plotting, and evaluating pre-depositional model ages from double-dating results and the outputs of inverse T–t modeling.

Figure 2 Overview of partial-reset-modeling strategies applied to a synthetic zircon $\mathrm{U}/\mathrm{Pb}$–(U-Th) $/$ He dataset. All information on this dataset is provided at the top panel. (a) True temperature-time (T–t) paths for each grain. The grey-shaded area marks the zircon (U-Th) $/$ He (ZHe) partial retention zone (grain size of the synthetic grains; kinetic parameters of Reiners et al., 2004) and (b) plot of (α-ejection-corrected) ZHe age vs. $\mathrm{U}/\mathrm{Pb}$ age. The panels below show the three different interpretation strategies of Reiners et al. (2005; c–e), Fosdick et al. (2015; f–h) and this study (i–k), each in three analytical steps. Detailed description of the specific steps for each strategy in the text.

2 Interpretation strategies for partially reset samples using combined zircon $\mathrm{U}/\mathrm{Pb}$–(U-Th) $/$ He dating

While it is clear when fully reset or unreset samples accumulated the bulk of their daughters, partially reset samples contain a mixture of pre- and post-depositional daughters. Interpreting partially reset zircon (U-Th) $/$ He (ZHe) data with standard inverse thermal-history modeling (e.g., Ketcham, 2005, 2024; Gallagher, 2012) starting at the depositional age suffers from two biases: (1) for samples with many ZHe ages pre-dating deposition, the models underestimate post-depositional temperature – to reach a good fit with the measured ages all T–t paths need to ensure accumulation and retention of enough He in a restricted time frame. (2) If most ZHe ages post-date deposition, matching the T–t paths with the measured data leads to overestimation of post-depositional temperature toward a full-reset scenario. A widely-used approach is forward modeling the thermal history with an “inheritance envelope” of pre-depositional residence near the surface or a pre-depositional cooling rate (e.g., Fillon et al., 2013; Cecil et al., 2014; Guenthner et al., 2015; Powell et al., 2016, 2020; Fox et al., 2019; Yonkee et al., 2019). Its main advantage is that it allows accumulating not only pre-depositional He but also radiation damage that affects the thermal sensitivity of the ZHe system (Hurley et al., 1956; Reiners et al., 2005; Guenthner et al., 2013; Gérard et al., 2022). This approach evaluates the model fit by comparing the predicted relationship between ages and effective uranium concentrations (eU) to the measured data. However, this approach requires at least endmember assumptions about the pre-depositional thermal history, which is usually unknown and varies from grain to grain.

Combining ZHe and $\mathrm{U}/\mathrm{Pb}$ dating on the same detrital grains imposes constraints on the pre- and post-depositional histories, as the $\mathrm{U}/\mathrm{Pb}$ of each grain marks the starting point of its thermal history. For example, the observed amount of He loss since crystallization is a proxy for the degree of reset (Reiners et al., 2005), and each grain's $\mathrm{U}/\mathrm{Pb}$ age is a high-temperature starting point for thermal-history modeling (Fosdick et al., 2015). Figure 2 illustrates three interpretation strategies on a synthetic $\mathrm{U}/\mathrm{Pb}$–ZHe dataset (Table S1 in the Supplement; n=10) with $\mathrm{U}/\mathrm{Pb}$ ages ranging from 50 to 205 Ma and grain sizes (equivalent spherical radii) between 38 and 122 µm. The grains underwent different thermal histories before deposition at 45 Ma (Fig. 2a; crosses: true pre-depositional ZHe ages) and were heated together to a maximum temperature of 168 °C at 25 Ma (thermal-history data in Table S2). The resulting (α-ejection-corrected) ZHe ages (31 to 124 Ma) span across the depositional age (dashed line) and correlate positively with the $\mathrm{U}/\mathrm{Pb}$ ages (Fig. 2b). This relationship reflects the ZHe-age variation due to different pre-depositional thermal histories.

Reiners et al. (2005) present a first-order interpretation approach to estimate the maximum-heating temperature for a simple post-depositional thermal history and a correction approach to obtain ZHe ages of pre-depositional cooling (Fig. 2c–e). Their first step is forward modeling the ZHe ages for a set of post-depositional thermal histories (Fig. 2c) assuming every grain stayed at surface temperatures between crystallization and deposition (grey horizontal line on the left). Figure 2d plots the measured ZHe ages against the $\mathrm{U}/\mathrm{Pb}$ ages with fields of modeled ages based on the maximum temperatures in Fig. 2c and the measured grain-size range. The further interpretation relies on the fraction of total He loss f_total:

$$\begin{array}{}\text{(1)}& f_{\mathrm{total}}=\mathrm{1}-{\displaystyle \frac{t_{\mathrm{ZHe}}-t_{\mathrm{reset}}}{t_{\mathrm{UPb}}-t_{\mathrm{reset}}}},\end{array}$$

where t_ZHe and t_UPb are the ZHe and $\mathrm{U}/\mathrm{Pb}$ ages, respectively, and t_reset is the time of post-depositional resetting (here: 25 Ma).

The grain plotting at the lowest temperature in Fig. 2d is the grain with the smallest total He loss (f_total=f_min). f_min also denotes the maximum possible He loss that can have occurred to all grains together after deposition, indicating the maximum post-depositional temperature. The grain with the lowest total He loss in Fig. 2d (arrow) suggest that heating reached between 160 and 180 °C, an interval containing the true maximum temperature of 168 °C.

f_min also allows correcting all ZHe ages for the maximum possible post-depositional He loss and determining a pre-depositional model ZHe age t_model:

$$\begin{array}{}\text{(2)}& t_{\mathrm{model}}=t_{\mathrm{reset}}+{\displaystyle \frac{t_{\mathrm{ZHe}}-t_{\mathrm{reset}}}{\mathrm{1}-f_{\min }}}.\end{array}$$

Figure 2e plots these model ages against the $\mathrm{U}/\mathrm{Pb}$ ages (grey circles) and compares them to the “true” pre-depositional ages calculated from forward modeling the paths in Fig. 2a up to the depositional age (pink crosses). While the pre-depositional ages are a first approximation, six out of ten grains deviate from the true pre-depositional ages by more than 10 Myr – one model age even post-dates the timing of deposition (dashed line). This discrepancy results from extrapolating the annealing kinetics of the least-reset grain to the other grains. If the least-reset grain is more retentive than the average grain, Eq. (2) underestimates the pre-depositional model ages (Fig. 2e). Reiners et al. (2005) discuss two additional limitations: first, the model assumes that the grains cooled rapidly after crystallization and stayed near the surface until deposition. This neglects pre-depositional residence in the PRZ. f_min is thus a maximum estimate of the post-depositional reset of the least-reset grain and the model age may still overestimate the time of pre-depositional cooling. This is particularly a problem for samples with many grains post-dating deposition. Second, depending on the relationship between the depositional age and the measured ZHe ages, the assumed value of t_reset may have a strong impact on the pre-depositional ages. Nevertheless, this interpretation approach is a simple method for narrowing down the maximum post-depositional temperature. It is also independent of a specific He-diffusion model, as long as the diffusion kinetics are sufficiently similar across grains.

In contrast, Fosdick et al. (2015) present an approach that models a combination of both pre- and post-depositional thermal histories. It uses forward modeling of T–t paths that consist of a pre-depositional segment (Fig. 2f, left of the depositional-age line), to which every grain rapidly cooled immediately after crystallization (grey vertical lines), and a post-depositional heating pulse (blue path, Fig. 2f). Fosdick et al. (2015) generate a large set of T–t paths with different maximum temperatures, timings and durations of residence in the PRZ and model the present-day single-grain ZHe ages for each path. They evaluate the modeled and observed ZHe age distributions with a two-sample Kolmogorov-Smirnov (K-S) test (cumulative-distribution plot in Fig. 2g), keeping only the best-fitting paths for model interpretation (K-S p value >68 %; blue and black lines). Finally, they compare the pre-and post-depositional histories of these best-fitting models (Fig. 2h). For our synthetic data, the seven best-fitting T–t paths cluster at maximum temperatures of 160 to 180 °C, which they reached between 40 and 15 Ma and held for 2–13 Myr. Their timing and maximum temperature of post-depositional heating agrees to a first degree with the input thermal history (pink line). The time of pre-depositional cooling through the PRZ for these models spans a large range from 200 to about 100 Ma containing most of the true pre-depositional ages (pink vertical ticks).

One of this method's drawbacks is that it does not evaluate pre-depositional cooling on a grain-by-grain basis, but only for the bulk sample. Thus, sources with different cooling histories contributing only few grains to the sample will be overlooked because the best-fitting paths are selected based on the bulk goodness-of-fit. Nevertheless, this approach is more flexible than that of Reiners et al. (2005) for testing different pre- and post-depositional histories and considers kinetic differences between grains.

The main limitation of both modeling approaches is that they assume either a specific pre-depositional history for every grain or at least uniform pre-depositional histories across grains. We therefore present an alternative strategy (Fig. 2i–k). Instead of making assumptions about the pre-depositional thermal history, it focuses exclusively on the post-depositional history shared by all grains. To this end, it makes use of two fundamental principles: (1) syn-depositional grains, i.e., grains that crystallized at or near the time of deposition, record only the post-depositional thermal history (Fosdick et al., 2015), and (2) cooling through the PRZ in the source areas happened after crystallization but before deposition. Our strategy combines the probabilistic evaluation of T–t paths of Fosdick et al. (2015) and the pre-depositional model ages of Reiners et al. (2005). Its first step consists in identifying syn-depositional or near syn-depositional grains based on the $\mathrm{U}/\mathrm{Pb}$ and the depositional ages (blue circles in Fig. 2b). We use these grains for inverse thermal-history modeling in HeFTy (Ketcham, 2005, 2024) to obtain candidate post-depositional T–t paths (Fig. 2i). In the second step, we evaluate each of these paths by calculating pre-depositional model ages. To this end, we determine the net amount of helium accumulated after deposition, He_post, and the remaining fraction ϕ_post of pre-depositional He by forward modeling the candidate post-depositional T–t paths (details in Appendix A). Analogous to Eq. (2), we calculate the pre-depositional He concentration He_pre from the present-day concentration He_obs, He_post and ϕ_post:

$$\begin{array}{}\text{(3)}& {\text{He}}_{\mathrm{pre}}={\displaystyle \frac{\left({\text{He}}_{\mathrm{obs}}-{\text{He}}_{\mathrm{post}}\right)}{{\mathit{\varphi }}_{\mathrm{post}}}}+{\text{He}}_{\mathrm{post}}.\end{array}$$

The pre-depositional ZHe model age then follows from a generic age equation (Härtel and Enkelmann, 2024) and adding the depositional age t_dep for scaling:

$$\begin{array}{}\text{(4)}& t_{\mathrm{model}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{238}}}}\ln \left(\mathrm{1}+{\displaystyle \frac{{\text{He}}_{\mathrm{pre}}}{{\text{eU}}_{\mathrm{pre}}}}{\displaystyle \frac{M_{\mathrm{238}}}{\mathrm{8}}}\right)+t_{\mathrm{dep}}.\end{array}$$

M₂₃₈, λ₂₃₈, and 8 are the molar mass of ²³⁸U, its decay constant, and the number of α-disintegrations in its decay chain, respectively. eU_pre is the pre-depositional effective uranium concentration, i.e., the concentrations of U and Th at the timing of deposition recalculated to the integrated decay rate of ²³⁸U between t_model and t_dep.

We evaluate each path based on the likelihood of the pre-depositional model ages falling into the expected envelopes between the individual $\mathrm{U}/\mathrm{Pb}$ ages and the depositional age (details in Appendix B). Figure 2j plots the model ages against the $\mathrm{U}/\mathrm{Pb}$ ages, color-coded by their likelihood. Each vertical string of squares represents the set of model ages of one grain for the different T–t paths. The higher-likelihood paths (dark squares) mostly fall between the equal-age (solid) and depositional-age (dashed) lines. For further interpretation, it is practical to focus on a subset of the highest-likelihood post-depositional thermal histories and their pre-depositional model ages. The likelihoods and their distribution vary from one dataset to another, depending on the number of dated grains, the age distribution and uncertainties, and the extent of reset. Finding a universally applicable algorithm for selecting such a subset based on a number or proportion of paths, or an absolute-likelihood threshold is therefore difficult. We apply a “knee-detection” criterion to determine a cutoff likelihood (dotted line in Fig. 2k; kneedle algorithm of Satopää et al., 2011). This criterion uses the point at which the likelihood difference between a path and the next-less-likely path drops significantly. The selection therefore contains the paths most consistent with crystallization and deposition constraints spanning a range of likelihoods, while the non-selected paths have low likelihoods and differ little in their fit with the constraints. Thereby, we select the highest-likelihood paths with the likelihood spread still allowing to distinguish how well they agree with the constraints, while we exclude less-likely paths that are more difficult to distinguish based on their likelihood. For the synthetic dataset, the knee-detection algorithm suggests four thermal histories (pink circles). Their post-depositional T–t paths (thick black lines in Fig. 2i) reach maximum temperatures of 160 to 190 °C between 25 and 20 Ma, close to the true thermal history (pink line). These best-fitting T–t paths also predict pre-depositional ages (black squares in Fig. 2j) matching the true pre-depositional ages within 10 Myr (pink crosses).

Building on the strategy of Fosdick et al. (2015) to evaluate a large range of possible T–t paths, our strategy uses Monte-Carlo-based inverse thermal-history modeling to find plausible post-depositional thermal histories that would be tedious to implement in pure forward-modeling approaches. In addition, our approach provides pre-depositional model ages on a grain-by-grain basis like that of Reiners et al. (2005). Interpretation of these pre-depositional ages can proceed in a similar way as combined $\mathrm{U}/\mathrm{Pb}$–ZHe dating results of unreset detrital samples to interpret source-cooling populations based on their age-age relationship (e.g., Reiners et al., 2005; Fosdick et al., 2024), for example as: (1) clusters of grains that overlap in both $\mathrm{U}/\mathrm{Pb}$ and ZHe age, suggesting that they crystallized and cooled together (e.g., Dunkl et al., 2024); (2) grains from shallow magmatic or volcanic sources, for which the ratio of $\mathrm{U}/\mathrm{Pb}$ to ZHe ages is approximately unity (e.g., Reiners et al., 2005; Saylor et al., 2012); (3) arrays of grains sharing the same crystallization age but with different ZHe ages, pointing to a similar source but having experienced different thermal histories (e.g., Reiners et al., 2005), or (4) grain populations showing the same ZHe age, but different $\mathrm{U}/\mathrm{Pb}$ ages, e.g., due to recycling from a fully reset (meta-)sedimentary source or erosion of grains from an igneous source and its contact aureole (e.g., Rahl et al., 2003; Saylor et al., 2012; Malusà et al., 2022). In addition, the pre-depositional model-age populations may show relationships to other single-grain provenance data such as grain size, mineral chemistry, or ϵHf values (e.g., von Eynatten and Dunkl, 2012; Sundell and Saylor, 2021). The good fit of our models to the synthetic dataset confirms the validity of this approach in theory. In the following section, we illustrate its use on actual thermochronological data from the Northern Canadian Cordillera.

Figure 3 Zircon $\mathrm{U}/\mathrm{Pb}$–(U-Th) $/$ He (ZHe) data and modeling results for an Imperial Formation sandstone from the Northern Canadian Cordillera. (a) Plot of measured ZHe age (corrected for α-ejection; see Härtel et al., 2026) vs. $\mathrm{U}/\mathrm{Pb}$ concordia ages with equal-age (solid), and depositional-age (dashed) lines. Blue symbols: syn-depositional grains. Bars show 2 s uncertainties. (b) Temperature–time (T–t) paths from inverse modeling of syn-depositional grains color-coded by log-likelihood. Light-blue boxes: Modeling constraints (Table S3). Black lines: best-fitting model paths. (c) Sorted log-likelihoods of model paths. The inset shows the region of the 20 highest-likelihood paths. Dotted line: cutoff for the selected best-fitting models (pink circles). (d) Plot of pre-depositional model ZHe ages vs. $\mathrm{U}/\mathrm{Pb}$ concordia ages for the selected paths color-coded by log-likelihood. Lines and uncertainty bars as in (a). (e) T–t paths of the selected paths (black lines) compared to envelopes of good-fitting paths from existing thermal-history models (Enkelmann et al., 2019a, blue lines; Powell et al., 2020, grey-shaded area; Spalding et al., 2023; red lines) and the ZHe partial-retention zone (PRZ; purple-shaded area; as in Fig. 1). N: number of samples.

3 Example: zircon $\mathrm{U}/\mathrm{Pb}$–(U-Th) $/$ He–Raman data from a Devonian sandstone, Northern Canadian Cordillera

We applied our new modeling strategy to zircon $\mathrm{U}/\mathrm{Pb}$–(U-Th) $/$ He–Raman multi-chronometer data from a Devonian sandstone (Imperial Formation) from the Mackenzie Plain, Northern Canadian Cordillera (Härtel et al., 2026; their sample NW-15). The data were acquired in three steps: Raman spectroscopic radiation-damage mapping of the ν₃ (∼ 1008 cm⁻¹) and external rotation (ER, ∼356 cm⁻¹) bands, He extraction from a small volume within the grains by laser-ablation noble-gas mass spectrometry (e.g., Boyce et al., 2006; Vermeesch et al., 2012; Evans et al., 2015), and LA-ICP-MS for determining the U and Th concentrations – to calculate the ZHe and zircon Raman (ZR) ages – and the $\mathrm{U}/\mathrm{Pb}$ age. For details on the analytical procedure and data treatment, see Härtel et al. (2026). Figure 3a shows the combined ZHe and concordant zircon $\mathrm{U}/\mathrm{Pb}$ data. $\mathrm{U}/\mathrm{Pb}$ ages range from 364 Ma to 2.9 Ga. ZHe ages span a range between 57 and 2030 Ma, across the depositional age of 364.3±5.7 Ma (dashed line in Fig. 3a; maximum-depositional age of Härtel et al., 2026); about 73 % of the grains post-date the depositional age. The ZHe ages also show a weak positive correlation with the $\mathrm{U}/\mathrm{Pb}$ ages. These observations confirm a high degree of partial reset by heating into the ZHe PRZ after deposition as suggested by thermal-history models on other Imperial-Formation samples from Northern Canada (Issler et al., 2005; Enkelmann et al., 2019a; Powell et al., 2020; Spalding et al., 2023).

For the inverse thermal-history modeling, we selected the three youngest grains with concordant $\mathrm{U}/\mathrm{Pb}$ ages <400 Ma (uncorrected ZHe ages: 106.1±20.1 Ma (2 s); 157.5±25.7 Ma; 266.7±40.3 Ma; blue circles in Fig. 3a), allowing for individual pre-depositional thermal histories of ≤36 Myr, less than 10 % of the post-depositional history. We performed inverse thermal-history modeling in HeFTy (Ketcham, 2005, 2024) with input uncertainties of 20 % to ensure finding T–t paths compatible with all three grains despite their ZHe-age differences. Appendix C describes how we adjusted the modeling procedures to account for differences between whole-grain He-diffusion kinetics and the diffusion kinetics of the volume that the He was extracted from by laser-ablation. We used the radiation-damage-independent diffusion parameters of Reiners et al. (2004) because the dose-equivalent α-radiation-damage densities of our grains (20 to 120×10¹⁶α g⁻¹, measured by Raman spectroscopy) are well within the range for radiation-damage-independent He diffusion (∼5 to 200×10¹⁶α g⁻¹, Guenthner et al., 2013). Table S3 lists the T–t constraints we set for the inverse thermal-history models, including the geologic observation of Lower Cretaceous strata unconformably overlaying the Imperial Formation. We set the stopping condition for the Monte-Carlo-based T–t inversion in HeFTy to 250 good-fit paths (goodness-of-fit >0.5).

For the resulting good-fit paths (Fig. 3b), we calculated the pre-depositional model ages and the path likelihoods. In theory, the syn-depositional grains permit full reset by heating to temperatures >190 °C before 200 Ma (yellow lines), but these paths have low log-likelihoods because they are incompatible with the crystallization ages of the older grains. Most of the other paths (blue and black lines) indicate heating during the Late Paleozoic to Early Mesozoic into the interval between 130 and 180 °C, followed by a second heating event to 50–160 °C in the Late Cretaceous to Paleogene. Out of the 250 paths, the knee-detection algorithm suggests six best-fitting paths (pink circles in Fig. 3c; thick black lines in Fig. 3b and e). Two of these paths show a first heating phase up to 140 °C in the Triassic, followed by ≥50 Myr residence in the ZHe PRZ (purple-shaded area) and cooling in the Jurassic to Cretaceous. Two paths reach their peak temperatures of 140 to 180 °C only in the Late Jurassic to Early Cretaceous before cooling nearly instantaneously. The two final paths do not show heating into the ZHe PRZ before the Late Cretaceous. For the Late-Cretaceous to Paleogene heating event, the best-fitting paths indicate peak temperatures of 100 to 160 °C between 90 and 70 Ma, or about 95 °C at 50 Ma.

Figure 3d plots the pre-depositional model ZHe ages for the six best-fitting paths against the $\mathrm{U}/\mathrm{Pb}$ ages of our grains. Overall, about 83 % of the grains (25 out of 30 grains) show the expected behavior, falling between the crystallization age (solid line) and the deposition age (dashed line). The main group of pre-depositional model ZHe ages falls in the Neoproterozoic to Paleozoic range, with a single grain straddling the equal-age line at 2.2 Ga. Two of the remaining pre-depositional ages are higher than the respective $\mathrm{U}/\mathrm{Pb}$ ages, indicating that the diffusion model overestimates their degree of reset relative to the measured ZHe age. One of them is a syn-depositional grain, whose misfit in Fig. 3d reflects the difference between its analytical uncertainty and the expanded uncertainty used for the inverse T–t modeling. Three pre-depositional ages post-date deposition because the diffusion model underestimates their degree of reset. The misfits between modeled and expected ages reflect several simplifications of the diffusion model, which we discuss in Sect. 4.1. Still, the 83 % of the grains overlapping with the expected envelope, confirm that a simple diffusion model can retrieve consistent T–t information from detrital data with a high degree of partial reset.

Figure 4 Zircon Raman (ZR) data and pre-depositional model ages for the Imperial-Formation sandstone. (a–b) ZR vs. $\mathrm{U}/\mathrm{Pb}$ age-age plots for the external rotation (ER; a) and ν₃ (b) Raman bands. (c–d) Plots of pre-depositional model ZR vs. $\mathrm{U}/\mathrm{Pb}$ ages for the ER (c) and ν₃ (d) bands color-coded by the log-likelihood. The lines denote equal ages (solid), and the depositional age (dashed). The bars show 2 s uncertainties.

To further validate our modeling results, we use the ZR data from the same grains dated with the $\mathrm{U}/\mathrm{Pb}$ and ZHe methods. Figure 4a and b plot the measured ZR ages from the ER and ν₃ Raman bands against the $\mathrm{U}/\mathrm{Pb}$ ages. The ZR ages spread across the depositional age, indicating that these dating systems were partially reset. The lower proportions of post-depositional ages in Fig. 4a and b (ER: 55 %; ν₃: 19 %) compared to the ZHe ages (73 %) agree with the ZR systems' lower susceptibilities to reset and the overlap of the ZHe PRZ with the ZR partial-annealing zones (ER: 180 to 450 °C; ν₃: 220 to 480 °C for 20 to 80 % annealing; Härtel et al., 2021). In addition, both ZR age distributions are bimodal, with a younger population (ER: 200 to 500 Ma; ν₃: 300 to 800 Ma) of grains that crystallized after 1.8 Ga, and an older population (ER: 700 to 1300 Ma, ν₃: 1.1 to 1.7 Ga) of grains with $\mathrm{U}/\mathrm{Pb}$ ages exceeding 1.8 Ga. We determined the pre-depositional model ZR ages based on the distributed-activation-energy kinetic model of Härtel et al. (2021; see Appendix A). Figure 4c and d plot these ages against the $\mathrm{U}/\mathrm{Pb}$ ages. In general, there is negligible age variation between T–t paths, which is explained by the low sensitivity of lattice damage to temperature differences at low degrees of damage annealing. Of the 29 pre-depositional ER ages, 24 (∼83 %) overlap with the expected age range; only three are lower than the depositional age and do not overlap within uncertainties, while two ages pre-date the $\mathrm{U}/\mathrm{Pb}$ ages. For the ν₃ band, all pre-depositional ages overlap with the expected age range within uncertainty. For both ZR ages, the pre-depositional ages preserve their bimodal distributions with the younger group ranging from the depositional age to 800 Ma (ER) and 900 Ma (ν₃), and the older group spanning from 800 to 1500 Ma (ER) and from 1.2 to 1.8 Ga (ν₃). The overall good fits of the pre-depositional model ZR ages to the expected age range confirm the consistency of the modeled T–t paths based on He-diffusion kinetics with the measured ZR ages and the radiation-damage-annealing kinetics.

While a geological interpretation of the T–t paths and pre-depositional ages for this Devonian sandstone is outside the scope of this paper, we compare our modeling results to existing thermal-history models from the Imperial Formation. Figure 3e plots our best-fitting T–t paths (black lines) together with the envelopes of good-fit T–t paths determined by Enkelmann et al. (2019a; blue lines; apatite (U-Th-Sm) $/$ He (AHe) and ZHe data), Powell et al. (2020; grey-shaded area; AHe data with T–t constraints from forward modeling of ZHe data and inverse modeling of apatite fission-track (AFT) data and stratigraphy), and Spalding et al. (2023; red lines; AFT data) – these thermal-history models used a total of five samples within 60 km from our sample location. In general, our best-fitting paths fall into the comparison envelopes, but they show three important deviations: (1) for the Carboniferous–Early Cretaceous part of the thermal history, our best-fitting paths mostly fall into the comparison envelopes, but show less temperature variation, likely due to the $\mathrm{U}/\mathrm{Pb}$ and depositional ages not allowing heating into the ZHe PRZ before the Triassic (Fig. 3b). (2) Some of our T–t paths indicate higher temperatures during the Early Cretaceous, however, this is likely due to different constraints used for each of the models. (3) Our temperatures for the Late Cretaceous–Paleogene heating event are higher than for the models of Powell et al. (2020) and Spalding et al. (2023), but only two of our paths exceed the temperatures suggested by the models of Enkelmann et al. (2019a). This is likely the result of Powell et al. (2020) and Spalding et al. (2023) mostly relying on lower-temperature thermochronometers, which are difficult to reconcile with some of the strongly reset ZHe data. Overall, while our T–t paths match previous modeling results from the Imperial Formation, our modeling approach could benefit further from using additional lower-temperature thermochronometers to better constrain post-reset cooling.

The pre-depositional ages also allow comparison to provenance and source-cooling interpretations of the Imperial Formation. The Imperial Formation is part of the Ellesmerian clastic wedge, representing a mixture of sources from the Canadian Shield and Proterozoic to Paleozoic terranes that collided with Laurentia in the Late-Devonian to Mississippian Ellesmerian orogeny (Hadlari et al., 2009; Beranek et al., 2010; Lemieux et al., 2011). Our pre-depositional ZHe ages (Fig. 3d) form a broad continuous distribution from the Neoproterozoic to the Paleozoic and thus do not contradict Archean to Devonian sources that were partially to fully reset during the Ellesmerian orogeny. The bimodal pre-depositional ZR ages (Fig. 4c and d) separate the cooling of the Canadian Shield in the Paleo- to Mesoproterozoic (older population) from Neoproterozoic to Paleozoic cooling in the younger source terranes (younger population), with a potential minor overprint during the Ellesmerian orogeny. The ZR ages of the younger grain population also overlap strongly with muscovite ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages of Powell and Schneider (2013) from the Imperial formation, which mostly range from 370 to 640 Ma. The similar closure temperatures (muscovite ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$: 350 to 400 °C, Harrison and Zeitler, 2005; ER: 260 to 310 °C, ν₃: 330 to 370 °C, Härtel et al., 2021), suggest comparable thermal histories for the zircon and muscovite source rocks. All in all, our modeling results match the existing thermochronological data from the Imperial Formation in the Mackenzie Plain.

4 Discussion

While the results show the overall viability of our modeling approach for partially reset samples, it also has some limitations. In the following section, we address the limitations and simplifying assumptions of our approach and their practical implications. Then we discuss the application of our modeling approach to other thermochronometers. Finally, we compare the possibilities of partial-reset modeling to strategies avoiding partially reset samples.

4.1 Limitations of the proposed modeling approach

Our modeling strategy successfully leverages the information contained in syn-depositional grains to disentangle pre- and post-depositional thermal histories, but there are four limitations that could hamper its use: first, the availability of syn-depositional grains is the method's Achilles heel. In general, the quality of the post-depositional thermal histories depends on the number and quality of the syn-depositional grains. As a simple heuristic, we suggest aiming for 3–10 syn-depositional grains as typically dated in conventional-ZHe bedrock studies (e.g., Reiners, 2005; Flowers et al., 2022) because the syn-depositional population experienced a single (post-depositional) thermal history similar to bedrock samples. Obtaining enough syn-depositional grains is not critical for strata that were sourced from regions with coeval magmatic activity (e.g. magmatic arcs, rift magmatism), but it complicates matters for strata with source regions where most of the datable zircon grains are eroded from significantly older rocks (e.g., cratons, metamorphic hinterland; Cawood et al., 2012). It is, however, possible to screen the sample for syn-depositional or near syn-depositional grains first. Daniels et al. (2018) applied $\mathrm{U}/\mathrm{Pb}$ analysis by LA-ICP-MS with short ablation pulses (e.g., 5 s) to a large number of detrital grains (n≥600) to find the youngest age population and define a maximum-depositional age. Such a screening approach can be carried out before polishing the grains to the full depth and performing double- or multi-dating. This procedure provides additional $\mathrm{U}/\mathrm{Pb}$ information from grain rims and offers the prospect of recognizing syn-depositional overgrowths on inherited cores from magmatic sources (Malusà et al., 2022). In this regard, laser-ablation-based ZHe dating has a decisive advantage over conventional whole-grain ZHe dating for addressing samples with few syn-depositional grains, as it allows acquiring dozens to hundreds of $\mathrm{U}/\mathrm{Pb}$-ZHe age pairs without biasing the selection toward syn-depositional grains. If no geo- and thermochronometric data from syn-depositional grains can be obtained, it is still possible to use other grains for inverse thermal-history modeling, given that their pre-depositional thermal history is constrained well enough. The most promising way to achieve this – although restricted to coarse clastic deposits – is by analyzing zircon from clasts instead of detrital grains (e.g., Rahl et al., 2011; Grabowski et al., 2013; Falkowski et al., 2016; Stevens Goddard et al., 2024). The main advantage is that other thermochronometers on different minerals or information on the source of the clast can be used together to constrain the pre-depositional thermal history well enough and model the post-depositional history. A different approach is to model the post-depositional thermal history with fully reset lower-temperature thermochronometers such as AFT or AHe dating and using the evaluation of the post-depositional model ages in the next step to weed out the paths not fitting the partially reset ZHe data. As a last resort in the absence of syn-depositional grains, it is also possible to construct simplified post-depositional T–t paths, e.g., in the way of Fosdick et al. (2015) by randomly sampling timing, temperature, and duration of heating for a large number of paths and evaluate their fit with the data with the likelihood criterion in the subsequent step.

Second, while our data show that the modeling strategy performs well at a high degree of partial ZHe reset, it becomes less useful once a sample approaches full reset. In this case, the inverse T–t modeling step will find many good-fit paths that lead to full reset. This increases the number of paths necessary for forward modeling to find solutions that are compatible with the crystallization and depositional ages. In addition, the uncertainty on the pre-depositional model ages increases while reset erases information on the pre-depositional history. This weakens the log-likelihood criterion for model evaluation and hampers the interpretation of source cooling.

Third, the syn-depositional grains may not be representative of the entire grain population in terms of He-retention kinetics. In this case, modeling is not infeasible but may require more inverse T–t modeling to find high-likelihood paths. If the syn-depositional grains are much smaller than the rest of the dated grains, they may have been fully reset after deposition. In this case, the inverse thermal-history modeling software will accept many paths because they reach temperatures just high enough to reset the syn-depositional grains despite being incompatible with the larger grains. Conversely, modeling large syn-depositional grains leads to T–t paths that do not fit with the thermochronological data from the smaller grains. We therefore recommend to select syn-depositional grains with variable size if possible.

Fourth, our modeling example in Sect. 3 relies on several simplifying assumptions that affect model accuracy, but can be dealt with in the future: (1) we assume a uniform U and Th concentration. Uranium and thorium zoning is common in zircon (e.g., Corfu et al., 2003). It determines He-concentration gradients and thus diffusion rates (e.g., Meesters and Dunai, 2002). If information on zoning is available, e.g., from U-Th mapping, it can easily be incorporated into the inverse- and forward-modeling routines by calculating the He accumulation profile based on the U and Th distributions. (2) We did not take into account effects of inherited radiation-damage. In our case study, this decision was supported by low radiation-damage densities directly measured by Raman spectroscopy. For samples containing grains with potentially high degrees of inherited radiation damage, it is possible to incorporate the damage-diffusivity relationship in the forward-modeling step (e.g., Guenthner et al., 2013; Guenthner, 2021). To avoid assumptions about the individual pre-depositional damage accumulation, we suggest using Raman spectroscopic measurements of present-day radiation damage and tracing it back through the post-depositional thermal history (Appendix A) to determine the He-diffusion parameters at every time step. In a subsequent step, it is then also possible to test the extent of radiation-damage annealing in the source regions (e.g., Resentini et al., 2026). (3) For modeling the retention of pre-depositional He during post-depositional heating, we assumed that the initial He-concentration profile within the grain was similar to an α-ejection profile. While this is not necessarily the case, we regard it as a reasonable estimate: if diffusion rounding of the He concentration profile occurred before deposition, it is at least in part diluted toward an α-ejection profile during pre-depositional cooling. Hence, the true He-concentration profile will only differ significantly from an α-ejection profile for a sample that first accumulated a large amount of He prior to deposition, was reheated to temperatures within the PRZ, and cooled immediately before deposition. Nevertheless, it is possible to adjust the pre-depositional He profile in the modeling procedure, e.g., to account for the effect of abrasion during transport. For abraded grains, the He-accumulation profile at the time of deposition falls between two endmembers (see Rahl et al., 2003): either the grains were abraded a long time before deposition and developed an α-ejection profile with respect to its present-day grain boundary, or their complete α-ejection profile was abraded very close to the time of deposition. It is therefore possible to determine a pre-depositional-age envelope for the abraded grains by using two different He-concentration profile, one with and one without α-ejection rim. (4) We assumed approximate ablation-pit positions within the grain for the laser-ablation ZHe data in Sect. 3 and placed them directly in the spherical transformations of our grains (Appendix C). Accurate measurements of the pit positions relative to the 3D grain boundaries would allow the translation of the pit volume to spherical coordinates. Sampling the modeled He concentrations from this translated volume would improve the accuracy of both, the inverse- and forward-modeling steps.

4.2 Other thermochronometers

The focus of our modeling approach is on zircon $\mathrm{U}/\mathrm{Pb}$–(U-Th) $/$ He double-dating, but the overall strategy can also be applied to other thermochronometers. For samples partially reset with respect to the AHe system, the combination with apatite $\mathrm{U}/\mathrm{Pb}$ dating can be applied in the same way. While the apatite $\mathrm{U}/\mathrm{Pb}$ system is more susceptible to Pb loss by thermally activated diffusion than the zircon $\mathrm{U}/\mathrm{Pb}$ system, its temperature-sensitivity window (closure temperature ≈350 to 570 °C; Cherniak, 2010) is much higher than that of the AHe PRZ (40 to 80 °C; Farley, 2000). Hence, the apatite $\mathrm{U}/\mathrm{Pb}$ age allows the identification of grains that formed or cooled from sufficiently high temperatures processes shortly before deposition. Compared to zircon, apatite is mechanically and chemically less stable (Morton, 1982; Jess et al., 2023), making it more useful for tracing the provenance and cooling histories of first-cycle sediment sources. In addition, its chemical composition is more diagnostic for source-rock lithologies (e.g., O'Sullivan et al., 2020).

Another double-dating scheme, for which partial-reset modeling holds great potential, is combined $\mathrm{U}/\mathrm{Pb}$ and fission-track dating, particularly in the light of the laser-ablation-based fission-track dating technique (e.g., Hasebe et al., 2004). Fission-track dating is popular for detrital thermochronology due to the time- and cost-efficient data acquisition for large samples and its compatibility with $\mathrm{U}/\mathrm{Pb}$ dating (e.g., Garver et al., 1999; Bernet and Garver, 2005; Chew and Donelick, 2012; Enkelmann et al., 2019b). Besides the ages, confined-track-length distributions in apatite are a powerful tool for thermal-history modeling (e.g., Green, 1988; Duddy et al., 1988; Ketcham, 2005). Still, for partially reset samples, the uncertainty about each grain's pre-depositional thermal history complicates the interpretation of AFT data (e.g., Fillon et al., 2013). A suggested solution to this problem uses the relative age of confined tracks with respect to their length for distinguishing pre- and post-depositional tracks and modeling the post-depositional thermal history (e.g., Keil et al., 1987; Jensen and Hansen, 2021; Jensen et al., 1992). However, this approach requires a large-enough number of post-depositional fission tracks, which is a potential problem for post-Paleozoic sedimentary samples or samples containing relatively U-poor apatite grains. In contrast, translating our strategy to $\mathrm{U}/\mathrm{Pb}$–fission-track double-dating would allow constraining the post-depositional thermal history from apatite or zircon fission-track ages of the syn-depositional grains. For AFT samples with large numbers of confined tracks, the track lengths of the syn-depositional grains could even provide a high-resolution post-depositional thermal history.

5 Partial-reset modeling vs. unreset or fully reset samples

The above discussion raises a question: (when) is modeling of partially reset detrital samples necessary? This depends first and foremost on the geological question to be addressed, e.g., if one aims to determine the cooling history of source rocks, the sediment-burial depth, or both. To determine cooling ages of source regions, relying on an undisturbed, higher-temperature thermochronometer is a viable strategy. However, this risks opening a large temporal gap (lag-time) between source-region cooling and deposition, potentially overlooking relevant lower-temperature processes. On the other hand, a fully reset low-temperature sample tracks cooling in the basin but only provides a minimum estimate for the burial temperature. Avoiding partial reset also requires additional temperature information, e.g., from stratigraphy, thermal maturation of the analyzed strata, or metamorphic grades in the source regions.

In contrast, the possibility of extracting thermal histories by partial-reset modeling allows analyzing samples with different degrees of reset, e.g., in a stratigraphic succession, across unconformities or faults, without changing chronometers. In addition, a partially reset thermochronometer is most sensitive to heating and cooling in the temperature window the sample reached after deposition. It thus offers the unique opportunity of determining the T–t history of post-depositional heating in more detail than a fully reset chronometer. It also records source cooling through relatively low temperatures, narrowing the gap between initial cooling and deposition. Still, incorporating additional temperature information from thermochronological, thermal-maturation, or sedimentological data is beneficial for partial-reset modeling, but not absolutely necessary to reconstruct the post-depositional thermal history. Finally, while mainly intended for interpreting the thermal history of partially reset samples, our modeling approach is still useful for unreset samples. Calculating pre-depositional ages based on T–t paths from lower-temperature thermochronometers allows testing the upper temperature limit of post-depositional heating: too high temperatures would move the pre-depositional model ages beyond their $\mathrm{U}/\mathrm{Pb}$ age.

6 PaRACAS – a jupyter notebook for partial-reset modeling

Together with the partial-reset-modeling strategy presented in the preceding sections, we introduce the software tool PaRACAS (Partial Reset Analysis for Cooling Ages of Sources). It is available at https://doi.org/10.5281/zenodo.18248224 (Härtel, 2026) as a jupyter notebook together with the underlying Python code and an input-data template with the synthetic data from Sect. 2.

PaRACAS takes as inputs inverse-modeling results in a HeFTy-output text file (.txt) and a comma-separated text file (.csv) with the ZHe and $\mathrm{U}/\mathrm{Pb}$ data, their uncertainties, the $\mathrm{Th}/\mathrm{U}$ ratios and the grain size. After the user specifies a selection of best-fitting T–t paths and a depositional age, PaRACAS carries out the forward modeling and calculates pre-depositional ages and likelihoods as described in Sect. 2 and Appendices A and B. It then plots the results as in Figs. 3 and 4. For ZHe data acquired by laser-ablation in-situ dating, PaRACAS allows adding He-extraction-pit sizes and their approximate positions to improve the match between measured ages and He-diffusion profiles and to determine an abrasion parameter for more accurate inverse modeling in HeFTy (see Sect. 3 and Appendix C). For samples without syn-depositional grains, PaRACAS can generate post-depositional T–t paths based on input ranges for the onset, duration, end and maximum temperature of a single heating pulse (see Sect. 4.1 and Fosdick et al., 2015). Besides ZHe data, PaRACAS also handles ZR thermochronological data, either using the selected paths from ZHe forward modeling as the input for determining pre-depositional ZR model ages as in Sect. 3, or modeling the ZR data on their own.

Plans for further development of the PaRACAS workflow include forward modeling of data from apatite $\mathrm{U}/\mathrm{Pb}$ dating combined with AFT or AHe dating (Sect. 4.2) and the incorporation of U-Th zoning in the diffusion models (Sect. 4.1). The jupyter notebook format and additional Python scripts also facilitate user-specific customization of the workflow.

7 Conclusions

In this paper, we describe a new modeling strategy for interpreting partially reset ZHe data from sedimentary samples using combined zircon $\mathrm{U}/\mathrm{Pb}$–(U-Th) $/$ He dating on individual grains. Its main improvement over existing approaches is that it does not assume the grains' pre-depositional thermal histories. It uses the $\mathrm{U}/\mathrm{Pb}$ ages to find grains that crystallized at or near the time of deposition and determines the post-depositional thermal history from inverse T–t modeling based on these grains. It then proceeds by forward modeling the post-depositional history and determining pre-depositional model ages for each grain. Finally, the likelihood of the pre-depositional model ages falling into the expected interval between $\mathrm{U}/\mathrm{Pb}$ and depositional age is a criterion to identify pairs of plausible T–t paths and pre-depositional ages. Applying this strategy to combined $\mathrm{U}/\mathrm{Pb}$ and ZHe data from a Devonian sandstone from the Northern Canadian Cordillera gives results that are internally consistent and agree with existing thermochronological data. Our approach is further validated by data from the higher-temperature and thus less-reset ZR dating system and pre-depositional ZR ages match existing muscovite ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ data from the Imperial Formation. While leveraging the thermal-history information from syn-depositional grains allows us to separate the pre- and post-depositional histories, it also limits our modeling approach to samples containing such grains. On the other hand, compared to other modeling approaches, our strategy does not require specific assumptions on the pre-depositional thermal history. In addition, calculating pre-depositional ages on a grain-by-grain basis allows to compare them to other single-grain data relevant for sediment provenance. Our approach uses a simple He-diffusion model, but other factors such as U and Th zoning, radiation damage, or grain abrasion can be integrated into the modeling procedure. While this paper focuses on zircon double-dating by the $\mathrm{U}/\mathrm{Pb}$ and ZHe methods or multi-dating including Raman thermochronology, the general principles of our modeling approach can be directly translated to combined apatite $\mathrm{U}/\mathrm{Pb}$ and AHe or AFT dating.

Appendix A: Model setup for pre-depositional age calculation

Our forward model for estimating the fraction of retention and the He accumulated after deposition used a radial diffusion model in spherical geometry. We first translated the dimensions of the analyzed grains to an equivalent spherical radius, using the equation of Ketcham et al. (2011) for a cylindrical shape and divided the radial profile into a discrete number of nodes. To determine the relative proportion of He produced at each radial distance from the center, we determined the α-ejection profile by Monte-Carlo simulation. The stopping distances were calculated in SRIM 2013 (Ziegler, 2013) based on disintegration energies of Haynes et al. (2014), excluding those energies occurring for <1 % of the parent nuclide's decay events. To estimate the respective contributions from the ²³⁸U, ²³⁵U and ²³²Th decay chains, we used the measured $\mathrm{Th}/\mathrm{U}$ ratios, assuming uniform U and Th concentrations throughout the grains.

The next step was to cut and simplify the post-depositional T–t paths provided by the inverse modeling in HeFTy (v. 2.17; Ketcham, 2005, 2024). First, we interpolated the T–t paths to short (e.g., 2 Myr) isothermal segments between deposition and the present day. The lower limit of the depositional-age constraint from HeFTy was used to exclude pre-depositional cooling from the forward models. We summarized consecutive steps with small temperature differences (e.g., 5 °C) in each T–t path by integrating the temperature-dependent diffusion coefficient D (e.g., Lovera et al., 1989):

$$\begin{array}{}\text{(A1)}& \underset{t_i}{\overset{t_{i+\mathrm{2}}}{\int }}D\left(T\left(t\right)\right)\mathrm{d}t=D_{i,i+\mathrm{1}}\left(t_{i+\mathrm{2}}-t_i\right),\end{array}$$

with t_i being the starting time of the ith time step.

The HeFTy modeling from the syn-depositional grains allows T–t paths leading to fully reset. As such solutions are at odds with the observation of a partially reset sample, we filtered out all fully reset paths based on an estimate for the total degree of post-depositional reset. We based this estimate on D_total each path's T–t conditions integrated into a single time step using Eq. (A1). The remaining pre-depositional He concentration depends on the normalized step length F:

$$\begin{array}{}\text{(A2)}& F={\displaystyle \frac{\sqrt{D_{\mathrm{total}}{t}_{\mathrm{total}}}}{R}},\end{array}$$

where t_total is the total time of the thermal history and R is the equivalent spherical grain radius. The analytical equation for the ratio of He retained after the thermal history, ϕ_post, is (Jost, 1952; Fechtig et al., 1961):

$$\begin{array}{}\text{(A3)}& {\mathit{\varphi }}_{\mathrm{post}}={\displaystyle \frac{\mathrm{6}}{{\mathit{\pi }}^{\mathrm{2}}}}\sum _{n=\mathrm{1}}^{\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}}{n^{\mathrm{2}}}}{e}^{-n^{\mathrm{2}}{\mathit{\pi }}^{\mathrm{2}}{F}^{\mathrm{2}}}.\end{array}$$

Figure A1 shows the decrease of ϕ_post with increasing F. We identified the paths with a median F of 0.75 (ϕ_post<0.3 %) among the analyzed grains as full-reset paths and set all their pre-depositional ages to zero.

Figure A1 Relationship between the post-depositional He retention ϕ_post and the normalized step length F.

We determined the degree of post-depositional He retention ϕ_post by forward modeling He diffusion using t_total and D_total. The He concentration was determined at every node in the radial profile using the Crank–Nicolson algorithm (Crank and Nicolson, 1947) with a Neumann boundary at the center (constant concentration gradient) and a Dirichlet boundary at the grain rim (infinite He sink).

We used the kinetic parameters of Reiners et al. (2004) for non-metamict zircon, however, it is possible to use other models such as the zircon-radiation-damage-and-annealing model (ZRDAAM, Guenthner et al., 2013, 2015; Guenthner, 2021). We calculated the fraction of post-depositional retention ϕ_post as the ratio of the remaining concentration C_loss only to the starting concentration C_initial. The whole-grain concentrations C are the radially weighted sums of the concentration C_r at the radial distance r from the grain center:

$$\begin{array}{}\text{(A4)}& C=\sum _{r=\mathrm{0}}^{r=R}{w}_{\mathrm{r}}{C}_{\mathrm{r}},\end{array}$$

The weight w_r at radius r is:

$$\begin{array}{}\text{(A5)}& w_{\mathrm{r}}=\mathrm{4}\mathit{\pi }{r}^{\mathrm{2}}\mathrm{d}r.\end{array}$$

dr is the radial distance between the He-concentration nodes.

We modeled the amount of He accumulated during the post-depositional thermal history with He production and diffusion at each step of the simplified T–t history. The He production used the relative concentrations of the α-ejection profile scaled by the He accumulation rate depending on the changes in U and Th concentrations over geologic time. We determined the pre-depositional He concentration from the observed He concentration and the modeling results using Eq. (3) and calculated the pre-depositional model age from Eq. (4) by iteration of the effective uranium concentration (eU; Härtel and Enkelmann, 2024; their Appendix A) based on the U and Th concentrations at the time of deposition. For α-ejection correction, we scaled the final He-concentration profile against the α-ejection profile. To approximate the model-age uncertainty, we divided the observed ZHe age uncertainty by ϕ_post. We tested the validity of our diffusion model against the analytical equation for volume diffusion from a spherical grain (Eq. A3) and the accumulation and age-calculation procedures against synthetic accumulation-only thermal histories.

For the ZR ages in Sect. 3, we first interpolated the HeFTy paths to isothermal segments as for the ZHe model. We applied an equivalent-time-based modeling approach (Duddy et al., 1988) similar to that of Guenthner (2021) for radiation-damage forward modeling. The principle of equivalent time assumes that future daughter loss only depends on the future T–t conditions and the current degree of reset, but not the specific past T–t history. This implies that any segment of a given thermal history can be translated to a different T–t pair that causes an equivalent degree of reset, which facilitates dealing with complicated thermal histories. Let the post-depositional history comprise k isothermal steps with temperature T_i and duration Δt_i. f_i is the cumulative degree of fractional annealing a portion of daughters produced in the ith step for all steps ≥i. The reset duration Δτ_res,i is then the time necessary to reach f_i by isothermal holding at T_i. f_i follows from the reset function R(t,T):

$$\begin{array}{}\text{(A6)}& f_i=R\left(\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{res},i},T_i\right).\end{array}$$

Δτ_res,i is easily defined for the last time step k, where it simply equals the duration Δt_k:

$$\begin{array}{}\text{(A7)}& \mathrm{\Delta }{\mathit{\tau }}_{\mathrm{res},k}=\mathrm{\Delta }{t}_k.\end{array}$$

The radiation-damage produced at the second-last step k−1 underwent the combined thermal history of steps k−1 and k. Its reset duration $\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{res},k-\mathrm{1}}$ follows from the duration of step k−1 and the equivalent time Δτ_eq:

$$\begin{array}{}\text{(A8)}& \mathrm{\Delta }{\mathit{\tau }}_{\mathrm{res},k-\mathrm{1}}=\mathrm{\Delta }{t}_{k-\mathrm{1}}+\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{eq}}\left(\mathrm{\Delta }{t}_k,T_k,T_{k-\mathrm{1}}\right).\end{array}$$

Δτ_eq is the time needed to reach the degree of reset of the kth step R(Δt_k,T_k) at temperature T_k−1:

$$\begin{array}{}\text{(A9)}& R\left(\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{eq}}\left(\mathrm{\Delta }{t}_k,T_k,T_{k-\mathrm{1}}\right),T_{k-\mathrm{1}}\right)=R\left(\mathrm{\Delta }{t}_k,T_k\right).\end{array}$$

f_i then follows from Eq. (A9) and Δτ_res,i of each previous step:

$$\begin{array}{}\text{(A10)}& f=\left[R\left(\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{res},\mathrm{1}},T_{\mathrm{1}}\right);\mathrm{\dots };R\left(\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{res},k},T_k\right)\right].\end{array}$$

We determined the post-depositional fraction of damage retention equals from the degree of annealing of the first time step. The post-depositional damage density D_postdep is the scalar product of the duration vector τ, the vector f of the fractions of reset at each time step, and ρ, the vector of damage-production rates at each time step:

$$\begin{array}{}\text{(A11)}& D_{\mathrm{postdep}}=\mathit{\tau }\cdot \mathit{f}\cdot \mathit{\rho }.\end{array}$$

We calculated the pre-depositional model ZR ages from Eq. (4). For the radiation-damage modeling, we assumed that most of the lattice damage present at the time of deposition behaved as unannealed, i.e. that no significant annealing took place immediately before deposition, analogous to the assumption of an α-ejection profile (Sect. 4.1). We used the damage-annealing kinetics following the distributed-activation-energy model of Härtel et al. (2021) but other radiation-damage annealing models (e.g., Geisler et al., 2001; Ginster et al., 2019) can also be used, as long as the activation energy of the annealing process is not a function of time (Duddy et al., 1988; Rufino et al., 2023).

Appendix B: Calculation of Likelihoods for each tested T–t path

To evaluate the T–t paths determined by HeFTy, we calculated the path likelihood L_i for the ith grain as the probability of the pre-depositional model age t_p,i falling between the $\mathrm{U}/\mathrm{Pb}$ age $t_{\mathrm{U}/\mathrm{Pb},i}$ and the depositional age t_dep:

$$\begin{array}{}\text{(B1)}& L_i=Pr\left(t_{\mathrm{dep}}<t_{p,i}<t_{\mathrm{U}/\mathrm{Pb},i}\right).\end{array}$$

This likelihood is approximately the integral of a Gaussian function with the model age as mean and the model-age uncertainty (1s) as standard deviation in the limits of the $\mathrm{U}/\mathrm{Pb}$ and depositional age (grey-shaded area in Fig. B1).

Figure B1 Probability-density plot illustrating the assignment of likelihoods to model paths on three grains with Gaussian bell curves describing the pre-depositional model ages (dotted lines) and uncertainties, and the dashed lines representing the $\mathrm{U}/\mathrm{Pb}$ and depositional ages. Shaded areas: model-age density in the expected age envelope.

To include uncertainties on the $\mathrm{U}/\mathrm{Pb}$ and depositional ages, we determined the likelihoods in a Monte Carlo approach. We sampled n values (here 1 000 000) for t_p,i, $t_{\mathrm{U}/\mathrm{Pb},i}$, and t_dep from normal distributions with ages as means and 1s uncertainties as standard deviations and calculated the likelihood by dividing the number of cases, in which t_p,i falls between the other two by n. For grain-thermal-history combinations, in which none of the iterations matched the expected envelope, we set the likelihood to that of half a match. This procedure for calculating the individual-grain likelihoods has several consequences: (1) it penalizes models producing pre-depositional ages outside the expected age envelope. (2) The single-grain likelihoods are relatively constant within the expected envelope, so that no pre-depositional cooling age within the allowed range is preferred over another. (3) The discrete lower limit of the single-grain likelihoods assigns a constant likelihood to model ages far from the expected age interval. Therefore, they lose influence on the model selection once they reach a certain distance from this envelope. We then determined the total log-likelihood LL of each T–t path as the sum of the logarithmized individual likelihoods.

Appendix C: Adapting the zircon (U-Th) $/$ He diffusion model to laser-ablation-based dating

To adapt our modeling strategy to the laser-ablation-based in-situ zircon (U-Th) $/$ He (ZHe) data in Sect. 3, we modified the inverse- and forward-modeling algorithms to only evaluate the part of the He-concentration profile that corresponds most closely to the position of the laser-ablation pit, from which He was extracted. In the absence of precise measurements, we determined pit positions using the following criteria: (1) we assumed spherical grains; (2) we assumed that the pits were centered on the centers of the horizontal polished grain sections; (3) for grains that were drilled through during the U and Th analysis by LA-ICP-MS, we determined the vertical position of the He-extraction pit from the distance of its floor to the bottom of the grain (see Härtel et al., 2026); (4) for the other grains, we estimated the vertical position of the pit based on the amount of material polished away during sample preparation (20 µm; Härtel et al., 2026). Based on the pit position, the grain size and the average pit shape and size (a truncated cone with a radius of 8 µm at the polished surface, pit walls at an angle of 22° to vertical, and a depth of 6.6 µm), we determined the volume distribution of the He-extraction pit in every grain as a function of the distance from the grain center by Monte-Carlo simulation.

We used this information for the inverse- and forward-modeling steps of our analysis. For the inverse modeling in HeFTy, we estimated an “abraded” parameter based on the volume distributions of the pits in the syn-depositional grains. This option in HeFTy assumes a spherical post-abrasion core. We therefore determined the size of a sphere surrounding the grain center, whose volume distribution relative to the grain center is as close to that of the ablation spot as possible. Figure C1a illustrates the volume distributions of an ablation pit (blue line; geometry as above, 10 µm distance from the grain bottom) for a grain with a spherical radius of 30 µm and an abraded core with 25 µm radius (orange line). We evaluated the misfit between the two distributions as the sum of the absolute areas between the curves (grey shaded areas). This value varies between 0 for a perfect match and 2 for two distributions that do not overlap. Figure C1b is a plot of the misfit of the ablation pit in Fig. C1a against different abraded depths for a grain with 30 µm radius. The global minimum of misfit is at about 10 µm and rapidly increases toward lower and higher values. We used the abraded portions for the minimum misfits as “abraded” values as inputs in HeFTy. For the forward modeling, we determined ϕ_post and He_post by modeling each grain and then sampling He from the concentration profiles based on the volume distributions of the He-extraction pits. We also used these volume distributions to correct the pre-depositional He concentrations for α-ejection.

Figure C1 Comparison of the volume-weighted radial positions of He-extraction pits from laser ablation and spherical cores of abraded grains. (a) Radial-volume-weight curves for the laser-ablation pit (blue) in a 30 µm radius zircon grain and the abraded core after 5 µm abrasion (orange). The shaded areas mark radial regions of misfit. (b) Plot of misfit between radial curves as in (a) as a function of the abraded portion.