Abstract

Despite playing a fundamental role in all models of Himalayan tectonics, minimal data constraining the structural evolution, metamorphic history, and offset magnitude of the South Tibetan detachment system (STDS) are available. Here, we integrate petrofabric, finite strain, and kinematic data with metamorphic and deformation temperatures to generate a structural model for the STDS in northwestern Bhutan. We divide the STDS into an ∼2-km-thick lower level that accommodated ∼6–13 km of thinning via ≥30–76 km of simple shear-dominant displacement within Greater Himalayan rocks, and an ∼3-km-thick upper level that accommodated ≥21 km of displacement via an upward decrease (from 44% to 2%) in transport-parallel lengthening within Tethyan Himalayan rocks. Peak metamorphic temperatures in the lower level are ∼650–750 °C, and two distinct intervals of telescoped isotherms in the upper level define a cumulative upward decrease from ∼700 to ∼325 °C. These intervals are separated by an abrupt upward increase from ∼450 to ∼620 °C, which we interpret as the result of post-STDS thrust repetition. Above the upper telescoped interval, temperatures gradually decrease upward from ∼325 to ∼250 °C through a 7-km-thick section of overlying Tethyan Himalayan rocks. Telescoped isotherms lie entirely above the high-strain lower level of the STDS zone, which we attribute to progressive elevation of isotherms during protracted intrusion of granite sills. This study demonstrates the utility of using gradients in fabric intensity and thin section-scale finite strain to delineate shear zone boundaries when field criteria for delineating strain gradients are not apparent.

INTRODUCTION

One of the most intriguing and enigmatic aspects of the structural architecture of the Himalayan orogen is the “metamorphic sandwich” of high-grade metasedimentary and igneous rocks of the Greater Himalayan (GH) package, which are situated between lower-grade metasedimentary rocks above and below (e.g., Yin, 2006; Kohn, 2014). At the top of the GH package, the top-to-north South Tibetan Detachment system (STDS) of normal faults and shear zones separates migmatitic GH rocks, which were deformed at mid-crustal levels, from overlying low-grade to un-metamorphosed, upper-crustal sedimentary rocks of the Tethyan Himalayan (TH) package (Fig. 1) (e.g., Burg and Chen, 1984; Herren, 1987; Burchfiel et al., 1992; Hodges et al., 1992; Godin et al., 2006; Kellett et al., 2018). Since its recognition, the STDS has become a key focus of Himalayan research, primarily for the kinematic challenge that it presents, requiring north-vergent displacement during construction of an overall south-vergent contractional orogenic belt. The STDS remains the most controversial first-order structure in the orogen, and plays an important (but widely varying) role in virtually every model of Himalayan orogenesis. For example, earlier studies attributed motion on the STDS to rotation of principal stress directions as a result of the topographic gradient between India and Tibet (Burchfiel and Royden, 1985), or as a consequence of surface uplift combined with cessation of slip on the underlying Main Central thrust (England and Molnar, 1993). More recently, the STDS has been interpreted as a structure that helped accommodate high-magnitude southward extrusion of a mid-crustal channel composed of GH rocks (e.g., Grujic et al., 1996; Beaumont et al., 2001; 2004; Jamieson et al., 2004), a normal fault system generated within an orogenic wedge governed by critical taper dynamics (e.g., Robinson et al., 2003, 2006; Kohn, 2008), a structure related to the dynamics of translation and structural thickening of the underlying GH section (Cottle et al., 2015; Larson et al., 2015), or as a backthrust system bounding the top of a southward-propagating tectonic wedge (Webb et al., 2007, 2011; He et al., 2015).

A critical prerequisite for interpreting the genesis of the STDS is an understanding of the structural and metamorphic field gradients that are the cumulative result of north-vergent shearing, and the implications of these gradients for offset magnitude. To date, few studies have documented the magnitude of isotherm telescoping across the full thickness of the STDS (Cottle et al., 2011; Law et al., 2011; Kellett and Grujic, 2012). In this study, we present data that illuminate trends in peak and deformation temperature, fabric development, finite strain, and kinematics along a transect through the STDS in northwestern Bhutan. These data facilitate definition of the boundaries and internal architecture of the STDS, allow estimation of offset magnitude, and are integrated with published geochronology to produce a detailed model of the structural evolution of the shear zone. We then use this case study to argue for the criteria that are most appropriate for definition of the boundaries of the STDS and other shear zones, and to demonstrate the utility of gradients in fabric intensity and finite strain when field criteria for delineating shear zone boundaries are not easily identifiable.

TECTONIC FRAMEWORK

Himalayan Geologic Background

The closure of the Neotethys Ocean culminated in collision of the Indian and Asian plates at ca. 50–55 Ma (e.g., Rowley, 1996; Leech et al., 2005; Najman et al., 2010). Following initial collision, deformation associated with continued India-Asia convergence has constructed the Himalayan-Tibetan orogen. South of the Indus-Tsangpo suture (Fig. 1A), rocks native to the Indian plate have been buried, metamorphosed, and stacked into a south-vergent thrust belt (e.g., Powell and Conaghan, 1973; LeFort, 1975; Mattauer, 1986). First-order structures in the thrust belt separate packages of rock with differing deformation histories and metamorphic grades. Immediately south of the suture, the TH package consists of low-grade to unmetamorphosed, Paleozoic–Mesozoic sedimentary rocks that were deposited on the southern Tethyan passive margin (e.g., Brookfield, 1993; Garzanti, 1999). To the south, the north-vergent, normal-sense STDS places TH rocks over high-grade (typically upper amphibolite-facies) metasedimentary and meta-igneous rocks of the GH package (Fig. 1A) (e.g., Burg, 1983; Herren, 1987; Burchfiel et al., 1992). GH rocks were buried to mid-crustal depths as early as the late Eocene–Oligocene (e.g., Hodges et al., 1996; Godin et al., 2001; Corrie and Kohn, 2011; Zeiger et al., 2015), and experienced partial melting and granite intrusion during the Miocene (e.g., Harrison et al., 1998). GH rocks overlie greenschist-facies metasedimentary rocks of the Lesser Himalayan package across the Main Central thrust, a major south-vergent shear zone that was active in the early Miocene (e.g., Gansser, 1964; Robinson et al., 2003; Kohn et al., 2004; Tobgay et al., 2012). The southern part of the orogen is defined by a brittle thrust belt that deforms Lesser Himalayan rocks into regional-scale duplex systems (e.g., Robinson et al., 2006; Long et al., 2011a; Webb, 2013). Near the modern deformation front, Miocene–Pliocene synorogenic rocks of the Subhimalayan package are involved in deformation (e.g., DeCelles et al., 1998; Ojha et al., 2000; Lavé and Avouac, 2000).

GH Rocks, TH Rocks, and the STDS in Northwestern Bhutan

The GH section in northwestern Bhutan consists of ∼8 km of migmatitic, garnet- and sillimanite-bearing paragneiss, with interlayered schist, quartzite, and orthogneiss (Gansser, 1983; Kellett et al., 2009; Long et al., 2011c; Warren et al., 2011; Tobgay et al., 2012). On the basis of youngest detrital zircon peaks and crystallization ages of granite intrusions, deposition of GH sedimentary protoliths in Bhutan is bracketed between the Neoproterozoic and Ordovician (Long and McQuarrie, 2010; Tobgay et al., 2012; McQuarrie et al., 2013). Above the GH rocks, an ∼50 km (E-W) by ∼40 km (N-S) exposure of TH rocks, often referred to as the Lingshi klippe or syncline (e.g., Kellett and Grujic, 2012; Tobgay et al., 2012), is preserved (Fig. 1C). At its base, the lithologically heterogeneous Chekha Formation is mapped as the lowest TH unit (e.g., Jangpangi, 1978; Gansser, 1983; Tangri and Pande, 1995; Grujic et al., 2002). In the western part of the Lingshi exposure, the Chekha Formation consists of biotite-garnet schist and calc-silicate gneiss, but marble is the dominant lithology in the eastern part (Gansser, 1983; Kellett and Grujic, 2012). In central Bhutan, deposition of the Chekha Formation is bracketed as Middle to Late Cambrian by ca. 515 Ma youngest detrital zircons (McQuarrie et al., 2013) combined with ca. 495 Ma trilobite fossils documented in overlying TH rocks (Hughes et al., 2011). Above the Chekha Formation, Paleozoic TH rocks in the Lingshi exposure include Cambrian and Carboniferous–Permian diamictite, quartzite, limestone, and slate (Gansser, 1983; Bhargava, 1995; Tangri and Pande, 1995). Higher in the section, Triassic–Jurassic and Cretaceous slate is exposed (Ganesan and Bose, 1982; Gansser, 1983).

At the Lingshi exposure, the STDS has been mapped as a north-vergent ductile shear zone that places the Chekha Formation over GH rocks (e.g., Grujic et al., 2002). Kellett et al. (2009) referred to this exposure, and to other STDS exposures mapped farther to the east in Bhutan, as the “outer STDS,” to contrast them with the main exposure of the STDS that has been mapped along the Bhutan-Tibet border (the “inner STDS”) (Fig. 1B) (Burchfiel et al., 1992; Edwards et al., 1996; 1999). Kellett et al. (2010) interpreted the outer STDS at the Lingshi exposure as a 2–3-km-thick high-strain zone defined by a dominantly top-to-north shear-sense, extending as much as ∼1 km below and ∼2 km above the GH-Chekha Formation contact.

GH rocks collected within the STDS zone attained peak pressure-temperature conditions of ∼6–9 kbar and ∼620–790 °C (obtained using THERMOCALC; Kellett et al., 2010). Raman spectroscopy of carbonaceous material (RSCM) thermometry defines peak temperatures of ∼300–450 °C at the top of the Chekha Formation, and ∼200–300 °C for Paleozoic and Mesozoic TH rocks (Kellett and Grujic, 2012). GH rocks experienced prograde monazite growth between ca. 26 and 21 Ma (Kellett et al., 2010). Monazite inclusions within garnet that grew prior to STDS-related deformation cluster at ca. 23 Ma, which is interpreted as the timing of partial melting of GH rocks (Kellett et al., 2010). Partial melting and peak metamorphism were followed by north-vergent shearing on the STDS, accompanied by retrograde growth of high-Y monazite rims in GH rocks between ca. 19 and 15 Ma (Kellett et al., 2010) and crystallization of granite sills within the Chekha Formation as young as ca. 16.5 Ma (Kellett et al., 2009). 40Ar/39Ar muscovite thermochronology from Triassic–Jurassic TH rocks near Lingshi and from a granite sill in the Chekha Formation indicate cooling through ∼425 °C at ca. 13–11 Ma (Kellett et al., 2009; Antolín et al., 2012).

On the northwestern margin of the Lingshi exposure, TH rocks are offset by the Lingshi fault (Fig. 1C), a SE-dipping brittle fault with normal and sinistral motion components (Gansser, 1983). The Lingshi fault has been correlated with the Yadong cross-structure, a series of NE-striking, active normal faults in southern Tibet (Fig. 1B) (Armijo et al., 1986; Wu et al., 1998). Motion on the Lingshi fault may have begun as early as ca. 14 Ma (Cooper et al., 2015).

TECTONOSTRATIGRAPHY OF THE DODENA-LINGSHI TRANSECT

In order to define trends in temperature, finite strain, and kinematics across the STDS in northwestern Bhutan, we mapped an ∼50-km-long transect between the towns of Dodena and Lingshi (Figs. 1C and 2). A total of 141 outcrops were examined, with samples collected at 54 of these localities (see Data Repository Table DR11). Rocks at all levels exhibit a dominant macroscopic foliation (Figs. 1 and 2), and mineral stretching lineations were observed in most outcrops and typically plunge shallowly NW (mean trend of 323°; Fig. 1C inset). Structural (i.e., foliation-normal) thicknesses were measured by projecting the apparent dip of foliation measurements onto a cross section (Fig. 3). Field stops and sample localities are projected onto a tectonostratigraphic column (Fig. 4), with structural distance listed relative to the base of the Chekha Formation. All structural distances in the following text are listed in meters below (negative numbers) or above (positive numbers) the base of the Chekha Formation. The transect encompasses a total thickness of 14,100 m, ranging from −3950 to +10,150 m. Lithologies and other features on the transect are described below from structurally low to high.

GH rocks: Between −3950 and −950 m below the base of the Chekha Formation, exposures are dominated by migmatitic (typically <5%–10% leucosomes by volume), sillimanite- and garnet-bearing paragneiss, which is interlayered with quartzite and schist. Between −950 and −650 m, deformed, foliation-parallel granite intrusions are interlayered with calc-silicate gneiss, quartzite, and paragneiss. Between −650 and 0 m, migmatitic orthogneiss with up to 30%–40% leucosomes by volume is dominant and is interlayered with migmatitic, sillimanite-garnet paragneiss. The structurally highest observed occurrence of sillimanite is at −25 m.

Chekha Formation (5550 m thick): We map the base of the Chekha Formation at an upward transition from orthogneiss to calc-silicate gneiss, which is similar to Gansser (1983) and Kellett and Grujic (2012). The basal 300 m of the Chekha Formation consists of calc-silicate gneiss interlayered with schist and quartzite. The structurally highest occurrence of leucosomes is at +75 m above the base of the Chekha Formation, and the structurally highest occurrence of garnet is at +250 m. The Chekha Formation is dominated by marble between +300 and +5550 m, with intervals of quartzite and phyllite between +1350 and +3200 m, and intervals of slate between +4000 and +5125 m. The structurally highest observed occurrences of biotite and muscovite porphyroblasts are at +1275 m and +4650 m, respectively.

Granite sills: Deformed granite sills are observed through the full examined thickness of GH rocks, and as high as +3200 m in the Chekha Formation (Fig. 4). However, they are most prevalent between −1800 and +800 m, where they locally account for up to 30%–40% of outcrops by volume (Figs. 5C, 5D). These sills typically range in thickness from ∼0.1 to 2 m, have a solid-state foliation, and often exhibit boudinage (Fig. 5A, 5B).

Paleozoic TH rocks (2900 m thick): Following Gansser (1983), we map the contact between the Chekha Formation and overlying Paleozoic TH rocks at an upward transition from marble to quartzite and slate. Paleozoic TH rocks are observed between +5550 m and +8450 m, and are dominated by slate (note: structural thicknesses were measured normal to slaty cleavage, and not to sedimentary bedding, which is locally preserved in TH slates), with a 200-m-thick basal quartzite, limestone intervals at +5825–5875 m and +7075–7350 m, and diamictite at +6925–7075 m (Fig. 4).

Mesozoic TH rocks (>1700 m thick): Mesozoic TH rocks are mapped after Gansser (1983), and are dominated by gray to black, graphitic slate. They exhibit a minimum thickness of 1700 m, between +8450 and +10,150 m above the base of the Chekha Formation. Near Lingshi, the axis of a NE-trending open syncline, which can be traced for at least ∼25 km (Fig. 1C), exposes the highest TH rocks on the transect.

Lingshi fault: The northern end of the transect crosses the NE-striking Lingshi fault, which places Cretaceous TH rocks in the hanging wall against Chekha Formation marble in the footwall. This relationship requires at least ∼6 km of down-to-SE, normal-sense offset along the transect (Fig. 3). However, map relationships show that offset on the Lingshi fault decreases significantly along strike to the NE and SW (Gansser, 1983; Kellett and Grujic, 2012; this study).

STRUCTURAL OBSERVATIONS

Shear-Sense Indicators

Outcrop- and thin section-scale shear-sense indicators on the Dodena-Lingshi transect (Fig. 6) include C′-type shear bands, SC fabrics, mica fish, oblique foliations in recrystallized quartz, asymmetric folds, asymmetric boudinage, and leucosomes, granite sills, and rigid porphyroclasts sheared into σ-objects. Shear-sense indicators were observed within a total of 29 outcrops and 15 thin sections (Fig. 4). Though both top-to-SE and top-to-NW indicators are widely distributed through the transect, top-to-NW shear is dominant between −1825 and +775 m (Fig. 4). Additional shear-sense information demonstrated by quartz c-axis plots and stretching directions in finite strain analyses are shown on Figure 4, and are discussed below.

Quartz Recrystallization Microstructures

Dynamically recrystallized quartz was observed in thin sections of all GH samples and the lowest 2700 m of the Chekha Formation (Fig. 7). Here, we used the morphology of recrystallized quartz to interpret dominant recrystallization mechanisms and to estimate deformation temperature (e.g., Stipp et al., 2002; also see review in Law, 2014). We utilized the approximate deformation temperature ranges for quartz recrystallization mechanisms calibrated by Law (2014; his fig. 21) from Himalayan rocks in proximity to the Main Central thrust (Fig. 4) (as a caveat, differences in strain rate between the Main Central thrust and the STDS may also affect the apparent deformation temperature ranges, as recrystallization mechanisms are also sensitive to strain rate; e.g., Law, 2014). Four different recrystallization microstructures were observed: (1) chessboard extinction (CBE), characterized by extinction domains that intersect at approximately right angles within large (≥∼0.25 mm) amoeboid quartz neoblasts (a neoblast is defined here as a visually distinct recrystallized grain) (Fig. 7A; e.g., Lister and Dornsiepen, 1982; Mainprice et al., 1986), which indicates deformation at ≥∼630 °C (Stipp et al., 2002); (2) grain boundary migration (GBM), which is indicated by ∼0.1–1.0 mm amoeboid quartz neoblasts (Fig. 7B; e.g., Guillope and Poirier, 1979; Urai et al., 1986), and represents deformation at ∼550–650 °C (Law, 2014); (3) subgrain rotation (SGR), which is defined by equigranular, polygonal, ∼10–100 µm quartz neoblasts (Fig. 7C; e.g., Poirier and Nicolas, 1975; White, 1977; Guillope and Poirier, 1979), and indicates deformation at ∼450–550 °C (Law, 2014); and (4) bulging (BLG), which is indicated by development of <∼10 µm subgrains at quartz porphyroclast boundaries (a porphyroclast is defined here as a relict, detrital grain) (Fig. 7D; Bailey and Hirsch, 1962; Drury et al., 1985), and represents deformation at ∼350–450 °C (Law, 2014).

On the Dodena-Lingshi transect, an overall pattern of recrystallization at progressively lower temperatures with increasing structural level is defined (Fig. 4). CBE was observed between −3950 and −1325 m (n = 5), GBM between −1825 and −25 m (n = 8), SGR between −1025 and +1975 m (n = 18), and BLG at +2700 m (n = 1). The average diameter of recrystallized quartz grains was measured for seven quartz-rich samples (see Fig. DR7 [footnote 1] for a graph of grain size versus structural height), which are distributed between −1550 and +150 m. Samples dominated by SGR recrystallization yielded diameters between ∼40 and ∼100 µm, and those dominated by GBM recrystallization yielded diameters between ∼100 and ∼230 µm. However, these data did not yield a consistent pattern of grain size with structural height.

RSCM THERMOMETRY

Carbonaceous material was analyzed from thin sections of GH metasedimentary rocks (n = 4), the Chekha Formation (n = 14), and Paleozoic-Mesozoic TH rocks (n = 3), in order to calculate peak metamorphic temperatures using the RSCM thermometer (e.g., Beyssac et al., 2002, 2003; Rahl et al., 2005). Measurements were made at Arizona State University, Tempe, Arizona, USA (see Data Repository [footnote 1] for analytical methods and supporting data). Carbonaceous material was typically present as ∼5–25 µm patches (Figs. 8A–8C), and sometimes as laminations (Fig. 8D). Examples of representative Raman spectra for each sample are shown on Figure 8E. We followed the analytical methods of Rahl et al. (2005), in which the RSCM thermometer for rocks that achieved peak temperatures between ∼100 and 740 °C was determined by measuring the height ratio (R1) and area ratio (R2) of four first-order Raman peaks (G, D1, D2, D3) in the relative wavenumber range of 1200–1800 cm–1. Mean peak temperatures of multiple measurements (varying between 9 and 16 for individual samples) are shown on Table 1, and are reported at a 2 standard error (SE) level, which takes into account internal uncertainty and the external uncertainty from the Rahl et al. (2005) calibration (e.g., Cooper et al., 2013; Long et al., 2016, 2017).

Four GH samples, which span structural heights between −3950 and −25 m, yielded overlapping temperatures between ∼650 and 700 °C (Fig. 4). The lowest seven Chekha Formation samples define an overall upward decrease from ∼650 to ∼450 °C between +175 and +1275 m. Above this, four Chekha Formation samples between +1975 and +2700 m define an abrupt increase to ∼620 °C, followed by an upward decrease to ∼480–510 °C. Moving structurally upward, three Chekha Formation samples between +3275 and +4425 m overlap between ∼320 and 360 °C, and three TH samples between +6750 and +8775 m overlap between ∼270 and 280 °C.

The 21 new RSCM temperatures from this study were combined with published peak temperatures from 12 samples from the Dodena-Lingshi transect and proximal areas (Table 2; Fig. 4). These include THERMOCALC temperatures determined from four GH samples (Kellett et al., 2010), which yielded temperatures of ∼740–790 °C between −1000 and −50 m, and ∼620 °C at −20 m. In addition, Kellett and Grujic (2012) presented RSCM temperatures from three Chekha Formation samples (+3050–5275 m), which yielded temperatures of ∼300–320 °C, two Paleozoic TH samples (+5700, +6075 m), which ranged from ∼250 to 290 °C, and three Mesozoic TH samples (+9550–9925 m), which ranged from ∼220 to 270 °C.

Combining the new and published temperature data sets yielded the following trends: Between −3950 and −50 m, peak temperatures are between ∼650 and 750 °C (n = 5), with one outlier of ∼790 °C (Fig. 4). Between −50 and +1275 m, temperatures decrease upward from ∼625 to 725 °C to ∼450 °C (n = 9, with one outlier (sample 113 at +450 m) that falls below this trend), which is best-fit by a telescoped, upright thermal field gradient of 163 ± 38 °C/km (R2 = 0.70). There is a data gap between samples at +1275 and +1975 m, across which temperatures increase upward from ∼450 to ∼620 °C. Above this, between +1975 and +3050 m, temperatures again decrease upward from ∼620 to ∼320 °C (n = 5), which is best-fit by a steep, upright field gradient of 257 ± 63 °C/km (R2 = 0.85). Above this, between +3050 and +9925 m, temperatures gradually decrease upward from ∼320–360 °C to ∼220–265 °C (n = 14), which is best-fit by a gentle, upright field gradient of 13 ± 3 °C/km (R2 = 0.70).

QUARTZ CRYSTALLOGRAPHIC FABRICS

In order to provide information on shear-sense, the 3-D strain field, deformation temperature, and kinematic vorticity (e.g., Lister and Hobbs, 1980; Schmid and Casey, 1986; Wallis, 1992, 1995; Kruhl, 1998; Morgan and Law, 2004), the orientations of quartz c-axes were measured in seven quartzite samples collected from structural levels between −1550 and +150 m (Table 3; Fig. 9). Quartz c-axis orientations were measured on foliation-normal, lineation-parallel thin sections using an automated crystal fabric analyzer (e.g., Wilson et al., 2009; Peternell et al., 2011; Larson et al., 2013; Larson and Cottle, 2014) (see Data Repository for details on analytical methods). The c-axis measurements are plotted on lower-hemisphere, equal-area stereonets in Figure 9.

All samples exhibit crystallographic preferred-orientations (CPO) characteristic of recrystallization at <∼650 °C (e.g., Lister and Hobbs, 1980; Behrmann and Platt, 1982; Bouchez et al., 1983). Samples 3 and 14 yielded weak type I crossed-girdle patterns and sample 82 exhibits a type II crossed-girdle pattern, which indicate that basal <a>, rhomb <a>, and prism <a> slip systems were active (e.g., Lister, 1977; Schmid and Casey, 1986). Sample 7 yielded a poorly developed crossed-girdle pattern, but its type (I or II) could not be visually determined. Samples 10 and 21A exhibit single-girdle patterns, which result from prism <a> and rhomb <a> dominant slip (e.g., Lister and Hobbs, 1980; Schmid and Casey, 1986). Sample 13 yielded a center point maxima, indicating prism <a> dominant slip (e.g., Mainprice et al., 1986; Passchier and Trouw, 2006). The fabric patterns from these seven samples are indicative of deformation within a 3-D strain field that approximates plane strain (e.g., Schmid and Casey, 1986).

Statistical parameters can be used to quantify the intensity (i.e., non-randomness) of fabric development, and have recently been applied as strain magnitude proxies to characterize strain patterns across Himalayan shear zones (Larson et al., 2017; Starnes et al., 2017). Here, we utilize the cylindricity index (B) of Vollmer (1990), which is based on eigenvector analysis of c-axis directions (e.g., Woodcock, 1977; Mainprice et al., 2014). Vollmer (1990) divided quartz fabrics into three end-member components: point (P), girdle (G), and random (R), where P + G + R = 1. The formation of a strong P or G fabric depends upon the active slip system during quartz recrystallization (Barth et al., 2010), and B is defined as the sum of the P and G values. B is a measure of the non-randomness of a fabric, where 0 represents a completely random fabric and 1 represents a completely non-random fabric. P, G, R, and B values (Table 3) were calculated from the quartz c-axis measurements using the program Orient (Vollmer, 2017 [Orient 3: Spherical projection and orientation data analysis program, www.frederickvollmer.com]) (see Fig. DR6 [footnote 1] for a ternary plot of P, G, and R values). B values on the Dodena-Lingshi transect (Fig. 9D) are low (0.21; sample 3) at −1550 m, increase to values between 0.48 and 0.87 between −1325 and 0 m (samples 7, 10, 13, 14, 21A), and decrease to 0.10 at +150 m (sample 82).

Assuming a consistent, critically resolved shear stress, negligible hydrolytic weakening, and an invariant strain rate, the opening angle exhibited in crossed-girdle quartz c-axis fabrics may be related to the temperature at which the fabrics were “locked in” during the final stages of dynamic recrystallization (e.g., Kruhl, 1998; Morgan and Law, 2004; Faleiros et al., 2016). Here, we use the pressure-independent calibration of Faleiros et al. (2016) to estimate opening angle deformation temperatures from four samples that yielded crossed-girdle fabrics. The samples exhibited a total range of opening angles between 61° and 83° (Fig. 9; Table 3), which correspond to deformation temperatures of 490 °C (sample 3), 579 °C (sample 7), 469–531 °C (sample 14), and 586–621 °C (sample 82).

Four samples exhibit asymmetric c-axis fabrics, which allow for interpretation of shear-sense. The central part of the fabric skeletons for samples 3, 10, 14, and 21A are inclined to the NW relative to foliation, defining a top-to-NW shear sense (Fig. 9B). This is consistent with top-to-NW microstructural shear-sense indicators observed in thin sections of all four samples.

Using the quartz shape-preferred orientation method of Wallis (1992, 1995), the angle between foliation and the line normal to the central girdle segment of a quartz c-axis fabric (β) can be combined with the mean elongation angle of recrystallized quartz neoblasts measured relative to foliation (θ′ISA1) to calculate the mean kinematic vorticity number (Wm). Wm is a ratio that quantifies the relative contributions of pure shear and simple shear (e.g., Means et al., 1980; Passchier, 1987; Tikoff and Fossen, 1993; Means, 1994); a value of 0 represents entirely pure shear, 1 represents entirely simple shear, and equal contributions occur at 0.71 (e.g., Law et al., 2004). Here, we assume the simplified case of steady-state plane strain for estimation of Wm (e.g., Fossen and Tikoff, 1993; Johnson et al., 2009), which is supported by fabric patterns indicative of approximate plane strain conditions observed between −1550 and +150 m. β and θ′ISA1 angles were measured for GH samples 3, 10, 14, and 21A (Fig. 9; Table 3; see Data Repository [footnote 1] for supporting data), which yielded Wm ranges of 0.74–0.83, 0.81–0.88, 0.53–0.64, and 0.74–0.83, respectively. This corresponds to a range of 55%–70% simple shear, with one outlier of 35%–45% (sample 14). β and θ′ISA1 angles were measured for Chekha Formation sample 82, which yielded a Wm range of 0.03–0.17 (5%–15% simple shear, 85%–95% pure shear).

FINITE STRAIN DATA

All GH samples exhibit complete recrystallization of quartz (e.g., Fig. 7). However, between +175 and +2700 m, several Chekha Formation samples contain non-recrystallized, elongated quartz porphyroclasts that are isolated within a mica- or calcite-rich matrix (Fig. 10A), and quartz within all Chekha Formation and TH samples above +2700 m is not recrystallized (Fig. 10B). We investigated the magnitude and orientation of 3-D finite strain in 19 Chekha Formation samples and 6 TH samples distributed between +175 and +9675 m by performing the Rf-ϕ method (e.g., Ramsay, 1967; Dunnet, 1969) on elongated, non-recrystallized quartz porphyroclasts (supporting data in the Data Repository). Two foliation-normal thin sections were cut from each sample, one parallel (denoted with “A”) and one normal (denoted with “B”) to stretching lineation, which are interpreted to approximate the XZ and YZ strain planes, respectively. 2-D strain ellipses from each thin section are plotted on Figure 4, and strain data are summarized in Table 4. For all samples, the 2-D strain ratio (Rs) in the “A” thin section was either greater than or equivalent within error to Rs in the “B” thin section, which supports the use of foliation and lineation to approximate the principal strain directions. In addition, the quartz c-axis fabrics on Figure 9, by analogy with experiments and numerical modeling of fabric development, also demonstrate that stretching lineation in the analyzed samples is parallel to the maximum finite stretch (X) direction (e.g., Lister et al., 1978; Lister and Hobbs, 1980).

Competency contrasts between matrix and clasts can result in heterogeneous strain partitioning (e.g., Sanderson, 1982; Holyoke and Tullis, 2006), and quartz clasts often exhibit lower elongation magnitudes than the surrounding mica-rich matrix (e.g., Tullis and Wenk, 1994; Treagus and Treagus, 2002; Yonkee et al., 2013). Therefore, we interpret that our strain magnitudes likely represent minimum values.

Our data show an overall decrease in strain magnitude with increasing structural level (Fig. 4), and can be divided into three vertical domains: (1) the lowest 5 Chekha Formation samples (+175 to +775 m), which yielded average Rs[X/Z] and Rs[Y/Z] ratios of 2.60 ± 0.10 (errors at 1 SE level) and 1.85 ± 0.20, respectively; (2) Chekha Formation samples between +850 and +4650 m (n = 12), which yielded average Rs[X/Z] and Rs[Y/Z] ratios of 1.80 ± 0.10 and 1.50 ± 0.05; and (3) Chekha Formation and TH samples between +4700 and +9675 m (n = 8), which yielded average Rs[X/Z] and Rs[Y/Z] ratios of 1.25 ± 0.05 and 1.20 ± 0.05.

On a Flinn diagram (Fig. 10C), most ellipsoids plot entirely in the flattening field (n = 15), with nine samples that exhibit uncertainty ranges that overlap the plane strain line, and one sample that plots entirely in the constrictional field. ϕ angles were measured relative to foliation and are therefore equivalent to the parameter θ′ of Ramsay and Huber (1983) (from this point on ϕ is referred to as θ′; see discussion DR5 (see footnote 1) for additional details on the orientation relationship between macroscopic foliation and shear direction). The sign convention used for θ′ is positive if the grain long axis is inclined toward the NW or NE relative to foliation, and negative if inclined toward the SE or SW. Thirty-nine of the 50 thin sections yielded θ′ values between ±10°. The remaining 11 thin sections yielded values between ±11–19°, with one outlying value of +42°. The Flinn diagram results and low θ′ values together indicate that most samples underwent layer-normal flattening strain, with some approaching plane strain conditions.

We estimated mean kinematic vorticity numbers (Wm) within “A” thin sections by comparing our data to Wm contours plotted on a graph of Rs versus θ′ (Fig. 10D) (referred to here as the Rs-θ′ method; e.g., Tikoff and Fossen, 1995) (refer to discussion DR5 for evidence supporting assumptions for estimating Wm using this method). This method assumes that quartz clasts were homogeneously elongated in the direction of maximum finite stretching, which is supported by low SE values for Rs (typically ±0.1–0.2) and θ′ (typically ±3–6°) for each analysis. For estimation of Wm, idealized plane strain is assumed (e.g., Fossen and Tikoff, 1993; Johnson et al., 2009), and since most of the analyzed rocks experienced flattening strain, Wm values are here interpreted to represent maxima (e.g., Tikoff and Fossen, 1995). Given the relatively low strain magnitudes of our samples (total Rs range of 1.1–2.9), the overestimation of Wm likely does not exceed ∼0.05 (Tikoff and Fossen, 1995). The range of Wm values reported for each sample was estimated from the ±1 SE θ′ error range (Fig. 10D; Table 4). Most samples (n = 21) yielded Wm values between 0.00 and 0.45 (70%–100% pure shear), with four outliers that yielded values as high as 0.65–0.85 (as little as 35%–55% pure shear).

The mean elongation direction of quartz porphyroclasts in “A” thin sections, as denoted by the sign of θ′, can also be used to determine shear-sense (e.g., Ramsay, 1967; Passchier and Trouw, 2006). Under our sign convention, negative θ′ values correspond to a top-to-NW shear-sense, and positive to top-to-SE (Fig. 10D). Most samples (n = 16) exhibit low (≤ ± 4°) θ′ values with error ranges that overlap with foliation, which we do not interpret as robust shear-sense indicators. However, the remaining nine samples exhibit θ′ values with error ranges that do not overlap with foliation. These include three Chekha Formation samples between +450 and +775 m (113, 89, 91) which define a top-to-NW shear sense, five Chekha Formation and TH samples between +1975 and +8775 m (97, 101, 46, 37, 36) that define a top-to-SE shear-sense, and a TH sample at +5700 m (45) that yielded a top-to-NW shear sense (Fig. 4).

Using methods outlined in Law (2010) and Xypolias et al. (2010), strain magnitude and Wm were integrated to calculate transport-parallel lengthening and flow plane-normal shortening (which is normal to the transport direction, and is from here on referred to as transport-normal shortening) for individual strain analyses (see Data Repository for additional information) (Table 4; Figs. 10E–10F). The data define an upward decrease in average lengthening and shortening for each strain domain. Domain 1 yielded an average lengthening value of 44 ± 5% (1 SE) and an average shortening value of 34 ± 2%. Domain 2 yielded average lengthening and shortening values of 20 ± 3% and 21 ± 2%, respectively. Domain 3 exhibits very low values, with 2 ± 1% lengthening and 3 ± 2% shortening. The average lengthening in the Y direction varies between 2% and 8% for all three strain domains.

DISCUSSION

Integration and Interpretation of Temperature and Structural Data Sets

The temperature, fabric, strain, and kinematic data sets allow generation of a structural model for the STDS in northwestern Bhutan. We interpret the base of the STDS zone, which we define as the section of rocks affected by deformation related to STDS shearing, at the onset of top-to-NW dominated shearing at −1825 m (Figs. 4 and 10G). This level coincides with the approximate base of prevalent (up to 30–40 vol%) granite sills. Quartz fabrics between −1550 and 0 m indicate approximately plane strain deformation, with well-developed CPO’s that increase in cylindricity from 0.21 at −1550 m to 0.48–0.87 between −1325 and 0 m. Therefore, though no finite strain data are available from this interval, the gradient in fabric intensity observed between quartz-rich samples of similar lithology implies an upward-increasing strain gradient (Fig. 11B) (e.g., Larson et al., 2017). Most Wm values in this interval range from 0.74 to 0.88, indicating a dominance (55%–70%) of simple shear. Published pressure data from GH samples at −1000 m (9.0 ± 0.8 kbar), −550 m (6.9 ± 1.5 kbar), and −50 m (6.2 ± 1.7 kbar) (Kellett et al., 2010) define a best-fit, upright field gradient of 2.92 ± 0.94 kbar/km (R2 = 0.91). This is ∼7–14 times steeper than a typical lithostatic gradient for mid-crustal rocks of 0.27 kbar/km, and corresponds to 86%–93% apparent vertical shortening (Fig. 11B). This shortening range is similar to other estimates across the STDS in the Everest region (90%–94%; Law et al., 2011) and in northwestern India (96%; Herren, 1987). The pressure data indicate that at least ∼6–13 km of apparent structural thinning was accommodated within the basal 1.8 km of the STDS zone. When integrated with Wm data, this thinning appears to have been accomplished dominantly by simple shear (55%–70%) with a less significant component of shortening normal to the shear zone boundaries. To accomplish this, the shear zone must have dipped northward while it was active (e.g., Law et al., 2011). Using a simple geometric model that attributes 55%–70% of the structural thinning to simple shear, offset ranges of 30–38 km and 59–76 km can account for the pressure difference between GH samples at −1000 and −50 m for a shear zone active at a dip angle of 10° and 5°, respectively. This corresponds to a range of shear strain values between 5.7 and 11.4. These offset ranges should be considered minima, as they were calculated to account for the pressure difference between samples that do not encompass the complete thickness of the basal part of the STDS zone.

Between −1100 and +700 m, many samples exhibit evidence for quartz recrystallization at temperatures up to ∼100 to ∼250 °C lower than peak temperatures (Fig. 4), as constrained by quartz fabric opening angles and semiquantitative bracketing using quartz recrystallization microstructures. Because of the difference between peak and deformation temperatures, we interpret these fabrics to have been generated along the exhumation path (e.g., Law, 2014), and therefore likely represent recrystallization during the later stages of shearing on the STDS. This interpretation is further supported by the dominance of top-to-NW shear-sense indicators and quartz fabrics within this interval, which also argue for recrystallization along the exhumation path (e.g., Law et al., 2011).

The peak metamorphic temperature data define two intervals with telescoped isotherms, one between −50 and +1275 m, and the other between +1975 and +3050 m (Fig. 4). We interpret the upward increase from ∼450 to ∼620 °C between +1275 and +1975 m to be the result of structural repetition of part of the Chekha Formation by an intraformational thrust fault or thrust-sense shear zone that post-dates STDS movement. Several lines of evidence support this interpretation, including: (1) the Chekha Formation is anomalously thick along this transect compared to other areas of the Lingshi exposure or farther east in Bhutan, where it typically ranges between ∼1000 and 2000 m thick (Gansser, 1983; Kellett et al., 2009, 2010; Long et al., 2011a, 2011c; Kellett and Grujic, 2012); (2) Kellett and Grujic (2012) inferred a S-vergent thrust fault in their cross section through the Linghsi exposure (their fig. 3B), which projects to an approximate structural height of +1650 m and therefore lies within the location range proposed here; (3) Gansser (1983) also speculated that thrust faults had thickened the Chekha Formation marble along this transect; and (4) south-vergent thrusting that is interpreted to post-date shearing on the STDS has been documented across northern Bhutan (Grujic et al., 2002, 2011; Kellett et al., 2009; Warren et al., 2011). Therefore, we suggest that the observed thermal pattern is most likely the result of thrust repetition of telescoped isotherms produced during earlier shearing on the STDS. By restoring sample 97 to the level of sample 89, which both have similar peak temperatures of ∼625 °C, structural repetition of ∼1300 m of the Chekha Formation is estimated.

We interpret the boundary between strain domains 2 and 3 at ∼4.7 km above the base of the Chekha Formation, which represents an upward transition to very low strain magnitudes (Rs ∼1.25), as the approximate top of the STDS zone (Fig. 11B). This defines a present-day thickness of 6.5 km for the STDS zone, which restores to an original thickness of ∼5.2 km after accounting for post-STDS thrust repetition. The cumulative temperature decrease across the STDS zone is ∼300–400 °C (from ∼625–725 to ∼325 °C). This is similar to temperature patterns documented across the STDS in the Everest region, where Cottle et al. (2011) measured an upward decrease from ∼650 to ∼340 °C across an ∼1-km-thick STDS zone, and Law et al. (2011) measured an upward decrease from ∼680 to ∼460 °C across an ∼0.5-km-thick STDS zone.

The finite strain data from TH rocks define a downward-increasing gradient in transport-parallel lengthening, which indicates that a component of top-to-NW displacement is broadly distributed through the upper part of the STDS zone (Fig. 11A, 11B). This is compatible with shear-sense observations, which indicate top-to-NW shearing as high as +3–7 km (Fig. 10G). Between the base of strain domain 3 (+4.7 km) and the top of domain 1 (+0.8 km), a 42% difference in lengthening is observed. A minimum NW-SE map length of ∼50 km of the Chekha Formation is preserved in northwestern Bhutan (Fig. 1C). Using this length, and assuming that strain magnitudes were similar across-strike, this corresponds to a minimum of ∼21 km of differential lengthening distributed through the upper part of the STDS zone. Therefore, the upper part of the STDS acted as a stretching fault (e.g., Means, 1989; Grujic et al., 2002), and pure shear-dominated lengthening contributed at least a few 10s of km to its overall offset magnitude.

The finite strain and vorticity data from TH rocks, when combined with the fabric, vorticity, and pressure data from GH rocks, can be used to divide the STDS zone into two distinct levels: (1) a lower level (−1825 to +175 m) dominated by simple shear, which accommodated significant apparent thinning (>6–13 km); and (2) an upper level (+175 to +4675 m) characterized by distributed, pure shear-dominated lengthening and shortening, which both progressively decrease upward (lengthening from 44% to 2%, and shortening from 34% to 3%) (Fig. 11B). Offset magnitude varies significantly between the two levels, with ∼20 km of distributed top-to-NW displacement accommodated in the upper level, and at least ∼30–75 km accommodated by simple shear in the lower level. These estimates are compatible with finite strain data in central Bhutan, which allow for up to ∼85 km of distributed N-vergent shearing that may have been fed northward into the STDS (Long et al., 2017). Studies in other areas of the Himalaya have estimated offset magnitudes of a similar order. In the Everest region, barometric gradients (Searle et al., 2002, 2003), telescoped isotherms (Law et al., 2011), and peak temperatures above and below the shear zone (Cottle et al., 2007, 2011) have yielded STDS offset ranges of ∼90–216 km, ∼25–170 km, and ∼15–67 km, respectively. In northwestern India, offset on the Zanskar shear zone, a western equivalent of the STDS, has been estimated between ∼25 and ∼120 km (Herren, 1987; Dèzes et al., 1999; Walker et al., 1999).

One interesting observation is that telescoped isotherms are located entirely within the gently shortened upper level of the STDS zone, and not within the highly thinned lower level, and thus appear to be spatially decoupled from strain magnitude. We suggest that isotherms were progressively elevated and compressed as a result of long-lived, distributed intrusion of granite sills during STDS shearing, which record crystallization spanning from ca. 24 to ca. 16 Ma (Kellett et al., 2009). A similar interpretation of significant heat contribution from granite sill injection during shearing has also been proposed for the STDS in the Everest region (Waters et al., 2018).

Evolution of the STDS in Northwestern Bhutan

The interpretations outlined above, when combined with timing constraints for metamorphism, granite intrusion, and shearing on the STDS from Kellett et al. (2009, 2010), define the following order of events (Fig. 12).

  • (1) 26–22 Ma: pre-STDS shearing (Fig. 12A): GH rocks that will eventually occupy the footwall and lower level of the STDS zone were buried to ∼9 kbar (∼35 km), and attained peak temperatures of ∼650–725 °C (Kellett et al., 2010; this study). Much of the burial of GH rocks can likely be attributed to construction of a south-vergent Paleogene thrust belt in overlying TH rocks, as documented in several places in southern Tibet (e.g., Ratschbacher et al., 1994; Wiesmayr and Grasemann, 2002; Murphy and Yin, 2003; Aikman et al., 2008). Prograde growth of monazite in GH rocks in northwestern Bhutan took place between ca. 26 and ca. 21 Ma, and peak metamorphism and partial melting are interpreted at ca. 22–23 Ma (Kellett et al., 2010). This brackets the earliest possible timing of shearing on the STDS. The Chekha Formation, which will eventually occupy the upper level of the STDS zone, attained peak temperatures of ∼300–350 °C, as indicated from the non-elevated “background” temperatures near the top of the section.

  • (2) 22–16 Ma: initial STDS shearing (Fig. 12B): Initial shearing was accompanied by retrograde growth of monazite in GH rocks from ca. 19 to ca. 15 Ma (Kellett et al., 2010), and intrusion and solid-state deformation of granite sills throughout the STDS zone, which crystallized between ca. 24 and ca. 16 Ma (Kellett et al., 2009). GH rocks in the lower level of the STDS zone accommodated high-magnitude, simple shear-dominated offset, which resulted in significant apparent thinning. Early, high-temperature quartz recrystallization microstructures (CBE, GBM) were generated while GH rocks were still >550 °C. In the upper level of the STDS zone, distributed top-to-NW offset was accommodated by an upward-decreasing gradient in pure shear-dominated lengthening. Protracted, spatially dense intrusion of granite sills (and also larger granite bodies, such as observed ∼20 km to the west of the studied transect; Fig. 1C) between ca. 24 and ca. 16 Ma enhanced temperatures throughout the STDS zone during shearing, eventually resulting in elevation and compression of isotherms into a relatively thin interval that lies entirely within the Chekha Formation.

  • (3) 16–>11 Ma: late-stage STDS shearing (Fig. 12C): Far along the exhumation path, quartz microstructures indicating SGR recrystallization were generated within the lower level of the STDS zone, while GH rocks were between ∼450 and 550 °C. A granite sill in the Chekha Formation cooled through ∼425 °C by ca. 11 Ma (Kellett et al., 2009), potentially bracketing the timing of development of SGR microstructures.

  • (4) 15–11 Ma: post-STDS thrusting (Fig. 12D): Motion on a top-to-S thrust fault that post-dated STDS shearing repeated a portion of the Chekha Formation section that contains the telescoped field gradient. The interpreted timing range is based on out-of-sequence thrusting documented across northern Bhutan (the Kakhtang and Laya thrusts), which is bracketed between ca. 15 and ca. 11 Ma (Grujic et al., 2002, 2011; Kellett et al., 2009; Warren et al., 2011).

Implications for Definition of the STDS Zone and Other Shear Zones

One interesting finding of our study is that the interval of telescoped isotherms lies entirely above the level of the STDS zone that experienced the highest offset magnitude and greatest apparent thinning (Figs. 11B, 12B, and 12C). This emphasizes the importance of defining shear zone boundaries on the basis of structural and/or kinematic criteria, and not exclusively on the spatial limits of metamorphic field gradients. This is in agreement with several studies that argue that defining the limits of a zone of high strain is the best way to delineate shear zone boundaries (e.g., Ramsay and Graham, 1970; Ramsay, 1980; White et al., 1980; Ramsay and Huber, 1983; Means, 1995; Searle et al., 2008).

As the STDS represents a top-to-N structure within an overall top-to-S orogen, it is tempting to use the limits of top-to-N shear-sense indicators to define its boundaries. However, we observed top-to-N indicators over a >10 km structural thickness (Fig. 10G), and top-to-N and top-to-S indicators are interspersed through an ∼11-km-thick section of GH and TH rocks in central Bhutan (Long et al., 2011b, 2017). Though we do map the base of the STDS zone at the onset of dominant top-to-N shearing, this is corroborated by an upward increase in fabric intensity that we interpret to represent a gradient in strain magnitude (e.g., Larson et al., 2017). We define the top of the STDS zone using a gradient in thin section-scale finite strain magnitude. Therefore, this case study demonstrates the utility of using a combination of quantitative fabric and finite strain analyses to delineate shear zone boundaries in situations where structural criteria indicative of strain gradients may not be easily identifiable in the field.

CONCLUSIONS

  1. We divide the STDS in NW Bhutan into an ∼2-km-thick lower level that accommodated ∼6–13 km of apparent structural thinning (86%–93% shortening) through high-magnitude (≥30–76 km), simple shear-dominant (Wm = 0.74–0.88), top-to-NW displacement, and an ∼3-km-thick upper level that accommodated ≥21 km of top-to-NW displacement via an upward decrease (from 44% to 2%) in pure shear-dominant (Wm = 0.00–0.45), transport-parallel lengthening.

  2. Peak temperatures within and beneath the lower level of the STDS zone are ∼650–750 °C. Two intervals of telescoped upright isotherms are observed in the upper level of the STDS zone, and define ∼160–260 °C/km field gradients. The lower interval is 1.3 km thick, and exhibits an upward decrease from ∼700 to ∼450 °C. Above this, temperatures increase abruptly from ∼450 to ∼620 °C, which is interpreted as the result of post-STDS thrust repletion. Moving structurally upward, temperatures decrease from ∼620 to ∼320 °C over a 1.1-km-thick interval, and temperatures decrease gradually from ∼320 to ∼250 °C through a 7-km-thick section of overlying TH rocks.

  3. Telescoped isotherms are located above the highly thinned lower level of the STDS zone, and therefore are spatially decoupled from strain magnitude. We interpret this to be a consequence of isotherm elevation during distributed granite intrusion that spanned the duration of STDS shearing. This work highlights the importance of defining shear zone boundaries on the basis of structural criteria, and not relying exclusively on the spatial limits of metamorphic field gradients. In addition, this study demonstrates the utility of using fabric and finite strain analyses to delineate shear zone boundaries when field criteria may not be available.

ACKNOWLEDGMENTS

We are indebted to Ugyen Wangda, Tashi Tenzin, and Yonten Phuntsho of the Bhutan Department of Geology and Mines. This manuscript was greatly improved by constructive reviews by Richard Law and Djordje Grujic. We would also like to acknowledge the pioneering work of Dawn Kellett on the South Tibetan detachment system, both along the studied transect and in other regions of Bhutan, which provides a foundation upon which our work builds. This work was funded by National Science Foundation EAR-1220300 awarded to S. Long and S. Gordon and by funds from Washington State University, Pullman, Washington, USA, awarded to S. Long.

1GSA Data Repository Item 2019160, Tables DR1–DR3, Figures DR1–DR7, and Discussions DR1–DR6, is available at http://www.geosociety.org/datarepository/2019, or on request from editing@geosociety.org.
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.