Metamorphic evolution of the Great Slave Lake shear zone

The Palaeoproterozoic Great Slave Lake shear zone (GSLsz) is a crustal‐scale strike‐slip structure, with a total length of >1,000 km and a width of ~25 km, that separates the Archean Rae and Slave cratons. The range of metamorphic rocks now exposed at the surface encompasses granulite facies mylonite through to lower‐greenschist facies ultramylonite and cataclasite, providing a potential type example of fault‐zone structure in the middle and lower crust. However, the metamorphic evolution of the units remains poorly quantified, hindering detailed structural and tectonic interpretations. Here, we use phase equilibria modelling and thermobarometry to determine the metamorphic conditions recorded by pelitic, mafic and felsic GSLsz mylonites. Samples from the entire range of granulite–greenschist facies units preserve evidence for nested clockwise pressure–temperature paths that are consistent with a single orogenic cycle. Our findings indicate that the northern Rae margin underwent pervasive crustal thickening with peak pressures in metasedimentary rocks reaching ~1.1 GPa. The crustal thermal gradient at the onset of thickening was ~650°C/GPa, whereas the final stages of equilibrium recorded by fine‐grained matrix minerals in all samples collectively define a metamorphic field gradient of ~1,000°C/GPa. Deformation microstructures are consistent with the main phase of dextral shear having been synchronous with or following peak metamorphism. The history of metamorphism and exhumation of the GSLsz is consistent with the Sibson–Scholz model for shear zones, with a narrowing of the deforming zone and the progressive overprinting of higher‐grade assemblages during exhumation through shallower crustal levels.


| INTRODUCTION
Characterization of the structural and petrological evolution of plate-boundary fault zones (sensu lato, including both upper crustal faults and their underlying shear zones) is central to modelling the occurrence (or absence) of plate tectonics. Observations of both active and ancient fault zones are typically interpreted in the context of the Sibson-Scholz model of fault-zone structure (Scholz, 1988;Sibson, 1977Sibson, , 1983. In this model, increases in pressure and temperature with depth induce a transition from dominantly pressuresensitive frictional failure in the upper crust to dominantly temperature-sensitive viscous flow in the lower crust, with a concomitant transition from narrow brittle faults to broad ductile shear zones. However, close examination of individual structures has revealed several complicating effects, including exhumation and/or burial during fault activity (e.g. Parsons et al., 2016;Toy et al., 2008), progressive fault-zone weakening (e.g. Holdsworth, 2004;Imber et al., 2008;Wallis et al., 2013) and reactivation (potentially with different kinematics) during distinct tectonic events (e.g. Holdsworth et al., 1997Holdsworth et al., , 2001. These effects are most prevalent in the continental lithosphere due to its longevity and structural complexity. As such, it is important to identify and characterize those continental fault zones that may pertain most closely to the basic Sibson-Scholz model to provide comparators against which to develop more complex models. The Great Slave Lake shear zone (GSLsz) is a potential type example of a deeply eroded continental transform plate boundary. Early studies of the GSLsz helped to establish the current paradigm for the structure of major fault zones in the continental crust (Hanmer, 1988(Hanmer, , 1997Hanmer & Lucas, 1985;Hoffman, 1987). Importantly, the central portion of the GSLsz runs through homogeneous felsic intrusive lithologies that allow the effects of deformation conditions on fault-zone structure to be isolated while mafic and pelitic enclaves are convenient for providing geothermobarometric constraints (Hanmer, 1988;Hanmer et al., 1992). Previous workers identified mylonitic units, with distinct metamorphic mineral assemblages, that exhibit broadly similar/compatible kinematic indicators and occur in an overprinting sequence of decreasing width, decreasing metamorphic grade and increasing intensity of deformation fabrics (e.g. grain size reduction, foliation development ;Hanmer, 1988;Hanmer et al., 1992). The mylonitic units with the lowest metamorphic grades are partially overprinted by brittle fault rocks (Hanmer, 1988). These characteristics are consistent with the Sibson-Scholz model of fault-zone structure that was initially developed for fault zones dominated by a reverse component, and thereby began to establish the relevance of the model to major strikeslip structures (Hanmer, 1988).
Despite their significance for models of continental faultzone structure, estimates of the metamorphic conditions recorded by each unit within the GSLsz have been based on sparse quantitative constraints (Hanmer, 1988). The lack of modern metamorphic studies on the GSLsz has left open questions regarding whether or not the laterally exposed metamorphic units reflect the vertical structure of a long-lived strike-slip fault zone. An alternative hypothesis is that the structure could have been reactivated at different crustal levels during temporally distinct deformation events. This latter model has been associated with several other major fault zones (e.g. Holdsworth et al., 1997Holdsworth et al., , 2001. In this contribution, we use phase equilibria modelling and thermobarometry to constrain pressure-temperature (P-T) paths for the four main metamorphic units of the GSLsz. The improved P-T constraints provided here are a significant development in discriminating between the competing hypotheses. The conditions of final equilibration, which fall on a single metamorphic field gradient despite spanning greenschist facies to granulite facies conditions, are consistent with all the mylonitic units having formed during a single orogenic cycle in a single overarching structure.

| GEOLOGICAL CONTEXT
Northern Laurentia is an amalgam of the Archean Slave, Rae, Hearne and Superior cratons, all of which are separated by Palaeoproterozoic belts that have unique magmatic, metamorphic and deformation histories (Hoffman, 1987). The tectonic boundary between the Slave and Rae cratons forms one of these Palaeoproterozoic belts and is itself divided, from north to south, into the Thelon Tectonic Zone (TTZ), GSLsz and the Taltson Magmatic Zone (TMZ; Figure 1a).
The TTZ and TMZ were initially thought to be continuous and have formed as a result of the collision and subsequent indention of the Slave craton into the Rae craton, with the GSLsz acting as a dextral strike-slip structure (Hoffman, 1987). However, recent geochronologic and isotopic data have cast doubt on the contiguous nature of the TTZ and TMZ. Instead, the north and south sections of the Slave-Rae boundary may have formed at a similar time but in distinctly different tectonic environments Card et al., 2014;Chacko et al., 2000;Whalen et al., 2018). The TTZ is interpreted as a continental arc that is related to the closure of an oceanic basin between the Slave craton and the northwestern part of the Rae craton at c. 1.97 Ga during the Thelon orogeny (Hoffman, 1987;Thériault, 1992;Whalen et al., 2018). In contrast, the TMZ is interpreted to have formed in a plate-interior rather than plate-margin setting (Chacko et al., 2000;Schultz et al., 2007). In this model, which is supported by the absence of a mantle signature in the chemistry of TMZ intrusions, the high-grade metamorphism and magmatism that characterizes the TMZ is the result of crustal thickening that itself was caused by far-field stresses associated with an active plate boundary located either to the west of the adjacent Buffalo Head terrane or to the east of the Rae craton (Chacko et al., 2000).
The TTZ records slightly older magmatic ages (2.07-1.96 Ga; Berman et al., 2018) than the TMZ (1.99-1.92 Ga; Bostock & Loveridge, 1988;Bostock, Van Breemen, & Loveridge, 1987, 1991Chacko et al., 2000). However, this trend is reversed for the metamorphic ages, with older metamorphic ages in the TMZ (1.94-1.92 Ga; Bethune et al., 2013;McDonough et al., 2000) and younger metamorphic ages in the TTZ (1.92-1.89 Ga; Berman et al., 2018). Further complicating matters, in the Union Island Group (part of the Great Slave Supergroup), which is located along the southeast margin of the Slave craton and adjacent to the GSLsz (Figure 2a), there is evidence for rift volcanism at c. 2.05 Ga (Sheen et al., 2019). Based on these findings, Sheen et al. (2019) suggest that the Slave and Rae cratons may have formed a contiguous crustal block prior to c. 2.05 Ga rifting.
The interpretation of Sheen et al. (2019) is consistent with models for the formation of Laurentia that include an earlier collision of Slave and Rae cratons during the c. 2.5-2.3 Ga Arrowsmith orogeny, that is now recognized along much of the western margin of the Rae craton Hartlaub et al., 2007;McNicoll et al., 2000;Schultz et al., 2007;Tersmette, 2012). Evidence in support of the Arrowsmith orogeny includes bimodal magmatism at 2.50-2.45 Ga (Schultz et al., 2007;Tersmette, 2012), silicic magmatism at 2.49-2.47 Ga  and the intrusion of post-orogenic granites at c. 2.3 Ga (Hartlaub et al., 2007;Tersmette, 2012). Metamorphism associated with the Arrowsmith orogeny ranges from upper-amphibolite to granulite facies between c. 2.5 and 2.3 Ga, and together with the magmatism, likely represents either an Andean-type margin with associated continental arc magmatism  or a Himalayan-style collision (Schultz et al., 2007;Tersmette, 2012).
With the available data, it is unclear how both the timing and tectonic setting of the GSLsz relate to the formation of the TTZ and TMZ. As the central segment of the Slave-Rae boundary, the GSLsz itself can be delineated over 1,300 km using aeromagnetic data, from the foothills of the Rocky Mountains in the west to the Thelon Basin in the east (Figure 1b). The mylonites of the GSLsz form a continuous structure that joins the north-trending linear magnetic highs of the TTZ and the TMZ (Figure 1b). Kinematic indicators within the GSLsz record dextral strike-slip motion and, based on reconstructions using the aeromagnetic data, the GSLsz accommodated over 700 km of dextral displacement (Hoffman, 1987). Geochronological constraints from the GSLsz are sparse, and the existing data were not collected with modern methods so we are unable to assess whether they represent mixing of multiple age domains. The best estimates of the time of deformation, based on zircon U-Pb data from variably deformed granitic dykes, indicate that the GSLsz first became active sometime between c. 2.00 and 1.86 Ga, possibly reactivating an earlier c. 2.56 Ga structure (Bowring et al., 1984;Hanmer et al., 1992). The cessation of ductile shear in the GSLsz is indirectly constrained by an unconformably overlying quartzite of the Great Slave Supergroup, which yields a maximum deposition age of 1,857 ± 11 Ma (Bowring et al., 1984).
Following the main ductile shearing event, the GSLsz and the adjacent southern margin of the Slave craton were dextrally offset by the McDonald fault (Bowring et al., 1984;Figure 2). Both the McDonald fault and the brittle component of the conjugate sinistral Bathurst fault (Gibb, 1978; Figure 1) are reported to have been active sometime after c. 1.84 Ga (Henderson & Loveridge, 1990;Ma et al., 2020) and prior to deposition of the Thelon basin at c. 1.72 Ga (Miller et al., 1989; Figure 1). F I G U R E 1 Regional maps. (a) Bedrock map for northern Laurentia presenting the locations of the Great Slave Lake shear zone (GSLsz) and conjugate Bathurst Fault (BF). Palaeoproterozoic belts: TMZ-Taltson Magmatic zone, TTZ-Thelon Tectonic Zone, WO-Wopmay Orogen, STZ-Snowbird Tectonic Zone, THO-Trans Hudson Orogen. Palaeoproterozoic basins: TB-Thelon Basin, ATH -Athabasca Basin. Modified after Ashton et al. (2013) and Thiessen et al. (2019). (b) Aeromagnetic map modified after Miles and Oneschuck (2016) [Colour figure can be viewed at wileyonlinelibrary.com] If the hypothesis that the TTZ and the TMZ do not define a continuous orogenic belt is correct, and/or the Slave and Rae cratons first collided during the Arrowsmith orogeny, then how does the GSLsz relate to the TTZ and TMZ in terms of both timing and tectonic affinity? In this paper, we provide new P-T estimates for representative samples from the GSLsz, which is the first step in elucidating the tectonometamorphic history of the central Slave-Rae boundary.

SAMPLE SELECTION
The width of individual map units within the GSLsz narrows with decreasing metamorphic grade (Hanmer et al., 1992), from ~20 km for the granulite, ~10 km for the upper-amphibolite mylonite, ~1 km for the epidote-amphibolite mylonite and ~100 m for the greenschist ultramylonite. Two major brittle faults, with unknown offsets, run parallel to the strike of the F I G U R E 2 (a) Inset indicating the location of our study area relative to the Slave-Rae boundary. (b) Map of metamorphic grades for the southwestern exposure of the Great Slave Lake shear zone. Shading corresponds to the protolith lithology as mapped by Hanmer (1988). (c) Closeup of the southwestern study area with the location of the samples collected for this study marked by white circles. The samples that were selected for petrological modelling are labelled on frames (b)  mylonitic foliation (Figure 2). To the north, the McDonald fault delimits the northernmost boundary of the GSLsz, juxtaposing greenschist ultramylonite against foliated plutonic rocks from the Slave craton and volcanic and sedimentary units from the Union Island Group. The second fault, which we term the Rocher River fault, runs through the centre of the GSLsz and separates epidote-amphibolite mylonite from greenschist ultramylonite. In addition to the metamorphic juxtaposition, the Rocher River fault divides the southwestern portion of the exposed GSLsz into two structural domains. South of the Rocher River fault, the mylonitic foliation strikes NE-SW, parallel to the strike of the fault. Whereas, to the north of the Rocher River fault, the mylonitic fabric strikes NNE-SSW and appears to be cut by the fault. Recent glaciation and the present-day arid climate contribute to near-continuous bedrock exposure in the southwestern portion of the exposed GSLsz. With the exception of the juxtaposition of units along late, minor brittle faults, the boundaries between the map units are diffuse and can be seen to be the result of lower metamorphic-grade assemblages superimposed on higher-grade rocks. For example, the boundary between the upper-amphibolite and epidote-amphibolite map unit occurs over ~100 m, with a progressive increase in the prevalence of the epidote-amphibolite assemblage overprint with distance into the area mapped as epidote-amphibolite facies. The overprinting of former, higher-grade metamorphic assemblages by lower-grade assemblages is associated with the strain related to pervasive right-lateral shear. Field evidence for the progressive nature of this overprint and its association with ductile deformation include the observation of upper-amphibolite mineral assemblages (garnet-present) preserved in the low-strain cores of mafic boudins that are flanked by Garnet-absent and epidote-bearing high-strain margins ( Figure S1). Given the progressive nature of the map-unit boundaries, we adopt the convention used by previous workers and mark the boundaries for each unit based on the first appearance of the lower-grade mineral assemblage.
In total, 60 samples and 50 cores were collected along two profiles that transect the southwestern exposure of the GSLsz (Figure 2c) with an additional eight samples and five cores collected along strike (~90 km to the northeast) from an area previously mapped as granulite facies metaplutonic and metasedimentary units (Figure 2b). Of all samples collected, six representative samples were selected to provide quantitative P-T estimates for the granulite, upper-amphibolite mylonite, epidote-amphibolite mylonite and greenschist ultramylonite map units.

| Analytical Methods
For the purpose of petrological modelling, effective bulk compositions were determined by combining mineral composition data with thin-section X-ray phase maps. Thinsection X-ray phase maps and mineral-area analyses were obtained at SGS, Burnaby, Canada using the QEMSCAN energy dispersive spectrum-based (EDS) automated mineralogy system mounted on a tungsten-source Tescan Vega scanning electron microscope. Entire thin sections were analysed at an excitation voltage of 20 kV and beam current of 10 nA, with a step size of 7 µm, providing high-resolution modal (area) analyses for all phases ( Table 1). The fully automated QEMSCAN system integrates backscattered electron (BSE) and EDS signals to characterize polished sample surfaces (Gottlieb et al., 2000). Maps of BSE brightness were calibrated against quartz, gold and copper standards. Spectra were processed using user-defined species identification protocol files that characterize phases on the basis of characteristic EDS spectra and BSE intensity. The high-resolution phase maps for all samples described in this study are available in Figure S2.
Mineral compositions were obtained at the University of British Columbia using a CAMECA SX-50 electron microprobe analyser (EMPA). Operating conditions include an excitation voltage of 15 kV, a beam current of 20 nA, a peak count time of 20 s with a background count time of 10 s and a spot size of 5 µm. The standards, X-ray lines and crystals used are provided in Appendix S1. Data reduction was performed using the method of Pouchou and Pichoir (1985). For all phases, locations for analysis were chosen to represent the full range of petrographic positions (e.g. matrix, corona, porphyroblast core and rim), and multiple spots were collected from select grains to assess inter-and intra-crystalline compositional variation. For all applicable phases, Fe 3+ was calculated on the basis of charge-balance using the software AX 6.2 (Holland, 2018).
While both thermodynamic data sets were tested for all samples, ds-55 was chosen for the pelitic and metagranodioritic rocks over the more recent ds-62 because it better predicted the observed mineral assemblages, and in particular, the composition and stability of Garnet in these samples (cf., Starr et al., 2020;Waters, 2019). Data set ds-62, with the associated activity-composition models, replicated the observed mineral assemblages for the mafic samples reasonably well with the one exception of augite, which was predicted to be stable over a wide range of conditions but is not present in either sample. The overstabilization of 1-2 modal% augite in petrological models at amphibolite facies conditions has been documented by multiple recent studies (Forshaw et al., 2019;Weller et al., 2019). As such, we did not rely upon the occurrence of augite in the calculated pseudosections for mafic samples for our interpretation of P-T conditions. Uncertainty on the absolute positions of phase assemblage boundaries in pressure-temperature space generally do not exceed ±20-30°C and ±0.05 GPa at 1 SD Powell & Holland, 2008), although the relative errors may be smaller when comparing results generated using the same data set and associated activity-composition models. Since the results we present here were calculated using two data sets, the absolute errors remain an important consideration when we compare the P-T estimates among samples.
The amount of H 2 O was set to minimally saturate (i.e. 0.05 mol.% free H 2 O) the assemblage in the immediate sub-solidus at a pressure of 0.7 GPa for all samples with peak mineral assemblages formed above the solidus (see Figure S3 for an example T-X H2O diagram). Bulk rock X Fe 3+ = Fe 3+ /(Fe total ) was set for each sample and is a function of the Fe 3+ determined for all ferric minerals by stoichiometric charge-balance ( Table 2). The effective bulk compositions for all phase diagrams discussed below are provided in Table 3.

| Thermobarometry
To complement the petrological modelling and provide independent constraints on the metamorphic conditions recorded in our samples, we applied two sets of independent thermometers, the Ti-in-biotite thermometer of (Henry et al., 2005) and the amphibole-plagioclase thermometers of (Holland & Blundy, 1994). Precision of the Ti-in-biotite thermometer varies with temperature, with the largest 1σ uncertainty reported as ±24°C for greenschist facies conditions, decreasing to ±12°C for amphibolite and granulite facies (Henry et al., 2005). For consistency, we report all Ti-in-biotite temperatures with the larger ±24°C uncertainty. Since all of the rocks we studied are quartz-bearing, we applied both the edenite-tremolite and edenite-richterite thermometers, which both have a 1σ uncertainty of ±40°C (Holland & Blundy, 1994). The application of the edenite-richterite thermometer was restricted to samples with X Ab < 90 because the equilibrium constant for this reaction requires an accurate estimate of the activity of anorthite. In addition to these independent thermometers, for all samples with suitable mineral assemblages, we used the average PT (avPT) and average T (avT) functions in thermocalc(versions tc340i with ds-62 and tc333 with ds-55), which performs a least-squares regression on a set of independent end-member reactions (Powell & Holland, 1994). For consistency with the petrological modelling, end-member activities for mafic samples were calculated using AX 6.2 (Holland, 2018) and using AX (Holland & Powell, 2009) for pelitic and metagranodioritic samples. The 1σ uncertainty for average PT calculations is expressed as an error ellipse, the dimension of which scales with the collective fit of the endmember reactions (Powell & Holland, 1994).

| GSL-18-BD25, Upper-amphibolite facies mylonitic schist
Sample GSL-18-BD25 is a pelitic and migmatitic schist (Figures 2c and 3c). The thin section from this sample comprises garnet (~2%), fibrolitic sillimanite (~1%), biotite (~14%), chlorite (~2%), quartz (~36%), plagioclase (~37%) and K-feldspar (~12%) with accessory rutile and ilmenite ( Table 2). Alignment of biotite, fibrolitic sillimanite, quartz ribbons and boudinaged K-feldspar defines a penetrative foliation ( Figure 3c) and oblique shear-band fabric. Quartz ribbons exhibit a crystallographic preferred orientation under the gypsum plate. The outcrop from which this sample and GSL-18-BD26 were collected has decametre-scale intercalations of mafic and pelitic rocks, with a distinct leucosome and melanosome developed in both lithologies (Figure 3d). On the thin-section scale, leucosomes comprise recrystallized quartz, K-feldspar and highly sericitized plagioclase (white-mica alteration). Melanosomes comprise biotite, sillimanite and retrograde chlorite (Figure 3c). Garnet occurs in the more aluminous melanosomes as porphyroclasts up to ~4 mm in diameter that are partially replaced by biotite along their margins. Analyses of garnet reveal a lack of internal zoning, with the exception of an increase in Mn and Fe in the outermost resorbed rim (Figure 4b). We interpret the core composition of Alm 0.73 Prp 0.18 Sps 0.02 Grs 0.05 And 0.02 to reflect the final stage of diffusional equilibration between garnet and the matrix minerals. There is no significant difference in composition between the matrix biotite and biotite grains adjacent to garnet rims, with Ti = 0.41 apfu and X Mg = 0.44, based on a 22-oxygen formula unit (Table 2). There are no detectable differences in the major element compositions of plagioclase or K-feldspar, independent of grain/porphyroclast morphology. Plagioclase has a mean X Ab = 0.71 and K-feldspar has a X Ab = 0.14 ( Table 2).
0.05 And 0.02 (Figure 4d, Table 2). The 'bell-shaped' grossular and spessartine profile is interpreted as prograde growth, with a local depletion in manganese during initial stages of growth (Kohn & Spear, 2000). Chlorite and biotite embay the outermost rim of garnet, and chlorite fills through-going fractures that were avoided when setting EMPA points (Figure 4e). An increase in Mn-content in the outermost 50 µm rim of garnet is interpreted as Mn-uptake during garnet resorption (Figure 4d). Millimetre-sized porphyroclasts of plagioclase and K-feldspar are rounded and wrapped by a very fine-grained matrix of biotite, chlorite and quartz (Figure 4e). Light-brown biotite occurs as both a very fine-grained matrix phase and intergrown with prismatic sillimanite in millimetre-sized schlieren (Figure 6b). There is no significant difference in biotite composition between the two petrographic positions, with Ti = 0.39 apfu and XMg = 0.40 based on a 22-oxygen formula unit ( extent of complete replacement by radial muscovite in some grains (Figure 6c). Rare, unaltered plagioclase grains have a mean XAb = 0.59, and XAb = 0.11 in K-feldspar.

AND THERMOBAROMETRY
The principle aim of petrological analyses is to determine the pressure and temperature conditions at which the GSLsz rocks underwent their final stage of ductile deformation and metamorphic recrystallization. However, as is commonly the case with garnet-bearing metamorphic rocks, portions of the prograde and/or retrograde metamorphic histories are also preserved (Spear & Selverstone, 1983;St-Onge, 1987). Here, we combine phase diagram-based petrological modelling with conventional multi-equilibria and single-phase thermobarometry to determine the metamorphic histories of the GSLsz rocks.

| Pressure-temperature estimates
6.1.1 | GSL-18-C47-Granulite facies paragneiss Figure 7 presents the calculated equilibria for GSL-18-C47. The diagram is contoured with isomodes of kyanite and garnet compositional isopleths (pyrope and grossular). As a consequence of setting H 2 O to minimally saturate the solidus, melt constitutes a minor phase (<1 modal%) in the immediate supra-solidus field. Only at temperatures above the vapour-absent muscovite dehydration reaction does any significant melting occur, which is accompanied by the growth of an aluminosilicate phase (Dyck et al., 2020). We interpret the aluminosilicate microstructures to have formed along a prograde path that started in the kyanite stability field with 3-4 modal% growth during muscovite melting. The flat compositional profiles in garnets from this sample (Figure 4a) indicate that garnet fully equilibrated with the matrix phases at some point in the rock's history. By comparing the measured garnet chemistry against the calculated isopleths of pyrope and grossular, we estimate that the final equilibration of garnet occurred at ~0.85 GPa and ~785°C. At these conditions, some biotite still remains in the predicted assemblage, and although the assemblage is not fully dehydrated, it straddles the boundary between upperamphibolite and granulite. With a two-pyroxene assemblage reported from nearby mafic lithologies (Hanmer, 1988) and confirmed in our own sample collection, we follow previous workers in mapping this unit as granulite facies. The intergrowth of muscovite and sillimanite indicates that the rock underwent both decompression and cooling following peak metamorphism, remaining in the sillimanite stability field during retrograde re-equilibration. When combined with the Ti-in-biotite results, we interpret that this final stage of deformation and (local) equilibration took place at ~0.7 GPa and ~710°C. We ran two sets of avPT calculations, one with the garnet core composition and a second using the outermost garnet rim composition. Uncertainty ellipses for results of both calculations are presented in Figure 7. Both avPT results are within error of temperature estimates given by Ti-in-biotite. However, we are cautious in interpreting either of the avPT results as being geologically meaningful due to the difficulty in assessing the extent to which matrix biotite and plagioclase re-equilibrated during garnet resorption. 6.1.2 | GSL-18-BD25-Upper-amphibolite facies mylonitic schist Figure 8a presents the pseudosection for GSL-18-BD25 calculated at 0.5-1.3 GPa and 600-900°C. Garnet is only stable at the higher end of the P-T conditions considered. As with the previous sample, we interpret the homogeneous compositional profile in garnet (Figure 4b) to reflect the final stage of diffusional equilibration between garnet and the matrix minerals, at ~1.1 GPa and ~760°C. At these conditions, garnet would have reached equilibrium with a matrix assemblage of melt-biotite-K-feldspar-plagioclase-kyanite. As sillimanite is the only aluminosilicate phase found in the matrix and the temperatures under which both garnet and sillimanite are stable are higher than those indicated by the pyrope and grossular isopleths, we interpret garnet as a porphyroclastic metastable phase. This interpretation is based on the simple assumption that once temperatures were high enough for garnet to effectively maintain equilibrium with matrix phases, the garnet composition would continue to re-equilibrate with increasing temperatures. If this were the case, and the rock had reached a temperature at which garnet was in equilibrium with sillimanite, we would expect to see a garnet composition that was 20%-30% more Mg rich.
Our best estimate for the conditions at which the mylonitic fabric was developed is given by the overlap of the Ti-in-biotite results with the stability field for the observed peak metamorphic assemblage of melt-biotite-K-feldsparplagioclase-sillimanite at ~0.7 GPa and ~710°C. The avPT results superficially appear to broadly overlap the conditions that we interpret as the final stage of equilibrium. However, there is poor agreement between the modelled isopleths and the avPT results generated with garnet core composition. If the garnet rim composition is used for the avPT calculation, the results are considerably less precise (Figure 8a). Leake et al., 1997Leake et al., , 2004 in terms of Si atoms per formula unit and Mg/(Fe 2+ +Mg). Accounting for the overstabilization of augite, the observed mineral assemblage of amphibole-garnet-plagioclasequartz (+melt) is found to be stable over a wide range of upper-amphibolite facies conditions. To further constrain the conditions of peak metamorphism, we plotted the observed garnet and amphibole chemistry on the modelled pyrope, grossular and Si-in-amphibole isopleths (Figure 8b). The pseudosection approach yields peak-metamorphic conditions for GSL-18-BD26 at ~0.9 GPa and ~800°C. There is a considerable overlap in the temperatures predicted by amphibole-plagioclase thermometers and those from the pseudosection results. Furthermore, the avPT results are within error of the results generated by the other methods, albeit shifted towards lower temperatures (Figure 8b). In these calculations, the fluid was considered to be pure H 2 O. Decreasing the activity of H 2 O would reduce the P-T estimates by ~0.01 GPa and ~10°C per 0.1 increment of X(H 2 O).

F I G U R E 5 Amphibole classification plots (after
6.1.4 | 6.1.4. GSL-18-C53-Epidoteamphibolite facies mylonitic schist In rocks for which the equilibrium volume (and the corresponding effective bulk composition) changed alongside the evolution of P-T conditions (i.e. with the chemical isolation of garnet cores), segments of prograde and/or retrograde history may also be preserved (Spear & Selverstone, 1983). In such cases, fractionation (modification) of the effective bulk composition may provide a more accurate estimate of the equilibration conditions (e.g. Evans, 2004;Konrad-Schmolke et al., 2008;Tinkham & Ghent, 2005). The sub-solidus region of the phase diagram constructed for GSL-18-C53 in Figure 9a was generated using the integrated composition and modal proportions of all phases, including the garnet core domain with water considered in excess. A second, fractionated, minimally water-saturated, phase diagram that excludes the garnet core domain was generated to constrain the conditions under which the rim of the garnet grew. For simplicity, we merge the two pseudosections along the water-saturated solidus (Figure 9a). Pyrope and grossular isopleths constrain garnet core growth to ~0.75°GPa and ~550°C, with a matrix assemblage F I G U R E 7 Calculated P-T pseudosection for sample GSL-18-C47. Relevant fields are contoured for pyrope and grossular isopleths and kyanite mode (legend, bottom right). Red polygon represents the garnet core composition (±10% of end-member mole fraction). The diagram is overlain by Ti-in-biotite temperature and avPT error ellipses (1 SD) for calculations including garnet core and rim compositions. The green polygon represents the range of conditions in which the Ti-in-biotite temperatures overlap with the observed stable mineral assemblage (bold text) [Colour figure can be viewed at wileyonlinelibrary.com] of muscovite-biotite-epidote-plagioclase-kyanite-quartzilmenite. The measured garnet rim composition is not replicated anywhere within the unfractionated P-T region of interest. Instead, the pyrope and grossular isopleths generated with the fractionated effective bulk composition constrain the conditions of garnet rim growth to ~0.8 GPa and ~725°C, in equilibrium with the observed peak metamorphic assemblage of biotite-plagioclase-quartz-sillimanitequartz-ilmenite (+melt).
AvPT results for GSL-18-C53, using garnet rim compositions, plot below the considered P-T range with a least square mean of 0.3 ± 0.1 GPa and 488 ± 13°C (correlation coefficient = 0.265 and significance of fit = 0.08). The presence of retrograde chlorite in the matrix shear bands and absence of retrograde staurolite is an indication that, despite continued deformation, equilibrium (on the thin-section scale) was not maintained as the rock cooled. As such, we do not interpret the avPT results as being geologically meaningful. Instead, the retrograde history is loosely constrained by the Ti-in-biotite temperatures and the upper-temperature stability of chlorite. Biotite last equilibrated at 686 ± 24°C in the sillimanite stability field. The alignment of chlorite within the shear fabric indicates the rock continued to deform as the rock cooled below ~600°C, at unknown pressures. 6.1.5 | 6.1.5. GSL-18-C50-Epidoteamphibolite facies mylonitic amphibolite We generated two pseudosections, presented in combination in Figure 9b, at P-T conditions of 0.5-1.1 GPa and 400-800°C for sample GSL-18-C50 because it exhibits petrological evidence for significant hydration following peak metamorphism. One set of equilibria relations was calculated with H 2 O as an excess phase and was used for evaluating sub-solidus phase stability (blue field shading). We used a second pseudosection (red field shading), constructed with H 2 O set to minimally saturated at the 0.7 GPa solidus, to evaluate the stability of garnet in melt-bearing assemblages. The results from both calculations were merged along the H 2 O-present solidus and are shown in Figure 9b.
If we consider the close proximity of this sample to GSL-18-C53, and assume similar peak metamorphic conditions for both samples, then the peak assemblage of GSL-18-C50 would have garnet in equilibrium with amphibole-plagioclase-Kfeldspar-quartz-ilmenite-titanite (+melt). While there is no garnet remaining in this sample, we interpret the rounded shape and size of the porphyroclastic microstructure (Figure 6e) as an indication that garnet was once stable and persisted during the initial development of mylonitic fabric. The results from the edenite-tremolite thermometer further indicate that this rock F I G U R E 8 Calculated P-T pseudosection for samples (a) GSL-18-BD25 and (b) GSL-18-BD26. Note pressure range differs between pseudosections. In both diagrams, relevant fields are contoured for pyrope and grossular isopleths and Si apfu in amphibole (legend, below each diagram). Bold text marks the observed stable mineral assemblage. Red polygons represent the garnet core composition (±10% of end-member mole fraction). Diagram (a) is overlain by Ti-in-biotite temperature and avPT error ellipses (1 SD). The green polygon represents the range of conditions in which thermobarometry estimates overlap with a peritectic K-feldspar-bearing field. Diagram (b) is overlain by the results of the Holland and Blundy (1994) amphibole-plagioclase thermometers and avPT calculation (1 SD). The yellow polygon represents the range of conditions in which the observed amphibole composition is bracketed by amphibole-plagioclase thermometry results [Colour figure can be viewed at wileyonlinelibrary.com] experienced temperatures in excess of ~650°C. However, the amphibole and plagioclase compositions may have equilibrated during retrograde metamorphism and deformation.
It was not possible to determine the extent to which this rock underwent decompression prior to the growth of the retrograde assemblage because the observed assemblage of amphibole-albite-clinozoisite-biotite-titanite-quartz-rutile is stable over a broad pressure range (0.4-1.0 GPa). Similarly, without garnet, there is an insufficient number of independent pressure-sensitive end-member reactions to calculate an avPT result. Instead, we made an avT calculation (excluding sanidine and anorthite), which yields a mean temperature estimate of 503 ± 16°C (significance of fit = 1.21) for reactions between amphibole-albite-clinozoisite-chlorite-ilmenite-quartz-H 2 O. 6.1.6 | GSL-18-BD7-Greenschist facies ultramylonite Figure 10 presents the phase equilibria for sample GSL-18-BD7 at P-T conditions of 0.2-0.7 GPa and 400-650°C. The overlapping stability of chlorite-epidote-K-feldspar±biotite, all found in the finely recrystallized matrix of GSL-18-BD7, constrains the final temperature of deformation to ~415-445°C. This temperature estimate is independently corroborated by an avPT calculation, based on equilibrium between chlorite-plagioclase-amphibole-epidote-quartz-H 2 O, which yields 0.4 ± 0.19 GPa and 442 ± 65°C (correlation coefficient = 0.986 and significance of fit = 1.13). When we combine the results from phase equilibria and avPT calculations, our best estimate for the conditions of greenschist facies metamorphism in the GSLsz is ~0.4 GPa and ~430°C.

| P-T synthesis and implications for
Slave-Rae tectonics Figure 11a summarizes the P-T results presented in this paper, overlain on a range of crustal thermal gradients typical of regional metamorphism (Brown, 2007; Thompson, 1984). Segments of clockwise P-T loops are recorded in samples from the granulite (GSL-18-C47), upperamphibolite (GSL-18-BD25) and epidote-amphibolite (GSL-18-53) units. Although there is a considerable range of final conditions of equilibrium for individual map units, spanning temperatures of ~430-800°C, the granulite, upperamphibolite and epidote-amphibolite samples all experienced peak metamorphism at similar P-T conditions (~0.85 GPa and ~750°C). The prograde record that is preserved in garnet chemistry and aluminosilicate microstructure indicates that metamorphism of the GSLsz rocks initially progressed under a thermal gradient of ~650°C/GPa. In four of the six samples that we studied, there is a record of decompression and retrograde mineral growth associated with the development of a penetrative shear fabric. The P-T conditions of the final stage of equilibrium recorded by matrix minerals in all six samples collectively define a metamorphic field gradient of ~1,000°C/GPa.
The nested nature of the clockwise prograde P-T paths preserved in the GSLsz rocks is consistent with a single metamorphic event involving crustal thickening and Barrovianstyle metamorphism (England & Thompson, 1984;Weller et al., 2013). The absolute timing of crustal thickening and the attendant development of metamorphic fabrics remain poorly constrained. However, in samples that we studied, it appears that crustal thickening either preceded or was coeval with ductile shear. This interpretation is, in part, supported by the absence of sigmoidal inclusion trails in garnet, which is consistent with crustal thickening preceding the development of a penetrative shear fabric. The one exception we note is the slightly sigmoidal garnet inclusion trails found in sample GSL-18-BD26 (garnet amphibolite). However, this sample only exhibits a weak and non-penetrative shear fabric and so it likely remained unaffected by the dextral shear. Moreover, at the regional scale, it is unlikely that the lower-crustal metamorphic and deformation features of the GSLsz would have been isostatically uplifted and exhumed to be present at the current surface erosion level without an earlier stage of crustal thickening.
Given the present configuration of the Slave and Rae cratons, the most apparent cause of the crustal thickening is convergence between Slave and Rae. However, without reasonable constraints on the timing of the prograde metamorphism and the provenance of sedimentary material, it is premature to speculate on the pre-GSLsz configuration of Slave and Rae. Moreover, the timing of the postulated transition from crustal thickening to lateral shear on the Slave-Rae margin is poorly documented. Although, the elevated F I G U R E 1 0 Calculated P-T pseudosection for sample GSL-18-BD7. The green polygon represents the predicted stability field for the observed mineral assemblage chlorite-plagioclase-epidote (±zoisite)-amphibole-titanite±biotite. Diagram overlain by avPT error ellipse (1 SD) [Colour figure can be viewed at wileyonlinelibrary.com] metamorphic field gradient and superposition of GSLsz fabrics on peak metamorphic assemblages implies that the ductile deformation associated with the GSLsz was active while the overthickened northern Rae margin was still relatively hot.

| Metamorphic evolution of continentalscale shear systems
In the Sibson-Scholz model of fault zone structure, illustrated schematically in Figure 11b, one or more zones of dominantly frictional deformation in the upper crust transition into zones of dominantly viscous deformation in the middle and lower crust, with the thickness of the deformed zones generally increasing with depth (Scholz, 1988;Sibson, 1983). Activity of a fault zone during exhumation by erosional unroofing commonly results in an overprinting sequence of down-temperature fault rocks (Hanmer, 1988;Jefferies et al., 2006;Norris & Toy, 2014;Phillips & Searle, 2007;Sibson, 1977). If deformation occurs over time-scales shorter than required for conductive relaxation of elevated crustal thermal gradients (<10 Ma for continental collision (England & Thompson, 1984), <100 Ma for margins with continental arcs; Xiao et al., 2005), then the conditions of final equilibrium in all levels should reflect a similar geothermal gradient.
Shear heating associated with ductile deformation in the GSLsz is an additional process that might have impacted the recorded temperatures. The magnitude of that contribution to the heat budget is difficult to assess because of the simultaneous exhumation that accompanied deformation, which resulted in cooling at an unknown rate. Previous estimates of shear heating in major ductile shear zones without exhumation suggest shearing may elevate temperatures by 10s (e.g. Thacher and England, 1998) to 100s (e.g. Nabelek et al., 2010) of degrees. However, models that estimate significant shear heating (e.g. Nabelek et al., 2010) require stresses during ductile deformation that are considerably higher than those estimated by palaeopiezometric analyses of similar shear zones (e.g. Chatzaras et al., 2015;Behr and Platt, 2014).
The improved P-T constraints provided by the present study are a significant development in discriminating between the GSLsz developing as a single structure active during one progressive deformation event, or a structure that was reactivated at different crustal levels during temporally distinct deformation events. The conditions of final equilibration, which fall on a single metamorphic field gradient despite spanning greenschist facies to granulite facies conditions (Figure 11a), are consistent with all the mylonitic units having formed during a single portion of the orogenic cycle, within a sufficiently short duration of time that there would have been little change to the crust thermal gradient. Furthermore, the retrograde path of granulite facies sample GSL-18-C47 indicates that individual rock packages did indeed follow this apparent geothermal gradient during exhumation and cooling from a maximum burial depth of ~30 km (Figure 11a).

| Quantifying P-T in shear zones
The application of equilibrium thermodynamics to rocks that have experienced a protracted history of shear and dynamic recrystallization poses some significant challenges that warrant further discussion. While a range of thermobarometry and petrological modelling tools can be employed to quantify portions of the P-T history for a rock, correctly interpreting the P-T results requires knowing how equilibration volume evolves with both metamorphic grade and degree of strain (Marsh et al., 2009). One of the main challenges we faced in this study was determining mineral assemblages that chemically equilibrated together during the final stage of deformation and recrystallization. This is particularly important in samples with porphyroclasts because (a) the partitioning of strain may trigger recrystallization and equilibration in the matrix phases but not the porphyroclastic phase (Koons et al., 1987) and (b) the length scale of intracrystalline diffusion are greater in the porphyroclasts than the finer-grained matrix, which may result in porphyroclast cores that do not equilibrate with the matrix (Florence & Spear, 1995).
In shear zones, in which deforming matrix phases can more efficiently maintain chemical equilibrium (Koons et al., 1987), the application of thermobarometric equations that include garnet endmembers, whether it be by avPT, TWQ (Berman, 1991), or another program, should be used with great caution on samples that show evidence of retrograde equilibration. For example, the presence of locally manganese-rich garnet rims likely indicated rim resorption processes and intragranular transport (Tracy, 1982) and should be avoided. This point is well illustrated in the samples that we studied, for which the disparity between the P-T results calculated by conventional thermobarometers and pseudosection modelling reflects disequilibrium between porphyroclastic garnet and matrix phases. We base this interpretation on the garnet compositional profiles for GSL-18-BD47, GSL-18-BD25, GSL-18-C53, all of which reveal an enrichment of manganese in the outermost garnet rims. In these cases, the issue with using garnet rim compositions is that garnet porphyroclasts did not efficiently diffusively equilibrate with the matrix minerals. Instead, garnet was consumed, resulting in the redistribution of its components (primarily almandine) into the matrix, leading to more iron-rich biotite. The effect of disequilibrium between garnet and the matrix is best illustrated in GSL-18-C53, in which the garnet isopleths (pseudosection approach) predict peak metamorphism at ~0.8 GPa and ~725°C, whereas multiphase thermobarometry (avPT) yielded a result of ~0.3 GPa and ~490°C. We consider the P-T results predicted by the pseudosection modelling for sample GSL-18-C53 to be more accurate than those predicted by avPT for two reasons. First, the pseudosection results are consistent with the observed mineral assemblages, notably the presence of sillimanite and melt (manifest as plagioclase-quartz-K-feldspar leucosome). Second, the thermal gradient from the pseudosection results (~900°C/GPa as opposed to ~1,625°C/GPa for the avPT result) is consistent with the record of crustal thermal gradient preserved in the other GSLsz samples.
Based on the findings of this study, we have the following recommendations for quantifying P-T in crustal-scale shear zones. Equilibria modelling remains the most robust approach for determining the prograde and peak metamorphic conditions. When garnet is present, comparing measured compositions with modelled isopleths of pyrope and grossular is particularly useful for constraining peak conditions. However, in rocks in which garnet is chemically homogenized by intracrystalline diffusion and there is little control on the extent of equilibration between garnet and matrix phases, determining absolute prograde and peak P-T conditions remains challenging. Determining the final conditions of shear and recrystallization is altogether more difficult. In samples with biotite-rich shear bands and folia, Ti-in-biotite thermometry appears to record the near-final conditions of ductile shear. We found Ti-in-biotite temperatures to be systematically lower than those predicted by phase equilibria modelling. We believe this discrepancy to be the result of strain-induced recrystallization of biotite as strain was localized to biotite-rich shear plans during exhumation and cooling. In amphibole-bearing compositions, we found that amphibole-plagioclase thermometry provides reasonable and consistent estimates for matrix equilibrium. At amphibolite facies conditions or lower, it is unlikely that significant diffusion would occur without recrystallization due to coupled Si-Al substitution in plagioclase (Grove et al., 1984). As such, amphibole-plagioclase temperatures in sheared rocks are most likely to record the final temperature at which plagioclase dynamically recrystallized. In general, unless retrograde metamorphism was accompanied by thorough recrystallization (as with GSL-18-BD7), it is challenging to accurately model retrograde assemblages because the effective bulk composition is unknown for this portion of the rock's metamorphic history. As such, when determining the retrograde conditions, results from phase equilibria modelling should only be used as a general guide to the sequence and broad conditions of retrograde reactions.
We were fortunate to have access to paired mafic and pelitic samples from the upper-and epidote-amphibolite facies units. In both cases, complementary segments of the GSLsz metamorphic history were preserved in the different lithologies. Given the consistency of values generated by the independent thermometers, avPT, and pseudosection modelling for sample GSL-18-BD26, we have confidence in our results. F I G U R E 1 1 (a) Summary of all P-T results reported in this study. Black circles represent the best estimates for the last recorded equilibrium conditions (with 1 SD error bars). Linear geotherms (dashed lines) intercept the surface at 25°C. The results suggest crustal thickening and Barrovian-style prograde metamorphism progressed under a geotherm of ~600°C/GPa. The last recorded equilibrium conditions yield a metamorphic field gradient of ~1,000°C/GPa. (b) Schematic cross-section of the GSLsz. Shear zone width is based on surface geology with the depth of each unit calculated assuming an average overburden density of 2,750 kg/m 3 [Colour figure can be viewed at wileyonlinelibrary.com] However, there is a ~40°C difference in the peak temperature recorded by GSL-18-BD26 (mafic) and the neighbouring GSL-18-BD25 (pelite). This discrepancy may be a consequence of using different data sets and activity-composition relations for the mafic and pelitic samples.
Finally, there is a wide range of independent single-phase thermometers that depend on the temperature-dependent substitution of minor or trace elements in metamorphic phases that may record additional information about the metamorphic evolution of shear zones. For example, Zrin-rutile (Ferry & Watson, 2007;Watson et al., 2006) and Zr-in-titanite (Hayden et al., 2008). When applying these thermometers to sheared rocks, it is important to assess whether or not the measured minerals recrystallized during shear. There are examples of titanite grains found in ductile shear zones that yield U-Pb ages that are associated with protolith crystallization rather than metamorphic recrystallization (Kohn & Corrie, 2011;Spencer et al., 2013). In these cases, if the isotopic systematics were not reset by shear (through the removal of incompatible Pb atoms), it seems unlikely that Zr would equilibrate with matrix phases and the Zr-in-titanite would reflect the temperature of shear. One thermometer that we did not use in this study that may be particularly useful for constraining the final temperatures of deformation in shear zones is the Tiin-quartz thermometer of Wark and Watson (2006). While we acknowledge the Ti-in-quartz thermometer has its own limitations (e.g. Kendrick & Indares, 2018), we intend to explore the application of Ti-in-quartz to the GSLsz rocks in a follow-up study.

Funding
Funding for this research was provided by a Natural Sciences and Engineering Research Council of Canada Discovery Grant to Dyck, Natural Environment Research Council Environmental Doctoral Training Grant to Goddard and Utrecht University start-up funds to Wallis. Additional fieldwork funding and support was offered by University College Oxford, The Geological Society and the Northwest Territories Geological Survey (NTGS contribution #0122). We thank Kyle Larson and Lorraine Tual for insightful reviews and Julia Baldwin for editorial handling.

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section. Appendix S1. EMPA operating conditions. Figure S1. Supporting field photos and photomicrographs. Figure S2. Phase maps. Figure S3. T-X H2O phase diagram for sample GSL-18-BD25.