Foreland-directed propagation of high-grade tectonism in the deep roots of a Paleoproterozoic collisional orogen, SW Montana, USA

Deeply eroded convergent belts offer exceptional natural laboratories in which to expand our understanding of crustal processes during continent formation. One such example is the recently recognized late Paleoproterozoic Big Sky orogen (e.g., O’Neill, 1998; Harms et al., 2004a), expressed as high-grade deep crust exposed within the Laramide basement-cored uplifts in southwest Montana, United States. Although a significant late Paleoproterozoic thermal disturbance was recognized along the northern margin of the Wyoming Province half a century ago (Giletti, 1966), it is only within the past 15 yr that substantial evidence emerged indicating that a significant collisional orogen transects central and southwest Montana (O’Neill, 1998; Mueller et al., 2002; Roberts et al., 2002; Berger et al., 2004; Harms et al., 2004a). However, major issues remain unresolved, including details of the scale and tempo of convergence, and how the internal dynamics of the Big Sky orogen evolved. Specific questions include the following. What is the along- and across-strike extent of the orogenic belt as defined by thermal and kinematic patterns, and what are the extents of the metamorphic core and foreland components? Can the timing of precollision convergence, initiation of collision, and duration of collision and postcollision relaxation processes be distinguished? These specific questions bear directly on larger questions regarding the assembly of Laurentia, including whether the rocks in southwest Montana represent one collision event or two (e.g., Mueller et al., 2005), or if there was any collision at all (Boerner et al., 1998).

New data are presented here from field and structural analysis, petrology, and geochronology in the Northern Madison Range, a key locality for clarifying the nature of the hinterland-foreland transition of the Big Sky orogen. Coupling high-spatial-resolution in situ geochronologic analysis of zircon and monazite with structural and petrologic observations allows recognition of links between accessory mineral growth, major phase reactions, and deformation textures. Results from this work show that Paleoproterozoic crustal depths of ≥25–30 km are exposed in the central portion of the Northern Madison Range, which considerably expands the across-strike extent of the orogen’s metamorphic core. Moreover, the age of thermotectonism in the study area appears younger than in other parts of the region affected by the Big Sky orogeny, suggesting propagation of peak metamorphism and deformation from NW to SE over 40–80 m.y. These results are similar to lateral growth and propagation patterns observed in other major collisional orogens (e.g., Jamieson et al., 2011; Staples et al., 2014).

The Wyoming Province is the southernmost component of the Archean core of North America and is bounded to the south, east, and northwest by Paleoproterozoic (ca. 1.8–1.7 Ga) orogenic belts (Cheyenne belt, Trans-Hudson orogen, and Paleoproterozoic component of the Great Falls tectonic zone, respectively; Gorman and Clowes, 2002; Mueller and Frost, 2006; O’Neill and Lopez, 1985). The northwestern part of the Wyoming Province is commonly referred to as the Montana metasedimentary province, one of several domains that are distinguished based on lithology and general age (Mogk et al., 1992; Mueller et al., 2005; Foster et al., 2006). The Montana metasedimentary province includes most Precambrian crystalline exposures in southwest Montana and is approximately bounded to the southeast by Late Archean shear zones in the northwestern Beartooth Mountains (Fig. 1A). The province is characterized by a relative abundance of supracrustal metamorphic rocks, including pelitic schist, quartzite, and carbonate rocks along with ca. 3.3–3.0 Ga quartzofeldspathic gneiss and other granitoids (Mogk et al., 1992; Mueller et al., 1996). To the southeast, the Beartooth-Bighorn magmatic zone includes a larger proportion of 2.8 Ga orthogneiss and associated trondhjemite-tonalite-granodiorite suites (Mueller and Frost, 2006). The Precambrian rocks in southwest Montana have been tectonically reworked multiple times and metamorphosed to granulite and upper-amphibolite facies (e.g., Giletti, 1966; Spencer and Kozak, 1975; James and Hedge, 1980; Mogk, 1992; Harms et al., 2004a, 2004b). The dominant structural grain across the region strikes NE-SW, defined by major foliation and fold axial surfaces (e.g., Spencer and Kozak, 1975; Harms et al., 2004b), shear zones (e.g., Erslev and Sutter, 1990; Kellogg and Mogk, 2009; Johnson et al., 2014), and the overall trend of a thermochronologic feature called Giletti’s line (defined in next paragraph; Giletti, 1966) and the Great Falls tectonic zone.

Early isotopic studies in the Montana metasedimentary province identified largely Archean protolith ages, but they also recognized a regional late Paleoproterozoic thermal overprint (Giletti, 1966, 1971; Mueller and Cordua, 1976; James and Hedge, 1980; Weyand, 1989; O’Neill et al., 1988; Mueller et al., 1993). Biotite and muscovite K-Ar and Rb-Sr dates are ca. 1.8–1.6 Ga northwest of a NE-trending transitional boundary that was originally drawn through the Southern and Northern Madison Ranges (Fig. 1B; Giletti, 1966, 1971; Reid et al., 1975; Mueller and Cordua, 1976). This feature, commonly referred to as “Giletti’s line,” was interpreted to represent the transition from warmer temperatures to the northwest, which induced complete resetting or new mica growth at 1.8–1.6 Ga, to cooler temperatures to the southeast, which were not high enough in Proterozoic time to significantly disturb the Ar system (Giletti, 1966). Subsequent work by Giletti (1971) revealed excess Ar in several biotite samples from the Gallatin Canyon localities, which called into question the exact position and nature of the thermal transition in the Northern Madison Range.

The spatial extent and magnitude of this late Paleoproterozoic thermotectonism remain unclear. Previous workers linked the Proterozoic thermal signature to relatively static low- to moderate-temperature overprinting, with associated localized greenschist-facies deformation, and assigned an Archean age to the regional upper-amphibolite- to granulite-facies tectonism (Mueller and Cordua, 1976; James and Hedge, 1980; Salt, 1987; Weyand, 1989; Erslev and Sutter, 1990; Mogk et al., 1992). More recent thermochronological (Harlan et al., 1996; Roberts et al., 2002; Brady et al., 2004a; Hames and Harms, 2013) and geochronological studies (Roberts et al., 2002; Mueller et al., 2004, 2005; Cheney et al., 2004; Ault et al., 2012; Alcock et al., 2013) documented high-grade metamorphism (up to granulite facies and 1.2 GPa) and deformation between ca. 1.8 and 1.7 Ga throughout the western ranges in the Montana metasedimentary province (Fig. 1). Harms et al. (2004a) called these tectonic events the Big Sky orogeny and suggested that they were caused by closure of an ocean basin and subsequent collision of the Archean Medicine Hat block with the Wyoming Province.

The extent of this overprint in the eastern part of the Montana metasedimentary province (Gravelly, Northern, and Southern Madison Ranges; Fig. 1B) is less certain. Much of the Northern Madison Range is located between the known Big Sky high-grade tectonism described here and a likely similarly aged, but localized Paleoproterozoic greenschist-facies shear zone, the Madison mylonite zone of Erslev and Sutter (1990), in the southern half of the Southern Madison Range (Fig. 1B). The low metamorphic grade of this shear zone and the apparent lack of warm Big Sky–aged temperatures (<300–400 °C) in much of the surrounding Southern Madison Range and Gravelly Range suggest that these ranges represent the Big Sky foreland. Consequently, the Northern Madison Range is an important area in which to further clarify the nature and location of the hinterland-foreland transition of the Big Sky orogen.

The Northern Madison Range (Fig. 2) is an ∼50-km-long, NW-SE–trending, Laramide-age, basement-cored uplift bounded to the south by the steeply dipping Spanish Peaks thrust fault (Garihan et al., 1983), which was last active during the Late Cretaceous (Kellogg and Harlan, 2007). The core of the range is predominantly Archean crystalline rock that has been multiply deformed and metamorphosed (Spencer and Kozak, 1975; Kellogg and Mogk, 2009). Recent work indicates that at least the northwest portion of the range shows evidence for Paleoproterozoic Big Sky tectonometamorphism (Ault et al., 2012).

The range consists of quartzofeldspathic orthogneiss, foliated granitoids, mafic amphibolite, supracrustal schist, and quartzite. Previous workers distinguished two separate terranes based on variations in lithologies and metamorphic assemblages (Salt, 1987; Weyand, 1989; Mogk et al., 1992). The Jerome Rocks Lake terrane contains predominantly trondhjemitic gneiss, sillimanite-bearing metapelite and migmatite, and it is separated by the Mirror Lake shear zone from orthogneiss with intercalated kyanite- and sillimanite-bearing metapelitic schist, quartzite, ultramafic schist, and amphibolite of the Gallatin Peak terrane to the south (Salt, 1987;Fig. 2). Salt (1987) and Mogk et al. (1992) interpreted the Gallatin Peak terrane as representing a ca. 3.2 Ga subduction-related calc-alkaline continental magmatic arc. U-Pb zircon dates from orthogneiss units in both terranes yielded a range of ca. 3.3–2.7 Ga crystallization ages (Weyand, 1989).

The dominant structural trend is a NE-SW–striking foliation that varies considerably in dip magnitude and direction, and that ranges in intensity, with several high-strain zones that are locally mylonitic. Spencer and Kozak (1975) interpreted two main stages of deformation in the range: (1) early gneissic foliation development and isoclinal folding of the foliation, and (2) subsequent ductile deformation and the development of broad folds. Discrete high-strain zones include the Mirror Lake and Big Brother shear zones (Salt, 1987; Mogk and Henry, 1988; Mogk et al., 1992), Crooked Creek mylonite (Kellogg and Mogk, 2009), the Hellroaring Creek shear zone, and the Spanish Creek mylonite (Johnson et al., 2014; Fig. 2). These ductile shear zones generally strike NE-SW and dip moderately to steeply to the SE or NW (Fig. 2). The timing of this deformation is not well constrained. Previous studies suggested a Late Archean age for the Mirror Lake shear zone and Crooked Creek mylonite (Weyand, 1989; Mogk et al., 1992; Kellogg and Mogk, 2009). However, the prevalence and proximity of Proterozoic high-grade metamorphism and deformation in nearby ranges to the south and west (e.g., Erslev and Sutter, 1990; Harms et al., 2004b) coupled with new data from this study suggest that some of the deformation is likely Paleoproterozoic.

The metamorphic grade observed within the Northern Madison Range varies from granulite to amphibolite facies (Spencer and Kozak, 1975; Salt, 1987; Mogk and Henry, 1988; Ault et al., 2012), with some greenschist-facies overprinting in the southeastern portion of the range in Gallatin Canyon (Condit et al., 2012). Previous thermobarometric estimates indicate relatively high pressure and temperature conditions of up to 1.0–0.7 GPa and up to 700 °C (Salt, 1987; Mogk and Henry, 1988; Kellogg and Mogk, 2009). The age of high-grade metamorphism was interpreted as Mesoarchean to Neoarchean (Mogk et al., 1992; Kellogg and Mogk, 2009; James and Hedge, 1980; Salt, 1987; Mueller and Frost, 2006). This interpretation is based on zircon crystallization ages for granitoids in the central part the Spanish Peaks area, interpreted to have been emplaced synchronously with metamorphism and deformation (Weyand, 1989). However, Ault et al. (2012) reported a 1753 ± 18 Ma metamorphic zircon age (ion probe, 2σ error) from a mafic dike that was metamorphosed and internally deformed at 1.2 GPa and 800 °C in the northwestern portion of the range (Fig. 2).

This study is focused in Bear Basin and its surrounding areas, which occupy the northwest portion of the Gallatin Peak terrane (Figs. 2 and 3). Previous detailed work in this region is limited to theses by Salt (1987) and Weyand (1989) and a recently published map at 1:24,000 scale by Vuke (2013). Salt (1987) mapped the major lithologies in the Bear Basin region, and Weyand acquired multigrain U-Pb thermal ionization mass spectrometry (TIMS) zircon data from these units. For this study, the Bear Basin region was remapped at 1:15,000 with a focus on structural geology and key units for relative and absolute timing of deformation and metamorphism (Fig. 3).

The oldest known rocks in the Bear Basin region are a series of ca. 3.2–3.0 Ga granitoid gneisses. The first major unit is a heterogeneous tonalitic gneiss with a zircon U-Pb crystallization age of 3244 ± 19 Ma (Weyand, 1989) that occurs in the vicinity of Gallatin Peak and south of the Mirror Lake shear zone. The unit consists of centimeter-scale layers of tonalite, diorite, and amphibolite, with variable proportions of Pl + Qz + Bt + Hbl (all mineral abbreviations after Whitney and Evans, 2010), and minor granitoid. South of the tonalitic gneiss exposures, a hornblende monzodiorite crops out (Fig. 3) that is similar to the tonalite in major mineralogy but is texturally more isotropic (non-layered), contains abundant millimeter-scale epidote, and yields a zircon U-Pb crystallization age of 3195 ± 43 Ma (Weyand, 1989). The third major unit is a porphyritic granodiorite exposed in the southern Bear Basin area (Fig. 3). It contains centimeter-scale relict K-feldspar phenocrysts and has a zircon U-Pb crystallization age of 3177 ± 36 Ma (Weyand, 1989). Salt (1987) documented inclusions of the previously described igneous units in this granodiorite and thus interpreted that the granodiorite is the youngest in the suite. Collectively, these units are the basis of Mogk et al.’s (1992) interpretation that the Gallatin Peak terrane represents part of a ca. 3.2 Ga calc-alkaline continental magmatic arc.

The fourth major unit is composed of intimately mixed and undifferentiated migmatitic biotite schist, granite, and granitic gneiss (Fig. 3, biotite gneiss unit). In the map area, this unit lacks sedimentary structures, and no aluminous metamorphic index minerals were observed other than local garnet. Weyand (1989) interpreted the protolith as igneous and reported a zircon U-Pb date of 2868 ± 24 Ma.

The youngest reported granitoid units within the Gallatin Peak terrane are small-volume foliated to unfoliated granites. Those we observed were 1–30-m-thick tabular sheets of Ms + Bt ± Grt leucogranite that are concordant with the foliated host rock. These leucogranites appear to be more abundant on the northern margin of the Hellroaring Creek shear zone (Fig. 2). Weyand (1989) reported a 2680 ± 130 Ma U-Pb zircon date for a weakly foliated granitic dike (Bt + Hbl + Ep + Pl) on the northern margin of the Gallatin Peak terrane.

Supracrustal rocks within the Gallatin Peak terrane include an intercalated package of aluminous schist, amphibolite, and minor quartzite. This unit crops out in central Bear Basin and in the adjacent upper drainage of Hellroaring Creek (Fig. 3). It is collectively referred to here as the Bear Basin schist. The schist contains a wide range of compositions, from aluminum-rich layers with varying amounts of kyanite (Fig. 4A), cordierite, sillimanite, staurolite, and garnet to less aluminous orthoamphibole-rich layers. Concordant garnet-hornblende amphibolite layers are interpreted as metavolcanics.

Salt (1987) interpreted a coarse Ged + Ky layer close to a contact with the Hbl-monzodiorite described in Bear Basin as a contact aureole, thus placing a minimum depositional age on the Bear Basin schist at ca. 3.2–3.1 Ga. However, the northwest contact between these two units is locally conglomeratic, with quartzite lenses and pebble-sized quartz and granitoid clasts (Fig. 4B) observed in at least two separate localities (Fig. 3). These observations suggest a depositional contact rather than an intrusive one, implying a maximum depositional age of ca. 3.1 Ga. Weyand (1989) also interpreted this unit as sedimentary in origin, based on whole-rock compositions, and implied that analyzed zircon, with ²⁰⁷Pb/²⁰⁶Pb dates ranging from 3025 to 3160 Ma, are detrital. This package is similar in lithology and appearance to the Spuhler Peak metamorphic suite found within the Tobacco Root Mountains, which has been interpreted to be a shallow-marine sequence of volcanics and intercalated sediments metamorphosed in a single event during the Big Sky orogeny (Burger et al., 2004).

Amphibolite dikes (Hbl + Pl ± Grt ± Cpx) appear to be some of the youngest intrusive units in the region (Fig. 5A), and they have similar deformational and crosscutting field relationships to the mafic dikes and sills observed in the Tobacco Root Mountains (Brady et al., 2004b) and in Bear Trap Canyon (Ault et al., 2012). The contacts of internally deformed dikes locally truncate early gneissic layering in tonalite units (Figs. 5B–5C). Although concordant amphibolite layers occur, no dikes with crosscutting contacts were observed in the Bear Basin schist, an observation similar to the Spuhler Peak metamorphic suite. Several undeformed postkinematic dikes thought to be late Proterozoic in age (Harlan et al., 1996) occur in the range, although none has been observed in the study area (Fig. 3).

At least three phases of deformation affected the rocks of the Gallatin Peak terrane (Table 1). The earliest structures likely represent more than one phase of deformation. These include early gneissic layering in the tonalite, early penetrative foliation in the Bear Basin schist, magmatic layering defined by variable concentrations of K-feldspar phenocrysts in the porphyritic granodiorite, and locally preserved inclusion trails within porphyroblastic garnet in deformed mafic dikes. Our observations indicate that these earliest fabrics cannot all be the same age. For instance, the gneissic layering in the tonalite is locally truncated by intrusive contacts of the deformed and metamorphosed mafic dikes (Figs. 5B–5C). However, for the purposes of this paper, all of these structures are collectively referred to as D₁, and all are variably transposed by D₂ structures (e.g., Fig. 4C). The early structures are locally preserved in low-strain regions such as the hinge of a map-scale F₂ fold in Bear Basin or in other low-strain zones in the upper drainage of Hellroaring Creek (Fig. 3). There, S₁ and F₂ enveloping surfaces tend to strike N or NW with moderate to steep dips, and they are distributed in a manner that is consistent with variable reorientation by later F₂ and probably F₃ folding (Fig. 6A).

D₂ deformation is the most pervasive in the range, and D₂ structures are ubiquitous within the rocks of the northern portion of the Gallatin Peak terrane. The penetrative S₂ foliation is NE striking and moderately SE dipping, with a mean orientation of 044°, 57° (Fig. 6B). S₂ surfaces contain a mineral lineation (L₂) defined by biotite, hornblende, or kyanite and/or a stretching lineation defined primarily by elongate feldspar grains or locally by stretched pebbles in conglomeratic horizons of the Bear Basin schist. L₂ has a mean orientation that is approximately downdip (60°→129°; Fig. 6B). F₂ folds are common and recognizable from outcrop to kilometer scales. These folds are tight to isoclinal (Fig. 4C), with axial surfaces subparallel to S₂ and with hinge lines that are subparallel to L₂ (Fig. 6C). The one exception is in a relatively low D₂ strain zone in the upper drainage of Hellroaring Creek, where F₂ hinge lines consistently plunge N to NE, possibly reflecting regional variability in the pre-D₂ orientation of the folded S₁ surfaces or that D₂ strain was not sufficiently high enough in this area to rotate the F₂ hinge lines to the direction of maximum stretching. The most prominent map-scale D₂ structure is a kilometer-scale, moderately inclined, and moderately southeast-plunging synformal isoclinal fold of the Bear Basin schist (Fig. 3). The hinge region of this fold is exposed in central Bear Basin (Fig. 4C).

The youngest deformation phase observed in the region (D₃) is limited to steeply inclined and moderately plunging, open to tight folds of the S₂ fabric (Fig. 4D). The mean axial surface of F₃ folds (001°, 68°) strikes more northerly than that of F₂ folds, but hinge line orientations are similar (mean trend and plunge 54°®152°). F₃ structures commonly have NE-vergent Z-asymmetry (Fig. 6C).

The Gallatin Peak terrane is bounded to the northwest by the Mirror Lake shear zone, described by Salt (1987) as an ∼500-m-wide, NE-striking, and moderately SE-dipping shear zone. We have not observed this structure, and it is not exposed in the study area. A second high-strain zone, here named the Hellroaring Creek shear zone, forms the approximate southeastern boundary of the study area (Fig. 2). This structure is a 1.5–2-km-wide zone of strongly foliated to mylonitic rocks that strikes southwest and dips steeply to the northwest (mean orientation 228°, 79°; Fig. 6B). Shear sense indicators record components of SE-vergent thrust displacement and right-lateral shearing, suggesting a transpressional environment. The dextral and reverse asymmetry of F₃ folds described here suggests that these folds may be a coeval and lower-strain expression of the Hellroaring Creek shear zone.

D₁ and D₂ structures described here both belong to the “first orogeny” of the more regional study by Spencer and Kozak (1975). Their “second orogeny” was based on the observation that the dominant foliation across the range varies in dip direction from SE to NW, and thus large-scale folds without an associated axial planar foliation were inferred. In the present study area, similar variation in the orientation of dominant penetrative fabrics is explained by overprinting of the S₂ foliation by the Hellroaring Creek shear zone in the eastern portion of the study area, which is tentatively related to D₃. Salt (1987) recognized two sets of fold styles in this region that appear to correlate with F₂ and F₃ as described here. It was concluded from that study that both fold sets formed at relatively high temperature, and thus probably developed close in time during one orogenic cycle.

Two samples of Bear Basin schist (GP7c and GP7f) and one mafic dike sample (AA09-61) were targeted to help constrain conditions and timing of metamorphism and deformation in the Gallatin Peak terrane. The Bear Basin schist was sampled in a pelitic horizon because the mineral assemblages are amenable to quantitative thermobarometry and monazite geochronology. The mafic dikes were targeted because they represent the youngest metamorphosed and deformed intrusive units in the study area. Pressure and temperature (P-T) equilibrium conditions calculated from all three samples constrain representative metamorphic conditions in the northwest portion of the Gallatin Peak terrane. Geochronology included electron microprobe (EMP) U-Th–total Pb analysis of monazite from the two Bear Basin schist samples and in situ secondary ion mass spectrometry (SIMS) U-Pb zircon dating of the mafic dike sample.

The two Bear Basin schist samples were collected near the hinge region of the kilometer-scale F₂ fold in Bear Basin (Fig. 3). Sample GP7c is from an ∼1–2-m-thick layer and contains Grt + Ky + Sil + St + Bt + Pl + Qz, with accessory zircon, monazite, and xenotime. The sample contains a dominant S₂ fabric, which is defined by (1) shape-preferred orientations of biotite, staurolite, kyanite, and sillimanite, and (2) alternating millimeter- to centimeter-scale layers of Pl + Bt + Qz and more aluminous layers. Sample GP7f, located several meters away, exhibits an S₁ fabric only partially transposed by F₂ folds. The sample is from a less-aluminous component of the schist, containing Grt + Bt + Pl + Qz and accessory monazite, zircon, and rutile. The GP7f thin section was cut perpendicular to the F₂ hinge line, so that traces of S₁ and S₂, both defined primarily by aligned biotite, are visible.

Garnet porphyroblasts in both Bear Basin schist samples are 0.5–3 mm in diameter, are subhedral to anhedral, and contain abundant inclusions of biotite, quartz, monazite, zircon, and fracture-filling late chlorite (Figs. 7A–7B). Kyanite is also locally included in garnet in sample GP7c (Fig. 7A), suggesting that it was the first aluminosilicate to stabilize. Both kyanite and sillimanite occur in the matrix as fabric-defining phases, with sillimanite texturally appearing to form after kyanite, (Fig. 7C). Although staurolite occurs in the matrix of sample GP7c (Fig. 7D), it is locally restricted as inclusions in garnet in other samples from the same outcrop (Fig. 7E), implying growth as an early prograde phase. As described in detail here, monazite occurs throughout the matrix in both schist samples and commonly as inclusions in several of the major mineral phases (Figs. 7A–7E). Xenotime was observed as an inclusion in staurolite (Fig. 7E), but not in the matrix in GP7c, and the phase was not observed in GP7f.

Sample AA09-61 was collected from an amphibolite dike that crops out on the flank of Gallatin Peak, ∼1 km north of Bear Basin (Fig. 3). Dike contacts are locally discordant to an early gneissic layering in the tonalite host (S₁ gneiss, Figs. 5B–5C). The dike contains a dominant penetrative S₂ fabric, defined by aligned hornblende and plagioclase, that appears axial planar to a meter-scale F₂ fold of the dike and host gneiss (Fig. 5A). The sample contains Hbl + Pl + Qz + Ilm ± Bt with ∼15–10-mm-diameter anhedral to subhedral garnet porphyroblasts (Figs. 5D and 7F). Garnet porphyroblasts contain abundant inclusions of hornblende, plagioclase, quartz, biotite, and ilmenite that define an earlier fabric oblique to the penetrative foliation within the dike. Garnet porphyroblasts are mantled by inner plagioclase and outer hornblende haloes. Plagioclase in the inner halo is evenly distributed around the garnet, whereas the outer hornblende haloes are commonly elongate in the S₂ fabric. Euhedral hornblende locally occurs within the plagioclase halo.

Elemental X-ray maps for Ca, Al, Mg, and Mn and mineral compositions were analyzed using a JEOL 8600 electron microprobe at the University of Colorado–Boulder. Beam conditions of 15 kV voltage, 50 nA current, 100 μs dwell time, and step size of ∼10–25 μm were used for X-ray maps. Quantitative mineral compositions were acquired with 15 kV voltage, 20 nA current, and a focused beam (∼1 μm diameter) for garnet and staurolite, and defocused beam (5–10 μm diameter) for micas, amphiboles, and feldspars, and count times ranging from 20 to 40 s. Select mineral compositions are reported in Table 2, and those used in thermobarometry calculations are noted.

Garnet compositions in both samples are dominantly almandine (X_Alm = 0.62–0.66 in GP7c; 0.72–0.76 in GP7f). Magnesium number (Mg#) is broadly homogeneous in the cores (0.30 in GP7c; 0.24 in GP7f) and sharply drops near the rims (0.24 in GP7c; 0.20 in GP7f). Spessartine content is also homogeneous in grain interiors, with a slight increase at the rim (from 0.09 to 0.11 in GP7c; 0.01 to 0.04 in GP7f), especially in smaller, isolated garnet islands (Figs. 8A–8B). Grossular garnet composition in sample GP7c is very low (X_Grs = 0.01), and the garnet is unzoned with respect to Ca. Calcium in garnet in sample GP7f is slightly zoned, with the core recording X_Grs = 0.04, a drop to 0.03 in a distinct moat, and an increase to 0.04 along the rims (Fig. 8B). These garnet compositional patterns suggest homogenization at high temperature with late modifications associated with limited resorption (X_Sps increase at margins) and diffusional Mg-Fe exchange (Tracy et al., 1976; Kohn and Spear, 2000).

Biotite is unzoned in both samples (Mg# = 0.60–0.62 in GP7c; 0.53–0.55 in GP7f). Staurolite in GP7c has Mg# = 0.24–0.25 and lacks zoning. Plagioclase in GP7f has minor zoning, with a higher-Ca core (An₂₂) that is surrounded by a volumetrically dominant low-Ca (An₁₉) moat, which in turn is surrounded by thin high-Ca rims with An₂₀ (Fig. 8B). Plagioclase variations in GP7c show a similar trend to GP7f, with slightly higher-Ca cores (An₀₉) and slightly lower-Ca rims (An₀₈).

Garnet in the mafic dike preserves little chemical variation (Fig. 8C). Grossular content varies from 0.27 to 0.29, with the highest X_Grs in the core. The Mg# is 0.11–0.12, with the exception of the rim, where it drops to 0.08. Spessartine content also increases at the rim of the garnet, with an average X_Sps of 0.07, rising to 0.11 within 10 μm of the grain edge.

Matrix hornblende and plagioclase have little zoning and are chemically homogeneous across much of the sample. The Mg# values of hornblende inclusions in garnet and in the matrix are 0.33–0.34 and 0.34–0.36, respectively. Hornblende grains in the haloes around garnet have similar compositions. Plagioclase compositions show minimal variation from An₃₂ to An₃₄, with plagioclase inclusions in garnet having the higher Ca content. The outer hornblende halo is interpreted to be a prograde feature developed by consumption of plagioclase during synkinematic garnet growth with respect to S₂ development. The inner plagioclase halo is interpreted to have formed during retrogression and partial garnet resorption, which is also supported by the distinct increases in Mn content at the extreme margins of the garnet porphyroblasts.

Thermobarometric calculations and phase assemblage modeling (pseudosection) were used to evaluate equilibrium peak P-T conditions and the relative timing of growth of the major metamorphic phases. X-ray maps of targeted garnet and matrix areas were made to evaluate mineral zoning and to ensure that mineral compositions most closely representing peak equilibrium conditions were chosen for P-T calculations (e.g., Fig. 8). First, calculations were made using the program TWQ 2.34 (Berman, 1991, http://gsc.nrcan.gc.ca/sw/twq_e.php) and the internally consistent database of Berman and Aranovich (1996; updated in 2007; Berman et al., 2007). Only well-calibrated reactions were used involving garnet, biotite, sillimanite, plagioclase, and quartz for the schist and garnet, and plagioclase, quartz, and hornblende for the mafic dike. Estimated uncertainties associated with these calculations are ±0.1 GPa and ±50 °C (Berman, 1991). All Fe was assumed to be Fe²⁺. This assumption likely has the greatest effect on the garnet-biotite temperatures, where the results may be overestimated by a few tens of degrees Celsius (Schumacher, 1991). Equilibrium P-T conditions for the three samples analyzed were 0.90 GPa, 730 °C for GP7c (Fig. 9A); 0.87 GPa, 700 °C for GP7f (Fig. 9A); and 0.92 GPa, 700 °C for AA09-61 (Fig. 9B). For the Bear Basin schist samples (GP7c and GP7f), duplicate TWQ 2.34 calculations where Fe³⁺ was considered via charge balance for biotite and garnet yielded temperatures ∼30 °C lower. Temperature estimates from AA09-61 showed no difference when considering Fe³⁺. Results from all three samples (∼0.9 GPa, 700–730 °C) are indistinguishable from one another within the estimated range of uncertainty, but the lower end of this temperature range is probably more accurate due to the possible presence of ferric iron in biotite.

Phase assemblage modeling of Bear Basin schist was performed to better understand the relative timing of growth of major mineral phases, the peak equilibrium metamorphic conditions, and to specifically evaluate the likelihood of metastable phases (e.g., staurolite) at higher P-T conditions. Modeling of sample GP7c (bulk composition: SiO₂ 71.32%, Al₂O₃ 14.48%, CaO 0.70%, FeO 4.44%, MgO 2.71%, K₂O 1.69%) in the CaKFMASH system was conducted using PerpleX 6.6.6 (Connolly and Petrini, 2002) with the updated thermodynamic database of Holland and Powell (2011). The results of this modeling indicate maximum thermal stability of staurolite in the range of 670–660 °C over pressures between 0.9 and 0.7 GPa, and stability for the interpreted peak assemblage at conditions consistent with those found using thermobarometric calculations (Fig. 9A).

High-spatial-resolution geochronological techniques are required to analyze accessory minerals because the zircon and monazite grains are small and/or compositionally zoned, with short dimensions of grains or compositional domains commonly ≤20 µm. Monazite and zircon were identified in thin section by automated scanning electron microscope analysis (QEMSCAN) at the Advanced Mineralogy Research Center at the Colorado School of Mines, Golden, Colorado. QEMSCAN instrument parameters are detailed in Ault et al. (2012).

Backscatter electron (BSE) images, acquired using a JEOL 8600 electron microprobe at the University of Colorado–Boulder, provided context and textures of select grains. All subsequent mapping and analyses were conducted at the University of Massachusetts–Amherst. X-ray maps for U Mβ, Th Mα, Y Lα, and Ca Kα, and in some cases for Nd Lα, were generated for each monazite grain of interest using a Cameca SX-50 electron microprobe. These maps documented compositional domains and identified areas of interest for quantitative analysis.

U-Th–total Pb monazite dates were acquired using the modified Cameca SX-100 (Ultrachron). Detailed analytical procedures and standard compositions used, including count times, standards, and a list of spectrometers, followed those presented in Appendix A of Dumond et al. (2008), Williams et al. (2006), and Jercinovic et al. (2008). Background collection and subtraction followed the multipoint background method of Allaz et al. (2011). For each chemically homogeneous domain, a weighted average of 3–8 analyses was calculated and reported with a 2σ uncertainty. The reported uncertainty represents the larger of either that calculated by propagating both analytical uncertainty on trace-element compositions through the age equation plus an estimated 1% uncertainty on background intensities (Williams et al., 2006), or two times the standard error of the mean (Table 3). By comparison with X-ray maps and evaluation of the analytical compositions, points that inadvertently sampled a compositional boundary or fell outside of the targeted domain were not included in the weighted mean calculation. The monazite consistency standard used is the Moacyr Brazilian pegmatite monazite with weighted mean isotope dilution (ID) TIMS dates of 506.4 ± 1.0 Ma (2σ, mean square of weighted deviates [MSWD] = 0.6) for ²⁰⁸Pb/²³²Th, 506.7 ± 0.8 Ma (MSWD = 0.83) for ²⁰⁷Pb/²³⁵U, and 515.2 ± 0.6 Ma (MSWD = 0.36) for ²⁰⁶Pb/²³⁸U (W.J. Davis [Geological Survey of Canada], 2007, personal commun.).

Monazite grains within the two Bear Basin schist samples, GP7c and GP7f, range in long-axis dimension from ∼20 to 140 μm and have shapes ranging from elongate to equant. Grains occur throughout the quartzofeldspathic matrix and are commonly included in, or in close proximity to, all major metamorphic phases (Grt, Ky, St, Sil, and Bt; Figs. 7A–7E). Three distinct monazite domain populations occur in the two samples, based primarily on composition, as well as supporting evidence from grain textures and zoning patterns. Therefore, some grains may contain more than one population. Population 1 consists of low-U core domains (≤0.5 wt% U; Table 3; Figs. 10A–10C). Population 2 consists of domains defining whole grains, inner rims (Fig. 10C), or outer rims with >1 wt% U (Fig. 10B). Both of these populations occur in grains located in the matrix, as inclusions in garnet in both samples, and as inclusions in kyanite and staurolite in sample GP7c. Grains with rim domains of population 2 are also locally intergrown and included within foliation-defining sillimanite in GP7c (Fig. 7C). Monazite with population 3 domains only occurs in the less-aluminous schist sample GP7f, where grains locally exhibit thin (∼5–10 μm), relatively high-Y rims (Figs. 10C–10D). These population 3 rim domains all contain 4–10 times the Y content of the core domains within the same grain. Grain GP7f-m3 contains all three populations (Fig. 10C). Both samples contain elongate monazite grains with long axes subparallel to a local S₁ or S₂ fabric (Figs. 10C and 10D). For example, grain GP7f-m3 contains a population 3 domain that forms an asymmetric tip aligned with the S₂ fabric (Fig. 10C).

Thirty-four grains were investigated through BSE and X-ray maps, and quantitative data were collected from 19 domains in 14 grains from the two schist samples (seven grains from each sample; Table 3; Fig. 11). This selection of grains represents the full range of textural settings and morphologies described earlier. Two population 1 domains from monazite GP7c-m12, a matrix grain occurring at the contact between biotite and quartz, yield Mesoarchean dates of 2995 ± 145 Ma and 3255 ± 14 Ma (all errors 2σ; Fig. 10A). All other population 1 domains, which represent the cores of several matrix grains, as well as the cores of two grains included in garnet (one from each sample), yield dates between 1752 ± 10 Ma and 1721 ± 15 Ma. Thus, population 1 was subsequently split into Archean population 1a and Paleoproterozoic population 1b, but no distinguishing characteristics other than age are apparent. Population 2 dates range from 1748 ± 8 Ma to 1720 ± 5 Ma, and they are mostly from rim domains of matrix grains, including two intergrown with sillimanite, but also include one whole grain and one rim domain included in garnet (Table 3). There is no apparent statistical distinction of dates from grains of different orientation (Figs. 10C–10D). The youngest dates are from population 3 rim domains in GP7f, ranging from 1706 ± 11 Ma to 1704 ± 10 Ma.

BSE and cathodoluminescence (CL) images of selected zircon targets were used to assess zoning and reveal internal grain structure. BSE images were acquired at the University of Colorado–Boulder, and CL images were taken at the University of Wyoming, Laramie, Wyoming. The thin section was trimmed into ≤7-mm-diameter chips containing the target zircon. The chips were mounted in epoxy, polished, cleaned, and coated with a nanometer-thick conductive coating of Au as detailed in Schmitt et al. (2010). Mounts also contained zircon standard AS3 (1099 ± 1 Ma; Paces and Miller, 1993). The U standard was NIST 91500 (Wiedenbeck et al., 1995).

Quantitative isotopic analyses were done on the high-sensitivity Cameca ims 1270 ion microprobe at the University of California–Los Angeles. The field aperture of the ims 1270 was adjusted to subsample secondary ions emitted from the interior of the primary beam pit in order to analyze grains with minimum dimensions as small as 4 μm at radiogenic yields typically >95% for radiogenic ²⁰⁶Pb (²⁰⁶Pb*). The reduction in the field aperture suppressed common Pb, and the surface was flooded with oxygen gas to increase the Pb* yield. The instrument operating conditions and procedures performed to optimize secondary ion yields for the present study were described in detail by Grove et al. (2003) and Schmitt et al. (2010). Zircon standard AS3 was measured at the beginning and end of each analytical session, during which a two-point linear drift correct was applied to the data. Th was not measured.

Of the 43 zircon grains evaluated in detail (optically and through BSE imaging), two zircon populations were identified in this sample based on morphology and texture. Although CL images were acquired, the internal structure of the grains is more clearly visible in BSE images. Zircon grains from AA09-61 range from 17 to 60 μm along their long axis. Similar to the monazite, zircon population distinctions are specific to domains, of which more than one commonly occur in a single zircon grain. Population 1 domains are distinctly inclusion rich and commonly contain small patchy variations in brightness in BSE images. Some zircons consist entirely of population 1 domains (Fig. 12A), but most contain population 1 domains as cores (Figs. 12B–12D). Population 2 domains are inclusion free and exhibit little to no zoning. They occur as rims of variable thickness around inclusion-rich population 1 domains (e.g., Figs. 12C and 12D) and as whole grains (Figs. 12E–12F). Texturally, zircons with population 1 domains occur as inclusions in garnet, in the hornblende and plagioclase haloes around garnet, and in the matrix. Zircon containing population 2 domains occurs in all of the same settings, but is more common in the matrix.

Twenty-four grains were selected for U-Pb analysis representing the full range of textural settings and grain sizes (Table 4). Of the 24 grains dated, eight grains were associated with garnet porphyroblasts, either as inclusions in garnet, inclusions within phases included in the garnet, or within the plagioclase and hornblende haloes around the porphyroblastic garnet. Sixteen zircon grains were analyzed from the matrix. Two of these 16 grains were within matrix hornblende, one was within matrix quartz, and the remaining one was are located along grain boundaries.

The U contents from the analyses provide additional distinction between populations 1 and 2. Zircon composed of dominantly population 1 had >500 ppm U, whereas analyses of whole grains of population 2 were <100 ppm U (Table 4). Analyses with U concentrations between 100 and 500 ppm U were those from grains with significant fractions of cores and rims where the beam likely sampled both domains (e.g., Figs. 12C–12D).

Two zircon analyses that sampled predominantly population 1 domains (z77, z79; Figs. 12A–12B) yielded concordant analyses (at 1σ) at ca. 2.45 Ga (²⁰⁷Pb/²⁰⁶Pb dates of 2458 ± 3 Ma and 2442 ± 7 Ma; Fig. 13A), with U concentrations in excess of 500 ppm (Table 3; Fig. 13B). Grain z79 is in garnet (Fig. 12B), and z77 is completely included in plagioclase, which is itself an inclusion in garnet (Figs. 12A and 13C). Analyses that sampled both cores and rims in multidomain grains were discordant (Figs. 12C–12D), with ²⁰⁷Pb/²⁰⁶Pb dates ranging between 1940 Ma ± 15 Ma (z41) and 2447 ± 4 Ma (z91). These analyses also had U concentrations >100 ppm. The array of discordant analyses likely reflects mixing as a result of sampling two age domains coupled with recent Pb loss. The recent Pb loss had a more noticeable effect on analyses sampling mixed domains of predominantly population 1, with the higher U concentrations (Table 4; Fig. 13B). Analyses of exclusively population 2 zircon were either concordant or reversely discordant, and together yielded a ²⁰⁷Pb/²⁰⁶Pb date of 1732 ± 21 Ma (2σ, MSWD = 1.5, n = 12). The seven analyses that overlapped concordia at 1σ have a weighted mean ²⁰⁷Pb/²⁰⁶Pb date of 1737 ± 28 Ma (2σ, MSWD = 1.8; Table 3; z71 and z50 are examples; Figs. 12E, 12F, and 13, inset). This date is our preferred timing of population 2 zircon growth.

At least three phases of deformation affected the northwest Gallatin Peak terrane. D₁ structures are likely composite, and zircon U-Pb data from the mafic dike point to an Archean age for at least some component of this deformation. Population 1 zircon is unlikely to be inherited because it would require the dike magma to have sampled an as-yet-unrecognized 2.45 Ga source while avoiding the widely recognized 3.2–3.0 Ga sources in the region. Thus, based on the most concordant data from population 1 zircon, the minimum age of the mafic dike that obliquely crosscuts S₁ fabric in the host tonalitic gneiss (Figs. 4C–4D) is ca. 2.45 Ga. Whether this population of zircon is magmatic or metamorphic is not yet clear. Support for the latter comes from a ca. 2.45 Ga metamorphic event recognized in the Tobacco Root Mountains and elsewhere in southwest Montana (Roberts et al., 2002; Cheney et al., 2004; Foster et al., 2006; Loehn, 2009; Krogh et al., 2011) and Utah (Mueller et al., 2011). The 2.45 Ga age may be synchronous with the deformation event that produced the foliation locally recorded as inclusion trails in garnet from this dike (also listed in the broad category of D₁ structures here). These data also imply that at least this dike in the Gallatin Peak terrane cannot be part of the suite from which a 2.06 Ga mafic dike was dated in the Tobacco Root Mountains (Mueller et al., 2004, 2005; Brady et al., 2004b).

D₂ structures are the most prominent deformation features in the terrane and studied samples. The prograde and peak mineral phases in the Bear Basin schist (e.g., Ky, St, Sil, Bt) and in the mafic dike (e.g., Hbl) are aligned with this fabric (including Ky locally defining L₂ mineral lineation), suggesting that peak metamorphism was synkinematic with respect to D₂.

Local kyanite and staurolite inclusions in garnet contrast with sillimanite occurrence solely in the matrix, suggesting that sillimanite stability followed that of kyanite and staurolite on a portion of a clockwise P-T path recorded by the Bear Basin schist. Thermobarometry from peak synkinematic assemblages in all samples yielded P-T conditions of ∼0.9 GPa and ∼700 °C. These conditions are near the Ky-Sil reaction line (Fig. 9), which is consistent with the presence of both aluminosilicates in the schist. Peak temperatures exceeded those predicted for staurolite stability from phase assemblage modeling. This is consistent with staurolite solely occurring as inclusions in garnet in some samples (Fig. 7C). These results are similar to undated conditions reported by Salt (1987) and Mogk et al. (1992) for the Gallatin Peak terrane. Clockwise P-T paths are inferred for other exhumed crystalline rocks in the Tobacco Root Mountains (Cheney et al., 2004) and the Ruby Range (Alcock et al., 2013).

The new data presented here indicate that D₂ deformation and prograde and peak metamorphism occurred in the northwest Gallatin Peak terrane ca. 1750–1720 Ma. This is supported both by zircon U-Pb data from the deformed mafic dike and by monazite U-Th–total Pb data from the Bear Basin schist. The morphology and lack of zoning in population 2 zircon domains are typical of metamorphic zircon (Corfu et al., 2003). These grains, with a weighted mean ²⁰⁷Pb/²⁰⁶Pb date of 1737 ± 28 Ma, occur as inclusions in garnet as well as in the matrix associated with S₂-aligned hornblende or plagioclase. Thus, garnet in the mafic dikes grew near or after ca. 1740 Ma. These observations also indicate that the S₂ fabric, which wraps around the garnet porphyroblasts, developed during or after garnet growth.

U-Pb analyses of zircon with population 1 and 2 domains are in some cases discordant, in part because the ion probe primary beam likely sampled both domains (Figs. 12C–12D). Results thus reflect both mixing between ca. 2.45 Ga cores and ca. 1.74 Ga rims and recent Pb loss. Grains with higher U are more susceptible to radiation damage-enhanced Pb loss, impacting the degree of discordance of these analyses (Fig. 13B).

Monazite populations 1b and 2 (1750–1720 Ma) are interpreted to have grown during prograde and peak metamorphism because they occur as inclusions in all major metamorphic phases (garnet, kyanite, staurolite) and in the matrix of the Bear Basin schist. In general, monazite in low-Ca pelitic bulk compositions, such as the Bear Basin schist, is commonly thought to become stable around the staurolite-in isograd (Kohn and Malloy, 2004) and to continue growing throughout the duration of a metamorphic event (Spear and Pyle, 2010).

An interpretation of episodic synmetamorphic monazite growth is further supported by distinct compositional variations between the populations. For example, Y content increases by a factor of two between populations 1b and 2 in sample GP7c, whereas the opposite trend is recorded in sample GP7f. However, both trends can be explained by prograde metamorphic growth, considering that monazite Y compositions are commonly dependent on the behavior of other Y-bearing phases such as xenotime and garnet (e.g., Spear and Pyle, 2010). In xenotime-bearing lithologies, monazite Y concentrations generally increase with increasing metamorphic grade (Pyle et al., 2001), tracking the progressive destabilization of xenotime. Observed xenotime inclusions in staurolite and total Y contents in monazite that are approximately an order of magnitude higher in GP7c than GP7f indicate that xenotime was a stable phase for a portion of the prograde sequence in GP7c. In contrast, xenotime was probably never stable in GP7f, and the decreasing trend in Y content in monazite between GP7f populations 1b and 2 probably reflects progressive depletion in available bulk Y as both garnet and monazite continued to grow. Thus, the range in population 1b and 2 monazite dates is interpreted to represent episodic growth through the main duration of late Paleoproterozoic metamorphism. Cheney et al. (2004) came to a similar conclusion for monazite data from an ion probe study of rocks in the Tobacco Root range.

Monazite population 3 has the highest Y content and dates that cluster around ca. 1705 Ma. The increase in Y is interpreted to reflect retrograde breakdown of garnet, which liberates Y into the system (Pyle et al., 2001; Gibson et al., 2004; Mahan et al., 2006). This is supported by evidence of late garnet resorption (Figs. 8A–8B), which may be due to fluid flow and/or exhumation after peak thermotectonism. Finally, the one Archean monazite grain identified in the Bear Basin schist (population 1a) could be either metamorphic or detrital. Weyand (1989) dated 3.1 Ga detrital zircon in schist from this unit, and a detrital origin for this similarly aged monazite is consistent with its textural location in a quartzofeldspathic layer (Bt + Pl + Qz only). Possible detrital sources are nearby, since the unit appears to sit in depositional contact with 3.2–3.1 Ga orthogneiss.

Our new results, combined with previously published data, suggest that exhumed tracts of rocks characterized by broadly decreasing Proterozoic paleodepth are exposed across the Northern Madison Range, and may extend farther to the southeast (Figs. 1 and 2). Of reported thermobarometric conditions across the range, two have timing constraints and are known to represent ca. 1.7 Ga paleodepths. These are ∼45–40 km near the Madison River (1.2 GPa metamorphic pressure in Bear Trap Canyon; Ault et al., 2012), which may be the deepest in southwest Montana, and ∼30 km in Bear Basin (∼0.9 GPa; this study). Thus, as much as a 15 km difference in paleodepth is recorded over a map distance of 30 km across strike of the dominant structural grain (Fig. 2).

K-Ar mica data were used to demarcate the transition in southwest Montana from complete to negligible Paleoproterozoic thermal resetting toward the interior of the Wyoming craton (Giletti, 1966). The original location of this boundary (Giletti’s line) passes through Gallatin Canyon, ∼15 km across strike and southeast from Bear Basin (Figs. 1 and 2). This would appear to require a rapid shallowing to the southeast of Proterozoic paleodepths to upper-crustal levels (i.e., ∼10 km) over 15 km from our farthest east paleodepth constraint of ∼30 km. However, problems with excess Ar brought into question the validity of the key Gallatin Canyon data (Giletti, 1971), and the true location of this thermal boundary in the Northern Madison range is still uncertain.

A tilted partial crustal section, stacked thrust sheets, or stacked fold nappes are among several viable options that could currently explain the constrained differences in paleodepth. Several previously identified but as yet poorly understood high-strain zones within the range, including the Hellroaring Creek shear zone (Fig. 2), may represent bounding structures of thrust sheets or nappes. More detailed work in the region is clearly justified given the general utility of partial crustal sections for better understanding the growth and evolution of orogenic crust (Percival et al., 1992; Miller and Snoke, 2009). Farther to the south, faulting at greenschist facies in the Southern Madison Range (Fig. 1; Madison Mylonite zone of Erslev and Sutter, 1990) and along-strike localities to the northeast (Mogk and Henry, 1988; Erslev, 1989) may represent upper-crustal deformation in the foreland region of the Big Sky orogeny (O’Neill, 1998; Harms et al., 2004a).

The spatial extent of the Big Sky orogen is poorly constrained, partly because it is only exposed in younger, isolated basement-cored uplifts (Fig. 1A). Previous studies invoke a ca. 1.80–1.71 Ga time interval for the Big Sky orogeny (e.g., Harms et al., 2004a). Some have suggested that it may represent the closure of a late Paleoproterozoic ocean basin (O’Neill, 1998; Roberts et al., 2002; Harms et al., 2004a; Vogl et al., 2004; Mueller et al., 2005; Foster et al., 2006; Alcock et al., 2013), recorded by ca. 1.86 −1.81 Ga arc-related granitoid rocks in the Little Belt Mountains (Fig. 1A; Mueller et al., 2002; Vogl et al., 2004).

New data presented here, along with previously published monazite and zircon geochronology, suggest a southeastward younging of high-grade metamorphism in SW Montana. The details of this observation provide important constraints on the pace, duration, and spatial patterns of hinterland growth in the Big Sky orogen. Figure 14A shows a compilation of all published zircon, monazite, and garnet geochronological data younger than 1900 Ma that we are aware of, plotted along a 150 km SE-trending profile centered among the basement exposures of SW Montana. The orientation of this profile was chosen to cross the structural grain of the dominant strike of major structures (e.g., foliation planes, fold axes, shear zones), the overall trend of Giletti’s line, and the Great Falls tectonic zone. This figure focuses on data that are linked to timing of peak metamorphic conditions equal to or exceeding 700 °C, although igneous zircon U-Pb data and ⁴⁰Ar-³⁹Ar data for mica and amphibole are also plotted for additional context. The southeastward younging trend of high-grade metamorphism is separately apparent from both monazite and zircon data sets and from isotopic (ion probe) and nonisotopic (U-Th–total Pb) techniques (Figs. 14B and 14C), suggesting no systematic bias from different mineral systems or analytical approaches. Additional description of the plotted data is given in the following paragraphs.

Within the Great Falls tectonic zone, there is evidence of late Paleoproterozoic arc magmatism in the Little Belt Mountains (Fig. 1A), Pioneer Mountains, and Biltmore anticline (Fig. 14A; Mueller et al., 2002; Vogl et al., 2004; Foster et al., 2006). Igneous zircon crystallization ages range from ca. 1890 to 1810 Ma (Vogl et al., 2004; Foster et al., 2006), while subsequent metamorphism and anatexis may have occurred in the Little Belt Mountains at ca. 1820−1790 Ma (Vogl et al., 2004).

In the most northwestern exposures of non-arc-related Precambrian rocks in the Highland Mountains (Fig. 1B), previous U-Pb zircon and monazite and Pb-Pb garnet studies concluded that peak metamorphism occurred ca. 1800–1770 Ma (O’Neill et al., 1988; Roberts et al., 2002; Mueller et al., 2005). Farther southeast across strike (∼40 km) in the Ruby Range, Pb-Pb garnet dates were similarly interpreted to record garnet growth and metamorphism from ca. 1800 to 1780 Ma (Roberts et al., 2002). This history is also supported by a more recent monazite U-Th–total Pb study of migmatitic metapelites (Alcock et al., 2013). These authors interpreted prograde metamorphism as early as ca. 1830 Ma, but culminating at ca. 1780 Ma. In the Tobacco Root Mountains, U-Pb dates from leucosomal zircon in migmatitic paragneiss and zircon in a metamorphosed mafic amphibolite dike are interpreted to record peak lower-granulite-facies metamorphism at ca. 1770 Ma (Mueller et al., 2004). A monazite U-Pb study from a variety of aluminous lithologies was alternatively interpreted to constrain prograde mineral growth from ca. 1780 Ma to peak conditions closer to ca. 1755 Ma (Cheney et al., 2004a). Continuing to the southeast in the Bear Trap Canyon area of the Northern Madison Range (Figs. 1 and 2), a deformed mafic dike records high-pressure granulite-facies metamorphism at ca. 1750 Ma (U-Pb zircon; Ault et al., 2012). Our new data from the central portion of the range are interpreted to record prograde to peak metamorphism from ca. 1750 to 1720 Ma. In summary, the younging of high-grade tectonometamorphism by at least 40 m.y., and perhaps as much as 80 m.y., extends across at least 80 km from the Highland Range to at least the central Northern Madison Range (Fig. 14). This suggests growth or migration of the hinterland of the orogen toward its foreland, where the latter is considered to be represented by substantially cooler Proterozoic temperatures (e.g., Giletti’s line or modern equivalent) and discrete greenschist-grade structures like the Madison mylonite zone (Fig. 1B).

Several explanations for this spatio-temporal pattern of Paleoproterozoic metamorphism appear plausible. First, the simplest interpretation is that the observed southeastward younging of tectonism may represent synconvergent lateral growth of the metamorphic core of the orogen by in-sequence thrusting and thickening (e.g., Willett et al., 1993), caused by either protracted or episodically accelerated collision. Similar propagation patterns representing >50 m.y. of continued convergence and growth of deformation and metamorphism into the foreland are observed in other well-known major collisional belts, such as the ca. 1090–995 Ma Grenville (e.g., Jamieson et al., 1995; Hynes and Rivers, 2010), the 55 Ma to present Himalayan (e.g., Hodges, 2000) orogens, and Jurassic to Early Cretaceous foreland-directed growth of the northern Canadian Cordilleran hinterland (Staples et al., 2014). Alternatively, synconvergent thickening of the core of the Big Sky orogen may have been followed by postconvergent shortening and associated metamorphism from thickening on the flanks of the mountain belt driven by gravitational spreading. This phenomenon is observed in numerical models and was invoked by Jamieson and Beaumont (2011) to explain the synconvergent Ottowan phase (1090–1020 Ma) of high-grade metamorphism in the core of the Grenville orogen versus the later Rigolet phase (1005–995 Ma) of thrusting and medium-grade metamorphism in the Grenville front tectonic zone. Other possibilities include one or a combination of additional tectonic processes that could have resulted in juxtaposition of domains with older and younger tectonometamorphic signatures, such as gneiss dome formation, channel flow, or later strike-slip displacements. Given our current understanding of the Big Sky orogeny, we lack clear constraints on which of these models is most plausible. Each of these scenarios could be tested through further structural and geochronological studies. For instance, expanding our current understanding of the poorly constrained high-strain zones within the Northern Madison Range may yield potential constraints on one or more of these models.

Our results suggest a systematic pattern in the age of peak metamorphism from northwest to southeast, but subsequent cooling histories and exhumation patterns are less well constrained. Several studies have reported ca. 1.8–1.7 Ga ⁴⁰Ar/³⁹Ar hornblende and biotite cooling dates, generally corroborating earlier work by Giletti (1966), although no ⁴⁰Ar/³⁹Ar data exist yet from the Northern Madison Range. Harlan et al. (1996) reported ca. 1.80 Ga ⁴⁰Ar/³⁹Ar dates from the Highland Mountains. A large ⁴⁰Ar/³⁹Ar data set from the Highland, Ruby, and Tobacco Roots Mountains yielded a mean of ca. 1.76 Ga (Roberts et al., 2002). Both of these results are older than the interpreted timing of peak metamorphism in the central Northern Madison Range (Figs. 1, 2, and 14), suggesting that these rocks had cooled through 300 °C and been exhumed prior to the time when other portions of the orogen experienced their peak conditions. However, a second and equally large thermochronology data set from the Ruby and Tobacco Root Mountains suggests significantly younger and rapid ca. 1.71 Ga cooling through the 500–300 °C temperature interval (Brady et al., 2004a). This is also supported by the youngest monazite (ca. 1730–1713 Ma) from the Tobacco Root Mountains, which is interpreted to have grown during retrograde decompression (Cheney et al., 2004). The youngest dates in that study (1713 Ma) are similar to those in this study (ca. 1705 Ma) in the Northern Madison Range, which is also interpreted to be related to garnet breakdown during exhumation. In summary, both differential and simultaneous exhumation may have occurred across different parts of the region, but more work is needed to clarify this history.

Significant uncertainty remains regarding how the kinematics of deformation varied spatially and temporally across the region. For example, Mueller et al. (2005) and Foster et al. (2006) pointed out that magmatism and associated metamorphism occurred almost 100 m.y. earlier in the Little Belt Mountains (ca. 1.86–1.81 Ga) than in southwest Montana, and they used this to suggest that the Big Sky orogeny may be a distinctly younger and separate collisional event. They envisioned an early orthogonal 1.9–1.8 Ga collision event, associated with suturing of the Archean Medicine Hat and Wyoming cratons (Great Falls tectonic zone), evolving to convergence with a stronger transpressional component associated with accretion of a separate Paleoproterozoic terrane during the Big Sky orogeny (e.g., Selway terrane of Foster et al., 2006). If so, a kinematic record of these transitions is likely preserved in the basement rocks of southwest Montana.

Spatial and temporal patterns of deformation and metamorphism yield insight into the kinematics, rheology, and overall tectonic significance of exhumed paleo-orogens, and they allow better understanding of the processes at work in active orogens of similar scales. Thermotectonism associated with the Big Sky orogeny modified much of southwest Montana from ca. 1.8 to 1.7 Ga. This contribution expands the known width of the high-grade metamorphic core of the orogen by ∼40 km farther southeast across strike from the northwest end of the Northern Madison Range. Structural relationships, metamorphic petrology, and U-Pb geochronology indicate that at least two episodes of late Paleoproterozoic deformation (D₂ and D₃), including the dominantly NE-striking and moderately SE-dipping foliation and associated isoclinal folds and lineation, affected rocks in the central Northern Madison Range. Synkinematic metamorphism, recorded in both mafic dikes and pelitic schist, reached peak conditions of 0.9–0.8 GPa and ∼700 °C at 1750–1720 Ma. Combined with other regional data, a pattern of southeastward younging of high-grade metamorphism by 80–40 m.y. is apparent across an ∼100 km transect from the Highland Mountains through the study area. Although additional work is needed to fully fingerprint the particular tectonic mechanism(s) at play here, this pattern suggests that lateral growth of the orogenic core toward the foreland during protracted collision is likely to be an important characteristic of one of North America’s most recently recognized major convergent belts.

We would like to thank Katherine Barnhart, Kevin Toeneboehn, Diana Rattanasith, and Omero Orlandini for assistance in the field, and Karl Kellogg for sharing his previous mapping work from the Spanish Peaks region. Financial support for this project was provided by National Science Foundation grant EAR-1252295 to Mahan, U.S. Geological Survey EDMAP grant G12AC20182 to Mahan and Condit, and student research grants from the Geological Society of America, Rocky Mountain Association of Geologist Stone-Holberg Scholarship, Tobacco Roots Geological Society, and Colorado Scientific Society to Condit. Axel Schmitt (formerly at University of California–Los Angeles), Kevin Chamberlain (University of Wyoming), and Mike Jercinovic (University of Massachusetts) are thanked for their assistance in their respective geochronology laboratories. Tekla Harms, Paul Mueller, David Foster, and Dan Gibson are thanked for constructive reviews of earlier versions of this manuscript.