Integration of petrologic, chronologic and petrophysical xenolith data with geophysical observations can offer fundamental insights into understanding the evolution of continental crust. We present the results of a deep crustal xenolith study from the northern Rocky Mountain region of the western U.S., where seismic experiments reveal an anomalously thick (10–30 km), high seismic velocity (compressional body wave, Vp > 7.0 km/s) lower crustal layer, herein referred to as the 7.x layer. Xenoliths exhumed by Eocene minettes from the Bearpaw Mountains of central Montana, within the Great Falls tectonic zone, include mafic and intermediate garnet granulites, mafic hornblende eclogite, and felsic granulites. Calculated pressures of 0.6–1.5 GPa are consistent with derivation from 23–54 km depths. Samples record diverse and commonly polymetamorphic pressure-temperature histories including prograde burial and episodes of decompression. Samples with barometrically determined depths consistent with residence within the seismically defined 7.x layer have calculated bulk P-wave velocities of 6.9–7.8 km/s, indicating heterogeneity in the layer. Shallower samples have markedly slower velocities consistent with seismic models. New monazite total U-Th-Pb data and a variety of additional published geochronology indicate a prolonged and episodic metamorphic history, beginning with protolith ages as old as Archean and followed by metamorphic and deep crustal fluid-flow events ca. 2.1 Ga, 1.8–1.7 Ga, and 1.5–1.3 Ga. We suggest that the 7.x layer in this region owes its character to a variety of processes, including magmatic underplating and intraplating, associated with multiple tectonic events from the Neoarchean to the Mesoproterozoic.

Continental crust provides an integrated record of crustal growth and differentiation through time. The lower crust is central to this evolution and its composition and degree of heterogeneity can have a profound influence on the rheology of the lithosphere, yet it is commonly the least constrained part of the lithospheric column, because data density and seismic constraints generally decrease with depth. Although the average composition of the lower crust is generally considered to be mafic, debate continues on this topic as a result of rare exposure and an incomplete and potentially biased record made available through lower crustal xenoliths (e.g., Rudnick and Taylor, 1987; cf. Hacker et al., 2011). Discrete layers of suspected mafic lower crust are identifiable as having high compressional wave seismic velocities (Vp) of >7.0 km/s. Such high-velocity lower crustal layers have been observed in a variety of settings, particularly in shields and platforms, continental arcs, and rift zones (Fig. 1) where the material is commonly considered a discrete product of mantle-derived magmatic underplating or intraplating (Warner, 1990; Kay and Mahlburg-Kay, 1991; Eaton, 2006; Karlstrom et al., 2005; Crowley et al. 2006; Cornwell et al., 2010; Ridley and Richards, 2010). Imbrication of oceanic crust has also been proposed in some tectonic settings (e.g., New Zealand; Okaya et al., 2007).

Although high seismic velocity layers are not uncommon, worldwide compilations show that they are generally relatively thin, averaging 3–8 km (Fig. 1; Christensen and Mooney, 1995; Rudnick and Fountain, 1995). In contrast, several seismic experiments in the Rocky Mountain region of Montana and Wyoming (USA), and Alberta (Canada) reveal a distinctly thick (10–30 km) lower crustal layer with Vp between 7.0 and 7.8 km/s (Henstock et al., 1998; Snelson et al., 1998; Clowes et al., 2002; Gorman et al., 2002; Eaton, 2006; Schutt et al., 2008). The unusual thickness of this high-velocity layer, which constitutes as much as half of the 40–60 km crustal column, makes it arguably one of the most anomalous regions of continental crust in North America. Initial observations indicate that Archean and/or Proterozoic magmatic underplating events may have played a role in contributing to the thickness of high-velocity material (Gorman et al., 2002: Chamberlain et al., 2003). However, the age and character of the lower crust in this region have until now remained largely unknown.

Xenoliths afford a unique opportunity to characterize the physical and chemical characteristics of the lower crust, and the data can be linked with modern geophysical studies (e.g., Rudnick and Taylor, 1987; Kay and Mahlburg-Kay, 1991). As extensive geophysical data sets become available for the western United States through large projects such as the National Science Foundation EarthScope program, it will be increasingly important to link these observations to compositional, petrologic, and chronologic information in order to both better refine the regional tectonic history and the range of composition and petrology corresponding to a particular set of seismic properties. Although observations from xenoliths are incomplete, perhaps unpredictably biased, and depth constraints have limitations, xenoliths commonly provide the best method of directly sampling relatively modern lower crust.

Here we report data from crustal xenoliths exhumed in Eocene minettes in the Great Falls tectonic zone in central Montana. The data include thermobarometry and calculated bulk seismic velocities as well as petrology and monazite geochronology documenting a dynamic and polyphase Paleoproterozoic history. We emphasize the heterogeneity in physical properties and history in the xenolith record, and explore the possible means of assembly for the high-velocity lower crustal layer in this region.

The Rocky Mountains of Montana and Wyoming record a protracted history of crustal growth and reactivation. This region consists of Archean cratons amalgamated in the Proterozoic and accreted Proterozoic terranes, both of which have been covered by Phanerozoic sedimentary rocks and late Mesozoic and Cenozoic volcanic rocks (Figs. 2A and 3). The principal features of this region that are relevant to this study are described in more detail in the following (and in the Discussion).

Medicine Hat Block, Wyoming Craton, and Great Falls Tectonic Zone

The Medicine Hat block is completely covered by the Western Canada Sedimentary Basin and has been studied exclusively by geophysical methods, xenoliths, and drill core samples (e.g., Ross, 2002). Drill core samples reveal that the upper crystalline crust of the Medicine Hat block is made of ca. 3.27–2.65 Ga metaplutonic gneiss (Villeneuve et al., 1993). This block is bounded to the north by the aeromagnetically defined Vulcan low, which is thought to represent a Proterozoic collisional boundary with the Hearne craton (Ross, 2002). Possible ages of collision are based on very sparse drill core data, including a single discordant ca. 2.1 Ga U-Pb analysis of titanite (Villeneuve et al., 1993) and ca. 1.8 Ga K-Ar biotite dates (Burwash et al., 1962).

The basement rocks of the Wyoming craton are observed in discrete Laramide-age uplifts and through xenolith studies, and contain three main domains that share a common isotopic fingerprint but have different Archean histories (Chamberlain et al., 2003; Frost, 1993; Mueller and Frost, 2006). The Montana metasedimentary province contains ca. 3.3–3.1 Ga igneous rocks and intercalated metasediments (Mueller et al., 1996). The Beartooth-Bighorn magmatic zone is primarily ca. 2.9–2.7 Ga metaplutonic tonalite-trondhjemite-granodiorite (Frost et al., 2006a). The southern accreted terranes were joined to the rest of the Wyoming craton by terrane accretion and arc magmatism at 2.68–2.50 Ga (Frost et al., 2006b). Dike swarms (ca. 2.1 Ga) are reported from exposed rocks at both the northern and southern margins of the Wyoming craton and are interpreted as rift related (Premo and Van Schmus, 1989; Cox et al., 2000; Mueller et al., 2004).

The Wyoming craton and Medicine Hat block are sutured together by the ca. 1.86–1.71 Ga Great Falls tectonic zone (Giletti, 1966; O’Neill and Lopez, 1985; Mueller et al., 2002; Harms et al., 2004). Some workers have speculated that the Great Falls tectonic zone may be an intracontinental tectonic feature (e.g., Boerner et al., 1998). However, more recent geological studies have identified ca. 1.85 Ga igneous rocks with arc affinities (Mueller et al., 2002), and high-grade contractional tectonism suggesting a major collisional boundary (Mueller et al., 2002; Harms et al., 2004). The timing of deformation and metamorphism attributed to activity along the Great Falls tectonic zone in southwestern Montana extends to as late as 1.71 Ga (Roberts et al., 2002; Harms et al., 2004). The craton is bounded to the south by the ca. 1.78–1.75 Ga Cheyenne belt (Karlstrom and Houston, 1984), and both cratonic blocks are bounded to the east by the ca. 1.83–1.72 Ga Trans-Hudson orogen (Bickford et al., 1990; Dahl et al., 1999).

Belt Basin Development

Mesoproterozoic sediments of the Belt Basin cover portions of the northwestern Wyoming craton and western Medicine Hat block. The Belt Basin reaches 15–20 km in thickness and defines a north-northwest–oriented rift basin near the western margin of North America (e.g., Harrison, 1972; Link et al., 1993). Rocks of the basin have also been identified in the subsurface of central Montana (Fig. 3). The sedimentary sequence is thought to have been deposited between 1.47 and 1.37 Ga (Evans et al., 2000). Mafic dikes emplaced across the central Wyoming craton between 1.47 and 1.45 Ga are interpreted to reflect an eastern rift arm of the Belt Basin (Chamberlain et al., 2003).

Montana Alkalic Province and Crustal Xenoliths

The projected trace of the Great Falls tectonic zone and the Wyoming craton are intruded by suites of potassic and alkali volcanics and kimberlites, including the Highwood Mountains, Sweet Grass Hills, Eagle Buttes, Little Rocky Mountains, Crazy Mountains, and the Williams and Homestead kimberlites. Many of these localities contain upper mantle and lower crustal xenoliths (Collerson et al., 1989; Hearn, 1989; Hearn et al., 1989; Joswiak, 1992; Carlson and Irving, 1994; Downes et al., 2004; Bolhar et al., 2007; Facer et al., 2009; Blackburn et al., 2010, 2011). Crustal xenoliths described in this study are from the Robinson Ranch and Little Sand Creek localities and were exhumed by Eocene minettes (54–50 Ma) of the Bearpaw Mountain volcanic field in north-central Montana (Marvin et al., 1980; Hearn, 1989; MacDonald et al., 1992) and within the Great Falls tectonic zone (Figs. 2A and 3).

Xenoliths from Little Sand Creek yielded U-Pb SHRIMP (sensitive high-resolution ion microprobe) dates on zircon that range from ca. 3.0 to 1.8 Ga (Bolhar et al., 2007). Northwest of the Bear Paw Mountains at the Sweet Grass Hills xenolith locality (Fig. 3), granulite xenoliths yield U-Pb zircon dates of ca. 1.8 Ga and Nd isochrons (garnet, clinopyroxene, and whole rock) ranging from 1.7 to 1.5 Ga and are interpreted to reflect modification of the lower crust (Davis et al., 1995). Davis et al. (1995) postulated that the mafic granulites represent either metamorphosed Archean crust or an addition to the lower crust ca. 1.8 Ga.

Modern Crustal Structure

Crustal thickness in Montana and Wyoming varies from 49 to 60 km, based on the Deep Probe active source seismic experiment (Snelson et al., 1998; Gorman et al., 2002), although more recent passive source seismic studies indicate somewhat lower crustal thickness estimates of 39–50 km for central Montana (e.g., Bensen et al., 2009; Gilbert, 2012). Much of the region is underlain by a layer of anomalously thick and seismically fast lower crust, referred to here as the 7.x layer (Fig. 2B; Gorman et al., 2002). Originally identified in the analysis of COCORP seismic lines (Morel-a-l’Huissier et al., 1987), the active source seismic refraction experiments Deep Probe, Lithoprobe, and SAREX (Southern Alberta Refraction Experiment) further refined the north-south extent and thickness variation of the 7.x layer (Snelson et al., 1998; Clowes et al., 2002; Gorman et al., 2002). It is a 10–30-km-thick layer that forms the lower crust of the northern two-thirds of the Wyoming craton, the Great Falls tectonic zone, and much of the Medicine Hat block (Gorman et al., 2002). Within the Wyoming craton, the 7.x layer has Vp values of 7.0–7.7 km/s that generally increase northward and an average thickness of 25 km. Within the Medicine Hat block, the layer has higher average Vp (7.6–7.9 km/s), appears thinner (16 km thick), and may display more upper surface topography.

Crust and upper mantle tomography from the Billings (Montana) seismic array is consistent with the presence of high-velocity lower crust as thick as 19 km east of the Deep Probe line (Schutt et al., 2008). To the west, the high-velocity lower crustal layer extends to the Yellowstone Caldera (Stachnik et al., 2008). In southern Wyoming, the high-velocity lower crustal layer is observed to end to the north of the Cheyenne belt (Rumpfhuber et al., 2009). Various mechanisms and ages have been proposed for the formation of the 7.x layer (Table 1), and these are explored further in the Discussion.

Analytical Methods

Modal mineralogy was quantified by automated scanning electron microscope analysis (Quantitative Evaluation of Minerals by SCANning, herein referred to as QEMSCAN) at the Advanced Mineralogy Research Center at the Colorado School of Mines. This technique employs a combination of simultaneously collected energy-dispersive spectroscopy (EDS) spectra and calibrated backscatter electron intensity to assign mineral identity to each mapped pixel (Pirrie et al., 2004; Hoal et al., 2009). The instrument is a Carl Zeiss EVO50 scanning electron microscope with four energy dispersive spectrometers. Instrument conditions for full thin section mapping were 25 kV accelerating voltage, 5 nA sample current, a beam size of 0.25–0.5 μm, pixel step sizes of 20–30 μm, and a variable dwell time necessary to accrue 1000 counts. The resulting field image maps (Fig. 4) are used to make general textural observations and to determine the textural setting of accessory phases. Because the xenoliths are too small for whole-rock measurements, mineral modes and compositions were used to estimate bulk major element chemistry for general characterization, petrological modeling, and bulk seismic velocity calculations. This method also allows avoidance of complications due to decompression during eruption and late-stage alteration effects on the bulk properties (e.g., Rudnick and Jackson, 1995) and allows consideration of effective bulk compositions for petrologic modeling (e.g., Stüwe, 1997).

Quantitative mineral compositions were determined using a JEOL JXA-8600 electron microprobe at the University of Colorado. Before quantitative analysis, X-ray maps were used to evaluate potential zoning in phases such as garnet and feldspar. For mapping, beam conditions were 15 kV and 100 nA. Dwell time was 40 ms per pixel. X-ray maps collected Mg Kα and Al Kα on TAP crystals, Ca Kα on a PET crystal, and Mn Kα on a LIF crystal. For quantitative analysis, beam size was ∼1 μm (focused) for all phases except feldspar (5 μm beam) and column conditions were 15 kV and 20 nA. A common suite of natural and synthetic standards was used. Full mineral compositions are reported in the Supplemental File1.

General Textural Observations, Bulk Compositions, and Mineral Assemblages

We analyzed 11 xenolith samples: 8 from the Robinson Ranch suite (ROB) and 3 from the Little Sand Creek suite (LSC). The samples display a diverse range of modal mineralogy, grain size, and textural characteristics (Table 2). All of the analyzed samples are garnet and feldspar bearing but with a wide range of abundances (e.g., garnet abundance ∼12%–46%, feldspar ∼5%–50%). Pyroxene and/or quartz are present in only about half the samples. Some samples contain deformation-related features such as gneissic layering defined by garnet-rich or quartzofeldspathic bands (Fig. 4; ROB5 and ROB7) and foliation defined by aligned biotite and/or amphibole (Fig. 4; ROB4 and ROB6), whereas others have generally granoblastic textures (Fig. 4; LSC04). Most samples preserve secondary products that are interpreted to have developed during volcanic transport from the lower crust to the surface. These include glassy kelyphitic rims and very fine grained and localized symplectite at grain boundaries (Padovani and Carter, 1977; Rudnick and Taylor, 1987; Messiga and Bettini, 1990). In addition, some samples contain porphyroblast-inclusion relationships or reaction textures that are interpreted to represent an important part of the pre-eruption evolution of the rocks (described herein). With one exception (LSC18; Fig. 4), late-stage surface alteration is minor.

Bulk major element compositions were calculated from modal proportions and major and minor mineral compositions (Table 3). For simplicity, the samples are divided into mafic (<52 wt% SiO2), intermediate (>52 wt% and <57 wt% SiO2), and felsic (>57 wt% SiO2) granulites. The mafic granulites are ROB3, ROB4, ROB8, ROB9, LSC04, and LSC06. Intermediate granulites are ROB1, ROB6, and LSC18, and the two felsic granulites are ROB5 and ROB7. On a total alkali versus silica diagram, the samples generally plot in the basalt-basaltic andesite or dacite-rhyolite fields (Fig. 5). The plotted range is similar to crustal xenoliths collected from the nearby Eagle Buttes (Joswiak, 1992) and from southern Wyoming and northern Colorado (Farmer et al., 2005; Mirnejad and Bell, 2008) with perhaps a tendency for lower alkali content in more mafic xenoliths.

Five of the six mafic granulites have >10 wt% bulk CaO (Table 3). The primary (peak) metamorphic mineral assemblage in these high-Ca mafic granulites, as well as in one of the intermediate granulites with similar CaO, FeO, and MgO (ROB1), is Grt + Cpx + Rt ± Pl ± Qz ± Ilm ± Hbl ± Bt ± Kfs (mineral abbreviations herein after Whitney and Evans, 2010). Quartz is a clear and significant component of the peak assemblage in only one of the five mafic granulites. However, rare small rounded quartz inclusions occur in the margins of garnet porphyroblasts in ROB8, suggesting that silica activity is near one in this sample as well. Plagioclase is interpreted as a stable component of the peak assemblage in all but one of the mafic granulites. In ROB9, plagioclase occurs as small rounded inclusions in some of the larger clinopyroxene grains (Fig. 6C) and as a fine-grained component of the recrystallized margins of clinopyroxene, but is otherwise absent from the matrix. The sixth mafic granulite (ROB3) has <2 wt% CaO and 22 wt% Al2O3 (sample is 46% garnet by mode) and its interpreted peak mineral assemblage is Grt + Bt + Pl + Kfs + Rt.

The remaining two intermediate granulites (other than ROB1) have distinctly different bulk compositions and mineral assemblages. ROB6 is a Grt + Opx + Oamph + Bt + Pl + Rt granulite gneiss. An extensive late reaction zone occurs between primary garnet and Opx1 that consists of Opx2 + Spl + Pl2 symplectite2 (Figs. 6B and 7B) (described in further detail herein). The primary assemblage in LSC18 is Grt + Pl + Kfs + Bt + Qz, but it is heavily altered such that ∼40% of the matrix is pectolite, apophyllite, and minor carbonate and zeolite, likely representing products of near-surface alteration. The two felsic granulites (ROB5 and ROB7) also have the peak assemblage Grt + Pl + Kfs + Bt + Qz ± Rt ± Ilm.

Ilmenite is present in all samples. However, in some samples where ilmenite is not considered part of the peak metamorphic assemblage, it occurs as either early inclusions in other minerals (e.g., garnet) or more commonly as secondary rims on rutile grains. Accessory zircon occurs in all samples, monazite is recognized in at least six samples (see Monazite Geochronology discussion), and apatite and titanite occur in most of the higher bulk Ca samples (Fig. 4; Table 2).

Mineral Compositions

Garnet is generally unzoned with the common exception of a slight increase in Ca and decrease in Mg# [XMg/(XMg + XFe) x 100] toward grain edges that is interpreted to represent diffusional exchange with matrix phases during isobaric cooling (Fig. 7A)3. Thus, core compositions were typically used for bulk composition estimation and pressure-temperature (P-T) calculations. Garnet compositions in the high CaO mafic granulites and in the high Ca intermediate granulite ROB1 are relatively rich in grossular (XGrs = 0.16–0.34) and Mg#s are variable among the samples (19–47) (Table A1 [see footnote 1]). All other samples have significantly lower grossular contents (XGrs = 0.01–0.06) and less variable Mg# (32–46).

Plagioclase grains are generally characterized by broad unzoned interiors with gradational decreases in anorthite content toward extreme margins, commonly spatially associated with garnet. This is interpreted to reflect diffusional exchange with garnet during isobaric cooling from high temperature. In most samples, plagioclase interiors have a relatively restricted composition of XAn between 0.23 and 0.41, where the range in each sample is generally <0.03 (Table A2 [see footnote 1]). The lowest anorthite content (XAn = 0.14) occurs in the low bulk Ca intermediate granulite ROB6 and as early Pl1 inclusions in coarse clinopyroxene in the hornblende eclogite ROB9. K-feldspar occurs in the primary assemblage of several samples and as a common secondary phase rimming corroded plagioclase in others, although little systematic compositional variation is observed. The latter may reflect high-temperature interaction with K-bearing fluids during a ca. 1.8 Ga metasomatic event identified in mantle xenoliths from this region (Carlson and Irving, 1994; Carlson et al., 2004; Rudnick et al., 1999;,Downes et al., 2004). K-feldspar in the hornblende eclogite (ROB9) contains ∼2.5% BaO.

All the high-Ca mafic granulites and one high-Ca intermediate granulite (ROB1) contain clinopyroxene (diopside or sodian diopside by the classification of Morimoto, 1988), whereas one intermediate granulite (ROB6) contains primary Opx. Pyroxene Mg#s range from 61 to 82 with generally little zoning or internal variation within each sample (Table A3 [see footnote 1]). Three exceptions where early pyroxene is extensively recrystallized and a later pyroxene has a distinctly different composition are described further in the following discussion. Amphibole occurs in two of the xenoliths (Table A4 [see footnote 1]). ROB6 contains an unzoned orthoamphibole that is gedrite (classification of Leake et al., 1997). ROB9 contains unzoned calcic amphibole that is ferropargasite and minor anthophyllite pseudomorphs that may have replaced orthopyroxene. All but two samples contain primary biotite, although with significantly variable abundances (Table 2). Mg#s range from 53 to 78 and TiO2 contents range from 3.3–6.3 wt% with little variation in each sample (Table A5 [see footnote 1]).

Polymetamorphic Samples

Of the 11 examined xenoliths, 4 contain either early or late assemblages and reaction textures that provide additional insight into the tectonometamorphic history of the deep crust in this region. ROB3 is a low-CaO (1.4 wt%), high-Al2O3 (22 wt%) quartz-free mafic granulite that contains 46% garnet by volume. The high modal proportions of garnet and biotite suggest that this xenolith may have undergone multiple episodes of melt extraction. For example, the SiO2-poor and [Al2O3 – (FeO + MgO + TiO2) – K2O] rich bulk composition is similar to diatexite residuum from extremely migmatized metasedimentary rocks in western Maine (Solar and Brown, 2001). The primary assemblage is Grt + Bt + Pl + Kfs + Rt, but inclusions in larger garnets preserve a distinct earlier assemblage. Two main garnet morphologies are present. Large (to 7 mm) anhedral grains contain abundant inclusions of ilmenite, biotite, and Zn-bearing hercynitic spinel (Fig. 6A). Smaller euhedral garnet grains, generally without inclusions, occur throughout the matrix. However, similar compositions of both garnet morphologies and the lack of significant major element zonation, other than that interpreted as having developed during retrograde diffusion, indicate homogenization of garnet compositions at relatively high temperature. The inclusion suite indicates an early metamorphic episode characterized by the assemblage Grt + Bt + Pl + Spl + Ilm. A distinct population of monazite is also associated with this early assemblage (see Monazite Geochronology discussion).

ROB6 is a low-CaO (1.2 wt%), high-MgO (13 wt%) quartz-free intermediate granulite gneiss. The primary assemblage is Grt + Opx + Oamph + Bt + Pl + Rt with a strong gneissic foliation defined by the distribution of garnet, pyroxene, and feldspar and a shape-preferred orientation of biotite and amphibole (gedrite). A symplectite intergrowth of Opx2, hercynitic spinel, and Pl2 is extensively developed along grain boundaries between Grt1 and Opx1 (Fig. 6B). The new Opx2 is distinctly lower in Mg# (51) and higher in total Al2O3 (9–10 wt%) than the peak Opx1 (Mg# 73, Al2O3 6–7 wt%; Table A3; see footnote 1). Development of these reaction domains prior to eruption in the volcanic host is suggested by (1) the large spatial extent and grain size of the reaction domains (several tens to hundreds of μm wide, grain sizes to several tens of μm) compared to the substantially narrower widths and grain sizes of kelyphite rims that are attributed to rapid heating and decompression during exhumation, (2) variations in Fe and Mg between reactants and products that correlate with proximity to grain boundaries and suggest late diffusional exchange of these elements, a process that is not expected to occur over significant distances during the short residence times (hours to days) in the erupting magma (Morin and Corriveau, 1996; Sparks et al., 2006; Rutherford, 2008), and (3) the presence of a distinct Proterozoic population of monazite spatially associated with the reaction domains (see Monazite Geochronology discussion).

ROB9 and LSC06 are high-Ca mafic granulites that contain Grt, Cpx (sodian diopside), and Rt as part of their peak assemblages. Plagioclase is not present as a matrix phase in ROB9. However, this sample contains plagioclase in two distinct pre-peak (M1) and post-peak (M3) metamorphism textural settings. First, small rounded inclusions of plagioclase (Pl1; XAn = 0.14) commonly occur in the cores of larger clinopyroxene grains (Cpx2; Fig. 6C). Second, the margins of peak Cpx2 (XJd = 0.17–0.18) are extensively recrystallized to symplectitic intergrowths of a less sodic Cpx3 (XJd = 0.05) and Pl3 (XAn = 0.24) (Figs. 6D and 7D). This sample also contains substantial hornblende (ferropargasite), biotite, and K-feldspar that occur in textural equilibrium with the peak assemblage, as well as anthophyllite pseudomorphs that may have replaced orthopyroxene (Table 2). Although plagioclase occurs in the matrix of LSC06 and is considered part of the peak assemblage, peak Cpx1 (XJd = 0.12) is similarly recrystallized to a less sodic Cpx2 (XJd = 0.05) and a new Pl2 (XAn = 0.26). Similar arguments to those here for ROB6 can be made for development of these reaction textures prior to exhumation in the volcanic host. For example, both recrystallized Cpx (locally whole grains are reconstituted) and very fine grained (glassy) kelyphite rims occur in sample LSC06, the latter locally appearing to have partly consumed the former.

In summary, some samples in this xenolith suite preserve evidence for a dynamic and polymetamorphic history in the lower crust. ROB3 contains evidence for an early spinel-bearing assemblage that may have initially developed at relatively low pressure prior to attaining peak conditions. Reaction textures and assemblages in ROB6 suggest that this sample underwent a relatively late reheating event. The high Ca mafic granulite ROB9, a hornblende eclogite (sensu lato), initially and ultimately resided in a P-T field where plagioclase was stable, but reached eclogite facies without stable plagioclase at its peak in the intervening period.

Pressure and temperature conditions were determined to further constrain the metamorphic history and potential residence depths of the xenoliths. Calculations used the program TWQ 2.34 (Berman, 1991; and the internally consistent thermodynamic database of Berman and Aranovich (1996, updated in 2007; Berman et al., 2007). In general, only well-calibrated reactions that are recommended in the TWQ program were used. Absolute errors are considered to be ∼±50 °C and ±0.1 GPa (Berman, 1991). Calculations for equilibrium P-T conditions were also made for three garnet granulite xenoliths from the Eagle Buttes (Fig. 3) using mineral compositions from Joswiak (1992) to allow for direct comparison with the new xenolith data presented here. Isochemical phase assemblage diagrams (pseudosections) were calculated for a subset of samples in order to better determine P-T stability limits for a particular phase (e.g., spinel stability in the garnet granulite ROB3; plagioclase stability in the hornblende eclogite ROB9). These calculations were made using PerpleX07 (Connolly and Petrini, 2002) and the 2007 updated version of the internally consistent database of Holland and Powell (1998). The calculation results, along with the specific reactions used and whether these data were supplemented with pseudosection calculations, are shown in Table 4. The Ti-in-biotite thermometer of Henry et al. (2005), which is calibrated for 0.4–0.6 GPa, was also used to estimate temperatures for samples with pressures <0.8 GPa. The absence of observed graphite may mean that these temperatures are minimum estimates (Henry et al., 2005). Pressures calculated with reactions involving quartz are maximum estimates for samples where quartz is not observed. Calculated temperatures may be minimum estimates due to some degree of diffusional exchange of Fe and Mg during cooling (Frost and Chacko, 1989; Spear and Florence, 1992).

Calculated peak pressures range from 0.6 to >1.7 GPa and calculated temperatures are generally in the range of 700–900 °C (Table 4; Fig. 8A). In addition, four of the samples preserve evidence for multiple assemblages where portions of a P-T path can be estimated or inferred. First, isochemical phase assemblage modeling indicates that spinel stability for the ROB3 bulk composition is limited to pressures below ∼0.6 GPa, which extends to lower pressures with increasing temperature (Fig. 8A). The Fe-Mg exchange reaction between garnet and spinel is calculated at 800 °C with TWQ, although the compositions may have been modified during younger metamorphic events. Despite this uncertainty in the temperature of early metamorphism, ROB3 apparently underwent relatively low pressures early in its history prior to being taken to pressures as high as 1.3 Ga during peak metamorphism (Fig. 8B). An increase in pressure also predicts a transition from ilmenite to rutile stability, which is consistent with the former as inclusions in garnet and the abundance of rutile in the matrix.

An increase in pressure is also inferred for ROB9 from the destabilization of plagioclase indicated by the presence of this phase in inclusions in Cpx and its absence in the peak assemblage. An isochemical phase diagram for this bulk composition places the plagioclase-out boundary at ∼1.7 GPa at 800 °C (Fig. 8A). From TWQ calculations, garnet-pyroxene and garnet-biotite temperatures are near 800 °C whereas a pressure of 2.0 GPa at that temperature is calculated using the Jd + Qtz = Ab barometer and the albite composition from coexisting K-feldspar. This is considered a maximum pressure due to the absence of free quartz in the assemblage. Restabilization of plagioclase, and thus some likely decompression (Fig. 8B), is indicated by pyroxene margins that are recrystallized to Cpx3 + Pl3. Estimated conditions for this late stage are 730 °C and 1.37 GPa. A similar decompression stage in LSC6 is suggested by extensively recrystallized Cpx2 + Pl2 domains. Peak and retrograde conditions in this sample are estimated as 730 °C, 1.42 GPa and 720 °C, 1.18 GPa (Table 4; Fig. 8B). Extensive reaction between peak garnet and orthopyroxene in ROB6 (estimated conditions 850 °C, 1.2 GPa) to a new Opx2 + Sp assemblage (estimated as >900 °C and 0.6–0.9 GPa) suggests a reheating event accompanied by some decompression.

The calculated peak pressures are interpreted to approximate the depths from which the samples were derived during exhumation, with the exceptions of the three polymetamorphic samples with evidence for late-stage decompression described above. Retrograde pressures for samples with late decompression textures are interpreted to indicate residence depths before exhumation. The calculated pressures range from 0.6 to 1.53 GPa, which correspond to ∼23 to ∼54 km depths using the average crustal density profile of Christensen and Mooney (1995). These data are consistent with 8 of the 14 samples shown in Figure 9 derived from within the modern seismically defined 7.x layer as projected from the Deep Probe profile. The remaining samples are midcrustal and derived from depths shallower than the 7.x layer.

Bulk seismic velocities of each sample were calculated using the physical properties spreadsheet of Hacker and Abers (2004) to establish the context of the samples within the seismic observations. Velocities were calculated for the equilibration pressure at three temperatures: the calculated equilibrium temperature, 25 °C for comparison with laboratory measurements commonly reported at room temperature, and 500 °C, which is a baseline temperature for modern lower crust in this region (Fig. 9; Blackburn et al., 2012b). The latter perhaps provides a better comparison to modern seismic observations. Aside from uncertainties regarding the modern geothermal gradient, the most readily quantifiable source of error in the calculated seismic properties is due to variations in modal composition, which was evaluated by making similar calculations over halves of the full thin sections (Table 5; Fig. 9). Typical variations in the mode of a major phase are in the range of 5%–10% of the full section quantity, which results in errors of <1% on the compressional velocities (Table 5). The largest modal variations are for gneissic samples (e.g., to ±30% variation from full section garnet mode in ROB5), which result in 1%–2% variation in calculated velocities. Table 5 presents both averages calculated by the Hacker and Abers (2004) spreadsheet, the Voigt-Reuss-Hill (VRH) and Hashin-Shtrikman (H-S); the former does not account for the geometric arrangement of mineral grains, but has nonetheless been shown to be a remarkably robust scheme for evaluating elastic properties and is thus the most commonly used in Earth science literature (e.g., Mainprice and Humbert, 1994), whereas the latter explicitly assumes a statistically random structure (Bunge et al., 2000). As the anisotropy of physical properties is not considered here, further consideration of the microstructural arrangement of minerals is beyond the scope of this contribution. However, the results of Naus-Thijssen et al. (2011) suggest that variations from the calculated VRH velocities due to microstructure are not likely to exceed those described above based on modal mineralogy variations.

The results suggest a distinct increase in bulk seismic velocities with residence depths at ∼35 km and deeper and are consistent with the seismically determined high-velocity layer (Fig. 9). At these depths, calculated velocities (VRH average) at 500 °C vary considerably from 6.97 to 7.85 km/s (7.16–8.02 km/s at 25 °C). These data compare well with the 7.0–7.7 km/s velocities modeled for the lower crust in central Montana from the Deep Probe study (Gorman et al., 2002). The interpreted residence depths are also consistent with estimated crustal thicknesses from the Deep Probe study, but do not preclude somewhat thinner crust, as suggested by recent passive source studies (e.g., Bensen et al., 2009; Gilbert, 2012), particularly since the highest residence pressure calculated (for LSC04) is a maximum estimate due to the absence of quartz. The calculated velocities also significantly exceed the 6.40–7.00 km/s range measured for garnet-free lower crustal xenoliths (at 25 °C and 1.0 GPa) from the Leucite Hills in southernmost Wyoming (Farmer et al., 2005), which is consistent with the absence of an imaged 7.x layer in that part of the craton (Gorman et al., 2002; Rumpfhuber et al., 2009).

Monazite geochronology was undertaken for five samples in order to further constrain the origin and metamorphic history of the Robinson Ranch xenolith suite. Slow diffusion of Pb in monazite makes it a robust chronometer for high-temperature events (Parrish, 1990; Cherniak et al., 2004), and the sensitivity of monazite to recrystallization makes it particularly useful for elucidating the history of polymetamorphic rocks (Williams et al., 2007). However, this same sensitivity to recrystallization commonly results in complexly zoned grains that may require high spatial resolution. Fortunately, the generally high Th and radiogenic Pb content of monazite allows utility of the high spatial resolution capability of the Th–U–total Pb method, which is capable of analyzing compositional domains as small as 2 μm (Jercinovic et al., 2008; Mahan et al., 2010).

Analytical Methods

Identification of monazite and its textural setting were performed using an automated mapping routine using QEMSCAN at the Colorado School of Mines similar to that described herein for major phases. The routine is specifically calibrated to rapidly target monazite at a step size of 4 μm. The subsequent protocol for monazite geochronology followed that of Williams et al. (2006). Select grains were X-ray mapped at the University of Massachusetts-Amherst using a Cameca SX-50; U Mβ, Th Mα, Nd Lα, and Ca Kα were measured on PET crystals and Y Lα was measured on a TAP crystal with an accelerating voltage of 15 kV, a beam current of 200 nA, and a pixel dwell time of 80 ms. These X-ray maps were used to define compositional domains and direct the analytical strategy.

U–Th–total Pb monazite dates were produced on a modified Cameca SX-100 (Ultrachron) at the University of Massachusetts (UMASS). Appendix A of Dumond et al. (2008) detailed current analytical procedures and standard compositions used at the UMASS microprobe facility, including a list of all spectrometers, column conditions, count times, standards, and overlap corrections. In addition to the overlap corrections listed in Dumond et al. (2008), an overlap of Nd Lβ on Eu Lα was made. Acceptance or rejection of points follows the criteria of Dumond et al. (2008). For each domain, a weighted average of accepted analyses was calculated after Williams et al., 2006. The 2σ error reported for each date (Table 6; Table B1 [see footnote 1]) is the larger of either the error calculated by propagating the analytical uncertainty on trace element compositions through the age equation plus an estimated 1% uncertainty in background intensities (Williams et al., 2006) or two times the standard error of the mean. The consistency standard used was the Moacyr Brazilian pegmatite monazite, which has the following weighted mean isotope dilution–thermal ionization mass spectrometry (ID-TIMS) dates: 506.4 ± 1.0 (2σ, mean square of weighted deviates, MSWD = 0.6) for 208Pb/232Th, 506.7 ± 0.8 Ma (MSWD = 0.83) for 207Pb/235U, and 515.2 ± 0.6 Ma (MSWD = 0.36) for 206Pb/238U (W.J. Davis [Geological Survey of Canada], 2007, personal commun.).

Samples and Results

All xenoliths in this study except the high-Ca mafic granulites contain monazite. We report data from five samples: the low-Ca mafic granulite ROB3, two intermediate granulites ROB1 and ROB6, and the two felsic granulites ROB5 and ROB7. In all, 48 compositional domains in 33 grains from 5 samples were analyzed (Table 6; full chemistry for each domain is reported in Table B1 [see footnote 1]).

ROB1 contains rare monazite that only occurs in secondary textural settings along grain boundaries. The grains are small (10–50 μm), unzoned, and commonly have high aspect ratios (to 1:10; Figs. 10A, 10B). All three dated grains contain relatively little Th, U, and Pb (<2.0 wt% Th; <1000 ppm U and Pb) and yield imprecise but consistent ca. 1.3 Ga dates (1312 ± 99 Ma to 1279 ± 44 Ma; Table 6). U-Pb TIMS dates for three zircon grains from this sample were reported in Blackburn et al. (2012b); two grains are <1% discordant with 207Pb/206Pb dates of 1683 ± 2 Ma and 1699 ± 1.3 Ma, and the third is 1.3% discordant with a 207Pb/206Pb date of 1761 ± 1 Ma.

The low-Ca, high-Al mafic granulite ROB3 contains abundant monazite in a variety of textural settings. At least three populations are recognized based on Th, U, and Y composition as well as textural setting. First, low Th (<5000 ppm), low U (<500 ppm), and low Y (<600 ppm) grains (15 μm) typically occur in early textural settings. Two grains in the population are recognized; one is included in a spinel, which is in turn included in the core of a large garnet porphyroblast (Figs. 10C, 10D), and the second is included in plagioclase. Three domains from these grains yield imprecise dates ranging from 2170 ± 128 Ma to 2039 ± 154 Ma. The second population (5–50 μm) also tends to occur in relatively early textural settings such as inclusions within garnet, and contains low Th (<1.0 wt% Th) and Y (<600 ppm) but relatively high U (1500–3300 ppm). Thus, Th/U ratios in this population are distinctly low (1–3; Table 6). Four domains range from 1890 ± 55 Ma to 1809 ± 29 Ma. The third population (50–270 μm) is the most abundant and occurs both as inclusions in garnet and commonly in the matrix. This population is characterized by generally high but variable Th concentrations (1–8 wt% Th). There are 17 domains that range from 1829 ± 15 Ma to 1738 ± 14 Ma. The large range of dates and the observation that some grains within this population have distinct cores and rims with older and younger dates suggest episodic monazite growth over an interval of several tens of millions of years. For example, ROB3 m6 contains inner and outer core domains with dates of 1794 ± 20 Ma and 1798 ± 16 Ma, but an outer rim with a date of 1747 ± 40 Ma (Figs. 10E and 11). Using a distinct 23 m.y. gap in the dates among these 17 domains, we distinguish two groups. The older group of 13 domains yield a weighted mean of 1793 ± 11 Ma (MSWD = 4.2) and the younger group of 4 domains yields a weighted mean of 1748 ± 7 Ma (MSWD = 0.96). The high MSWD of the former group may indicate that it still comprises multiple monazite growth episodes that are indistinguishable without additional work.

ROB5 also contains evidence for multiple populations of monazite (25–200 μm). Similar to ROB3, the first population characteristically occurs in early textural settings such as inclusions in garnet and yields imprecise 2.0–2.2 Ga dates (Fig. 10H). However, the composition of these domains is highly variable and spans the full range of observed compositions in this sample. Two grains that occur in matrix quartz and feldspar yield similar dates of 1795 ± 18 Ma and 1785 ± 26 Ma (Fig. 10J). One grain, which occurs along a garnet-plagioclase grain boundary, has a sharply zoned core and rim that yielded distinctly younger dates of 1662 ± 18 Ma and 1691 ± 15 Ma, respectively.

ROB6 and ROB7 both contain abundant monazite (grain sizes from 10 to 250 μm), commonly with complex zoning patterns suggesting multiple episodes of dissolution and reprecipitation (Figs. 10G and 10I–10L). In ROB6, a common pattern is for monazite grains to have high-Th (∼6–7 wt% Th), low-Y (<1000 ppm) resorbed cores with distinctly higher Y (>5000 ppm) overgrowths. Such grains occur both in the matrix and as inclusions in orthopyroxene. Other rare texturally early domains (e.g., cores or grains included in garnet) contain either very low Th (<3000 ppm) or very high U and Y (both >2.5 wt%). All texturally early domains described here yielded dates between 1.79 and 1.78 Ga. A weighted mean of all three analyzed domains is 1784 ± 10 Ma (MSWD = 0.36). The texturally youngest monazite population is restricted to grains that touch or occur within late reaction zones containing Opx2 + Sp2 (Figs. 7B and 10G), suggesting that they grew in association with this reaction. These grains have small overgrowths (10–20 μm) that occur outside of the high Y overgrowths described here (Figs. 10I, 10K). They have similarly high Y (5000–8000 ppm) but distinctly lower U (<2200 ppm). All three analyzed domains give a weighted mean date of 1715 ± 35 Ma (MSWD = 1.3).

The felsic granulite ROB7 contains three compositional and textural populations (Fig. 10L), all of which occur in the matrix but are also locally included in garnet. The oldest population occurs as innermost resorbed cores with the highest Y (∼1.8 wt%) and U (>6000 ppm) and the lowest Th concentrations (∼6 wt%). One domain from this population was dated as 1776 ± 15 Ma; the second undated population with lower U and Y concentrations overgrows the former but is itself resorbed. The third and texturally youngest population occurs as rims on most grains and has the lowest Y and U (<3000 ppm and <2000 ppm, respectively) and the highest Th concentrations (to 10.6 wt%). Two of these rim domains give dates of 1754 ± 12 Ma and 1740 ± 15 Ma. Several monazite grains occur within late veins of Kfs, where they are resorbed and partially overgrown by apatite.

Heterogeneous High-Velocity Lower Crust

Heterogeneity is a notable characteristic of the entire suite of xenoliths studied from central Montana. The suite displays a diverse range of modal mineralogy and textural characteristics, bulk major element chemical compositions (although the deepest xenoliths have generally basaltic compositions), and bulk seismic velocities, even within the seismically defined 7.x layer. In addition, there is significant variation in metamorphic histories and geochronological signatures from the xenoliths. Multiple samples preserve evidence for prograde burial with at least one recording evidence for an early mid-crustal metamorphic event, and others record evidence for retrograde tectonic decompression and/or later reheating events. Along with monazite U–Th–total Pb geochronologic data from this study, published U-Pb data from zircon from xenoliths in this locality and others in the region (Davis et al., 1995; Scherer et al., 2000; Bolhar et al., 2007; Blackburn et al., 2011) indicate igneous, metamorphic, and/or fluid flow events before 2.6 Ga, and at 2.1 Ga, 1.83–1.68 Ga, and 1.5–1.3 Ga. Collectively, these data suggest that the 7.x layer beneath Montana is a composite feature.

Incremental Assembly of the Central Montana Lower Crustal 7.x Layer

We propose a history of incremental assembly of the high-velocity lower crust in this region from Archean to Mesoproterozoic time motivated by the heterogeneity in composition and chronology observed in xenolith samples (Fig. 12). We acknowledge that any model that specifies amounts or mechanisms for material added to the lower crust based on currently available data would be highly speculative. Instead, our intent here is to provide a discussion of what events and processes that may have resulted in the addition of material to the 7.x layer in this region through time. We present a compilation based upon regional studies of surface geology, geophysical observations, xenolith studies, and new data presented here.

Archean Components

Several mechanisms for Archean-aged components of the 7.x layer in Montana and Wyoming have been proposed. The generally small volume of exposed post-Archean magmatic rocks in the Wyoming craton has led some to suggest that the high-velocity lower crust has a significant if not an entirely Archean origin (Snelson et al., 1998; Chamberlain et al., 2003). Proposed mechanisms include a restite associated with widespread magmatism at 2.9–2.75 Ga in the Bighorn subprovince or an underplate during the formation of the 2.7 Ga Stillwater layered mafic intrusion (Chamberlain et al., 2003, and references therein). Zircon data from felsic xenoliths, broadly granitic or tonalitic in composition, from the Sweet Grass Hills indicate Archean protolith ages, although they have pressures suggesting derivation from above the 7.x layer in that region (Davis et al., 1995; Blackburn et al., 2011). Bolhar et al. (2007) used whole-rock and feldspar Pb isotopic compositions from Bearpaw Mountains crustal xenoliths to argue for 2.8–4.0 Ga protolith ages, although no depth information for the xenoliths is provided, and their model dates are not confirmed with zircon data. While no direct geochronological evidence for an Archean origin to components of the 7.x layer is yet published, it seems likely that some of the layer is Archean.

Rift-Related Underplating and Intraplating at 2.2–2.0 Ga

Major episodes of continental rifting during the 2.2–2.0 Ga interval are suggested for both the northern and southern margins of the Wyoming craton based on surface geology. In southwestern Montana, these include ca. 2.06 Ga mafic dikes (Mueller et al., 2004) and 2.1 Ga contact metamorphism associated with ultramafic intrusions (Alcock and Muller, 2010). Mafic dikes and sills of similar age also occur farther south in central (Harlan et al., 2003) and southern Wyoming (Premo and Van Schmus, 1989; Cox et al., 2000).

Montana xenoliths also record deep crustal thermal events that we associate with this rifting episode. Two of the xenoliths from this study contain texturally early 2.2–2.0 Ga populations of metamorphic monazite, whereas a crustal xenolith from nearby Biebinger Ranch in the Bearpaw Mountains contains zircon rims with a concordant 207Pb/206Pb date of 2167 ± 26 Ma (Bolhar et al., 2007). We suggest that this metamorphism may have been associated with the passage of rift-related mafic magmas through the crust and likely some degree of magmatic emplacement in the lower crust.

Possible Addition ca. 1.8–1.7 Ga During Great Falls Tectonic Zone Activity

A ca. 1.8–1.7 Ga high-velocity magmatic underplate beneath the Wyoming craton and Medicine Hat block was proposed (Gorman et al., 2002; Clowes et al., 2002; Eaton, 2006), primarily based on limited reports of geochronology on lower crustal xenoliths from the Sweet Grass Hills of southern Alberta and northern Montana (Davis et al., 1995). These data indicate multiple zircon growth events between 1.85 and 1.69 Ga. However, Davis et al. (1995) interpretation is primarily one of granulite facies reworking of lower crust during this interval rather than favoring a clear signal of new magmatic addition to the lower crust. This time frame also coincides with available constraints for the age of convergent tectonic activity along the Great Falls tectonic zone. Arc-related plutonism was ongoing at 1.86 Ga in central Montana (Mueller et al., 2002), and several ranges in southwestern Montana contain 1.82–1.71 Ga records of high-grade tectonism that are associated with the Great Falls tectonic zone (Roberts et al., 2002; Big Sky orogeny of Harms et al., 2004).

The dominant monazite populations in four of the five xenoliths in this study also reflect metamorphism during the interval 1.83–1.66 Ga. This is also the interval within which deformation-related foliation and gneissic layering in some xenoliths are likely to have developed. The 1.89–1.81 Ga population of monazite in ROB3 may coincide with arc-related magmatism described by Mueller et al. (2002). The most prominent populations in this and several other xenoliths are 1.79–1.78 Ga or 1.75 Ga and likely represent the timing of peak metamorphic events. Increases and late-stage decreases in pressure recorded by some xenoliths in this time frame, such as ROB3 and ROB6, respectively, reflect crustal thickening and thinning, possibly associated with Great Falls tectonic zone collision. Similar processes may also explain the record of up-pressure destabilization and post-peak stability of plagioclase in the hornblende eclogite, although no geochronological data are yet available for this sample. Additional published data from Bearpaw Mountains crustal xenoliths include 1.76–1.68 Ga zircon from mafic granulite ROB1 (same sample as in this study; Blackburn et al., 2012b), 1.85–1.75 Ga zircon from a felsic granulite (Bolhar et al., 2007), and 1.71 Ga zircon from a garnet pyroxenite (Scherer et al., 2000). Carlson and Irving (1994) reported 1.78 Ga monazite from mantle xenoliths from the Highwood Mountains, interpreted as evidence for fluid or melt interactions with the lithospheric mantle.

Multiple mechanisms for 1.8–1.7 Ga addition of high-velocity and high-density lower crust in the Wyoming craton and Medicine Hat block have been proposed (e.g., Gorman et al., 2002; Clowes et al., 2002; Chamberlain et al., 2003; Eaton, 2006); however, the xenolith record points to granulite facies and high-temperature eclogite facies metamorphism of lower crust. It remains unclear whether a significant volume of new magmatic material was added at that time. Postorogenic collapse may have been accompanied by extension-related mafic magmas that could have added some new material to the lower crust and provided local sources for reheating of older host material. Chamberlain et al. (2003) suggested that a component might reflect mechanically imbricated oceanic material added during suturing. All of these ideas, and likely others unstated, are speculative until future work provides additional information.

Rift-Related Addition ca. 1.5–1.3 Ga

A previously proposed Mesoproterozoic origin for the 7.x layer is based on the widespread occurrence of 1.47–1.45 Ga mafic dikes across the Wyoming craton that are interpreted to reflect an eastern extension of the Belt Basin rift system (Chamberlain et al., 2003). Deposition of Belt Basin sediments, which extend in the subsurface to as far east as 107.5°W in central Montana (Fig. 3), and associated magmatism were episodic, with major pulses of rifting at 1470–1400 Ma (Zartman et al., 1982; Höy, 1989; Anderson and Davis, 1995; Sears et al., 1998; Evans et al., 2000) and ca. 1370 Ma (Doughty and Chamberlain, 1996). Postdepositional magmatism and associated contact metamorphism of sediments that were probably part of the rifting cycle are documented to be as young as 1320 Ma (Schandl et al., 1993; Zartman and Smith, 1995).

Crustal xenoliths from central Montana appear to also record cryptic deep crustal thermal and/or fluid-related processes that may be associated with the end of the Belt Basin rift formation. Monazite grains from the mafic granulite xenolith ROB1 yield dates of ca. 1.3 Ga and have textural settings and morphologies that indicate growth after peak metamorphism. Rutile that initially grew as part of the primary assemblage in a range of deep crustal xenoliths from the Bearpaw Mountains (including several xenoliths from this study) records the onset of Pb retention just after 1.50 Ga (Blackburn et al., 2012a). A garnet-clinopyroxenite xenolith exhumed in the northern Bearpaw Mountain field yielded Sm-Nd whole-rock and multimineral dates of 1.36–1.35 Ga, and zircon U-Pb intercepts of 1707 ± 15 and 1268 ± 34 Ma (MSWD = 0.49) (Scherer et al., 2000). The Sm-Nd dates indicate either the time the sample cooled through closure, which would not be inconsistent with cooling from 1.47 to 1.45 Ga magmatism, or the most recent time at which the Sm-Nd system was reset.

Magmatic addition to the crust at 1.47–1.45 Ga is clear from shallow geologic observation, and a similar magmatic contribution to the 7.x layer seems reasonable (Chamberlain et al., 2003). The ca. 1.3 Ga chemical and isotopic signals may have developed in response to the addition of heat and/or fluids from a magmatic underplate near the end of Belt Basin rifting, but the lack of evidence for near-surface magmatism at that time makes this possibility less compelling. Although the nature of this disturbance is unclear, the vein- and fracture-filling morphology of the ROB1 monazite suggests a fluid circulation event.

Simplified Consideration of Incremental Assembly

The high-velocity lower crust at the northern margin of the Wyoming craton makes up almost half of the ∼55-km-thick crust (Gorman et al., 2002). This is a much greater proportion than that seen in averaged worldwide compilations of shields and platforms (Fig. 1; Christensen and Mooney, 1995). Using a simple plane strain calculation, we tracked the thickness of a hypothetical column of initially 35-km-thick crust that has a generalized history similar to that described in the preceding discussion. During two rift-related thinning and magmatic underplating episodes (at 2.1 Ga and 1.4 Ga) the crust is thinned by 10 km, and 10 km of mafic lower crust is added. This is generally consistent with observations from the Rio Grande Rift (West et al., 2004; Wilson et al., 2005) and the main Ethiopian Rift (Cornwell et al., 2010). During one intervening 1.8 Ga collisional episode, the crust is homogeneously thickened by 20 km. The result is that an initial normal thickness (35 km) column of crust that starts with a 5 km (Archean) mafic lower crustal layer thickens to 55 km, the lower half of which is a mafic layer similar to that observed in the Wyoming craton. This calculation does not include the effects of other processes that likely occurred (e.g., magmatic differentiation, mechanical underplating, Cenozoic evolution), and the effects of the processes considered likely varied spatially (e.g., Chamberlain et al., 2003). Nevertheless, we conclude that incremental assembly is a reasonable way to achieve the modern crustal configuration considering the current state of knowledge of the regional geologic history.

Petrological, geochronological, and bulk seismic velocity data for crustal xenoliths from central Montana provide new insight into the nature and modes of formation of the lower crust in this region. Similar to some other past studies of deep crust from xenoliths (Rudnick and Taylor, 1987) and exhumed lower crust (e.g., Williams and Hanmer, 2006), we emphasize the heterogeneity in geologic history and physical properties of the deep crust from this xenolith record. We suggest that the modes of formation of the high seismic velocity lower crustal layer were also likely heterogeneous. From our xenolith data, earlier suggestions for the origins of the lower crust in this region, and considering the regional tectonic history, we present a model for incremental assembly of the 7.x layer involving episodic magmatic and possibly mechanical underplating and intraplating (and accompanying fractionation processes) associated with multiple regional tectonic events from Archean to Mesoproterozoic time.

The thickness of the high-velocity lower crustal layer in this region is significantly greater than most similar layers observed elsewhere. However, the observed crustal structure appears consistent with a history of incremental assembly based on the surface and xenolith observations available. In this context, perhaps it is the processes that have allowed preservation of the thickness of this lower crustal layer rather than the thickness that is anomalous. This study provides an example of how integration of xenolith data, surface observations, and geophysical studies can elucidate the formation, evolution, and present-day structure of the continental lithosphere.

We thank M.L. Williams and M.J. Jercinovic at the University of Massachusetts Ultrachron facility, J. Drexler at the University of Colorado microprobe facility, and S. Appleby at the Colorado School of Mines QEMSCAN facility. We thank C. Gerbi and D. Schutt for constructive reviews that improved the manuscript. The research presented here was funded in part by a Geological Society of America Graduate Student Research Grant to Barnhart and National Science Foundation EarthScope grants EAR-07464246 to Mahan and EAR-0746205 to Bowring.

1Supplemental File. Excel file of seven supplemental mineral composition tables and one monazite chemistry table. If you are viewing the PDF of this paper or reading it offline, please visit or the full-text article on to view the Supplemental File.
2We use a numeral as a subscript to a mineral abbreviation to denote different generations of mineral growth in polymetamorphic samples. Thus, an early generation of plagioclase associated with the earliest observed metamorphic event M1 is denoted as Pl1 and a subsequent generation of plagioclase is denoted Pl2.
3We use X with a subscript to denote the proportion of mineral end-member components. For example, XMg = 0.25 indicates that 25% of the mineral is the Mg-endmember of the solid-solution. The subscript may also be an abbreviation for the named end-member mineral (e.g., XGrs = 0.25 indicates that 25% of the garnet composition is grossular, which is the Ca-endmember of the series).