Whitehorse Trough records Late Triassic–Cretaceous accretionary orogenesis in the Northern Canadian Cordillera via detrital mineral thermochronometry

The northern Canadian Cordillera comprises a complex assemblage of arc, oceanic, and microcontinental terranes that have accreted to the western margin of Laurentia (ancestral North America) since the Mesozoic (Fig. 1A; Nelson et al. 2013). The arc and oceanic “Intermontane” terranes that occupy the core of the orogen were emplaced along North America's western passive margin before or during the Early Jurassic (e.g., Murphy et al. 1995; Colpron et al. 2005, 2015; Nixon et al. 2020, 2022; Nelson et al. 2022). Their current configuration has been imaged geophysically as an apparent broad synformal (the Northern Intermontane synform) stack, with the oceanic Cache Creek terrane overlying Stikinia–Quesnellia, which in turn overlies the Yukon–Tanana terrane (Calvert 2016; Calvert et al. 2017). However, interpretations of Intermontane terrane assembly and accretion history vary, ranging from Paleozoic assembly of Cache Creek and Stikinia–Quesnellia arc terranes offshore of Laurentia (Zagorevski et al. 2021) to Late Triassic to Early Jurassic oroclinal bending of Yukon–Tanana and Stikinia–Quesnellia arc terranes with oceanic Cache Creek terrane entrapment (Mihalynuk et al. 1994; Colpron et al. 2015). Recent models involve end-on collision of the Stikinia arc with the Yukon–Tanana terrane, analogous to the modern interaction of the Aleutian arc with Kamchatka (George et al. 2021; Nelson et al. 2022), followed by entrapment of the Cache Creek terrane and development of a sinistral transform fault system linked to the southward retreat of an Early Jurassic arc (Colpron et al. 2022).

Syn-orogenic strata assigned to the Laberge Group were derived from rapidly eroded arc and basement rocks of the Intermontane terranes and shed into the Whitehorse Trough, a >600 km long basin that formed on Stikinia and the Cache Creek terrane in southern Yukon and northern British Columbia during the Early to Middle Jurassic (Colpron et al. 2015; Kellett and Zagorevski 2021, 2022). The Whitehorse Trough developed within a sinistral transtensional setting in the arc-collision zone from early Sinemurian to Toarcian (van Drecht et al. 2022). Middle Jurassic imbrication of the Cache Creek terrane and Stikinia and inversion of the basin into a fold and thrust belt (Figs. 1B and 1C; Wheeler 1961; Hart 1997; English et al. 2005) were followed by lithospheric delamination, surface uplift, and emplacement of 174–168 Ma post-collisional granitoids in Yukon (Colpron et al. 2022). Whitehorse Trough rocks underwent Cretaceous dextral transpression because of ongoing plate convergence along western North America (White et al. 2012; Monger and Gibson 2019), and to the south, the younger Bowser Basin was inverted to form the northeast-vergent Skeena fold and thrust belt in northern British Columbia (Nelson et al. 2013). Thus, the Whitehorse Trough rocks have been sensitive recorders of upper crustal processes during a key period in the tectonic evolution of the orogen. The precise timing of burial, shortening, exhumation, and reburial has not been well defined, but would better reveal the orogenic events contributing to the evolution of the Intermontane portion of the Canadian Cordillera.

Detrital mineral thermochronology of multiple mineral-decay systems with varying closure temperatures provides a unique and powerful dataset in unravelling the development of syn-orogenic sedimentary basins (Fig. 2). Higher temperature thermochronometers that retain their detrital ages provide a means to fingerprint sediment sources, calculate source exhumation rates (via lag times across multiple chronometers), and constrain maximum depositional ages (MDA), including for source rocks that have been otherwise erased from the geological record (Garver et al. 1999; Kellett et al. 2018). During stratigraphic and structural burial of sediments, low temperature thermochronometers may undergo open system behavior due to heating, providing information about minimum rates and peak thermal conditions of burial, potentially to a greater degree of confidence than thermal maturity data such as vitrinite reflectance and RockEval, for which the presence of second cycle organic matter can be difficult to diagnose (e.g., Nzoussi-Mbassani et al. 2005). The subsequent progressive passage through specific partial retention or annealing zones during final cooling directly relates to latest exhumation. Sedimentary basins in orogenic settings may be subjected to shortening (i.e., inversion) during or immediately following deposition, or may experience multiple cycles of shortening/tectonic burial—exhumation—subsidence/renewed deposition, as evidenced by unconformities. Comparison of thermal histories for detrital minerals above and below such unconformities can provide further constraints on these multiple burial cycles. Fault displacements cause differential burial and/or differential tectonic exhumation that may also produce contrasting syn- to post-depositional thermal histories within basin rocks in different structural positions (e.g., Malusá and Fitzgerald 2019 and references therein).

Geological evidence suggests that Whitehorse Trough strata experienced two structurally controlled burial episodes, separated by an intervening angular unconformity. Here we test those geological constraints against a range of detrital mineral geo- and thermochronometers including published zircon U–Pb, rutile U–Pb, muscovite ⁴⁰Ar/³⁹Ar, and new zircon (U–Th)/He (hereafter ZHe), apatite fission track (AFT), and apatite (U–Th)/He (AHe) thermochronometers, sampled from the Laberge Group and representative under- and overlying units along the length of the Whitehorse Trough. Available thermal maturity data provide peak temperature estimates reached in the basin, but do not provide any timing information, making it difficult to distinguish, date, and quantify the timing and magnitude of the two burial/heating episodes. We apply forward and inverse modeling approaches, focused particularly on the ZHe thermochronometer, to constrain timing and rates of each heating/burial episode and relate these to the Mesozoic accretionary orogenic cycles experienced by the Intermontane terranes.

The Whitehorse Trough, first named and characterized by Wheeler (1961), is an elongate, narrow, orogen-parallel basin in southern Yukon and northern British Columbia (BC), Canada, that contains fine- to coarse-grained siliciclastic strata and minor felsic pyroclastic and epiclastic rocks (Fig. 1). The historical designation of “trough” refers to the narrow and linear geometry of the basin. The term “Whitehorse Trough” was originally assigned to carbonate and siliciclastic upper parts of the Triassic Stuhini/Lewes River Group (BC/Yukon local terminology) and Lower to Middle Jurassic Laberge Group strata and intervening volcaniclastic rocks (Fig. 1B; e.g., Lowey 2008). Both the Stuhini/Lewes River Group and Tantalus Formation represent different depositional environments in terms of basin extent and sedimentation type from Laberge Group. They are separated from these Lower to Middle Jurassic strata by local angular unconformities, and we use the term Whitehorse Trough here to refer exclusively to the Laberge Group (White et al. 2012; Calvert et al. 2017; Hutchison 2017; Kellett and Zagorevski 2021; Colpron et al. 2022; van Drecht et al. 2022).

The Laberge Group is underlain by siliciclastic (including Mandanna member, Aksala formation) and carbonate (including Hancock member of the Aksala formation, Sinwa formation) rock units that mark the waning stages of the Stuhini/Lewes River Group, a major Late Triassic volcano–plutonic complex developed on Stikinia (Fig. 1B; English and Johnston 2005; Bordet et al. 2019). Similar siliciclastic rocks deposited on Cache Creek terrane and Quesnellia (Kedahda and Shonektaw formations, respectively) have common origins and represent an overlap sequence across all three terranes, implying Late Triassic continuity between these terranes (Colpron et al. 2015; Zagorevski et al. 2018, 2021). The Late Triassic depositional environment on Stikinia indicates an arc-proximal, marginal marine setting (Shirmohammad et al. 2011; Long 2015). The oldest Laberge Group strata are separated from Stuhini/Lewes River Group by an (locally angular) unconformity (Dickie and Hein 1995; Shirmohammad et al. 2011; Bordet et al. 2019; van Drecht et al. 2022).

Laberge Group strata comprise a >3-4 km-thick succession deposited during the Early to early Middle Jurassic (e.g.,Shirmohammad et al. 2011; White et al. 2012). The Laberge Group has been sub-divided into broadly coeval facies/units: the Tanglefoot/Takwahoni formations (Yukon/BC terminology) consist of terrestrial to marginal marine conglomerate, greywacke, shale and volcaniclastic rocks; the Richthofen/Inklin formations (Yukon/BC terminology) include submarine fans and mass flows with siltstone, sandstone, and conglomerate; and the Nordenskiold facies consist of local Pliensbachian pyroclastic flows (Johannson et al. 1997; Lowey and Long 2006; Lowey 2008; Colpron and Friedman 2008; Hutchison 2017, and references therein). We preferentially collected Tanglefoot/Takwahoni facies medium to coarse-grained sandstone units as these contain abundant detrital zircon and apatite. Coeval Lower Jurassic chert to argillite/greywacke strata overlie Cache Creek terrane (Fig. 1B), and potentially represent a continuous overlap sequence of different facies across Stikinia and Cache Creek (Colpron et al. 2015; Kellett and Zagorevski 2021; Zagorevski et al. 2021).

The abundance of conglomerate and lithic wacke throughout the Laberge Group, including throughout the well-exposed, complete, and fossiliferous stratigraphic section at Lisadele Lake, BC (Fig. 1A; Shirmohammad et al. 2011), shows an evolution of sediment sources during the Early to Middle Jurassic (Dickie and Hein 1995; Hart et al. 1995; Hart 1996; Johannson et al. 1997). Sedimentary and volcanic clasts dominate basal Sinemurian(?) strata at both Lisadele Lake and Atlin Lake, BC, and are likely derived from the underlying Stuhini Group. Clasts in lower Pliensbachian strata are predominantly volcanic and hypabyssal in type, with sedimentary clasts being rare in comparison to underlying strata. The proportion of plutonic clasts steadily increases to ∼60% through upper Pliensbachian to lower Toarcian strata, with a concomitant decrease in the proportion of volcanic clasts. This progression has been interpreted as tracking the evolving tectonic setting of sediment sources from Sinemurian flank uplift (sedimentary clast-dominated) to lower Pliensbachian transitional arc (volcanic and hypabyssal clast-dominated) to upper Pliensbachian arc dissection (intrusive rock clast-dominated) of (probably) the nearby Stuhini/Lewes River and Hazelton groups, ultimately exposing their plutonic roots (Dickie and Hein 1995). More recent work has shown that Early Jurassic detrital zircon in the Whitehorse Trough have evolved Hf isotopic signatures rather indicative of syn-collisional plutons (Colpron et al. 2022; van Drecht et al. 2022).

Detrital zircon U–Pb studies have shown that most Whitehorse Trough strata were sourced from proximal Late Triassic to Early Jurassic igneous rocks, with secondary Paleozoic and rare Precambrian sources (Shirmohammad et al. 2011; Colpron et al. 2015; Kellett et al. 2018; Kellett and Iraheta-Muniz 2019; Bordet et al. 2019; Kellett and Zagorevski 2021; van Drecht et al. 2022). Most Laberge Group samples are characterized by major or unimodal zircon age peaks that are within ∼10 myr of biostratigraphic or MDA in the BC (i.e., southern) portion of the Whitehorse Trough, and within ∼5 myr of MDAs in the Yukon (i.e., northern) portion (Kellett and Zagorevski 2021; van Drecht et al. 2022). In combination with the clast compositions, this relatively short lag time supports a proximal, volcanic to shallow plutonic arc origin for most detrital zircon grains in Laberge Group sandstones (Kellett and Zagorevski 2021).

A discrete interval of metamorphic detritus, including garnet, eclogite, granulite, amphibolite, mica schist, and ultra-high pressure garnet peridotite, occurs in Upper Pliensbachian strata at Eclogite Ridge, northern BC, with abundant detrital garnet in coeval strata at Atlin Lake likely representing a common stratigraphic horizon (Fig. 1A; Canil et al. 2006; Kellett et al. 2018). Detrital rutile, dated in situ in eclogite clasts (2.2–2.9 GPa, ∼800 °C) from Eclogite Ridge, yielded a U–Pb cooling age of 182 ± 15 Ma, indicative of rapid Late Triassic or Early Jurassic exhumation along a subduction channel (Kellett et al. 2018). Upper Toarcian strata at Lisadele Lake contain clasts of quartzofeldspathic schist and gneiss from which detrital white mica yielded an age of ca. 195 Ma (Fig. 1A; Shirmohammad 2006; Kellett and Iraheta-Muniz 2019), whereas Stikinia and Cache Creek terrane lithologies generally lack muscovite; metamorphic muscovite is common in some Yukon–Tanana terrane rocks, with the majority of published bedrock muscovite ⁴⁰Ar/³⁹Ar cooling ages being between 200 and 190 Ma (Breitsprescher et al. 2004; Joyce et al. 2015; Kellett and Zagorevski 2021; Yukon Geological Survey 2022). These thermochronometric data and additional evidence for an Early Jurassic Barrovian-type metamorphic event (e.g., Currie and Parrish 1993; Dusel-Bacon et al. 2002; Berman et al. 2007; Gaidies et al. 2021) suggest a Yukon–Tanana source for metamorphic detritus deposited into the Whitehorse Trough (Canil et al. 2006; Kellett et al. 2018).

Ammonite biostratigraphy locally provides precise depositional age constraints for Laberge Group strata (e.g., Johannson et al. 1997), with the well-exposed Lisadele Lake section ranging from Pliensbachian to Bajocian (ca. 190.8–168.3 Ma; Shirmohammad et al. 2011; Walker et al. 2018). Detrital zircon U–Pb datasets for Laberge Group strata span the length of the Whitehorse Trough (Shirmohammad et al. 2011; Colpron et al. 2015; Kellett et al. 2018; Kellett and Iraheta-Muniz 2019; Bordet et al. 2019; van Drecht et al. 2022; van Straaten et al. 2022). Various methods are available to determine MDAs from detrital zircon U–Pb datasets. Here we use the MDAs reported in Kellett and Zagorevski (2021): Laberge Group strata range from ca. 208–174 Ma (Kellett and Zagorevski 2021 and references therein; van Drecht et al. 2022). Since the youngest underlying Lewes River Group MDA from Kellett and Zagorevski (2021) is ca. 202 Ma, we consider 202 Ma rather than 208 Ma as an upper MDA value for the Laberge Group. These depositional constraints are valuable for incorporating the timing of deposition into thermal history models.

Laberge Group and underlying Stikinia strata were deformed into an SW-verging fold and thrust belt with minimum shortening estimates ranging from 15% to 35%, which included the formation of the King Salmon, Kehlechoa, and Nahlin faults (Fig. 1C; Mihalynuk et al. 1995; Gabrielse 1998; English et al. 2002, 2005; Mihalynuk et al. 2017; White et al. 2012; van Straaten and Gibson 2017; van Straaten and Bichlmaier 2018a). The fold and thrust belt was exhumed and overprinted by Cretaceous dextral strike-slip faults, including the Teslin fault (Figs. 1A and 1C; Gabrielse et al. 2006; Colpron 2011; White et al. 2012).

The exposed top of the Laberge Group is a disconformity upon which Middle Jurassic to Lower Cretaceous siliciclastic strata were locally deposited (Fig. 1B). In Yukon, Bathonian(?) to Lower Cretaceous conglomeratic strata of the Tantalus Formation were deposited in confined, mountain valley-type fluvial and lacustrine beds (Bostock 1936; Colpron et al. 2015; Long 2015) and record surface uplift as a result of lithospheric delamination (Colpron et al. 2022). Sedimentary rock sources for the Tantalus Formation are indicated by an abundance of chert, quartz, silicified mudstone clasts (Long 2015), and recycled detrital zircon grains from the Laberge Group (Colpron et al. 2015). In BC, siliciclastic rocks locally unconformably overlie the Laberge Group in the footwall blocks of the King Salmon and Kehlechoa faults at the SW—most extent of the Whitehorse Trough (e.g., Fig. 1C). These strata have been correlated with the Bowser Lake Group, based on their MDAs, clast compositions, and lithologies (Shirmohammad et al. 2011; van Straaten and Gibson 2017; van Straaten and Bichlmaier 2018a): the Bowser Lake Group is the main unit of the Late Jurassic to mid-Cretaceous Bowser Basin, a major sedimentary basin located SW of the Whitehorse Trough. Bowser Lake Group conglomeratic strata have abundant chert clasts like that of the Tantalus Formation. Chert clasts with early Permian to Triassic fossils in the Bowser Lake Group are generally interpreted to have Cache Creek terrane provenance (Cordey et al. 1991; Cordey 1992; Mihalynuk et al. 2004; Evenchick and Thorkelson 2005; Colpron et al. 2015), but other proximal sources of roughly coeval chert are exposed in the Stikine assemblage (Logan et al. 2000; Mihalynuk et al. 2012), Tsabayhe Group (Read 1984), and Sinwa Formation (Long 2015; Mihalynuk et al. 2017). The Bowser Basin was inverted in the Cretaceous to form the northeast-vergent Skeena fold and thrust belt (e.g., Evenchick 1991).

Laberge Group and under- and overlying strata yield VRo and equivalent values that range from 0.4 to 2.4, corresponding to peak rock temperatures ranging from ∼80 to 200 °C (Nielsen et al. 2017), with few samples yielding values > 3.5 (Fig. 3A). The highest and lowest rock temperatures in Yukon are reported for Laberge Group and Tantalus Formation samples along the western margin of the Whitehorse Trough, with intermediate values in the central and eastern portions of the basin (Fig. 3B). Laberge Group samples from the Atlin Lake area in BC show a basin margin-perpendicular gradient from low thermal maturity in the NE to high thermal maturity in the SW (Fig. 3C). In the context of this thermochronology study, we therefore expect most or all detrital apatite to preserve post-depositional AHe and fission track ages, with closure/annealing temperatures for these systems likely reached or exceeded during basin evolution (e.g., Fig. 2). However, the ZHe closure temperature range overlaps significantly with the range of peak temperatures experienced in the basin, and we expect only local resetting of He in detrital zircon. Regions recording the highest rock temperatures, such as the SW-most edge of the Laberge Group in the Atlin area, are more likely to preserve post-depositional ZHe ages, whereas regions recording the lowest rock temperatures, such as the NW extent of the Tantalus Formation (e.g., in the Tantalus Butte area), are more likely to preserve pre-depositional ZHe ages. There are no thermal maturity data available for significant areas of the Whitehorse Trough footprint.

We sampled a wide range of geographic locations (samples span the full length of the basin), structural positions (across-strike sampling), and stratigraphic levels (including representative samples from underlying and overlying stratigraphic units) that are at least 1 km away from syn- and post-depositional intrusions (Table 1, Fig. 1A). Units that were likely to yield detrital zircon and detrital apatite grains were prioritized. Accordingly, we targeted mostly fine- to medium-grained sandstone and greywacke. Rock samples were subjected to electrical pulse disaggregation and density separations including water table and heavy liquid methods, followed by magnetic separation using a Frantz Magnetic Separator, to isolate apatite- and zircon-rich mineral fractions. Where possible (12 samples), all three thermochronometric methods of ZHe, AFT, and AHe were applied. Together, 43 samples yielded 171 new ZHe single grain ages, 12 new AFT central ages, and 158 new AHe single grain ages.

Apatite and zircon grains were manually selected using a Nikon SMZ1500 stereoscopic microscope at ×106 magnification under plane- and cross-polarized light and photographed and measured (width and length) using NIS Elements v4.10 imaging software. These measurements were used to calculate the α-ejection correction (Farley et al. 1996). Preference was given to euhedral, transparent, inclusion- and crack-free grains with a smallest dimension being not less than 70 µm. All selected grains were packed individually in Nb tubes. For a small subset of samples, we employed a double-dating approach (Reiners et al. 2005): zircon grains that had been mounted, polished down to their mid-sections, and U–Pb dated (Kellett and Iraheta-Muniz 2019) were then plucked from their mounts and dated with our (U–Th)/He methodology, but applying a half α-ejection correction, since half the grains’ volume have been removed by polishing. This approximation for the α-ejection correction necessarily introduces additional uncertainty in the calculated ages compared to analysis of full grains. They should be regarded with caution and not used for thermal history modeling. For those samples, the grain number relates to the dated grain from the previous U–Pb analyses, and the published ²⁰⁶Pb–²³⁸ U ages of those grains are provided in supplementary data. Note that a preliminary summary of some of these data was presented in Kellett and Zagorevski (2022). That presentation did not include methods, any raw or calculated single crystal age results, or the full data tables, all presented here.

He measurements were completed on an in-house built He extraction line equipped with a 45 W diode laser and a Pfeiffer Vacuum Prisma quadrupole mass-spectrometer at Dalhousie University. Each Nb packet was heated with the laser to release He. Apatite grains were heated to 1050 °C for 5 min, whereas zircon grains were heated to 1250 °C for 15 min. After He extraction, a precisely measured aliquot of ³He was added to the sample and the ³He/⁴He ratio was measured using the quadrupole mass spectrometer. This procedure was repeated once for apatite grains to ensure no He was left in the grain. For zircon grains, the He extraction step was repeated up to 3–4 or more times until the amount of He in the last re-extraction was less than 1% of total He extracted from the grain. Typical errors are in range of 1.5%–2% (1σ) for both zircon and apatite. Samples were analyzed in groups of 36, and each group included two Durango apatite standards for apatite analyses (32.0 ± 1 Ma; Farley 2002) or two Fish Canyon tuff zircon standards in the case of zircon analyses (28.3 ± 2.6 Ma; Reiners 2005).

After He extraction, apatite and zircon grains were combined with a mixed ²³⁵U, ²³⁰Th, and ¹⁴⁹Sm spike and dissolved. Apatite grains were dissolved in 7 N HNO₃ at 80 °C for 1.5 h, whereas zircon grains were dissolved in high-pressure dissolution vessels in a mixture of concentrated HF and HNO₃ at 200 °C for 96 h. Isotopic ratios were measured using an iCAP Q ICP-MS at Dalhousie University. Raw data were reduced using the Helios software package developed by R. Kislitsyn and D. Stockli for (U–Th)/He data reduction. The α-ejection correction was calculated from grain surface-to-volume ratios (Farley et al. 1996). Standards were subjected to the same analytical procedures as unknown samples to ensure data accuracy, reproducibility, and reliability. Dalhousie University lab results for Durango apatite and Fish Canyon tuff zircon at the time of the analyses were 31.4 ± 1.81 Ma (1σ, n = 85) and 28.3 ± 2.6 Ma (1σ, n = 81), respectively.

AFT samples were processed at the Thermochronology laboratory at Dalhousie University. Apatite aliquots were mounted in Araldite epoxy on glass slides, ground and polished to expose internal grain surfaces, then etched for 20 s in 5.5 M HNO₃ at 21 °C to reveal spontaneous fission tracks. All mounts were prepared using the external-detector method (Hurford and Green 1983). Samples and CN₅ glass standards were irradiated with thermal neutrons in the Oregon State University reactor. After irradiation, the low-U muscovite detectors that covered apatite grain mounts and glass dosimeter were etched in 40% HF for 45 min at 21 °C to reveal induced fission tracks. Samples were analyzed using a Kinetek computer-controlled stage driven by the FTStage software (Dumitru 1993) attached to a Zeiss Axioplan microscope. Dry counting was done at a magnification of ×1000 and 5 to 21 grains were analyzed per sample. Fission track ages were calculated using a weighted mean Zeta calibration factor (Hurford and Green 1983) based on IUGS ages standards (Durango, Fish Canyon, and Mount Dromedary apatites) (Miller et al. 1985; Hurford 1990). The Zeta calibration factor for the operator (I. Coutand) is 370.6 ± 5.0.

The complete datasets including raw data are available in the supplementary data, and AFT results are shown in Table 2. Single crystal ZHe and AHe ages are shown with AFT central ages in kernel density estimate plots in Fig. 4B, and in along- and across-strike age profiles of the Whitehorse Trough in Fig. 5. Age results are summarized in Table 1 and map view in Fig. 4A as mean (ZHe, AHe) or central (AFT) ages, except ZHe ages are reported in Table 1 as ranges for samples in which ZHe ages are older than MDAs and hence include age data from potentially multiple source rocks. Errors in He ages are considered to be at the level of the long-term reproducibility of the Durango apatite and Fish Canyon zircon standards for the Dalhousie University Noble Gas laboratory at the time the analyses were completed, 6.4% for ZHe and 4.8% for AHe. Double dating results are summarized in Table 1 and Fig. A1. For reporting and displaying mean ZHe and AHe ages, samples with outlier individual grain ages have been identified as having intra-sample dispersion. The outlier(s) were removed from the average age calculations (Table S1, supplementary data). We considered samples with standard error of the mean (SEM) values >8 to show dispersion and removed outliers to reduce the SEM to <8, except for small (<3) n datasets. This is an arbitrary cut-off, but we emphasize that this approach is used here for data display and summary only, not for modeling or interpreting thermal histories. In two cases (ZE09006, ZE09079) where 7 of 8 grain ages were older than the MDA of the sample, the anomalously young grain (significantly younger than the MDA) was also considered an outlier and removed, although the SEM for those samples was <8.

171 single crystal ZHe ages range from ca. 222 to 42 Ma. Three Laberge Group samples (15KNA010, ZE09006, and ZE09016), one Lewes River Group sample (ZE09079), and two Tantalus Formation samples (15KNA031 and 15KNA034) yield ZHe single crystal ages equivalent to or older than their MDA ranges. These data (colored black in Fig. 4A) are from rock units located in the northernmost part of the basin and in the low thermal maturity region of the Atlin Lake area (Fig. 3C) and form a Triassic–Early Jurassic age peak, here designated as Group 1 (Fig. 4B). The remaining 29 ZHe samples yield single crystal ages that are generally younger than their MDA ages or age ranges, and form two age peaks (Fig. 4B). The older of these is Late Jurassic and includes Laberge Group samples 15KNA032, 15KNA035C, 15KNA027, 16EB384, 16RGI45307, 16BvS-256, Cache Creek terrane sample 18KNA067, and Bowser Lake Group sample 17BvS25-242, designated as Group 2 (Fig. 4B). The younger age peak is Cretaceous–Paleogene, ca. 120–50 Ma, and includes Laberge Group samples 18KNA016, 18KNA021, 18KNA022, 18KNA027, 18KNA030, 15KNA030, 15KNA029, 15KNA028, 15KNA005, 15KNA012, 15KNA014, S3, 16EB257, 18KNA024, 18KNA036, 18KNA046, and ZE111001, Lewes River Group samples 17EB225 and 16EB421, Cache Creek terrane sample 18KNA065, and Bowser Lake Group sample 15KNA023, designated as Group 3 (Fig. 4B). Both the Late Jurassic and Cretaceous–Paleogene ZHe cooling age groups are distributed along the length of the basin, with the youngest (Paleogene) ages of the latter group apparently confined to a central portion of the basin in southern Yukon, as well as sample ZE111001 from the footwall of the King Salmon thrust in the Atlin Lake area (Fig. 5). There is a gap in sampling between southern Yukon and the Atlin Lake area, such that the full extent of the Paleogene ZHe age pattern is undefined.

Twelve AFT central ages range from ca. 95–30 Ma, with all but the oldest (Cache Creek sample 18KNA067) and youngest (Laberge Group sample 18KNA046) ages forming a narrow population with dominant age peak at ca. 50 Ma (10 ages ranging ca. 62–45 Ma; Fig. 4). Radial plots are provided in the supplementary data.

157 AHe ages range from 228 to 13 Ma (one additional grain yielded 510 Ma), forming a broad age peak at ca. 50–40 Ma (Fig. 4). In general, age data from all samples are significantly younger than MDA ranges for those units: only four single grains from four different samples were equivalent or older.

The along-strike distribution of thermochronological data for the Whitehorse Trough (Fig. 5) indicates that all the identified ZHe age groups are represented across the basin, suggesting that to a first order, the rocks may have experienced a common thermal evolution. In detail, there are likely differences for different portions of the basin, as suggested in preliminary review of a subset of this data (see Kellett and Zagorevski 2021). Gaps in representation (e.g., no Cretaceous ZHe ages in the southern portion of the basin, no Late Jurassic ZHe in the central, i.e., Atlin Lake, portion of the basin) may be real or a result of limited sampling.

Elevation swaths along the length (Fig. 5A) and width (Fig. 5B) of the basin show that elevation ranges between sea level and 2500 m, with our samples from ∼560 to 1900 m (Table 1, Fig. 6C). Within our range of sampling elevations, ZHe age Groups 1–3 are represented at both low and high elevation, except for the youngest (<50 Ma) data, which are only present in low elevation samples (Fig. 5C). AFT and AHe data are also present among both high and low elevation samples, with no apparent correlation between age and elevation at this scale and sample density.

The goal of this study is to constrain, in time and magnitude, the two structurally controlled burial events that contributed to the tectonic and thermal evolution of the Whitehorse Trough. To do this, we apply forward and inverse modeling approaches to better understand the distribution of the low-temperature thermochronology data presented above, especially the broad but punctuated age record preserved in the ZHe system. Our working hypothesis is that the different age groups observed in the ZHe data are the result of different magnitudes of heating among the samples during these two burial events. We applied the complementary modeling tools, HeFTy (v2.1.4, © 2022, Richard Ketcham) and TcPlotter (Whipp et al. 2022) for forward modeling.

We constructed three bounding, representative T-t paths in HeFTy that capture the known constraints of the geologically defined history of the Laberge Group (Paths A–C, Fig. 6A). Geological constraints are: c1-Source rock cooling/exhumation; c2-Early Jurassic deposition into the basin; c3-Middle Jurassic heating during sedimentary and structural burial; c4-Middle–Late Jurassic exhumation that produced the top Laberge angular unconformity, followed by c5-deposition of Bowser Lake Group and Tantalus Formation strata on top of the Laberge Group, and heating of unknown degree/duration. Constraints c1, c2, c4, and the model end point are fixed, while c3 and c5 are varied between end member values. Model constraint c1, “initial” igneous cooling conditions, represents rapid cooling of Triassic to Early Jurassic igneous detritus that are the dominant sources of Laberge Group strata. Evidence for rapid cooling of detrital sources includes short lag times between dominant detrital zircon age peaks and MDAs, short lag times between metamorphic cooling ages and MDAs, and abundant volcanic, hypabyssal, and plutonic clasts indicating collapse of an active volcanic arc adjacent to the basin (Kellett and Zagorevski 2022). Model constraint c2, deposition in the basin, is informed by the range of MDAs reported for Laberge Group (202–178 Ma; Kellett and Zagorevski 2022). Model constraint c3, the first burial event (B1 in Fig. 6), allows burial-related heating during deposition of Laberge Group and shortening into a fold and thrust belt. We varied this constraint between minimum and maximum heating conditions of 80 °C and 200 °C, corresponding to the range of peak thermal maturity conditions in the basin (Fig. 3). Model constraint c4 was informed broadly by the top-Laberge unconformity and depositional age range of Bowser Lake Group and Tantalus Formation (Fig. 1). Model constraint c5 is not well defined in time or temperature by geological constraints. The timing of peak temperature has been arbitrarily fixed at 110 Ma, and as for c3, we varied this constraint between minimum and maximum heating conditions of 80 °C and 200 °C, corresponding to the range of peak thermal maturity conditions in the basin (Fig. 3).

Forward modeled end-member Paths A–C (Fig. 6A) are identical except for peak T conditions reached during c3 and c5. They include the following nodes: 400 °C at 200 Ma (constraint c1), 0 °C at 180 Ma (c2), peak burial B1 T at 170 Ma (c3), 0 °C at 150 Ma (c4), peak burial B2 T at 110 Ma (c5), and 0 °C at 0 Ma (Fig. 6B). These paths are intended to capture the overall shape of the Laberge Group T-t history based on existing geological information, to explore whether such a set of T-t paths could produce the pattern and range of ZHe ages measured in this study, by varying the magnitude of heating during B1 and B2. We applied a zircon equivalent spherical radius (ESR) of 50 µm (range in our data are 31–89 µm), and effective U (eU) concentration of 100 ppm (range in our data are 0.04–3537 ppm), and the zircon radiation and annealing model of Guenthner et al. (2013). We ran 27 models in which first c3 peak T was increased in increments of 10 °C from 80 to 200 °C, and then held at 200 °C, while c5 peak T was increased in increments of 10 °C from 80 to 200 °C.

ZHe ages calculated from the forward models are shown in Figs. 6B and 6 reported in the supplementary data. Path A, in which peak temperatures during B1 and B2 reach only 80 °C, yielded a ZHe age of 187 Ma. There is an incremental decrease in ZHe age for increasingly higher peak temperature during B1, especially between 120 and 150 °C (Fig. 6A), eventually producing a ZHe age of 159 Ma in Path B, for which the peak temperature during B1 is 200 °C, and that for B2 remains fixed at 80 °C. The ZHe age incrementally decreases from 159 Ma to 83 Ma as the peak temperature for B2 is raised progressively from 80 to 200 °C (and B1 peak T remains fixed at 200 °C; Fig. 6A): the largest shift in ZHe age occurs across the B2 peak T range of 120–150 °C. ZHe age ranges from 187 to 83 Ma across the 27 models.

The ZHe thermochronometer can be highly sensitive to grain size (Reiners and Farley 2001) and eU (Flowers et al. 2009; Guenthner etal. 2013; Guenthner 2021), or it can be relatively insensitive to these parameters, depending on the particular T-t history (Whipp et al. 2022). We fixed ESR and eU at representative values in the above forward models. With TcPlotter, we explored the sensitivity of the end-member T-t Paths A, B, and C under the full range of ESR and eU observed in our data. We used TcPlotter version 0.3.3 (doi:10.5281/zenodo.7684503) to overlay our single crystal ZHe data on contoured ZHe age plots in which ESR and eU form the plot axes, scaled to the observed ranges of our dataset, for each path (Fig. 7). We plotted data from Group 1 against Path A, Group 2 data against Path B, and Group 3 data against Path C. Tantalus Formation and Bowser Lake Group data were omitted from the plots because they experienced a single heating/burial event and are therefore not represented by any of the three forward models.

The combination of all three thermochronometers applied to an individual sample provides the most detailed information regarding rock thermal histories, so the 12 samples for which ZHe, AFT, and AHe ages are available were prioritized for inverse thermal history modeling using HeFTy, as guides to interpret the broader dataset. This subset of Laberge Group and one Cache Creek terrane sample spans a large extent of the Whitehorse Trough, from near Carmacks (Yukon) to Atlin (BC) (Table 1, Figs. 1 and 4). We lack full multi-thermochronometry data from the Lisadele Lake and Dease Lake areas in the southernmost regions of the basin (Fig. 4). To explore thermal histories for these regions, we ran inverse models based on ZHe and AHe data only. From the Dease Lake area we modeled 16RGI45307, a Laberge Group sample in the footwall of the King Salmon thrust/hanging wall of the Kehlechoa fault. From the Lisadele Lake area, we modeled ZHe and AHe data from nearby Laberge Group samples 15-KNA-023 and 15-KNA-016, respectively, which exhibit stratigraphic continuity (footwall of the King Salmon thrust). We also selected Group 1 sample 15KNA010 for inverse modeling, to explore Path A-type thermal histories, although we lack AFT data for this sample.

HeFTy inverse controlled random search models were run using the zircon and apatite radiation damage and annealing models of Guenthner et al. (2013) and Flowers et al. (2009), respectively. He stopping distances reported in Ketcham et al. (2011), and the apatite fission track annealing model of Ketcham et al. (2007). HeFTy allows up to seven individual sample data. The Jurassic–Cretaceous thermal history under investigation is most comprehensively represented by the ZHe data, with AFT and AHe ages predominantly defining latest cooling in the early Paleogene. Consequently, our modeling approach prioritizes the single crystal ZHe data, with the AFT age, and one AHe single grain age included to constrain the final cooling stage. Each model run in HeFTy includes 2–5 single crystal ZHe ages, AFT age, and one single crystal AHe age, with the selected AHe grain most closely representing the mean AHe age for the sample. Where not all ZHe data could be included (i.e., samples with > 5 single crystal ages such as 15KNA005), grains were selected to span the range in eU and grain size. Model inputs and constraints summarized here are detailed in the supplementary data on a sample-by-sample basis, including any available depositional age constraints from published detrital zircon U–Pb and biostratigraphy. He diffusion in zircon may be significantly affected by accumulated alpha decay damage, which is a function of the crystal eU and the time the crystal has spent under conditions too cool to anneal the damage (Guenthner et al. 2013; Guenthner 2021; Whipp et al. 2022). Consequently, and to compare with the forward models, each sedimentary sample is modeled from an initial, “fully annealed” high-temperature zircon condition, which is intended to represent the time at which the detrital zircon was exhumed from its source region (i.e., constraint c1 from the forward models).

Inverse thermal modeling in HeFTy can be used to assess whether hypothesized geological and thermal histories are compatible with thermochronological datasets. In particular, hypotheses including heating require user-defined constraints. Based on the insights gathered from forward modeled thermal histories that explored different peak heating conditions during two simulated stages of heating/burial, with intervening exhumation, we applied as broad constraints for c1-5 as possible, allowing us to quantify maximum allowable heating/burial conditions during B1 and B2 on a sample-by-sample basis, as well as cooling/exhumation rates and timing during the two post-depositional cooling episodes. We developed inverse thermal models that guide T-t paths for all modeled samples through the same constraints as the forward models: c1-Source rock cooling/exhumation, c2-Early Jurassic deposition into the basin, c3-Middle Jurassic heating (simulating structural and sedimentary burial, analogous to B1 of the forward models), c4-Middle–Late Jurassic cooling (simulating exhumation producing the top Laberge angular unconformity and allowing deposition of Bowser Lake Group and Tantalus Formation strata), and c5-Early–Late Cretaceous heating events of varying degree and duration (analogous to B2 of the forward models). Thus, each inverse model has two heating stages and three cooling stages. The constraints used to guide heating achieved in each burial phase are intentionally broad (80–200 °C) to allow for different sample-specific peak temperature conditions during the two heating events, where necessary. Two models (Lisadele Lake and 18KNA046) allow for even higher temperatures during B1 and B2, to 250 °C, based on the observation that rocks in the footwall of the King Salmon thrust record the highest thermal maturity values (Fig. 3). Also, the c4 temperature window is broader (0–40 °C) than allowed during c2, initial Laberge Group deposition, since not all samples necessarily reached the surface during exhumation and development of the Laberge Group unconformity. For two samples (18KNA036 and 18KNA046), it was necessary to guide final cooling to ranges that include sample AFT and AHe ages within a broad temperature window (140–40 °C) to avoid restricting possible cooling paths. The resulting “acceptable” and “good” thermal histories are shown in Fig. 8, and the mean paths from all models are coloured by Group (1–3) and overlaid for comparison in Fig. 9. “Acceptable" fits correspond to a goodness of fit value of 0.05–0.49, while “good” fits correspond to a value of 0.50 or higher.

Inverse models are presented and discussed in the context of ZHe Groups 1–3. Thermal history modeling was applied to one sample from Group 1: 15KNA010. Thermal histories for five Group 2 samples were modeled: Laberge Group samples 15KNA032, 15KNA035C, 15KNA027, 16RGI45307, Cache Creek terrane sample 18KNA067, and the combined Lisadele Lake samples (Fig. 8A). These models each show mean thermal history paths consistent with reaching peak thermal conditions during B1. Thermal histories for eight samples from Group 3 were modeled. Group 3 models with Paleogene ZHe ages are presented together because of the common Paleogene rapid cooling event they define (Fig. 8C). Despite this difference, all the Group 3 models show much greater B2 temperatures compared to the Group 2 samples, up to and exceeding 200 °C. Whereas the magnitude of peak temperature is relatively well-defined, the timing of peak temperature is less well defined, with both Early (e.g., 15KNA005) and Late (e.g.,18KNA024) Cretaceous times favored by our models. 18KNA036 and 18KNA046 models show long burial residence times during the Cretaceous, between ca. 100 and 50 Ma. Some ZHe individual crystal data prevented any thermal history solutions within the parameters we allowed. These are indicated in the model log (supplementary data) and highlight an additional source of uncertainty in these data. We suggest the most likely explanation is unaccounted for zoning in U and Th, resulting in inaccuracy in α-ejection corrections. We assumed uniform U and Th distribution in the dated zircon, and in calculating the α-ejection corrections. Applying this approach on grains with high U/Th rims would result in an underestimation of the single grain U–Th/He age. In contrast, applying this approach on grains with high U/Th cores would result in an overestimation of U–Th/He age.

As first B1 and then B2 peak burial temperature was varied in the HeFTy forward models, the ZHe age incrementally decreased from 187 Ma (older than the c1 deposition constraint of 180 Ma) down to 83 Ma (younger than the last burial/heating stage). These results demonstrate that ZHe ages spanning 100 myr, a comparable range to our observed ZHe dataset, can be produced from a single T-t history representative of Laberge Group only by modifying the magnitude of heating during the two burial/heating stages. The sensitivity of ZHe age during 120–150° C heating indicates that for these conditions, the ZHe partial retention window is approximately within this temperature range (Fig. 6B). Further, the distribution of ages from this set of forward models shows that as peak temperatures are varied within this T-t geometry, Early Jurassic, Late Jurassic, and mid-Cretaceous ages (the three cooling stages) are overall favored (Fig. 6A).

The Path A plot developed with TcPlotter shows that ZHe age under such a thermal history should be relatively insensitive to ESR, and differences in eU within the natural ranges captured by our sample set could produce up to 20 myr difference in ZHe age (Fig. 7A). That is, under a Path A-type thermal history (rapid exhumation of Late Triassic/Early Jurassic source material, low T heating during B1 and B2) variations in single crystal eU may result in up to 20 myr of within-sample age differences. Most of Group 1 ZHe single crystal ages fall within this Path A age range (Fig. 7A), though there is dispersion towards younger ages (i.e., the data points are more commonly darker/younger than the colour surrounding them in the plot, which represents the model age for those ESR and eU values). The Path B plot, in which the ZHe system opens during B1, similarly shows insensitivity of ZHe age to our natural range of ESR and ∼16 myr of age variation due to eU differences (Fig. 7B). We compared Group 2 ZHe single crystal ages with this Path B, and while many overlap, they also skew to younger than model ages. These younger data could reflect intermediate samples between Path B and Path C (i.e., heating during B2 was > 80 °C but < 200 °C), among many other possible scenarios. The Path C scenario (Fig. 7C) shows that the spread in ESR and eU of our Group 3 ZHe single crystal data fall mostly within the 80 Ma to 96 Ma age contours, a ∼16 myr difference in ZHe age. However, it is apparent from the overlay of Group 3 ZHe single crystal ages that they show much greater variability than the range predicted by varying ESR and eU, and they do not correspond well with the predicted ages for Path C, with most ages being older or younger than modeled. This suggests that B2 timing and duration in our forward models are not a particularly good fit for our data. Accordingly, there may be considerable variation in B2 timing, duration, and/or peak temperature across the Whitehorse Trough such that one path alone is not representative.

The inverse model for Group 1 sample 15KNA010 shows the lowest maximum temperatures of any of the models, for both burial/heating stages (Figs. 8A and 9). We infer that Group 1 samples, while buried along with the rest of the sedimentary rocks, were only heated to temperatures intermediate between the ZHe retention zone and AFT annealing zone (Fig. 8A) during B1 and B2, allowing them to preserve ZHe dates that exceed their depositional ages. Path A thus is a good analogue for Group 1 samples. Group 2 inverse models (Figs. 8B and 9) show maximum temperatures during B1 heating of Laberge Group and Jurassic Cache Creek terrane strata, with lower temperatures reached during B2 heating, similar to forward model Path B. These models best define the magnitude and timing of the B1 burial/heating event. All six models show an envelope of paths compatible with 20–30 myr of post-depositional progressive heating, reaching maximum temperatures of 150–200 °C in the Middle to Late Jurassic, and cooling over ∼15 myr during the Late Jurassic. Maximum heating and cooling rates for B1 are thus approximately 10 °C/myr and −13 °C/myr, respectively. Group 3 models (Figs. 8C, 8D, and 9) produced T-t envelopes in which peak temperatures during B2 are as high as or higher than during B1 (again, ∼150–200 °C), allowing for partial to full resetting of the ZHe ages during this second heating event. The duration of heating to reach B2 is between ∼30 and 50 myr, such that the maximum heating rate for B2 is ∼7–4 °C/myr. The Group 3 samples with Paleogene ZHe ages are similar except they show longer duration at high temperature during B2, and all three show rapid and late (Paleogene) final cooling, resulting in near coincidence of all three thermochronometers.

We find that these inverse thermal histories satisfy both the available geological constraints and our low temperature thermochronometer dataset. Though all samples contribute to our understanding of the Whitehorse Trough thermal history, Group 1 samples most clearly define pre-depositional cooling of Stikinia and Yukon–Tanana terrane source rocks, Group 2 samples best define B1 heating rates and magnitudes, and Group 3 samples best define B2 heating rates and magnitudes. Consequently, we can explore the spatial distribution of established high and low temperature conditions during B1 and B2 (Fig. 10). The B1 pattern suggests that central parts of the basin experienced higher temperatures than the northern edge and parts of the western flank. This may be because these regions correspond to the deepest regions of the trough (van Drecht et al. 2022), or they occupied a deep (e.g., footwall) structural position in the Late Jurassic fold and thrust belt (Fig. 5B), or both. The B2 pattern is complex and shows a domanial pattern of high versus low heating that could be explored further with detailed sample transects across regions such as the Atlin Lake area, or the northern Yukon section of the basin, with accompanying vertical sample transects within those domains.

The heating of sediments and sedimentary rocks occurs via burial (depositional, structural, or both) or proximity to magmatic heat sources and/or hydrothermal circulation cells. In framing our forward and inverse models, we have ascribed heating and cooling to burial and exhumation, respectively. However, the Intermontane terranes host extensive plutonic suites (Fig. 1), including syn-collisional, Late Triassic to Middle Jurassic plutons that overlapped with the timing of Whitehorse Trough deposition (e.g., Sack et al. 2020; Colpron et al. 2022), and it is therefore necessary to ascertain to what degree they may have influenced the thermal history of proximal samples. ZHe data are plotted against the shortest mapped distance to Jurassic, Cretaceous, and Paleogene plutons to accomplish this examination (Fig. 11).

The closest samples to Jurassic intrusions are 16BvS35256 (ZHe_mean = 144 Ma) and 16RGI45307 (mean ZHe = 158 Ma), which are 1.5 km and 2.4 km away (Fig. 11A), respectively, from the Snowdrift Creek pluton, which has a zircon U–Pb crystallization age of 160.43 ± 0.16 Ma, and K–Ar cooling ages between ca. 161 and 147 Ma (van Straaten et al. 2017). The proximity and slightly younger ZHe ages could reflect heating in a contact aureole around the pluton rather than sedimentary or structural burial. We note that the sampled outcrops are also in the footwall of the King Salmon thrust, which was active during Late Jurassic (Mihalynuk et al. 1995). The next closest sample with a Jurassic ZHe age is >8.5 km away from its nearest pluton neighbor. The modeled Early Jurassic T-t history for 16RGI45307 is similar to the other T-t histories for Group 2 samples. We conclude that it is challenging in this case to resolve the competing contributions of depositional and structural burial versus pluton heating from thermochronological data or T-t histories alone. In the case of the samples with Jurassic ZHe ages, it seems most plausible that the modeled Late Jurassic rapid heating and cooling event was produced via burial by overlying sedimentation, structural burial during development of the W-verging fold belt, and perhaps locally, Jurassic magmatism, followed by exhumation and local pluton cooling.

The closest sample to Cretaceous intrusions is 18KNA065 (ZHe_mean = 83 Ma), which is 1.1 km from a Teslin suite (123–115 Ma) pluton (Colpron 2011; Fig. 11B). This sample was not modeled, but the ZHe and AHe ages are similar to data from Group 3 modeled samples such as 18KNA022, which are situated several kilometers away from Cretaceous plutons. The next closest sample with Cretaceous or older ZHe age, at 2.6 km distance from a Cretaceous pluton, is ZE09079, which yielded Jurassic ZHe ages (Fig. 11B). We conclude that heating and cooling during the Late Cretaceous must largely be the product of burial (depositional and possibly structural) followed by exhumation, with magmatic heating playing a minor, local role.

The nearest sample to a Cenozoic intrusion is 15KNA027 (ZHe_mean = 149 Ma), which is 2.3 km from a Nisling Range suite pluton with a zircon U–Pb crystallization age of 53.6 Ma (Hart 1997). We conclude that the ZHe system for this sample was not sensitive to magmatic heating. However, the AFT and AHe ages for this sample are 54 and 42 Ma, respectively, and may have been reset during pluton emplacement. Thus, we conclude that the ZHe Cenozoic ages in our sample set are not affected by Cenozoic magmatic heating, although some AHe and AFT data may be affected. As with the older magmatic suites, despite the proximity of some samples to magmatic bodies with slightly older crystallization ages, there is no systematic correlation of ZHe thermochronometric age versus proximity to the plutonic suites (Fig. 11). Thus, although we cannot fully rule out a contribution of magmatic heating to the resulting age pattern, it is likely to have played a minor role in the thermal history of the Whitehorse Trough.

While we can compare our sample locations to mapped plutonic suites, we cannot rule out the possibility of unexposed plutons that may have contributed to Jurassic to Cenozoic thermal history of the Whitehorse Trough. Similarly, we cannot rule out the possibility that the various Mesozoic plutonic and associated volcanic suites generated large hydrothermal circulation cells that may have extended tens of kilometers beyond the igneous bodies themselves and may have accessed structural pathways. For example, it has been proposed that southern Yukon, including the Whitehorse Trough region, experienced an extensive hydrothermal event during ca. 70 Ma emplacement of the Carmacks Group volcanic succession, thought to be a product of mantle upwelling (Wynne et al. 1998). Our data do not resolve a specific heating event ca. 70 Ma, but depending on the duration of inferred mantle upwelling, the “final cooling” event marked by similar Cenozoic AFT and AHe ages across our entire study area could reflect the terminal stage or waning of that geodynamic/thermal event. More recently, it has been proposed that collision between the Stikinia arc and Yukon–Tanana terrane was followed by lithospheric delamination, asthenospheric upwelling, and emplacement of the 174–168 Ma Bryde and 4th of July instrusive suites in Yukon and British Columbia (Colpron et al. 2022). These plutons locally intrude the Whitehorse Trough, Stikinia, and Cache Creek terrane and cross-cut the SW-verging thrust belt that imbricates these tectonic elements (e.g., Bickerton et al. 2020; Colpron et al. 2022). The juxtaposition of hot asthenosphere against the base of the crust after delamination may result in a short-lived heating and thermal re-equilibration of the crust (Göğüş and Ueda 2018). Isostatic uplift due to thermal buoyancy could have triggered extensional deformation and increased surface erosion. These competing processes may have contributed to the relatively rapid heating and cooling cycle recorded in B1, in addition to depositional and structural burial followed by basin inversion and shortening.

The long-term tectonic evolution of the Whitehorse Trough is reinterpreted below based on our new and published thermochronological constraints from detrital minerals and thermal history models in the context of the Intermontane terranes (Fig. 12).

Ca. 210–200 Ma: the latest Triassic was marked by carbonate shelf deposition along the eastern flank of the Lewes River arc as subduction beneath Stikinia waned. The Yukon–Tanana terrane, Stikinia/Quesnellia, Atlin terrane (including the Kutcho arc), and Cache Creek terrane were amalgamated (Nelson et al. 2022; Colpron et al. 2022; see Atlin and Cache Creek terrane as redefined in Zagorevski et al. 2021) before or during this time interval. Latest Triassic to Early Jurassic subduction to the south formed the Hazelton arc on Stikinia in northern BC (modern coordinates, not shown in Fig. 12), while in Yukon, Stikinia terrane collided with Yukon–Tanana terrane, possibly “end-on,” i.e., at a high angle of collision as proposed by George et al. (2021), which resulted in Barrovian and high-pressure metamorphism (Currie and Parrish 1993; Dusel-Bacon et al. 2002; Canil et al. 2006; Berman et al. 2007; Joyce et al. 2015; Kellett et al. 2018; Gaidies et al. 2021; Soucy La Roche et al. 2022), emplacement of the peraluminous Minto suite plutonic rocks at mid-crustal depths (Sack et al. 2020; Colpron et al. 2022), and sinistral strike-slip opening of the Whitehorse Trough (van Drecht et al. 2022).

Ca. 200–175 Ma: the Whitehorse Trough opened and broadened through trans-tension along basin-bounding faults such as the proto-Big Salmon–Teslin and Llewellyn faults (van Drecht et al. 2022; Fig. 12B). Early Jurassic arc collision and crustal thickening were followed by the exhumation of Stikinia and Yukon–Tanana terrane basement rocks (e.g., Willow Lake fault; Knight et al. 2013), and syn-collisional plutons. Laberge Group strata were deposited into the Whitehorse Trough, fed by detritus from the collision zone to the north. High pressure metamorphic rocks were rapidly exhumed along the Yukon Tanana/Stikinia suture and medium pressure metamorphic rocks of Yukon–Tanana terrane were also exhumed (Berman et al. 2007; Staples et al. 2014; Kellett et al. 2018). Exhumation (Colpron et al. 2015, 2022) and erosion of this collision zone supplied limited Yukon–Tanana terrane metamorphic detritus to the Whitehorse Trough (Canil et al. 2006; Kellett et al. 2018).

Ca. 175–160 Ma: During the Middle–Late Jurassic, with amalgamation of Yukon–Tanana/Stikinia/Quesnellia/Cache Creek/Atlin terranes completed, the Whitehorse Trough underwent inversion via shortening, producing a W-vergent thin-skinned fold and thrust belt, involving the initiation of the southwest-vergent King Salmon thrust fault, and the juxtaposition of Atlin and Cache Creek terranes over Stikinia along the Nahlin and other southwest vergent thrust faults (Figs. 1C, 9C, and 12C). Early Jurassic sinistral strike-slip faults in Yukon, such as the Braeburn, Goddard, and Laurier faults (van Drecht et al. 2022), may have accommodated some Late Jurassic shortening (e.g., White et al. 2012). The vergence of the fold-and-thrust belt support westward advance of North America during the early break-up of supercontinent Pangea (Monger and Gibson 2019). Laberge Group strata were buried to depths up to 6–10 km, heating at 10 °C/myr during regional shortening, and ZHe ages were reset in the most deeply buried strata (this study). Mid-crustal rocks were metamorphosed (Staples et al. 2016), consistent with crustal thickening in the hinterland of the growing Cordilleran orogen. Laberge Group strata formed part of an orogenic topography after shortening and were likely subjected to orographic erosion via surface erosive processes (e.g., Willett 1999), with cooling of deeply buried Laberge Group rocks occurring up to −13 °C/myr. The top Laberge Group unconformity probably developed because of this regional exhumation. This could be when the northern Intermontane terranes developed into a synformal, “flower” type structure (e.g., Calvert et al. 2017). Proposed lithospheric delamination during this period (Colpron et al. 2022) might have contributed to the short cycle of heating and cooling, by: pulling down crust in the “drip” phase (contributing basin accommodation space and increasing burial); heating the crust during the delamination phase (accelerating heating of the buried sediments); and isostatic rebound of the crust (increasing elevation of sediments that could result in more aggressive erosion of the fold and thrust belt).

Ca. 160–140 Ma: Following Late Jurassic exhumation of the Laberge Group and Cache Creek Jurassic strata, the Whitehorse Trough region underwent renewed subsidence, perhaps related to a backarc setting as east-dipping subduction initiated beneath the northern Insular terranes (Colpron et al. 2015; Beranek et al. 2017; Fig. 12C), or to a period of transtension within the Intermontane terranes, as the North American craton shifted from a westward to north-northwestward trajectory (Monger and Gibson 2019). Subsidence allowed for local Late Jurassic to Early Cretaceous deposition of Bowser Lake Group strata in northern BC and Tantalus Formation strata in Yukon, reburying and reheating the exhumed Laberge Group rocks. Elevated topography created by shortening of Cache Creek terrane strata during the Late Jurassic basin inversion became subject to erosion, delivering chert into these successor basins from the east.

Ca. 140–120 Ma: Laberge Group and overlying Tantalus Formation and Bowser Lake Group strata underwent burial during the Early Cretaceous. The variable reheating of Laberge Group rocks occurred at a rate of up to 4–7 °C/myr, locally reaching temperatures of up to 200 °C and burial depths of ∼6-10 km, subject to apparent structural controls. For example, greater heating of footwall strata compared to hanging wall strata across the King Salmon fault, including resetting of detrital ZHe in Bowser Lake Group strata in the Dease Lake and Lisadele Lake regions, appears to require reactivation of Jurassic structures in the vicinity of the King Salmon fault (Figs. 4, 6B, and 11E). Early Cretaceous metamorphism of mid-crustal rocks to the north of the Whitehorse Trough (Staples et al. 2014, 2016) further indicates an episode of crustal thickening and shortening in the Intermontane terranes. The shift in trajectory of the North American plate from NW or N during 160–140 Ma to SW during 140–120 Ma could have favored this shortening and thickening phase (Fig. 12; Monger and Gibson 2019).

Ca. 120–70 Ma—Whitehorse Trough strata resided in the upper crust, with local exhumation likely controlled by dextral and transpressional structures such as the Braeburn and Teslin faults in central Yukon (e.g., White et al. 2012). There is evidence for lithosphere-scale extension within the northern Canadian Cordillera during the mid-Cretaceous, including the development of metamorphic core complexes like the Australia Mountain domain in central Yukon (Pavlis 1989; Pavlis et al. 1993; Staples et al. 2016). During this period, the Whitehorse Trough reached its present-day structural geometry, bisected by Jurassic thrusts and Cretaceous strike-slip faults (Fig. 11F). Paleomagnetic evidence for widespread hydrothermal circulation in south and central Yukon (Wynne et al. 1998), latest Cretaceous emplacement of the Carmacks volcanic suite (Grond et al. 1984), and evidence for post 120 Ma low pressure–high temperature metamorphism in Yukon–Tanana terrane west of Whitehorse Trough in the Atlin, BC area (Soucy La Roche et al. 2022) suggest a heightened geothermal gradient during this period.

Ca. 65–50 Ma: Paleocene to Eocene rapid cooling and exhumation is the final tectono-thermal event recognized by Whitehorse Trough strata. AFT and AHe ages across the >600 km strike-length of the basin show little variation, with peaks at ca. 50 Ma and 50–40 Ma, respectively, being closely spaced in time. The Whitehorse Trough lies between the Tintina and Denali faults, which are post-Cretaceous dextral strike-slip faults with >400 km displacement each (e.g.,Gabrielse et al. 2006). The Tintina fault extends to the Moho (Cook et al. 2004). It is possible that slip along those structures allowed for block uplift and rapid exhumation of Whitehorse Trough and surrounding rocks. However, several competing factors during Paleocene–Eocene, including changes in plate boundaries and trajectories (Vaes et al. 2019), and thermal relaxation after a magmatic flare-up within the Coast Mountains batholith (Cecil et al. 2021) could have contributed to this cooling pattern as well.

A new low-temperature thermochronology dataset is presented for the Whitehorse Trough of the northern Canadian Cordillera. ZHe ages from detrital zircon in 35 samples fall into three main groups: Triassic–Early Jurassic (ca. 200–180 Ma), Late Jurassic (ca. 160–150 Ma), and Cretaceous–Paleogene (ca. 120–60 Ma), while AFT ages are closely distributed around 50 Ma and AHe ages around 50–40 Ma, both from detrital apatite. The distribution of samples representing these different age groups and known geological constraints on the evolution of the Whitehorse Trough suggest that the broader basin experienced a thermal history involving two phases of burial and heating, with structural and stratigraphic positions controlling the degree of reheating. Temperature conditions during both burial stages locally exceeded 150 °C. We propose that both events are attributable to sedimentary and structural burial during accretion of the Intermontane terranes. Regional Paleogene AFT and AHe ages, and local Paleogene ZHe ages warrant future investigation.

We gratefully acknowledge lab support from Roman Kislitsyn, field support from Catherine Mottram and Mathilde Banjan, access to samples and discussions with Maurice Colpron and Bram van Straaten, project support from Steve Irwin, and Discovery and Capital helicopter support for safe and efficient field sampling. This work benefitted from critical feedback by anonymous reviewers. This is Natural Resources Canada contribution no. 20220220.

All data presented in this study are available in the supplementary data.

Conceptualization: DAK, IC, AZ

Data curation: DAK, IC, KD

Formal analysis: DAK, IC, DG, KD

Funding acquisition: DAK

Investigation: DAK, AZ

Methodology: IC, DG

Project administration: DAK

Resources: DAK, AZ

Software: DAK

Validation: DAK, KD

Visualization: DAK, IC

Writing – original draft: DAK

Writing – review & editing: DAK, IC, AZ, DG, KD, LB

This work was supported by the Geological Survey of Canada's Geomapping for Energy and Minerals II and GeoNorth programs.

Supplementary data are available with the article at https://doi.org/10.1139/cjes-2023-0082.