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 age dispersion 10 in detrital AFT samples. We present an interpretation and modelling strategy 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 the Yukon, Canada, both with complicated geological histories and long residence times in the AFT partial annealing zone. Effective Cl values (eCl; converted from rmr0 values), derived from detailed apatite elemental data, are used to define AFT statistical kinetic populations with significantly different total 15 annealing temperatures (~110–245 °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 Dpar. Elemental analyses and measured Dpar for Phanerozoic samples from the Yukon and Northwest Territories confirm that Dpar has low precision and that Cl content alone cannot account for the compositional and associated kinetic variability observed in natural samples. An inverse multikinetic AFT model, AFTINV, is used to obtain thermal history 20 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 AFT data. The high-quality Devonian sample (39 single grain ages and 202 track 25 lengths) has two kinetic populations that require three cycles of heating and cooling (each subsequent event of lower intensity) to obtain close-fitting solutions. The younger and more westerly Permian sample with three kinetic populations only records the latter two heating events. These results are compatible with known stratigraphic and thermal maturity constraints and the QTQt software produces similar results. Model results for these and other samples suggest that elementalderived eCl values are accurate within the range, 0–0.25 apfu (rmr0 values of 0.73–0.84), which encompasses most of the data 30 https://doi.org/10.5194/gchron-2021-22 Preprint. Discussion started: 2 August 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction 35
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 40 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 (Green et al., 1989;Ketcham et al., 1999) 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 45 temperature range using a thermal model with appropriate annealing kinetics (e.g., Gallagher, 1995Gallagher, , 2012Green et al., 1989;Issler et al., 2005;Ketcham, 2005;Ketcham et al., 2000Ketcham et al., , 2007Willett, 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 50 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 55 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 65 et al., 1999Ketcham 65 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 70 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, 2008bPowell et al., 2018Powell et al., , 2020Schneider and Issler, 2019) where the conventional parameters, Dpar and Cl content, have failed.

75
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 80 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 245 °C) than typical fluorapatite (~110 °C) and therefore they can resolve details of the thermal history beyond the range of a single fluorapatite 85 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. https://doi.org/10.5194/gchron-2021-22 Preprint. Discussion started: 2 August 2021 c Author(s) 2021. CC BY 4.0 License.

Multikinetic AFT methodology 90
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. 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.  Figure 1 is a flowchart outlining the key steps for multikinetic AFT analysis and data interpretation. 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 105 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. 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 110 requires iterative modelling to refine implicit model parameters, explicit boundary conditions, and investigate different thermal history styles.

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., 115 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). 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 proxy -thus 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 120 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. 125 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 c 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 130 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;135 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 140 counted age grains naturally increases the statistical probability of c 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 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 145 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 grain mounting and etching procedures, spontaneous tracks are counted, Dpar is measured for individual apatite age grains, and grain x-y coordinates are recorded so that subsequent measurements can be linked to the age grains. The sample is analysed using LA-ICP-MS to obtain U, Th, Sm and U-Pb age data for the AFT age grains, ensuring that the laser spot coincides with the track count area 155 to minimize any potential problems with inhomogeneous U distributions. 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 option, U, Th, Sm 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 160 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. 165 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. We recommend that elemental data be acquired using electron probe microanalysis (EPMA) rather than by LA-ICP-MS even if the latter method may be more convenient to integrate in the 170 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). In our experience, the following suite of elements are useful for identifying kinetic populations 175 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 180 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). 185 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). Elemental weight % oxide values 190 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 stochiometric 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 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. 195

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): 200 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 205 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 210 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. 215 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 220 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) fig. 1) using an equation that relates rmr0 to measured Cl (given in figure 7 of Ketcham et al., 1999): (2) 225 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: 230 Equations (2) and (3)  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 240 retentivity beyond what is expected from Cl alone. eCl represents the Cl concentration that would be required 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 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 245 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 and Dpar analyses from separate measurements on grains with both age and length data (step 2, represent the actual Cl and Dpar measurements that would be required to produce the same rmr0 value as derived from 250 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. These results are encouraging because they suggest that chemical 260 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 ( 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. 275 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) to 280 investigate if the sample 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 285 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 290 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. 295  Table 1. 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 Table 2 (see Issler  300 et al. (2021) for complete AFT data). All kinetic populations pass the c 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 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.  305 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.

310
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 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. 315 (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).  Table 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 325 and shallower Cretaceous AFT samples from the well bear no resemblance to it. Not all of the grains could be probed (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 330 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 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 335 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. 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 340 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).
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 345 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 350 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, 360 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 365 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. 370 Table 4 shows 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 375 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. The static temperature limits (column 5, Table 4) were estimated based on maximum temperatures inferred from organic maturity data, the general degree of annealing among the different kinetic populations, and regional 380 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). Further details on model boundary conditions and 385 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 390 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 395 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 Triassic -but it could occur anywhere 400 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. fitting, 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 and 7 display the model and c-axis projected measured track length 415 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 the age of the oldest surviving track in each population (assumed to be ~ 2 µm; Ketcham et al., 2000) and they

AFTINV
give an estimate of how much of the thermal history is preserved in each population. 420 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 425 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 Ma but two higher temperature peaks at 192 Ma are associated with higher older temperature peaks at 70 Ma in the second 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 430 thermal history for each kinetic population (coloured stars). Modelling suggests that population three has retained tracks from the initial phase of cooling at very high temperatures after approximately 540 Ma whereas population two records predepositional 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 435 thermal peak one and it provides some constraint on the second thermal event. Joint modelling of all kinetic populations improves model resolution. 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. 440 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 445 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 450 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 455 population two suggest that it is most sensitive to the pre-depositional cooling history and the first high-temperature peak.
Most of the retention ages for population one are older than, or similar to, the second thermal peak -indicating that it has 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. 460  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.

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 465 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 data -therefore unnecessary complexity was 470 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 475 and underwent resampling during QTQt inversions.
The history for Permian sample P013-12 exhibits a two-pulse heating history in 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). The QTQt results for the Devonian Imperial Fm. sample are notable, since preliminary AFTINV test 480 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 three-peak history is the more likely scenario that provides the best fit to the AFT data (Fig. 9). The EX and ML models 485 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. The AFTINV and QTQt results are very similar, even with the subtle tradeoffs between the different thermal minima/maxima inflection points and preferred model population kinetic parameters. 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 model -simplest, 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 495 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 500 °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.

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 520 detrital samples from broad geographic areas of northern Canada based on published Powell et al., 2018Powell et al., , 2020Schneider 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 525 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 530 behaviour is underrepresented and underutilized in thermal history analysis.
Well characterized multikinetic samples can contain 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 535 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 540 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 545 measurements represent accurate kinetic parameter values over the full (kinetic) model range.
An important point is that successful modelling of multikinetic samples relies on the relative annealing behaviour that is implicit in the rmr0 model(s) of Ketcham et al. (1999Ketcham et al. ( , 2007. The basic assumption is that the same annealing mechanism applies to all apatite, but that composition controls the temperature at which annealing occurs. This reasonable assumption 550 was used successfully by Ketcham et al. (1999Ketcham et al. ( , 2007 to account for experimental annealing data for apatite of variable composition. We extend this approach to natural multikinetic AFT samples and we are able to show that a common thermal history can reproduce AFT ages and lengths for the different kinetic populations. Ketcham et al. (1999) advised against general use of the Carlson et al. (1999) elemental-based multivariate equation for predicting rmr0 values for unknown apatites due to limited calibration data. We agree that this equation (or the equivalent equation in Ketcham et al., 2007) cannot 555 predict rmr0 values accurately over the full range of apatite compositions that are likely to be encountered in nature, and accurate prediction for poorly represented apatite may well be an unattainable goal. However, this should not be a deterrent to using a method that has the potential to yield better results than conventional approaches, keeping in mind that the conventional parameters are no more accurate and involve the same annealing model. In order to advance the field of thermochronology, it is essential to pursue the logical consequences of the annealing experiments and show that it is possible 560 to use existing techniques to better understand complicated AFT samples.
The examples in this paper, other published studies (Powell et al., 2018(Powell et al., , 2020Schneider and Issler, 2019) and unpublished results (Fig. 2) suggest that kinetic parameters can be accurately determined for eCl values of approximately 0 -0.25 apfu (rmr0 values of 0.73 -0.84), a range that encompasses most of the data for the annealing experiments (Ketcham et al., 1999) 565 and represents more typical apatite compositions. In our experience, the Carlson et al. (1999) multivariate equation overestimates track retentivity for endmember fluorapatite with rmr0 values > 0.84 (requiring more negative eCl values than predicted), and underestimates retentivity for more exotic higher retentivity apatite with rmr0 values < 0.73 (requires higher eCl values than predicted). A successful tactic we have used is to anchor model calculations on the kinetic population having the best calibrated annealing behaviour within the above kinetic parameter range. In this case, the model eCl value can be 570 fixed at the average eCl value for the anchor population and eCl values for atypical populations can be adjusted as required to obtain successful solutions that fit all populations. Issler et al. (2005;their Fig. 17) used this approach to investigate how the objective function value of the exponential mean solution changes as the kinetic parameter of one population is adjusted with respect to a population with a fixed kinetic value. It is possible to determine the optimal kinetic parameter offset that yields the lowest minimum object function value (closest fit). Optimization generally occurs in a kinetic parameter interval 575 where changing the kinetics has little effect on the thermal history solution (i.e., minor changes in annealing sensitivity), yet model misfit is minimized. For software like QTQt, manual 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.  580   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 585 '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 590 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 595 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 600 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 < 605 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.
Our modelling of natural and synthetic (McDannell and Issler, 2021) AFT data demonstrates that multikinetic samples can retain a record of multiple heating and cooling events under suitable geological conditions due to the relative annealing 610 behaviour of kinetic populations that are sensitive to different parts of the thermal history. The ability to recover multicycle histories depends on many factors including sufficient and well distributed AFT and elemental data to constrain kinetic population interpretations, a favourable thermal history (decreasing thermal intensity with time) that caused significant annealing of multiple kinetic populations, and appropriate application of modelling strategies. Even if the first two conditions are met, careful modelling is needed to obtain successful solutions. Nondirected Monte Carlo modelling software 615 such as AFTINV and HeFTy depend on user-defined model boundary conditions and kinetic parameter values -often requiring iterative modelling to refine model parameters to obtain solutions. Overly simplistic thermal histories or unsuitable kinetic parameters may result in no or few solutions and this may be attributed incorrectly to problems with the data. For example, a two-cycle heating/cooling model failed to yield solutions for our high-quality Devonian sample. When boundary conditions were adjusted to allow for the possibility of an extra heating event, the model converged on three well resolved 620 heating cycles, a scenario that is compatible with available geological constraints and is independently supported by QTQt model results without user-imposed boundary conditions. The QTQt model allows for looser boundary conditions and greater uncertainty in the style of thermal history because it implements a reversible jump Markov chain Monte Carlo optimisation algorithm that automatically modifies kinetic parameters (within specified ranges) and the number of heating cycles to try to fit the input data. 625 Clearly not all samples will be multikinetic, and not all multikinetic samples will share a common thermal history if there has been insufficient annealing to eliminate differences in provenance (i.e., inherited pre-depositional histories).
Nevertheless, we have seen many cases where a common thermal history works for multikinetic samples. Possible explanations for this include rapid exhumation of heterogeneous source areas, strong annealing/thermal resetting of some or 630 all AFT populations, and mixing of detrital and syn-depositional volcanic components. If interpreted populations are incompatible, then solutions will not be obtained by adjusting kinetic parameters. Ultimately, the value of the approach will be judged on its ability to generate spatially coherent thermal histories over different stratigraphic intervals across study regions. An encouraging sign is that we have been able to use essentially the same kinetic parameters for samples from the same stratigraphic units in the Yukon that have experienced different degrees of heating. We strongly recommend that 635 elemental data be collected for detrital AFT samples from areas with complicated geological histories and for other samples with unexplained age dispersion. Age dispersion is a desirable characteristic of samples when viewed in a multikinetic framework due to the potential for enhanced thermal-history resolution.

Conclusions
It is common for detrital AFT samples of variable cation and anion composition to have significant age dispersion that is 640 caused by multikinetic annealing. Under these conditions, AFT age and length data can be sorted into discrete kinetic populations with different annealing temperatures using eCl values (derived from rmr0 values obtained using multi-elemental https://doi.org/10.5194/gchron-2021-22 Preprint. Discussion started: 2 August 2021 c Author(s) 2021. CC BY 4.0 License. data). In general, these kinetic populations are unresolved or poorly resolved using the conventional single kinetic parameters, Dpar or Cl content; Dpar has low precision and Cl alone neglects how other elements influence track retentivity.
Modelling of dozens of samples from northern Yukon, two of which are presented here, indicates that a complicated record 645 of multiple heating and cooling cycles can be retained in multikinetic samples under certain geological conditions (heating cycles of decreasing intensity through time) due to the different relative annealing behaviour of the kinetic populations.
Accurate prediction of elemental-derived kinetic parameters is unlikely and may not be attainable for all natural apatite populations, however this is not a requirement for using the method and the same problem exists with conventional kinetic parameters. Absolute kinetic parameters are best constrained for eCl values within the range of 0-0.25 apfu (rmr0 values of 650 0.73-0.84) which represents the more commonly encountered apatite compositions from published annealing experiments.
Effective Cl values in this range can be fixed in the model and inaccurately predicted parameters for more exotic apatite compositions can be adjusted by exploiting the relative annealing inherent in the rmr0 model. As expected, the uncertainty range on the less constrained higher eCl values increases as the degree of AFT annealing decreases. Overall, the model is tolerant of these uncertainties and a range of eCl values for higher retentivity populations can still produce similar solutions 655 at the 0.5 significance level. Considering these results, age dispersion in multikinetic samples should be viewed as desirable for enhancing thermal history resolution rather than as hindrance to data interpretation and modelling.

Supplement
The Supplement contains the data for figure 2.

Author Contributions 660
DRI developed the method, interpreted and modelled the data using AFTINV, and wrote the paper. KTM performed QTQt modelling, was involved in conceptual discussions, and contributed to writing and editing the paper. PBO did the AFT analysis (first at A to Z Inc. and then at GeoSep Services) and coordinated the elemental analyses at Washington State University. LSL collected the samples, arranged for the analyses, and provided geological context for the study.

Competing interests 665
The authors declare that they have no conflict of interest.

Acknowledgements
These ideas were developed over many years of working with detrital AFT samples and benefitted from many influences.
Early on, Ray Donelick impressed upon DRI how widespread age dispersion is in detrital samples and his pioneering work with the LA-ICPMS AFT method enabled us to acquire the data for this study. Richard Ketcham generously helped DRI to 670 better understand his annealing model so it could be incorporated into AFTINV. Careful AFT and multi-elemental analyses by Sandy Grist (formerly of Dalhousie University) for samples from the Mackenzie Delta region provided the first clear evidence for distinct kinetic populations in Canadian samples. We thank Chance Oil and Gas Ltd. (formerly Northern Cross (Yukon) Ltd.) for providing the well sample used in this study. Owen Neill (formerly of the Peter Hooper GeoAnalytical Laboratory, Washington State University) and Ken Severin (Advanced Instrumentation Laboratory, University of Alaska 675 Fairbanks) are recognized for EPMA analyses and their assistance in organizing and processing the data. This study has been supported through the Natural Resources Canada Geo-Mapping for Energy and Minerals Program Yukon Basins Project.