150 Myr of Episodic Metamorphism Recorded in the Yukon-Tanana Terrane, Northern Canadian Cordillera: Evidence from Monazite and Xenotime Petrochronology

The in situ analysis of isotopic and trace element contents in petrochronometers combined with detailed characterization of their textural relationships with metamorphic index minerals is a powerful tool to link time constraints with pressure (P) and temperature (T) conditions [1, 2]. This analytical approach is particularly useful in polymetamorphosed domains where petrographic relationships can be ambiguous or difficult to interpret. Aluminous metamorphic rocks are excellent targets for these studies because they commonly contain petrochronometers such as monazite or xenotime, and a wide range of metamorphic index minerals that can provide first-order constraints on P–T conditions (e.g., [3]). Al₂SiO₅ polymorphs (andalusite, kyanite, and sillimanite) are key minerals to reconstruct the P–T evolution of metamorphic domains around the world (e.g., [4–6]) because the P–T conditions of the transitions between these polymorphs are independent of the geochemical composition of the equilibration volume and therefore less affected by geochemical variations due to open system behavior.

The Yukon-Tanana terrane of the northern Canadian Cordillera is characterized by a polymetamorphic evolution that results from its complex tectonic interactions with the Laurentian margin and other adjacent microcontinents, arc-back arc systems, and oceanic domains. The basement of this terrane locally contains abundant aluminous metamorphic rocks that provide an excellent opportunity to reconstruct its polyphase tectonometamorphic evolution. In this contribution, we investigated a locality of the Yukon-Tanana terrane near Atlin Lake in northwestern British Columbia, Canada (Figure 1) that experienced several episodes of metamorphism recorded by all three Al₂SiO₅ polymorphs [7, 8]. There, two contrasting Al₂SiO₅ crystallization sequences were proposed based on petrographic observations (And→Ky→Sil, [7]; Ky→And→Sil, [8]; mineral abbreviations after Whitney and Evans [9]), and the lack of robust geochronological constraints linked to P–T conditions during each metamorphic phase hindered their interpretation in a tectonic context. We used detailed in situ monazite and xenotime petrochronology to constrain crystallization of Al₂SiO₅ polymorphs and other metamorphic index minerals and to assess the qualitative P–T–t evolution of the region. Our monazite and xenotime petrochronological investigation indicates that the P–T–t path successively reached the stability fields of kyanite, sillimanite and andalusite during three distinct metamorphic episodes over 150 Myr from the Permian to the Early Cretaceous. These data enable evaluation of local sedimentation models and regional tectonic models involving the interaction of the Yukon-Tanana terrane with adjacent terranes. Moreover, the current study demonstrates the usefulness of combining detailed textural analysis of petrochronometers in metamorphic rocks with their in situ isotopic and trace element analysis as a tool to understand the tectonic evolution of the lithosphere.

The Canadian Cordillera comprises a collage of Paleozoic to Cenozoic terranes that experienced a complex magmatic, tectonic, and metamorphic history prior to, during, and following their accretion to Laurentia [10–12] [13, 14]. The Yukon-Tanana terrane is the largest peri-Laurentian terrane (i.e., originated along the western margin of Laurentia) exposed in eastern Alaska, central Yukon, and northwestern British Columbia (Figure 1). The Yukon-Tanana terrane rifted from the Laurentian margin in the Late Devonian and underwent several cycles of magmatism, deformation and metamorphism prior to returning to the Laurentian margin [12, 13, 15, 16]. Late Devonian and older siliciclastic rocks of the Snowcap assemblage are the oldest rocks in the Yukon-Tanana terrane and form the basement to Paleozoic and Mesozoic magmatic rocks [12, 14, 17, 18].

The Yukon-Tanana terrane records several episodes of eclogite facies, Barrovian and low-pressure contact metamorphism over ca. 250 Myr (Table 1; [7, 19–29]). These metamorphic episodes have been used to support often-contradictory models that attempt to reconcile the collisional history of Yukon-Tanana terrane with adjacent Stikine, Slide Mountain and Cache Creek terranes and/or with the Laurentian margin (e.g., [12, 24, 25, 30–35]). The distinct metamorphic conditions experienced by different localities at different times, however, imply that the Yukon-Tanana terrane is likely a composite terrane or records a heterogeneous tectonic evolution (e.g., [32, 33, 36, 37]). The Yukon-Tanana terrane of northwestern British Columbia represents an important piece in the northern Cordilleran tectonic puzzle. There, it has been interpreted as a stratigraphic equivalent to the basement of the Stikine terrane (e.g., [38]), which is a long-lived (Early Devonian-Middle Jurassic) composite island arc terrane [39], or as a separate terrane that was tectonically juxtaposed with the Stikine terrane in the Late Triassic [35, 40] to early Middle Jurassic [41]. Receding glaciers in the Atlin Lake area expose low variance metamorphic assemblages in metapelitic rocks, including coexisting metastable kyanite, sillimanite, and andalusite. As such, this area provides the opportunity to constrain the metamorphic evolution of the Yukon-Tanana terrane exposed in northwestern British Columbia and to compare this evolution to the existing tectonic framework (Figure 1).

The Yukon-Tanana terrane in northwestern British Columbia is separated from the Stikine terrane to the east by the Llewellyn fault zone and is intruded by the Cretaceous to Eocene Coast Plutonic Complex in the west (Figures 1 and 2). The Llewellyn fault zone is a subvertical, brittle-ductile deformation zone that was episodically active from the Cretaceous to the Eocene [42, 43], though older, Triassic motion has been inferred [8]. The Yukon-Tanana terrane in the Atlin Lake area has been subdivided into the Florence Range and the Boundary Ranges metamorphic suites (Figure 2; see review in [8]). These two metamorphic suites are separated by the folded, 1 to 4 km thick Wann River shear zone [7, 41]. The Wann River shear zone is subvertical near Llewellyn glacier, north dipping in the Wann River headwaters, and subhorizontal on White Moose Mountain where a small klippe of Florence Range metamorphic rocks is preserved in the hanging wall (Figure 2; see also cross-sections in Figure 3 of Currie and Parrish [44]). This shear zone is interpreted to have been active between 185 and 170 Ma based on sheared Early Jurassic plutons and U–Pb titanite and rutile dates that were interpreted as postmetamorphic cooling ages [7, 41, 45]. A common retrograde metamorphic overprint (e.g., garnet replaced by chlorite and biotite) hindered accurate pressure-temperature estimates using conventional geothermobarometry in the Florence Range and Boundary Ranges suites [7]. Furthermore, abundant Cretaceous to Eocene plutons locally overprinted metamorphic assemblages related to Paleozoic and early Mesozoic tectonism and metamorphism [8].

The Florence Range suite [46, 47] comprises biotite-muscovite quartzofeldspathic schist and gneiss with minor layers of metapelite, quartzite, marble, calc-silicate gneiss, and amphibolite. These were interpreted to represent a metamorphosed passive margin sedimentary sequence [7] and correlated to the Snowcap assemblage of the Yukon-Tanana terrane in Yukon [18]. Previous workers inferred two different sequences of crystallization of Al₂SiO₅ polymorphs based on petrographic observations [7, 8]. Currie [7] suggested an early medium T and low P event in the andalusite stability field, followed by a medium T and medium P event in the kyanite stability field. A third high T and medium P suprasolidus event in the sillimanite stability field was inferred from Early Jurassic multigrain discordant U–Pb ID-TIMS dates on zircon from a leucosome and monazite from a garnet-muscovite-biotite schist (⁠$177±7Maand177.5±1.5Ma$⁠, respectively). In contrast, Mihalynuk et al. [8] suggested a kyanite to andalusite to sillimanite crystallization sequence without absolute timing constraints.

The Boundary Ranges suite [48–50] is composed of muscovite-$biotite±garnetquartzofeldspathic$ schist, chlorite-$actinolite±garnet$ schist, and rare graphitic schist. The provenance of the Boundary Ranges suite is ambiguous. Protoliths have been interpreted to represent island arc volcano-sedimentary rocks correlative to the basement of the Stikine terrane (Stikine Assemblage; [44]), or to the Yukon-Tanana terrane [8]. In general, metamorphism in the Boundary Ranges suite occurred at lower temperature than in the Florence Range suite, though rare andalusite and kyanite have been reported [8]. The age of metamorphism in the Boundary Ranges suite is unconstrained but postdates the intrusion of an Early Mississippian granitic orthogneiss [44].

Several suites of variably metamorphosed and deformed felsic to intermediate intrusive rocks occur in the area (Figure 2). The Late Devonian Bighorn Creek orthogneiss (⁠$366±9Ma$ igneous crystallization age, U–Pb zircon, [44]) intrudes the Boundary Ranges suite northwest of White Moose Mountain. The Wann River gneiss (⁠$270±5Ma$ igneous crystallization age, U–Pb zircon; [45]) is a pervasively sheared hornblende-plagioclase banded orthogneiss that is tectonically interlayered at the map and outcrop scales with the Florence Range suite within the Wann River shear zone [44]. The variably deformed Early Jurassic Tagish Lake suite, including the Hale Mountain Granodiorite (⁠$187±1.1Ma$ igneous crystallization age, U–Pb zircon, [7, 45]) and the Mt. Caplice quartz monzodiorite (⁠$180.9±3.1Ma$ igneous crystallization age; [44]), intrudes the Boundary Ranges suite and the Wann River shear zone. The Tagish Lake suite has been correlated to the Long Lake plutonic suite in Yukon ([8, 51]). Cretaceous to Eocene felsic to intermediate plutonic and volcanic rocks of the Coast Plutonic Complex and Sloko suites form an overlap assemblage in the Florence Range and Boundary Ranges suites (130–55 Ma igneous crystallization ages; [8, 52] and references therein).

Isochemical phase diagrams were calculated using THERIAK-DOMINO [53] and thermodynamic data of Holland and Powell [54], with solution models as in Gaidies et al. [55] and silicate melt models from White et al. [56]. Calculations were done in the MnO-Na₂O-CaO-K₂O-FeO-MgO-Al₂O₃-SiO₂-H₂O (MnNCKFMASH) system. Rutile and ilmenite are excluded from the calculated diagrams because the position of their field of stability is highly dependent on the chemical composition of the system and fO₂ and could not be reliably constrained in our interpretations. Excluding these accessory phases does not have a significant effect on the stability fields of the major phases used in our interpretations but simplifies phase diagrams and improves readability (e.g., [57–60]). TiO₂ and P₂O₅ were thus omitted from the input chemical composition, such that THERIAK-DOMINO recalculated proportions of other elements. H₂O is considered in excess at subsolidus conditions, whereas we used a fixed H₂O proportion at suprasolidus conditions based on the amount of water necessary in each sample to saturate the system just below the solidus at 600 MPa. We did not assess the effects of varying fO₂ because there is no evidence of elevated fO₂, such as the presence of magnetite, hematite and epidote in the low variance samples (25B, 31, 46B, and 49B) that are critical to constrain P–T conditions. The chemical composition used for phase equilibrium modelling was determined with X-ray fluorescence (XRF) on fist-sized sample fragments free of obvious alteration. Samples were crushed, split, pulverized with mild steel, fused using lithium metaborate/tetraborate in Pt crucibles, and analyzed on a Panalytical Axios Advanced wavelength dispersive XRF at ActLabs (Ancaster, Ontario, Canada).

Thin sections were scanned with a Tescan Mira3 Field Emission Scanning Electron Microscope (SEM) equipped with an X-MAX 80 Silicon Drift Detector energy dispersive spectrometer (EDS) at the Geological Survey of Canada (Ottawa, Canada). Automated backscattered electron (BSE) Feature Identification (Oxford Instruments AZtec 4.2) of grains with a BSE brightness higher than that of garnet and manual EDS were used to detect monazite and xenotime, identify matrix minerals and inclusions in porphyroblasts, and complement petrographic observations of mineral textures. Monazite and xenotime as small as 5 μm are detected with this identification routine, ensuring that every exposed crystal large enough to be dated was identified. Monazite and xenotime used for further chemical characterization and petrochronology (see following) were selected based on size, texture, lack of cracks and inclusions, and textural relationship to metamorphic minerals.

Selected monazite and xenotime were chemically mapped for Y, U, Th, Ca, and Si (Y, U, Th, Dy, and Si for xenotime) at the Geological Survey of Canada (Ottawa) on a JEOL JXA-8230 electron microprobe equipped with 5 vertical wavelength-dispersive spectrometers running ProbeforEPMA software. X-ray maps of monazite and xenotime were used to identify zonation, detect inclusions, and select laser spot locations. One garnet grain was chemically mapped for Ca, Mn, and Mg with the same instrument to assess chemical zoning. For all X-ray maps, electron microprobe conditions were as follows: an acceleration voltage of 20 kV, beam current of 200 nA, and dwell time of 100 ms. Step size was 0.1-0.5 μm for monazite and xenotime maps depending on grain size and 10 μm for garnet maps.

U-Th-Pb isotopes and trace element concentrations were measured using an 8 μm diameter spot by laser-ablation split stream inductively-coupled plasma mass spectrometry at the University of California (Santa Barbara). U-Th-Pb isotopic ratios were measured on a Nu Plasma 3D multicollector ICP-MS, and trace element concentrations (U, Th, Pb, Ca, Y, and La–Lu except Pm) were measured on an Agilent 7700S quadrupole ICP-MS. Analytical methods are described by Kylander-Clark et al. [61] with modifications as outlined by McKinney et al. [62]. Mass bias and Pb/U and Pb/Th downhole fractionation were monitored and corrected for using a primary reference monazite, “44069” (424 Ma, ²⁰⁶Pb/²³⁸U isotope dilution-thermal ionisation mass spectrometry age (ID-TIMS) [63]). Normalization of xenotime unknowns to “44069” monazite primary reference material does not cause any resolvable difference in accuracy and is therefore appropriate [62]. A common Pb correction was not applied because of the general absence of measurable ²⁰⁴Pb in monazite and xenotime [61, 62]. The relatively small uncertainty in date due to ²⁰⁴Pb is therefore incorporated in the external reproducibility of the secondary reference materials. Two secondary reference monazites “FC-1” (55.7 Ma ²⁰⁶Pb/²³⁸U ID-TIMS age [64]) and “Trebilcock” (273 Ma ²⁰⁶Pb/²³⁸U ID-TIMS age [65]) were analyzed repeatedly throughout the analytical session and treated as unknowns to monitor isotopic data accuracy. Trace element concentrations were normalized to an in-house reference “Bananeira” monazite and, based on the long-term reproducibility of secondary reference monazites, are accurate to 3–5%. During the analytical period, repeat analyses of FC-1 gave a weighted mean ²⁰⁶Pb/²³⁸U date of $55.7±0.5$ Ma and $meansquareofweighteddeviatesMSWD=1.2n=16$⁠, a weighted mean ²⁰⁷Pb/²³⁵U date of $55.7±0.9MaandMSWD=2.4n=16$⁠, and a weighted mean ²⁰⁸Pb/²³²Th date of $57.8±1.4MaandMSWD=5.4n=16$⁠. Repeat analyses of Trebilcock gave a weighted mean ²⁰⁶Pb/²³⁸U date of $276.5±2.7MaandMSWD=1.3n=16$⁠, a weighted mean ²⁰⁷Pb/²³⁵U date of $274.4±2.5MaandMSWD=0.61n=16$⁠, and a weighted mean ²⁰⁸Pb/²³²Th date of $260.6±3.4MaandMSWD=2.3n=16$⁠. Data reduction, including corrections for baseline, instrumental drift, mass bias, downhole fractionation, and uncorrected date calculations were carried out using Iolite version 2.5 [66]. Diagrams were constructed with Igpet 2018 and Microsoft Excel 2016. Uncertainties include external reproducibility of primary reference materials for ²⁰⁷Pb/²⁰⁶Pb, ²⁰⁶Pb/²³⁸U, and ²⁰⁸Pb/²³²Th ratios and are quoted at $2σ$⁠.

Several precautions were taken to ensure that individual dates represent unique and homogeneous age domains. Time-resolved isotopic signals, chemical maps, and photomicrographs were used to identify analyses that may have sampled multiple chemical or age domains, inclusions, microfractures, or matrix material during ablation in the pit (3-4 μm depth). Trace element content of each spot was examined to detect the influence of inclusions (e.g., anomalous high Y from xenotime in monazite, Ca from apatite or allanite, and Th from thorite). Analyses affected by the abovementioned criteria were excluded from geochronological interpretations. Analyses with >5% uncertainty on ²⁰⁷Pb/²³⁵U or ²⁰⁶Pb/²³⁸U were also excluded from the interpretations. In total, 186 of the 841 analyses (~22%) were rejected for analytical reasons. All remaining analyses used for interpretations are <15% discordant (²⁰⁷Pb/²³⁵U vs. ²⁰⁶Pb/²³⁸U dates). From the 655 analyses included in geochronological interpretations, 446 of the 476 monazite dates are <5% discordant (²⁰⁸Pb/²³²Th vs. ²⁰⁶Pb/²³⁸U dates) and 145 of the 179 xenotime dates are <5% discordant (²⁰⁷Pb/²³⁵U vs. ²⁰⁶Pb/²³⁸U dates). Age interpretations are based on the ²⁰⁶Pb/²³⁸U monazite and xenotime dates because of low accumulated radiogenic ²⁰⁷Pb. The terms “date” and “age” refer to results and their interpretation in a tectonometamorphic context, respectively.

Thirty-seven metapelite samples were collected during targeted fieldwork in the Atlin Lake area (Figure 2 and Figure S1) to constrain the tectonometamorphic evolution of the area. Samples were collected away from exposed rocks of the Coast Plutonic Complex to limit, as much as possible, the overprinting effect of Cretaceous to Eocene contact metamorphism. Samples selected for petrochronology were chosen based on their metamorphic assemblages, prioritizing low variance samples with index minerals such as Al₂SiO₅ polymorphs, cordierite, K-feldspar, or staurolite that are more likely to reveal information about P–T conditions of metamorphism. The Boundary Ranges suite was not dated because none of the thirteen samples scanned with the SEM AZtec Feature tool contain monazite or xenotime sufficiently large to be analyzed with laser-ablation split-stream (>10 μm). The absence of monazite and xenotime in samples from the Boundary Ranges suite could be due to unfavourable bulk compositions and/or low metamorphic P–T conditions. The absence of medium to high P–T metamorphic index minerals other than garnet (all samples) and andalusite (one sample) suggests that the metamorphic conditions were too low for monazite and xenotime crystallization in the samples from the Boundary Ranges suite (e.g., [67, 68]).

Florence Range suite samples contain abundant monazite and xenotime to support petrochronological investigations. Six samples containing Al₂SiO₅ polymorphs, cordierite, K-feldspar, or staurolite that encompass the entire range of metamorphic index minerals observed in the Florence Range suite were selected for petrochronology (samples 25B, 31, 42, 46B, 83, and 49B; see Figure 2 and Figure S1). As the petrochronological record of monazite and xenotime can be incomplete in rock samples that share a common P–T–t evolution (e.g., [69]), two additional samples lacking index minerals were analyzed to further support petrochronological results (samples 67 and 77A). Overall, the eight analyzed samples cover a wide geographic and structural area and are representative of all metamorphic assemblages we and previous workers (e.g., [7, 8]) identified in the Florence Range suite. Samples 25B, 31, and 67 are from Willison Glacier area; sample 42 is from White Moose Mountain; and samples 46B, 83, 49B, and 77A are from the Eastern Wann River shear zone (Figure 2, Table 2). A summary of petrographic observations critical for interpretations is reported herein, and detailed descriptions of individual samples are provided in the supplementary Text S1. Mineral occurrences in each sample are summarized in Table 2.

All samples are coarse-grained, mica-rich, porphyroblastic gneisses with a well-developed foliation. Leucocratic (⁠$Qz+Pl±Kfs±Grt$⁠) segregations ranging from a few millimeters to a few centimeters are common in metapelitic outcrops (Figure 3(a)). They are interpreted as in situ leucosome based on their mineralogical composition, which is consistent with the presence of index metamorphic minerals characteristic of upper amphibolite facies metamorphism, such as prismatic sillimanite. Microtextural evidence indicates the former presence of partial melt, such as euhedral peritectic phases (e.g., Kfs) and biotite-quartz-$plagioclase±sillimanite$ symplectites potentially indicating a reaction between a peritectic phase and cooling melt (samples 46B, 25B, and 49B below) [70, 71]. Delicate features such as cuspate volumes and thin films of quartz around refractory minerals, veinlets of quartz and plagioclase, or microgranophyric intergrowths of quartz and alkali feldspar (e.g., [70–72]) are difficult to recognize unequivocally. These may have been obliterated by recrystallization, either during deformation along the Wann River shear zone, where half of our samples were collected, or during subsequent, lower-grade metamorphism.

In most samples, garnet is partially replaced by coronae of $Pl+Bt+Fe‐sulfide±whitemicaWm±Sil±And±Crd±Chl$⁠. In the replacement corona around garnet in sample 25B, cordierite occurs in the inner corona and andalusite in the outer corona (Figure 3(b)). In sample 31, circular patches of $Bt+Sil+Qz+Pl+Ilm+Fe$-sulfide are likely pseudomorphs after garnet, which has been completely replaced (Figure 3(c)). Garnet in Florence Range samples is typically characterized by gradually decreasing Ca and Mn from core to rim, consistent with a single phase of garnet growth (Figure 4(a); see also Dyer [73]). In contrast, garnet in samples 49B and 83 is characterized by distinct inclusion assemblages and/or chemical compositions in the cores and the rims. Staurolite, and more rarely kyanite, are only observed as inclusions in garnet in sample 49B (Figure 3(d)). Garnet in sample 49B is characterized by distinct compositions in the core and the rim separated by a sharp discontinuity that points to multiple episodes of garnet growth (a central section of the largest garnet of the sample was chemically mapped in another study by Dyer [73]; Figures 4(b) and 4(c)). Although the chemical composition of garnet grains in the thin section used for monazite and xenotime petrochronology was not characterized, the chemical zoning is expected to be similar to the garnet mapped in Dyer [73] on the basis of the unique occurrence of staurolite inclusions in all garnet cores. In sample 83, garnet is characterized by increasing Mg and decreasing Mn and Ca from core to rim, except for several Ca-rich patches separated from Ca-poor garnet by sharp discontinuities locally present at the rim (Figures 4(d) and 4(e) and Figure S2). This garnet is resorbed significantly; therefore, it is difficult to assess the preresorption distribution of the Ca-rich areas. We interpret the Ca-rich garnet patches at the rim to represent a second phase of garnet growth, as in sample 49B that shows a similar but better-preserved Ca-rich rim.

Kyanite forms foliation-parallel porphyroblasts partially replaced by coarse flakes of white mica in samples 31 and 46B and by fine grained white mica (margarite and muscovite) in sample 42 (Figure 3(e)), whereas it only occurs as small inclusions in garnet in sample 49B (Figure 3(d)). Where present, sillimanite occurs as foliation-parallel fibrolite, rare prisms, and as randomly oriented mats in garnet strain shadows and in coronae around resorbed garnet (Figure 3(f)). Andalusite is only present in samples from the toe of Willison Glacier (samples 25B and 31, Figure 2). There, andalusite forms randomly oriented euhedral to subhedral porphyroblasts that overprint the matrix foliation and locally overgrew kyanite (Figures 3(f) and 3(g). Andalusite is not observed in contact with sillimanite. Anhedral cordierite poikiloblasts overprint the sillimanite, biotite and white mica that define the foliation in samples 25B and 49B (Figure 3(g)). In sample 25B, cordierite commonly encloses andalusite crystals (Figure 5(a)). Anhedral to euhedral K-feldspar, identified with EDS-SEM, is rare and <1 mm in diameter in samples 46B and 83 (Figure 3(h)). Allanite is only abundant in sample 67, in which it is commonly rimmed by epidote in the matrix and forms inclusions in garnet. In other samples, allanite is absent or occurs as $rare<100μm$ in diameter grains in the matrix or as inclusions in garnet cores and rims. All samples from the Eastern Wann River shear zones (samples 46B, 83, 49B, and 77A) are characterized by a sinistral C-S fabric defined by white mica, biotite, and sillimanite (Figure 3(i)).

Variably abundant monazite and xenotime are present in all selected samples. Monazite is present in the matrix of all samples and commonly forms inclusions in porphyroblasts. Monazite occurs as single inclusions in the rim of garnet in samples 46B and 77A (Figure 5(b)) and in the cores of garnet in samples 83 and 49B (Figures 5(c) and 4(e)). Monazite is also part of polymineralic inclusions in garnet rims of samples 42 and 77A and rarely occurs in coronae around garnet in samples 25B and 46B. Monazite inclusions in kyanite occur in samples 31 and 42, but these inclusions intersect thin cracks and are separated from kyanite by thin (<5 μm) rims of white mica (Figure 3(e)). Monazite in textural contact and intergrown with sillimanite is common in samples 31, 46B, and 49B (Figure 5(d)). Cordierite and andalusite contain inclusions of monazite in samples 25B and 31 (Figures 5(a) and 5(e)). Xenotime is typically found in the matrix, in particular in coronae around resorbed garnet (Figure 3(c)). In sample 83, xenotime occurs as a single inclusion near the edge of a partially resorbed garnet porphyroblast. In sample 49B, inclusion-rich xenotime with irregular boundaries is intergrown with sillimanite. In sample 25B, xenotime commonly forms inclusions in cordierite and andalusite (Figures 5(e) and 5(f)).

The samples investigated herein are not ideal to determine precise P–T conditions of metamorphism with phase equilibrium modelling, in particular with composition and volume isopleth thermobarometry, because they contain evidence for polymetamorphism, suprasolidus conditions, and zoned minerals. As such, evaluation of the effects of element sequestration in zoned minerals (e.g., [74]), melt loss (e.g., [75]), and metasomatism (e.g., [76]) would be required to assess the chemical composition of the equilibration volume during each metamorphic event to obtain reliable quantitative P–T constraints from mineral composition or volume isopleths. Nonetheless, phase equilibrium modelling is useful to illustrate the predicted stability in P–T space of metamorphic index minerals, which variably depends on the chemical composition of the rocks. Then, the least variable mineral stability fields can be used to estimate P–T conditions semiquantitatively during monazite and xenotime crystallizing events. We calculated isochemical phase diagrams using the analyzed bulk-rock chemical composition of all dated samples (Table S1) to provide a petrogenetic framework encompassing a wide range of geochemical compositions (Figure 6(a)). Most samples plot within or close to the average field of low-Al pelite composition, but sample 31 plots near the field of high-Al pelites of Spear [77].

Individual phase diagrams calculated with the composition of each dated sample are included in Figure S3, whereas Figure 6(b) illustrates the minimum and maximum extents of the stability fields of metamorphic index minerals predicted for all samples. The positions of the Al₂SiO₅ polymorph transitions do not depend on the chemical composition of the model system given that these polymorphs hardly act as solid solution phases. The maximum extent of the predicted stability fields of Al₂SiO₅ polymorphs is consistent with constraints from experiments, thermodynamic modeling, and natural mineral assemblages (e.g., [78]). As such, the variable presence and textural relations among the Al₂SiO₅ polymorphs are especially useful guides to the metamorphic evolution of the rocks. Accordingly, kyanite indicates conditions over 0.6 GPa at 600°C and over 1.0 GPa at 800°C. Andalusite indicates conditions between 500 and 700°C below 0.4 GPa. Sillimanite crystallizes above 550°C on the low-P side of the Ky–Sil transition and on the high-P side of the And–Sil transition. In our calculations, cordierite is indicative of high T (>600°C) and low to medium P (<0.5 GPa) metamorphism. K-feldspar is predicted to be stable across the range of pressure considered (0.2 to 1.0 GPa) above 700°C but is restricted to $P<0.7$ GPa between 600 and 700°C. However, the extent of the stability fields of cordierite and K-feldspar varies moderately with the chemical composition of the model system. The predicted stability field of staurolite is very sensitive to the chemical composition of the model system and indicates moderate T (~500–625°C) but cannot be used confidently to constrain P.

Monazite and xenotime displaying diverse textural relationships with metamorphic index minerals were selected for in situ laser-ablation split-stream U-Th-Pb and trace-element analysis. Wherever possible, inclusions in garnet, kyanite, andalusite, and cordierite or grains associated with metamorphic reactions (e.g., garnet partial replacement) were analyzed. Dates and trace element compositions are presented in Figures 7 and 8, and selected examples of monazite and xenotime grains with date spots are presented in Figure 9. Complete results are presented in Table S2, and locations of monazite and xenotime in thin section, BSE images, and compositional X-ray maps are presented in Figure S4.

Monazite and xenotime dates in sample 25B range from $187±7to161±6Ma$ and $177±6Mato121±4Ma$⁠, respectively (Figure 7). Monazite inclusions in cordierite yield dates between 187 and 161 Ma (21 spots on 4 grains, $n=21$⁠; 4, Figure 9(a)). Monazite inclusions in andalusite yield dates between 177 and 161 Ma (⁠$n=11$⁠; 3). Monazite in the matrix yields dates between 185 and 179 Ma (⁠$n=26$⁠; 3). A monazite grain located within a corona around garnet yields dates between 178 and 172 Ma (⁠$n=5$⁠; 1). Xenotime inclusions in cordierite yield dates between 169 and 121 Ma (⁠$n=15$⁠; 3, Figure 9(b)). A xenotime inclusion in andalusite yields dates between 170 and 145 Ma (⁠$n=7$⁠; 1). Xenotime in the matrix yields dates between 177 and 141 Ma (⁠$n=13$⁠; 3). Xenotime located within a corona around garnet yields dates between 172 and 151 Ma (⁠$n=9$⁠; 2). Yttrium and Gd/Yb in monazite and xenotime do not vary systematically with date. Two monazite grains (one inclusion in cordierite, one in a corona around garnet) have lower Y and higher Gd/Yb compared to other monazite of the same age.

Monazite and xenotime dates in sample 31 range from $185±6to164±5Ma$ and $185±7Mato157±5Ma$⁠, respectively, with two outliers at $522±17$ Ma (monazite) and $228±8$ Ma (xenotime) (Figure 7). A monazite inclusion in kyanite yields dates between 183 and 171 Ma (⁠$n=7$⁠; 1), whereas an eighth spot from the core of the grain yields a date of 522 Ma. Monazite inclusions in andalusite yield dates between 179 and 164 Ma (⁠$n=15$⁠; 3, Figure 9(c)). Monazite in the matrix yields dates between 185 and 169 Ma (⁠$n=26$⁠; 4). Monazite in textural contact and more rarely intergrown with sillimanite in pseudomorphed garnet yields dates between 179 and 167 Ma (⁠$n=22$⁠; 5). Xenotime located around pseudomorphed garnet yields dates between 185 and 157 Ma (⁠$n=51$⁠; 5), and one spot from the core of a 100 μm large xenotime grain yields a 228 Ma date. Yttrium and Gd/Yb in monazite do not vary with date. However, there is a variation with position in the thin section: three monazite grains (one included in kyanite, one included in andalusite, and one in the matrix) from the same compositional layer away from pseudomorphed garnet have lower Y and higher Gd/Yb compared to monazite within or near pseudomorphed garnet, suggesting local variations in the trace element budget (Figure S5). Yttrium and Gd/Yb in xenotime do not vary with date.

Monazite and xenotime dates in sample 67 range from $182±6$ to $168±6$ Ma and $181±6$ Ma to $166±6$ Ma, respectively (Figure 7). Monazite in the matrix yields dates between 182 and 168 Ma (⁠$n=11$⁠; 4). Xenotime in the matrix yields dates between 179 and 166 Ma (⁠$n=7$⁠; 2). Xenotime located within a corona around garnet yields dates between 181 and 166 Ma (⁠$n=12$⁠; 5). Yttrium and Gd/Yb do not vary with date or textural position.

Monazite and xenotime dates in sample 42 range from $187±6to169±6Ma$ and $171±6Mato168±6Ma$⁠, respectively (Figure 7). One monazite grain at the edge of a millimeter-sized inclusion of quartz and chloritized biotite within garnet yields two dates at 185 and 182 Ma. Monazite inclusions in kyanite yield dates between 186 and 169 Ma (⁠$n=24$⁠; 3, Figure 9(d)). Monazite in the matrix yields dates between 187 and 175 Ma (⁠$n=39$⁠; 5). Two spots on one xenotime grain located in a corona around garnet yield dates of 171 and 168 Ma. In monazite, Y and Gd/Yb display a weak increase and decrease, respectively, with decreasing date, especially for $spots<180$ Ma. Monazite inclusions in garnet and kyanite have higher Y and lower Gd/Yb compared to contemporaneous monazite in the matrix.

Monazite and xenotime dates in sample 46B range from $191±7to175±6Ma$ and $183±7Mato170±6Ma$⁠, respectively (Figure 8). Monazite inclusions in garnet rims yield dates between 191 and 184 Ma (⁠$n=14$⁠; 3, Figure 9(e)). Monazite in the matrix yields dates between 189 and 176 Ma (⁠$n=28$⁠; 3). Monazite in textural contact with sillimanite yields dates between 188 and 175 Ma (⁠$n=24$⁠; 3). A monazite grain located within a corona around garnet yields dates between 184 and 181 Ma (⁠$n=4$⁠; 1). Xenotime located within a corona around garnet yields dates between 183 and 170 Ma (⁠$n=12$⁠; 6). Yttrium is lower and Gd/Yb higher in 191-184 Ma monazite included in garnet and in the core of a monazite in textural contact with sillimanite compared to younger monazite. $Monazite<184$ Ma has variable, but generally higher Y and lower Gd/Yb. Yttrium and Gd/Yb in xenotime do not vary with date or textural context.

Monazite and xenotime dates in sample 83 range from $182±6to171±6Ma$ and $180±6Mato168±7Ma$⁠, respectively, with two outliers at $253±9Ma$ (monazite) and $243±8Ma$ (xenotime) (Figure 8). A monazite grain fully included in the core of a garnet porphyroblast yields a 253 Ma date (⁠$n=1$⁠; 1; Figure 4(e)), whereas a monazite partly included in the rim of a second garnet grain next to an embayment yields two dates at 179 and 178 Ma (⁠$n=2$⁠; 1). Monazite in the matrix yields dates between 182 and 171 Ma (⁠$n=48$⁠; 5). A xenotime inclusion near the edge of a third partially resorbed garnet porphyroblast and, therefore originally located in the core of garnet, yields a date of 243 Ma. A xenotime grain in the matrix yields a 170 Ma date, and xenotime located within a corona around garnet yields dates between 180 and 168 Ma (⁠$n=8$⁠; 4). Yttrium and Gd/Yb in monazite and xenotime are homogenous and do not vary with date, except for the 250-240 Ma monazite and xenotime inclusions in garnet that have higher Y and lower Gd/Yb compared to <185 Ma grains in the matrix. Two spots on a monazite partly included in garnet have higher Y and similar Gd/Yb compared to contemporaneous monazite in the matrix.

Monazite and xenotime dates in sample 49B range from $270±9to167±5Ma$ and $186±7Mato170±5Ma$⁠, respectively (Figure 8). Monazite inclusions were analyzed in half of the garnet porphyroblast cores (3 of 6) exposed in the thin section and yield dates between 270 and 217 Ma (⁠$n=18$⁠; 5, Figure 9(f)). Monazite in the matrix yields dates between 190 and 175 Ma (⁠$n=40$⁠; 3). Monazite in textural contact with sillimanite yields dates between 188 and 167 Ma (⁠$n=20$⁠; 3). Xenotime in the matrix yields dates between 185 and 173 Ma (⁠$n=5$⁠; 2). Xenotime intergrown with sillimanite yields dates between 182 and 170 Ma (⁠$n=21$⁠; 4, Figure 9(g)). Xenotime located within a corona around garnet yields dates between 186 and 174 Ma (⁠$n=7$⁠; 3). Yttrium and Gd/Yb in monazite included in garnet vary with date and grain; four monazite inclusions in two garnet grains have low Y and are generally older (270-243 Ma) than a fifth inclusion in a third garnet grain that is characterized by high Y and younger dates (244-217 Ma). No clear variation in Gd/Yb is observed in monazite included in garnet. Monazite outside of garnet shows a faint increase in Y and decrease in Gd/Yb in $spots<185$ Ma compared to the older dates. Yttrium and Gd/Yb in xenotime do not vary with date or textural context.

Monazite and xenotime dates in sample 77A range from $192±7to176±6Ma$ and $181±6Mato172±6Ma$⁠, respectively (Figure 8). Monazite in polymineralic inclusions in garnet rims yields dates between 184 and 178 Ma (⁠$n=4$⁠; 2). Monazite in the matrix yields dates between 192 and 176 Ma (⁠$n=62$⁠; 7). Xenotime located within a corona around garnet yields dates between 181 and 172 Ma (⁠$n=8$⁠; 5). Yttrium and Gd/Yb in monazite are higher and lower, respectively, in $spots<185$ Ma compared to $spots>185$ Ma (Figures 9(h) and 9(i)), irrespective of textural position (inclusion in garnet or matrix). Yttrium and Gd/Yb in xenotime do not vary with date or textural context.

Textural characterization of metamorphic index minerals is useful to identify metamorphic events and reconstruct their relative timing. In addition, several textural observations indicate that monazite and xenotime are likely to provide robust constraints on the absolute timing of metamorphic mineral growth and resorption. In particular, monazite inclusions in the two garnet populations characterized by different inclusions assemblages and chemical compositions are expected to constrain the timing of these two phases of garnet growth (Figures 3(d), 4, 5(b), 5(c), 9(e), and 9(f)). Monazite and xenotime may also provide constraints on the age of crystallization (or maximum age for petrochronometer inclusions) of other metamorphic index minerals such as sillimanite, with which monazite and xenotime are commonly in textural contact or intergrown (Figures 5(d) and 9(g)), or andalusite and cordierite that engulfed several inclusions of the two petrochronometers (Figures 5(a), 5(e), 5(f), 9(a), 9(b), and 9(c)). Monazite and xenotime petrochronology also has potential to date mineral resorption; coronae around resorbed garnet commonly contain abundant xenotime likely associated with this reaction (Figure 3(c)). Section 7.1 of our discussion explains how monazite and xenotime data are linked to metamorphic textures to provide absolute age constraints on the metamorphic evolution of the Florence Range suite.

Acquisition of trace element data and U-Th-Pb isotope ratios combined with petrographic observations and textural context interpretation (petrochronology) enables direct correlation between dates, chemical composition, and growth conditions (e.g., [79–83]). In the current study, there are clear links between dates and textural context, which we use primarily to interpret the monazite and xenotime dates. In addition, variations in trace element contents in some samples are used to support our interpretations. The Y and HREE contents of monazite, xenotime, and garnet are controlled by their crystallization and resorption histories. In general, monazite grown during or after garnet or xenotime growth yields high Gd/Yb (an inverse proxy for HREE abundance) and low Y content, whereas monazite grown during or after garnet or xenotime breakdown should have a low Gd/Yb and high Y content [67, 84–86]. In metapelitic rocks in general and in our samples specifically, garnet grows with increasing temperature and/or pressure (except during staurolite crystallization) and may be consumed during decompression, cooling, or melt crystallization (Figure 6(b) and Figure S3), whereas xenotime reacts out during garnet crystallization and may be produced following garnet breakdown [67]. Trace element patterns in monazite must be carefully evaluated along with microtextural context to attribute trace element content variations to the appropriate metamorphic reactions, especially those involving garnet and/or xenotime. The Y and HREE contents in monazite coexisting with xenotime may also increase with temperature of crystallization (e.g., [87]). However, the Y and HREE trends identified in our samples correlate with changes in textural settings involving garnet growth or consumption, suggesting that the stability of this phase exerted a dominant control on Y and HREE distribution in monazite. The links between dates, textural context, and trace element compositions of monazite in the samples investigated herein indicate that the trace element data can be used to assist reconstructing the metamorphic reaction history and evaluate the effects of local variations in trace element availability.

A single spot on a Th-rich monazite core from sample 31 in Willison Glacier area yields a $522±17$ Ma date. We interpret this anomalously old date as a potential inherited detrital core, which would provide a maximum age of deposition for the Florence Range suite. Although monazite typically recrystallizes during greenschist to lower amphibolite prograde metamorphism (e.g., [88, 89]), rare detrital monazite cores or inherited metamorphic grains can be preserved up to upper amphibolite facies (e.g., [90–93]). A Cambrian maximum age of deposition is consistent with other age constraints from the correlative Devonian and older Snowcap assemblage in Yukon and Alaska [23, 94] and considering that Late Devonian to Early Mississippian siliciclastic and volcanic rocks unconformably overlie the Snowcap assemblage [15, 95].

Chemical maps of garnet in samples 49B and 83 from the Florence Range suite indicate chemically distinct cores and rims, which points to two phases of garnet growth (Figure 4; [73]). All garnet cores in sample 49B are also characterized by the unique presence of staurolite inclusions (e.g., Figures 3(d) and 5(c)), suggesting that the mapped garnet is representative of other garnet grains in the sample. The monazite and xenotime dates from inclusions in garnet in samples 49B and 83 are also unique: they are older than monazite and xenotime from the matrix and from other samples. In sample 49B, monazite inclusions in three garnet porphyroblasts yield dates between $270±9$ Ma and $217±8$ Ma (Figure 8). These inclusions can be split into two groups based on their age, geochemical composition and textural context. The first four inclusions (270-243 Ma) have low Y and are fully enclosed in the core of two garnet grains (Figures 5(c) and 9(f) and Figure S4). These 270-243 Ma monazite inclusions are interpreted to have crystallized during garnet growth because the low Y content is indicative of growth synchronous with a high-Y phase such as garnet. The fifth inclusion (244-217 Ma) has higher Y, is located near the edge of another, half resorbed garnet, and intersects a large crack (Figure S4). We interpret this inclusion to have (re)crystallized following garnet partial breakdown. In sample 83, one monazite inclusion and one xenotime inclusion in two garnet grains yield additional Permian to Early Triassic dates (⁠$253±9$ Ma and $243±8$ Ma, respectively; Figure 8). The monazite inclusion (⁠$253±9$ Ma) is located in the chemically distinct core of the mapped garnet grain (Figures 4(d) and 4(e)), and both inclusions are interpreted to have grown before or during garnet growth.

The preservation of Permian-Triassic monazite in the garnet cores and their absence in the matrix of samples 49B and 83 suggest that all >215 Ma monazite grains not shielded by garnet were resorbed or recrystallized by the earliest Jurassic (<190 and<183 Ma monazite in the matrix in samples 49B and 83, respectively). The absence of early Jurassic monazite or xenotime inclusions in garnet, despite their abundance in the matrix (Figure 8), suggests that the growth of the garnet cores that contain Permian-Triassic monazite and xenotime inclusions occurred before the crystallization of Early Jurassic monazite. Therefore, Permian-Triassic monazite almost certainly could not have been included in garnet that crystallized during Early Jurassic metamorphism (see below). Another alternative, i.e., garnet growth between 215 and 190 Ma, fails to explain the clear presence of Permian-Triassic monazite, is inconsistent with the lack of 215-190 Ma monazite and xenotime, and would require a metamorphic event at a time during which metamorphism is not recognized in the Yukon-Tanana terrane (Table 1).

Garnet growth during the Permian to Middle Triassic (270–240 Ma) occurred along a prograde to peak segment of a P–T path that passed through the stability fields of staurolite and kyanite because these minerals are included in garnet cores of sample 49B (Figure 10(a)). This interpretation relies on the assumption that staurolite and kyanite are approximately the same age as the monazite inclusions and did not crystallize during an older metamorphic event for which there is no petrographic or geochronological evidence. The presence of a Y-rich monazite inclusion intersecting a crack near the edge of a half resorbed garnet in sample 49B (Figures 10(a) and S4) suggests that garnet breakdown likely occurred during retrograde metamorphism during the Late Triassic (240–215 Ma). This interpretation is also supported by a $228±8$ Ma date measured on the core of a 100 μm large xenotime grain in sample 31 (Figure 7), which we attribute to xenotime growth facilitated by Y liberation during garnet resorption.

Overall, evidence of Permian-Triassic metamorphism is rare in the Florence Range suite. With a single exception, Permian-Triassic monazite and xenotime are only preserved as inclusions in garnet, suggesting that they must have reacted out subsequently in the matrix. Shielding by garnet is critical to the preservation of Permian-Triassic metamorphic evidence, but it is not clear why Permian-Triassic garnet has not grown, or was not preserved, in every sample. One possibility is that the measured bulk composition of samples 49B and 83 is different from the original bulk composition, which could have been richer in MnO and have facilitated garnet growth (e.g., [96]). Alternatively, early Jurassic metamorphism (see below) may have largely obliterated the mineral assemblages of past metamorphic events, as has been recorded in other areas. For example, in the Canadian Cordillera, Vice et al. [97] analyzed monazite in a sample from the Snowcap assemblage that likely underwent early Jurassic metamorphism, based on the results of Clark [98] approximately 130 km to the north of their study area, but that only contains monazite formed during a Late Cretaceous kyanite-grade overprint. In the Himalaya, there is ample geochronological evidence that upper amphibolite facies metamorphism is the result of the Cenozoic India-Asia collision (e.g., [99, 100] and references therein), but detailed investigations have found evidence for a cryptic, earlier Cambro-Ordovician metamorphic event such as monazite and allanite inclusions in garnet [92, 93, 101, 102] or early Cambrian garnet porphyroblasts [103]. These studies demonstrate that an upper amphibolite facies metamorphic overprint can eliminate most textural and geochronological evidence of earlier garnet-grade or higher metamorphic events. In the current study, the systematic in situ analysis of petrochronometers preserved in metamorphic porphyroblasts was thus critical to obtain the maximum information about the P–T–t evolution of the metamorphic suite.

The earliest Jurassic marks a second episode of garnet growth during prograde metamorphism (Figure 10(a)). A garnet rim in sample 46B contains three low-Y and high Gd/Yb monazite inclusions that yield dates between $191±7Maand184±6Ma$ and constrain the maximum age of garnet rim growth to ca. 185 Ma (Figures 5(b), 8, and 9(e)). The timing of garnet core growth remains unconstrained. In addition, matrix monazite dated at 195-185 Ma in samples 49B and 77A is characterized by low-Y and high Gd/Yb, indicating crystallization postdating or synchronous with a Y- and HREE-rich mineral such as garnet (Figures 8 and 9(i)). The absence of >185 Ma monazite in other dated samples may be due to the stability of monazite relative to allanite. Allanite is preserved in the core and rim of garnet (samples 42 and 67; Figure S4) indicating that garnet grew under allanite-stable conditions in these samples (e.g., [68, 89, 104–106]). In sample 46B, monazite inclusions are preserved only in the rim of garnet (Figure 5(b) and Figure S4), whereas allanite inclusions are preserved in the core and rim, indicating that the garnet core grew under allanite-stable conditions and only the garnet rim grew under monazite-stable conditions.

A significant change in petrochronometer abundance, geochemical composition, and textural setting occurs at ca. 185 Ma, suggesting a change in metamorphic evolution. Approximately 75% of all monazite spots yield dates between 185-175 Ma and are generally characterized by higher Y and lower Gd/Yb compared to >185 Ma monazite (Figures 7, 8, 9(i), and 10(a)). Xenotime started to crystallize in the matrix and more commonly in coronae around resorbed garnet (Figure 3(c)) as early as ca. 185 Ma and became abundant at ca. 180 Ma (Figures 7 and 8). Garnet breakdown starting at ca. 185 Ma, supported by partial to complete resorption textures, likely increased the available Y and HREE budget, which resulted in a relative increase in Y and HREE in monazite and abundant xenotime crystallization (Figure 9(i)). Monazite and xenotime <185 Ma, and particularly <180 Ma, are commonly in textural contact or intergrown with sillimanite mats near resorbed garnet (Figures 5(d) and 9(g) and Figure S4), indicating that the Florence Range suite was in the sillimanite stability field when these monazite and xenotime grains crystallized. Because garnet is unlikely to break down significantly under increasing pressure or isobaric heating (except during staurolite crystallization, for which there is no evidence in the late Early Jurassic), the addition of sillimanite to the metamorphic assemblage at the expense of garnet is best explained by decompression across the sillimanite-in reaction that has a positive slope (Figure 10(b); see garnet volume isopleths in Figure 6(b) and Figure S3; e.g., [107, 108]). The absence of sillimanite inclusions in earliest Jurassic garnet further suggests that garnet did not grow in the presence of sillimanite, although the possibility that some of the sillimanite is prograde cannot be excluded. In addition, because garnet can break down during cooling (Figure 6(b) and Figure S3), it is possible that some of the 185–170 Ma monazite records a cooling segment of the P–T path (Figure 10(b)). Field and microtextural evidence for partial melting suggests that the rocks reached suprasolidus conditions (Figure 3(a)), whereas accessory K-feldspar in samples 46B and 83 (Figure 3(h)) indicates these rocks barely crossed the K-feldspar-in reaction (Figure 6(b) and Figure S3). Metamorphic reactions consistent with decompression and/or cooling thus indicate that rocks in the Florence Range suite may have started exhuming during the late Early Jurassic. Sinistral deformation (present day orientation) along the Wann River shear zone occurred during or after sillimanite crystallization (Figure 3(i)) bracketed by associated monazite growth at 185–170 Ma (samples 46B, 49B, and 83), consistent with previous interpretations (e.g., [41]). Deformation along the Wann River shear zone synchronous with decompression and/or cooling suggests that it may be one of the structures that accommodated tectonic exhumation, although the initial geometry and shear sense of the shear zone are not well constrained (e.g., [41, 44, 109]).

Monazite inclusions in kyanite (samples 31 and 42) and some inclusions in garnet (samples 42, 83, and 77A) yield dates between 185 and 170 Ma that are difficult to reconcile with other petrochronological interpretations if they represent the maximum age of porphyroblast growth (Figures 7 and 8). However, all monazite inclusions in kyanite intersect cracks filled with white mica (Figures 3(e) and 9(d) and Figure S4), suggesting that monazite may have formed during the partial replacement of kyanite by white mica, not before or during kyanite crystallization (e.g., [97, 110]). In sample 42, monazite inclusions in kyanite have higher Y and lower Gd/Yb compared to synchronous monazite in the matrix, compatible with (re)crystallization during garnet resorption. An inverse trend is observed in sample 31. However, the other two monazite grains from the same compositional layer located away from pseudomorphed garnet also have lower Y and higher Gd/Yb compared to monazite in the matrix, which suggests that the REE budget in sample 31 was controlled by local availability of these elements in the layer rather than being linked to a reaction involving kyanite and white mica (Figure S5). Thus, the age of kyanite porphyroblasts in the matrix of these samples cannot be determined confidently; they crystallized prior to 185 Ma but could be either Permian-Triassic or earliest Jurassic. The 185–175 Ma inclusions in garnet in samples 42, 83, and 77A are part of polymineralic inclusions near embayments (Figure S4) and thus could be in contact with the outside of garnet but appear enclosed due to the 2D sectioning of a complex 3D geometry. Monazite (re)crystallization in garnet hosts is also possible via cracks that could serve as pathways for fluids (e.g., [111]). The higher Y and lower Gd/Yb from these monazite inclusions in garnet compared to synchronous monazite in the matrix suggest a local increase in Y and HREE budget, likely related to the resorption of garnet. (Re)crystallization of 185–170 Ma monazite inclusions in kyanite (samples 31 and 42) and garnet (samples 42, 83, and 77A) after porphyroblast crystallization and during garnet breakdown in the sillimanite stability field (Figure 10(a)) is compatible with abundant 185–175 Ma xenotime near resorbed garnet (Figure 3(c)) and 185–175 Ma monazite and xenotime in textural contact or intergrown with sillimanite (Figures 5(d) and 9(g).

Two samples from the Willison Glacier area (25B and 31) contain abundant randomly oriented porphyroblasts of cordierite and andalusite (Figure 3(g)), providing new constraints on low-P metamorphism in this area (Figure 10). Furthermore, garnet in sample 25B is partially replaced by coronae of $Pl+Bt+Fe‐sulfide+Crdinnercorona+Andoutercorona$⁠, suggesting that andalusite crystallized prior to cordierite (Figure 3(b)). The youngest monazite spots are $161±6Maand164±5Ma$ in samples 25B and 31, respectively (Figures 7, 9(a), and 9(c)). These monazite grains are included in andalusite (both samples) and in cordierite (sample 25B), indicating that andalusite and cordierite grew after ca. 160 Ma. In sample 25B, xenotime (re)crystallized between ca. 175 and ca. 120 Ma and thus provides a more complete dataset (Figure 7). The youngest inclusions in andalusite (⁠$145±6$ Ma) and cordierite (⁠$121±4$ Ma) provide an upper limit on the crystallization age of these minerals (Figure 9(b)). The ca. 25 Myr older xenotime dates from inclusions in andalusite are consistent with crystallization of andalusite before cordierite, as indicated by petrographic observations. However, this 25 Myr interval may not be representative of the true time interval between andalusite and cordierite crystallization because it is based on a limited number of inclusions of variable ages, which constrain the maximum age of porphyroblasts rather than the timing of their growth. The spread of xenotime dates between 170 and 120 Ma indicates that P–T conditions remained sufficiently high (i.e., chlorite zone or higher; [67]) for xenotime (re)crystallization throughout this period. Sample 49B from the Eastern Wann River shear zone also contains cordierite, but the lack of monazite or xenotime inclusions hinders the assessment of its timing of crystallization. Because cordierite overprints the foliation in this sample, it is probable that it crystallized during a posttectonic overprint, similar to the one observed in Willison Glacier area.

The Florence Range suite records a protracted metamorphic evolution spanning over 150 Myr from the Permian to the Early Cretaceous. Linking monazite and xenotime data with the crystallization sequence of Al₂SiO₅ polymorphs and additional metamorphic index minerals (K-feldspar and cordierite) provides first-order constraints on the metamorphic P–T conditions of each metamorphic event (Figures 6(b) and 10 and Figure S3). Kyanite crystallized first during the earliest preserved metamorphic event at 270–240 Ma. P–T conditions during Permian to Middle Triassic metamorphism are constrained to >600°C and >0.6 GPa based on kyanite inclusions in garnet interpreted to have crystallized at that time (Figure 3(d)). The presence of staurolite inclusions in garnet also supports a prograde path that eventually reached kyanite grade (Figures 3(d) and 5(c)). However, the high sensitivity of staurolite stability to the geochemical composition of the model system (Figure 6(b) and Figure S3) precludes additional inferences on the prograde P–T conditions. Retrograde metamorphism to unconstrained P–T conditions likely followed during the Late Triassic (240–215 Ma). The age of kyanite porphyroblasts in the matrix cannot be determined with available data, but they may have crystallized during the Permian to Middle Triassic metamorphic event and/or during a second episode of prograde metamorphism documented during the earliest Jurassic (195–185 Ma). Abundant fibrolite and prismatic sillimanite and rare K-feldspar crystallized between 185 and 170 Ma at peak T under suprasolidus conditions (>675°C and< 0.8 GPa). Significant garnet breakdown during this metamorphic phase is best explained by decompression and/or cooling along this segment of the P–T–t path. Finally, andalusite and cordierite crystallized below 0.4 GPa after ca. 120 Ma. The garnet replacement texture that resulted in sequential andalusite and cordierite crystallization observed in one sample (Figure 3(b)) is most consistent with a third heating event that reached at least 600°C based on the lowest temperature at which cordierite is predicted to be stable in the chemical compositions of our model systems.

The polymetamorphic Al₂SiO₅ crystallization sequence that we identified (Ky→Sil→And) based on petrographic observations and monazite and xenotime petrochronology is different than those from previous studies that were based solely on petrographic observations or indirect evidence (And→Ky→Sil, [7]; Ky→And→Sil, [8]). Currie [7] does not describe replacement textures of andalusite by kyanite, or the relationship between andalusite and sillimanite, and notes that the relationship between kyanite and tectonic fabrics is uncertain because kyanite is isolated by a corona of $whitemica±sillimanite$ (fibrolite). Partial replacement of kyanite by sillimanite is compatible with our petrographic observations and petrochronological interpretations. Mihalynuk et al. [8] describes andalusite directly overgrowing kyanite, as we observed in our sample 31 (Figure 3(f)), and swirling sillimanite overgrowing andalusite. In contrast, we found clear evidence that sillimanite crystallized before andalusite (i.e., sillimanite intergrown with Early Jurassic monazite and andalusite overgrowing Early Cretaceous xenotime inclusions). Nonetheless, we cannot rule out the possibility for a second generation of sillimanite postdating andalusite growth (e.g., [8]), during the Early Cretaceous contact metamorphism. This possibility is not shown in Figure 10 because we did not observe unambiguous evidence for this second generation of sillimanite. Our study thus demonstrates the utility of acquiring absolute age constraints linked texturally to Al₂SiO₅ polymorphs crystallization to decipher polyphase metamorphic evolution.

At the scale of this study, petrochronology data are consistent between samples and do not indicate that individual metamorphic phases were diachronous across the Atlin area. However, some metamorphic phases are only recorded in certain samples, likely due to either poor preservation of petrochronometers, or the limited influence of a metamorphic phase on the (re)crystallization of petrochronometers. With the exception of the Early Mississippian metamorphism (eclogite facies metamorphism, [20, 23]; unconstrained P–T conditions, [22, 24]), the Florence Range suite records similar periods of metamorphism as the Yukon-Tanana terrane in Yukon and northern British Columbia (Table 1). However, the P–T conditions recorded in the Florence Range suite do not necessarily match the P–T conditions recorded elsewhere in the Yukon-Tanana terrane (see the following).

The Laberge Group east of the Llewellyn fault zone (Figure 1) records deep erosion of the Late Triassic Stuhini-Lewes River arc indicated by a change in provenance from volcanic- to plutonic-dominated clasts [112–115]. The Laberge Group is locally characterized by a brief influx of metamorphic detritus at ca. 180 Ma (e.g., [30, 40, 114, 116–118]), including 202–182 Ma high and ultra-high-P clasts (#11 in Figure 1; [30, 117, 118]). The source of these clasts remains enigmatic, though Florence Range and the adjacent Boundary Ranges suites have been previously suggested as the source of petrographically similar clasts and consistent with north-easterly paleoflow directions commonly preserved in the late Pliensbachian Laberge Group [113]. We did not identify any evidence of earliest Jurassic high-P metamorphic rocks that could be the source of the high and ultra-high-P clasts in the Laberge Group [30]. Furthermore, our data indicate that the Florence Range suite was in the sillimanite stability field during the deposition of metamorphic detritus in the Laberge Group (ca. 180 Ma) and was not exhumed to the surface until at least 120 Ma, the age of andalusite- and cordierite-grade metamorphism. Consequently, the presently exposed Florence Range suite is unlikely to have been the direct source of the metamorphic clasts in the Laberge Group. It is possible that structurally higher, now eroded metamorphic units contributed lower grade metamorphic detritus into the Laberge Group basin as garnet resorption (and potentially decompression and exhumation) in the Florence Range suite started at ca. 185 Ma. Alternatively, earliest Jurassic metamorphic detritus may have come from elsewhere in the Yukon-Tanana terrane, which experienced Barrovian metamorphism during the Early Jurassic (Table 1; see Section 7.3.2). The modern proximity of metamorphic detritus-containing Laberge Group and Florence Range suite is in part controlled by strike-slip displacement along structures such as the Llewellyn fault zone [42, 43]; however, along-basin sediment transport in the Laberge Group throughout the Jurassic [113] may have sourced metamorphic detritus from the Yukon portion of the Yukon-Tanana terrane.

The Yukon-Tanana terrane preserves a complex tectonometamorphic history that has been attributed to its evolution from a late Paleozoic arc to late Paleozoic or early Mesozoic collision with Laurentia (Table 1). The middle Permian to Late Triassic Barrovian metamorphism is the earliest tectonometamorphic event preserved in the Florence Range suite. This episode overlaps the middle to late Permian eclogite facies (e.g., [19, 21, 119]), late Permian to Middle Triassic amphibolite facies [24], and Middle Triassic low pressure, pyroxene hornfels facies [96] metamorphism in Yukon. Eclogite facies metamorphism is only recognized in the eastern part of the Yukon-Tanana terrane (#5–7 in Figure 1). It has been attributed to subduction erosion of the Yukon-Tanana terrane in an overriding plate setting [120], the result of a collision with the downgoing Laurentia [25], or the result of a collision with an intra-oceanic arc in the overriding plate (Slide Mountain terrane, Figure 1; [32, 33]). The Middle Triassic amphibolite facies event (ca. 260-240 Ma) in the Stewart River area (#1 in Figure 1) is interpreted to have been the result of intra-arc thickening induced by migration of the Yukon-Tanana terrane towards Laurentia [24]. The low-pressure metamorphic event (pyroxene hornfels facies) documented in the Stewart River area ~100 km to the south at $245±1$ Ma [96] may be related to middle to late Permian magmatism (Klondike cycle of Nelson et al. [18]) or extensional tectonics in a complex tectonic configuration [32, 33, 121–123].

Middle Permian to Late Triassic Barrovian metamorphism in the Florence Range suite is distinctly different from several other localities of the YTT. The presence of kyanite and staurolite inclusions in chemically distinct garnet suggests similar or deeper burial during the middle Permian to Middle Triassic than in the Stewart River area [24, 96], but not at as high pressure as documented in the eastern part of the Yukon-Tanana terrane (e.g., [19, 21, 119]). Permian igneous rocks in the Florence Range suite (i.e., $270±5$ Ma Wann River gneiss; [45]) are restricted to the Wann River shear zone, would not have contributed to significant burial, and do not occur in proximity to all samples that record evidence for a Permian-Triassic metamorphic event (e.g., sample 31, location in Figure 1). Furthermore, Permian-Triassic monazite and xenotime crystallization lasted approximately 55 Myr after the crystallization of the Wann River gneiss, an interval far too long to be associated with contact metamorphism around a unit that is spatially restricted. Permian igneous rocks are therefore unlikely to be the direct cause of this metamorphic episode. Middle Permian to Middle Triassic burial of the Florence Range suite to mid-crustal depth (i.e., >0.6 GPa) is consistent with contraction-related burial; however, the cause of this contractional event cannot be discerned in this area. The Permian-Triassic tectonometamorphic evolution of the Yukon-Tanana terrane thus appears to vary significantly across the terrane, which suggests that this terrane should not be viewed, at least from a tectonometamorphic perspective, as a homogeneous entity during that time.

The pervasive Early Jurassic middle to upper amphibolite facies metamorphic event is the most recognizable event throughout the Yukon-Tanana terrane (Table 1; e.g., [7, 24, 26, 96, 98]) and is also recorded in the Florence Range suite. P–T conditions in the Florence Range suite (⁠$suprasolidussillimanite+K$-feldspar grade) are broadly similar to the rest of the Yukon-Tanana terrane, though coeval Jurassic plutonism (e.g., Tagish Lake Suite, [8]; Minto and Long Lake suites, [51]) likely affected peak temperature conditions in close proximity to these plutons. The overall homogeneity in P–T conditions suggests that a large portion of the Yukon-Tanana terrane underwent a similar tectonic evolution during the Early Jurassic, in contrast to the heterogeneous metamorphism recorded during the middle Permian to Late Triassic. The Early Jurassic tectonometamorphism has been attributed to the collision of the Yukon-Tanana terrane with Laurentia (e.g., [24, 31, 98]) or with the Stikine terrane (e.g., [30, 35, 41, 44]).

The Yukon-Tanana terrane in the study area is separated from the Laurentian margin by other Intermontane terranes (Figure 1). Hence, it is not feasible to assess models of Yukon-Tanana–Laurentia collision here. However, the proximity of the study area to the Stikine terrane does allow comparison of the tectonic history even though the Stikine terrane was metamorphosed at much lower grade (zeolite to greenschist facies; [8, 39, 124]). Similar to Florence Range, there is ample evidence of Late Triassic–earliest Jurassic deformation within the northern Stikine terrane. This evidence includes Late Triassic burial and metamorphism of the Late Triassic Lewes River Group (equivalent to Stuhini Group in British Columbia) volcanic rocks near the contact with the Yukon-Tanana terrane in central Yukon (Carmacks copper deposit; [34]). Late Triassic Lewes River–Stuhini arc underwent magmatic shut off (e.g., [125]), followed by development of an extensive latest Triassic fold-and-thrust belt (e.g., [125–127]), uplift and erosional incision into the plutonic middle crust in the Early Jurassic [113, 114, 128], and Early to Middle Jurassic Hazelton Group magmatism and sedimentation [35, 125, 126, 128–130]. The coeval contractional deformation of the Stikine and Yukon-Tanana terranes may indicate that these terranes either collided in the Late Triassic or collided prior to that and underwent another contractional episode together. The former interpretation is more consistent with the juvenile isotopic compositions and lack of zircon inheritance in ca. 215 Ma and older plutonic suites that intrudes the Stikine terrane, suggesting no interaction with evolved continental crust until the latest Triassic [131]. The Stikine and Yukon-Tanana terranes are accepted to be together by Late Triassic to Early Jurassic because they share magmatic (Minto, Long Lake and Tagish suites; [8, 51]) and sedimentary histories [40, 114, 115]. The initiation of the contractional event in the Yukon-Tanana terrane and Stikine terrane predates the Early Jurassic crystallization of metamorphic monazite (ca. 195 Ma) in the Florence Range suite, reflecting the time lapse necessary to burry and heat the rocks to metamorphic conditions sufficient for monazite crystallization.

Overall, the Permian-Triassic and Early Jurassic tectonometamorphism in the Yukon-Tanana terrane suggest two contraction-related episodes of deformation and metamorphism. However, there is significant uncertainty with respect to what Yukon-Tanana terrane interacted with. For example, Jurassic collision of Yukon-Tanana terrane with the Laurentian margin is inconsistent with the presence of detritus derived from the Yukon-Tanana terrane in the Triassic sediments on the Laurentian margin [25]. Similarly, Late Triassic to Jurassic collision of Yukon-Tanana and Stikine terranes is inconsistent with the widely accepted basement-cover relationship of the Yukon-Tanana terrane below the Stikine terrane (e.g., [12–14, 38, 132, 133]). These uncertainties cannot be solved without reevaluation of the tectonometamorphic history of all Intermontane terranes.

The Early Cretaceous in the Florence Range suite is characterized by posttectonic low-P metamorphism that mainly affected samples in the Willison Glacier area. There, several meter-scale undeformed granitic dykes and sills crosscut the metamorphic units and are probably part of the Cretaceous to Eocene Coast Plutonic Complex ([8, 52, 134] and references therein). These and related intrusions at depth likely drove low-P contact metamorphism in parts of the Florence Range suite. Similar contact metamorphism is also documented in the Stewart River area to the north (#1 in Figure 1; [24]). Further east near the tectonic boundary with Laurentia, contrasting metamorphic conditions indicating Barrovian metamorphism during the Middle Jurassic to Early Cretaceous in the Finlayson Lake district (#8 in Figure 1) have been interpreted as the result of foreland-ward propagation of regional deformation followed by mid-crustal exhumation [29]. Barrovian metamorphism also occurred during the Early Cretaceous in the Australia Mountain area (#2 in Figure 1), which is interpreted as a metamorphic core complex exposing the Cretaceous infrastructure of the Yukon-Tanana terrane [28]). However, this metamorphic domain lacks post-Early Mississippian arc magmatism and pre-Middle Jurassic metamorphism, characteristic of the Yukon-Tanana terrane, and is more likely part of parautochthonous North America (Cassiar terrane in Figure 1; [31]). In southwest Yukon, along the tectonic boundary between the Yukon-Tanana terrane and the Wrangellia terrane to the west (#4 in Figure 1), the Snowcap assemblage records a Late Cretaceous amphibolite facies overprint attributed to the final collision of the Insular terranes with the Yukon-Tanana terrane [97]. Altogether, these data demonstrate that the tectonometamorphic evolution of the Yukon-Tanana terrane following the main phases of collisional orogenesis was highly variable depending on E-W positions, the across-strike position in the greater Cordilleran orogen (present coordinates). Localities in the middle of the Yukon-Tanana terrane, such as Stewart River, Aishihik Lake, and Atlin areas (# 1, 3 and 10 in Figure 1), may have a better potential to preserve the Early Jurassic and older tectonometamorphic evolution of the terrane, as they appear to have escaped tectonic overprint during the final assembly of the Canadian Cordillera. Post-Jurassic contact metamorphism may overprint previous events, so it is advisable to avoid exposures close to known Cretaceous to Eocene plutons.

The in situ analysis of petrochronometers and the evaluation of their textural relationships with metamorphic index minerals, particularly Al₂SiO₅ polymorphs, is key to reconstruct the first order metamorphic evolution of a metamorphic domain. Laser-ablation split-stream petrochronology on monazite and xenotime records complementary segments of the P–T–t path, illustrating the usefulness of xenotime even though this petrochronometer is not as widely used as monazite. Our new data constrain the following sequence of metamorphic phases that have affected all or parts of the Florence Range suite:

• (1)

  Permian to Middle Triassic (270–240 Ma): first episode of cryptic prograde metamorphism to the kyanite stability field (>600°C, >0.6 GPa)

• (2)

  Late Triassic (240–215 Ma): retrograde metamorphism to unconstrained P–T conditions

• (3)

  Earliest Jurassic (195–185 Ma): second episode of prograde metamorphism to unconstrained P–T conditions

• (4)

  Late Early Jurassic (185–170 Ma): pervasively recorded metamorphism in the suprasolidus stability field of sillimanite and K-feldspar, during decompression and/or cooling (>675°C, <0.8 GPa)

• (5)

  Middle Jurassic to Early Cretaceous (170–120 Ma): retrograde metamorphism without complete exhumation

• (6)

  Late Early Cretaceous (<120 Ma): locally recorded third episode of prograde contact metamorphism through the stability fields of andalusite and cordierite (>600°C, <0.3–0.4 GPa)

This sequence of metamorphic phases indicates that two separate contractional events affected the Yukon-Tanana terrane and were overprinted by local contact metamorphism. Equally detailed investigations of the P–T–t evolution of metamorphic units in adjacent terranes of the Northern Cordillera will be required to resolve the uncertainties relative to what Yukon-Tanana terrane interacted with during each of these contractional events.

Datasets related to this article are included in the supplementary material (Tables S1 and S2).

An earlier version of this manuscript has benefited from comments by D. Kellett, M. Smit, D. Moynihan, D. Gibson, and D. Pattison. This manuscript is NRCan contribution 20210423.

The authors declare no conflict of interest with respect to this publication.

We thank Pat Hunt and Daniele Regis for assistance on the scanning electron microscope, Katherine Venance for electron microprobe mapping, and Lisel Currie for supplying unpublished first-hand documentation about the geology of Atlin area. Norm Graham of Discovery Helicopters provided reliable transportation to the field area. This research was funded by the Geological Survey of Canada, GEM-II Cordillera project.