Geophysical extent of the Wyoming Province, western USA: Insights into ancient subduction and craton stability

Defining the boundaries and structural framework of blocks of crustal lithosphere of varying age is fundamental to understanding numerous geological processes, from the onset of plate tectonics and the tectonic evolution of the continents to the distribution of mineral resources and the interpretation of xenolith data. In terms of fundamental Earth processes, the formation of the Archean Wyoming Province (WP) and the orogens along its margins span a time window from the Paleoarchean to the Paleoproterozoic, during which time the continental nuclei composing Laurentia formed and assembled. This time period also is marked by major transitions in Earth history, including changes in ocean and atmosphere chemistry, the rise of biogenic carbon, increased continental weathering, and the onset of modern plate tectonics (e.g., Kump et al., 2011; Bindeman et al., 2018; Hawkesworth et al., 2020; Swanner et al., 2020). Thus, the WP provides an opportunity to investigate Earth evolution and plate tectonics through this critical time window.

Save for exposures along the Cheyenne belt in southern Wyoming, USA (Fig. 1), the boundaries of the Archean WP are concealed (Mueller and Frost, 2006). In addition, younger tectonics, in particular Basin and Range extension and Cenozoic magmatism along the Snake River Plain (SRP), overprint the western province margin. Recent studies have called into question the eastern margin of the WP (e.g., Worthington et al., 2016; Kilian et al., 2016; Hrncir et al., 2017; Chamberlain and Mueller, 2019) while questions remain as to the stability and degree of lithospheric modification the WP lithosphere has undergone (Kusky et al., 2014; Dave and Li, 2016; Snyder et al., 2017; Bedle et al., 2021).

In this paper we present a 3-D electrical resistivity model that provides new constraints on the extent of the WP, its surrounding Proterozoic orogens, and the internal architecture interpreted to reflect its protracted Archean assembly. Together with a new geochronologic compilation and incorporating independent geological and geophysical constraints, we refine the margins of the WP, explore changes in subduction from the Neoarchean to the Paleoproterozoic in a global context, and discuss the stability of the WP in light of recent and ongoing tectonism.

Originally called the “Wyoming Archean province” (Condie, 1976), the Wyoming Province is defined as “the region in Wyoming and adjacent states underlain by rocks of Archean age” (Houston et al., 1993). The WP is the most southwesterly (present-day coordinates) of several large Archean blocks within the core of North America (Fig. 1, inset) and is located south of the Medicine Hat block (MHB) and Hearne Province and west of the Superior Province (SP). Although not explicit in the original definition, the boundaries of the WP have come to be associated with Proterozoic orogens (Foster et al., 2006; Mueller and Frost, 2006). Because a majority of the intervening basement is covered by younger, mainly Phanerozoic, sedimentary rocks, crust within the limits of these bounding orogens is assumed to be contiguous and mostly, if not exclusively, Archean. The major bounding Paleoproterozoic orogens are the 1.78–1.74 Ga Cheyenne belt along the south, the ca. 1.7 Ga Farmington zone to the west, the 1.86–1.72 Ga Great Falls tectonic zone to the north, and the 1.79–1.72 Ga southern Trans-Hudson orogen to the east (Fig. 1) (Barnett et al., 1993; Chamberlain, 1998; Dahl et al., 1999; Harms et al., 2004; Mueller et al., 2005). These bounding orogens include both reworked Archean crust and juvenile Paleoproterozoic crust (Ball and Farmer, 1991; Dahl et al., 1999; Foster et al., 2006; Condit et al., 2015; Mueller et al., 2016). Precambrian exposures are limited to the cores of Late Cretaceous to Eocene basement-involved uplifts of the Laramide orogeny (<10% of the area of underlying Archean crust) and the province is overprinted by extensive Cenozoic volcanism in the northwest, making determination of the precise extent of the province difficult. The boundaries suggested by Foster et al. (2006), shown as green dashed lines on Figure 1, are widely but not universally accepted. For example, extrapolations of the southern boundary west of the Cheyenne belt include those that run along the Wyoming-Colorado state line (Mueller et al., 2011) and those that place the southern boundary farther north (Nelson et al., 2011). In the north, the boundary is associated with the Great Falls tectonic zone, a 150–250 km wide, NE-SW trending zone originally defined on the basis of high-angle faults and lineations in sedimentary cover as well as magnetic and gravity anomalies (O’Neill and Lopez, 1985). However, interpretation of the precise location and width of this zone varies (c.f. Holm and Schneider, 2002; Carlson et al., 2004; Gu et al., 2018). The location of the eastern boundary, typically placed to include the Black Hills, has alternatively been suggested to lie some 250 km to the west near the Bighorn Mountains based on seismic and magnetic data (Worthington et al., 2016; blue dashed line on Fig. 1). Additional geophysical data are needed to resolve uncertainties in the placement of province boundaries.

The geologic history of the Archean Wyoming province is interpreted from studies of the limited exposures provided by basement-cored uplifts in Wyoming, Montana, and South Dakota (USA). Combined geochemical, petrologic, and structural studies have led to the identification of three sub-provinces: the Montana metasedimentary terrane (MMT), the Beartooth-Bighorn magmatic zone (BBMZ), and the southern accreted terranes (SAT) (Mogk et al., 1992; Mueller and Frost, 2006) (Fig. 1). The limited exposures within these sub-provinces make determination of their extents and boundaries problematic. Nonetheless, discrete histories have become apparent through observations and interpretations, as described below.

The Montana metasedimentary terrane (MMT) is confined to the northwestern corner of the WP and is exposed in the Tobacco Root Mountains and Ruby, Gravelly, Gallatin, and Madison ranges of southern Montana (Mogk et al., 1992). The subprovince is composed of tonalitic and trondhjemitic gneisses ranging in age from ca. 3.5 to 2.7 Ga that are intercalated with metasupracrustal rocks including orthoquartzite, iron formation, marble and calcareous gneiss, pelitic to mafic schist, as well as mafic and ultramafic rocks (e.g., Mogk et al., 2004; Mueller et al., 2004). The greenschist- to amphibolite-facies metasupracrustal rocks suggest a miogeoclinal or passive-margin environment; it is not clear that they constitute a single sequence. Sparse volcanic rocks and cross-cutting intrusive rocks combined with the youngest detrital zircons suggest deposition no earlier than mid-Neoarchean for much of the supracrustal assemblage(s) (Mueller et al., 1993).

The BBMZ is dominated by Meso- and Neoarchean granitoids and gneisses in the Beartooth, Bighorn, Wind River, and Owl Creek mountains of Wyoming (Frost et al., 2006a; Mueller et al., 2010), with older, mainly metasupracrustal rocks, located in the eastern Beartooth Mountains of Wyoming and Montana (Henry and Mueller, 1982; Mueller et al., 2014) and Paleoarchean orthogneisses in the Granite Mountains of central Wyoming (Frost et al., 2017). The ca. 2.85 Ga Barlow Springs Formation of Langstaff (1995) in the southern BBMZ, composed of quartzite, banded iron formation, metapelite, and amphibolite, is interpreted as a passive-margin sequence. The boundary between the MMT and BBMZ is inferred to lie north of both the Teton Range and the South Snowy block of the Beartooth Mountains. In the Beartooth Mountains, the boundary zone lies within the North Snowy block and Yankee Jim Canyon of the Beartooth Mountains (Mogk et al., 1992). Although the MMT and BBMZ are defined by different lithologic associations, different age modes defined by magmatic and detrital zircons, and distinct metamorphic, depositional, and magmatic histories, they share a distinctive geochemical trait: evolution in a long-term high U/Pb environment (Wooden and Mueller, 1988; Frost et al., 1998). Juxtaposition of the MMT and BBMZ in the Neoarchean is suggested by a ca. 2.55 Ga leucogranite stitching pluton that intrudes a major NE-SW-trending shear zone between the MMT and BBMZ that is exposed in the North Snowy Block of the Beartooth Mountains (Mogk et al., 1988, 1992).

The SAT include blocks of Neoarchean supracrustal rocks dominated by juvenile graywacke and mafic rocks, such as at South Pass in the southern Wind River Range, the Rattlesnake Hills and Bradley Peak in central Wyoming, and the northern Sierra Madre of southern Wyoming (Frost et al., 2006b; Souders and Frost, 2006). These blocks are interpreted to have accreted to the BBMZ and deformed prior to the intrusion of the undeformed ca. 2.63–2.62 Ga Louis Lake and Wyoming batholiths in the Wind River Range and Granite and Laramie Mountains (Frost et al., 1998; Bagdonas et al., 2016). The boundary between the BBMZ and the SAT is exposed only in the South Pass area, where the juvenile, Neoarchean Miners Delight Formation (Frost et al., 2006b) was thrust onto ca. 3.0 Ga quartzofeldspathic gneiss (new U-Pb zircon data presented in online materials) and remnants of a passive-margin sequence at ca. 2.65–2.63 Ga (Frost et al., 2006b). Farther east in central Wyoming, the boundary is obscured by intrusion of the 2.62 Ga Wyoming batholith, which extends across the Granite, Pedro, Shirley, and Laramie mountains (Bagdonas et al., 2016). Because the Nd isotopic composition of the BBMZ is distinctly less radiogenic than the more juvenile southern accreted terranes, Nd isotopic ratios can be used to estimate the location of the boundary (Frost et al., 2006b). BBMZ basement and evolved metasedimentary rocks deposited on that basement have ε_Nd values at 2.6 Ga of <–4 (Frost et al., 2017; Frost and Da Prat, 2021), therefore the southernmost samples with ε_Nd <–4 may approximate the limits of the BBMZ. The “ε_Nd = –4” corridor runs through the Granite Mountains and across the southern Shirley Mountains to the Laramie Mountains (Fig. 1).

Geophysical data provide constraints on the crust and upper mantle beneath the WP. Domains of positive and negative aeromagnetic anomalies form broad arcs extending across northern Wyoming and southeast Montana (Fig. 2A); the continuity of these geophysical anomalies across areas of little to no basement exposure provides some of the strongest evidence for the interpreted extent and continuity of the BBMZ. In areas where bedrock is exposed, granitic rocks correspond to positive anomalies while gneisses and supracrustal rocks are associated with magnetic lows (Sims et al., 2001). Ground-based gravity data over the WP (Fig. 2B) are dominated by a series of intense highs and lows associated with basement-cored uplifts and basins, respectively. The gravity highs in particular reveal the presence of several concealed basement uplifts, such as a prominent high northwest of the Black Hills uplift.

The Deep Probe seismic refraction/wide-angle reflection experiment imaged the crust below a north-south profile (Fig. S1¹) crossing the WP from the western Bighorn Basin to the Granite Mountains and from the southern Wind River Mountains into Colorado (Henstock et al., 1998; Snelson et al., 1998; Gorman et al., 2002). Their models show an interpreted step in Moho depth as well as a >20-km-thick lower-crustal high-velocity layer beneath the WP to the north of the Granite Mountains. This layer is absent farther south beneath the SAT, a finding in accordance with seismic imaging beneath the Sierra Madre Mountains (Morozova et al., 2005). A change in upper-mantle velocity also occurs north and south of the Granite Mountains (Snelson et al., 1998). The thick lower crust is characterized by high velocities (V_P) of 7.0–7.3 km/s and densities >3000 kg/ m³, suggesting a mafic composition. Gravity modeling along a north-south profile through the Powder River Basin further suggests a high-density, lower-crustal layer is present beneath northeast Wyoming (Snelson et al., 1998). However, the Bighorn Arch Seismic Experiment found the 7.x layer to be absent or discontinuous beneath the Bighorn uplift (Worthington et al., 2016) yet present both east and west of it. A lower-crustal layer with distinctly higher velocities than the WP is imaged beneath the MHB and has been interpreted in terms of Proterozoic underplating following accretion of the WP to the MHB (Gorman et al., 2002).

Archean cratons are characterized by their thick, buoyant, and viscous lithospheric roots and high relative integrated yield strength (Bedle et al., 2021). Global tomographic models, including CAM2016 (Priestley et al., 2018) and SGlobe_rani (Chang et al., 2015) suggest that WP lithosphere is thinner than many other Archean cratons. This observation has led to suggestions that the WP mantle lithosphere has been destabilized and destroyed subsequent to its formation (Kusky et al., 2014; Dave and Li, 2016; Snyder et al., 2017). Several events may have modified and weakened the cratonic lithosphere by introduction of fluids, heat, and/ or magma. A Paleoproterozoic metasomatic event has been identified within the Eocene Absaroka Volcanic Supergroup in northwestern Wyoming (Feeley, 2003) and from mantle peridotite xenoliths from the diatreme field in the Missouri Breaks area in central Montana (Carlson et al., 2004). A Late Cretaceous–Paleocene metasomatic event attributed to shallow subduction of the Farallon plate affected Absaroka volcanic rocks and peridotites and peridotite xenoliths from the Grassrange diatreme field (Feeley, 2003; Carlson et al., 2004). Late Cretaceous shallowing of the Farallon plate may also have eroded and displaced the base of the WP below ~150 km (Snyder et al., 2017). In addition, the intrusion of magma associated with the Yellowstone hotspot may modify and erode the cratonic lithosphere (Kusky et al., 2014; Dave and Li, 2016). Relative to the rest of the WP, this northwest corner is geophysically distinct, with reduced seismic velocity (e.g., Yang et al., 2011; Porritt et al., 2014) (Figs. 2D and 2E), elevated conductivity (Kelbert et al., 2012; Meqbel et al., 2014), and elevated heat flow (Fig. 2C) (Blackwell et al., 2011) in addition to changes in the amplitude, trend, and wavelength of potential-field anomalies (Figs. 2A and 2B). Dave and Li (2016) have further suggested that small-scale convection has modified the cratonic root beneath the eastern WP based upon a region of low shear-wave velocity below 150 km depth (centered upon 45.5°N, 105.5°W); this anomaly is not present, however, in other regional seismic models.

More detailed information about the lithospheric structure of the WP is provided by regional tomographic models, such as the Western United States (WUS) model of Schmandt and Humphreys (2010a, 2010b). Bedle et al. (2021) observe that a number of geologic features correlate with variations in compressional- and shearwave velocities imaged by this model, including seismically slow lithosphere trending NE-SW along the Snake River Plain-Yellowstone hotspot track, and slow velocities associated with the Basin and Range and Rio Grande Rift (Figs. 2D and 2E). They also note a difference based on dV_S data in the thickness of cratonic lithosphere beneath the northern and southern parts of the WP, with thicknesses of >250 km beneath the BBMZ, and thinner, ~200-km-thick cratonic lithosphere beneath the SAT.

A 3-D resistivity model has been estimated from inversion of full magnetotelluric (MT) impedance and vertical magnetic-field transfer functions across the western United States (west of 101°W longitude). Data from 1525 MT sites were inverted using the ModEM inversion code (Egbert and Kelbert, 2012; Kelbert et al., 2014). This paper focuses on the northeast quarter of the model area, a region covered by ~500 MT stations (Fig. S1).

Inverted data consist of long-period soundings collected during the EarthScope program (Kerr, 2013) and more than a dozen local data sets. All data are publicly available through the Incorporated Research Institutions for Seismology Data Management Center (Kelbert et al., 2011). Both full impedance (Z) and vertical magnetic-field (T) transfer functions (tippers) were inverted for the combined data set at 19 periods from 5 to 20,000 seconds. The model mesh has a uniform 10 km horizontal cell size; the vertical mesh is non-uniform, with logarithmically spaced cell thickness (0.04 in log units) starting from 25 m at the surface. The 100 Ωm start model contains over 3.9 million model cells. The electrically conductive Pacific Ocean, though far from the model areas of interest, does have an impact on long-period Z and T responses up tô1000 km from the coastline. We have included the conductive ocean a priori in the start model using 1 arc-minute resolution bathymetry (NOAA National Geophysical Data Center, 1985). Data were inverted with statistically determined errors but subject to error floors defined as 5% of the square root of |Zxy·Zyx| and 0.03 in T.

A sequential inversion approach (Bedrosian et al., 2018) was applied to balance the fit between impedance and tipper data and to build structure progressively within the model while simultaneously reducing error floors. A 100 Ωm prior model was applied in all inversion steps and serves to damp unnecessary structure in those parts of the model domain with little data sensitivity. A final normalized root-mean-square (nRMS) data misfit of 3.96 was obtained, representing a 90% reduction in misfit relative to the starting halfspace model (nRMS = 38.27). The final resistivity model is available in netCDF format from the IRIS Earth Model Collaboration (Bedrosian, 2021).

In map view (Fig. S1), impedance data misfits are uniformly distributed, although elevated misfits are concentrated in areas of dense station coverage, where resistivity variations occur over short length scales that cannot be adequately modeled by the 10 km model cell size. Tipper data are, with few exceptions, well fit across the modeling area (Fig. S2), reflecting the sensitivity of these data to longer-wavelength structure. Examination of data misfit by station, period and component (Fig. S3) reveals a uniform fit with no systematic patterns in misfit, save for a general increase in misfit at shorter periods, again reflecting short-wavelength resistivity variations that cannot be reproduced in this regional model. The diffusive nature of MT fields in the Earth serves to degrade resolution with increasing depth; however, a minimum horizontal resolution of 30 km is estimated based upon the applied model covariance (which defines the degree of smoothing regularization) and the horizontal cell size.

Published U-Pb zircon, monazite, titanite, and baddeleyite age determinations of igneous and metamorphic rocks from the WP and adjacent areas were compiled to aid in determining the location of Archean province boundaries. In addition to a comprehensive list of published Archean U-Pb dates from the WP, the compilation also includes age determinations from the MHB and Paleoproterozoic ages determined for rocks from WP margins and its bounding orogens. Seven U-Pb dates on Paleoproterozoic mafic dike samples from within the province are also compiled. Because no published age determinations were available for basement gneisses at the southern end of the Wind River Range, we report a new U-Pb zircon age for quartzofeldspathic gneiss sample 04SP-QF of 3014 ± 2 Ma. These data, along with a spreadsheet and shapefile for the compilation of 243 U-Pb age determinations, are available in the Supplemental Material (see ^{footnote 1}).

Several additional data sources are considered as part of our investigation. Magnetic data (Fig. 2A) are extracted from a North American compilation analytically continued to a height of 305 m (Bankey et al., 2002; https://mrdata.usgs.gov/magnetic/). A merged gravity grid (Fig. 2B) was corrected for crustal thickness assuming an Airy-Heiskanen isostatic model (https://mrdata.usgs.gov/gravity/ and http://gdr.agg.nrcan.gc.ca/gdrdap/dap/index-eng.php). The heat flow map (Fig. 2C) is from (Blackwell et al., 2011). The body-wave seismic tomography models (Figs. 2D and 2E) were downloaded from the IRIS Earth Model Collaboration (Schmandt and Humphreys, 2010b).

Within the resistivity model, upper-crustal structure (Fig. 3A; 2 km depth) is dominated by the conductive sedimentary section of the Western Interior Seaway in addition to a series of conductive basins and resistive uplifts formed during the Laramide orogeny and subsequent Basin and Range extension (e.g., Teton Range). Comparison of the resistivity model to areas of Precambrian outcrop suggests that uplifted areas extend well beyond their mapped exposure. Several resistive uplifts are additionally offset from Precambrian outcrop. The Bighorns uplift, for example, is shifted east relative to an imaged resistor, consistent with NE-vergence of the master thrust on the eastern margin of the uplift (Sims et al., 2001). A number of high resistivity areas are located where Precambrian outcrop is absent. For example, a resistive region located northwest of the Black Hills uplift most likely reflects an Archean granitic uplift concealed beneath the eastern Powder River basin, consistent with gravity and magnetic highs (Figs. 2A and 2B), high P-wave velocities (Worthington et al., 2016) and a pronounced shallowing of the basin at this location (Blackstone, 1993). Whereas the Laramide uplifts are known to be thick-skinned and bounded by faults interpreted to extend >25 km (Lynn et al., 1983; Yonkee and Weil, 2015), there appears to be little expression of them in the resistivity model below ~10 km depth (Fig. 3B). We interpret this pattern to reflect a lack of resistivity contrast across these deep-seated faults.

The tectonic framework of the region is more clearly imaged at mid- to lower-crustal depths (Figs. 3B and 2C). At these depths, a first-order contrast in resistivity (Fig. S1) defines the boundary between stable North America and the active magmatic and tectonic domains of the western United States (Meqbel et al., 2014; Bedrosian and Feucht, 2014). The latter are characterized by elevated conductivity interpreted to reflect small amounts of lower-crustal melt and saline fluids beneath the Basin and Range (Wannamaker et al., 2008), the Yellowstone-Snake River Plain (Kelbert et al., 2012), and the Rio Grande Rift (Feucht et al., 2017, 2019). These three regions are further marked by an abundance of geothermal areas (shown on Fig. 3), reduced crustal- and upper-mantle seismic velocities (e.g., Schmandt and Humphreys, 2010a), and elevated heat flow values (e.g., Blackwell et al., 2011). In contrast, the geochemically depleted WP is imaged as a more coherent resistive domain, particularly at lower-crustal and upper-mantle depths (Figs. 3C and 4). The WP is characterized by faster seismic velocities, lower heat flow values, a general absence of geothermal prospects, and an abundance of Archean U-Pb zircon ages (Figs. 2C–2E).

A series of high-conductivity belts (HCBs) mark the northern, eastern, and southern margins of the WP (Figs. 3C and 3D). Their geometry, depth, and peak conductivity values (>1 S/m) are consistent with interconnected graphite and sulfides in metasedimentary rocks that are commonly emplaced during the subduction of carbonaceous passive-margin sediments that were deposited within an anoxic water column (Boerner et al., 1996; Bedrosian and Finn, 2021). This genetic model is discussed in greater detail in the discussion. For the purpose of what follows, these HCBs are interpreted to mark locations of past subduction and, as such, are considered proxies for the margins of the WP.

The ca. 1.78–1.74 Ga Cheyenne belt (CB) is characterized in the resistivity model as a prominent mid- to lower-crustal HCB (Figs. 3B–3D and 4) coincident with a magnetic low and offset to the southeast of a linear gravity high (Figs. 2A and 2B). This HCB is largely confined to >15 km depth but is imaged as shallow as 7 km adjacent to the Hartville uplift and between the Uinta and Sierra Madre mountains. The conductor coincides with Paleoproterozoic metasedimentary rocks mapped in both the Medicine Bow Mountains and the Hartville uplift of southern Wyoming (Karlstrom et al., 1983). There is a sharp transition from Archean to Paleoproterozoic U-Pb age determinations within the Medicine Bow Mountains, with the conductor located beneath the Paleoproterozoic rocks and with no Archean ages within or south of it. Both Paleoproterozoic and Archean age determinations within the Hartville uplift occur in close proximity to one another and are located along the edge of the HCB (Figs. 3B–3D). The CB conductor is less prominent to the west of 108°W but is interpreted to project south of the Uinta Mountains as far west as ~112°W, beyond which its conductive signature is overprinted by a broad region of elevated crustal conductivity beneath the Basin and Range province (Fig. 3B–3D).

Tectonic models of the CB are unresolved, with geological and seismic studies arguing for southeast (Karlstrom and Houston, 1984; Chamberlain, 1998; Tyson et al., 2002) and northwest- (Jones et al., 2010) directed subduction, sinistral transpression (Sullivan and Beane, 2013), the wedging of Paleoproterozoic crust within the WP (Morozova et al., 2002), and the imbrication of Paleoproterozoic crust beneath the southern edge of the WP (Hansen and Dueker, 2009). The resistivity model does not resolve this question, as there is no discernable dip on the lowercrustal HCB. However, its location beneath the Paleoproterozoic Green Mountain arc and to the southeast of uplifted Archean crust (suggested in part by the linear gravity high to the northwest) is most consistent with models invoking southeast-directed subduction of continental passive-margin sediments accompanied by uplift and deformation of the Archean province margin immediately to the northwest (Patel et al., 1999; Chamberlain et al., 2003).

The Trans-Hudson orogen (THO) is marked by two parallel HCBs separated by a 200-km thick section of Paleoproterozoic crust (Fig. 4) (Bedrosian and Finn, 2021). Based on prior modeling, the southern THO has been interpreted to be a composite margin, with a western belt delineating the edge of the WP and the eastern belt marking the edge of the Superior Province (SP) (Bedrosian and Finn, 2021). Only the westernmost belt falls within the domain of the resistivity model presented here.

The western HCB skirts the eastern edge of the Black Hills uplift and trends due north (Figs. 3B and 2C), coincident with several gravity highs and a magnetic quiet zone along much of its length (Figs. 2A and 2B). At ~46°N, this HCB is split into a western branch which ends abruptly at 49°N, and an eastern belt that continues north into Canada (Bedrosian and Finn, 2021). The latter can be traced geophysically for more than 1500 km from the Black Hills to 56°N; this conductive lineament coincides with the North American Central Plains conductor that was roughly located using magnetic-variometer arrays (Camfield et al., 1971; Alabi et al., 1975). Similar to the CB, the western HCB extends deep (15–40 km) within the crust, has peak conductivities in excess of 1 S/m, and is <50 km wide. As described in Bedrosian and Finn (2021), this HCB lies down-dip of a package of west-dipping seismic reflectors at 48.5°N (Nelson et al., 1993) which together are interpreted to mark west-directed subduction of juvenile Paleoproterozoic arc terranes against the Archean WP. The eastern HCB (not shown) is similarly correlated with east-dipping reflectors and taken to mark east-directed subduction beneath the SC (Bedrosian and Finn, 2021). Given the lack of crystalline basement exposure, there are few age determinations along the THO. However, Dahl et al. (1999) bracket deformation and collisional magmatism between ca. 1.79 and 1.72 Ga in the Black Hills. The uplift preserves both Archean and Proterozoic magmatic ages, consistent with its position at the edge of the HCB. Proterozoic orogenic gold deposits in the Black Hills lie just west of the HCB (Figs. 3C and 3D), also consistent with their formation during Paleoproterozoic orogenesis above a west-directed subduction zone (Bedrosian and Finn, 2021).

The northern margin of the WP is marked by a linear HCB that can be traced for more than 600 km from the Tobacco Root Mountains in the west tô105°W where it appears cut by the conductor along the THO (Fig. 3D). This HCB is coincident with a linear low-velocity zone imaged in the upper mantle (Fig. 2D) (Schmandt and Humphreys, 2010a). Although somewhat less conductive (~0.1 S/m) than the belts marking the THO and CB, this belt is of similar depth and width to both. It is tempting to attribute this conductor to the Great Falls tectonic zone (O’Neill and Lopez, 1985; Boerner et al., 1998; Gifford et al., 2018), however, its location is nearly 150 km southeast of both mapped exposures of the Great Falls tectonic zone (GFTZ) and potential-field lineaments to the northeast (Figs. 2A and 2B) that have been interpreted to mark its extension beneath cover. A sub-parallel HCB to the north, however, does overlap with the GFTZ as originally mapped (O’Neill and Lopez, 1985) and extends northeast almost to the Little Rocky Mountains, which are Archean in age (gray solid line on Fig. 1). A resistive region centered at 48°N, 108°W (Fig. 3C) appears to interrupt the northern HCB, which in contrast to its southern counterpart, does not extend as far east as the THO. The Archean age of the Little Rocky Mountains and their isotopic similarity to rocks of the MHB (Gifford et al., 2018), lead us to suggest that the resistive domain in northeast Montana is a promontory of the MHB extending as far south as 47°N.

The age, extent, and tectonic significance of the GFTZ is poorly understood due to concealing Phanerozoic cover, with Precambrian rocks limited to exposures in the Little Belt and Little Rocky Mountains and crustal xenoliths brought up in Eocene volcanic rocks in the Sweetgrass Hills, Missouri Breaks, and Grassrange areas (Fig. 1) (Gorman et al., 2002; Carlson et al., 2004). The western part of the region (west of 110°W) has been overprinted by Cenozoic volcanic rocks, although Foster et al. (2012) suggested on the basis of isotopic signatures that the WP boundary runs through the Pioneer batholith. The GFTZ has been alternately interpreted as an intracratonic tectonic zone (O’Neill and Lopez, 1985), a Proterozoic suture between the Medicine Hat Block and the WP (Mueller et al., 2002, 2005), or an Archean shear zone or suture (Boerner et al., 1998; Buhlmann et al., 2000; Gorman et al., 2002).

The resistivity structure imaged between the WP and the MHB (Figs. 3 and 4) may be the result of protracted orogenesis occurring over two distinct phases marked by the north and south HCBs. Arguments for a two-phase diachronous orogen have been advanced by others (Harms et al., 2004; Mueller et al., 2005, 2016; Condit et al., 2015; Gifford et al., 2018) and are similar to that proposed for the THO (Bedrosian and Finn, 2021). The earlier phase of the Great Falls orogeny, which we refer to as the Medicine Hat phase, involved subduction of oceanic crust beneath the Little Belt arc (Fig. 5A). This phase is identified mainly from the Little Belt Mountains, although correlative rocks also are present in the Pioneer Mountains (Foster et al., 2006). Ca. 1.86–1.85 Ga gneisses of the Little Belt Mountains have major and trace element compositions indicative of a subduction-related magmatic arc, and Nd isotopic compositions suggest they were derived primarily from a depleted mantle source with little to no involvement of older crust (Mueller et al., 2002). An initially southeast-dipping subduction zone is postulated to account for the lack of crustal material incorporated into the Little Belt arc during its construction. North-dipping mantle reflectors beneath the MHB suggest a polarity reversal to northwest-directed subduction along the southern margin of the MHB (Gorman et al., 2002; Mueller et al., 2002; Gifford et al., 2018). As subduction continued, the postulated Medicine Hat ocean was consumed, likely by ca. 1.79 Ga (Mueller et al., 2016). Lower crustal xenoliths from the Sweetgrass Hills in the southern MHB yield Paleoproterozoic ages, leading Gorman et al. (2002) to speculate that Medicine Hat crust may have thickened by underplating at this time, a process that appears to have accompanied subduction and closure of the Medicine Hat ocean.

Collision of the MHB and WP followed closure of the Medicine Hat ocean in which the composite Medicine Hat-Little Belt margin overrode the northern Wyoming province, resulting in deformation, metamorphism, and partial melting within the Montana metasedimentary terrane (Mueller et al., 2002, 2005, 2016; Condit et al., 2015). These collisional events are known as the Big Sky phase of the Great Falls orogeny (Fig. 5B) (Harms et al., 2004; Mueller et al., 2016). The age of peak metamorphism migrated southward, from ca. 1800 to 1770 Ma in the Highland Mountains and Ruby Range to ca. 1750–1720 Ma in the Northern Madison Range, a progression consistent with migration of the orogen toward the foreland during protracted collision (Condit et al., 2015). The southward younging of peak metamorphism during the Big Sky phase may reflect the foreland growth of the orogen as thrust sheets of MHB and LB arc rocks advanced over the northern edge of the WP (Condit et al., 2015), analogous to what has been argued for the Cheyenne belt (Houston et al., 1989; Morozova et al., 2005) and along the distant Penokean orogen (Schulz and Cannon, 2007). A similar set of thrust sheets may have advanced to the northwest during the earlier Medicine Hat phase; remnants of these now eroded thrust sheets may be reflected in the linear magnetic-field anomalies near 48°N, 108°W (Fig. 2A).

In north-central Utah, the Farmington Canyon Complex coincides with an abrupt north-trending transition in crustal resistivity and conductance (Fig. 3). In sharp contrast to high resistivity within the WP (>1000 Ωm), elevated midto lower-crustal conductivity to the west is characteristic of the Great Basin writ large, where small amounts of fluids and melt within the deep crust are interpreted to drive geothermal and mineral systems (Wannamaker et al., 2008). This modern-day tectonic boundary, which coincides with the Wasatch fault zone (Fig. 1), is further reflected by gradients in heat flow and seismic velocity (Fig. 2) as well as in the distribution of geothermal resources, which cluster along the Wasatch front and to the west of it (Figs. 2 and 3).

This transition marks the modern boundary between actively extending Basin and Range Province and tectonically stable North America, however, it may further mark the western edge of the Archean WP. As viewed in cross-section (Fig. 4), the Farmington zone is distinct from the CB, THO, and GFTZ in that it lacks a discrete HCB. This, however, is not unexpected since the modern overprint of conductive crustal fluids (aqueous or magmatic) is expected to mask any lower-crustal conductors of mineralogic origin. Age data from the region suggest the lateral contrast in resistivity marks the Archean province margin, with ages >2.4 Ga just east of the boundary (Mueller et al., 2011) and a lack of ages >2.0 Ga farther west (Nelson et al., 2011). Paleoproterozoic passive-margin sequences and metamorphism at 1.74 Ga (Nelson et al., 2002; Mueller et al., 2011), coupled with nearby Paleoproterozoic arc rocks, together support this boundary marking a relict Archean continental margin similar to those imaged along the southern, eastern, and northern province margins.

North of 42°30′N, where the Wasatch fault system terminates, the western margin of the WP is speculative. The lack of exposed basement hinders our ability to precisely define this boundary. Most geophysical data in this region (seismic, gravity, magnetic, and heat flow) are dominated by the Neogene to modern overprints of Basin and Range extension and Yellowstone-Snake River Plain magmatism. Stanciu et al. (2016) note a step in crustal thickness near 44.5°N, 113.5°W between the Farmington zone and Grouse Creek block (Fig. 1)—this we take to represent the maximum western extent of the WP. Similarly, Archean ages within the Teton Range constrain the margin to lie west of it (111°W). Farther north, Archean and Paleoproterozoic age dates constrain the province boundary to project just east of the Tendoy Mountains in southern Montana (Fig. 1) (Mueller et al., 2016).

We use the resistivity model to estimate the location of the western margin in this region. At upper-crustal depths, Laramide structures obscure the margin whereas the model at lower-crustal depths reflects active tectonics; it is only at mid-crustal depths (Fig. 3B) that the model preserves a conductive-to-resistive transition that could be interpreted as a remnant of the province margin as (from 43°N to 44.5°N). Unlike farther south along the Farmington zone, where edge of Basin and Range extension appears to coincide with the Archean province margin, here the imprint of Cenozoic tectonics (e.g., Miocene initial uplift of the Teton Range, Eocene volcanism in the Absaroka Range) penetrates 200–300 km east into the Wyoming province (Figs. 1–4). The effect of such events on the stability of the WP is discussed in a subsequent section.

Neoarchean ages for basement granitoids within the Black Hills uplift, together with isotopic studies linking these same rocks to the core of the WP, are traditionally interpreted as evidence that the Black Hills remained connected to the core of the WP throughout the Proterozoic (Gosselin et al., 1988; McCombs et al., 2004). The Black Hills also record (1) evidence of one or more rifting events between ca. 2.48 and ca. 1.96 Ga inferred to represent the breakup of Kenorland, as well as (2) regional folding and metamorphism by ca. 1.85 Ga, (3) deformation events related to accretion along the Cheyenne belt and closure of the Trans-Hudson orogen, and (4) post-orogenic magmatism at ca. 1.72 Ga (Gosselin et al., 1988; Redden et al., 1990; Dahl et al., 1999; McCombs et al., 2004; Hrncir et al., 2017). The superposition of these latter three events together suggests a fundamental tectonic boundary lies within or just east of the Black Hills.

More recently, Worthington et al. (2016) argued that changes in the geophysical character of the crust east of the Bighorn uplift suggest a diminished WP, with a proposed margin just east of the Bighorn Mountains (blue line, Fig. 1). This alternate interpretation implies that the Black Hills are not rooted in Archean lithosphere but are rather part of an Archean crustal remnant afloat within a broad Proterozoic orogen (Chamberlain and Mueller, 2019). This proposed margin and the implied WP extent have been incorporated in some plate reconstructions for the amalgamation of Laurentia (Kilian et al., 2016), although aspects of this model have been called into question (Hrncir et al., 2017).

The evidence for interpreting the WP margin immediately east of the Bighorn uplift include a west-to-east change in seismic reflectivity, the absence of a 7.x layer beneath the Bighorn uplift, and a high-angle magnetic boundary at the edge of the Bighorn gneissic domain (Worthington et al., 2016). With the addition of our resistivity model, we argue that the totality of evidence does not support a Proterozoic suture adjacent to the Bighorn uplift but instead suggests the Archean WP extends east to the Black Hills as originally proposed. The presence of thick resistive lithosphere (Figs. 3, 4, and 6) to the east of the Bighorn uplift is most consistent with Archean craton and in sharp contrast to surrounding Proterozoic lithosphere (compare for example the Colorado Province and WP in Fig. 4). In addition, the absence of a high-conductivity belt along the proposed province margin is at odds with imaged HCBs accompanying other Paleoproterozoic sutures (e.g., CB, THO, and GFTZ).

In considering seismic evidence for the eastern province margin, there is little imprint on the Moho across the Bighorn uplift, which varies smoothly throughout the region as part of a broad crustal upwarp (Yeck et al., 2014). Seismic anisotropy, though dominated by asthenospheric fabric, also shows no change in either direction or intensity of shear-wave splitting across the Bighorn Mountains (Hongsresawat et al., 2015). In contrast, receiver functions image a west-to-east increase in crustal thickness of ~10 km near the Black Hills (Thurner et al., 2015). Furthermore, the presence of a 7.x layer east of the Bighorn Mountains and beneath the Black Hill uplift (Schulte-Pelkum et al., 2017) is consistent with other parts of the BBMZ (Snelson et al., 1998) that are known to be Archean and also exhibit 7.x layers. The reflectivity boundary imaged by Worthington et al. (2016) just east of the Bighorn uplift may reflect the juxtaposition of different Archean domains as opposed to Archean against Proterozoic crust. Finally, we note that at least one of the arcuate magnetic highs characteristic of the BBMZ (Fig. 2A), more than 200 km in length and attributed to Archean granitic rocks (Sims et al., 2001), is continuous across the proposed Proterozoic boundary of Worthington et al. (2016).

We thus consider it most probable that west-to-east variations in seismic reflectivity and 7.x-layer thickness reflect structural domains internal to the Archean WP and that the eastern edge of the Archean WP more likely lies beneath or just east of the Black Hills uplift. Our preferred eastern province margin (Fig. 1) is demarcated by the western HCB of the Trans-Hudson orogen (Bedrosian and Finn, 2021) and supported by geologic evidence, Archean age determinations in the Black Hills and Hartville uplifts, and a change in Moho depth across the Black Hills.

The resistivity model, in combination with metamorphic age determinations, can be used to place constraints on the time of formation of bounding HCBs and the Paleoproterozoic construction of Laurentia. Most definitive is the relative age of the HCB along the Big Sky phase of the Great Falls orogen, which is truncated by the Trans-Hudson HCB, suggesting that final collision of the MHB and WP by ca. 1.79 Ga preceded subduction along the western edge of the Trans-Hudson orogen (Mueller et al., 2016). In the Black Hills, north-south shortening dated at ca. 1.79–1.78 Ga (Dahl et al., 1999) is similar to timing and shortening directions identified along the Cheyenne belt, Hartville uplift, and Laramie Mountains that are interpreted to reflect accretion from the south (Allard and Portis, 2013). Slightly younger, ca. 1.77–1.74 Ga shortening recorded in the Black Hills is oriented east-west and is interpreted to record collision of the WP with the Superior Province along the southern Trans-Hudson orogen. Final suturing appears to have been complete prior to the intrusion of the ca. 1.72 Ga Harney Peak Granite in the Black Hills (Allard and Portis, 2013). Closure along the southern Trans-Hudson orogen after collision along the Cheyenne belt may help account for the lack of a continuous conductor between the Cheyenne belt and its extension farther east (Bedrosian and Finn, 2021). Taken together, a remarkable series of orogenic events affected the southern Laurentian margin, including all sides of the WP, during this relatively short interval in geologic time.

The Archean WP, as defined in Figure 1, has an estimated area of more than 350,000 km². With the exception of the northwest corner, the entirety of this region is resistive (>500 Ωm) to depths in excess of 200 km (Figs. 3, 4, 6, and 7A). This same domain is characterized by low heat flow (Blackwell et al., 2011) and moderate to high seismic velocities (e.g., Schmandt and Humphreys, 2010a; Shen and Ritzwoller, 2016). As will be discussed below, the resistivity model beneath the WP is most consistent with dry, refractory lithosphere to depths in excess of 200 km.

In sharp contrast to the HCBs bounding three sides of the WP, the boundaries between internal domains of the WP are marked by a noted lack of high conductivity. We focus our attention here on the ε_Nd = -4 isotopic corridor separating the SAT and BBMZ (Fig. 1), as this boundary is not affected by the overprint of Basin and Range extension and Yellowstone-Snake River Plain magmatism (as is the MMT-BBMZ boundary). As discussed previously, we interpret the isotopic and geochronologic changes across this boundary, together with the presence of supra-crustal belts of juvenile material, as indications of Neoarchean juxtaposition of these two distinct sub-terranes of the WP. In addition to the lack of a 7.x lower crustal layer in the SAT (Snelson et al., 1998; Morozova et al., 2002), its mantle lithosphere is characterized by lower V_P and V_S perturbations compared to the BBMZ (Figs. 7B–7E).

Along the ~300 km length of the ε_Nd = –4 corridor, there is virtually no elevated conductivity in the lithospheric column (Figs. 3, 4, and 6). The Laramide uplifts visible within the upper crust of the resistivity model (Fig. 3A) show some alignment with the ε_Nd = –4 corridor, suggesting a degree of structural inheritance between this Archean boundary and deformation during the Laramide orogeny. At mid- and lower-crustal depths, a resistive “corridor” follows this same boundary and may mark the ca. 2.63 Ga magmatic arc rocks that stitch the SAT to the BBMZ. This is further reflected in crustal conductance (Fig. 3D), where a dark band of low conductance follows Neoarchean granite and granodiorite within the Wind River Range and Granite Mountains (Stuckless et al., 1985; Frost et al., 1998; Bagdonas et al., 2016). The implications of the dramatically different resistivity expression of this Neoarchean suture compared to that of the Paleoproterozoic CB, THO, and GFTZ are discussed below.

The ε_Nd = –4 corridor marks the Neoproterozoic collision of juvenile terranes against the southern margin of the BBMZ. Ca. 2.64 Ga syncollisional granite intrusions north of the suture (Frost et al., 2006b; Frost and Da Prat, 2021) suggest subduction was directed beneath the BBMZ. Following accretion, voluminous subduction-related continental-arc magmatism continued and the suture at the surface was largely overprinted by intrusion of voluminous granitic batholiths (Bagdonas et al., 2016). Despite this evidence for modern-style plate tectonic processes within the Wyoming Province from at least ca. 2.8 Ga (Mogk et al., 2022), the geophysical signature of the sub-province boundary between the BBMZ and SAT is dramatically different from that of Paleoproterozoic orogens that mark the WP boundaries. Crustal conductance is 100–1000 times lower in the Archean suture than along the surrounding Paleoproterozoic sutures (Fig. 3D). Does this dichotomy reflect changes in subduction processes between the Neoarchean and the Paleoproterozoic?

The source of high conductivity beneath ancient sutures is considered to be organic carbon and metallic sulfides deposited, most commonly as black shales, along long-lived passive margins within an anoxic or euxinic water column (Boerner et al., 1996). Sediment subduction transports these rocks to lower-crustal depths where the carbon is metamorphosed to graphite. The resulting graphitic and sulfidic metasedimentary rocks can remain stable over geologic time scales to depths of 100–150 km. In considering the absence of HCBs elsewhere, Murphy et al. (2022) advanced three possible explanations: (1) magmatic or structural dismemberment of an original HCB by later tectonics, (2) incompatible sediment composition or redox state needed to preserve carbon and/or sulfides, or (3) insufficient volume of metasedimentary material.

Is it possible that granite magmatism along the ε_Nd = –4 corridor intruded and overprinted an earlier HCB between the BBMZ and SAT? Syn- and post orogenic magmatism are abundant along the ε_Nd = –4 corridor, however we consider it unlikely that this magmatism was extensive enough to erase all remnant conductivity along this 300-km-long corridor. As a point of comparison, extensive 1.4 Ga magmatism within the southern Laramie Mountains (Sims et al., 2001) occurred along the Cheyenne belt, however, our models image a continuous HCB across the region of 1.6 Ga magmatism (Figs. 3B–3D).

Alternately, passive-margin assemblages in the Neoarchean may not have been enriched in carbon or sulfides, and thus unable to produce a conductivity anomaly. For example, a purely quartzofeldspathic sedimentary pile is not expected to yield a conductivity anomaly. There is, however, evidence that these conductive mineral phases were present on Earth at the time subduction was occurring along the ε_Nd = –4 corridor. Maxima in the global abundance of black shales occur from ca. 2.7 to 2.5 Ga and again from ca. 2.0 to 1.7 Ga (Condie, 2001); these abundance peaks cover both the periods of subduction we compare. The lack of a euxinic water column in the Neoarchean (Canfield, 1998; Poulton et al., 2004, 2010) would have limited sulfide precipitation from seawater, thus removing one mechanism for generating conductive phases. However, prior to the rise of atmospheric oxygen, subaerial sedimentary s uccessions (≥2.45 Ga) include abundant reduced-facies minerals, for example pyrite pebbles within placer deposits (Rasmussen and Buick, 1999; England et al., 2002). Thus, even prior to the ca. 1.8 Ga transition from a ferruginous to euxinic ocean, which led to a pulse of seafloor sulfide deposition, we expect a baseline of organic carbon and sulfides preserved within black shales beneath the anoxic Archean ocean (Condie, 2001).

A final possibility is insufficient sediment input. Passive-margin supracrustal rocks along the ε_Nd = –4 corridor are sparse, restricted to the ca. 2.85 Ga Barlow Springs Formation in the Granite Mountains (Langstaff, 1995) and ca. 2.64 Ga pelitic schists of the Bluegrass Creek area in the Laramie Mountains (Patel et al., 1999). In contrast, an abundance of sediments deposited in the Paleoproterozoic, shed from the recently emerged continental landmass, resulted in the emplacement of HCBs along the CB, THO, and GFTZ. Much thicker sequences of Paleoproterozoic supracrustal rocks are present along the WP margins, including the Snowy Pass Supergroup along the Cheyenne belt (Karlstrom et al., 1983), and quartzite, metapelite, and meta-graywacke platform to deep-marine sequences in the Black Hills (Redden et al., 1990). We consider a relative lack of sedimentary input in the Archean to be the most plausible explanation for the lack of an HCB along the ε_Nd = –4 corridor.

Around 2.5 Ga continental freeboard is thought to have increased, as recorded in O₂ isotope data (Bindeman et al., 2018), in response to mantle cooling, decreased radiogenic heat production and the strengthening of continental lithosphere (Flament et al., 2008). This likely led to a dramatic rise in continental weathering and in turn to the generation and burial of organic carbon from 2.22 to 2.06 Ga as evidenced by a positive δ¹³C excursion (Karhu and Holland, 1996; Kump et al., 2011). These changes appear to have been temporally linked to the rise in atmospheric O₂ on Earth (Bekker, 2003; Bekker et al., 2004). Molybdenum enrichments in black shales provide an independent record of continental weathering (Scott et al., 2008), also exhibiting a pronounced increase from 2.2 to 2.0 Ga coincident with organic carbon burial and the rise in atmospheric O₂ content. These global changes account for the appearance of HCBs in the Paleoproterozoic.

Is the lack of HCBs prior to the Paleoproterozoic unique to the WP? To examine this further we compile MT studies over known sutures ranging in age from the Archean through the Cenozoic (Table S1). This compilation suggests the geophysical distinction between Archean and Paleoproterozoic sutures may be global. Considering only those studies employing 3-D inversion of array data (where the geometry and along-strike extent of conductors can reliably be estimated), the compilation reveals that: (1) HCBs of 100 km to 1000 + km in length are imaged along sutures as young as Cenozoic in age, (2) there exists a relative abundance of Paleoproterozoic sutures accompanied by HCBs, and (3) there are no HCBs imaged beneath Archean sutures.

This final observation reflects both a lack of unequivocally Archean sutures and the limited number of 3-D MT studies in Archean terranes. Nevertheless, when considering less reliable studies (2-D inversion and/or 2-D data sets), only three instances of elevated conductivity beneath recognized Archean sutures are identified—each with caveats (Table S1). We therefore speculate that increased continental weathering, driven by the emergence of continents from the Neoarchean ocean, led to a dramatic increase in sedimentary input to subduction zones and the formation of HCBs beginning in the Proterozoic and continuing through to today. Superimposed on this step change is the relative abundance of Paleoproterozoic HCBs, reflecting a pulse of seawater pyrite precipitation ca. 1.8 Ga, driven by a rise in euxinic conditions along continental margins in response to an increased flux of organic carbon and accompanied by the disappearance of dissolved iron in the oceans (Canfield, 1998; Poulton et al., 2004, 2010).

Studies invoking modification or destruction of part or all of the WP (Kusky et al., 2014; Dave and Li, 2016; Snyder et al., 2017; Bedle et al., 2021) cite as evidence for: (1) the absence of thick lithosphere, (2) petrologic and geochemical data from xenoliths suggesting metasomatism of the deep lithosphere or wholesale replacement of the Archean subcontinental lithospheric mantle (SCLM), and (3) reduced seismic velocities beneath the SAT and in northeast Wyoming. We consider each line of evidence below, but we first examine the resistivity model and the constraints it provides in terms of geometry and resistivity.

Characteristic 1-D profiles through the BBMZ and SAT, the locations of which are indicated on Figure 6A, are shown on Figure 7A, as well as at a second location within the BBMZ that overlaps with the eastern Snake River Plain. For comparison, we also include 1-D resistivity profiles from the Basin and Range Province as well as the southern Superior Province (46°N, 99°W). Since the SP falls outside the model domain of this study, the latter profile is extracted from the model of Bedrosian and Finn (2021). At each location, resistivity curves are averaged within a 0.5 degree box; seismic velocity curves show V_P and V_S perturbations extracted at these same locations from the model of Schmandt and Humphreys (2010b) which has a 0.25 degree resolution (Figs. 7B and 7C). We focus our attention on the base of the lithosphere, where the effect of metasomatism, such as from flat-slab subduction during the Laramide (Coney and Reynolds, 1977; Saleeby, 2003), is expected to be most pronounced.

Examination of the 1-D resistivity curves shows the BBMZ and SAT to be in lockstep within the mantle lithosphere to depths of 150 km and with resistivities of 500–1000 Ωm (Fig. 7A). Both V_P and V_S perturbations are similarly aligned to 125 km depth, below which the BBMZ and SAT differ by ~2% (Figs. 7B and 7C). Based upon synthetic resolution tests (Schmandt and Humphreys, 2010a, their fig. 3), we estimate the resolution of these models to be on the order of 1%–2%, comparable to the imaged velocity difference between the BBMZ and SAT. The resistivity profile through the southern SP is more resistive than the WP but they overlap near the base of the lithosphere. Some of this difference stems from the screening effect of overlying sedimentary basins (Green River and Powder River basins), which act to depress the inverted resistivity model beneath the WP. In velocity, the SP is indistinguishable from the WP (Figs. 7B and 7C) whereas profiles through the tectonically active Basin and Range and SRP provinces are considerably slower (>4%) and more conductive (10×) than the BBMZ, SAT, and SP.

The resistivity of continental lithosphere is to first order controlled by temperature. The sharp contrast between the active and inactive provinces (Figs. 3C, 3D, 4, and 6) thus reflects dramatic differences in geotherms and surface heat flow values (30–50 mw/m² versus 80–110 mW/ m2; Fig. 2C). Both the Yellowstone-Snake River Plain (YSRP) and BBMZ soundings are located within Archean lithosphere of the WP (Fig. 6A), however, the geophysical signature of the former bears little relation to the latter due to the combined effects of temperature, melt, and fluids. While all indications (V_P, V_S, heat flow, resistivity, volcanism) point to modification of this northwest corner of the WP, it is uncertain whether the long-term stability of the SCLM is at risk. The majority of Archean provinces on Earth maintain thick and cool lithospheric mantle, protecting the overlying crust from subsequent reworking and recycling (Kusky et al., 2014; Pearson et al., 2021; Bedle et al., 2021). Geodynamic modeling by Wang et al. (2015) suggests that mantle plumes cannot on their own erode depleted cratonic mantle without precursory metasomatism. Whether the northwest corner of the WP remains stable after passage of the Yellowstone hot spot may thus depend upon the extent of metasomatism.

What can be said of the degree and extent of metasomatism beneath the WP? Turning to those inactive parts of the province (characterized by low heat flow and a lack of post-Archean internal deformation), variations in lithospheric resistivity cannot be explained in terms of temperature or partial melt. Compositional variations are also unlikely to produce significant resistivity variations without appealing to specific minerals such as phlogopite or garnet that are unlikely to be present in the quantities needed to dramatically lower bulk resistivity. Jones (1999) and Selway (2014) both demonstrate the minor effect of composition on lithospheric mantle conductivity. The most likely explanation for observed conductivity variations in the SCLM is dissolved water in nominally anhydrous minerals. Olivine, the major component of mantle peridotite, can accommodate water at concentrations up to a few hundred parts per million (ppm) (Bell and Rossman, 1992) and laboratory data show significant enhancement in conductivity for olivine, orthopyroxene, and garnet at low temperatures and modest water content (Dai and Karato, 2009; Yoshino et al., 2009; Poe et al., 2010). Thus, these hydrated minerals can considerably reduce bulk resistivity relative to their dry mineral counterpart.

The resistivity model provides evidence for intact Archean lithosphere beneath both the BBMZ and the SAT. Cross-sections and depth slices (Figs. 4 and 6) image thick resistive lithosphere over a conductive asthenosphere, with the electrical lithosphere-asthenosphere boundary (LAB) at depths of 200–250 km and thickening to the southwest beneath the SAT. This picture contrasts with surrounding Proterozoic lithosphere, such as the Colorado Province (Fig. 4), which is considerably thinner (120–150 km) from both electrical (Fig. 4) and seismic (West, 2004) LAB estimates. Beneath the BBMZ and SAT, this thick lithospheric mantle has likely persisted since the Neoarchean. The stability of SCLM is dependent upon its state of dehydration, making it viscous relative to the underlying asthenosphere and thus resistant to delamination. Extensive melting during craton stabilization is considered to have depleted the SCLM of most of its water. Given that increased water content in olivine results in a corresponding decrease in viscosity, it is reasonable to expect the olivine water content of SCLM to be less than that of the asthenosphere, which is estimated to be ~100 wt. ppm (Dai and Karato, 2009). Studies of SCLM xenoliths by Peslier et al. (2010) also suggest the water content of SCLM is typically <100 wt. ppm, although it may vary throughout the lithospheric column.

Resistivity profiles for both the BBMZ and the SAT (Fig. 7A) show a general decrease in resistivity with depth from the Moho to the LAB. We compare these inverted resistivity-depth models to the unified hydrous conductivity model of Gardés et al. (2014), a model of olivine conductivity as a function of temperature and water content that is calibrated against laboratory measurements. As input to this model we assume a conductive geotherm for the region (red line, Fig. 7A) that incorporates variable heat production and an estimated surface heat flow of 40 mW/m² (Hasterok and Chapman, 2011). The LAB implied by this geotherm is at 190 km depth, broadly consistent with both our electrical and seismic LAB estimates described below.

Figure 7A shows the Gardés et al. (2014) model for a dry lithospheric mantle and for mantle with water concentrations of 100 wt. ppm olivine, which as described above we consider a “breaking point” for lithospheric stability. The resistivity profiles beneath the BBMZ and SAT, specifically near the base of the lithosphere where widespread metasomatism is expected to be most likely, are consistent with a dry olivine mineralogy. Considering the resistivity contrast between the WP and neighboring Proterozoic provinces (Figs. 3, 4, and 6), the comparable thickness of WP and SP lithosphere (Fig. 6D), and the thinner lithosphere beneath the Proterozoic Colorado Province (Fig. 4), we conclude that the core of the WP maintains a thick stable Archean lithosphere and that the effects of metasomatism, whether during the Laramide or otherwise, have not considerably affected the stability of the WP. Previous MT modeling by Meqbel et al. (2014) also conclude that the average resistivity structure beneath the WP is consistent with the olivine dry solidus.

Whereas WP lithosphere appears to be largely intact, the resistivity model provides additional insight into tectonic processes affecting its margins. Figure 7D and 7E estimates the extent of the WP at various depths as imaged in the resistivity model. Contours range from 300 Ωm at 60 km depth to 100 Ωm at 200 km, with the latter approximately corresponding to the dry olivine resistivity model near the base of the lithosphere (Fig. 7A). As such, these contours can loosely be viewed as bounding that portion of the WP consistent with intact Archean lithosphere. Areas outside of these contours but within the interpreted WP surface extent (black line) are then areas where the WP may have been modified through some combination of melt, fluids, or metasomatism.

Area 1 (Figs. 7D and 7E) is most clear cut: starting at mid-crustal depths and continuing through the entire lithospheric mantle is a prominent “bite” out of the resistive WP (Fig. 3, 4, 6, and 7) that encompasses all of the MMT and a portion of the BBMZ. Low seismic velocities and elevated heat flow also characterize this region. The affected region includes Yellowstone and the eastern Snake River Plain but notably continues at least 100 km east of Yellowstone and nearly 200 km south of it.

Area 2, in the southwest WP, is a promontory of Archean crust extending west to the Wasatch Fault but with more conductive lithosphere beneath it than is typical of the SAT. We interpret this ~100 km ingress of conductive lithosphere to reflect elevated geotherms due to its adjacency to the Basin and Range. This may in turn be amplified by small-scale convection at what is presumed to be a step in lithospheric thickness between the eastern Great Basin and the SAT (Wannamaker et al., 2008; Hopper and Fischer, 2015). The reduced lithospheric velocities beneath the Basin and Range province also penetrate into this region (Figs. 7D and 7E).

Area 3 is an enigmatic conductive “tongue” (Figs. 3C and 6A), less than 50 km in width, projecting into the center of the WP from the SRP. We note that the resistivity model images this feature in the lower crust and upper mantle, but only to ~100 km depth. The geometry of this feature is further reflected in higher heat flow (Fig. 2D), in reduced V_P at 60–100 km depth (Fig. 2D), and in the distribution of Cenozoic volcanic rocks, specifically the Eocene Absaroka volcanic province (Fig. 1). The spatial correlation with Eocene magmatism and elevated heat flow suggests a pre-Tertiary structural corridor that has served as a locus for deformation and guided heat, melt, and fluids through the crust.

We now consider the xenolith and seismic data that have been interpreted as evidence for modification or destruction of part or all of the WP. A primary argument for loss of the Archean SCLM (e.g., Humphreys et al., 2015; Snyder et al., 2017) comes from the Williams and Homestead peridotites in the Missouri Breaks and Grassrange areas, respectively (Fig. 1), that sample depths of 130–170 km and yield Re-Os model ages of <200 Ma (Carlson et al., 2004). Given the revised boundaries of the WP (Fig. 1), both localities now fall outside of the WP and are instead located within the MHB or Great Falls orogen. The only evidence for metasomatism that falls squarely within the WP is from small-volume lamproites within the Leucite Hills, which record two generations of metasomatism in veins, one at >1 Ga and a younger event <100 Ma (Mirnejad and Bell, 2006). The nature of these infiltration events are important: metasomatism that is localized in veins is unlikely to affect the average composition and stability of the SCLM (Carlson et al., 2005). The presence of thick resistive lithosphere to 250 km depth throughout the SAT (Figs. 4 and 6) as well as the limited extent of the Leucite Hills lamproites suggest that neither the ancient nor recent metasomatism beneath the Leucite Hills have severely affected the stability of the WP.

The second line of evidence for craton modification, the seismic LAB, is notoriously difficult to constrain beneath cratonic regions, and to our knowledge there exists no compelling evidence to support claims that the WP lithosphere is thinner than that of typical cratons (200–250 km). While global LAB models (e.g., Priestley et al., 2018) indicate an LAB ranging from 80 to 170 km beneath the WP, these models are of low spatial resolution (2°), with the WP falling along a broad lateral gradient between thin lithosphere of the Basin and Range and thick lithosphere beneath the Great Plains. The highest resolution LAB estimates come from S_P-phase receiver functions, which capture large and rapid velocity gradients due to an anomalously hot asthenosphere (Fischer et al., 2010; Abt et al., 2010; Hopper et al., 2014); the absence of such phases beneath the WP (and cratons in general) indicates that the LAB beneath the WP is marked by a diffuse velocity gradient extending over 50–70 km in depth (Fischer et al., 2010; Hopper and Fischer, 2015). Any seismically determined LAB depth estimates for the WP bear this degree of uncertainty. Surface-wave tomography is one of the only means of identifying the seismic LAB beneath cratons, defined alternately as the maximum negative velocity gradient, the depth to an absolute velocity, a velocity anomaly contour, or as a depth below which lateral velocity variations cease (Eaton et al., 2009; Fischer et al., 2010). Hopper and Fischer (2015) estimate the LAB beneath the WP, applying the negative gradient approach to the model of Porritt et al. (2014). They estimate an LAB of ~175 km depth beneath the BBMZ and between 170 and 200 km beneath the SAT; both are indistinguishable from their LAB estimates for the SP. Given the uncertainty in defining the seismic LAB beneath the WP, and acknowledging differences between seismic and electrical LAB determinations (Eaton et al., 2009), we consider seismic LAB estimates to be in line with the 200–250 km electrical LAB beneath the WP (Fig. 4) and indistinguishable from the LAB depth beneath the SP (Fig. 6D).

A final argument for modification of the WP is based upon reduced velocities in the SAT lithosphere below ~125 km depth (Bedle et al., 2021) (Figs. 7C and 7E). The Deep Probe experiment (Fig. S1) identified reduced P-wave velocities in the upper mantle beneath the SAT relative to the BBMZ (Henstock et al., 1998; Snelson et al., 1998). The model of Schmandt and Humphreys (2010b) also find V_P and V_S to be lower beneath the SAT than the BBMZ (Figs. 7B–7E). In both studies, velocity differences between the two sub-provinces are on the order of 1%–2%, a level which can be explained by differences in composition or geotherm without necessarily invoking metasomatism. Given their distinct Archean histories prior to 2.65 Ga, it is not unexpected that the BBMZ and SAT have distinct geophysical signatures. Snelson et al. (1998) concluded that upper-mantle velocities beneath the WP, though slightly lower than that of typical cratons, the mantle lithosphere beneath the SAT is warmed, but otherwise intact.

We also note the reduction in upper-mantle resistivity under the westernmost SAT (Figs. 6A and 6B), which we attribute to elevated temperatures adjacent to the eastern Basin and Range Province. Although the SAT lacks the Paleoarchean history and the lower-crustal 7.x-layer characteristic of the BBMZ, we find no evidence of widespread lithospheric modification or destruction that cannot be alternatively attributed to an elevated geotherm. In summary, none of the evidence for widespread craton modification is definitive; to the contrary, the WP lithosphere appears to remain thick, buoyant, and stable to this day.

Our results are not, however, at odds with models of Laramide-flat-slab subduction nor with some degree of attendant metasomatism resulting from a period when oceanic Farallon lithosphere extended beneath the province. The totality of evidence across the western United States, including the shutoff of arc magmatism, widespread foreland deformation, and time-progressive volcanism during slab rollback, is elegantly explained by flat-slab subduction models (Humphreys, 1995), while modern analogues (e.g., Andes; Yakutat slab, Alaska, USA) offer supporting comparisons (Ramos and Folguera, 2009; Abers, 2013). The Cenozoic magmatic record is also consistent with metasomatism beneath the WP—the question is to what degree. Outside of the SRP, volcanism is focused along the periphery of the WP (Fig. 1). One exception to this is the Eocene Absaroka volcanic province, derived from an ancient depleted mantle source but with isotopic and trace-element ratios recording mantle metasomatism during shallow subduction of the Farallon plate (Feeley, 2003). This region, we note, is spatially coincident with the conductive “tongue” imaged at <100 km depth within the SCLM (Fig. 6A). Other volcanism within the WP is volumetrically minor (e.g., Leucite Hills, Rattlesnake Hills), and as described above, shows indications of metasomatism, though not pervasive. We thus conclude that the effects of metasomatism are highly variable, and that metasomatism due to flat-slab subduction may modify, but not necessarily destroy, cratonic lithosphere (Bedle et al., 2021).

The Archean Wyoming Province and its bounding Paleoproterozoic orogens span an interval of Earth history that experienced some of the planet’s most fundamental changes, including the onset of plate tectonics, subaerial emergence of the continents and increased continental weathering, assembly of the first super-continents, shifts in ocean chemistry, the rise of atmospheric oxygen, and increased biogenic carbon. Many of these events are recorded in the resistivity structure of the lithosphere, as illustrated by this study.

• A 3-D electrical resistivity model, combined with a new geochronology compilation, images high conductivity belts surrounding Archean crust that identify the extent of the Wyoming Province, particularly along its unexposed northern, eastern, and southern boundaries. The model indicates that the northwest boundary is 150 km south of prior estimates and supports an eastern boundary that encompasses the Black Hills.

• Internal boundaries between Archean subprovinces are not marked by deep-crustal conductive belts, suggesting an absence of conductive metasedimentary rocks at depth along these sutures. The onset of conductive sutures in the Paleoproterozoic appears to be a global phenomenon and may reflect an increase in continental weathering between 2.6 and 1.8 Ga as continents emerged from the Neoarchean ocean.

• Although the Wyoming Province has been identified as a craton with modified lithosphere, evidence for modification is limited to its western margin, where the Snake River Plain-Yellowstone hotspot, Basin and Range extension, and extensive Eocene magmatism have introduced heat and fluids. Metasomatism related to shallow subduction of the Farallon plate may cause enduring compositional changes but outside of the northwest corner, the WP maintains a thick and presumed stable lithospheric keel.

Resistivity models offer a powerful means with which to image Archean lithosphere and the orogens that bound them. The corresponding insensitivity of resistivity to upper-mantle composition, when coupled with seismic models, permits an understanding of continental architecture even in areas impacted by modern tectonics, where the effects of temperature, melt and metasomatism are otherwise difficult to distinguish. The Archean WP is encircled by Paleoproterozoic orogens, suggesting that cratonic lithosphere has not diminished since its formation in the Archean. The 3-D extent of cratonic lithosphere on Earth today may in fact remain close to its original extent, and modification may lead to craton destruction rarely and only under exceptional circumstances.

This paper grew out of a 2019 EarthScope synthesis workshop held in Bozeman, Montana, USA. This work was supported by the U.S. Geological Survey (USGS) Mineral Resources Program. The magnetotelluric data used in this study are available through the Incorporated Research Institutions for Seismology Data Management Center at http://ds.iris.edu/spud/emtf/. Both the 3-D resistivity model and the seismic model of Schmandt and Humphreys (2010b), to which we compare, are available at http://ds.iris.edu/ds/products/emc/. Isostatic gravity grids are available from https://mrdata.usgs.gov/gravity/ and http://gdr.agg.nrcan.gc.ca/gdrdap/dap/index-eng.php. The magnetic data are available from https://mrdata.usgs.gov/magnetic/). This research used resources provided by the CoreScience Analytics, Synthesis, and Libraries Advanced Research Computing group at the USGS. We thank Jamey Jones and two anonymous reviewers for their helpful comments, and Ben Drenth, Carol Finn, Heather Ford, Dave Mogk, Paul Mueller, Ben Murphy, and Jared Peacock for discussions on aspects of this work. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.