The Deep Structure and Rheology of a Plate Boundary-Scale Shear Zone: Constraints from an Exhumed Caledonian Shear Zone, NW Scotland

Contractional fault systems that cut the crust and upper mantle are necessary consequences of plate convergence and continental collision [1]. At any one time, much of the relative plate motion is localized onto relatively narrow zones of high strain [2]. In the upper crust, these zones are manifested as discrete brittle faults or systems of faults that exhibit a frictional rheology. Below the brittle-to-ductile transition (BDT), which delineates the shift from dominantly brittle and frictional behavior to deformation dominated by crystal-plastic and diffusional processes, they are thought to widen into broad zones of distributed strain commonly referred to as ductile shear zones. In plate boundary-scale settings, these shear zones may be upwards of 20 km wide (depending on fault regime and geothermal gradient) for the quartz-rich crystalline continental crust [3–8] and even wider in the upper mantle [9].

Studies on the San Andreas system in California suggest that some strike-slip faults continue as discrete, narrow fault zones down to the Moho [10–12], whereas SKS-splitting and xenolith data suggest that the San Andreas transform system may form a ~100 km wide shear zone in the upper mantle [13]. A 1-2 km thick zone of mylonites, derived from mid-crustal depths, has been described along the Alpine Fault in New Zealand, suggesting a mid-crustal shear zone of at least these dimensions [14]. These two systems exemplify the significant uncertainty concerning the deep structure and geometry of plate boundary-scale fault systems.

Despite these and other studies that have made significant progress in understanding the deep roots of fault systems below the BDT, key information is still missing. This includes (1) the variation of stress as a function of depth in contractional and strike-slip fault systems, although considerable progress has been made on normal fault systems (e.g., [15]); (2) the geometry, thickness, and internal structure of shear zones as they transect the lithosphere; and (3) the mechanical properties (i.e., rheology) of these zones, including variability in space and time.

Developing a naturally-constrained model for ductile shear zone geometry and rheology has broad implications for a better understanding of the processes that control plate interaction and faulting. Lithosphere-scale faults likely behave as systems with decoupling or complex feedbacks between regions where behavior is dominated by transient, stick-slip, frictional events in the upper crust and steady-state-dominated, ductile behavior in the mid and lower crust. Rocks around the BDT probably retain the majority of crustal strength and therefore play a significant role in the loading and activation of faults near the surface [16, 17]. As such, it is critical to characterize what happens at depth, within and below the BDT, to better understand how the crust behaves within the brittle realm. Ultimately, field-based observations on natural systems are required to construct models and provide validation for experimental research on the rheology of crustal rocks.

Large displacements on plate boundary shear zones, and the complexities of the processes by which they are exhumed, can make it difficult to reconstruct the structure of such shear zones with depth. In the examples outlined above, both are strike-slip systems, meaning that barring oblique slip or postfaulting uplift and erosion, rocks deformed at mid- to lower-crustal levels are unlikely to become exposed. In a normal-sense fault system, footwall rocks are cooled and exhumed with motion along the fault system, often preserving a zonation in dynamic microstructures that records strain localization and the spatial and temporal evolution of the shear zone. Unlike normal-sense systems, reverse-sense systems generally bury rocks, advecting cold rocks downwards in the footwall, and thickening the crust. For rock microstructures that record a reverse-sense evolution to be preserved, hanging wall rocks must undergo syn-deformational exhumation to prevent static overprinting of dynamic microstructures.

Here, we present an integrated field, microstructural, and analytical study from the Caledonian frontal thrust system in northwestern Scotland to better understand the geometry, structure, and rheology of a lithosphere-scale contractional shear zone from the BDT to the lower crust. To do this, we present estimates of differential stress (σ_d), pressure (P), and temperature (T) of deformation, along with field and microstructural observations from a continuous shear zone profile to reconstruct shear zone rheology, internal structure, and geometry.

The Caledonian orogen extends from western Ireland north to Svalbard and records a collisional event that was likely of Alpine or Himalayan dimensions [18, 19], spanning a period of nearly 200 Ma from the late Cambrian to latest Devonian [20]. Broadly speaking, the Caledonian orogeny is the result of the closure of the Lower Paleozoic Iapetus Ocean accommodated along several intraoceanic subduction zones, and subsequent collision of Laurentia, Baltica, Avalonia, as well as oceanic arcs and other minor terranes [18, 21–23].

In northwestern Scotland, two main phases of the Caledonian orogeny are expressed: the Grampian (475–460 Ma) and the Scandian (445–425 Ma) ([24] and references therein). The Grampian phase, attributed to the collision of an oceanic arc, marked the end of the eastern Laurentian passive margin and deposition of Cambrian to Ordovician shelf facies rocks, culminating in the emplacement of ophiolites and Barrovian-style metamorphism [21]. The Scandian phase, during which many of the ductile thrusts of the orogenic prowedge were formed, was a result of the collision of Baltica with Laurentia (c. 435-420 Ma), followed closely by Avalonia (c. 425 Ma) [18, 22–24].

The foreland to the Caledonian orogenic belt comprises an Archean to Mesoproterozoic basement gneiss complex unconformably overlain by Proterozoic and Cambro-Ordovician sedimentary rocks. The basement Lewisian complex includes multiply-deformed gneisses comprising 3000–2700 Ma meta-granitoids, basites, and sediments that were affected by the Badcallian (c. 2450 Ma) and later Laxfordian (c. 1800 Ma) orogenic events [25]. Rocks of the Lewisian complex are unconformably overlain by Proterozoic (c. 1200–1050 Ma) continental red beds of the Torridonian Group and Cambro-Ordovician Laurentian shelf sequence rocks including the Eriboll Formation (quartzites), An t-Sron Formation (fine to medium-grained clastic rocks), and carbonate rocks of the Durness Group. Lewisian gneisses and Eriboll Formation quartzites are the only foreland rocks that crop out within the study area. Lewisian gneisses comprise an assemblage of feldspar (orthoclase + minor albite) + quartz + white mica + epidote group minerals + chlorite + opaques. Eriboll Formation quartzites that are unaffected by ductile strain have equant quartz grains with grain sizes up to >1 mm. Minor components of feldspar + white mica are also present.

The retrowedge metamorphic hinterland of the Caledonian orogenic belt is dominated by rocks of the Moine Supergroup, tectonically interleaved with basement gneisses correlative with the Lewisian, and intruded by Caledonian granitoids (Figure 1). In northwestern Scotland, hinterland rocks make up four major ductile thrust sheets; from structurally lowest to structurally highest these are the Moine, Ben Hope, Naver, and Skinsdale Nappes. These ductile thrust sheets preserve internal gradients in metamorphic grade that increase structurally up from lower greenschist facies in the structurally lowest parts of the Moine Nappe to migmatitic upper amphibolite facies in the upper Ben Hope, Naver, and Skinsdale Nappes [26–28]; this field gradient is interpreted to be related to west-directed (present orientation) Scandian thrusting [29].

The Moine Supergroup comprises strongly deformed Proterozoic metasediments consisting chiefly of psammites and pelites with minor marble and metaconglomerate [30, 31]. Rocks of the Moine Supergroup are likely more distal, time-equivalents of Torridonian Group rocks that crop out in the Moine Thrust Zone and Caledonian foreland [32]. The broad metasedimentary Moine Supergroup (herein referred to as Moine schist) is traditionally subdivided into three distinct groups: the Morar, Glenfinnan, and Loch Eil, although in the present study area, all Moine schist is included within the Morar group. Rocks of the Morar group are dominantly psammitic in composition with a metamorphic assemblage of quartz + white mica + feldspar (albite + orthoclase) + epidote group minerals ± biotite ± garnet ± opaques with accessory phases including titanite, zircon, apatite, and monazite. Pelitic layers, richer in micas, feldspar, garnet, epidote, and rare staurolite crop out in lenses, although none are exposed in the study area. Chlorite occurs throughout the sequence, generally as a retrograde breakdown product from biotite and/or garnet. Plagioclase is present as relict sedimentary porphyroclasts and smaller grains incorporated into the matrix, whereas potassium feldspar is chiefly incorporated into the matrix but occurs as porphyroclasts at lower structural levels. Garnets, where present, display varying morphologies from the euhedral to anhedral [33] and skeletal. Meter-scale thick lenses of amphibolite (hornblende + quartz + feldspar ± garnet ± white mica ± biotite ± titanite ± epidote group minerals) crop out at the base of the Ben Hope Nappe and as rare inliers within the Moine schist.

Scandian deformation in NW Scotland can be traced from the Moine Thrust Zone (MTZ), an imbricated zone delineating the transition from the foreland to the hinterland of the Caledonian orogenic belt, eastwards for up to 40 km into the hinterland where deeper, mid- to lower-crustal structural levels are exposed. Exposure of the MTZ and overlying hinterland rocks extends approximately N-S for at least 200 km (Figure 1) and possibly greater than 450 km when considering suggested offshore continuations to the north and south [34]. The MTZ forms a belt of imbricated foreland rocks bound by the structurally lowest Sole Thrust below and Moine Thrust above.

Timing of thrusting is based on dating of intrusive rocks and deformational fabrics; K-Ar and Rb-Sr dating of syn-kinematic micas in mylonites [35] and U-Pb dating of zircons in syn-kinematic intrusions [29, 36] constrains activity to between ~435–425 Ma. This timing is consistent with reported Rb-Sr ages on muscovite and ⁴⁰Ar/³⁹Ar cooling ages from hornblende and muscovite [37]. Thrusting was roughly WNW-directed, recorded by the pervasive stretching and mineral lineation and tight to isoclinal folds rotated into parallelism with the lineation [38–40].

This study examines rocks exposed along a transect that extends from the southern head of Loch Eriboll at the settlement of Strabeg, in the footwall of the Moine Thrust (MT), SSE across the southern continuation of the Eriboll Peninsula, and the Ben Hope Thrust to the summit of Ben Hope (hereon referred to as the Strabeg transect; Figures 2(a) and 2(b)). This area was chosen based on the clear lithologic expression of the MT (Figure 3), the relatively simple tectonic stratigraphy (Figure 4), and extensive, continuous exposure along much of the transect. The transect is taken along a line trending 110°-290°, parallel to the average trend of the stretching lineation and inferred direction of transport [38–40]. In total, the transect spans 9 km on the ground, corresponding to a structural thickness of over 3 km. The lower structural levels of the transect record greenschist grade metamorphism, while the structurally higher levels record metamorphism at lower- to mid-amphibolite facies [43–45].

Three primary ductile thrusts crop out along this transect (Figures 1, 2, and 4(a)). The structurally lowest Lochan Riabhach Thrust (LRT) emplaces Lewisian gneiss over the foreland Eriboll Formation quartzite. Structurally above, the Moine Thrust (MT) emplaces Moine schist of the Moine Nappe onto variably mylonitized Lewisian gneiss and Eriboll Formation quartzite. Exposed along the west flank of Ben Hope, the Ben Hope Thrust (BHT) emplaces Moine schist, along with thin layers of basement gneiss (likely Lewisian) and amphibolite, which make up the Ben Hope Nappe, over Moine schist of the Moine Nappe (Figure 2).

Normal faults have been previously mapped striking N-S from the head of Loch Hope (Figures 1 and 2 and references therein). Due to a lack of exposure in this section of the transect, these cryptic structures have an unknown orientation and displacement. Because there is no distinct gap or repetition of lithology or microstructure across these structures, we conclude that the offset along these faults and any resulting discrepancy in structural distance across them is minimal.

The sequence of thrusting generally propagated towards the foreland with earlier, hinterland thrust sheets transported “piggy-back” on younger thrusts [34, 38, 46]. In detail, and especially at the orogen-scale, the system likely had a more complex evolution with out-of-sequence thrusting and reactivation of older thrusts [47, 48]. Within the study area, however, fold geometries and interactions are consistent with a locally simple foreland-propagating model [49]. Additionally, the MT is the primary tectonic boundary between the Moine schist and foreland Cambro-Ordovician and Lewisian rocks, and it is therefore interpreted to be the approximate basal detachment of this broader crustal-scale shear zone system at depth, although foot wall rocks are involved in deformation [26]. Lastly, geochronological evidence indicates that the Scottish Caledonides were undergoing active erosion and cooling throughout much of the Scandian phase (syn-deformational cooling and erosion is suggested ca. post-425 Ma; [50] and references therein). We support these ideas and proceed with the interpretation that faults exposed along the Strabeg transect are part of the same shear zone at depth, which, within the study area, repeatedly propagated towards the foreland. We therefore treat the BHT, MT, and LRT as representing different structural levels of the Scandian shear zone (Figure 5). The hinterland thrust sheets that were uplifted piggy-back on these thrusts were simultaneously exhumed, preserving the rock microstructure and exposing deep levels of the same shear zone (Figures 5 and 6).

Structural distance is taken as the foliation-perpendicular distance from the projected MT plane to the point of interest. Present-day structural distances along the Strabeg transect are calculated in order to place our samples in the context of the shear zone as a whole, and to provide a basis for reconstructing the geometry of the shear zone when it was active. Distances are calculated using the lithologic expression of the MT (i.e., the contact between Moine schist and underlying Lewisian gneiss) as the reference plane (Figure 4(a)). Although the lithologic contact likely does not coincide with the most localized, highest strain rocks, it provides an appropriate passive boundary within the shear zone.

For distance measurements, the MT is assumed to extend as a planar feature at a constant dip beneath the area of interest. This is supported by the DRUM, MOIST, and LISPB seismic studies, which image reflectors interpreted to be thrust planes cutting through the upper and middle crust at a constant angle [51–53]. We assume that the thrust plane dips at 15°ESE, based on measurements in the field and published British Geological Survey mapping (Loch Eriboll sheet 114 W). We also include a correction for the change in elevation between the point of interest and the trace of the MT.

For structural distance measurements to be valid, we assume that the foliation is parallel or subparallel to the shear plane. This is based on the fact that the foliation in a shear zone rotates towards the shear plane with progressive strain. We measured strain intensities on two outcrops of basal Moine metaconglomerate (Figure 1) which resulted in values of $D=2.5$ and $D=2.3$⁠, where $D$ is a dimensionless measure of strain magnitude defined on a logarithmic Flinn diagram (Figure S3). Additional attempts to quantify strain using the nearest neighbor technique on feldspar porphyroclasts within the Moine schist were unsuccessful, likely due to the invalid initial condition assumption of a random distribution of clasts. Despite the lack of widespread quantitative finite strain estimates within the Moine schist, we posit that finite shear strains are generally large enough to approximate the foliation as being parallel to the shear plane, thus validating our calculations of structural distance (see [26, 54] for further discussion).

Based on the examination of >50 thin sections from the Strabeg transect, we have identified four distinct yet gradational microstructural domains (Figures 4(a) and 4(c)–4(g)). These domains record progressive shear zone evolution and strain localization and are subdivided based on quartz recrystallization mechanism, rheological behavior (e.g., phases accommodating strain) and to a lesser extent, lithology. In the following section, we describe microstructures and defining characteristics for each of these domains, with emphasis on quartz microstructure and recrystallization mechanism.

Across the transect, quartz recrystallizes by bulge nucleation (BLG), subgrain rotation (SGR), and high-temperature grain boundary migration (high-T GBM). BLG, SGR, and GBM, as outlined in Stipp et al. [55], broadly correspond to the experimental regimes 1, 2, and 3 of Hirth and Tullis [56], respectively. The reader should take note however that regime 3 of Hirth and Tullis [56] is better related to a transitional SGR/GBM zone of Stipp et al. [55]. For the remainder of this manuscript, we will adopt the nomenclature outlined in Stipp et al. [55].

Within the Eriboll quartzite and the tectonically overlying Lewisian gneiss, Domain 1 records deformation structurally below the MT and in rocks affected by movement along the Lochan Riabhach Thrust (LRT). Domain 1 spans a structural distance of approximately 250 m, extending from the lower limit of crystalplastic deformation in the footwall of the LRT to approximately 100 m structurally below the MT (Figure 4(a)). Strain is accommodated chiefly by crystal plastic deformation of quartz and mica as well as subsidiary brittle fracture in quartz, mica, and feldspar. Dynamic recrystallization in quartz is dominated by bulge nucleation, indicated by serrated grain boundaries and bulges, as well as very limited subgrain rotation, evidenced by the development of core and mantle microstructures with mantling subgrains roughly equal in size and dimension to recrystallized grains and an abundance of low-angle misorientations (Figures 4(g) and 7(a)–7(c), S2). Percent of recrystallized quartz increases structurally up from the undeformed Eriboll Fm. in the footwall of the LRT, where detrital grains are preserved at the lowest structural levels, to complete recrystallization at <20 m below the lithologic MT. Where present, feldspar generally forms large (500 to >1000 μm) porphyroclasts which show variable intensity of brittle fracture typically increasing structurally up towards the MT. Fractures in quartz and feldspar are commonly filled with precipitated quartz and/or calcite, indicating solution-precipitation creep involving these phases (Figure 7(c)). Phyllosilicates, including white mica and chlorite, are the primary phases defining the foliation and generally wrap around larger feldspar porphyroclasts.

The characteristics differentiating Domain 2 from Domain 1 are the near-complete dynamic recrystallization of quartz, chiefly by SGR, and extensive reduction in grain size. We subdivide Domain 2 into two subdomains based on foliation intensity, grain size, and interconnection of phyllosilicate layers (Figures 7(d) and 7(e)). Domain 2a (Figures 4(f) and 7(d)), characterized by finer grain sizes and strong lamellar folia of alternating quartz-rich and homogenized, fine-grained, polyphase (quartz + feldspar + white mica ± chlorite ± minor accessory) material, extends from approximately 100 m structurally below the MT to >100 m structurally above the MT (Figure 4(a)). Domain 2b, extending structurally up to approximately 200 m above the MT, is characterized by increased grain size and anastomosing phyllosilicate networks that commonly wrap around feldspar porphyroclasts (Figures 4(e) and 7(e)). Matrix material in Domain 2b records less mixing (homogenization) of fine-grained quartz + feldspar + mica compared to Domain 2a (e.g., Figures 7(e) and 7(f)). Recrystallized quartz-dominated lamellae in Domain 2b also show variable buckle folding whereas quartz-dominated lamellae in Domain 2a show no evidence of folding and are mostly parallel to subparallel to the macroscopic foliation (Figures 7(d) and 7(e)). In all, Domains 2a and 2b extend for a structural distance of approximately 300 m above the MT (Figure 4(a)).

Within Domain 2, quartz textures record pervasive dislocation-driven dynamic recrystallization chiefly by SGR, indicated by (1) development of core and mantle microstructures with mantling subgrains roughly equal in size and dimension to recrystallized grains; (2) development of a pervasive oblique shape preferred orientation (SPO) of recrystallized grains; and (3) misorientation angle distributions within dynamically recrystallized quartz (Figure S2). Limited evidence in the form of quartz-filled pressure shadows in porphyroclasts and rare quartz fill in fractured garnets and feldspar porphyroclasts record a relatively minor component of dissolution-reprecipitation creep (i.e., pressure solution; Figure 7(e)). Where quartz is mixed with other phases, generally mica and feldspar, the resulting mixture is significantly finer-grained due to the pinning of grain boundaries (Figures 7(d) and 7(e)). Because of the small grain size and abundance of phase boundaries that may act to promote quartz dissolution rates [55, 56], quartz in regions dominated by fine-grained, polyphase material (Figures 7(d) and 7(e)) likely underwent increased deformation by pressure solution, although we lack direct evidence to support this. In pure quartz domains where grain boundaries are unaffected by nonquartz phases, grains show a ubiquitous SPO generally inclined at 10-50° to the foliation in the direction of shear (Figures 4(e) and 7(f)).

Feldspars occur as equant to elongate clasts with diameters of up to >1000 μm. Feldspar long axes are parallel to the foliation and clasts are commonly wrapped by micas and other phyllosilicates (e.g., chlorite). There is a correlation between feldspar shape and surrounding material; whereas feldspar grains surrounded primarily by quartz tend to fracture, feldspar grains surrounded by a polyphase matrix are elongate (Figure S1). The presence of elongate clasts in a polyphase matrix is suggestive of crystalplastic deformation or spallation of the feldspar; however, there is only limited direct evidence (e.g., core and mantle structure, bulging, or sutured grain boundaries) to support crystal plastic deformation.

Discontinuous microscale shear bands, inclined to the primary foliation, are common in Domain 2b (Figure 8(a)). These appear to be coeval with the primary foliation: there is no difference in microstructure in the bands, and they are defined primarily by the deflection of the mica foliation. The relative abundance and intensity increase structurally up from very high strain rocks in which shear bands have been obliterated to a maximum at a structural distance of 200-300 m above the MT. Structurally above this, both intensity and abundance decrease until the transition to Domain 4, above which shear bands are generally absent. Most well-developed shear bands indicate a top-WNW sense of shear although rare, less developed bands indicate top-ESE, consistent with a strain geometry dominated by simple shear but with a minor component of flattening.

The structurally higher part of Domain 2 contains a zone of Lewisian gneiss <4 m thick, in contact with Moine schist structurally above and below. This Lewisian inlier is likely either an imbricate thrust slice, in which case the lower contact is a fault and the upper contact is a nonconformity, or the core of a tight to isoclinal fold which incorporated basement rock, in which case both the upper and lower contacts are nonconformable. We observe no exposure of a brittle fault or high-strain zone. At the microscale, the grain size, shape, and recrystallization mechanism in recrystallized quartz do not vary drastically across this region. Below, we address the rheologic implications of this Lewisian inlier with respect to effective shear zone width.

Domain 3 is marked by the transition from quartz recrystallizing by subgrain rotation to high-T GBM. Quartz grains in this domain are significantly coarser grained (Table 1) and preserve an irregular, amoeboid texture characteristic of high-T GBM (Figures 4(d), 7(f), and 7(g)). These amoeboid grains include minor subgrain development and commonly show undulose extinction, indicating limited buildup of dislocations in the crystal lattice. At higher structural levels within Domain 3, straight grain boundaries and stable 120° grain triple junctions [59] mark recovery and likely grain growth driven by minimization of surface energy (γ-GBM of [60]). Based on the relative structural position of dominant recrystallization mechanisms, we interpret γ-GBM to have overprinted dynamic high-T GBM microstructures (i.e., grains were annealed), although it is plausible that the two modes coexisted.

Feldspar and mica are also coarser grained than in Domains 1 and 2; single mica grains >500 μm are common, as are feldspar porphyroclasts >1 mm in diameter. Feldspars occur in two distinct populations: plagioclase (generally sodic) occurs primarily as porphyroclasts in variable proportions whereas orthoclase occurs as both porphyroclasts and as subordinate grains incorporated into the matrix. Within Domains 3 and 4, some plagioclase feldspar grains are irregularly zoned with more sodic cores (~An_40-50) and irregular calcic (~An_60-70) overgrowths; these observations are consistent with other chemical analyses of feldspars from within the Moine Nappe [35, 38]. Feldspars show variable degrees of dynamic recrystallization, where recrystallized, new grains tend to be a calcic plagioclase, regardless of the composition of the grain they originated from. Within the basement gneiss slice that crops out immediately above the BHT, recrystallized feldspars have distinctive lozenge shapes, with grain long axes parallel to the foliation and grain sizes on the order of 100 to >500 μm.

Domain 4 occupies the structurally highest levels of the Moine Nappe and continues up into the Ben Hope Nappe. Domain 4 is similar to Domain 3 although quartz in Domain 4 has undergone pervasive annealing and probable grain growth, resulting in a larger grain size (Figures 4(c) and 7(h), Table 1). Quartz preserves a variable to strong CPO although the dynamically recrystallized microstructure has been overprinted by annealing of grains (indicative of recovery and γ-GBM). The most abundant phases besides quartz are feldspar and mica, which have similar compositions, grain sizes, and textures to those in Domain 3.

The Stipp and Tullis [62] and Cross et al. [61] piezometers are calibrated to grain sizes up to ~50 μm, where microstructures indicate dynamic recrystallization primarily by BLG and SGR. The calibration has been extrapolated to grain sizes up to >120 μm, corresponding to deformation primarily by SGR and high-T GBM, but yields only minimum differential stress estimates for these large grain sizes [60, 68]. We also calculate stresses for microstructures which clearly show evidence of recovery (e.g., γ-GBM in Domain 4) and accept that these calculations are only crude estimates of minimum differential stress.

We measured dynamically recrystallized grain sizes of quartz in 19 samples using both electron backscatter diffraction (EBSD) and Fiji optical analysis software, an extended toolbox for ImageJ (http://www.fiji.sc; [69]). Thin sections were cut and analyzed parallel to the lineation and perpendicular to the foliation. EBSD acquisition details and data cleaning routine are detailed in Appendix A.

For coarser-grained samples (>100 μm), optical images were acquired and grain boundaries defined in Fiji to quantify grain size distribution. Multiple photomicrographs were taken with differently angled polarizers to more accurately distinguish grain boundaries.

All grain size measurements are expressed as the diameter of a circle with an equivalent area to calculated grain polygons. A minimum of 350 grains were used to calculate the mean grain size within a sample. We calculated the grain size distribution frequency peak, as outlined by Lopez-Sanchez and Llana-Fúnez [70]. However, we found that our EBSD-based grain size analyses contained an abundance of small grains which result in a positively skewed distribution and a frequency peak significantly lower than the arithmetic mean or root mean square (RMS) grain size. For this reason and for consistency with the piezometric calibrations, final recrystallized grain size for each sample is taken as the RMS of all measurements after cleaning, with errors reported as 1σ (Table 1).

We used EBSD data from all 19 samples to determine crystallographic preferred orientations (CPO) of recrystallized quartz as evidence for recrystallization by dislocation creep, qualitative conditions of deformation, and active slip systems throughout the structural section of the Strabeg transect. $c$-, $a$-, and $m$-axis pole figures (PF) and X-direction inverse pole figures (IPF-X) were constructed from a calculated orientation distribution function using MTEX (http://www.mtex-toolbox.github.io) based on a one point-per-grain statistical calculation. PFs are oriented in the X-Z plane with the foliation oriented parallel to the X-direction, thus interpreted to lie close to the shear plane in the direction of shear. PFs were plotted as lower hemisphere, equal-area projections, whereas IPFs were plotted as X-direction upper hemisphere projections. Individual scales were used for each PF and IPF to bring out CPO intricacies in low-intensity fabrics that would otherwise be lost if a common scale were applied. CPO intensity was measured using the ODF-based m-index of Skemer et al. [71] and j-index [72]. Each PF and IPF was constructed using only dynamically recrystallized grains; host grains are identified following the grain orientation spread discrimination technique and are excluded [61]. The quartz crystallographic data we present are derived from both monophase quartz and polyphase regions and were constructed from a number of grains greater than or equal to # grains (EBSD) from Table 1.

At all structural levels within the Moine Nappe, quartz grains exhibit evidence of grain boundary pinning by both phyllosilicates and fine-grained feldspars. These polyphase aggregates systematically contain demonstrably finer grain sizes compared to areas of pure quartz, (Figures 7(d) and 7(e)) consistent with observations elsewhere (e.g., [73–75] and references therein; [76]); for this reason, we attempt to record grain size measurements from areas free of nonquartz phases (e.g., recrystallized quartz veins, quartz-rich lenses). In some cases, chiefly in the coarser-grained rocks at higher structural levels in Domains 3 and 4, it is difficult to obtain a statistically significant number of grains from monophase regions. In this case, we include some grains that have pinned boundaries and accept that errors associated with these measurements may be significantly higher than those in monophase regions.

At lower structural levels within Domain 2a where total strain is likely higher, phyllosilicate-rich or fine-grained polyphase interconnected layers alternate with monophase quartz or feldspar dominated layers, defining the primary foliation. Polyphase layers show evidence of higher shear strains compared to the pure quartz or quartz + feldspar layers (Figure 8(b)). Similarly, the presence of buckle-style folds in quartz layers at higher structural levels within Domain 2b indicates a viscosity contrast between these layers and the surrounding polyphase matrix (Figure 7(e)).

Based on these two primary lines of evidence, we interpret the quartz and quartz + feldspar layers as more competent, with more strain accommodated in the weaker polyphase layers. These weak polyphase layers may deform by a combination of dislocation creep and grain size sensitive creep, including grain boundary sliding (GBS), especially parallel to the long axis of phyllosilicates. Because the areas where we collected grain size measurements are possibly stronger, our stress determinations may overestimate the average stress in the system.

We determined deformation temperatures (T) and pressures (P) for 11 samples using the Ti-in-quartz (TitaniQ) thermobarometer of Wark and Watson [77], as calibrated by Thomas et al. [78, 79], in conjunction with thermodynamic modeling of Si-in-phengite [80] and TiO₂ activity ([81]; Figure 9). Temperature and pressure calculations can be attributed to a specific quartz- and mica-dominated microstructure, thereby determining the conditions at which the microstructure was formed. Both quartz and white mica record progressive recrystallization and reduction in grain size along the transect, suggesting that both phases recrystallized and reequilibrated to ambient conditions together. As such, the temperatures and pressures we calculate represent the last stage of deformation before the rock left the active shear zone (Figure 6), and may not be the same those calculated using traditional thermobarometry on the metamorphic assemblage as a whole. For 9 samples that record crossed-girdle $c$-axis topologies, we use the quartz $c$-axis opening angle thermometer calibration of Faleiros et al. (2016) for additional temperature estimates, chiefly to compare to our TitaniQ-based temperature estimates as well as previous opening angle temperatures from the area ([83–86] and references therein).

For each sample, we calculate T, P, and a(TiO₂) using a combination of Ti-in-quartz (Figure 10(a)), Si-in-phengite (Figure 10(b)), and titania activity (a(TiO₂)) pseudosection modeling (Figures 9(a) and 9(b)); details of these analytical techniques and methods are outlined in Appendix B. a(TiO₂) is contoured in P-T space based on the bulk composition of each sample (see [81]). Due to the temperature dependence of titanium activity in our samples, modeled a(TiO₂) generally increases with temperature (Figures 9(b) and 10(c)). From the solubility equation for Ti in quartz (e.g., Equation (B.1); [79]), the position of the equilibrium line in P-T space for a given Ti concentration in quartz ([Ti]) will shift to lower temperatures with increasing a(TiO₂) in the system. By plotting the positions of the equilibrium line for the measured [Ti] as a function of a(TiO₂), and contouring for a(TiO₂) from 0 to 1, we can determine the points in P-T space where the equilibrium lines at different a(TiO₂) intersect the corresponding contours. We then graphically determine where this array intersects the modeled Si-in-phengite isopleth (corresponding to the measured Si PFU). Errors are reported as 1σ of the measured analytical values for Si-in-phengite PFU and Ti-in-quartz concentration (Figure 9(b)).

Grain sizes along the Strabeg transect range from $~9±4μm$ at the structurally lowest levels, where recrystallization by BLG dominates, to $~127±55μm$ just below the transition from Domain 3 to 4 (sample LS-105), to $~175±87μm$ at the structurally highest levels, where recrystallization occurs chiefly by high-T GBM (Table 1, Figure 11(a)). These values correspond to differential stresses ranging from ~119 to 15 MPa, respectively, using the 1 μm step size calibration of Cross et al. [61] (Figure 11(b)). For structurally higher samples in which grain sizes were measured both optically and by EBSD, we observe a marked decrease in average recrystallized grain size measured by EBSD compared to optical measurements from the same sample. White [75] also reported different grain sizes determined by electron and optical microscopy and suggested that such measurements not be interposed. Optically measured recrystallized grain sizes range from $~68±25μm$ at a structural distance of 866 m, to $~143±73μm$ just below the transition from Domain 3 to 4 (sample LS-105), to $~255±121μm$ at the structurally highest levels (Figure 11(a)). These values correspond to differential stresses ranging from ~23 to 8 MPa using the piezometer of Stipp and Tullis [62] (Figure 11(b)).

Our EBSD-based grain size measurements from the lowest structural levels are similar to previous optical determinations of recrystallized grain size from the region. White [75, 87] measured grain sizes of 14.6 μm from the footwall of the MT at Eriboll, Ord and Christie [88] recorded grain sizes as small as 12.7 μm in the footwall of the MT in the Assynt Culmination, and Weathers et al. [89] determined grain sizes of 10-20 μm from the Assynt Culmination and Eriboll (Figures 1 and 2). These measurements were taken chiefly from within quartzites of the Eriboll Fm. in the footwall of the MT. Additionally, Francsis [83] measured recrystallized grain sizes optically 5 km to the south of the Strabeg transect; these measurements ranged from $26.2±10μm$ at 80 m below the MT, up to $203.5±96.1μm$ at 2294 m above the MT, in the hanging wall of the BHT. These grain size measurements are larger than those we measured by EBSD, but are close to our optical measurements.

All samples show some degree of a CPO in recrystallized quartz, consistent with dislocation creep as the primary deformation mechanism. CPO intensities, measured by the $m$-index [71] and j-index [72], show significant variability (Figure 12). All CPOs are inclined in an orientation consistent with top-WNW shearing. CPOs along the Strabeg transect are dominated by Y-maxima, crossed-, and single girdles, consistent with the activity of basal <a>, rhomb <a>, and prism <a> slip systems; we see no evidence for prism [c] slip at these structural levels. Slip systems are likely at least partially temperature-dependent with basal <a>, rhomb <a>, and prism <a> representing progressively higher temperatures. They may also depend on finite strain, with prism <a> favored by higher strains (Figure 12, [14, 90] and references therein; [91]). Temperatures, based on active slip systems, generally increase structurally upward from the LRT to the BHT although the effect of temperature in isolation on CPO topology remains unclear [14].

A general trend emerges in the shape of the CPO through the structural section (Figure 12). At the lowest structural levels around the LRT (Domain 1), CPOs are defined by weak $c$-axis crossed-girdles suggesting activity of basal <a>, rhomb <a>, and prism <a> slip systems. At structural levels immediately above and below the MT (Domain 2a, b), $c$-axes show somewhat stronger crossed-girdles and single girdles. Single girdles are generally stronger and may show weak complementary crossed-girdle legs. CPOs in the lower structural levels of Domain 3 are characterized by strong maxima normal to the foliation (Y-maxima), which are indicative of prism <a> slip. Structurally higher (i.e., higher structural levels of Domain 3 and Domain 4), textures transition to single girdles with elongate Y-maxima parallel to the long axis of the girdle, indicative of prism <a> and rhomb <a> slip. Secondary maxima generally expressed along the pole figure periphery indicates smaller contributions from the basal <a> slip system. As in Domain 2, generally weak maxima exist in some samples along the pole figure periphery that represent crossed-girdle legs. Strong Y-maxima in the immediate hanging wall of the BHT and at the structurally lower levels of Domain 3 may result from large shear strains [90, 91] or enhanced hydrolytic weakening (e.g. [84, 92]).

As illustrated in Figure 12, many of the crystallographic textures plot as single girdles with absent or very weak crossed-girdle legs that may not be expressed in the contouring. These patterns differ from previously published crystallographic data from northwest Scotland (e.g., [45, 84–86] and references therein) which commonly show crossed-girdle $c$-axis patterns, used to determine temperatures based on crossed-girdle opening angle. The paucity of crossed-girdle patterns in our EBSD-based $c$-axis data can likely be attributed to the techniques used for fabric measurement and contouring (R. Law, personal communication); whereas EBSD-based crystallographic data are produced by one point-per-grain analyses including all grains within a large area, measurements by universal stage, utilized in many previous studies, require individual measurements of grains that occupy the entire section thickness, and may introduce subjectivity.

Like other studies that have employed the TitaniQ thermobarometer in deformed terranes (e.g., [93–96]), our analyses reveal low Ti concentrations in dynamically recrystallized quartz. Despite low Ti concentrations, the data follow a generally well-defined trend, increasing structurally up from 0.58 ppm near the MT, to 2.57 ppm in the hanging wall of the BHT (Figure 10(a)), whereas Si-in-phengite shows a broadly increasing trend over the same structural section from 3.23–3.29 PFU (Figure 10(b)). Utilizing the modeling approach detailed in Appendix B, Ti concentrations correspond to temperatures and pressures ranging from approximately 350–450°C and 270–560 MPa (2.7–5.6 kbar; Figures 4(b) and 10(c), Table 1). Quartz $c$-axis opening angle temperatures broadly increase structurally up, albeit with significantly more scatter, from 334–555° C (Figure 12, Table 1). Both TitaniQ and $c$-axis opening angle temperatures agree with trends identified by previous workers, who have found structurally upward-increasing field gradients within individual nappes [45, 97]. Our calculated TitaniQ temperatures are also generally lower than temperatures estimated by quartz recrystallization mechanism; recrystallization by BLG, SGR, and GBM are suggested to occur within temperature ranges of ~300-400°C, 400-525°C, and >525°C, respectively [55, 98]. Along the Strabeg transect, rocks deforming by BLG have calculated temperatures of ~350°C, recrystallization by SGR corresponds to temperatures of ~350–380° C, and recrystallization by high-T GBM range from ~370–450°C. It is worth noting, however, that the recrystallization mechanism temperature estimates of Stipp et al. [55] are based on syn-kinematic mineral assemblages, representing “close to peak” temperatures in a system that continued deforming while cooling, and therefore may be capturing a different part of the P-T-t path than TitaniQ estimates.

TitaniQ-derived deformation temperatures calculated here tend to be significantly lower than published garnet-biotite Fe-Mg exchange (GARB) metamorphic temperatures, quartz $c$-axis opening angle deformation temperatures, and multisystem P-T analyses [43–45, 48, 97]. We offer the following possible explanations for this discrepancy.

• (1)

  Temperatures derived from TitaniQ are strongly dependent on a(TiO₂), and small changes in a(TiO₂), especially at low values (i.e., a(TiO₂) <0.3), can result in significant differences in P-T conditions. Although the method we follow yields a simultaneous, unique solution for P, T, and a(TiO₂), it is possible that analytical errors arising from XRF, SIMS, or EPMA may result in a(TiO₂) values that are systematically too high, and hence temperatures that are too low. For our TitaniQ-based temperatures to come within the error of multisystem and GARB values, a value of a(TiO₂) of <0.1 is necessary. However, based on the temperature dependence of a(TiO₂), the presence of syn-deformational titanite, and the homogeneous distribution of Ti at both the grain- and sample-scales, we consider a(TiO₂) values this low to be unlikely.

• (2)

  A more fundamental flaw in our understanding of the mechanics of Ti mobility and substitution into quartz affected by crystal-plastic deformation may also be possible. Ashley et al. [99] suggest that low [Ti] in dynamically recrystallized quartz could be a result of local reequilibration from subgrain boundaries and dislocation arrays migrating through quartz grains. The reequilibration in this case is suggested to be buffered or regulated by the composition of the intergranular medium which is typically Ti-undersaturated with respect to the overall assemblage and may not represent the actual a(TiO₂). More recent experimental work however indicates that Ti concentrations are reequilibrated during recrystallization, reflecting bulk a(TiO₂) [100, 101].

• (3)

  Another possible explanation for the discrepancy in temperatures stems from what portion of the P-T-t history a thermometer or thermobarometer captures. Metamorphic temperatures calculated from metamorphic thermobarometers or multisystem equilibria (GARB included) likely record “peak” or “near-peak” prograde metamorphic temperatures. Within rocks of the Moine Supergroup, garnets are present as porphyroblasts of variable size, morphology, and chemistry and are commonly chemically zoned with relative proportions of Mn higher in the cores and Fe and Ca higher in the rims (see [33, 43, 45, 48, 97]). Optically, zoning is commonly defined by linear or spiral inclusion trails either in the cores or rims that are discordant to the primary foliation. This physical and chemical zoning indicates that garnets record growth over multiple episodes, likely at variable P-T-X. At structurally lower levels along the Strabeg transect garnets tend to be anhedral or skeletal, suggesting that they could be out of equilibrium. In contrast, if [Ti] reequilibrates with dynamic recrystallization, TitaniQ temperatures should record the conditions when plastic deformation and dynamic recrystallization ceased (i.e., when the rock left the actively-deforming shear zone).

• (4)

  Based on our interpretation of shear zone evolution, deformation temperatures recorded in the quartz microstructures should be well below peak conditions. $c$-axis opening angle temperatures are inconsistent with this interpretation whereas our TitaniQ temperatures, albeit low, are consistent. Although the similarity between opening angle temperatures and petrologic thermobarometers (e.g., [84] their Figure 24) is compelling, considerable uncertainty remains regarding the effects of strain geometry, recrystallization mechanism, strain rate, and water content on $c$-axis opening angles (see [84] for further discussion). Furthermore, existing $c$-axis opening angle calibrations are based almost exclusively on temperature estimates derived from peak or near-peak metamorphic assemblages ([84] and references therein; Faleiros et al., 2016); although these may represent conditions during deformation, they are not necessarily representative of the last deformation conditions the rocks last experienced before being “locked in”. As such, we are hesitant to assign temperatures calculated from $c$-axis opening angles as being representative of any specific (or consistent) time in the deformation history (e.g., [84]).

• (5)

  Lastly, based on microstructural observations, the trace of the MT (i.e., base of the Moine Nappe) exposes different structural levels along strike ranging from very shallow, mylonite-on-foreland (e.g., Knockan Crag – Figure 1) to high-temperature GBM microstructures within mylonitic Moine schist lying on mylonitic Oystershell rock (a phyllonite thought to be derived from Lewisian basement gneisses that commonly crops out within the MTZ) along the north coast [34, 38, 46, 102]. So, although the Moine Nappe is a coherent tectonic unit, pressure and temperature conditions vary both along and across strike, as do conditions at the base of the nappe. This is plausibly a result of differential syn- and postdeformation erosion and exhumation. In this case, it is likely that previously published estimates are measuring temperatures and pressures attained at different structural levels and at different times within the Moine Nappe

For the structural thickness of a shear zone to remain constant through time, the shear zone must not undergo significant thinning or volume loss. We argue here that the shear zone deformed primarily by plane strain and simple shear (noncoaxial) and did not undergo significant volume loss. In principle, this can be tested by a variety of techniques that allow the kinematic vorticity number (⁠$Wk$⁠) to be estimated. Thigpen et al. [85], for example, use rigid grain analysis to calculate $Wk$ values between ~0.6–0.7 for rocks within the Moine Nappe on the Eriboll Peninsula, corresponding to ~60–50% pure shear (coaxial). However, in numerous deformed terranes, calculation of $Wk$ by rigid grain analysis appears to underestimate the component of noncoaxial strain relative to other methods [103–105] and uncertainties associated with these methods may be significant [106]. Law (2010) reports values of $Wk$ from within 100 m structurally above and below the MT at the Stack of Glencoul, using rigid grain analysis for the Moine mylonites, and two additional techniques involving the geometry of the CPO for the mylonitic Cambrian quartzites below the Moine thrust. The values range from 0.75–0.65 (45–55% coaxial shear) to 0.99–0.90 (10–30% coaxial shear), and generally increase (i.e., a larger component of simple shear) with proximity to the MT. A significant component of coaxial shortening in such high-strain rocks presents serious compatibility issues, however. X-Z finite strain ratios in from detrital and recrystallized grains in the Eriboll Fm. commonly range between 10 and 19 (Law, 2010). Assuming no volume loss, and $Wk≈0.77$⁠, an X-Z ratio of 19 requires 65% shortening normal to the shear plane, and hence a stretch of 2.1 in the direction of shear (Law, 2010).

The situation in the Moine mylonites is much more extreme. Strain data are lacking, but if we take a very conservative estimate of 1 km displacement across the basal 100 m of mylonite, giving a component of simple shear strain γ of 10, take $Wk≈0.75$ as estimated by Law (2010), and assume plane strain and no volume loss, then shortening normal to the shear plane is 88%, and the pure shear related stretch in the direction of shear is 8.2. Unless the entire thrust pile above the mylonite zone experiences the same amount of stretch, this requires extrusion of many kilometers of mylonite at the thrust front. This type of extrusion model has been invoked to explain the emplacement of the Greater Himalayan Sequence [107–109], but it requires a reversal of the shear sense across the mylonite zone, and we see no evidence for this. Furthermore, there is no evidence for thinning of the Caledonian orogen as a whole by a factor of 10 or more, which would be required to avoid extrusion of the shear zone. Lastly, although kinematic vorticity analyses commonly indicate general shear, quartz $c$-axis fabrics (e.g., this study, [85, 86, 110]) show Type 1 crossed-girdles consistent with plane strain and primarily simple shear [111, 112]. As such, we assume plane strain and simple shear dominated, acknowledging that there may be subsidiary components of general shear.

The structural continuity, distribution of microstructures, and a metamorphic field gradient that increases structurally up collectively suggest that our samples from the Strabeg transect record the progressive evolution of a single Scandian-aged shear zone. The structurally lowest rocks cropping out immediately below the lithologic MT record the narrowest, most localized, highest stress, and lowest temperature sections of the shear zone. Conversely, structurally higher rocks from near the top of the Moine Nappe, which record higher deformation temperatures and lower stresses, record the deeper, wider portions of the shear zone (Figure 6). Rocks in the lower Ben Hope Nappe record conditions of deformation of the Scandian shear zone at greater depths than rocks from the top of the Moine Nappe. Therefore, structural distances calculated from samples within the Ben Hope Nappe represent cumulative thicknesses of the active parts of the Moine and Ben Hope Nappes. We include these rocks in this analysis acknowledging that the calculation of shear zone width at these highest structural levels yields only an approximate value, as we cannot constrain the displacement on the BHT.

As strain localizes up dip in an idealized and simplified reverse-sense fault system, hanging-wall rocks near the upper margin of the shear zone will leave the shear zone as they move up dip, and the shear zone narrows. Their microstructures will therefore be “locked in”, recording the conditions of the shear zone at the depth they left it (Figure 6). If deformation is dominated by simple shear in which strain is homogeneously distributed (i.e., the shear zone at any given depth deforms at the same stress, temperature, and strain rate), the width of the shear zone for a given microstructure and conditions of deformation will be recorded as the structural distance from the projected fault plane to the sample of interest. With these assumptions, we make a space-for-time substitution to reconstruct shear zone geometry, internal structure, and rheology.

We present a model for shear zone geometry as a function of depth based on an assumed fault dip (15°), plus calculated shear zone thickness and depth. Depths are estimated from pressure calculations (Table 1, Figure 10(c)), based on a density of 2.75 g cm^-3 for rocks of the Moine Supergroup (Rollin, 1994). We neglect the effect of tectonic overpressure, as our stress measurements indicate that this did not exceed 100 MPa (1 kbar), which is within the uncertainties on our pressure estimates. Our reconstruction, based on thickness-depth data from 10 samples along the Strabeg transect, is illustrated in Figure 13. Because of the large errors associated with calculated pressures, not all samples plot in a sequence of increasing structural thickness with increasing depth, as would be expected in an ideal scenario. Additional uncertainty in total shear zone thickness arises from the thickness of the footwall material incorporated into the shear zone (teal-colored zone in Figure 13). The lower limit of shearing is poorly constrained, but this uncertainty is likely to be only a fraction of the total shear zone thickness. We base this on the relatively minor (i.e., 10 m-scale) thickness of the Lewisian basement that is incorporated along the MT and BHT (Figures 1, 2(a), and 2(b)). Similar to what is observed in the hanging wall Moine schist, the amount of basement incorporated into the shear zone likely increases with increasing depth.

Even with these uncertainties, this model illustrates a general trend of a wider shear zone at depth. We estimate a thickness of ~2.5 km at ~20 km depth, and extrapolating downwards along the projected fault dip results in a shear zone structural thickness >5 km at ~25 km depth. We base our synoptic shear zone (Figure 6) on this modeled profile.

Experimental rock studies have been paramount in our understanding of rock strength and rheology in regions of both brittle and ductile deformation. However, the difference between experimental deformation conditions and geologic conditions for temperature, stress, and strain rate is large, and we must therefore rely on scaling relationships to apply experimental results to geological conditions. These relationships require validation from field-based studies on naturally-deformed rocks to determine how well experimental constraints approximate deformation at natural conditions. In the following section, we use calculated deformation temperatures, pressures, and differential stresses to compare our natural data to currently published flow laws for dislocation creep in quartz. We then use strain rates predicted from published flow laws and compare them to the geometrical strain rates calculated independently from field data.

Flow laws are plotted in T-σ space with strain rate contouring from 10^-10 to 10^-16 s^-1 (Figure 14). Data from the Strabeg transect generally lie between strain rates of 10^-15 and 10^-13 s^-1 for Hirth et al. [113] and between 10^-16 and 10^-14 s^-1 for Tokle et al. [114]. Data show a trend of increasing strain rates at lower temperatures and higher stresses, reflecting the shear zone narrowing due to localization of strain and resulting tendency towards higher strain rates.

The Tokle et al. [114] flow law is proposed specifically for prism <a> slip whereas rocks along the Strabeg transect clearly show evidence for activation of multiple slip systems (Figure 12). For quartz that shows contributions from basal <a>, rhomb <a>, and prism <a> slip systems (i.e., single girdle CPOs), Tokle et al. [114] suggest grain boundary sliding (GBS) as a deformation mechanism that connects the prism <a> limiting and basal <a> limiting dislocation creep regimes. We see no evidence for GBS in the quartz microstructure of samples with single girdle CPOs (e.g., parts of Domain 3 and 4). Dynamically recrystallized grain sizes in these samples are significantly larger than those where GBS has been proposed in other naturally-deformed quartz-rich rocks [119–122].

A value for imposed plate velocity is calculated from the total displacement (shortening) and duration of thrusting on two faults below the MT in the Assynt Culmination, located ~30 km south of the Strabeg transect (Figure 1). This calculated value is then compared to generalized tectonic velocities for the region. Shear zone thickness is taken as the structural (orthogonal) distance to the projected fault plane, plus the thickness we determine of the shear zone in the footwall. These thicknesses are based on four detailed structural and microstructural transects parallel to the direction of displacement, and we justify the thickness calculations below.

In the above sections, we have discussed, in detail, the calculation of shear zone thickness relative to the lithologic MT. Examination and quantification of rock microstructures indicates that this lithologic contact is of little rheological significance and acts primarily as a passive marker from which we measure the structural distance. We have presented evidence for moderate to high strains preserved in rocks from the footwall of the lithologic MT. Therefore, for calculations which concern the total thickness of the active shear zone, we add the distance from the lithologic MT (defined as 0 m) to LS-26 (172 m below the MT) to each thickness measurement at and above the MT (for LS-27, at 58 m below the lithologic MT, we calculate a shear zone $thickness=172–58m$⁠). LS-26 shows moderate strains compared to structurally higher rocks indicating that it was likely at the lower shear zone margin (Figure 7(c)). There is no exposure of the brittle MT along the Strabeg transect, and we estimate that the minimum thickness of the shear zone at the deformational depths recorded by the structurally lowest samples was likely on the order of 100–150 m.

Displacement rate (⁠$V$⁠) is calculated by dividing total thrust convergence (km) by duration of thrusting (Myr). Elliott and Johnson [34] estimate 20–25 km and 28 km of displacement on the Glencoul and Ben More thrusts in the Assynt Window, respectively. The timing of deformation is well-constrained based on U-Pb geochronology on a suite of alkali intrusions. Syn-kinematic (Loch Ailsh Pluton, early parts of the Loch Borralan Pluton, and the Canisp Porphyry sills) and postkinematic (Loch Borralan Pluton) intrusions constrain thrusting to have been active between $430.6±0.3Ma$ and $429.2±0.5Ma$⁠, although it could possibly have initiated earlier [36]. Taking an estimate for shortening (50 km) and time span (0.6–2.2 Myr) yields a calculated displacement rate of 23–83 mm yr^-1. As a conservative estimate, we use the minimum displacement rate of 23 mm yr^-1 but also calculate an upper limit based on 35 mm yr^-1 (~50% increase in displacement rate). These timing constraints are from the Assynt Window, which is 40–50 km south of the Strabeg transect, but the continuity of structures along strike strongly suggest contemporaneity.

The tectonic rate of orogen-normal motion for the Scandian phase of the Caledonian Orogeny in NW Scotland is constrained between 30–60 mm yr^-1 between Laurentia and Baltica based on plate tectonic reconstruction models [18]. Our local calculated rates are in general agreement with these tectonic-scale estimates.

To use this displacement rate to calculate strain rates, we make several assumptions. First, we assume that the displacement rate did not vary significantly in time or along strike. This is especially important since the displacement rate is calculated based on thrusting below the MT, requiring the displacement rate to be constant over a period of c. 5–10 Myr prior to motion along the Glencoul and Ben More thrusts. Second, we assume that strain was accommodated along a single active fault (shear) zone at any given time. It is, however, likely that minor components were partitioned along subsidiary structures, which is why we elect to use the minimum value for the estimated displacement rate. Although assumptions are made to independently calculate strain rate and we accept that these are first-order estimations, it is important to realize that order-of-magnitude strain rates will not be significantly affected by even moderate errors in the values for velocity and/or shear zone thickness. We are confident that we can constrain the displacement rate to within ± ~50% and shear zone thickness to ± ~20%, which will not drastically affect calculated strain rates.

As discussed above, microstructures at any given structural distance from the fault plane likely preserve the conditions present in the shear zone at that given thickness (and time). Therefore, we assign temperature, pressure, and stress to a specific structural thickness, which we then attribute to a calculated strain rate.

Although some estimates are within error, flow law-based strain rates (from [110, 111]) are consistently lower for given temperature, pressure, and stress conditions compared to independent field-based calculations, sometimes by >1 order-of-magnitude (Figure 15, Table 2). Note that if differential stresses were lower than those calculated (due to strain partitioning between quartz-rich and polyphase regions), or calculated displacement rates were higher (i.e., not the minimum estimate we use), the discrepancy in experimental versus natural strain rates would become even larger. These differences may indicate (a) that dislocation creep in quartz is not the primary control on rheology, (b) inaccuracies in estimates of field-based strain rate, and/or (c) selected flow laws do not accurately model rock rheology at geological conditions. We address each of these possibilities below.

• (a)

  Based on our optical and textural analyses, we conclude that quartz is the primary phase controlling viscous rheology in most parts of the shear zone. Psammitic Moine schist is a quartz-rich rock and it is probable that a rock with a higher proportion of feldspar (e.g., a granitic or granodioritic composition) will be stronger, reflecting the rheology of the feldspar, or of a polyphase mixture of feldspar and quartz [123]. We do however observe differences in rock rheology related to composition and grain size variability at the thin section-scale (Figure 8(b)). Interconnectivity of phyllosilicates can result in weaker rock strength due to easy glide between phyllosilicate basal planes [124–127]. Additionally, reduction in grain size due to pinning in a polyphase material may lead to a switch to grain size sensitive creep in quartz [128, 129].

• (b)

  As discussed above, even moderate uncertainties in shear zone width or displacement rate do not substantially affect our order-of-magnitude estimates in strain rate. Furthermore, our field-based strain rates (Figure 15) are within the range of other independent estimates of strain rates from plate boundary-scale fault systems (Sassier et al., 2009; [129, 130]).

• (c)

  Significant variability exists between published calibrations of flow laws for dislocation creep in quartz [113, 114, 131–135]. Variability in values of H, n, and r in equation (2), which may stem from differences in starting materials, experimental conditions (including deformation apparatus), poorly constrained fluid content, among other uncertainties, can result in order-of-magnitude differences in predicted strain rate. More recent work also introduces a pressure sensitivity for activation enthalpy (H), adding further complexity [133]. We contend that experimental-based constitutive laws describing rock rheology should be thoroughly tested, and if necessary altered, to better fit well-constrained field-based data, as suggested by Hirth et al. [113].

The most significant difference between predicted and field-based strain rate is in samples LS-72 and LS-74, which are derived from the zone above the Lewisian inlier in Domain 2. Both samples record temperatures and stresses lower than those from samples structurally above or below (Table 1). If this inlier is an imbricate slice and acted as a mechanically strong block, then the effective shear zone width represented by both LS-72 and LS-74 would decrease by ~25-30%, resulting in faster strain rates closer to those predicted based on sample temperature, pressure, and stress. The effect of decreasing shear zone width on structurally higher samples would be minor as the thickness of the imbricate slice is a smaller proportion of the total shear zone width. The width of this secondary zone is likely only a fraction of the wider primary shear zone, which would rectify lower strain rate estimates (for given stress, temperature, and pressure conditions) calculated from flow laws. A secondary localization is also consistent with lower calculated temperatures; we do not however observe a smaller recrystallized grain size, which is also predicted by this model (Table 1).

The strength of the upper crust is thought to be controlled primarily by frictional processes along active faults following a Mohr-Coulomb criterion for frictional sliding. Empirical and experimental estimates of the coefficient of friction (μ) generally range between 0.6–0.85 for undamaged crustal materials [136]. More realistic estimates indicate that the upper crust is much weaker than predicted by Byerlee’s Law because of low values of the effective coefficient of friction of fault gouges and clay-rich fault rocks (e.g., $μ<0.1$⁠). These low values of the effective coefficient of friction are likely temperature dependent due to the thermal stability of clay minerals, and increase closer to the BDT as clays become unstable and/or frictionally stronger phases are precipitated ([16] and references therein). Because of the transition from frictionally weak to frictionally strong materials close to the BDT, Byerlee-type friction can still be used to approximate rock strength in this region [16].

Assuming a simplified system where rock strength is controlled primarily by μ, frictional strength is depth-dependent due to the increase in normal stress (pressure) with increasing depth. At temperatures of ~300°C for quartz-rich rocks, thermally-activated creep processes become the dominant deformation mechanism, following a power-law relationship between stress and strain rate (equation (2)). Viscous creep processes are highly temperature-dependent, resulting in an exponential decrease of crustal strength with depth. In these models, the strongest region of the crust is predicted to be the transition between frictional and viscous rheology (see [16, 17]). In the following section, we construct a first-order, naturally-constrained, crustal strength profile for a contractional fault system by plotting shear stress and deformation depth data from rocks along the Strabeg transect in conjunction with frictional strength as predicted by Byerlee’s Law. Note that stresses measured by paleopiezometry (Table 1) are calibrated in uniaxial compression (⁠$σ2=σ3$⁠) and are expressed in terms of differential stress (⁠$σd=σ1–σ3$⁠). To compare these stresses to those driving brittle faulting, we convert them to shear stresses in a plane stress environment by dividing the values of differential stress by √3 as outlined in Paterson and Olgaard [137].

Rocks at and below the trace of the MT at the Strabeg transect record quartz microstructures indicative of viscous deformation, but lack evidence for extensive frictional or brittle behavior, likely because the structural level of the BDT on the MT has been eroded along much of the fault trace. Faults closer to the foreland (west), however, record deformation conditions at the BDT. Therefore, we quantify paleostress on a sample (GT-3) from the immediate hanging wall of the Glencoul Thrust (GT) at Loch Glencoul (red star in Figure 1) which places Lewisian gneiss on a thin veneer of An t-Sron siltstones and sandstones, which in turn conformably overlie Eriboll quartzite [47]. Coexisting brittle (e.g., quartz fracture and offset along brittle shear bands) and ductile (e.g., incipient dynamic recrystallization of quartz by BLG) structures clearly place this sample around the BDT when the microstructure formed (Figures 8(c) and 8(d)). Although the precise depth that sample GT-3 records is difficult to determine, we calculate an approximate depth of ~11 (+4/-2) km based on a temperature of $~280±20°C$ for the lower limit of crystal plastic behavior in quartz and a geothermal gradient of $25±5°Ckm−1$⁠. This depth estimate is consistent with previous estimates of 200–250°C and 200–300 MPa in the immediate footwall of the MTZ (i.e., to the west of GT-3; [39]) and is also consistent with depths for samples from the Strabeg transect which approach the BDT. EBSD analysis on sample GT-3 yields an average recrystallized grain size of 3.9 ± 1.5, resulting in $σd=225$ (+92/-46) MPa based on the Cross et al. [61] piezometer, which is equivalent to a shear stress in a plane stress environment of $σs=130$ (+53/-26) MPa.

Following the Andersonian model for fault classification, principal stresses driving deformation in a thrust regime are oriented such that the least compressive principal axis, $σ3$⁠, is vertical and is taken as the effective pressure (equal to the lithostatic pressure less the pore fluid pressure). $σ1$⁠, the largest compressive stress, is horizontal, acting in the direction of shortening. Assuming plane stress, $σ2=σ1+σ3/2$⁠. Because of the orientation of these stresses relative to the Earth’s surface, shear stress along a fault, an analog for crustal strength, is predicted to be larger for a reverse regime compared to corresponding shear stresses in a normal or strike-slip regime for a given depth [138]. In a dry system where $σ3$ is equal to the lithostatic pressure, rock strength in a contractional regime without preexisting faults (i.e., $μ=0.85$⁠) is predicted to exceed 350 MPa. As discussed above, this is an unreasonable estimate of bulk rock strength as thrust faults have been shown to have very low (e.g., <0.1) effective coefficients of friction ([139]; see also discussion in [16]). This estimate was also made assuming a dry system with no effect of pore fluid pressure. The effect of pore fluid pressure in rocks that follow a Byerlee-type frictional rheology is to lower the effective normal stress; indeed, high pore fluid pressures have long been invoked as a mechanism to explain apparently low values of effective friction along thrust faults. The presence of syn-deformational fluids along the GT and MT is evidenced by precipitation of quartz and pervasive chloritization and sericitization of nominally anhydrous minerals.

The effect of pore fluids in ductile rocks is more difficult to predict as thermally-activated processes (e.g., dislocation creep and dynamic recrystallization) likely result in a change in the area of frictional contact and therefore the effect of pore fluid pressure (P_f; see discussion in [67]). We therefore modify our estimates for Byerlee-type frictional behavior through the introduction of the term α (⁠$1≥α≥0$ as a function of asperity yield stress, modifying the P_f) following Hirth and Beeler [67]. We use the same values to calculate an asperity yield stress but use a strain rate of 10^-13 s^-1 as calculated from Hirth et al. [113] for the structurally lowest rocks in the Strabeg transect.

Based on paleopiezometric-derived estimates of shear stress on rocks deformed within the BDT, maximum shear stress is unlikely to have exceeded ~180 MPa (maximum shear stress within error for GT-3 is 183 MPa). For modeled shear stresses of Byerlee-controlled rocks to be consistent with values that we observe, $μ=0.50$⁠, and a ratio of pore fluid pressure to lithostatic pressure, $λ=0.36$⁠, are required (Figure 16).

This study presents a reconstruction of the geometry, internal structure, and rheology of a plate boundary-scale shear zone from ~10–20 km depth. We draw the following conclusions:

• (1)

  Temperatures and pressures of deformation are consistent with a tectonically-inverted field gradient (in which isograds dip to the east) that preserves a record of shear zone evolution accompanied by exhumation.

• (2)

  Quartz appears to be the primary phase controlling shear zone rheology, deforming chiefly by dislocation creep. Quartz recrystallization mechanisms within the active shear zone evolve through time from high-T GBM to BLG as temperature decreases and strain localizes. The change in recrystallization mechanism corresponds roughly with a change in active slip systems from prism <a> dominated at higher temperatures and lower stresses to a combination of prism <a>, rhomb <a>, and basal <a> at lower temperatures and higher stresses.

• (3)

  Shear zone geometry in a contractional, plate boundary-scale fault system has been independently constrained and shown to broaden with depth, to a thickness of ~2.7 km at ~20 km depth.

• (4)

  Independent, field-based estimates of strain rate decrease from approximately 10^-11–10^-12 from the BDT to 20 km depth, consistent with a shear zone model that localizes strain at lower temperatures and pressures, producing higher strain rates in narrower portions of the shear zone.

• (5)

  Rock microstructures indicate that dislocation creep in quartz is the primary deformation mechanism controlling rock rheology across the Strabeg transect and likely throughout most of the Scandian shear zone. However, independent, field-based estimates of strain rate are generally higher than strain rates predicted by existing constitutive relationships for dislocation creep in quartz, often by greater than one order-of-magnitude.

• (6)

  Calculated depths and stresses are plotted in conjunction with frictional rock strength to construct a naturally-constrained stress profile through the middle- to lower-crust in a thrust environment. We estimate that shear stress decreases from 130 to <10 MPa from the BDT to 20 km depth. Shear stress measurements compare well with those determined from a normal-sense shear zone in the Whipple Mountains, California, by Behr and Platt [15]. Note that we modified the stress estimates from Behr and Platt [15] so as to not reflect the proposed Holyoke and Kronenberg [66] stress recalibration (discussed above) that the authors originally used. Our estimated peak shear stress of 130 MPa at the BDT is somewhat higher than that in the Whipple Mountains (107 MPa shear stress), as is expected given the difference in tectonic regime. The profile shows that for values of $μ<0.50$⁠, ductile rocks from immediately below the BDT were deformed at higher stresses than rocks deforming by frictional mechanisms at higher structural levels. This suggests that the strength of the crust is controlled by this shallow ductile region of the middle crust, as suggested by Behr and Platt [16].

• (7)

  We present a method for application of TitaniQ to psammitic, semipelitic, or otherwise low a(TiO₂) mylonites by simultaneous solution of deformation temperature, pressure, and a(TiO₂). This method requires thermodynamic modeling of bulk rock composition, measurements of Ti-in-quartz, and an independent thermobarometer.

EBSD analyses were performed on a JEOL-7001F scanning electron microscope equipped with an EDAX Hikari EBSD detector at the University of Southern California Center for Nanoimaging. We used an accelerating voltage of 20 kV, probe current of 14 nA, working distance of 15 mm, and step sizes from 500 nm to 8 μm, where step size is at most ¼ the diameter of the smallest grain. Map coverage was variable, based on recrystallized grain size and step size; in coarse-grained samples, multiple maps were stitched together for a total map size of $~3×15mm$⁠, whereas in finer-grained samples, maps were as small as $~2×2mm$⁠. Map data were collected, cleaned, and plotted using OIM 7 software package by EDAX. Cleaning in OIM 7 consisted of a grain confidence index standardization, followed by a neighbor point confidence index correlation and a neighbor point orientation correlation. No grain dilation procedure was applied. Grains were defined using the OIM grain reconstruction routine; high-angle grain boundaries were defined by a minimum misorientation of 10° between neighboring grains and subgrain boundaries allowed to complete down to a misorientation of 2° based on EBSD angular precisions and community conventions [137–139]. We followed the routine of Cross et al. [58] to discriminate recrystallized and relict grains based on the grain orientation spread (GOS) trade-off curve threshold within each sample. Although the GOS has been suggested to not be an accurate metric to discriminate recrystallized and relict grains since newly recrystallized grains continually accumulate intragranular strain [140], we applied this filter for consistency with the calibration of Cross et al. [58].

In all analyses, quartz was the only phase indexed. In cases where other phases were present in the map, and to eliminate misindexing due to damage to the sample surface, we removed indexed points with an average confidence index <0.4. Because of our chosen step size of ¼ the diameter of the smallest grain, grains with <16 pixels were also removed as they are likely artefacts of the cleanup routine. EBSD analyses regularly revealed mottled or spotted patterns of pixels within grains; these patterns have misorientations of $60±5°$ around <0001>, likely recording Dauphiné twin boundaries as well as regions of systematic misindexing. These regions were removed in postprocessing to maintain consistency with the paleopiezometer calibration [58].

The TitaniQ thermobarometer utilizes the temperature and pressure dependence of the ⁴⁸Ti for ²⁸Si substitution in quartz to calculate P–T conditions. Higher temperatures or lower pressures allow for increased Ti concentration within the quartz crystal structure. Kohn and Northup [141] and Spear and Wark [142] demonstrated that Ti content in quartz reequilibrates as a result of dynamic recrystallization, suggesting that TitaniQ will actually yield deformation temperatures. Later work, however, suggested that recrystallization by GBM is necessary to facilitate Ti reequilibration, and that at temperatures below ~500° C, where recrystallization by BLG and SGR likely dominate, Ti concentrations are not reset [91]. On the basis of lower Ti concentrations in recrystallized grains compared to host porphyroclasts and low overall Ti concentrations, Ashley et al. [96] claimed that TitaniQ cannot be applied to determine the deformation temperature of quartz recrystallized by SGR. In contrast, however, more recent experimental work has demonstrated that dynamic recrystallization in quartz enhances the kinetics of Ti equilibration, regardless of recrystallization mechanism, confirming that TitaniQ does indeed record deformation temperatures [98].

The degree to which Ti substitutes for Si is also dependent on the activity of TiO₂ in the system; at titanium activity (a(TiO₂)) less than unity, the availability of Ti to substitute decreases, thus reducing the concentration of Ti in quartz. Ghent and Stout [143] estimate the following a(TiO₂) for assemblages in which Ti-bearing phases are in equilibrium: rutile $TiO2=1$⁠; ilmenite $FeTiO3=0.8$⁠; titanite $CaTiSiO5=0.6$⁠. However, using these values for a(TiO₂), numerous workers have reported anomalously low temperatures (e.g., [87–90]).

To address this issue, we implement a(TiO₂) modeling across P-T space combined with contouring of TitaniQ $XTiO2quartz$and Si-in-phengite isopleths to accurately constrain activities for each sample (Figure 9(b); see also [78] for a complete discussion of the modeling approach). Using this approach, we calculate values of a(TiO₂) for Moine psammites in the range 0.20-0.35 (Table 2), considerably lower than the broad estimates suggested by Ghent and Stout [143]. They are however in agreement with thermodynamic modeling by Ashley and Law [78] who calculate a(TiO₂) for greywacke (i.e., psammitic) compositions of <0.3 and by Kidder et al. [92] who suggest a value of $aTiO2=0.1$ for mylonites from the Alpine Fault. These lower activities result in higher temperatures for a given Ti concentration.

Due to the relatively low Ti concentrations, samples are best measured on a Secondary Ion Mass Spectrometer (SIMS). Prior to analysis, samples were examined by cathodoluminescence (CL) in the panchromatic and blue-UV wavelengths to look for Ti heterogeneity and/or zoning in quartz (e.g., Wark and Spear, 2005; [142, 145]). Within the Moine Schist, CL revealed no detectable heterogeneity, suggesting homogeneous Ti distribution in quartz. Within the Lewisian gneiss, recrystallized quartz is generally darker than the larger gneissic grains, indicating that Ti concentration in recrystallized grains is lower. This is consistent with observations that dynamic recrystallization resets Ti concentration in quartz [93–95].

Analyses were conducted on a Cameca 6f SIMS at Arizona State University (ASU). A ¹⁶O^- beam was accelerated to -12.5 kV with a primary beam current of ~15 nA and a stationary beam with a diameter of 10–30 μm. For standards, three silica glasses with concentrations of 0, 100, and 500 μg/g (ppm) were synthesized and characterized at the University of Edinburgh (https://www.ed.ac.uk/geosciences/facilities/ionprobe/standard-materials-available/tiquartzstandards). Study of these samples using the ASU Cameca 6f SIMS at the conditions described above provided a calibration curve to correlate Ti⁺/Si⁺ ion ratios with TiO₂ concentrations (Figure A1). For unknown samples, analyses started with a 300-second sputter to remove the conductive Au coating, ensuring the ion beam was contacting the rock sample. This was followed by measurements of ²⁷Al, ⁴⁰Ca, ³⁰Si, ⁴⁸Ti, and ⁴⁹Ti over 30 cycles per site. Depending on sample complexity, between 7–10 sites were chosen per sample; site locations were based on CL imaging, proximity to Ti-bearing phases, and quartz microstructure. In each measurement, we monitored two titanium isotopes ⁴⁸Ti and ⁴⁹Ti because of the unresolvable interference of ⁴⁸Ca (a minor isotope) on ⁴⁸Ti. If a Ca-rich area is unintentionally analyzed (e.g., through accidental overlap of the analyzed area on a plagioclase crystal), the ⁴⁹Ti signal will suggest a lower concentration. While the lower-intensity from this minor titanium isotope will result in a larger uncertainty, the accuracy of the analysis will be better than using the sum of ⁴⁸Ti and ⁴⁸Ca. ²⁷Al and ⁴⁰Ca isotopes were used to monitor for microinclusions within the quartz as described by Kohn and Northrup [141]. Data cleaning and post-processing consisted of removal of errant values and spikes or inconsistencies in ²⁷Al and ⁴⁰Ca counts that would indicate contact with a microinclusion or grain boundary.

After data cleaning, a final, representative sample Ti concentration was calculated by taking the mean of each site; individual sites with Ti concentrations >1σ and which failed the Student’s $T$-test were considered outliers, disregarded, and not included in the final mean calculation. Analytical results are tabulated in Table 1 and Figure 10(a).

Because the TitaniQ calibration we use is sensitive to both pressure and temperature, we model Si-in-phengite isopleths to independently constrain P-T conditions based on the composition of white mica in the muscovite-celadonite solid solution [77]. The Si content per formula unit (PFU) of phengite, a high celadonite content white mica favored by high P and low T, increases with increasing pressure from 3.0 (pure muscovite) to ~3.9 (high-Si phengite).

Si-in-phengite contents were measured for each sample that has a Ti-in-quartz analysis (Table S2). Analyses were performed on a JEOL JXA-8200 electron probe microanalyzer (EPMA) at the University of California, Los Angeles. An accelerating voltage of 15 kV, beam current of 15 nA, and spot diameter of 5 μm were used throughout. Neither high-contrast backscatter imaging nor multiple EPMA line transects analyses per grain revealed compositional zoning within white mica grains. Analytical results are tabulated in Table 1 and Figure 10(b).

X-ray fluorescence (XRF) bulk rock analyses were performed for each sample that has Ti-in-quartz and Si-in-phengite measurements. Lithologically homogeneous, unweathered portions of samples were powdered in an agate ball mill. From this powder, ~2 g was separated and loss on ignition calculated. XRF analyses of major elements were measured on a Bruker M4 Tornado Micro-XRF Spectrometer at the California Institute of Technology. All P is assumed to be confined predominantly to apatite; P₂O₅ was removed from the bulk chemistry and remaining oxide ratios normalized to 100%. All Fe is assumed to be FeO (2+ oxidation state) and corresponding XRF data were corrected for this assumption; no correction was made for Fe²⁺/Fe³⁺ ratio. The results are tabulated in Table S3.

Thermodynamic modeling is carried out for the system MnO-Na₂O-CaO-K₂O-FeO-MgO–Al₂O₃-SiO₂-H₂O-TiO₂ (MnNCKFMASHT) in Perple_X [146], which calculates stable phases in P-T space based on minimization of Gibb’s free energy. We use the 2004 update to the Holland and Powell [147] thermodynamic database and the following solution models: garnet, Gt(HP); chlorite, Chl(HP); phengite, Pheng(HP) for potassic phengite; chloritoid, Ctd(HP); staurolite, St(HP) from Holland and Powell (2004); biotite, Bio(TCC) allowing for Tschermak exchange [148]; ideal solution models for hydrous cordierite (hCrd) and ilmenite-geikielite-pyrophanite (IlGkPy) from Holland and Powell [147]; high structural state feldspar (feldspar) from Fuhrman and Lindsley [149]. The system is modeled as saturated in SiO₂ and H₂O and excludes hedenbergite (hed), stilpnomelane (stlp and mnsp) as none of these phases are expected in rocks of this composition nor are they observed in the higher-grade Moine rocks. Melt is not included as P-T conditions are mostly below the minimum melting temperature, and we do not observe evidence for melting at these structural levels (Figure 9(a)).

Data used in this study are not currently available online but are available by request to anyone who wishes to access them.

The authors declare that they have no conflicts of interest.

We would like to sincerely thank Sarah Roeske for editorial handling. Our thanks to Rick Law, Jeffrey Rahl, and 2 anonymous reviewers for thorough and insightful reviews that greatly improved the content and organization of this manuscript. Special thanks are extended to Rick Law, for especially helpful reviews and discussions during the whole period of this research project. We are grateful to Richard Hervig and Lynda Williams at the SIMS facility at Arizona State University for assistance with the TitaniQ measurements. This research was funded in part by NSF grant EAR-1650173 to J. Platt and a Dornsife Doctoral Fellowship to A. Lusk.