The Bridger Range of southwest Montana, USA, preserves one of the most temporally extensive sedimentary sections in North America, with strata ranging from Mesoproterozoic to Cretaceous in age. This study presents new detrital zircon geochronologic data from eight samples collected across this mountain range. Multidimensional scaling and non-negative matrix factorization statistical analyses are used to quantitatively unmix potential sediment sources from these and 54 samples compiled from previous studies on regional correlative strata. We interpret these sources based on reference data from preserved strata with detrital zircon signatures likely representative of ancient sediment sources. We link these sources to their sinks along sediment dispersal pathways interpreted using available paleogeographic constraints. Our results show that Mesoproterozoic strata in southwest Montana contain detritus derived from the nearby craton exposed along the southern margin of the fault-bounded Helena Embayment. Middle Cambrian strata were dominated by the recycling of local sources eroded during the development of the Great Unconformity. In Devonian–Pennsylvanian time, provenance in southwest Montana shifted to more distal sources along the northeastern to southeastern margins of Laurentia, but more western basins received detritus from outboard sources along a tectonically complicated margin. By the Late Jurassic, provenance in the developing retroarc foreland basin system was dominated by Cordilleran magmatic arcs and fold-thrust belt sources to the west. Eastward propagation of the fold-thrust belt caused recycling of Paleozoic and Jurassic detritus into the foreland basin to dominate by the Early Cretaceous.

Sedimentary rocks in the Bridger Range of southwest Montana, USA, record the Mesoproterozoic–Mesozoic history of sedimentation in this region of North America, with Mesoproterozoic to Cretaceous stratigraphy exposed as a generally eastward dipping homoclinal panel across the range (Fig. 1). The >8 km of stratigraphy collectively record critical information about source-to-sink relationships in this portion of Laurentia during much of the Mesoproterozoic, Paleozoic, and Mesozoic. During this time interval, the western margin of North America underwent rifting associated with the breakup of Proterozoic supercontinents and subsequent thermal relaxation and early Paleozoic passive margin sedimentation (e.g., Devlin and Bond, 1988; Sears et al., 1998). The western margin of North America transitioned to an active margin likely beginning in the Devonian (e.g., Dorobek et al., 1991; Crafford, 2008; Beranek et al., 2016). This margin became tectonically consolidated with eastward subduction of the Farallon plate underneath North America by Late Jurassic time, culminating in Cordilleran orogenesis and the creation of a ∼3000 km fold-thrust belt and foreland basin system (e.g., DeCelles, 2004; Yonkee and Weil, 2015). Throughout this evolution, southwestern Montana occupied a central position near the western margin of Laurentia. Consequently, the strata deposited there provide a crucial record of this tectonic history.

Siliciclastic intervals in these rocks contain detrital zircons, which can be radiometrically dated to determine their crystallization age (e.g., Gehrels et al., 2006). Sampling large populations of zircons from strata spanning a broad temporal and spatial interval can reveal changes in regional sediment routing systems through comparison of the distributions of detrital zircon ages to potential ancient sediment sources (e.g., Gehrels et al., 1996, 2011; Gehrels and Ross, 1998; Park et al., 2010; Laskowski et al., 2013). These ancient sources may consist of first-cycle grains derived from Precambrian basement provinces or more recent magmatism. However, zircons are more commonly sourced from erosion and recycling of previously deposited sedimentary rock (e.g., Schwartz et al., 2019).

The evolution of sedimentary provenance is strongly influenced by tectonic factors (Gehrels, 2014). Tectonic events, such as rifting and mountain building, trigger drainage reorganization and uplift new sediment sources. These changes are recorded in the detrital zircon record, which can be used to assess the tectonic and paleogeographic evolution of the region (e.g., Fuentes et al., 2011; Doe et al., 2012; Laskowski et al., 2013; Beranek et al., 2016). While similar analyses have been conducted on some correlative strata across the North American continent (e.g., Fuentes et al., 2011; Gehrels et al., 2011; Laskowski et al., 2013; May et al., 2013; Gehrels and Pecha, 2014), southwestern Montana is comparatively less studied. Our new data from this region provide new constraints on the evolution of the provenance of sediment feeding basins occupying this region and the coupled tectonic evolution driving these changes. Specifically, we test or constrain: (1) the hypothesis that Middle Cambrian strata deposited during the Sauk transgression sourced most of their detritus from the underlying basement (e.g., Malone et al., 2017; Matthews et al., 2018); (2) the timing and nature of the regional transition from northern/northeastern-dominated sedimentary provenance (i.e., Gehrels and Pecha, 2014) to eastern/southeastern-dominated provenance (i.e., Chapman and Laskowski, 2019) during the Paleozoic; (3) the hypothesis that exotic terranes off the west coast of Laurentia provided sediment to Paleozoic strata preserved in the Northern Rocky Mountains (e.g., Beranek et al., 2016); (4) the timing of the shift from eastern- to western-dominated provenance in the Jurassic foreland basin system (e.g., Fuentes et al., 2011); and (5) the hypothesis that Early Cretaceous intraforeland uplifts controlled regional provenance patterns (DeCelles, 1986).

In this paper, we present a detrital zircon data set sampled from a stratigraphic section in the Bridger Range spanning Mesoproterozoic through Cretaceous time. This data set, consisting of 2133 analyses from eight original samples, is coupled with 8691 analyses from 54 samples (Table 1) compiled from the Northern Rocky Mountain region (Northern Rockies), herein considered from southern Wyoming, USA to southern Alberta/British Columbia, Canada. This study couples descriptive statistical analysis (e.g., Vermeesch, 2013; Saylor and Sundell, 2016) with a newly developed inferential statistical modeling technique, non-negative matrix factorization (Saylor et al., 2019), to quantitatively unmix source age distributions from this data set. This technique honors that sediment sinks are combinations of end-member sources, which are themselves mixtures of detrital zircons from crystalline and recycled sedimentary sources. This method answers calls for statistically robust analysis from Schwartz et al. (2019), who show that North American basement ages are spatiotemporally nonunique in the detrital zircon record of the continent and, therefore, incapable of directly tying source areas to sediment sinks on the basis of their presence or absence in a detrital zircon sample alone. Through this, we hope to (1) present a workflow for provenance analysis of spatially and temporally extensive detrital zircon data sets using modern statistical techniques, (2) utilize detrital zircon provenance interpretations coupled with existing petrologic and paleocurrent constraints, where available and applicable, to better constrain sediment-dispersal routes delivering detritus to the Northern Rockies, and (3) reconstruct the paleogeographic evolution of this region through time and test existing hypotheses for the tectonic processes driving this evolution.

The north-trending Bridger Range, adjacent to Bozeman in southwest Montana, exposes Archean metamorphic and Mesoproterozoic through Cretaceous sedimentary rocks (Fig. 1). This range is situated in a region of overlap among four major tectonic provinces each with its own distinctive structural style (Fig. 1A). These are: (1) a series of extensional structures in the cratonic basement associated with the Mesoproterozoic development of the Belt Basin (Winston, 1986); (2) a fold-thrust belt active in the Jurassic–Paleocene, regionally referred to as the Sevier Belt (DeCelles, 2004); (3) a region of intraforeland basement-involved thrusts active in the Late Cretaceous(?)–Eocene, regionally referred to as the Laramide province (DeCelles, 2004; Carrapa et al., 2019); and (4) a zone of Paleocene–recent extensional deformation associated with the extensional collapse of the North American Cordillera, herein considered to include Paleogene low-angle detachments superimposed on the fold-thrust belt and later high-angle normal faults regionally considered to comprise the Basin and Range Province (Constenius, 1996).

The Bridger Range preserves Mesoproterozoic rocks deposited in the Belt Basin (Fig. 1A). The tectonic setting of this basin is the subject of much debate. Previous interpretations of the Belt Basin include a rift basin (Sears and Price, 1978), intracratonic basin (Winston, 1986), and back-arc or transtensional basin (Ross and Villeneuve, 2003). Most models, with the exception of an intracratonic basin model (e.g., Winston, 1986), require some association between rifting of the supercontinent Nuna, also referred to as Columbia (e.g., Rogers and Santosh, 2002), and the development of this basin.

The Helena Embayment formed an asymmetric east-northeast–trending branch of the Belt Basin extending along the central Montana trough and encompassing the location of the Bridger Range (Fig. 1A; Winston, 1986). This block was down-dropped in Mesoproterozoic time along bounding fault zones. The east-west–striking Perry and Garnet lines bound the southern and northern limits of the basin, respectively, and the north-south–striking Townsend line bounds the basin on its western margin (Winston, 1986). The Helena Embayment may represent a failed arm of a rift that reactivated the suture between the Wyoming and Medicine hat blocks (Harrison et al., 1974; Winston, 1986; Price and Sears, 2000; Sears et al., 2004).

In Neoproterozoic time, rifting of the supercontinent Rodinia led to the development of a passive margin along the western margin of Laurentia (Bogdanova et al., 2009). In Middle Cambrian time, thermal relaxation of this rifted margin along with rising eustatic sea levels led to a period of transgressive sedimentation referred to as the Sauk transgression (Sloss, 1963). The oldest Phanerozoic rocks in the Bridger Range were depositing during this period (McMannis, 1955). Throughout Paleozoic time, this dominantly passive margin-style marine sedimentation on the stable shelf of the North American craton gradually transitioned to an active collisional margin, marked by terrane accretion and enigmatic mountain building events on the western margin of Laurentia (Dickinson, 2006). In the Late Devonian–early Mississippian, thrusting and emplacement of allochthons occurred, mostly in Nevada (USA) and nearby regions, in an event referred to as the Antler orogeny (Crafford, 2008). In the Idaho-Montana region, there is little evidence of thrusting during this time (Dickinson, 2004). Nevertheless, the Laurentian platform subsided in Montana and Idaho, forming the elongate Antler foreland basin system (Savoy and Mountjoy, 1995). In Pennsylvanian time, the development of intracratonic block uplifts in Colorado (USA) and surrounding regions resulted from collision of Laurentia’s southern margin with South America and Africa in an event termed the Ancestral Rocky Mountain orogeny (Kluth and Coney, 1981, May et al., 2013). In Montana, far-field stresses from this collision are interpreted to have caused subsidence of the central Montana trough and uplift of a lowland in the Beartooth shelf region (Fig. 1; Maughan, 1990).

In Middle Triassic–Middle Jurassic time, a nascent magmatic arc system developed as the western margin of North America evolved from a region of fringing arcs and interarc oceanic basins to a tectonically consolidated margin associated with east-dipping subduction of the Farallon Plate beneath North America by Late Jurassic time (Dickinson, 2004; Colpron et al., 2007). This resulted in the establishment of the North American Cordillera, a retroarc orogenic system extending more than 6000 km from southern Mexico to Alaska, USA (e.g., DeCelles, 2004).

A complex and extensive foreland basin system developed during Cordilleran orogenesis as a result of lithospheric flexure (Fig. 1A; Kauffman and Caldwell, 1993; DeCelles, 2004). This basin, referred to as the Cordilleran retroarc foreland basin system or Western Interior Basin, paralleled the developing orogen to the west throughout Late Jurassic into Eocene time and eventually occupied an area of >5 × 106 km2 (Fig. 1A; DeCelles, 2004). Eastward propagation of the fold-and-thrust belt resulted in the uplift and erosion of Proterozoic–Paleozoic strata and Mesozoic foreland basin sediments, recycling this detritus into the foreland basin system (e.g., Fuentes et al., 2011; Laskowski et al., 2013).

During Late Cretaceous–Eocene time, the foreland basin system was partitioned by uplifts exposing Archean–Proterozoic basement along moderately dipping reverse faults (Fig. 1A; DeCelles, 2004). The Bridger Range is one of the northernmost instances of these uplifts. During the Paleocene, northwest-trending basement structures were reactivated with a reverse sense of offset, creating a fault-propagation fold. Proterozoic through synorogenic strata were uplifted in this ancestral Bridger anticline (Lageson, 1989). Beginning in the Paleocene, the Cordillera began to gravitationally collapse along normal faults mainly reactivating fold-thrust belt structures (Constenius, 1996). The west limb of the Bridger anticline was down-dropped along a steep north-south–striking normal fault, forming the Gallatin Valley (Lageson, 1989). The east limb of this anticline now forms the Bridger Range.

Eight samples were collected from siliciclastic intervals of Mesoproterozoic–Cretaceous strata along a roughly east-west strike-perpendicular transect across Sacagawea Peak (Proterozoic–Paleozoic samples) and near Fairy Lake (Mesozoic samples) in the northern Bridger Range, Montana (Fig. 1B). Here, we briefly describe the sedimentary facies and depositional environments of the clastic intervals sampled in this study (Fig. 2).

In our study area, the Belt Supergroup assemblage consists solely of the Mesoproterozoic LaHood Formation (Fig. 2). The age of the LaHood Formation is bracketed by the underlying 1.79 Ga basement in the Little Belt Mountains and the intrusion of the Moyie and Perma sills at ca. 1.45 Ga (Anderson and Davis, 1995; Sears et al., 1998; Ross and Villeneuve, 2003). The LaHood Formation was deposited atop Proterozoic and possibly Archean rock of the Great Falls Tectonic Zone (GFTZ) of central and southwestern Montana (Mueller et al., 2016; Fox, 2017). Thickness estimates for the LaHood Formation range from 2400 to 3500 m (Schmidt et al., 1988; Nilsen, 1991), with ∼2000 m exposed in the Bridger Range (McMannis, 1955, 1963). The LaHood Formation is divided into two members, the Corbly Canyon Conglomerate and the Reese Creek Sandstone. The latter interfingers with the Johnson Canyon Carbonate member of the lower Newland Formation to the north (Fox, 2017). Sample “Yla” was collected from the Reese Creek Sandstone, which locally consists of a reddish-purple slightly pebbly coarse arkose. This sandstone is poorly sorted with visible lithic fragments and abundant feldspar. Bedding is generally massive (>2 m). Down-section of the sample location, beds <2 m thick are interbedded with thinly bedded coarse sandy laminated shale. Graded bedding, soft sediment deformation, and flute casts were documented where these two lithofacies are interbedded by Fox (2017), who inferred these sediments were deposited by the distal frontal lobes of debris flows. Elsewhere, facies within the LaHood are highly variable and have been interpreted to comprise alluvial fan to submarine basin-plain facies (Nilsen, 1991).

Paleozoic stratigraphy in the Bridger Range is Middle Cambrian to lower Pennsylvanian in age and is dominated by thick carbonate sequences punctuated by deep-marine to eolian siliciclastics (Poole and Sandberg, 1977; Maughan, 1993). The Middle Cambrian Flathead Sandstone was deposited in a shallow-marine to marginal-marine environment atop the LaHood Formation in the northern portion of the Bridger Range and atop Proterozoic crust in the southern portion following a global temporally extensive hiatus (Sloss, 1963). Fluvial facies of the Flathead Sandstone have been documented in northern Wyoming, where it is divided into two units with inferred tidal flat–fluvial and shallow marine depositional settings, respectively (e.g., Bell, 1970). This unit thickens to the northwest of our study area and thins to the east as strata onlap the craton and to the south across a SW-NE–trending topographic high in southern Montana (Forhan, 1949). South of this topographic high, it thickens again into the Big Horn Basin of northern Wyoming (Bell and Middleton, 1978). In the northern Bridger Range, the Flathead Sandstone is a white feldspathic sandstone containing sparse cross-beds and quartz pebbles outcropping on the western side of the range (McMannis, 1955). Together with correlative Neoproterozoic to Cambrian sandstone units, the Flathead Sandstone forms the basal Sauk transgressive sequence that marks the onset of passive margin sedimentation in western North America (Sloss, 1963; Yonkee et al., 2014). Sample “Cf” was collected from basal beds of unit.

The Devonian Sappington Formation is stratigraphically situated between the unconformably underlying Devonian Three Forks Formation and the overlying Mississippian Lodgepole Formation (Fig. 2). It was deposited in a small equatorial basin covered by a shallow epicontinental sea in southern Montana during the Late Devonian and early Mississippian (Peterson, 1981). This basin likely formed in response to the emerging Antler orogen west of the study area during the Late Devonian (Dorobek et al., 1991). The Sappington Formation is regionally subdivided into an upper organic-rich mudstone member, a middle dolomitic siltstone member, and a lower organic-rich mudstone member, recording deposition in a restricted offshore marine to tidal flat-deltaic environment (Phelps et al., 2018). Based on stratigraphic similarities and correlation of these facies, it has been suggested as an outcrop analogue to the oil producing Bakken Shale of the Williston Basin to the east in North Dakota, USA (Phelps et al., 2018). Sample “Ds” was collected from the Sappington Formation and analytical results were published in Chapman and Laskowski (2019). The sample is included here due to its location along the present transect.

In the Bridger Range, the Big Snowy Group unconformably overlies a karst surface with local relief of as much as ∼30 m developed at the top of ∼500 m of carbonate Madison Group strata (Fig. 2; Shean, 1947). To the east, this unit conformably overlies the Charles Formation or its equivalent in the Madison Group (Maughan and Roberts, 1967). The Big Snowy Group was deposited in a trough that extended from central Montana to the Williston Basin in North Dakota (Gardner, 1959; Maughan and Roberts, 1967). Deposition in the Bridger Range took place on the edge of this trough, and consequentially this group is absent south of the Pass Fault to the south of our study area (McMannis, 1955). The Kibbey Sandstone is the most widespread member of the Big Snowy Group, locally consisting of a brick red sandstone deposited in a shallow marine environment with a maximum thickness of ∼75 m in the Bridger Range (McMannis, 1955). This formation has been dated as Chesterian on the basis of fossil evidence (Easton, 1955). Sample “Mk” was collected from this unit.

The Amsden Formation unconformably overlies the Big Snowy Group in the Bridger Range (Fig. 2; McMannis, 1955). Here, we adapt the nomenclature most commonly used in the Bridger Range of “Amsden Formation,” following McMannis (1955), Guthrie (1984), and Vuke et al., (2007), but acknowledge that other regions of Montana have elevated the Amsden to group status (e.g., Bottjer, 2017; Stanton and Silverman, 1989). The Amsden Formation consists of a basal sandstone overlain by red siltstones and upper interbedded dolomite and quartz sandstone. The Amsden Formation and lithostratigraphically correlative rocks were deposited along a laterally extensive portion of the passive margin, including Idaho, southwestern and central Montana, and northern Wyoming during a marine transgression linking these basins (Sando et al., 1975). The upper carbonate unit of this formation has yielded Morrowan conodonts in the Bridger Range (Guthrie, 1984). The overlying Quadrant Formation is separated from the underlying Amsden Formation by a major sequence boundary (Maughan and Roberts, 1967). Beds of the upper Quadrant Formation sampled near Three Forks, Montana (∼45 km west of our study area) have yielded Desmoinesian fusulinids (Thompson and Scott, 1941), suggesting this unit is locally lower(?)–middle Pennsylvanian in age. This unit constitutes a progradational sequence consisting of a lower dolostone and sandstone unit and an upper cross-stratified quartz arenite unit, interpreted to represent erg deposits (Maughan and Roberts, 1967). This unit thickens substantially eastward into the Snowcrest trough extending into central Montana (Saperstone, 1986). The Quadrant Formation was deposited in a subtidal, prograding to eolian, environment and is lithostratigraphically correlative with the Tensleep Sandstone of the Bighorn Basin (Saperstone and Ethridge, 1984; Saperstone, 1986). Sample “IPMa” was collected from sandstone interbedded in the Morrowan upper carbonate member of the Amsden Formation. Sample “IPq” was collected from lower(?)–middle Pennsylvanian basal sandstone beds of the Quadrant Formation.

The Mesozoic stratigraphy in the Bridger Range was deposited during Upper Jurassic to Cretaceous time (Fig. 2). In southwestern Montana, an unconformity of ∼150 m.y. separates the Pennsylvanian Quadrant Formation from the Upper Jurassic Ellis Group. The Ellis Group consists of ∼30–110 m of marine strata subdivided into three members—the Sawtooth Formation, Rierdon Formation, and Swift Sandstone (McMannis, 1955). The Swift Sandstone consists of a basal conglomeratic zone followed by locally fossiliferous yellow calcareous cross-bedded sandstones—deposited in a shoreface environment—marking a marine regression (McMannis, 1955; Fuentes et al., 2011). The Ellis Group is irregularly distributed and displays abrupt lateral facies changes and local unconformities (Parcell and Williams, 2005). Reported ages for the Ellis Group range from Bajocian–Bathonian to late Oxfordian–Kimmeridgian and the temporal intervals represented by internal unconformities have been disputed within the literature (e.g., Cobban, 1945; Imlay, 1980; Porter, 1989; Parcell and Williams, 2005; Fuentes et al., 2009, 2011). Recent detrital zircon and palynology work by Fuentes et al. (2011) constrains the age of the Ellis Group in northern Montana to Bajocian–mid-Oxfordian. The Ellis Group is temporally correlative with the early foreland basin deposits of the Fernie Formation of southwestern Canada (Poulton et al., 1994) and the Gypsum Spring and Sundance formations of the Bighorn Basin (Parcell and Williams, 2005). Sample “Je” was collected from the Swift Sandstone member of the Ellis Group.

The Late Jurassic Morrison Formation overlies the Ellis Group. In Montana, the Morrison Formation is conformable with the underlying Swift Sandstone (Demko et al., 2004). In the Bridger Range, the Morrison Formation consists of ∼35–135 m of variegated shale and mudstone with interbedded calcareous sandstones likely representing low-energy fluvial to lacustrine deposition (McMannis, 1955). The Morrison Formation in northern Montana has been dated as Oxfordian to Kimmeridgian using palynology (Fuentes et al., 2011), which agrees with previous work in Utah, Colorado, and Wyoming (Kowallis et al., 1998; Litwin et al., 1998; Turner and Peterson, 2004). The Morrison Formation is correlative with the upper Fernie Formation and lower Mist Mountain Formations of southwestern Canada (Poulton et al., 1994). Sample “Jm” was collected from the Morrison Formation.

A regional unconformity separates the Lower Cretaceous Kootenai Formation from the underlying Morrison Formation. Locally, this unconformity is irregular and channelized (McMannis, 1955). The Kootenai Formation was deposited in a fluvial-to-lacustrine environment within the foredeep depozone of the foreland basin (DeCelles, 1986). In the Bridger Range, the basal member of the Kootenai Formation consists of trough cross-bedded coarse salt-and-pepper sandstone composed of quartz and black chert clasts overlying chert and quartzite pebble conglomerate (DeCelles, 1986). This member, from which sample “Kk” was collected, locally contains chert pebbles and wood fragments and is overlain by quartz arenite. The age of this member as well as its relationship to the upper members of the Kootenai Formation have been the subject of debate, with researchers working in central Utah suggesting continuous deposition between the Morrison Formation and the overlying conglomeratic beds (Roca and Nadon, 2007). However, recent detrital zircon data from northern Montana yielded Hauterivian ages, suggesting a large unconformity separates the Morrison and basal Kootenai formations (Fuentes et al., 2011). The upper Kootenai Formation records a transition from high-energy stream environments to lower-energy fluviolacustrine settings (Suttner, 1969). Detrital zircon samples from these beds have yielded euhedral zircons of early Albian age, inferred by Fuentes et al. (2011) to reflect syndepositional volcanism. The overlying basal Blackleaf Formation has been dated as late Albian–early Cenomanian (Cobban and Kennedy, 1989), constraining the maximum age of the Kootenai Formation.

U-Pb Detrital Zircon Geochronology

Zircon crystals from the eight samples collected in the Bridger Range were extracted by crushing and grinding, followed by separation with water, heavy liquids, and a Frantz magnetic separator. Zircon grains were retained in the final heavy mineral fraction. A large split of recovered zircon grains was poured or hand mounted on a 1” epoxy mount. Fragments of Sri Lanka, FC-1, and R33 zircon crystals were mounted to serve as standards. Approximately 20 μm of epoxy was sanded from the pucks to expose zircon grains. The pucks were then polished. Backscattered electron images of mounts were generated with a Hitachi 3400N scanning electron microscope at the Arizona LaserChron Center at the University of Arizona (Tucson, Arizona, USA) to provide a guide for locating analysis pits in the core of zircon crystals free from excessive zoning or radiation damage. Analysis spots were placed using automated spot-picking software to reduce human-introduced bias.

Zircon U-Pb ages were obtained using laser ablation–inductively coupled plasma–mass spectrometry at the Arizona LaserChron Center (Gehrels et al., 2006, 2008; Gehrels and Pecha, 2014; Pullen et al., 2018). These analyses consisted of ablation of zircon with a Photon Machines Analyte G2 excimer laser equipped with HelEx ablation cell using a spot diameter of 20 μm. The ablated material is carried in helium into the plasma source of an Element2 high-resolution–inductively coupled plasma–mass spectrometer, which sequences rapidly through U, Th, and Pb isotopes. With the laser set to an energy density of ∼5 J/cm2, a repetition rate of 8 Hz, and an ablation time of 10 seconds, ablation pits are ∼12 μm in depth. Each analysis consists of 5 s on peaks with the laser off (for backgrounds), 10 s with the laser firing (for peak intensities), and a 20 second delay to purge the previous sample and save files.

Following analysis, data reduction was performed with an in-house Python decoding routine and an Excel spreadsheet (E2agecalc). 204Hg was subtracted from the 204 signal to yield 204Pb intensity. Common Pb was corrected for using the measured 206Pb/204Pb and the assumed composition of common Pb (Stacey and Kramers, 1975). In-run fractionation was corrected for by comparing measured and known ratios for zircon standards then applying these fractionation factors to unknowns using a sliding-window average. Internal uncertainty was calculated by summing 206Pb/238U, 208Pb/232Th, 206Pb/207Pb, and 206Pb/204Pb measurement uncertainties with any overdispersion factor. External uncertainty was calculated from the scatter of standard analyses, uncertainty in ages of standards, uncertainties in the composition of common Pb, and uncertainties in the decay constants for 235U and 238U. Ages <900 Ma reflect the 206/238 ratio while ages >900 Ma reflect the 206/207 ratio.

Approximately 300 grains were analyzed from most samples. One-hundred and five grains were analyzed from the Kibbey Sandstone sample due to poor zircon yield. All new analytical data and a detailed list of concordance filters are reported in Supplemental item 11. The resulting interpreted ages were prepared as normalized age-probability kernel density estimate (KDE) diagrams (Fig. 3) using the DetritalPy python script (Sharman et al., 2018). These diagrams show each age as a normal distribution with a fixed “kernel” with a bandwidth of 10 m.y. and sum all kernels from a sample into a single curve. This bandwidth was chosen based on qualitative comparisons of KDEs with a bandwidth of 5, 10, and 20 m.y. A bandwidth of 10 m.y. best balances oversmoothing and undersmoothing and most intuitively represents the detrital zircon age distribution. A comparison of KDEs using these bandwidths is available in Supplemental item 2 (see footnote 1). Curves are normalized according to the number of constituent analyses, such that each curve contains the same area, and then stacked to facilitate the comparison of age distributions between samples composed of variable numbers of analyses. KDEs are used in lieu of probability density plots, another common visualization technique within detrital zircon geochronology, to avoid counter-intuitive results from variable precision and number of analyses (Vermeesch, 2012). Maximum depositional ages (MDAs) were calculated using the DetritalPy python script (Sharman et al., 2018) based on the methods of Dickinson and Gehrels (2009) and Coutts et al. (2019) (Table 2).

Multidimensional Scaling (MDS) Analysis

To quantitatively compare detrital zircon populations, multidimensional scaling (MDS) was implemented using a MATLAB program, DZmds, from Saylor et al. (2018) (Fig. 4). This approach is a variety of principle component analysis adapted to facilitate the comparison of multiple detrital zircon data sets (Vermeesch, 2013; Che and Li, 2013; Sharman and Johnstone, 2017). This method allows for the visualization of the relative dissimilarity between detrital populations by converting intersample dissimilarity into disparity via linear transformation and rearranging this disparity in Cartesian space, while seeking to minimize stress, or the misfit between disparity and Cartesian distance. We input kernel density estimates with a kernel of 10 m.y. generated from nine samples collected from the Bridger Range (8 from this study, 1 from Chapman and Laskowski, 2019), as well as from a source compilation of 53 samples from Middle Cambrian–Lower Cretaceous strata in the Northern Rockies region; a total of 10,824 individual detrital zircon ages composed this input data set. Samples were grouped with unit(s) deposited during a similar time interval in the Bridger Range. Cross-correlation coefficients were calculated from these distributions and used as a metric of dissimilarity. This statistic was used because it is more sensitive to the number and relative proportion of age modes within a detrital age distribution compared to other popular statistics, such as the Kolmogorov-Smirnov (K-S) test distance statistic (Vermeesch, 2013; Saylor et al., 2018; Sundell et al., 2018). The ideal number of dimensions for the MDS plot output was selected based on the inflection point in a scree plot, which plots stress against number of dimensions. The goodness of fit of the model output was assessed using the stress statistic. The grouping of samples in MDS space suggests a similar detrital zircon age distribution (Sharman and Johnstone, 2017) and can be used to infer common provenance regime between samples (e.g., Quinn et al., 2018).

Non-Negative Matrix Factorization (NMF) Analysis

Inverse modeling was implemented to characterize sediment source age distributions. This MATLAB program—DZnmf, from Saylor et al. (2019)—utilizes non-negative matrix factorization (NMF) to unmix the detrital zircon age distribution of input samples into “synthetic” source populations, whose age distributions are statistically factorized. This program implements an optimum factorization rank routine to determine the optimum number of synthetic sources capable of reconstructing the age distributions of a set of sink samples. NMF has proven its ability to help unravel complex sediment mixing patterns by leveraging differences in the age spectra between samples to highlight geologically meaningful age populations whose abundance differs in a correlated fashion between samples—these ages are interpreted to compose a source (Sharman and Johnstone, 2017; Saylor et al., 2019; Chapman and Laskowski, 2019; Laskowski et al., 2019, Leary et al., 2020).

The algorithm takes sink sample age distributions as its input. We input kernel density estimates from our compilation of 62 samples (8 new, 54 compiled; Table 1; analytic results of compiled samples in Supplemental item 3; see footnote 1). The NMF algorithm outputs potential source age distributions (Fig. 5) along with a weighting function for each sink sample, describing the relative contribution from each synthetic source required to reconstruct the age distribution of the sink sample (Supplemental item 4 [see footnote 1]; pie charts in Figs. 68). During modeling, sink samples were grouped by provenance regimes, as determined by clustering in MDS space (Fig. 4). Grouping samples into common provenance regimes produced the most intuitively geologically meaningful reconstructed sources, best accounting for recycling of sedimentary rocks and increasing complexity of age distributions up section (Schwartz et al., 2019). This grouping achieved a higher cross-correlation coefficient—a measure of goodness of fit between reconstructed and input sources—versus modeling all samples together.

The optimum number of factorized sources is determined by the NMF algorithm as a point of inflection in the decrease of the final residual versus the number of sources modeled (i.e., factorization rank, Saylor et al., 2019). To identify the optimum number of sources, we modeled 2–20 sources for each sample group and identified this point of inflection for inputs of 10 m.y. and 20 m.y. kernel bandwidths.

Zircons with ages representing most major North American basement provinces have been shown to persist up section following their arrival in the Northern Rocky Mountains due to sedimentary recycling (Schwartz et al., 2019). Thus, the presence of zircon having an age coincident with the age of a specific region of North American basement alone is not diagnostic of unique sediment source areas or dispersal pathways linking surface exposures of these basement provinces to depocenters. This has been recognized by previous authors, who have exploited a number of statistical tools, such as K-S testing and/or MDS analysis (e.g., Laskowski et al., 2013; Saylor and Sundell, 2016; Sundell et al., 2018; Quinn et al., 2018) and traditional provenance indicators, such as petrographic modal analysis (e.g., Fuentes et al., 2011; Sundell et al., 2018) to compliment traditional detrital zircon techniques and group samples of similar provenance. However, this categorical assignment of provenance interpretation groups assumes that provenance behaves as a discrete variable. The weighing factors assigned by NMF provide a means by which to honor the integrated nature of provenance—existing as a mixture of several end-member sources, each itself a mixture of igneous and detrital zircon age populations—and provide a more rigorously quantitative basis for interpreting regional sediment-dispersal patterns.

This method is less biased by a priori assumptions of the multitude or age distribution of potential sources that affect forward modeling approaches similarly utilized to statistically determine relative contributions by sources to sediment sinks (e.g., Sundell and Saylor, 2017). Forward modeling uses age distributions from modern rocks as model inputs to approximate ancient sediment sources, which assumes that age distributions of ancient sources are reflected in modern rocks. This is problematic given that many sources are inherently poorly preserved in the geologic record and sedimentary sources are frequently heterogeneous through space and time. However, forward modeling approaches have been shown to be effective in instances where the potential sources are well known and characterized (e.g., Sundell and Saylor, 2017).

NMF is effective in the present study because our data set is large, dissimilarity among sink samples is high, and sink samples are well characterized (average n = 176). These criteria are necessary for effective use of this approach (Saylor et al., 2019). However, several considerations may affect the results: (1) geologic context is required to differentiate primary, mixed sources from recycling of (meta-) sedimentary sources; (2) the algorithm will factorize the most basic elements of sink samples—therefore, mixed sources may be subdivided into their constituent parts; and (3) while the algorithm will converge on the most statistically likely solution, NMF inversions are inherently non-unique (Saylor et al., 2019). While results should be evaluated in light of independent geologic evidence, analysis of age distributions of synthetic sources and spatial patterns of source weights provides a powerful new tool from which to glean insight into complex sediment dispersal systems.

Isotopic Results

Detrital zircon ages for individual samples collected from the Bridger Range are plotted on normalized stacked kernel density estimate diagrams with a kernel width of 10 m.y. (Fig. 3). Age peaks discussed herein form a graphical peak consisting of four or more grains. Peaks were located with the DetritalPy peak finder routine (Sharman et al., 2018). MDAs are reported in Table 2.

The sample from the Mesoproterozoic LaHood Formation is dominantly composed of Archean grains—centered around a peak at 3276 Ma (Fig. 3). Additionally, Paleoproterozoic grains are present—clustered around a peak at 1794 Ma. The Middle Cambrian Flathead Sandstone yielded an age spectrum partially mirroring that of the underlying LaHood Formation, with prominent age-probability peaks centered around 3251 Ma and 1782 Ma similar to those of the LaHood (Fig. 3). However, the Late Archean (ca. 2.6–3.0 Ga) age group is much more pronounced within the Cambrian. Later Paleozoic samples can be broadly characterized by a shift to dominance of ca. 400–460 Ma and ca. 990–1900 Ma zircon ages, with minor additional Archean-aged grains (Fig. 3). The youngest group of ages within the Devonian Sappington Formation forms bimodal peaks at 403 Ma and 459 Ma with an additional cluster of ages around 1056 Ma. A spread of ages between ca. 1.3–1.85 Ga and minor peaks of ca. 2.6–2.8 Ga grains are also present. Though separated by more than 500 m of carbonate, the Mississippian Kibbey Sandstone has a zircon age spectrum very similar to that of the Devonian Sappington Formation. The Archean peak present within the Sappington Formation is subdued within the Kibbey Sandstone. The lower Pennsylvanian Amsden and lower(?)–middle Pennsylvanian Quadrant formations have generally similar age distributions which broadly reflect the patterns of the previously discussed Paleozoic samples; these formations can be differentiated from the Sappington Formation and Kibbey Sandstone by their sharp ca. 430 Ma peak, as well as the presence of ca. 1.8–2.0 Ga grains that are sparse to absent in underlying Paleozoic units.

Mesozoic samples from the Bridger Range are characterized by the appearance of detrital zircon ages interpreted to be very close to depositional ages (Fig. 3). The Upper Jurassic Swift Sandstone is dominated by two prominent age-probability peaks centered at 165 Ma and 233 Ma, respectively. Subsidiary detrital zircon age populations in the Swift Sandstone generally reflect ages represented in the underlying Paleozoic strata; however, there is a notable increase of ca. 500–700 Ma and 1.4 Ga grains. The Upper Jurassic Morrison Formation, while less dominated by Triassic–Jurassic ages than the Swift Sandstone, also contains abundant ages centered around 157 Ma and 270 Ma. Grains of 500–700 Ma ages are abundant, with a broad peak centered at 596 Ma. The distribution of >700 Ma ages generally follows the trends of the Devonian–Pennsylvanian passive margin strata. The Cretaceous Kootenai Formation has a very similar age distribution to the underlying Morrison Formation, with abundant Triassic–Jurassic arc grains with peaks at 158 and 232 Ma. The Proterozoic–Archean age distribution resembles that of the Morrison Formation, with the notable addition of a pronounced peak centered at 1829 Ma and a higher abundance of ca. 2.5–3.0 Ga grains.

MDS Analysis Results

Multidimensional scaling (MDS) analysis of Northern Rockies sink samples compiled from each period of interest are shown on a 3-D MDS plot in Figure 4. Inferred trends in provenance within MDS space are plotted as arrows. A stress value of 0.12 was obtained using a 3-D MDS plot and cross correlation values as a metric of dissimilarity. Mesozoic samples form one tight cluster while Pennsylvanian and Devonian samples form a tight pair. Mississippian samples are nearest to Devonian samples but don’t group tightly. Mesoproterozoic and Cambrian samples plot closest to each other, but by themselves.

NMF Analysis Results

Non-negative matrix factorization analysis was performed on Cambrian, Devonian–Pennsylvanian, and Mesozoic sample compilations separately based on commonality of provenance regime reflected by grouping of samples in MDS space (Fig. 4). This helps account for the recycling of Paleozoic passive margin strata into the Mesozoic foreland basin system. Our grouping approach achieved the highest cross-correlation, used as a metric of goodness-of-fit, between the age distribution of the input sink samples and those reconstructed by the algorithm. The Mesoproterozoic LaHood Formation was excluded from NMF modeling due to the paucity of analyses in the literature to use as an input compilation of sink samples.

NMF analysis for inputs of 10 and 20 m.y. bandwidths produced similar synthetic sources but differed in their identification of the optimum number of sources (Supplemental item 4). 10 and 20 m.y. bandwidths yielded: either five or four optimum sources for Mesozoic samples, respectively; either five or three optimum sources for Devonian–Pennsylvanian samples, respectively; either four or three optimum sources for Cambrian samples, respectively. However, we interpret the 10 m.y. bandwidth results to be more reasonable when considering geologic context and prior provenance interpretations. For example, a 5-source solution for Devonian to Pennsylvanian (D–IP) samples generates a source with a dominant sharp age peak of ca. 426 Ma (D–IP source 4) and a source with a broader age peak centered at ca. 440 Ma and a very broad peak centered at ca. 1070 Ma (D–IP source 3). These sources are grouped in the 4-source solution generated with an input of 20 m.y. kernel bandwidths. We interpret that splitting these sources is informative and provides information about the switch between provenance from the Franklinian orogen and the Appalachian orogen (see Provenance and Paleogeographic Evolution section below). Therefore, we focus our discussion on the results produced by the model with an input of 10 m.y. kernel bandwidths.

Synthetic sources are plotted as stacked kernel density estimates in Figure 5; crustal provinces of North America, which likely originally supplied the zircon ages represented by these synthetic sources, are plotted as vertical bars. Reference curves from modern rocks which provide an analogue to the detrital zircon age distribution of various hypothesized source regions, compiled from previous studies (references in figure caption) are overlaid as unfilled KDE curves in Figure 5 to facilitate comparison to synthetic sources reconstructed by NMF. Source regions for each synthetic source are discussed below. Weighing functions describing the relative contribution from each synthetic source (Fig. 5) required to best construct the age distribution of each sink sample are shown as pie diagrams in Figures 68. Raw NMF synthetic source age distributions, goodness-of-fit metrics, and weighing functions are provided in Supplemental item 4.

We discuss the synthetic sources generated by NMF modeling of strata of the Northern Rocky Mountains below. Synthetic sources are interpreted to have been major sources of sediment to the Northern Rocky Mountain region but not all are interpreted to have contributed to samples from the Bridger Range. A detailed discussion of the regional provenance evolution incorporating these synthetic sources is found in the Provenance and Paleogeographic Evolution section below.

Cambrian Sources

NMF analysis predicts four synthetic sources for Cambrian samples (Fig. 5), which we interpret to dominantly reflect the recycling of local rock units (Fig. 6). Cambrian source 1 represents detritus sourced from the Lemhi Arch region (Fig. 1A), based on the prominent 500 Ma age-probability peak. These grains are nearly unique in the detrital zircon record of North America and were derived from the Beaverhead plutonic suite east-central Idaho, as shown by Lund et al. (2010) and Link et al. (2017). Additionally, the subsidiary ca. 1.7 Ga population within this source is common within strata of the Belt Supergroup that the Beaverhead plutonic suite intrudes (Lewis et al., 2010). Cambrian source 2 is composed of detritus derived from the Wyoming craton, evidenced by the dominant ca. 2.88 Ga age population, characteristic of the ca. 2.82–2.95 Ga Beartooth-Bighorn Magmatic Zone (Mogk et al., 1992; Frost et al., 2006), as well as minor ca. 3.3–3.1 Ga grains, associated with Archean gneisses of the Wyoming Province (Chamberlain et al., 2003).

Source 3 represents detritus from the Cheyenne Belt of southern Wyoming. This source likely requires a major igneous source of zircons based on the unimodality of the ca. 1.78 Ga age peak. Grains of the ca. 1.78 Ga peak are present in the 1.74–1.79 Ga Green Mountain magmatic arc of southeastern Wyoming, northern Colorado, and northern Utah (Fig. 1A; Hills et al., 1968; Premo and Loucks, 2000). The subsidiary ca. 1.86 Ga age population is present within the proximal reworked Archean crust of northern Utah (Whitmeyer and Karlstrom, 2007). While grains of similar ages have been identified in the Big Sky orogeny (Ault et al., 2012; Condit et al., 2015) and Little Belt arc (e.g., Mueller et al., 2002; Foster et al., 2006; Gifford et al., 2014) in southwest-central Montana, Paleoproterozoic zircons in these sources do not form a well-defined peak at ca. 1.78 Ga and are commonly associated with >2 Ga zircons. Additionally, the shoreline transgression southeastward onto the craton largely buried this region by the early Middle Cambrian making this source unavailable during sedimentation of the upper Flathead Sandstone in Wyoming, such as in the Big Horn Basin (May et al., 2013). For these reasons, we favor an interpretation in which the Cheyenne Belt supplied the zircons represented by this synthetic source.

Source 4 represents either detritus recycled out of upper strata of the Belt Supergroup, in the case of samples overlying these rocks, or primary input from the Yavapai/Mazatzal and mid-continent granite provinces to the south, in the case of late Middle to Upper Cambrian samples in Idaho and Wyoming. The bimodal Paleoproterozoic age population that dominates this source represents grains originally derived from the GFTZ (Foster et al., 2006; Harms et al., 2006; Labadie, 2006; Gifford et al., 2014) and sourced from reworked Archean crust and Yavapai/Mazatzal arcs (Hoffman, 1988; Whitmeyer and Karlstrom, 2007). A reference curve showing the average age distribution of Belt Basin strata—compiled from a data set of 4314 ages from 61 samples—is provided for comparison (Fig. 5). The offset of the Paleoproterozoic peak of the Belt Basin strata versus the synthetic source may reflect a transition toward greater abundance of ca. 1.8 Ga ages up-section, as documented in LaHood Formation (Fox, 2017) and discussed below. Vast amounts of strata of the Belt Supergroup were likely eroded during the development of the Great Unconformity. These strata may have had a larger proportion of these ca. 1.8 Ga ages, which were then incorporated into Cambrian units. However, there is little evidence for this up-section transition in the Belt Basin east of the Helena Embayment.

Devonian–Pennsylvanian Sources

NMF analysis reconstructs five potential synthetic source age distributions for Devonian–Pennsylvanian samples. Overall, we interpret these sources to reflect a shift to the dominance of more distal northern, northeastern, and southeastern sources (Figs. 5 and 7); additionally, recycling of previously deposited Neoproterozoic–Paleozoic rocks likely contributed to source age distributions. Source 1 reflects detritus derived dominantly from recycling of Ordovician strata. The age distribution of this synthetic source matches closely with that of preserved North American Ordovician strata (Fig. 5B; Gehrels and Pecha, 2014). The dominant ca. 1.8–2.0 Ga and ca. 2.5–2.78 Ga ages were originally sourced from the Peace River Arch of the northern Canadian Shield, which exposed Archean–Paleoproterozoic crust in the Ordovician; during this time, this detritus dominated the western passive margin of Laurentia from Canada through Sonora (Gehrels and Pecha, 2014). The Ordovician inundation and subsequent later Paleozoic draining of the Canadian shield may have facilitated the erosion and recycling of this strata.

Source 2 represents detritus derived from Proterozoic terranes associated with mid-Paleozoic orogenesis. This likely reflects arcs built on Proterozoic crust outboard of the Laurentian cratonic margin. The similarity between the age distribution of strata of the eastern Klamath terrane (Figs. 1A and 5; Grove et al., 2008) and this source indicates this or analogous terranes, sutured to the western coast of Laurentia beginning in Devonian time, may be good candidates. Much of these terranes have been overprinted by subsequent Cordilleran magmatism and deformation, but the Klamath or an analogous terrane may have been juxtaposed to the western margin of Laurentia near the latitude of Montana in the middle Paleozoic (Beranek et al., 2016 and references therein). Alternatively, in Pennsylvanian time, this source may reflect detritus from rocks of the Yavapai/Mazatzal Province exposed in uplifts of the Ancestral Rockies orogen and routed northward across mid-continent granites. Zircon Hf analysis, such as that conducted by Beranek et al. (2016), would assist in differentiating between these two potential source regions likely represented by this synthetic source age distribution.

Source 3 is dominated by Appalachian and Grenvillian ages, with subsidiary populations of Mesoproterozoic mid-continent granite/rhyolite province ages (Fig. 5B). The age distribution of this synthetic source closely resembles the age distribution of Devonian through Permian strata that accumulated in the Appalachian orogen (Fig. 5B; Gray and Zeitler, 1997; McLennan et al., 2001; Eriksson et al., 2004; Thomas et al., 2004; Becker et al., 2005, 2006; Park et al., 2010; Gehrels et al., 2011). This suggests routing of sediment from the high Appalachian/Grenville orogens along the east or southeast margins of Laurentia through the U.S. midcontinent.

Source 4 has a very prominent Appalachian peak with subdued subsidiary age populations at ca. 1.1 Ga and 1.6–1.8 Ga. We interpret this to represent sediment sourced from the Franklinian and Caledonian orogens along the northeastern margin of Laurentia; there, rocks of Appalachian age are directly juxtaposed against the Archean craton (Whitmeyer and Karlstrom, 2007), resulting in Appalachian ages present in a larger abundance than Grenvillian ages. This source approximates the Mesoproterozoic and younger portion of the age spectrum of Devonian through Cretaceous strata which accumulated along the Arctic margin in the Franklinian-Caledonian clastic wedge (Fig. 5B; McNicoll et al., 1995; Gehrels et al., 1999, 2011; Røhr et al., 2010; Beranek et al., 2010). The Paleoproterozoic ages present in the Franklinian wedge were likely derived from recycling Ordovician strata (Beranek et al., 2010). These ages are accounted for by source 1, which is generally coupled to this source during primary northern input in the Devonian (Fig. 7A). Source 5 (Fig. 5B) was likely derived from the clastic wedge of the Antler orogeny. The age distribution of this synthetic source is a close match to that of strata that accumulated on top or adjacent to allochthons of the Antler orogeny in northern Nevada (Fig. 5; Gehrels and Dickinson, 2000). The erratic peaks of the reference curve and mismatch in the Mesoproterozoic portion of the age distribution may be an artifact of the small sample size composing the reference curve (n = 52).

Jurassic–Cretaceous Sources

A shift to the dominance of western source terranes is interpreted from the five synthetic source age distributions generated for the Mesozoic (Figs. 5C and 8). The arrival of depositional-aged grains signifies the importance of Cordilleran arcs and terranes while heavily mixed Paleozoic and older age populations within sources reflect recycling of previously deposited passive margin strata. Source 1 (Fig. 5C) dominantly represents the latter, with an age distribution mirroring a composite of Paleozoic strata of Wyoming and Montana (Fig. 5C, this study and references in Table 1) and Mesozoic Eolianites of southwestern North America (Fig. 5C; Laskowski et al., 2013 and references therein). This reflects uplift and subsequent erosion of these units by eastward propagation of the fold-thrust belt (e.g., Fuentes et al. 2011). Source 2 (Fig. 5C) similarly shows recycling of passive margin strata. The dominance of a broad ca. 1.8 Ga peak is characteristic of the passive margin strata of Canada, again evident in the compilation of Gehrels and Pecha (2014) (Fig. 5C). Mismatches in the Mesoproterozoic–Paleozoic portion of the age spectrum may result from these ages being grouped with source 1, which is generally coupled with this source. Additionally, a younger population of grains coincident with the Coast Range arc is present (Gehrels et al., 2009; Cecil et al., 2018). This suggests routing of detritus from this arc through the uplifted and eroding strata of the Canadian passive margin.

Source 3 (Fig. 5C) dominantly represents detritus from the Triassic Wallowa and Olds Ferry terranes of the Blue Mountains province of eastern Oregon and west-central Idaho (ca. 255–215 Ma) (Fig. 1A; Armstrong et al., 1977; LaMaskin et al., 2011; Schwartz et al., 2011; Gaschnig et al., 2017) and Triassic Intermontane terranes of northern Washington (USA)/southern British Columbia (Canada) (ca. 180–350 Ma) (Coney and Evenchick, 1994; Roback and Walker, 1995). Input from the Triassic Mojave arc (Busby-Spera et al., 1990) during Morrison time is also likely represented by this source (e.g., sample 16; Laskowski et al. 2013). Source 4 (Fig. 5C) is composed of detritus sourced directly from volcanism or reworked from tuffs of Cordilleran arcs. This source is almost entirely composed of ages defining a peak at ca. 160 Ma. This age coincides with the Jurassic flare-up experienced by the Sierra Nevada and Coast Range arcs (Fig. 1A), which had a strong component of extrusive volcanism relative to intrusive plutonism within the Sierra Nevada arc (Paterson and Ducea, 2015). This source likely represents ash-fall detritus in depositional-aged samples (e.g., Oxfordian Swift Sandstone). In Kimmeridgian–Cretaceous time, this source may represent a drainage system linking the actively eroding Jurassic plutons of the Sierra Nevada arc to the foreland basin system or reworking of earlier tuffs. Source 5 (Fig. 5C) has a similar Jurassic age peak, correlative with the Coast Range and Sierra Nevada arcs. However, this source includes ages that reflect detritus recycled from the eroding North American passive margin strata. This represents a drainage system routing arc detritus, likely derived from plutons and volcanic deposits proximal to the arc, through the eroding passive margin into the foreland basin. The temporal overlap makes differentiating arc sources difficult, but spatial patterns, as well as the presence of ages associated with the U.S. passive margin, suggest primary derivation from the southern Coast Range arc or northern Sierra Nevada arc.


NMF sources centrally guide our interpretation of Cambrian–Cretaceous provenance. However, due to the paucity of data on the LaHood Formation, we qualitatively interpret Mesoproterozoic provenance based on comparisons to the zircon age distributions of modern rock exposures that are likely representative of ancient sediment sources. Our results from the Reese Creek Sandstone member of the Mesoproterozoic LaHood Formation (Fig. 3) further support the interpretations of previous studies of the LaHood Formation in the Bridger Range (Guerrero et al., 2016; Fox, 2017), Jefferson Canyon (Ross and Villeneuve, 2003; Mueller et al., 2016), and Horseshoe Hills (Guerrero et al., 2016; Mueller et al., 2016). Archean ages are dominant in our sample, with a major peak defined by >100 analyses centered at 3.27 Ga. These zircon ages correspond with the age of both the Wyoming Province to the south and Medicine Hat Block to the north.

Additionally, a population of Eoarchean grains is present, centered around two minor peaks at 3.82 Ga and 3.87 Ga, reflecting recycling of Archean sedimentary rocks. A likely analogue to these sedimentary sources exists in Archean quartzites from the Beartooth, Ruby, and Tobacco Root uplifts of the Northern Wyoming Province. Because these quartzites, whose Archean age distribution (Mueller et al., 1998) mirrors that of the LaHood Formation, were deposited before the suture of the Wyoming and Medicine Hat blocks, we favor provenance within the proximal Wyoming Block. This interpretation is supported by the presence of a subsidiary population of ca. 2.8 Ga grains, likely reflecting zircons of the Beartooth-Bighorn Magmatic Zone (i.e., ca. 2.81 Ga) of the northern Wyoming Province (Mogk et al., 1992). This detritus was likely supplied by rapid down-dropping along the Perry Line, proving an abundance of detritus derived proximally from the south.

The second major population of detrital zircon ages within the LaHood Formation forms a broad peak centered around 1.79 Ga. Grains of this age were likely sourced from Paleoproterozoic sources along the Great Falls Tectonic Zone (GFTZ). Metamorphic zircons from the Big Sky orogeny associated with the GFTZ in southwestern Montana have yielded ages of 1.73–1.75 Ga (Ault et al., 2012; Condit et al., 2015) along with dominantly Archean ages (Weyand, 1989; Mogk, 1990). Though these zircons may be slightly too young to account for the ca. 1.79 Ga age peak in our sample from the LaHood Formation, high-grade tectonism occurred in the northwestern portion of this orogeny, exposed now in the Highland Mountains, at ca. 1.81–1.78 Ma. Therefore, we interpret that the Big Sky orogeny cannot be ruled out as a source for the ca. 1.79 Ga zircons in our sample from the LaHood Formation. The Little Belt magmatic arc of the GFTZ may have also supplied zircons of these ages. Zircons in Paleoproterozoic crust of the Little Belt Mountains, Montana (Mueller et al., 2002) and granitoid xenoliths from proximal to the Grass Range, Montana (Gifford et al., 2014) have yielded ages from ca. 1.73–1.94 Ga. Though exposure of this arc is limited and the ages of its zircons have not been robustly characterized, it was likely active at ca. 1.79 Ga. The ca. 1.79 Ga age population is absent or subdued in several samples of the LaHood from previous studies in the Bridger Range (Guerrero et al., 2016; Fox, 2017) and Jefferson Canyon (Mueller et al., 2016). When present in samples from these studies, it is centered at ca. 1.86–1.84 Ga. This may be an artifact of the small number of grains analyzed per sample in these studies (n = 6–51) or reflect a shift toward Paleoproterozoic provenance up-section, as suggested by Fox (2017).

A minor peak of ca. 2.55 Ga grains present within our sample supports a dominantly southern source. This age population is uncommon in the detrital record of Laurentia, but is present in the Prichard Formation and Ravalli Group of the western Belt Basin. Non-Laurentian sources of these zircons have been proposed in the western Belt Basin, including Northern China, India, Antarctica, and Australia (Anderson, 2017). The Prichard Formation and Ravalli Group also contain an abundance of grains in the range 1.49–1.61 Ga, which are interpreted to be non-Laurentian due to the well documented North American Magmatic Gap of that time interval (Ross and Villeneuve, 2003). However, these ages are absent in our sample from the LaHood Formation. This absence, coupled with whole-rock Nd results from the LaHood Formation suggestive of Wyoming/Hearn sources (Frost and Winston, 1987), suggest that an entirely Laurentian source for the LaHood Formation is likely.

Grains of this ca. 2.55 Ga age are associated with several regions of the Wyoming Province, including the North Snowy Block (Mogk et al., 1988), anatectic melts of the Tobacco Root Mountains (Mueller et al., 2004), and low-volume magmatism in Gallatin and Teton ranges (Mogk et al., 1992; Zartman and Reed, 1998). Additionally, these ages are present in higher abundances within the southern terranes of the Wyoming Province (Frost et al., 1998). Rare earth element patterns and potassium feldspar abundances within the LaHood provide further evidence for a component of more distal detritus (Fox, 2017).

High topography was likely present in the southern reaches of the Wyoming Province during deposition of the LaHood Formation, providing input of detritus, including ca. 2.55 Ga zircons, from more distal regions of the craton. Grains of this age are dominant in the Newland Formation, which interfingers with the upper LaHood Formation (Anderson, 2017). This shift to carbonate sedimentation coupled with the relative increase in far-traveled detritus indicates that rapid down-dropping along the Perry Line may have ceased toward the end of deposition of the LaHood Formation, resulting in a cessation of the input of coarse detritus composed of sediment from proximal Archean blocks. This indicates a subsidence pattern typical of an aulacogen (Einsele, 2013)—fast subsidence at the onset of basin formation, reflected by rapid down-dropping along the Perry Line delivering coarse detritus from the proximal craton, followed by slowing subsidence, reflected by the cessation of proximal detritus supply.

In consideration of the paucity of published large-n U-Pb data sets from the LaHood Formation and difficulty correlating it to other units of the Belt Basin, basin-wide paleogeographic reconstruction is not attempted. However, trends within the LaHood Formation are evident upon comparison with previous studies—a higher of abundance of Paleoproterozoic grains is present within eastern exposures of the LaHood Formation (Guerrero et al., 2016; Fox, 2017; this study) compared to western exposures (Ross and Villeneuve, 2003; Mueller et al., 2016). This may reflect the proximity of eastern exposures of the LaHood Formation to the Paleoproterozoic crust of the Little Belt arc of the GFTZ. This difference is especially pronounced within lower members of the formation, which has been inferred by Fox (2017) to imply that early sedimentation occurred in a series of compartmentalized sub-basins, later integrated into one large basin.


The Cambrian Flathead Sandstone reflects a shift in provenance following the Great Unconformity, locally representing a period of approximately one billion years for which no sedimentary record is preserved. The Flathead Sandstone in the Bridger Range is best reconstructed with 49% input from the proximal Archean craton (Fig. 5A, Cambrian (C) source 2), 47% input from underlying rocks of the Belt Supergroup (Fig. 5A, C source 4), and 4% input from southern sources in the Cheyenne Belt (Fig. 5A, C source 3), suggesting that detritus was predominantly derived from recycling of proximal rocks (Fig. 6).

Much of the detrital zircon age distribution of the Flathead Sandstone reflects that of the underlying Mesoproterozoic LaHood Formation (Fig. 3), despite the major tectonic events occurring in the intervening ∼1 b.y. separating the deposition of these two units. We interpret this to reflect local recycling of the underlying strata into the Flathead Sandstone. Detrital zircon age distributions reflecting ages of underlying rocks are nearly ubiquitous within basal Sauk strata of the North American Neoproterozoic–Cambrian passive margin (e.g., Warren et al., 2008; Malone et al., 2017; Matthews et al., 2018). Regional provenance patterns and paleogeographic reconstruction highlight this trend (Fig. 6). Cambrian units directly overlying the Great Unconformity generally reflect the detrital zircon age distribution of proximal underlying strata (Fig. 6, samples 43–48, sample Cf). NMF results for samples that overlie strata of the Belt Supergroup (Fig. 6, samples 43–46) are weighted heavily in favor of recycled Belt Basin detritus (Fig. 5A, C source 4).

The vast period of subaerial erosion preceding the Sauk transgression (Sloss, 1963) likely provided an abundance of locally derived sediment. The abruptly rising sea levels of this transgression inundated the region, reworking some of this previously eroded detritus into the Flathead Sandstone. This is supported by a widespread conglomeratic zone at the base of the unit (e.g., McMannis, 1955; Bell, 1970; Malone et al., 2017). This conglomerate has been observed to be transitional to a highly weathered zone up to 30 cm in thickness at the top of the underlying Precambrian basement (Bell and Middleton, 1978). This zone has been interpreted as a regolith developed prior to the sedimentation of the basal Flathead Sandstone (Bell and Middleton, 1978; Beebee, 1997; Beebee and Cox, 1998). Additionally, the basal Flathead Sandstone is rich in coarse feldspar clasts and granitic lithics similar in composition to nearby basement exposures, supporting a local provenance source (e.g., Beebee and Cox, 1998).

Local highs on the pre-Flathead Sandstone topography likely provided a significant amount of detritus to the lower Flathead Sandstone as well. This is supported by the presence of zones in the basal Flathead rich in texturally immature feldspar grains which are found proximal to basement highs of either feldspar-rich crystalline basement (Hubert, 1962) or arkosic strata of the Belt Supergroup (Graham and Suttner, 1974). Regionally, the thickness of the lower Flathead Sandstone is strongly controlled by the topography of the underlying basement, with thicknesses of zero observed over basement highs (Graham and Suttner, 1974; Bell and Middleton, 1978). These highs likely provided a more enduring source of locally derived detritus following the incorporation of existing regolith.

The Neo- to Mesoarchean age distribution of the Flathead Sandstone sample from the Bridger Range contains several populations not present in the underlying LaHood Formation. Prominent peaks in the Flathead Sandstone centered at 2.68 Ga and 2.95 Ga are absent in underlying Belt stratigraphy. Grains of these ages were likely derived from the proximal Wyoming craton (Fig. 5A, C source 2), rather than sources farther south. Approximately 6 km to the south of the sample location is the Perry Line, south of which Belt Supergroup strata are absent and the Flathead Sandstone directly overlies Archean basement of the Wyoming Province (Berg et al., 2000). Additionally, isopach data reveal an ESE-WSW–trending topographic high composed of Archean rocks of the Wyoming Province present in southern Montana during deposition of the Flathead Sandstone (Forhan, 1949). This topographic high, in combination with locally exposed basement, likely supplied this Neo–Mesoarchean detritus to the Flathead Sandstone at the location of the Bridger Range.

Samples from higher in the Flathead Sandstone stratigraphy (Fig. 6, samples 49–53) do not strongly display correlation with underlying basement. This suggests that once available detritus liberated during the development of the Great Unconformity was recycled into the basal Flathead Sandstone and correlative strata and local basement sources were blanketed by sediment, a more complex provenance regime developed with further-traveled sources playing an important role. Cambrian sources 3 and 4 (Fig. 5A) are weighted heavily in these units (Fig. 6). This source may reflect primary detritus shed from uplifted Yavapai/Mazatzal and Cheyenne Belt strata of the Transcontinental Arch to the southeast. Early Cambrian uplift of the Transcontinental Arch has been inferred based on the disappearance of Grenville-age (1250–950 Ma) zircons in the Northern Rockies (Link et al., 2017). These zircons, inferred to be derived from the east coast of Laurentia, are found in the Ediacaran lower Brigham Group of Idaho, but are absent in all Cambrian samples in our compilation, as well as the Lower Cambrian upper Brigham Group in Idaho (Link et al., 2017). This topographic high cut off western basins from eastern sources and likely uplifted Yavapai/Mazatzal and Cheyenne Belt basement rocks to the southeast of our study area. This source area is consistent with dominant E-W directed paleocurrent measurements from Flathead Sandstone exposures in this region (Malone et al., 2017 and references therein). This interpretation is also consistent with increasing compositional maturity of the Flathead Sandstone up-section—with the disappearance of coarse feldspar grain common in basal beds (Bell and Middleton, 1978)—which is indicative of a transition to more distal sources.

Source 4 may also reflect recycled Belt strata. Belt Supergroup rocks, being highly resistant to erosion, likely formed topographic “islands” in much of southwestern Montana during upper Flathead Sandstone sedimentation, as interpreted by Graham and Suttner (1974) on the basis of regions where the Flathead Sandstone is absent. This facilitated the shedding of Belt Supergroup detritus to upper Flathead Sandstone and Upper Cambrian units once paleotopographic lows were filled in during the initial transgression. However, the absence of feldspar grains in most exposures of the upper Flathead Sandstone (e.g., Beebee and Cox, 1998) suggests these topographic highs of Belt strata were likely only a locally significant source of sediment for upper Flathead Sandstone strata.

In contrast to the more distal sediment dispersal pathways supplying detritus to southwestern Montana in Late Cambrian time, sub-basins in Idaho continued to incorporate proximal detritus into the Upper Cambrian. These units, including the Furongian Worm Creek Formation of southeastern Idaho Formation (Link et al., 2017), display a strong affinity for 500–490 Ma zircons. These zircons have been interpreted by (Link et al., 2017) to have been sourced from the nearby Big Creek-Beaverhead magmatic belt exposed in the Lemhi Arch of north-central Idaho, a northwest-trending region of thinned or missing Neoproterozoic to Cambrian strata (Sloss, 1954). There, the Beaverhead plutons (500–485 Ma) intrude strata of the Mesoproterozoic Swauger and Apple Creek formations (Fig. 5A, C source 1). Samples from the Upper Cambrian Pilgrim Formation in the Highland Mountains of southwestern Montana, Du Noir Member of the Gallatin Limestone in Wind River Canyon, Wyoming, and upper Wilbert Formation of the central Idaho thrust belt contain these 500–490 Ma zircons as well (Link et al., 2017), indicating this source may have been routed east of the Lemhi Arch into southwestern Montana in Late Cambrian time. Zircons of the Beaverhead pluton are younger than the depositional age of the Flathead Sandstone (Lund et al., 2010; Link et al., 2017) and are consequently absent in the detrital zircon population of Flathead Sandstone samples. However, the Mesoproterozoic rocks exposed in the Lemhi Arch are in contact with Ordovician strata across an angular unconformity in much of the region, so it is possible that eroding Mesoproterozoic rocks of the uplifted Lemhi Arch supplied sediments along similar dispersal routes prior to the intrusion of these plutons at 490–500 Ma.

Our results are broadly consistent with previous interpretations of Cambrian paleogeography, wherein sediment deposition was on a west-facing shelf of a passive margin, with deposition grading eastward from subtidal through fluvial environments—with the fluvial systems tapping the highs of the Transcontinental Arch (Sloss, 1963; Gehrels et al., 2011; Matthews et al., 2018). Additionally, the presence of 500–490 Ma zircons related to the Big Creek-Beaverhead plutons in Upper Cambrian rocks of southeastern and central Idaho, western Montana, and northern Wyoming suggest the Lemhi Arch was a regionally important source of sediment in the Northern Rockies by at least Late Cambrian time.


Units deposited in Late Devonian to Pennsylvanian time reflect a major change in provenance from local recycling during Cambrian time to more distal sources associated with collision on nearly all margins of Laurentia. NMF analysis predicts five sources (Fig. 5B) supplying sediment to the Devonian–Pennsylvanian system: (1) recycled Ordovician strata, widespread but dominantly recycled off the Canadian shield; (2) Proterozoic crust from terranes along the western Laurentian margin and basement uplifts of the Yavapai/Mazatzal Province; (3) the Appalachian clastic wedge distally to the east; (4) the clastic wedge of Arctic orogens to the north; and (5) the Antler clastic wedge/allochthons to the south. The Upper Devonian Sappington Formation of the Bridger Range (Fig. 7A) is best reconstructed with ∼68% input from distal eastern sources (Fig. 5B, D–IP source 3), ∼17% input from Proterozoic terranes (Fig. 5B, D–IP source 2), ∼8% input from sources with headwaters in the arctic orogens (Fig. 5B, D–IP source 4), and ∼7% input from recycled sources of the Antler orogen (Fig. 5B, D–IP source 5).

This interpretation underscores the increasing importance of distally derived detritus in the sediment budget of Paleozoic depositional systems of the Northern Rockies. However, it is unlikely that this detritus was derived directly from these long-traveled sources. The Late Devonian Sappington Basin was surrounded by the topographic high of the Antler forebulge to the west, central Montana uplift to the north, and the Beartooth shelf to the south and east (Fig. 7A; Phelps et al., 2018). These highs at least partially cut the Sappington Basin off from the Antler foredeep to the west, Alberta Basin to the north, and Bighorn Basin to the south. However, the provenance signature of the Sappington Formation is very similar to that of the Bakken Formation of the Williston Basin to the northeast (Fig. 7A). These basins were likely joined by a seaway occupying the central Montana trough during sea level highstands (Phelps et al., 2015). The Beartooth shelf has been suggested as a sediment source for the Sappington Basin (Phelps et al., 2018). Based on the similarity of provenance between these two basins, we propose that long-traveled sediment routed to the Williston Basin may have been deposited on the eastern high of the Beartooth shelf during sea level highstands. During periods of low sea level, this sediment was eroded and passed into the Sappington Basin (Fig. 7A). This sediment was then distributed across the Sappington Basin likely by counterclockwise currents (Phelps et al., 2018).

Regional provenance patterns show a dominance of eastern-derived detritus in the Williston and Sappington basins. The presence of the Transcontinental Arch (Carlson, 1999) as well as several major intracratonic basins to the east and southeast of the Northern Rockies, such as the Anadarko, Michigan, and Illinois basins (Miall, 2019), posed a barrier to routing this sediment directly from eastern orogens. However, a northern route directing sediment from the northern Appalachian and Franklinian-Caledonian clastic wedge in northern Canada and the Arctic may have been viable by Late Devonian time, as suggested to explain provenance of Upper Devonian strata of southern British Columbia by Gehrels and Pecha (2014). The dominance of D–IP source 3 in the Sappington and Williston basins is consistent with a sediment routing system with headwaters in the northeastern Canada. Large fluvial systems likely brought detritus from the Caledonian-Appalachian clastic wedge along routes north of the Transcontinental Arch to the Williston Basin. This interpretation is consistent with the absence of eastern-derived detrital zircons in Devonian strata of the Grand Canyon, which are interpreted to have derived detritus predominantly from basement rocks of the southwestern United States (Gehrels et al., 2011). The southwestern United States was likely cutoff from eastern sources by the Transcontinental Arch to the east and Williston and Sappington basins to the north.

Strata of the Devonian Antler foredeep (Fig. 7A, samples 39–40) are dominated by a source with a strong Proterozoic peak at ca. 1650 Ma and subsidiary ca. 440 Ma peak (Fig. 5B, D–IP source 2; Fig. 7A). This fits with the hypothesis of Beranek et al. (2016) that sediment was routed from outboard Paleozoic arcs built on Proterozoic crust. Beranek et al. (2016) infer a western source for the Milligen and Jefferson formations of central Idaho, which have a strong affinity for our D–IP source 2, on the basis of detrital zircon U-Pb and Lu-Hf results. While no such arcs are found accreted at the latitude of Montana, strata of the eastern Klamath terrane of California (Fig. 1A) have a very similar detrital zircon signature (Grove et al., 2008). This may suggest that a transcurrent tectonic setting broadly consistent with the hypothesis of Eisbacher (1983) was present during this time in the western United States, connecting the Ellesmerian orogen of the Arctic to the Antler orogen of Nevada. Sinistral southward displacement of exotic crustal fragments, such as the Yreka subterrane of the Klamath mountains of California and the Okanagan subterrane of the Kootenay terrane in southeastern British Columbia (Fig. 1A), at the rates proposed by Colpron and Nelson (2011) places such terranes near the latitude of the Northern Rockies during the Late Devonian. Transpressional regions along this margin may have created uplifted regions which shed detritus onto the Laurentian shelf. This is supported by Middle Devonian deformation in the Purcell Mountains of British Columbia, near where the exotic Okanagan subterrane was sutured to North American (Colpron and Nelson, 2009) and by the presence of mid-Paleozoic conglomeratic strata of the Milford group in southern British Columbia, which contain 418 Ma and 431 Ma granitoid boulders with an inferred western source (Thompson et al., 2006).

In contrast, Lemieux et al. (2007) and Gehrels and Pecha (2014) suggest predominantly Laurentian sources for contemporaneous strata of the Chase and Sassenach formations of southern British Columbia, with possible input from the Alexander Terrane—hypothesized to have been located along the paleo-Arctic margin in Devonian time (Soja, 1994; Gehrels et al., 1996; Colpron et al., 2007; Grove et al., 2008). We largely favor western sources for samples with strong affinity for D–IP source 2 (Fig. 7A, samples 39–40) and Laurentian sources for samples farther east, with a higher proportion of D–IP source 2 and 3 (Fig. 7A, samples 37, Ds). Mixing of these two sources may have occurred on the western extents of the continental margin, such as in the Sassenach Formation, (Fig. 7A, sample 38) where a hybrid of these two source affinities is present.


During Mississippian time, the differences between samples within central Montana and the Antler foredeep were magnified (Fig. 7B). The Kibbey Sandstone sample from the Bridger Range is weighted in favor of eastern and arctic sources (Fig. 5B, D–IP sources 3 and 4) at 64% and 27%, respectively, with the remaining 7% of detritus favoring Proterozoic terranes (Fig. 5B, D–IP source 2). Eastern derived detritus may have been sourced directly from eastern orogens via oceanic currents routed through the inundated Colorado sag—a region between the persistent Paleozoic highs of the Uncompahgre and Front Range tectonic cores of the Transcontinental Arch (Carlson, 1999)—into Montana during early Mississippian time given the similarity of D–IP source 3 to the Appalachian clastic wedge (Fig. 5B; Chapman and Laskowski, 2019). Erosion of Late Devonian–early Mississippian proximal shelf and Williston Basin strata likely also supplied recycled sediment to the Kibbey Sandstone. Northern routes draining the Franklinian-Caledonian clastic wedge were likely still active during Mississippian time, as suggested by the abundance of D–IP source 4 (Fig. 5B). This provenance is consistent between samples within Montana (Fig. 7B, samples 35, Mk) but contrasts strongly with those deposited in the Antler foredeep, where almost no eastern-derived detritus is present.

Samples deposited in western basins favor sediment recycled from Ordovician strata (Fig. 5B, D–IP source 1; Gehrels and Pecha, 2014). Beranek et al. (2016) proposed these strata may have been uplifted along a transpressional restraining bend along the Saint Mary-Moyie transform zone of northern Idaho/northeastern Washington and a releasing bend forming the Copper Basin of central Idaho. Within the Thompson and Davis assemblage of British Columbia, the basal member requires a Proterozoic source like that derived from the Roberts Mountain allochthon of Nevada (D–IP source 5) as well as dominant input from recycled Ordovician strata (D–IP source 1; Cole et al., 2015). Uplift along this transpressional zone may have exhumed Ordovician strata atop accreted Proterozoic crust by early Mississippian time, possibly in the uplifted block termed the Okanagan high (Fig. 7B; Thompson et al., 2006; Colpron and Nelson, 2009; Beranek et al., 2016). Additional detritus from the Roberts Mountain allochthon may have been routed northward through the Antler foredeep to central Idaho—an interpretation broadly consistent with paleocurrent trends (e.g., Di Bona, 1982).


By Pennsylvanian time, provenance within the Northern Rockies system was well mixed, reflecting recycling of local sedimentary rocks combining with primary input from eastern orogens and basement rocks of the northern Colorado and Utah region. At the location of the Bridger Range, NMF analysis of samples from the Amsden Formation and overlying Quadrant Formation yield similar source weights (Fig. 7C). Every synthetic source population is represented with dominance of those sourced from peri-Laurentian orogens (Fig. 5B, D–IP sources 3 and 4), accounting for ∼45% of the age distribution of the Amsden Formation and ∼74% of that of the Quadrant Formation. Appalachian-derived detritus (Fig. 5B, D–IP source 3) was ubiquitous within the Northern Rockies system in Pennsylvanian time (Fig. 7C). The weight of Appalachian sources generally decreases in abundance away from the Wyoming sag—a region between the persistent Paleozoic highs of the Cambridge and Front Range tectonic cores of the Transcontinental Arch—which likely provided a marine corridor between the Northern Rockies and Appalachian systems during early Pennsylvanian time (Carlson, 1999). This suggests that detritus from this source may have been transported by northwest-directed ocean currents driven by equatorial trade winds, which is generally reflected in dominantly E-W paleocurrent trends (Leckie, 1964). Alternatively, a network of westward prograding deltas, large river systems with outflows in the Great Plains, and SW-migrating dune fields may have been the dominant drivers of this sediment dispersal, as suggested by Chapman and Laskowski (2019) and references therein.

Sediment routed from Proterozoic terranes (Fig. 5B, D–IP source 2) is also widespread in Pennsylvanian strata. It is most abundant in southern samples, decreasing in abundance to the north as depocenters become more distal to basement sources of Yavapai/Mazatzal affinity exposed in southeastern Wyoming and the Four Corners region (Fig. 7C; Mallory, 1972). These blocks were exhumed by the Ancestral Rockies orogeny (Kluth and Coney, 1981; Maughan, 1990), which resulted in basement uplift of Proterozoic crust in the four-corners region through Wyoming due to collision with South America/Africa along the southern margin of Laurentia.

Broadly homogenous provenance within the Northern Rockies system was probably driven by recycling of previous Paleozoic sediment—this effectively averaged the age distributions resulting from the prior ∼200 m.y. of Paleozoic provenance evolution. Recycling was facilitated by eustatic sea level lowstands in the early Pennsylvanian at levels unprecedented since Precambrian time (Haq and Schutter, 2008). Regions on the Beartooth shelf and in north-central Montana were subaerially exposed by low sea level in the early Pennsylvanian and likely shed detritus recycled from underlying Paleozoic sedimentary rocks such as the Big Snowy Group, from which more than 300 m was likely eroded in south-central Montana (Maughan, 1989). Recycled Ordovician strata uplifted in northwestern Montana, Idaho, and Alberta has been proposed as a significant source of sediment for Pennsylvanian erg deposits of the Quadrant Formation and correlative strata (Maughan, 1989). The abundance of D–IP source 1 in the Quadrant Formation of western Montana (Chapman and Laskowski, 2019) supports this interpretation. However, this source decreases in abundance to the east in the thicker Quadrant deposits of central-Montana. These northern sources may have incorporated detritus shed from the Franklinian clastic wedge (D–IP source 4) by middle Pennsylvanian time, as suggested by Leary et al., 2020. This is reflected in the relative abundance of this source in our Quadrant Formation sample from the Bridger Range and supported by a shift to these ages in middle Pennsylvanian strata farther south in the Rocky Mountains (Leary et al., 2020).

Late Jurassic

In Mesozoic time, the western margin of the North American continent transitioned to a predominantly two-plate system, with Farallon plate subduction beneath North America (DeCelles, 2004; Evenchick et al., 2007). Regional subsidence in the foreland began in Bajocian time, with marine sedimentation of the Ellis Group in the distal retroarc region (DeCelles and Currie, 1996; Fuentes et al., 2009, 2012). The onset of orogenesis is recorded in the detrital zircon record of the Northern Rockies, with provenance shifting from eastern to dominantly western potentially by Bajocian time (Stronach, 1984; Fuentes et al., 2011) and certainly by Oxfordian time (Syzdek et al., 2019, this study). Broadly, Mesozoic strata of the Northern Rockies record a basin-wide unroofing sequence recycling Paleozoic passive margin cover into the foreland basin of the developing orogen and the mobilization of detritus derived from newly intruded/accreted arcs and terranes (Fig. 8A).

The Middle Jurassic strata of the Ellis Group represent the earliest strata that have been hypothesized to have been deposited in a basin driven by flexural subsidence in the foreland (DeCelles and Currie, 1996; Fuentes et al., 2009, 2011). A western source was hypothesized for the Middle Jurassic Sawtooth Formation of northern Montana based on a single Mesozoic zircon (Fuentes et al., 2009). However, the onset of arc-derived detritus reaching the foreland basin system was evidently delayed until at least the Callovian in the Powder River Basin of Wyoming (Fig. 1A), where sediment from the lower Callovian sandstone of the Jurassic Sundance Formation is devoid of Mesozoic grains while the upper sandstone contains them in abundance (Syzdek et al., 2019). Additionally, the MDA of the Swift Sandstone of the Ellis Group in the Bridger Range is 154.8 ± 0.58 Ma (Table 2), which supports a younger depositional age compared to the Ellis Group in northern Montana (Fuentes et al., 2011). This may suggest diachronous shifts from cratonic sources to Cordilleran sources within the foreland basin. Alternatively, the Sawtooth Formation may have been largely tapping eastern cratonic sources with sparse input of ashfall arc zircons. More detrital zircon data on the Sawtooth Formation would likely clarify this point. In either case, the Northern Rockies system was unequivocally dominated by western sources associated with the developing fold-thrust belt and cordilleran arcs/terranes by late Oxfordian time.

At the location of the Bridger Range, NMF analysis weights detritus shed from western arc sources and routed through the Sevier fold-thrust belt (Fig. 5C, 68.3% Mesozoic (Mz) source 5) as the dominant source in the Swift Sandstone. Additionally, Triassic grains derived from the Wallowa/Olds Ferry/Intermontane terranes are present in abundance (Fig. 5C, 25.7% Mz source 3). Given the lack of underlying strata containing these young ages, it is likely these represent early-cycle grains routed along generally east-directed fluvial systems draining western arcs and accreted terranes to the marginal marine environment of Swift Sandstone deposition (Fig. 8A). The abundance of arc detritus in the Bridger Range sample relative to other Oxfordian strata in the Northern Rockies may indicate this region received input directly from fluvial systems with their origins within the magmatic arc.

Regionally, samples from northern Wyoming (May et al., 2013; Syzdek et al., 2019) show an increased abundance of grains derived from recycling of the U.S. passive margin while samples in northern Montana display a well-mixed source signature (Fig. 8A). This likely reflects northwest directed basin-axial sediment transport within the foreland basin system, likely via longshore drift within a marginal marine environment. This interpretation is supported by dominantly S-N directed paleocurrent trends (e.g., Ballard, 1961) and previous provenance studies (Hamblin and Walker, 1979; Miles et al., 2012; Raines et al., 2013; Quinn et al., 2018). East- to southeast-directed routing of sediment directly from the fold-thrust belt to the foreland basin in Montana was also likely of importance, shown by the addition of Canadian passive margin detritus in northern Montanan samples (Quinn et al., 2018).

By Kimmeridgian–Tithonian time, recycled U.S. passive margin detritus (Fig. 5C, Mesozoic (Mz) source 1) dominated within most southern samples (Fig. 8B). This trend is consistent at the location of the Bridger Range where the Morrison Formation received detritus from the recycled passive margin (Fig. 5C, ∼70% Mz source 1), as well as western arcs and terranes (Fig. 5C, ∼30% Mz source 1, 2, 4). This reflects the continued eastward propagation and development of Late Jurassic hinterland thrust systems in the Northern Rockies as a topographic barrier to arc detritus (Fuentes et al., 2012). This is evident in the lack of zircons derived from the Intermontane Belt within samples in the Canadian foreland basin (Fig. 8B). Sediment from Cordilleran arcs/terranes is concentrated in samples 13 and 16, which were likely along major fluvial systems traversing the fold-thrust belt. Basin axial transport was present during deposition of the Morrison Formation, but may have at times had a southward-directed component, as indicated by transport of detritus from the Canadian passive margin (Fig. 5C, Mz source 2) to samples in northern Wyoming.

These provenance interpretations are broadly consistent with the deposition of the Swift Sandstone and Morrison Formation in a distal back-bulge depozone of the developing foreland basin, as suggested by DeCelles (2004) and Fuentes et al. (2011) based on slow subsidence rates and tabular geometry of the Late Jurassic deposits. This region may have benefitted from increased subsidence due to dynamic coupling with the upper mantle above mantle lithosphere subducting at a low angle (Currie, 1998) and/or mantle processes related to the subsidence of the Williston Basin (Fuentes et al., 2011). This back-bulge depozone was likely complicated by the activation of basement structures such as the Sweetgrass Arch, Kevin-Sunburst Dome, and South Arch (Fig. 1A), explaining intraformational unconformities of Ellis Group strata in Montana (Cobban, 1945; Suttner, 1969; Peterson, 1981; Suttner et al., 1981; Parcell and Williams, 2005; Fuentes et al., 2011).


By Aptian time, the fold-thrust belt had propagated into Idaho, recycling uplifted passive margin strata into the foreland basin (DeCelles, 2004 and references therein). The abundance of chert in the basal Kootenai Formation and correlative units of the Northern Rockies reflects recycling of Pennsylvanian–Permian rocks, such as the chert-rich Phosphoria Formation, to the west (Rapson, 1965; Suttner et al., 1981; DeCelles, 1986; Schwartz and DeCelles, 1988; Fuentes et al., 2011). The Kootenai Formation and correlative Aptian strata in Alberta and south through Utah represent the first unequivocal record of foredeep sedimentation within the Northern Rockies system (Suttner, 1969; DeCelles, 1986; Schwartz and DeCelles, 1988). These units thicken abruptly to the west, from <70 m at the longitude of the Bridger Range to 400 m in western Montana, reaching maximum thicknesses of greater than a kilometer in central Utah (DeCelles, 2004 and references therein). This wedge shape, as well as flexural modeling (e.g., Jordan, 1981) and subsidence curves (e.g., Fuentes et al., 2011) provide abundant evidence for deposition in a foredeep setting.

In the Bridger Range, the Kootenai Formation incorporates detritus from all five Mesozoic synthetic source populations in near equal measure (Figs. 5 and 8C). This likely reflects recycling of Late Jurassic foredeep deposits equivalent to the Swift Sandstone and Morrison formations, which are not preserved in the rock record. These strata were situated to the west and subsequently thrust over structurally lower elements of the frontal thrust belt during Early Cretaceous propagation of the fold-thrust belt and largely eroded, creating a missing or “phantom foredeep” (Royse, 1993). The existence of this “phantom foredeep” is supported by 200–300 °C burial temperatures—suggesting ∼6 km of burial—indicated by conodont alteration indices from Mississippian strata underlying the proposed foredeep (Sandberg and Gutschick, 1984). Additionally, relicts of this foredeep, such as the Late Jurassic Fernie Formation, are preserved in the hanging walls of major thrust systems of the Canadian Rockies (Stronach, 1984). This recycled detritus likely overwhelmed first-cycle arc grains within the southwestern Montana system, which is supported by the lack of syn-depositional aged zircons in the basal Kootenai Formation in the Bridger Range (Table 2). The MDA for our Kootenai Formation sample indicates a Kimmeridgian depositional age. However, given robust palynologic data suggesting Hauterivian onset of Kootenai Formation deposition (Fuentes et al., 2011), our MDA is likely skewed by the relative abundance of recycled grains, largely from the underlying Morrison Formation and its foredeep equivalent, versus first-cycle grains. In contrast, samples from the Kootenai Formation northern Montana contain depositional-age arc zircons (Fig. 8, sample 5), suggesting routing of arc sources to the foreland basin, potentially in a structural gap between the Hawley Creek and Moyie thrust sheets (Fig. 8C; Fuentes et al., 2011).

Samples from coeval strata in northern Utah (sample 10, Fig. 8C), show a heavy contribution by an ca. 120 Ma undiluted arc source (Fig. 5C, Mz source 4), with contribution by recycled Mesozoic eolianites (Laskowski et al., 2013). This suggests a relatively direct dispersal pathway between the magmatic arc, likely the Sierra Nevada arc, and foreland basin of southern Wyoming (Fig. 8C, sample 10) in the Early Cretaceous. In contrast, Canadian Albian samples (Fig. 8C, samples 1–4) have a much smaller component of arc-derived grains and are dominated by detritus recycled from the Canadian passive margin by the eastward propagation of the thrust sheets within the Intermontane Belt (Fig. 8C).

Tectonic partitioning of the foreland basin during deposition of the Kootenai Formation has been argued on the basis of evidence for regionally complex paleocurrent directions and the abrupt deflections of these directions, intra- and interformational unconformities, intraformational recycling, and dramatic base-level alterations within this formation (DeCelles, 1986). Provenance data support this interpretation—Lower Cretaceous strata in Montana exhibit largely heterogenous provenance in small distances east-west (Fig. 8). This may indicate compartmentalized, largely basin-axial fluvial to lacustrine systems each tapping different source regions to the west and experiencing little mixing perpendicular to strike. This is illustrated by the relatively similar provenance between samples at a similar longitude (i.e., Fig. 8C, samples 6 and 7, sample Kk), but stark differences with slight longitudinal variation (i.e., Fig. 8C, sample 5 versus sample 6 and 7, sample Kk versus sample 8), potentially indicating a tectonically partitioned basin system. This partitioning was likely most intense in southwestern Montana, where the basal conglomeratic member of the Kootenai Formation is reworked over the present Tobacco Root, Madison, Gravelly, and Beartooth Archean blocks (Fig. 1A; DeCelles, 1986). This suggests that these blocks were being uplifted during Kootenai Formation deposition, creating intraforeland drainage divides and producing differences in provenance signatures across these features. This intraforeland partitioning and subsequent diversification of provenance signatures within the foreland basin has been documented in Paleogene strata of the Northern Rockies (Schwartz et al., 2019)—we hypothesize this effect may have been occurring, albeit on a much smaller scale, as early as the Early Cretaceous. Recent radiometric thermochronological data from these ranges, as well as the Highland Mountains and Ruby Range (Fig. 1A), places the inception of uplift of these ranges earlier than previously thought (Carrapa et al., 2019). Exhumation beginning as early as 120 Ma in the Beartooth and Ruby ranges—with exhumation in all studied ranges under way by at least 80 Ma—was inferred by Carrapa et al. (2019). This contrasts with previous studies that suggest exhumation in southwestern Montana began at ca. 80 Ma or later (Peyton and Carrapa, 2013 and references within). Exhumation of Archean basement in the region during the Late Cretaceous is supported by sandstone petrography that suggests Blacktail-Snowcrest Arch in southwest Montana was actively forming and sourcing the synorogenic Beaverhead Conglomerate during this time (Haley, 1986; Haley and Perry 1991; Garber et al., 2020).

An extensive record of the tectonic evolution of the Northern Rocky Mountains and the coupled evolution of its sedimentary basins spanning from the Mesoproterozoic through the Cretaceous is preserved in the stratigraphy of the Bridger Range, southwest Montana. The detrital zircon record within this range and correlative strata of the Northern Rockies elucidates a complex story of the evolution of sediment-dispersal pathways and the basins they supply. The subtleties of this coupled tectonic and provenance evolution are deconvolved using inverse statistical modeling, NMF. This work demonstrates that NMF is capable of providing a statistically robust basis for interpreting complex sediment-dispersal patterns that involve first-order and recycled sources and provides a work-flow by which to do so. The evolution of the provenance and paleogeography of the Northern Rockies system is as follows:

  • (1) The Mesoproterozoic Helena Embayment was likely an aulacogen, which began as a series of compartmentalized basins deriving sediment from local sources such as the proximal Archean craton and the arcs and reworked crust of the Great Falls Tectonic Zone exposed along the down-dropping faults of the Perry and Garnet Line. These sub-basins became integrated through time. As subsidence slowed, the basin derived an increasing proportion of sediment from sources in the southern Wyoming Province.

  • (2) Following development of the Great Unconformity, Middle Cambrian sedimentation in the Northern Rockies predominantly involved recycling of local Proterozoic–Archean sources in a transgressive marine shelf setting. As topographic lows were infilled, sediment sourced from eroding strata of the Proterozoic Belt Basin and the Transcontinental and Lemhi arches became dominant.

  • (3) Provenance shifted from local to distal Laurentian sources in the middle Paleozoic, with northern sediment-dispersal routes tapping the northern Appalachian and Franklinian clastic wedges supplying sediment to the eastern Northern Rockies system by Devonian time. By the early Mississippian, there was likely the addition of a component of southeastern sources, which breached the Transcontinental Arch. Mid-Paleozoic tectonism had a strong influence on western depozones—a complex plate setting along the western margin of Laurentia facilitated sediment routing from terranes and uplifted Ordovician strata to the west. By the Pennsylvanian, detritus from the Ancestral Rockies and Appalachian orogens became more important. Regional recycling contributed to a well-integrated provenance signal across the Northern Rockies during this time.

  • (4) The establishment of a coherent North America Cordilleran system with a retroarc fold-thrust belt by Late Jurassic time led to the development of a retroarc foreland basin. This is recorded by a shift to the dominance of western provenance by Oxfordian time, evident in the arrival of syn-depositional arc-derived grains. Early foreland basin deposition is recorded in back-bulge deposits of the Swift Sandstone and Morrison formations and correlative units. Basin-axial transport supplied detritus recycled out of the fold-thrust belt mixed with first-order arc grains to this depozone. By the Cretaceous, detritus recycled from Paleozoic and foreland basin strata dominated foreland basin provenance as the fold-thrust belt propagated into Idaho. Within southwestern Montana, this system was likely partitioned by early intraforeland uplifts centered above modern basement-cored uplifts.

Fruitful discussions with Andrew Laskowski, Caden Howlett, and Stuart Parker were instrumental throughout this process. Constructive reviews by Tim Lawton and Marc Hendrix greatly improved the quality of this manuscript. We would like to thank Editor Brad Singer and Associate Editor Rajat Mazumder for the handling of this manuscript. We would also like to thank Colorado Plateau Geosystems Inc. for academic usage of their paleogeographic maps. We would like to gratefully acknowledge the assistance of GeoSep Services for help with mineral separation and the Arizona Laserchron Center for assistance with detrital zircon analyses. Funding was provided by Montana State University, the Montana State University Undergraduate Scholars Program, and the Montana Academy of Sciences.

1Supplemental Material. Analytical techniques, data reduction methods, ICP-MS zircon U-Pb datasets, comparison of age-probability curves generated using PDPs and various KDE bandwidths, sample information, and NMF details and results. Please visit to access the supplemental material, and contact with any questions.
Science Editor: Brad S. Singer
Associate Editor: Rajat Mazumder
Gold Open Access: This paper is published under the terms of the CC-BY license.