Simulating sedimentary burial cycles – Part 2: Elemental-based multikinetic apatite fission-track interpretation and modelling techniques illustrated using examples from northern Yukon

Compositionally dependent apatite fission track (AFT) annealing is a common but underappreciated cause for AFT age 10 dispersion in detrital AFT sedimentary samples. We present an interpretation and modelling strategy for samples with variable apatite composition that exploits multikinetic AFT annealing to obtain thermal histories that can provide more detail and better resolution compared to conventional methods. We illustrate our method using a Permian and a Devonian sample from thenorthern Yukon, Canada, both with complicated geological histories and long residence times in the AFT partial annealing zone. Effective Cl values (eCl; converted from r mr0 values), derived from detailed apatite elemental data, are used 15 to define AFT statistical kinetic populations with significantly different total annealing temperatures (~110–185245 °C) and ages that agree closely with the results of age mixture modelling. These AFT populations are well-resolved using eCl values but exhibit significant overlap with respect to the conventional parameters, Cl content or D par . Elemental analyses and measured D par for Phanerozoic samples from the Yukon and Northwest Territories confirm that D par has low precision and that Cl content alone cannot account for the compositional and associated kinetic variability observed in natural samples. An 20 inverse multikinetic AFT model, AFTINV, is used to obtain thermal history information by simultaneously modelling multiple kinetic populations as distinct thermochronometers with different temperature sensitivities. A nondirected Monte Carlo scheme generates a set of statistically acceptable solutions at the 0.05 significance level and then these solutions are updated to the 0.5 level using a controlled random search (CRS) learning algorithm. The smoother, closer-fitting CRS solutions allow for a more consistent assessment of the eCl values and thermal history styles that are needed to satisfy the


Introduction
Apatite fission track (AFT) thermochronology is a well-established method for constraining low-temperature (< 150 ⁰C) thermal histories for a broad range of rock types in a wide variety of geological settings (e.g., Gallagher et al., 1998;Gleadow et al., 2002;Lisker et al., 2009;Malusà and Fitzgerald, 2019;Wagner and Van den Haute, 1992).Fission tracks (FTs) are linear damage zones within apatite crystals that form continuously through time by the spontaneous fission decay of 238 U.An AFT 'age' is a function of the measured track density, which in turn depends on the uranium concentration and the orientation and length of the observed tracks.AFTs form with an initial length of ~16 m but undergo temperaturedependent length reduction (thermal annealing), yielding a track length distribution that reflects the style and rate of heating and cooling of the apatite-bearing sample (Gleadow et al., 1983).For typical fluorapatite, FTs show variable partial annealing between ~20 ⁰C and ~110 ⁰C (Donelick et al., 1990;Green et al., 1989;Ketcham et al., 1999;Tamer and Ketcham, 2020) and are totally annealed at higher temperatures.AFT age and length data can be used to constrain the time-temperature history of a sample over a broad temperature range using a thermal model with appropriate annealing kinetics (e.g., Gallagher, 1995Gallagher, , 2012;;Green et al., 1989;Issler et al., 2005;Ketcham, 2005;Ketcham et al., 2000Ketcham et al., , 2007;;Willett, 1997).
In general, application of the AFT method can be straightforward for crystalline rocks with apatite grains of similar composition and a common formation age.In contrast, sedimentary rocks can pose special challenges for the interpretation and modelling of AFT data because they may contain mixed populations of detrital apatite of variable composition from source areas with significantly different thermal histories.This can lead to discordant AFT grain-age distributions where the variance of the grain ages is greater than expected for the analytical error and therefore the grains cannot be treated as a single age population (Brandon, 2002;Galbraith and Laslett, 1993;Vermeesch, 2019).Overdispersed AFT grain ages are quite common for sedimentary rocks and may result from mixed provenance (e.g., Carter and Gallagher, 2004;Coutand et al., 2006;Garver et al., 1999) and/or variable composition-dependent annealing of fission tracks (Barbarand et al., 2003;Carlson et al., 1999;Crowley et al., 1991;Green et al., 1986;Ketcham et al., 1999;Ravenhurst et al., 2003).The sedimentary samples used in this study have discordant AFT grain-age distributions that we attribute to differential thermal annealing related to apatite composition.
Chlorine content (e.g., Green, 1995;Green et al., 1985Green et al., , 1986) ) and Dpar, the mean length of fission track etch figures on the polished mineral surface parallel to the crystallographic c-axis (Burtner et al., 1994;Donelick, 1993), are two widely used parameters for constraining AFT annealing kinetics (e.g., Barbarand et al., 2003;Donelick et al., 2005;Ketcham et al., 1999Ketcham et al., , 2007)).Unfortunately, these single-parameter methods may only work for simple cases, which may be one of the reasons why there are relatively few published multikinetic AFT studies.The empirical rmr0 parameter (Carlson et al., 1999;Ketcham et al., 1999Ketcham et al., , 2007) ) is potentially more useful for characterizing multikinetic AFT annealing behaviour because it can account for the effects of variable elemental composition, but it has largely been overlooked in the scientific literature.Individual rmr0 values can be calculated for each analysed grain using a multivariate equation with elemental data obtained from electron probe microanalysis (Carlson et al., 1999;Ketcham et al., 2007).The rmr0 parameter was calibrated based on the results of AFT annealing experiments for a range of apatite compositions and it is a measure of the relative resistance to annealing for a given apatite compared with the most retentive apatite used in the experiments (Ketcham et al., 1999(Ketcham et al., , 2007)).It has been applied successfully to constrain multikinetic AFT thermal histories for a number of areas in northern Canada (Issler, 2011;Issler and Grist, 2008a, 2008b, 2014;Powell et al., 2018Powell et al., , 2020;;Schneider and Issler, 2019) where the conventional parameters, Dpar and Cl content, have failed.
In a previous paper (McDannell and Issler, 2021), we used synthetic data to show in principle that it is possible to recover a thermal history with multiple heating events from a single high-quality multikinetic AFT sample using inverse modelling techniques.The obvious implication, discussed in McDannell and Issler (2021), is that multikinetic AFT interpretation can provide more thermal history information than conventional approaches.The purpose of this paper is to demonstrate how compositionally variable detrital apatite grains within natural samples can be grouped into separate statistical kinetic populations that behave as distinct thermochronometers using elemental data and other supporting information.Our method allows for virtually all the sample AFT data to be used; only grains with obviously poor quality AFT or chemical analyses or, more rarely, single grains with unique chemistry or a different provenance may need to be excluded in some cases.The AFT kinetic populations of this study have a wider range of total annealing temperatures (~110 °C to 18245 °C based on model thermal histories) than typical fluorapatite (~110 °C) and therefore they can resolve details of the thermal history beyond the range of a single fluorapatite population.We use an inverse thermal-history model (AFTINV; Issler, 1996;Issler et al., 2005;Schneider and Issler, 2019) to show that it is possible to extract records of multiple heating events preserved in natural multikinetic AFT samples with long residence times in the partial annealing zone.We conclude that age dispersion is a desirable feature of multikinetic AFT samples that can be exploited to resolve thermal histories in unprecedented detail compared with conventional approaches.

Multikinetic AFT methodology
We use two samples from northern Yukon, Canada, to illustrate how we interpret and model multikinetic AFT data.Table 1 summarizes basic sample location and stratigraphic information.and a geological map with plotted sample locations is available in Issler et al. (2021).Percent vitrinite reflectance (%Ro) data (Issler et al., 2021) indicate that both samples experienced paleotemperatures that were high enough (~135 -175 °C) to cause substantial AFT annealing.The geological setting and implications of model results for these samples will be discussed in more detail elsewhere as part of a larger regional study of northern Yukon

AFT and elemental data acquisition
This section discusses the type of data required for multikinetic AFT thermochronology; more details on sample analysis are in Issler et al. (2021).Our AFT data were acquired using the LA-ICP-MS method (Chew and Donelick, 2012;Cogné et al., 2020;Donelick et al., 2005;Hasebe et al., 2004) although the technique works equally well using the older external detector method (EDM; Hurford and Green, 1982).A key difference is that the U data needed for AFT age determination are acquired for the spontaneous track count area after counting is completed (LA-ICP-MS method) whereas counting of spontaneous and induced (proxy measure for U from sample irradiation) tracks are done at the same time (EDM).The LA-ICP-MS method has some distinct advantages compared with EDM.Sample throughput is significantly greater and costs are lower because there is no need to irradiate the sample with thermal neutrons to determine U concentration by proxythus avoiding the need for cool down waiting periods and the counting of induced tracks in mica detectors.Instead, U is measured directly via laser analysis within spontaneous track count areas.The LA-ICP-MS method is well suited for the routine collection of AFT data from a large number of apatite grains because it involves less counting and faster processing than EDM.Typically, 40 single grain AFT ages and 100-200 track lengths are obtained per sample, depending on apatite yield.
Generally, this amount of data is sufficient for most multikinetic samples with two or three kinetic populations, but more data may be required for samples with unevenly distributed populations or with more than three populations.In contrast, many EDM AFT studies have used a lower number of age grains per sample (usually ≤ 20) for thermal history studies.The greater number of counted age grains naturally increases the statistical probability of  2 failure that may complicate mixture model interpretation (McDannell, 2020;Vermeesch, 2019).
It has been suggested that the LA-ICP-MS method may be less reproducible in certain circumstances than EDM because it tends to yield AFT data with larger age dispersion and more  2 failures (Ketcham et al., 2018;Vermeesch, 2017).This has been more recently attributed to U zonation (i.e., Cogné and Gallagher, 2021) but we do not see evidence for this in our samples.Unlike many of the examples shown in Cogné and Gallagher (2021) where grain count areas were large and multiple U spots were measured in high, low, and zero Ns regions, our samples typically have U ablation spots located within smaller count areas to mitigate potentially unrecognizable zoning.We think that both the LA-ICP-MS and EDM methods can yield equally good results for the age grains that are measured, which has been repeatedly demonstrated in published results when EDM and ICP-MS methods were compared or ICP-MS AFT dates were referenced to samples with well-determined absolute ages (e.g., Ansberque et al., 2021;Cogné et al., 2020;Hasebe et al., 2004;Iwano et al., 2019;Seiler et al., 2013;Soares et al., 2014).Arguably, the separate collection of spontaneous count data and U concentration make the LA-ICP-MS method more objective than EDM because there is less opportunity for observer bias to influence results (Donelick et al., 2005).In the former method, spontaneous track areas are counted prior to obtaining U data so there is no way to assess relative age differences among grains during counting.For EDM, there are two issues to consider.First, for many studies, the number of age grains counted is generally much less than 40 (usually ≤ 20).The greater number of counted age grains naturally increases the statistical probability of  2 failure that may complicate mixture model interpretation (McDannell, 2020;Vermeesch, 2019).Second, extra information is available from induced tracks to influence (consciously or unconsciously) how areas are selected for counting (e.g., O'Sullivan, 2018).If only 20 grains are counted 155 and detailed kinetic parameter data are lacking, it makes sense to avoid grains with strongly divergent spontaneous to induced track ratios.The question becomes, "what is not being counted and how does that factor into sample age dispersion?" Figure 1 illustrates the key steps in our multikinetic workflow.Step 1 summarizes the procedures needed for acquiring AFT and related data using the LA-ICP-MS method.Following standard mineral separation, and grinding, polishing and etching of apatite crystals to expose spontaneous tracks, grain mounts are typically 252 Cf -irradiated to increase the number of confined tracks for length measurement (this may not be necessary for samples with high track densities such as Precambrian samples).ing and etching procedures, Then, spontaneous tracks are counted, Dpar is measured for individual apatite age grains (average of four Dpar measurements where possible), and grain x-y coordinates are recorded so that subsequent measurements can be linked to the age grains (Step 1, Figure 1).The sample is re-etched to reveal horizontal confined tracks and their lengths, angles with respect to the mineral c-axis, and Dpar are measured and x-y coordinates are recorded for the measured grains.Finally, Tthe sample is analysed using LA-ICP-MS to obtain U, Th, Sm, and U-Pb age, and trace element (as an option) data for the AFT age grains, ensuring that the laser spot coincides with the track count area to minimize any potential problems with inhomogeneous U distributions. (Jepson et al., (2021) discuss how U-Pb age, trace element, and AFT data can be used to enhance thermal history interpretations.Normally the sample mount is then irradiated with 252 Cf and re-etched to increase the number of confined tracks for length measurement, but this may not be necessary for samples with high track densities such as Precambrian samples.Track lengths, their angles with respect to the mineral c-axis, and Dpar are measured and x-y coordinates are recorded for the measured grains.As an additional option, U-Pb age, U, Th, Sm, and other trace element and U-Pb age data can be acquired for the length grains as well.Although apatite U-Pb ages and Dpar measurements are low precision, the data can be useful for the qualitative assessment of AFT kinetic populations because it is not always possible to obtain elemental data for every grain.In particular, U-Pb ages can be very useful for distinguishing between detrital and volcanic components in a sample.Routine application of the above procedures means that ages and lengths are measured separately and only a subset of grains may have both age and length measurements.There are significant advantages to obtaining length measurements from all the age grains, but the analysis is more time-consuming and is not done routinely.
Step 2 (Fig. 1) involves the acquisition of detailed elemental data for sample apatite grains having AFT age and length measurements.These data are critical for constraining the annealing parameters that are required to recognize and group the data into different statistical kinetic populations.Currently, wWe recommend that elemental data be acquired using electron probe microanalysis (EPMA) rather than by LA-ICP-MS even if the latter method allows for elemental data to be acquired at the point of age measurement and it ismay be more convenient to integrate in the workflow.We have used both methods and they give accurate information for the elements analysed.However, at present, F content cannot always be determined accurately using LA-ICP-MS and therefore OH content cannot be estimated.Without OH information, track retentivity can be underestimated, causing significant overlap of populations in kinetic parameter space, a result consistent with previous studies that emphasize the influence of OH on kinetic parameters (Ketcham et al., 1999;Powell et al., 2018).Analytical Formatted: Subscript conditions for elemental analysis are summarized in the elemental data files included with the Issler et al. (2021) sample report.For the sake of efficient sample processing, a single setup was used for elemental analysis and time-dependent corrections were used to deal with halogen migration (similar to (Nielsen and Sigurdsson, 1981), with the knowledge that crystals oriented with their c-axis parallel to the electron beam could yield some inaccurate results.
In our experience, the following suite of elements are useful for identifying kinetic populations for most multikinetic samples: Fe, Mn, Mg, Na, Sr, La, Y, Ce, F, Cl, Ca, P, Si and S. Fe is very important because current kinetic model calibrations suggest it has a stronger influence on annealing than other elements and it is common in apatite recovered from Phanerozoic rocks of northern Canada.In addition, multikinetic detrital samples can have significant concentrations of elements such as Mn, Mg, Na and Sr. Anions (F, Cl, OH) are known to have a strong effect on annealing behaviour and both F and Cl are essential for estimating OH (e.g., Barbarand et al., 2003;Carlson et al., 1999).Ca and P are useful for assessing the quality of probe analyses and ensuring stoichiometric calculations are accurate.Finally, Si and S have been observed in high abundance for some samples although their contributions to annealing cannot be properly accounted for with existing models.These data may still be useful for assessing relative track retentivity among grains, particularly Si since there is experimental evidence to suggest it greatly enhances track retentivity (in wt % quantities) with respect to the Durango apatite age standard (Tello et al., 2006).The raw analysis time for our samples (Peter Hooper Geoanalytical Lab, WSU) is approximately 3.7 minutes per spot (16 grains per hour), excluding the time for setup (1-2 hours), standardization (8 hours) and point picking (150 grains per hour).Hourly billing is capped at 12 hours per day so it is advantageous to run samples in large batches for 48-72 hours straight which can reduce the hourly rate by 30-40 %.EPMA increases the average cost of AFT analysis by approximately 20% for our samples.
We try to obtain EPMA measurements on a smooth "clean" surface but this is not always possible for small grains with many etched tracks and other imperfections.Missing elements and track void space in the electron beam excitation volume can result in elemental totals that are less than the 97-100 wt% expected for good analyses (see Issler et al. (2021) data tables) and this can happen for samples with AFT data acquired by LA-ICP-MS or EDM.Fortunately, we have observed that suboptimal analyses still yield an elemental signature that allows for discrimination of different kinetic populations and replicate elemental analyses with good and lower elemental totals yield similar results for the study samples.EPMA is undertaken in two passes per grain mount using the x-y coordinates to identify grains with age and length measurements (step 2, Fig. 1).As a result, there may be sets of replicate elemental analyses corresponding to grains with both types of measurements.Replicate analyses (both Dpar and elemental) are very useful for assessing the reproducibility and relative precision of conventional and elemental-based kinetic parameters (see below).Results indicate that compositional zoning is not a common problem for the samples we have studied which is important because the grains are not probed at the exact point of age measurement (laser ablation precedes EMPA).This could contribute to the occasional compositional outlier in kinetic parameter space if elemental zoning is present but kinetic populations are still better resolved compared to when conventional parameters Dpar or Cl are used (see below).Changing the order of steps so that EPMA is done before laser ablation may help alleviate this problem and reduce the number of AFT grains without elemental data but it could delay sample analysis time by up to several weeks because samples must be transferred back and forth between labs and schedules need to be coordinated.Our current method is efficient and works well for the majority of our samples but it can be modified as needed to deal with more problematic samples.For example, if compositional zoning is a significant issue, it might be better to do laser ablation after EPMA and obtain track lengths from age grains only.
Elemental weight % oxide values are converted to atoms per formula unit (apfu) based on an ideal apatite formula (step 3, Fig. 1).We used in house software (Probecal) that incorporates the stoichiometric model of Ketcham (2015) to calculate apfu values (including estimation of OH content) for the Permian sample (Table 1).Elemental data were collected in 2011 for the Devonian sample and provided as apfu values directly from the fission-track laboratory.The apfu values are needed for calculating the empirical annealing kinetic parameter, rmr0 (Carlson et al., 1999;Ketcham et al., 1999Ketcham et al., , 2007)).

Calculation of rmr0, eCl and eDpar values
Composition-dependent AFT annealing was investigated using laboratory experiments (Carlson et al., 1999) and the observed difference in annealing behaviour between different compositional groups was approximated using the equation (Ketcham et al., 1999): where rlr and rmr are the reduced fission-track lengths corresponding to apatite that is less resistant and more resistant to thermal annealing, respectively, and k and rmr0 are fitted parameters.The most resistant apatite in the experimental data set is B2 apatite from Bamble, Norway (highly enriched in Cl and OH) and it is the reference apatite used in equation ( 1).The parameter, rmr0, is the reduced length of the more resistant apatite at the point in the thermal history when the less resistant apatite is totally annealed.From a computational standpoint, this means that annealing calculations are only required for the most resistant reference apatite and the degree of thermal annealing can be determined for any number of less resistant apatite populations using equation (1).Elemental apfu values can be used with an empirical multivariate equation (equation ( 6) in Carlson et al., 1999) to calculate rmr0 values for each analysed apatite grain per sample (step 4, Fig. 1).Ketcham et al. (2007) proposed an alternate equation for rmr0 that was derived from the analysis of the combined annealing experimental Formatted: Subscript data of Carlson et al. (1999) and Barbarand et al. (2003).Either equation can be used but we prefer the original Carlson et al. (1999) equation because it resolves kinetic populations better (less grain overlap) for the samples we have studied.Small values of rmr0 represent highly track-retentive apatite (rmr0 = 0 for B2 apatite) whereas retentivity decreases with increasing rmr0.Although endmember fluorapatite has a nominal rmr0 value of 0.84 for the Ketcham et al. (1999) annealing model, we have encountered less retentive apatite with higher rmr0 values and the model can be extrapolated beyond this limit.
Conversely, the model will not work properly in the unlikely event that an apatite is more retentive than the reference B2 apatite.
For most AFT studies, detailed elemental data are unavailable and Dpar or Cl content are used as kinetic parameters that are converted to rmr0 values in thermal-history models for the annealing calculations.Here, we do the opposite and convert element-based rmr0 values for each apatite grain into "effective Cl" (eCl) values (apfu) (e.g., Issler et al., 2018;McDannell et al., 2019;Schneider and Issler, 2019;McDannell and Issler, 2021) (step 4, fig. 1) using an equation that relates rmr0 to measured Cl (given in figure 7 of Ketcham et al., 1999): where Cl* = Abs(Cl-1).Similarly, we can convert rmr0 values to "effective Dpar" (eDpar) values using the Ketcham et al.
(1999) expression that relates rmr0 to Dpar: Equations ( 2 to yield an equivalent rmr0 value for the Ketcham et al. (1999) annealing model.The data transformation is temporary because eCl values revert back to rmr0 values when used for thermal history modelling.

Relative precision of rmr0 versus conventional kinetic parameters
In order to investigate the relative precision of elemental-based rmr0 values with respect to Dpar and Cl content, elemental data were compiled for fifty-two Phanerozoic AFT samples collected from northern Yukon and the Northwest Territories of northern Canada and then calculated rmr0 values were converted to eCl and eDpar values (step 4, Fig. 1) using equations ( 2) and (3).Replicate elemental (a single EPMA spot per AFT analysis) and Dpar analyses from separate measurements on grains with both age and length data (step 2, Fig. 1) are very important for assessing the reproducibility of kinetic parameter values (Fig. 2).The eCl and eDpar values represent the actual Cl and Dpar measurements that would be required to produce the same rmr0 value as derived from elemental data.Significant differences between these calculated and measured values imply that there are incompatibilities between rmr0 and these conventional kinetic parameters.The data used for figure 2 are in the Supplement; these data will be published in more detail elsewhere after interpretation and modelling have been completed for these samples.Approximately 95 % of the eCl values are within ± 0.06 apfu and a small number of grains show differences as high as 0.24 apfu, likely representing chemical zoning within a sample.A similar plot of 322 replicate results for F is included in the supplement with Table S1 and it shows that 92% of the measurements are within ± 0.2 apfu with larger variations associated with zoning and nonstoichiometric F values.Overall, the accuracy of the halogen measurements is sufficient for estimating OH contents and calculating rmr0 values.These results are encouraging because they suggest that chemical zoning is not a widespread problem and that single-spot probe analyses may be adequate to generate consistent eCl values for most apatite grains.In contrast, replicate Dpar analyses show much larger variation that represents uncertainties of up to ± 20-50 % of the typical measured range in a sample (Fig. 2c).Such large variations in Dpar could make it difficult to resolve different kinetic populations if they are present.Despite the relatively high precision of eCl and Cl measurements, there are strong systematic differences between them (Fig. 2d).The eCl values are skewed to the right of the 1:1 line, indicating higher track retentivity than predicted by Cl alone due to the contributions of OH and elevated cation concentrations.Apatite with low Cl values can still have high eCl values as a result.Dpar is a function of apatite solubility, which in turn is related to mineral composition.
The influence of apatite composition on Dpar may be the reason why data points are more evenly distributed on the Dpar Formatted: Subscript versus eDpar plot (fig.2e).Large differences between eDpar and Dpar seem to be explained mainly by the imprecision of the Dpar measurements (fig.2c).
In summary, the elemental-based rmr0 parameter has significantly higher precision than Dpar and, unlike the single parameter Cl, it incorporates experimentally determined contributions of multiple elements to AFT annealing behaviour.Analysis of Phanerozoic rocks in northern Canada indicates that heterogeneous apatite composition is widespread in sedimentary rocks.
Therefore, rmr0 is superior to Dpar or Cl for recognizing and characterizing multikinetic AFT populations in samples with variable apatite composition.

Kinetic population interpretation
Preliminary interpretation of the data gathered in steps 1 to 4 (Fig. 1) involves visual display of single grain AFT ages on a radial plot (Galbraith, 1990) and age mixture modelling (Galbraith and Green, 1990;Sambridge and Compston, 1994) 3 summarizes the range of measured elemental concentrations for the apatite grains with age and length measurements (see Issler et al. (2021) for elemental data) and it is clear that there are elevated cation and OH concentrations that can increase track retentivity.In all cases, the average Cl value for each population is less than the average eCl values calculated using the Ketcham et al. (1999) or Ketcham et al. (2007) rmr0-Cl equations (equivalent rmr0 values are shown also).The last two columns in Table 2 show the eCl values (and equivalent rmr0 values) that were used for thermal history modelling (see below).2.  The Permian well cuttings sample illustrates some of the issues that can arise when dealing with natural multikinetic samples.It had modest apatite recovery and was processed in two aliquots to obtain 64 single grain ages and 94 track lengths (Fig. 4).Sample drilling contamination does not appear to be a problem; the sandstone interval is overlain by a shale section and shallower Cretaceous AFT samples from the well bear no resemblance to it.Not all of the grains could be probed  4a.A single grain age in each of population two and three crosses the boundary dividing the populations but values are within the expected error range (Fig. 2a).Two relatively high precision ages (one with a single track length; Fig. 4b) plot within population two space but are assigned to population one based on their age.We have more confidence in the age (rather than eCl) because U concentration was determined at the point where the tracks were counted (see section 2.1).Although elemental measurements can be made with high precision, they are the last step after the grains have been subject to multiple treatments (e.g., etching, laser ablation, Cf-irradiation and etching) and it may not be possible to find a clean spot for EPMA that is close to where the age was measured.Therefore, if compositional zoning is present, the grain may appear as an rmr0 outlier.Other lower precision young ages were kept with population two because they had little effect on the population age even if they increase the age dispersion (Fig. 4a).In our experience, the LA-ICP-MS method produces consistent and reliable AFT ages and therefore we use as much of the data as possible.If a higher precision AFT age matches an existing age population but plots as an outlier in kinetic space, we prefer to reassign it to the matching population rather than omitting the grain.
The Devonian outcrop sample is of high quality with 39 single grain ages and 202 track lengths that clearly define two robust kinetic populations in eCl-space (Fig. 5a, b); two single-grain ages from population two are inferred to overlap with population one but this is within measurement uncertainty (Fig. 2a).Both populations overlap completely with respect to Cl (Fig. 5c, d) and Dpar (Fig. 5e, f), indicating that multikinetic behaviour would be undetected or incorrectly interpreted for this sample using conventional kinetic parameters.Well-resolved populations with low age dispersion (Fig. 5a, b), close agreement between radial plot and kinetic population ages (Fig. 3 and Table 2), high thermal maturity, and abundant AFT data make this an excellent and unambiguous example for illustrating compositionally controlled multikinetic AFT annealing.The next step is to obtain thermal history information from the multikinetic samples using inverse modelling techniques (step 8, Fig. 1).

Thermal history modelling of multikinetic AFT data
We use a newer version of the inverse thermal model AFTINV v. 6.1 (Issler, 1996;Issler et al., 2005) to illustrate how detailed thermal history information can be extracted from multikinetic AFT data.We refer the reader to McDannell and Issler (2021) for a recent application of the software to a synthetic AFT dataset.AFTINV uses the Ketcham et al. (1999) multikinetic annealing formulation and can model up to four kinetic populations.Similar to the HeFTy model (Ketcham, 2005), thermal histories are generated randomly and fits between calculated and observed AFT ages and lengths are assessed using p-values.Temperatures are calculated at fixed, user-specified time points using randomly selected heating and cooling rates.Different thermal history styles (Tstyles) can be combined to create complicated thermal histories with multiple phases of heating and cooling.For example, Tstyles 4 (random heating, then cooling), 5 (random cooling, then heating) and 10 (cooling then two cycles of heating and cooling with randomly selected thermal minima) are used for modelling the two samples in this study (Table 4).Monte Carlo calculations terminate when a specific number of solutions (usually 300) exceed the 0.05 level of significance.Unlike HeFTy, a controlled random search (CRS; Price, 1977;Willett, 1997) learning algorithm is then used to try and improve the initial 0.05 level solution set to a higher significance level (typically 0.5 level).
Calculations proceed iteratively until either all members of the solution set exceed the new significance threshold or no further improvements can be made.
Table 4 shows  4) were estimated based on maximum temperatures inferred from organic maturity data, the general degree of annealing among the different kinetic populations, and regional maturity and stratigraphic trends.Dynamic limits are given for thermal minimum (column 6, Table 4) and thermal maximum (column 7, Table 4) inflection points and they are set equal to the static temperature limits where geological constraints are lacking.Heating and cooling rate limits are estimated based on regional geological information and model sensitivity tests.Vitrinite reflectance (%Ro) values are calculated for the entire post-depositional thermal history and for the last phase of heating and cooling using the basin%Ro model (Nielsen et al., 2017) which provides better fits to observed maturity profiles in northern Canada than the (Sweeney and Burnham, (1990) EASY%Ro model (e.g., (Issler et al., 2016;Powell et al., 2020).Further details on model boundary conditions and supporting geological data are the subject of a future paper that deals with the regional thermal history and its geological implications.Model sensitivity runs were undertaken to determine the style of thermal history and the suite of kinetic parameters that were needed to obtain model solutions that closely fit the AFT data.Table 2 lists the model eCl values that yield a set of 300 solutions at the 0.5 significance level.For the Permian sample, the average eCl values of 0.05 and 0.24 apfu could be used for kinetic populations one and two, respectively, but the eCl value for population three had to be increased from the average of 0.37 to 0.55 apfu to obtain solutions that closely fit all three populations.For the Devonian sample, the average eCl value (0.03 apfu) was used for population one but the model eCl value had to be increased substantially from the average value of 0.20 to 0.50 apfu for population two.A three-cycle heating model was attempted for the Permian sample but it cannot be constrained by the data because the first thermal peak had a lower temperature than the second peak.Models using only two cycles of heating (Paleozoic-Mesozoic burial/erosion and Late Cretaceous-Cenozoic burial/erosion) were unsuccessful in finding solutions for the Devonian sample.Thin remnants (20 -30 m) of Upper Triassic (Carnian-Norian) strata occur in isolated places in the northern part of the study region (Norris, 1981) suggesting a burial event of unknown significance.The three-cycle history was set up to allow for a reburial event starting in the upper Triassicbut it could occur anywhere between 0 °C and 150 °C (Table 4; i.e., no requirement to cool to near-surface temperatures prior to reburial).If a Triassic-Early Cretaceous burial/exhumation event were not required by the data, then only a minor inflection should appear on the thermal history over this time interval.

AFTINV
Figure 6 and 7 show AFTINV model results for the Permian and Devonian samples, respectively.The upper panels of both figures show acceptable solution space (≥ 0.05 significance level) defined by the light grey Monte Carlo solutions.The dark grey curves represent the "good" 0.5 level solutions obtained from updating the light grey solutions using the CRS algorithm.These results indicate that the different kinetic populations are mutually compatible and can be modelled with the same thermal history.The blue curve is the exponential mean of the 300 0.5 level solutions and it provides an excellentfitting, smoothed thermal history.The green curve is the closest fitting minimum objective function CRS solution.The initial temperature search space (yellow shaded area) was estimated using regional geological information and model sensitivity analyses.It is larger than the acceptable solution space but small enough to limit the time spent interrogating unproductive regions of solution space.The initial temperature limits collapse to the bounds enveloping the acceptable thermal solutions when the model converges.The lower panels in Figure 6  the shortest track ever measured; Ketcham et al., 2000) and provide an uppermost temperature and time limit for track survival.However, very short tracks are rarely observed and maximum temperatures constrained by the AFT data may be significantly lower.and they give an estimate of how much of the thermal history is preserved in each population.
Two cycles of heating and cooling were used to model the post-depositional thermal history of the Permian sample based on preserved stratigraphy in the well (Fig. 6).A broad range of solutions is permitted at the 0.05 significance level due to the generally low number of track lengths for this sample.Most of the lengths are in population two and therefore it has the strongest influence on the thermal history.The CRS algorithm defines a much narrower region of "good" solution space at the 0.5 level, but it allows for different modes.Most of the peaks for the first heating cycle occur between 160 Ma and 180 heating cycle.The 70 Ma peaks are generated by higher Cretaceous heating rates that are permitted by the sparse length data for population one.Model retention ages were used to estimate total annealing temperatures of 115 °C, 160 °C and 245 °C for population one, two and three, respectively (Fig. 6).Retention ages for the exponential mean solution are plotted on the thermal history for each kinetic population (coloured stars).Modelling suggests that population three has retained tracks > 2 m from the initial phase of cooling at very high temperatures after approximately 540 Ma whereas population two records pre-depositional cooling at lower temperatures after 420 Ma.Any record of pre-depositional cooling is lost for population one due to thermal resetting by the first heating event; it records cooling after 140 Ma and therefore provides no constraint on the magnitude and time of thermal peak one but it is sensitive to the second heating event.Population 2 is most sensitive to thermal peak one and it provides some constraint on the second thermal event.Joint modelling of all kinetic populations improves model resolution.The dashed line (upper panel, Fig. 6) coincides with a change to a steady cooling rate (~1.2 °C/My) below 185 °C at ~440 Ma and marks the upper temperature limit that can be resolved from modelling the AFT data.
The most robust result is the predicted > 30 °C difference in burial temperatures between the two heating cycles.The highly retentive population three has little influence on the post-depositional thermal history.Effective Cl values had to be adjusted higher (than calculated) for population three in order to obtain good solutions that improved the model fit to population three AFT data without significantly changing model temperatures.
The Devonian outcrop sample is of much better quality in terms of data abundance and, although it has only two kinetic populations, it constrains a more complicated and well resolved thermal history with three cycles of heating and cooling of decreasing intensity with time (Fig. 7).The first phase of rapid Paleozoic heating is consistent with regional organic maturity and sedimentological evidence for rapid burial.Results show substantial cooling (> 100 °C) prior to reheating starting in the Triassic.This scenario is permitted but not forced by the model which only requires that an inflection point occur somewhere within the initial temperature envelope over the interval, 240-210 Ma.The third heating/cooling cycle is consistent with the widespread occurrence of Cretaceous strata of generally low organic maturity across the region (Link and Bustin, 1989;Reyes et al., 2013).The three-cycle model works so well that the initial set of Monte Carlo solutions found thirteen 0.5 level solutions (irregular dark grey curves) and the CRS algorithm had no problem updating the remaining solutions to the 0.5 level.Although a Triassic to Early Cretaceous burial/exhumation event was not expected for the area, newer results from the Permian sample to the west (Fig. 6) indicate that heating was sufficient at this time to strongly anneal AFT parameters and overprint Paleozoic thermal maturity at that location, providing further evidence that high-quality multikinetic data can preserve information on multiple heating events.Estimated total annealing temperatures for population one and two are 110 °C and 230 °C, respectively.The very old retention ages (generally > 500 Ma; Fig. 7) and high annealing temperature for population two suggest that it is most sensitive to the pre-depositional cooling history and the first high-temperature peak.
Formatted: Font: Symbol been strongly annealed and retains no information from the first heating event and is most sensitive to the second thermal event.The larger uncertainty on the timing of the third lower temperature peak is consistent with the lower level of AFT annealing.Model results suggest that the AFT data can resolve the thermal history below 175 °C after ~480 Ma (dashed line in upper panel, Fig. 7).
The P013-12 and LHA003 samples were modelled assuming a common inherited, pre-depositional history for each kinetic population in order to better resolve the post-depositional thermal history by taking advantage of relative annealing.We cannot determine whether the kinetic populations within each sample have a shared inheritance because this information has been degraded by thermal annealing of the less track-retentive populations.For LHA003, the pre-depositional thermal record for the lower retentivity population one was erased completely by thermal annealing and therefore pre-depositional cooling is only constrained by population two.For P013-12, pre-depositional thermal history has been erased for population one and population two experienced significant post-depositional annealing.Pre-depositional cooling is dominated by population three which can easily overlap with any residual cooling record for population two.Overall, the post-depositional thermal histories are not very sensitive to the pre-depositional cooling for these samples.

QTQt
Model results from the Bayesian QTQt software (Gallagher, 2012) are shown for comparison with model output from AFTINV.General model setup was the same as the approach discussed in McDannell and Issler (2021), except here the same geologic constraints and %Ro data that were used in the AFTINV models were used in QTQt as well.The differences between the QTQt representative output models (i.e., maximum likelihood (ML), maximum posterior (MP), maximum mode (MM), and expected (EX); see Gallagher and Ketcham, 2020) are discussed in McDannell and Issler (2021) and in Figure 8.
We primarily discuss the ML model since it provides the best fit to the observed data.We prevented more complex models from being accepted during simulations unless they provided a better fit to the datatherefore unnecessary complexity was prohibited and the ML and MP models often end up being similar, with LHA003 being an exception.The overall boundary conditions and heating-style assumptions applied to the AFTINV inversions were absent for QTQt modelling; since the latter relies on the data to directly inform the level of model complexity (i.e., t-T history style).We fixed the eCl kinetic parameter for the well-determined kinetic population most similar to fluorapatite, which was the best constrained group in the annealing experiments of Carlson et al. (1999).The remaining kinetic populations were allowed to vary within uncertainty and underwent resampling during QTQt inversions.
The history for Permian sample P013-12 exhibits a two-pulse heating history in general agreement with the AFTINV results.
Maximum temperatures of ~150 °C were achieved at ca. 135 Ma, followed by a second heating event to ~110 °C at ca. 60 Ma (Fig. 8).These temperatures are at the upper end of the range defined by the CRS solutions in Figure 6 (132)(133)(134)(135)(136)(137)(138)(139)(140)(141)(142)(143)(144)(145)(146)(147) but the 95% confidence region overlaps with the AFTINV results.The QTQt results for the Devonian Imperial Fm. sample are notable, since preliminary AFTINV test simulations investigated the possibility of a two or three-peak thermal history and the ability of the data to resolve the latter scenario within the regional geologic context of preserved Mesozoic outliers.Other than the depositional age of the Imperial Fm., a large constraint box was placed between 55 ± 55 °C to allow for a potential thermal minimum of unknown magnitude at 107 ± 10 Ma (the approximate depositional age of the overlying Cretaceous rocks).The QTQt results demonstrate that a threepeak history is the more likely scenario that provides the best fit to the AFT data (Fig. 9).The EX and ML models suggest the first thermal maximum is at ca. 300 Ma and ~155-190 °C, the second peak is at ca. 155 Ma and ~100-115 °C, and the third occurs near 60 Ma and 70-80 °C.QTQt predicts a younger time for the first thermal peak than the AFTINV CRS solutions (Fig. 7; ~320-360 Ma) but the temperatures and times for the other peaks show good agreement and both models overlap in the 95% confidence region.AFTINV uses larger eCl values and more closely fits the AFT parameters for the most retentive population in both samples.minima/maxima inflection points and preferred model population kinetic parameters.Model results differ in detail for a number of reasons.Compared with QTQt, AFTINV uses more model points, constructs thermal histories differently, allows for manual fine tuning of kinetic parameters, and generates a much larger set of "acceptable" and "good" solutions.QTQt 565 generally prefers simpler histories and there is a trade-off between the number of time-temperature points and data fit.QTQt converts LA-ICP-MS AFT data to EDM AFT data and uses Ns and Ni count data rather than ages for modelling whereas AFTINV models either EDM or LA-ICP-MS AFT data using central or pooled ages, depending on  2 and age dispersion statistics.This difference is most evident for sample P013-12 where QTQt uses younger observed population ages for model input than AFTINV (compare Figs. 6 and 8).570 Figure 8: QTQt model results for Permian detrital AFT sample P013-12 (Jungle Creek Formation).(a) Time-temperature plot of the histories retained post burn-in, coloured by relative probability, with warmer colours denoting higher probability.Individual models include the Maximum Likelihood (highest likelihood and best fit; red curve), Max.Posterior (preferred Bayesian modelsimplest, balancing fit with complexity; gold curve), Max.Mode (peak of the marginal distribution at 1 m.y.intervals; white curve), Expected (weighted mean of the marginal distribution ± 95% credible interval; black curves).The latter two models are summaries of the posterior distribution and not directly sampled during the inversion.(b) Max.likelihood fit to the track length distribution and predictions for AFT age and MTL for population one using a fixed eCl value of 0.05 apfu.(c) Max.likelihood fit to the track length distribution and predictions for AFT age and MTL for population two, allowing resampling of the calculated eCl value.(d) Max.likelihood fit to the track length distribution and predictions for AFT age and MTL for population three, allowing resampling of the calculated eCl value.General information: Green box is depositional age of 290 ± 5 Ma (15 ± 15 °C).Yellow box is geologic constraint at 108.75 ± 6.25 Ma (20 ± 20 °C).Note: All populations utilize the central AFT age in QTQt, whereas the pooled age is used in AFTINV for samples with < 10% age dispersion.Model setup: 500,000 iterations (burn-in) and 500,000 iterations (post burn-in; shown).Prior model space: 300 ± 300 Ma and 125 ± 125 °C.Modern surface temperature of 25 ± 10 °C and maximum allowed ∂T/∂t of 20 °C/My.Proposal moves were rejected if proposed outside of the prior and more complex models were rejected for equivalent likelihood.See text for further discussion.Ma (55 ± 55 °C).Model setup: 500,000 iterations (burn-in) and 500,000 iterations (post burn-in; shown).Prior model space: 350 ± 350 Ma and 125 ± 125 °C.Modern surface temperature of 2.5 ± 2.5 °C and maximum allowed ∂T/∂t of 20 °C/My.Proposal moves were rejected if proposed outside of the prior and more complex models were rejected for equivalent likelihood.

Discussion
Our multikinetic data interpretation and modelling techniques are designed to improve thermal history resolution by exploiting compositionally controlled AFT annealing in samples with high age dispersion.Although numerous factors can contribute to age dispersion, we conclude that multikinetic annealing is the dominant cause of dispersion for Phanerozoic detrital samples from broad geographic areas of northern Canada based on published (Issler et al., 2018;Powell et al., 2018Powell et al., , 2020;;Schneider and Issler, 2019) and unpublished (e.g., Fig. 2) results.This result should apply to other areas that have experienced similar amounts of burial and exhumation and therefore it is of global significance.The method is not restricted to sedimentary rocks and has been applied to Precambrian basement and Proterozoic metasedimentary rocks as well (McDannell et al., 2019;McDannell et al., in press).Heterogeneous apatite compositions are common for multikinetic detrital AFT samples and therefore kinetic populations are much better resolved using the multi-elemental rmr0 parameter than the conventional kinetic parameters, Cl content or Dpar.Although Cl content can be measured with sufficient accuracy and precision, it ignores the documented effects of cation and OH concentrations on track retentivity (e.g., Barbarand et al., 2003;Ketcham et al., 1999).Dpar is influenced by apatite composition but replicate Dpar analyses show that it has low accuracy and precision (Fig. 2).The dearth of published multi-elemental data for AFT studies suggests that true multikinetic behaviour is underrepresented and underutilized in thermal history analysis.
Well characterized multikinetic samples can contain may yield significantly more information than samples with a single AFT population.Therefore, more data and effort are required to interpret and model multikinetic samples, especially if they come from areas with complicated tectonic histories.Generally, 40 age and 100-200 track length measurements are sufficient for typical multikinetic samples with two or three populations, depending on how the data are distributed among the populations.Samples with more populations are less common and may require additional processing to obtain sufficient data to better resolve each population.We can use most or all of the AFT measurements because our interpretations are constrained by multiple parameters (elemental, U-Pb age and Dpar data).In our experience, multikinetic detrital samples are best interpreted as having discrete kinetic populations that are defined by grouping age and length data using eCl values.The close correspondence between population ages inferred from kinetic parameters and peak ages derived from age mixture modelling supports the discrete model approach and provides compelling evidence that differential AFT annealing is controlling age populations.A discrete population model is the simplest interpretation that is consistent with our data and requires less assumptions than a more continuous model that subdivides the data into finer groups using pre-determined incremental kinetic parameter values.The latter model ignores the results of age mixture modelling and assumes that all measurements represent accurate kinetic parameter values over the full (kinetic) model range.
adjustment is not required because the model resamples kinetic parameters in an attempt to fit the data.In addition to the uncertainty with rmr0 calculations, high retentivity populations can span a broad range of kinetic parameter space and it can be difficult to obtain representative average eCl values for undersampled populations (Fig. 4) or for populations with unevenly distributed data that cluster at one end of the range (Fig. 5).For the Permian sample, average eCl values for kinetic populations one and two are within the range of more typical apatite and require no adjustment whereas the population three eCl value had to be increased to obtain 300 CRS solutions that fit all three populations.For the Devonian sample, the eCl value was fixed for population one and the eCl value for population two was adjusted relative to population one until 300 solutions at the 0.5 significance level were obtained.Preliminary investigations of resampling the kinetic parameter for all populations within QTQt, suggests that the 'fluorapatite' population often remains relatively stationary, whereas 'exotic' endmember kinetic populations are sampled outside of their calculated eCl range.This aligns with results from the annealing experiments and lends support for our AFTINV approach of relative kinetic adjustment during interpretation and modelling.
The allowable range for estimated eCl values for high retentivity populations depends on the degree of annealing.At low levels of annealing, there is lower sensitivity to the thermal history and a broad range of eCl values can yield good solutions.
Numerous models were run for the Permian sample, but results are only presented for the model with an eCl value of 0.55 apfu for kinetic population three (Fig. 6) because it gave the highest number of solutions at the 0.5 significance level.Models with eCl values between 0.45 apfu and 0.55 apfu also gave a significant number of 0.5 solutions.The effect of changing eCl was to improve the fit to AFT parameters while having a negligible effect on the thermal history which is largely constrained by populations one and two.For the Devonian sample, a range of eCl values is permissible as well for the higher retentivity population two but we only show results for an eCl value of 0.5 apfu (Fig. 7) because it gave 300 solutions at the 0.5 level with the broadest temperature envelope.The number of 0.5 level solutions decreased when eCl was increased to 0.55 apfu or decreased to 0.45 apfu.When eCl was increased from 0.5 apfu to 0.55 apfu, the region defined by the 0.5 level solutions narrowed and the corresponding average peak temperature shifted from 173 °C at 341 Ma to 178 °C at 351 Ma with calculated average vitrinite reflectance increasing from 1.55 %Ro to 1.72 %Ro.The variation was much less (1 °C and 1 Ma) when eCl was decreased to 0.45 apfu.We do not consider these uncertainties in eCl values to be a problem for multikinetic modelling because the effects on the thermal history are minor and well within available geological constraints.
In our experience, relative annealing can significantly limit the allowable offset in eCl values (incremental changes < 0.05 apfu) between lower retentivity populations that experienced strong annealing.In some cases, incremental changes must be < 0.05 apfu or solutions cannot be obtained at the 0.5 level.For our Permian sample, good solutions were obtained without adjusting the average eCl values for populations one and two.
. The early Permian (~295-285 Ma) sandstone sample is from well cuttings collected approximately 70 m below a pre-Cretaceous unconformity in the Eagle Plain west of the Richardson Mountains.Approximately 1 km of Cretaceous strata (≤ 0.6 %Ro) overlie the unconformity, which represents approximately 170 million years of missing geological record.The second sandstone sample is from an Upper Devonian (~375-365 Ma) outcrop northeast of the first sample on the western flank of the Richardson Mountains.These samples provide an opportunity to investigate whether multikinetic samples can retain important details of the thermal history in areas with a complicated tectonic history and much of the stratigraphic record missing.

Figure 1 :
Figure 1: Flowchart outlining the key steps involved in the acquisition, interpretation, and modelling of multikinetic AFT data.
) and (3) can be used to transform measured kinetic parameters (i.e., Dpar and Cl) to rmr0 values or vice versa by rearranging the equations in terms of Cl and Dpar (see equation (1) of McDannell and Issler (2021) for eCl).The Ketcham et al. (2007) multikinetic model has similar equations with slightly different coefficients.For example, the rmr0 value for fluorapatite in this model is 0.83, which translates to an eCl value of ~0.03 apfu and an eDpar of ~1.85 µm.The corresponding rmr0, eCl and eDpar values for the Ketcham et al. (1999) model are 0.84, 0.0 apfu and 1.75 µm, respectively.There are several advantages to converting nonlinear rmr0 values to "linear" eCl values for data interpretation purposes.First, eCl values are more evenly distributed on linear x-y plots which enhances the visual display and interpretation of the data.Second, arithmetic averages of eCl values can be used to estimate representative eCl values for different kinetic populations.Third, eCl values can be compared with the conventional parameter, Cl content, to show how other elements enhance AFT retentivity beyond what is expected from Cl alone.The eCl parameter represents the Cl concentration that would be required Formatted: Not Superscript/ Subscript

Figure
Figure 2a and 2b show a comparison of replicate eCl and Cl values, respectively, from single spot elemental analyses on apatite grains from fifty different samples.Both eCl and Cl are reproducible within ± 0.03 apfu for the majority of the data, representing a 5-10 % variation over the range of measured values.As expected, a somewhat higher percentage (91 % versus 84 %) of Cl values show closer agreement than the eCl values which were determined using multiple elements.
used.Age and length data are well resolved for each kinetic population when plotted with respect to eCl (Fig.4a, b and Fig. 5a, b) but show significant to complete overlap when plotted with respect to Cl content (Fig. 4c, d and Fig. 5c, d) and Dpar (Fig. 4e, f and Fig. 5e, f).Population boundaries were set at eCl values of 0.134 and 0.312 apfu for the Permian sample (Fig. 4a, b), and 0.725 apfu for the Devonian sample (Fig. 5 a, b).Percent vitrinite reflectance (%Ro; Table 2) measurements indicate that both samples exceed the threshold for full organic maturity (0.6 %Ro), implying significant AFT annealing.Average Dpar and Cl values are listed for each kinetic population in Table 2.Although average Dpar increases with population age, these values are not viable kinetic parameters due to the severe population overlap in Dpar space.Table

Figure 3 :
Figure 3: Radial plots of single grain ages for the multikinetic Permian (left) and Devonian (right) AFT samples.Points are colourcoded according to the kinetic populations determined in figure 4 and 5. Peak ages are from age mixture modelling (estimated percentage of grains per population in brackets) and show good correspondence with kinetic population ages summarized in Table2.

Figure 4 :
Figure 4: Plots of (a) single grain ages and (b) track lengths grouped into different colour-coded kinetic populations using eCl values for the Permian well sample.Boundary between populations is indicated by vertical black lines.Similar plots of population ages and lengths with respect to Cl concentration (c and d) and Dpar (e and f) using the same colour-coding as in (a) and (b).

Figure 5 :
Figure 5: Plots of (a) single grain ages and (b) track lengths grouped into different coloured-coded kinetic populations using eCl values for the Devonian outcrop sample.Boundary between populations is indicated by the vertical black line.Similar plots of population ages and lengths with respect to Cl concentration (c and d) and Dpar (e and f) with same colour scheme as (a) and (b).
the boundary conditions used to model the Permian (P013-12) and Devonian (2009LHA003) samples.The Permian sample was modelled with an initial pre-depositional cooling event followed by two cycles of random heating and cooling (Tstyle 10).The Devonian sample has an extra cycle of heating and cooling that was accommodated seamlessly by adding an additional random cooling/heating (Tstyle 5) and heating/cooling (Tstyle 4) cycle.Each model starts at a high temperature (245-250 °C) and cools to surface temperatures (0-30 °C) within a time range for deposition that captures uncertainties in the stratigraphic age.The model uses static and dynamic temperature limits.Static limits define the entire model search space whereas dynamic limits are applied only at model inflection points to focus calculations into favourable regions of solution space to improve model efficiency.The static temperature limits (column 5, Table and 7 display the model and c-axis projected measured track length distributions (coloured histograms) and model retention age distributions for each kinetic population.Model track length distributions are shown for the exponential mean (blue) and minimum objective function (green) solutions along with the regions encompassing track length distributions for the 0.05 (light grey) and 0.5 (dark grey) level solutions.Model retention ages represent a theoretical age forof the oldest (shortest) surviving track in each population (assumed to be ~ 2 µm based on

Figure 6 :
Figure 6: Upper panel shows AFTINV thermal history results for the Permian sample.Light grey lines are statistically acceptable Monte Carlo solutions (≥ 0.05 significance level); dark grey lines are good solutions (≥ 0.5 level).The black curves bounding model solutions are not valid solutions.The blue curve is the exponential mean (EM) of the 300 good solutions; the green curve is the closest fitting minimum objective function (MOF) solution.Lower panels show model and observed track length distributions, and the distribution of model retention ages (age of oldest track) for the different kinetic populations.The goodness of fit (GOF) probability for the age and length data is given for the exponential mean solution.Tanneal is the estimated total annealing temperature for each population.Uncertainties on average retention age, average peak temperature and average peak time are two standard deviations.

Figure 7 :
Figure 7: Upper panel shows AFTINV thermal history results for the Devonian sample.Light grey lines are statistically acceptable Monte Carlo solutions (≥ 0.05 significance level); dark grey lines are good solutions (≥ 0.5 level).The black curves bounding model solutions are not valid solutions.The blue curve is the exponential mean (EM) of the 300 good solutions; the green curve is the closest fitting minimum objective function (MOF) solution.Lower panels show model and observed track length distributions, and the distribution of model retention ages (age of oldest track) for the different kinetic populations.The goodness of fit (GOF) probability for the age and length data is given for the exponential mean solution.Tanneal is the estimated total annealing temperature for each population.Uncertainties on average retention age, average peak temperature and average peak time are two standard deviations.

Figure 9 :
Figure 9: QTQt model results for Devonian detrital AFT sample LHA003 (Imperial Formation).(a) Time-temperature plot of the histories retained post burn-in, coloured by relative probability, with warmer colours denoting higher probability.Individual model descriptions are 590

Table 1 .
Sample location information Figure1is a flowchart outlining the key steps for multikinetic AFT analysis and data interpretation.The order of steps is based on (1) efficiency and speed of analysis, (2) maximizing the number of track lengths, (3) minimizing selection bias for age grains, and (4) obtaining replicate elemental data.The method can be modified to optimize for other factors or to deal with particular sample conditions but this may increase the cost or the time for analysis.Steps 1 and 2 involve the acquisition of AFT and apatite elemental data.Steps 3 and 4 concern the processing of data to obtain kinetic parameter values that are associated with AFT age and length measurements.Steps 5 and 6 deal with visual displays of the data to aid in the interpretation of kinetic populations.Step 7 involves assessing the interpretation by considering all available data in the context of measurement uncertainty and missing information.The goal here is to try to use all the available data except for obviously poor analyses that provide no useful information (i.e., inaccurate U measurements).Exclusion of any outlier data should be justified on reasonable analytical or geological grounds.Inverse modelling of the interpreted data is carried out in step 8.The objective is to obtain thermal solutions that provide close fits to the data within acceptable statistical precision while satisfying available geological constraints.In general, this requires iterative modelling to refine implicit model parameters, explicit boundary conditions, and investigate different thermal history styles.

Table 1 .
Issler et al. (2021)ch, 2012)may have more than one age population (step 5; Fig.1).We use the DensityPlotter software(version 8.4;Vermeesch, 2012)to plot single grain AFT ages and estimate age populations.Single grain ages and track length measurements are then plotted with respect to eCl values to determine whether discrete age and associated length populations can be identified as statistical kinetic populations (step 6; Fig.1).Although boundaries separating interpreted kinetic populations are set at a fixed eCl value, overlap of some population age and length data across boundaries is Age mixture modelling yields three and two age populations for the Permian and Devonian samples, respectively.Single grain ages are colour-coded according to the interpreted kinetic populations in figures 4 and 5 and the model peak ages (Fig.3) are within one standard deviation of the kinetic population ages summarized in Table2(seeIssler et al. (2021)for complete AFT data).All kinetic populations pass the  2 test but the complete samples fail with high age dispersion (Fig.3).Pooled ages are used for the kinetic populations if age dispersion is < 10 %; otherwise, the central age is expected, given the documented uncertainties in eCl values (Fig.2a).Other complications may include missing elemental data for some grains, compositional zoning, outlier ages associated with a different provenance, and poor-quality data.Some ambiguity in the interpretation can be reduced by including other data such as apatite U-Pb age data, Dpar values, replicate elemental analyses, and quality control procedures to refine the interpretations (step 7, Fig.1).Kinetic population interpretation is facilitated if all length measurements come from age grains; if not, length data must be sorted using eCl values and information from replicate measurements.Single grain ages can be "coloured-coded" with respect to the kinetic populations on a radial plot for comparison with the radial arms representing model age populations.There is strong evidence to infer multikinetic behaviour if age populations based on eCl values correspond closely with those determined from age mixture modelling, especially if there is other evidence (e.g.organic maturity) to indicate that burial temperatures were hot enough for substantial AFT annealing.Figure3shows radial plot results for the Permian (Jungle Creek Fm.) and Upper Devonian (Imperial Fm.) AFT samples listed in

Table 3 .
Summary elemental data for AFT age and length measurements (atoms per formula unit) approximately 30 % of the age and length data have no associated elemental data) but there is enough information to demonstrate that three populations are present.Data from the unprobed grains were sorted using age, Dpar and information from replicate analyses that link length data to age grains.Missing probe data is not a serious concern for this sample; population ages are similar (within a few million years) if age grains without probe data are excluded.Sixteen of the length measurements without elemental data are associated with age grains that were sorted into various populations.Eleven lengths with intermediate Dpar values were included with the dominant population two; the mean track length changes by < 0.3 m if these data are excluded.Four age grains are interpreted to overlap with other populations in figure

Table 4 .
Model boundary conditions/constraints for two AFT samples