Foreland basin response to middle Cretaceous thrust belt evolution, southwestern Montana, USA

In Cordilleran orogenic systems, subduction of an oceanic plate beneath a continental plate creates a volcanic arc and fold-thrust belt on the over riding plate. This results in topography that serves as a sediment source and creates flexural deflection of the lithosphere, allowing for deposition in the adjacent foreland basin. Whereas the arc and foldthrust belt are often only partially preserved due to subsequent deformation and erosion, foreland basin strata can record the entire evolution of the hinterland system and are essential to interpreting the tectonic histories of ancient margins. The North American retroforeland basin is an archetypal foreland basin that records Cordilleran orogenesis, generally considered to have initiated in the Jurassic, with subduction-driven retroarc contraction during the Sevier orogeny (e.g., DeCelles, 2004). Following Early Cretaceous sinistral oblique plate convergence, substantial dextral oblique subduction (e.g., Seton et al., 2012) and terrane translation initiated in middle Cretaceous time (e.g., Enkin et al., 2006; Wyld et al., 2006). High rates of middle Cretaceous plate convergence were concomitant with a marked increase in retroarc fold-thrust belt deformation at ca. 105 Ma (e.g., Pană and van der Pluijm, 2015; Yonkee et al., 2019). Though many basementinvolved uplifts of the Laramide province are thought to have initiated after ca. 70 Ma (e.g., Lawton, 2019; Peyton et al., 2012), previous work in southwestern Montana has demonstrated that certain Laramidestyle basement uplifts were exhumed down to upper Paleozoic strata by Coniacian–Campanian time (ca. 88–81 Ma; Nichols et al., 1985; Garber et al., 2020). Preliminary lowtemperature thermochronometric data suggest that these structures may have begun exhuming as early as ca. 100 Ma (Carrapa et al., 2019). Southwestern Montana contains a nearly continuous and well-preserved record of sedimentation that allows for testing models of middle Cretaceous tectonism. Middle Cretaceous strata (ca. 106–90 Ma), however, have not been investigated in the context of regional provenance and are critical for constraining the timing and character of Sevier fold-thrust belt progression and inception of basementinvolved thrusting.

A longstanding problem in foreland provenance studies is that sedimentary rocks deposited during orogenesis are often composed of detritus from older rocks; therefore, detrital zircon U-Pb ages extracted from foreland strata necessarily predate the age of the recycled source rock. The robustness of zircon leads to excellent preservation across many sedimentary cycles, which results in recurring, continent-scale distributions of dominant age populations (Schwartz et al., 2019). Graphical comparisons of individual detrital zircon populations among basin samples and potential sedimentary source units may reveal some provenance patterns, but with small-n analyses (n = 100), spectra do not accurately characterize diverse populations (Saylor and Sundell, 2016). In contrast, large-n (n = 600) detrital analyses—coupled with robust detrital zircon characterization of potential eroded source strata—and statistical models linking both sample and source signatures, are vastly more useful to accurately infer provenance and relative source contributions (Sundell and Saylor, 2017).

Our work presented herein focuses on the provenance of the Albian–Cenomanian (middle Cretaceous) upper Blackleaf Formation in southwestern Montana and the subsequent tectonic implications. We aim to identify and accurately quantify relative proportions of various sources of recycled sediment to develop a detailed provenance model for the upper Blackleaf Formation using large-n data sets and inverse Monte Carlo mixture modeling (Sundell and Saylor, 2017). Basin-wide sampling of two successive lithostratigraphic units allows us to discriminate shifts in provenance and sediment transport both spatially and temporally. These shifts can be linked to hinterland tectonics and thrust-front progradation and offer new insights into the overall middle Cretaceous tectonic development of western North America.

Oblique subduction of the Farallon plate beneath the North America plate drove terrane accretion and translation, arc volcanism, and development of the Sevier fold-thrust belt from the late Jurassic through the Eocene (e.g., Dewey and Bird, 1970; Coney and Evenchick, 1994; DeCelles, 2004; Dickinson, 2004). At the latitude of Idaho and Montana (~44°–46.5° N), the Salmon River Suture Zone (SRSZ) marks the boundary between accreted terranes to the west and the Laurentian craton to the east (Fig. 1; Manduca et al., 1993; Tikoff et al., 2001). Middle Cretaceous arc magmatism, beginning at ca. 120 Ma and continuing to ca. 90 Ma (Gaschnig et al., 2017), was concentrated near the SRSZ and is today preserved as plutonic remnants that crop out in a >200-km-long and ~5-km-wide band, although it is estimated to have been ~85–100 km wide prior to Late Cretaceous shearing (Giorgis et al., 2005). Transpression along the Western Idaho Shear Zone (WISZ) between ca. 104–90 Ma deformed both the SRSZ and the middle Cretaceous plutons (Manduca et al., 1993; McClelland et al., 2000; Giorgis et al., 2005; Giorgis et al., 2008; Braudy et al., 2017; Gaschnig et al., 2017; Giorgis et al., 2017). Subsequent evolution of the magmatic arc resulted in widespread emplacement of the Idaho batholith across central Idaho through the Eocene (Gaschnig et al., 2010; Gaschnig et al., 2017).

Inboard of the arc region, horizontal crustal shortening resulted in eastward propagation of the Sevier fold-thrust belt from central Idaho toward southwestern Montana (Skipp, 1988; Tysdal, 2002). The deformed stratigraphic package in central Idaho consists of the Mesoproterozoic Belt Supergroup (e.g., Link et al., 2016), northeastward-thinning Neoproterozoic–Ordovician rift and passive margin strata (e.g., Brennan et al., 2020), eastwardthinning Devonian–Mississippian mixed siliciclastic and carbonate strata, a Mississippian–Pennsylvanian westward-thickening carbonate bank with minor siliciclastic intervals (e.g., Link et al., 1996; Skipp et al., 1979b; Beranek et al., 2016), poorly preserved Permian sandstone and carbonate (e.g., Mahoney et al., 1991; Rankey, 1997; Link et al., 2014), and Triassic carbonate and mudstone only present today in southeastern Idaho and near the Idaho-Montana border (Peterson, 1988).

The Idaho-Montana fold-thrust belt consists of two overlapping components (Kulik and Schmidt, 1988; Parker and Pearson, 2021): (1) an upper, Sevier-type thrust system that links the more inboard, western part of the thrust belt in central Idaho with similar thrusts in east-central Idaho and southwestern Montana (Fig. 1); and (2) a deeper thrust system that constitutes the Laramide belt of southwestern Montana but also projects beneath the shallower thrust system in east-central Idaho. Relative timing constraints indicate that although the upper and more inboard thrust system is generally older, early thrusting on the deeper system was coeval with the overlying fold-thrust belt and eventually truncated it (Kulik and Schmidt, 1988; Tysdal, 1988; McDowell, 1997; Parker and Pearson, 2021).

The Sevier thrust system that constitutes the fold-thrust belt in central Idaho accommodated shortening of a relatively thick succession of Neoproterozoic through upper Paleozoic strata (Dover, 1981; Rodgers et al., 1995; Montoya, 2019; Porter, 2021). The western part of this thrust system is unique in that it exhumed Neoproterozoic through Cambrian sedimentary rocks that are absent to the northeast above the Lemhi Arch in eastcentral Idaho (Fig. 1), where Ordovician strata sit unconformably upon Mesoproterozoic quartzite (Brennan et al., 2020). Several significant contractional structures within the western portion of the fold-thrust belt were crosscut by 97–91 Ma plutons, which gives a younger bound on fold-thrust belt timing (Montoya, 2019; Brennan et al., 2020; Porter, 2021). These results suggest cessation of deformation within this western portion of the upper thrust system by middle Cretaceous time. Based upon sandstone petrography and detrital zircon U-Pb geochronology for foreland basin strata, Rosenblume et al. (2021b) suggested unroofing to lowermiddle Paleozoic levels within this portion of the fold-thrust belt by Aptian–early Albian time (ca. 125–110 Ma). To the northeast, along the western margin of the Lemhi Arch, shallowly detached thrusting was described by Anastasio et al. (1997) and Hedlund et al. (1994) in the Lost River Range where Mississippian and younger rocks are duplexed and folded. Farther northeast in east-central Idaho and into southwestern Montana, this upper thrust system—which here consists primarily of the Thompson Gulch, Medicine Lodge, and McKenzie thrusts—is confined to Ordovician and younger strata that unconformably overlie Mesoproterozoic quartzite, with few constraints on absolute timing of thrust activity (Parker and Pearson, 2021).

Near the Idaho-Montana border, the deeper thrust system deformed Mesoproterozoic quartzite of the Belt Supergroup, Cambro-Ordovician plutons (Lucchitta, 1966; Parker and Pearson, 2021), and Archean and Paleoproterozoic crystalline basement rock (Skipp, 1988). Some major thrusts here likely reactivated early Paleozoic normal faults and involve large-wavelength folds that suggest relatively deep detachment levels (Parker and Pearson, 2021). The major thrusts near the Idaho-Montana border consist of the Hawley Creek, Baby Joe Gulch, and Poison Creek thrusts, which may be along-strike equivalents (Skipp, 1988; see Lund, 2018, for an alternative view); we refer to them collectively as the Hawley Creek thrust system. Additional major thrusts within this part of the deeper thrust system include the Cabin thrust (Skipp, 1988) and Freeman–North Fork thrusts (Lonn et al., 2016). In the footwall of the Hawley Creek thrust system within and east of the Lemhi Arch, Lund (2018) and Parker and Pearson (2020, 2021) showed that displacement within the Thompson Gulch thrust system, which involves Devonian and younger strata, is crosscut by the deeper Hawley Creek thrust system in the Beaverhead Range. Approximately 65 km to the northeast in southwestern Montana, Tysdal (1988) showed that early thrusts detached in Mississippian limestones were subsequently folded and crosscut by the basement-involved Jake Canyon thrust. The observation that early thrusting within and east of the Lemhi Arch had a basal detachment in Devonian and younger rocks that predated deformation along deeper Mesoproterozoiclevel detachments (Parker and Pearson, 2021) provides an opportunity to evaluate the timing of the transition between the two thrust systems.

Exhumation of folded Mesoproterozoic Belt quartzites in the Lemhi Range in the hanging wall of the Hawley Creek thrust system is constrained by (U-Th)/He zircon thermochronology at 68–57 Ma (Hansen and Pearson, 2016). Garber et al. (2020) provide detrital evidence from the foreland basin that unroofing of ca. 1380 Ma plutons that are currently exposed within the Hawley Creek system (Lemhi Range and Salmon River Mountains) and in the hanging walls of the Freeman–North Fork thrusts (northern Beaverhead Mountains) occurred at ca. 83 Ma; motion on the Hawley Creek thrust in the Beaverhead Range was interpreted to result in erosion of ca. 500 Ma plutons near the Idaho-Montana border by ca. 68 Ma. These results suggest a general eastward progression of thrustrelated exhumation ranging from Early Cretaceous thrusting in central Idaho to Late Cretaceous thrusting in eastcentral Idaho and southwestern Montana.

The potentially oldest preserved foreland basin strata in southwestern Montana belong to the Morrison Formation, which represents muddy fluvial systems inferred to record deposition east of the flexural forebulge (DeCelles and Giles, 1996) in the backbulge depozone (DeCelles, 2004; Fuentes et al., 2009; Fuentes et al., 2011; Quinn et al., 2018). Eastward migration of the forebulge is hypothesized to explain a regionally extensive Late Jurassic–Early Cretaceous unconformity that is overlain by unequivocal foredeep deposition of the Aptian–Albian Kootenai Formation (DeCelles, 1986; Schwartz and DeCelles, 1988; Fuentes et al., 2011). The Kootenai Formation in southwestern Montana is dominantly nonmarine, with fluvial conglomerate, sandstone, and mudstone interbedded with lacustrine carbonate. Transverse fluvial systems delivered sediment from hinterland terranes and the fold-thrust belt located to the west of the basin, whereas axial systems transported eroded Jurassic continental strata from south of the basin toward the Boreal Sea in Canada (Rosenblume et al., 2021b).

The Blackleaf Formation overlies the Kootenai Formation in southwestern Montana (Fig. 2); the contact is largely conformable. In accordance with previous studies (e.g., Suttner et al., 1981; Schwartz and DeCelles, 1988), we follow the lithostratigraphic nomenclature established by Schwartz (1972). In successive chronostratigraphic order, Blackleaf Formation members are simply named A, B, C, and D. In northwestern and eastern Montana and Wyoming, coeval strata are assigned other names (Fig. 2). Published sedimentological summaries for the Blackleaf Formation are provided by Schwartz (1982) and Dyman and Nichols (1988) based on detailed work from Schwartz (1972) and Dyman (1985).

The Blackleaf Formation represents continued foredeep and forebulge deposition with more rapid basin filling and accelerated subsidence compared to the underlying Kootenai Formation (Fuentes et al., 2011). The lowermost member, Blackleaf-A (13–87 m thick; Dyman and Nichols, 1988), is characterized by lithic and quartz-arenitic sandstone with interbedded mudstone, whereas the overlying Blackleaf-B (17–70 m thick) is dominantly black shale with an increasing sand fraction from east to west. Overall, the A-B series records fluvio-deltaic sedimentation followed by anoxic marine conditions as the Western Interior Seaway flooded the basin (Suttner et al., 1981). Rosenblume et al. (2021a) provide detrital zircon U-Pb data from Blackleaf-A that indicate continued sedimentary recycling from the Idaho fold-thrust belt into the western basin but a primarily Appalachian contribution from transcontinental river systems in eastern basin sandstones.

Blackleaf-C (4–96 m thick) contains fine- to coarsegrained lithic sandstone with sparse interbedded siltstone and mudstone (Dyman and Nichols, 1988; Schwartz, 1972). Sandstones exhibit low-angle cross-bedding and channel geometries. In the west, the contact between B and C is gradational and characterized by increasing sand content. Local sharp-based deltaic sand bodies comprise the lower C sandstone transition upward into paralic and alluvial sandstone, representing overall progradation and a transition from marine to nonmarine deposition. Marine deposition for the lower Blackleaf-C is bolstered by the presence of marine invertebrates and glauconite, whereas lenticular channelized sandstone and sparse woody debris indicate alluvial deposition in the upper Blackleaf-C. A general eastward thinning is observed across the basin; however, local thickness variations are attributed to focused deposition at delta fronts and variable placement of gradational contacts in different sections (Schwartz, 1972). Overall, Blackleaf-C is more texturally and compositionally immature than Blackleaf-A with a greater variety of lithic fragments (Schwartz, 1972; Suttner et al., 1981).

Blackleaf-D conformably overlies Blackleaf-C and is the thickest of all Blackleaf units (270–380 m); it is largely typified by fine-grained, muddy lithologies with interbedded thin lenticular lithic and volcaniclastic sandstone (Schwartz, 1972; Dyman and Nichols, 1988). Pebble-cobble conglomeratic lenses are sometimes present, and fossils are rare. Scattered, thin coal beds indicate organicrich, low-energy terrestrial deposition in the west, and thin limestone beds are occasionally present in the east. Lower Blackleaf-D sediments are attributed to muddominated, low-energy fluvial systems that grade upsection into more coastal plain and/or brackish settings in the west and shallow marine settings in the east, reflective of minor transgression. Blackleaf-D is texturally immature, with widely variable quartz, chert, and volcanic fractions documented across the basin. Volcaniclastic detritus is common along with an increase in finegrained plutonic detritus in the east. Based on compositional data, Suttner et al. (1981) suggest sediment sourcing from recycled Paleozoic–Mesozoic strata in Idaho, specifically citing pebbles that indicate erosion of upper Paleozoic limestone and quartz arenite. Other models infer volcaniclastic sources, specifically from the Boulder batholith (e.g., Schwartz, 1972; Schwartz, 1982), however, that igneous body is now known to be younger than the Blackleaf Formation (81–76 Ma; Scarberry et al., 2021), precluding it as a sediment source.

The upper Blackleaf Formation and its stratigraphic equivalents in northern Montana and Wyoming (Fig. 2) are considered to be Albian–Cenomanian–Turonian in age (ca. 113–89 Ma). Zartman et al. (1995) dated volcanogenic zircons from bentonites in the Vaughn member (Blackleaf-D equivalent) in southwestern Montana and assigned a pooled age range of 97–95 Ma. More recently, Singer et al. (2020) re-evaluated several bentonite ages in central Montana and Wyoming using high precision ⁴⁰Ar/39Ar analysis of sanidine crystals. Bentonites within the Taft Hill (Blackleaf-C equivalent) and lowest Vaughn Members in central Montana are dated to 103.08 ± 0.11 Ma and 102.68 ± 0.07 Ma, respectively. The Muddy Sandstone in Wyoming returns an age of 101.23 ± 0.09 Ma, and two Mowry Shale samples in Wyoming come in at 98.17 ± 0.11 Ma and 97.52 ± 0.09 Ma (Fig. 2).

Detrital zircon studies provide similar maximum depositional ages (MDAs), suggesting that MDAs in the Blackleaf Formation could approximate its true depositional age. Gentry et al. (2018) report MDAs from the Cokeville, Sage Junction, and Aspen formations (Blackleaf-D equivalent) in western Wyoming of 101.5 ± 0.7, 101.5 ± 0.5, and ca. 98.9 ± 0.5 Ma, respectively. Fuentes et al. (2011) provide an MDA near the base of the Vaughn Member in northern Montana of ca. 97.39 +1.49 −0.46 Ma, and Quinn et al. (2018) reports an MDA for the upper sandstone of the Flood Member (probably equivalent to Blackleaf-C based on their stratigraphic section) of 103 ± 1.7 Ma. Rosenblume et al. (2021a) provide multiple MDAs for the underlying Blackleaf-A in southwestern Montana, some from the same stratigraphic sections as this study, with ages ranging from ca. 113 Ma to ca. 106 Ma. Rosenblume (2021) gives ages for the overlying Frontier Formation, with deposition initiating at ca. 91 Ma. All together, these detrital ages consistently indicate Blackleaf-C deposition initiated ca. 106 Ma, was ongoing at ca. 103 Ma, and Blackleaf-D deposition commenced shortly after but by at least ca. 102.7 Ma and continued through at least 92 Ma.

Localities in southwestern Montana were selected to account for west-to-east variations in depositional environment, grain size, and lithology. Samples for detrital zircon geochronology were taken from the lowermost sandstones in Blackleaf-C and -D (Figs. 1 and 2). For sampling context, stratigraphic sections were measured and described, and samples are identified with abbreviations corresponding to their sample locality and the sampled stratigraphic level in meters. Measured sections and GPS coordinates are given in Figure S1 and Table S1¹ in the Supplemental Material. Blackleaf-C sandstone was identified in the field as the first sandstone above the mudstone-rich Blackleaf-B, whereas Blackleaf-D samples were taken at the first appearance of a thin lithic or volcaniclastic sandstone above Blackleaf-C. Blackleaf-C was not identified at sample site LPB, and Blackleaf-D was not identified at locality RRB (Fig. 1).

Five detrital zircon samples were collected, each from Blackleaf-C and -D. Six of the samples were split, and one half was crushed using a jaw crusher and disk mill at the University of Iowa; the other half was disaggregated using selective fragmentation (SELFRAG) at the U.S. Geological Survey (USGS) Mineral Separation Lab in Denver, Colorado, as part of an independent comparative methods test. Mineral separation procedures and descriptions of test results are presented in the Supplemental Material (see ^{footnote 1}).

Detrital zircon ages were measured via rapid laser ablationinductively coupled plasma mass spectrometry (LA-ICPMS) at the University of Arizona’s LaserChron Center using a NuPlasma multicollector ICPMS (Sundell et al., 2020) with n = 600 grains analyzed per sample and a spot size of 20 microns. Detrital zircon cores were targeted and exposed via polishing. Duluth Gabbro (FC) zircons (ca. 1099 Ma) and R33 (ca. 420 Ma) zircons were used as standards. Common Pb corrections were made using Stacey and Kramers (1975). Analyses with greater than 10% discordance or 5% reverse discordance were filtered out only for grains older than 600 Ma due to difficulties in distinguishing lead loss in younger grains. A threshold of 900 Ma was used to select best ages between the ²³⁸U/²⁰⁶Pb (<900 Ma) and ²⁰⁶Pb/²⁰⁷Pb (>900 Ma) systems. Data reduction was completed using AgeCalcML (Sundell et al., 2020). Samples that displayed an obvious negatively skewed tail in their youngest grains (e.g., Spencer et al., 2016) were reevaluated to avoid those with potential lead loss. Care was taken to only eliminate grains that had high error and little to no temporal overlap with the adjacent grain so as to minimize impact on maximum depositional age calculations. Analyzed grain ages, including rejected ages for negatively skewed tails, are reported in the supplemental tables (see ^{footnote 1}).

Quantitative determination of relative proportions of recycled sediment sources was attained through inverse Monte Carlo mixture modeling. Existing detrital zircon data were compiled from clastic sedimentary rocks representing the entire stratigraphic sequence beneath the Blackleaf Formation in Montana and Idaho. Modeled strata were grouped by similar detrital zircon U-Pb age distributions into ten groups: Mesoproterozoic Belt Supergroup, Neoproterozoic, Cambrian, lower Paleozoic (Ordovician, Devonian, and Mississippian), upper Paleozoic (Pennsylvanian–Permian), Triassic, Jurassic, and three Cretaceous groups (Kootenai types I, II, and Blackleaf A-West). The DZMix modeling software runs an inverse Monte Carlo simulation, mixing the potential source zircon spectra into random proportions during the process of 10,000 iterations to find a combination that best fits the probability density plot (PDP) of the desired Blackleaf sample (Sundell and Saylor, 2017). The best-fit model reports the possible sources in their relative proportions and is evaluated using a cross-correlation coefficient (R²) for similarity between the Blackleaf PDP and the modeled PDP. An R² > ~0.70 was accepted as indicating all potential zircon sources are reasonably represented, whereas an R² < ~0.70 indicated that some potential zircon sources may be missing (Sundell and Saylor, 2017). Most samples displayed an R² > ~0.80.

Ultimately, we eliminated Cretaceous and Jurassic strata from the source database along with Cretaceous and Jurassic zircon grains from our samples due to low statistical correlation (R² values) from resulting outputs. Sufficient evidence exists for sourcing of Jurassic and Cretaceous grains directly from igneous rocks related to the Cordilleran magmatic arc (e.g., Gaschnig et al., 2017; Rosenblume et al., 2021b), rendering them inappropriate for mixture modeling, which is designed to evaluate recycled sedimentary sources. A more detailed discussion of arc-related sources is presented below. Mixture models consistently show the best R2 values when only including grain ages and potential sources that are Triassic and older (>201 Ma; Table S2 [^{footnote 1}]). A lack of matching age modes for Jurassic and younger grains between our samples and most underlying Cretaceous and Jurassic strata further supports a 201 Ma cutoff in the models. As expected, eliminating younger grains and source strata greatly increased R² values and identified recycled sedimentary source components that better align with visual inspection of probability density plot age modes. Potential biases common to all detrital zircon provenance analyses, such as variable zircon fertility in the source rocks, hydraulic sorting during sediment transport and deposition, and mineral separation or sampling biases, may be present but no specific evidence for bias has been identified in the data.

Thin sections for upper Blackleaf Formation sandstones (N = 12) were prepared by Spectrum Petrographics and analyzed using the Gazzi-Dickinson point-counting method (n = 400) to calculate grain abundances and approximate sandstone composition (e.g., Ingersoll et al., 1984). Petrographic analysis was conducted using a Nikon Eclipse 50i POL polarizing microscope, PetrologLite x64 software, and an automated stepping stage. Grain abundance data are presented in Table S3 (^{footnote 1}). Chert is considered a lithic fragment, and volcaniclastic lithic fragments were identified based on the presence of thin, randomly oriented microlites in a rounded clast. Very fine–scale chalcedonic silica, here classified as chert, also bears a similar appearance, sometimes leading to difficulty in classification.

Sandstone point-count results corresponding to each detrital zircon sample are plotted on ternary diagrams (Fig. 3), as quartzfeldsparlithics (QFL), sedimentary (LS), volcanic (LV) and metamorphic (LM) lithics as fractions of the total lithics, and limestone (Lsl), chert (Lsc), and sandstone and mudstone (Lss) lithics as fractions of the total sedimentary lithics. Overall, our point-count data align well with previous compositional data (Suttner et al., 1981). Detrital zircon age distributions for each sample are visualized and discussed in terms of specific age modes from relative probability plots (Fig. 4) and DZMix model outputs of relative contributions from potential source units (Fig. 5).

At the westernmost Blackleaf-C locality, MM1-168 is a lithic sandstone (41% lithics), of which 91% of lithic grains are sedimentary and 9% are volcanic. Chert makes up 49% of the sedimentary lithics, with sandstone and mudstone fragments comprising 42% and limestone 9% (Fig. 3). Twenty-seven meters up-section, MM1–195 is a fine-grained sublithic sandstone (25% lithics). Of those, 82% are sedimentary and 18% are volcanic, and 63% of the sedimentary lithics are sandstone and mudstone with 37% chert fragments. The detrital zircon probability plot for MM1-195 (n = 578) has prominent age modes of 106, 168, 1794, and 1850 Ma, with minor modes at 268, 1060, ~1400, 1927, 2087, and 2709 Ma (Fig. 4). 20.0% of grains are derived from the Cordilleran arc. Mixture modeling identifies an ~55% recycled contribution from lower Paleozoic (Ordovician, Devonian, and Mississippian) strata, 28% from Cambrian strata, and minor input from upper Paleozoic (Pennsylvanian–Permian) strata (8%) with R² = 0.87 (Fig. 5).

In the northern part of the study area, NMB-175, also a fine-grained lithic sandstone, is composed of 56% lithic fragments (Fig. 3). Only 16% are volcanic, with the remaining 84% sedimentary composed of primarily sandstone and mudstone (56%) with minor chert (24%) and limestone (15%). A detrital zircon probability plot (n = 542) displays a prominent age mode of 104 Ma, a secondary mode at 1787 Ma, and minor age modes at 114, 167, 434, 555, 1172, 1468, 1713, 1873, 2115, 2586, and 2728 Ma (Fig. 4). 43.7% of grains are derived from the Cordilleran arc. Recycled sedimentary source contributions identified via modeling in decreasing order are Cambrian (~35%), Triassic (~23%), upper Paleozoic (Pennsylvanian–Permian; ~20%), lower Paleozoic (Ordovician, Devonian, and Mississippian; ~12%), and Neoproterozoic (~10%) with R² = 0.69 (Fig. 5).

In the southeastern portion of the basin, LINCMB-68 is a lithic sandstone (37% lithics) with 63% mudstone and 23% limestone lithics (Fig. 3). Detrital zircon geochronology (n = 507) shows a prominent age mode at 102 Ma, with minor modes at 162, 172, 393, 408, 427, 533, 553, 622, 1066, 1501, 1651, 1783, 1875, and 2709 Ma (Fig. 4). 7.7% of grains are derived from the Cordilleran arc. Primary recycled sedimentary sources identified via mixture modeling are Triassic (~52%) and upper Paleozoic (~23%) strata with R² = 0.80 (Fig. 5).

The easternmost sample, CMB-85, is a lithic sandstone (35% lithics), of which 70% are limestone grains (Fig. 3). Detrital zircon geochronology (n = 475) shows a prominent Cretaceous age mode at 105 Ma, followed by minor 165, 420, 596, 1065, 1387, 1505, 1666, 1780, and 2776 Ma age modes (Fig. 4). 5.1% of grains are derived from the Cordilleran arc. Mixture modeling identifies significant recycled contribution from Triassic strata (~50%), with minor input from upper Paleozoic (Pennsylvanian–Permian; ~25%) and Neoproterozoic (~10%) strata with R² = 0.79 (Fig. 5).

In the central portion of the basin, RRB-65 is a lithic sandstone (80% lithics) of which 84% are sedimentary and 16% are volcanic. Forty percent of the sedimentary lithics are sandstone and mudstone and 51% are chert (Fig. 3). Detrital zircon ages (n = 569) have prominent modes at 105 and 1832 Ma, with secondary modes at 116, 169, 1050, 1779, 1832, 1917, 2072, 2580, and 2691 Ma (Fig. 4). Fiftytwo percent of grains are derived from the Cordilleran arc. Mixture modeling identifies recycled sediment source contributions of ~50% lower Paleozoic (Ordovician, Devonian, and Mississippian) and ~30% Cambrian strata into the basin with R² = 0.91 (Fig 5).

At the westernmost locality, MM1-240 is a lithic sandstone (33% lithics), with sedimentary grains making up 95% of all lithics, of which 37% are chert, 34% are sandstone and mudstone, and 28% are limestone fragments (Fig. 3). MM1-278, collected 38 m upsection, is also a lithic sandstone (70%), of which 90% are sedimentary and 7% are volcanic. Sedimentary lithics are 43% chert, 35% sandstone and mudstone, and 22% limestone. Detrital zircon geochronology of MM1-240 (n = 515) displays prominent age modes of 169, 427, 441, 1652, 1778, and 1842 Ma, and minor 119, 551, 589, 620, 1055, 1508, 2023, and 2726 Ma modes (Fig. 4). 2.5% of grains are derived from the Cordilleran arc. Mixture modeling identifies roughly even recycled contributions from Triassic and upper Paleozoic (Pennsylvanian–Permian) strata (45% and 40%, respectively) with little other input with R² = 0.82 (Fig. 5).

In the northern region, NMB-250 is a volcaniclastic sandstone with 70% lithic fragments that are mainly sedimentary (54%) and volcanic (43%; Fig. 3). Sedimentary lithics are dominantly sandstone and mudstone (58%) with minor chert (37%). Detrital zircon geochronology (n = 575) reveals three prominent age modes at 102, 166, and 1781 Ma (Fig. 4). Minor modes at 115, 1431, 1837, 2588, and 2954 Ma are also present. 32.5% of grains are derived from the Cordilleran arc. Mixture modeling identifies recycling of exclusively Cambrian strata (95%) as the sediment source with R² = 0.82. Triassic source strata and Triassic zircons from this sample were not included in the preferred model as they produced a sub-par R² value of 0.61 (Fig 5).

In the southeastern part of the study area, LINCMB-85 is a volcaniclastic sandstone with 46% lithic fragments that are mainly sedimentary (77%) and volcanic (36%; Fig. 3). Sedimentary lithics are sandstone and mudstone (55%), chert (28%), and limestone (17%). Detrital zircon geochronology (n = 590) shows a single prominent age mode of 103 Ma (Fig. 4). 91.7% of grains are derived from the Cordilleran arc. Because the number of grains that are >201 Ma is so low in this sample (n = 49), mixture models were unreliable with R² = <0.51. Visual inspection of age modes (344, 431, 627, 1498, 1782, and 1835 Ma) suggests a recycled Cambrian, lower and upper Paleozoic, and Triassic sedimentary source. This result is consistent with mixture modeling results from nearby sample CMB-98.

In the easternmost basin, CMB-98 is a lithic sandstone (32% lithics), of which 94% are sedimentary and only 6% are volcanic (Fig. 3). Sedimentary lithics are notably chert (42%) and sandstone and mudstone (42%), with minor limestone (16%). Detrital zircon geochronology (n = 533) displays a single prominent age mode of 104 Ma with minor 115, 165, 418, 1063, 1437, 1783, 2686 Ma modes (Fig. 4). 12.4% of grains are derived from the Cordilleran arc. Mixture modeling identifies recycling from almost all modeled strata, including upper Paleozoic (Pennsylvanian–Permian; ~30%), Triassic (~23%), Cambrian (20%), Neoproterozoic (~18%), lower Paleozoic (Ordovician, Devonian, and Mississippian; ~12%) with R² = 0.84 (Fig. 5).

In the southern part of the study area, LPB-145 is a volcaniclastic sandstone that contains 30% lithic fragments with 64% sedimentary and 36% volcanic. Of the sedimentary lithics, 54% are chert, and 44% are sandstone and mudstone (Fig. 3). Detrital zircon geochronology (n = 588) reveals prominent age modes of 92 and 103 Ma, with minor modes at 157 and 1835 Ma (Fig. 4). 72.0% of grains are derived from the Cordilleran arc. Mixture modeling identifies recycled lower Paleozoic strata (Ordovician–Devonian–Mississippian; ~40%) as the primary sediment source, with minor influx from Neoproterozoic (~18%), Cambrian (~16%), and Triassic (15%) strata with R² = 0.81 (Fig. 5).

The most appropriate method for calculating maximum depositional ages (MDAs) from detrital zircon U-Pb data is still a topic of debate; however, consensus is trending toward using multiple overlapping grain ages from high-n data sets (Coutts et al., 2019; Herriott et al., 2019; Sharman and Malkowski, 2020). Detrital zircon MDAs for Blackleaf-C and -D are presented in Table 1, with preferred MDA in bold, and are calculated on the basis of youngest single grain (YSG), youngest single concordant grain (YSCG), youngest graphical peak (YGP), and youngest single population (YSP) (Dickinson and Gehrels, 2009). YSP is the weighted average of the youngest cluster of grains with a mean square weighted deviation (MSWD) of ~1 (Coutts et al., 2019). This method was preferred for most samples due to (1) an abundance of near-depositional grains, (2) younging upsection of YSP MDAs indicating they likely approximate true depositional age (TDA), and (3) the young (ca. 100 Ma) nature of the grains necessitating a higher-n MDA calculation to reduce the impacts of potential lead loss that is difficult to account for in young (<600 Ma) grains. High-resolution zircon sampling throughout the basin between this study and the underlying Blackleaf-A and Kootenai formations (Rosenblume et al., 2021a; Rosenblume et al., 2021b) and the overlying Frontier Formation (Rosenblume et al., 2018) requires a method of calculating MDAs that prioritizes the preservation of stratigraphic relationships in sampled sections. Honoring age constraints from below and above in the same stratigraphic section also ruled out methods such as YPP and YSCG that occasionally produced ages significantly older than MDAs from other methods or from underlying strata.

Blackleaf-C and -D YSP ages display a consistent ~.02–1.1 m.y. difference within the same stratigraphic section. Upsection younging in detrital zircon MDAs has been shown to be consistent with arc-basin connectivity and to closely approximate TDA (Schwartz et al., 2017; Daniels et al., 2018; Johnstone et al., 2019). Blackleaf-C MDAs cluster ca. 106–101 Ma and Blackleaf-D ca. 105–92 Ma. MDAs also young away from the fold-thrust belt, consistent with clastic wedge progradation sourced from exhumed strata due to active thrusting (e.g., Mars and Thomas, 1999), which was originally inferred from detailed sedimentologic studies of the Blackleaf Formation (Schwartz, 1972). Because all Blackleaf-C and -D samples, except LPB-145, are from the base of each interval, the onset of Blackleaf-C and -D deposition in southwestern Montana is constrained to 106 Ma and 105 Ma in the west to 105 Ma to 100 Ma in the east, respectively. Blackleaf-D sedimentation continued until at least 92 Ma, determined from sample LPB-145 taken near the top of Blackleaf-D. These MDAs are consistent with previous geochronological work and detrital studies in Wyoming and Montana (e.g., Zartman et al., 1995; Quinn et al., 2018; Singer et al., 2020).

To evaluate the spatial component of the mixture model results, tighter geographic constraints must be placed on the modeled source rocks. Based upon prior work on the timing of deformation in the Sevier belt, we assume sources for the Blackleaf Formation to be confined to Idaho; however, the ability to parse out distinct recycling signatures via mixture modeling allows for even stronger spatial constraints on provenance regions. Spatial constraints on the modern geographic extent of potential recycled sedimentary sources are described here from oldest to youngest.

Mesoproterozoic rocks from the Belt Basin are exposed throughout northeastern Idaho and western Montana (Fig. 1), and in Idaho are typified by a strong detrital zircon age mode between 1710 and 1740 Ma and subsidiary age peak at ca. 1440 Ma (Link et al., 2016). Because these age peaks are not prevalent in any Blackleaf sample (Fig. 4) and none of the models require contributions from Belt strata, we conclude that thrust sheets bearing Belt strata were either in the early stages of unroofing shallower stratigraphy or not unroofing at all during deposition of the Blackleaf Formation.

Neoproterozoic and Cambrian units in Idaho are the oldest strata exhumed by the fold-thrust belt west of the Lemhi Arch (Brennan et al., 2020), with scattered exposures between the northern end of the Lost River Range and the WISZ (Fig. 1; Lund et al., 2003; Brennan et al., 2020). These rocks were formerly more widespread in central Idaho, but were overprinted and buried by Late Cretaceous and Paleogene extrusive and intrusive magmatism of the Idaho batholith and Challis Volcanic Group (Ma et al., 2016; Stewart et al., 2017). Correlative Neoproterozoic and Cambrian rocks are also exposed in southeastern Idaho on the Paris-Putnam thrust sheet (Yonkee et al., 2014; Link et al., 2017), but detrital zircon data from Blackleaf-equivalent and younger foreland strata from southeastern Idaho and western Wyoming indicate that those rocks were still buried during Blackleaf time (Malone et al., 2017; Gentry et al., 2018); this rules out southeastern Idaho as a source of recycled Cambrian detritus in southwestern Montana during Blackleaf time.

Neoproterozoic units in central Idaho are dominated by a Ectasian-Stenian “Grenville” (ca. 900–1300 Ma) signature that is limited in our samples—especially those bearing Cambrianderived grains—suggesting that zirconbearing Neoproterozoic units were not a source during Blackleaf time. However, both margin-wide (Yonkee et al., 2014; Linde et al., 2017; Matthews et al., 2018) and central Idaho studies (Brennan et al., 2020) reveal a regional Middle Cambrian detrital zircon signature typified by a 1780 Ma peak, which is present in several of our samples; we thus consider recycled Middle Cambrian rocks a likely source for some Blackleaf Formation detritus.

Lower Paleozoic (Ordovician, Devonian, and Mississippian) strata are exposed across the central Idaho fold-thrust belt (Fig. 1) and are typified by a detrital zircon signature with a prominent ca. 1840 Ma age mode and minor ca. Ma 1920 and Ma ca. 2080 modes (Baar, 2009; Beranek et al., 2016; Ma et al., 2016; Brennan et al., 2020). Lower Paleozoic strata were a significant sediment source for the Kootenai Formation and Blackleaf-A (Rosenblume et al., 2021a; Rosenblume et al., 2021b). Notably, Mississippian strata grade from carbonatesiliciclastic Antler-derived turbidite flysch in the western foldthrust belt (Link et al., 1996) to a carbonate bank on the Lemhi Arch (Sandberg et al., 1975; Skipp et al., 1979b; Link et al., 1988), with Upper Mississippian units comprising almost 1000 m of outer carbonate bank in eastern Idaho (Skipp et al., 1979b).

Upper Paleozoic (Pennsylvanian and Permian) strata are presently exposed in the hanging wall of the Pioneer thrust and along the southern and western flanks of the Lemhi Arch (Fig. 1) but likely were ubiquitous atop lower Paleozoic strata based on exposures across the fold-thrust belt (Skipp et al., 1979a; Skipp et al., 1979b). Anastasio et al. (2004) summarize several lines of evidence for burial of Mississippian strata in the Lost River Range to ~7 km depth during Cretaceous synorogenic deformation, suggesting significant overlying upper Paleozoic and younger strata. Lithologies include both interbedded carbonate and siliciclastic sediments (Skipp et al., 1979b; Mahoney et al., 1991; Rankey, 1997). Detrital zircon signatures for these units are characterized by strong ca. 425 Ma, ca. 1070 Ma, and ca. 1640 Ma age modes (Link et al., 2014; Garber et al., 2020; Leary et al., 2020; Rosenblume et al., 2021b).

Triassic strata are sparsely exposed in central Idaho, and existing detrital zircon data are scarce (Fig. 1). It is reasonable, however, to assume that Triassic rocks were once present atop thrust sheets currently eroded to Paleozoic levels, because the Triassic Dinwoody marine basin likely extended across south-central Idaho (Paull and Paull, 1994; Hofmann et al., 2013) with additional nonmarine Triassic material deposited outside of its mapped extent. Isopach mapping of Triassic strata across Idaho, Montana, Wyoming, and adjacent states suggests that while they are not preserved in central Idaho, contours of >1 km are truncated along the erosional limit in eastern Idaho and imply Triassic deposition continued farther to the west (Peterson, 1988). The central Idaho fold-thrust belt is currently unroofed to Paleozoic levels (Fig. 1); therefore, a lack of Triassic strata can be explained by erosion. Unroofing of Triassic strata is documented to have occurred in southeastern Idaho during Blackleaf time in the Wyoming foreland (Gentry et al., 2018). Detrital zircon age modes in the Triassic strata are typified by multiple 400–600 Ma age peaks and a 1640 Ma peak seen in both Triassic zircon data from this study area (Rosenblume et al., 2021a) and in southern British Columbia (Gehrels and Pecha, 2014).

For reasons discussed earlier, the high percentage of Late Jurassic and middle Cretaceous grains, particularly <110 Ma, in most upper Blackleaf samples are interpreted to be mostly derived from Cordilleran arc-related rocks and not sedimentary recycling. The most proximal preserved Cordilleran arc is the Idaho batholith in central and western Idaho. The oldest part of the batholith is the ca. 110–100 Ma suture zone suite (Gaschnig et al., 2010, 2017), which was spatially coincident with and deformed by the WISZ from 104 to 92 Ma (Giorgis et al., 2008; Braudy et al., 2017), partially coincident with deposition of the upper Blackleaf Formation. Zircon ages from the suture zone suite match the upper Blackleaf Formation (Giorgis et al., 2008; Unruh et al., 2008; McKay et al., 2017) with prominent 113–118 Ma, 100–110 Ma, and ca. 92 Ma ages present, notably with the 110–100 Ma suites occupying the easternmost portion of the suite and the inferred locus of active deformation during Blackleaf deposition (Braudy et al., 2017; Gaschnig et al., 2017; McKay et al., 2017). While these complexes were emplaced and deformed at depth, transpression was occurring in the shear zone during and after emplacement, and contemporaneous surficial volcanism and/or erosion of overlying plutonic bodies and earlier volcaniclastic deposits linked to SRSZ magmatism likely contributed significant zircons to the upper Blackleaf Formation. In addition to the intrusive component of the Idaho batholith, Hannon et al. (2021) recently identified bentonites in the distal foreland basin that are attributed to volcanic activity along the SRSZ spanning 130–100 Ma.

Jurassic zircon grains in the Blackleaf Formation and equivalent units occur at a higher percentage in strata to the north than to the south. For example, Jurassic zircon ages (largely ca. 160 Ma) comprise 4% of upper Blackleaf grains in southwestern Montana, whereas the northernmost sample in Blackleaf-D has 10% Jurassic ages. An upper Blackleaf equivalent sandstone in the Mowry Formation to the south in Wyoming has 0% Jurassic grains (May et al., 2013), whereas detrital zircons from an upper Blackleaf sandstone in northern Montana have 14% Jurassic grains (Fuentes et al., 2011). Rosenblume et al. (2021b) suggest that a dominant age mode of ca. 160 Ma zircons in the underlying Kootenai Formation is a fingerprint of the southern Omineca Belt (e.g., Webster and Pattison, 2018), which was migrating northward and likely located just west of Idaho during the middle Cretaceous according to previous paleogeographic reconstructions (Wyld et al., 2006). Jurassic grain fractions in southwestern Montana decrease through time from the Kootenai Formation to the upper Blackleaf Formation yet increase spatially in the upper Blackleaf toward the north. This suggests that some Omineca-derived sediment may have still been making its way into the basin during upper Blackleaf time but was largely overshadowed by more proximal arc activity in the SRSZ.

Multidimensional scaling (MDS) plots are a useful way to visualize sample groupings. Blackleaf-C and -D detrital zircon PDPs are crosscorrelated against one another and plotted on a unitless spatial plane with respect to how similar or dissimilar they are. When plotted along with potential sources, groups of similar provenances can be ascertained (Fig. 6). Like the mixture models, samples and sources plotted are evaluated only based on Triassic and older grains; therefore, the mixture models and MDS plot only address sedimentary recycling. As such, the MDS plot is used to visualize similarities and dissimilarities among Blackleaf samples themselves as well as exposed units in the Idaho fold-thrust belt and southwestern Montana foreland basin. An MDS plot that includes all grain ages is in the Supplemental Material (see ^{footnote 1}) and demonstrates that, when using the entire age spectrum, the relative abundance of arc-derived grains (<200 Ma) controls the positions of the samples, negating the effects of recycled sedimentary provenance.

Provenance Group 1 (Fig. 6) includes samples MM1-195 and RRB-65 from Blackleaf-C and LPB-145 from Blackleaf-D. A composite plot shows dominant peaks at ca. 1840 Ma, ca. 1910 Ma, and ca. 1780 Ma. Averages from mixture modeling indicate primary sourcing from lower Paleozoic strata (~50%), minor input from Cambrian strata (~25%), and minimal input (<10% each) from various other sources (Fig. 5). Group 1 also shows the most similarity with the lower Paleozoic strata on the MDS plot (Fig. 6). Point-count data reveal that these samples generally have a higher proportion of chert lithics (35%–60%) and lower proportion of limestone lithics (<10%) relative to other Blackleaf samples (Fig. 3).

Sourcing from lower Paleozoic and Cambrian strata is geographically tied to central and/or north-central Idaho west of the Lemhi Arch, where the Idaho batholith and shear zone plutons intruded Cambrian and lower Paleozoic units (e.g., Ma et al., 2016). While the remaining exposures of the fold-thrust belt in northcentral Idaho currently bear both Cambrian and lower Paleozoic strata, igneous overprinting throughout central Idaho hinders our ability to draw definitive source-to-sink pathways between mapped units on thrust faults and upper Blackleaf sample localities. However, Group 1 samples all bear significant (20%–73%) arc grains (Fig. 4), implying a catchment that included recycled sediment sources adjacent to the arc.

Provenance Group 2 (Fig. 6) includes LINCMB-68 and CMB-85 from Blackleaf-C and MM1-240 from Blackleaf-D. Dominant age peaks range between 400 and 600 Ma, ca. 1060 Ma, ca. 1500 Ma, ca. 1660 Ma, and ca. 1780 Ma. Mixture modeling indicates primary sourcing from Triassic (~50%) and upper Paleozoic (~30%; Pennsylvanian–Permian) strata (Fig. 5). Point-count data indicate a composition that is consistently ~60%–66% quartz and 33%–40% lithics, with higher percentages of limestone sedimentary lithics (>30%) relative to other samples (Fig. 3). Low (~5%) concentrations of arc grains in Group 2 suggest sourcing away from the Cordilleran arc in western Idaho. Given that the central Idaho portion of the Idaho-Montana fold-thrust belt had already been exhumed down to lower Paleozoic levels (Group 1; Rosenblume et al., 2021a; Rosenblume et al., 2021b), the interpretation of recycling from younger strata for Group 2 indicates an alternate sediment source exhuming younger strata. This is consistent with eastward thrust propagation into east-central Idaho in regions presently eroded to upper and lower Paleozoic levels or overprinted by the eastern Snake River Plain to the south (Fig. 1).

Provenance Group 3 (Fig. 6) only includes NMB-250 from Blackleaf-D and is dominated by a ca. 1780 Ma peak. Mixture modeling indicates almost exclusive recycling of Cambrian strata (95%) with little to no lower Paleozoic influence. This 1780 Ma signature is consistent with regional Cambrian detrital zircon spectra from Brennan et al. (2020) and Matthews et al. (2018). The absence of Cambrian strata on the Lemhi Arch and restriction of the fold-thrust belt progression to Idaho during Blackleaf deposition limit provenance to a central Idaho source, similar to that of Group 1; this source region now exposes Mesoproterozoic units or is overprinted by the Idaho batholith (Fig. 1). A relatively high arc grain fraction (33%) implies a similar catchment to Group 1.

Provenance Group 4 (Fig. 6) includes NMB-175 from Blackleaf-C and LINCMB-85 and CMB-98 from Blackleaf-D. Samples in this group are typified by age peaks seen both in Groups 2 and 3. Mixture modeling results indicate various levels of contributions from primarily Cambrian, upper Paleozoic, and Triassic units, with a relatively higher Cambrian component in NMB-175 (~35%) and higher upper Paleozoic/Triassic input in CMB-98 (~55%; Fig. 5). The Cambrian zircon signature is consistent with sourcing from central Idaho (Groups 1 and 3), and the upper Paleozoic and Triassic signatures match a younger more inboard source in eastern Idaho (Group 2; Fig. 1). Therefore, mixing of separate catchments is implied. LINCMB-85 from Blackleaf-D did not have sufficient grains that are >201 Ma to produce a successful mixture model; however, similarity in age peaks with the other Group 4 samples in the MDS plot and geographic and temporal proximity all indicate a similar sediment source as CMB-98. The similarity between combined Groups 2 and 3 with Group 4 is visualized on Figure 6 where the PDP of combined Group 2 and 3 signatures (purple line) neatly overlaps Group 4.

Group 1 samples on the MDS plot show similarities to underlying foreland basin strata, including Type II Kootenai (Rosenblume et. al, 2021b) and Blackleaf-A Group 2 (Fig. 6; Rosenblume et al., 2021a), although Kootenai Type II lies slightly closer to the lower Paleozoic strata. Both units were interpreted to be primarily sourced from lower Paleozoic and/or Cambrian units exhumed in thrust sheets in central Idaho. Ascertaining whether our samples represent continued sourcing from those same central Idaho thrust sheets or initial unroofing of Blackleaf-A and Kootenai sediments atop younger thrusts in eastern Idaho is critical to our reconstructions.

Here we outline several lines of evidence against recycling of underlying Cretaceous strata as the primary source for Group 1. First, the progression of lower Paleozoic sources during Kootenai Type II deposition to a larger Cambrian contribution during deposition of Blackleaf-C and -D suggests continued unroofing of older strata rather than recycling of the underlying Cretaceous strata. Second, the overwhelming presence of <110 Ma zircon grains inferred to have been sourced from the Cordilleran arc in western Idaho, which are not seen in either Kootenai or Blackleaf-A zircon signatures (Rosenblume et al., 2021a; Rosenblume et al., 2021b), indicates that a direct connection likely existed between hinterland source regions in western Idaho and the foreland basin. Airfall was also not a major contributor of arc grains as maximum depositional ages from our samples (Table 1), and bentonite ages in the upper Blackleaf (e.g., Singer et al., 2020) are distinctly younger than the prominent Cretaceous age peaks in all samples (Fig. 4). Group 2 (Triassic–Upper Paleozoic) recycled signatures also consistently have much lower <110 Ma zircon fractions, associated with nascent thrusts in eastern Idaho away from the Cordilleran arc, whereas Groups 1 and 3, sourced from lower Paleozoic–Cambrian strata in central Idaho, consistently have higher >110 Ma fractions, suggesting catchments drawing directly from the arc. Direct airfall transport from the arc to the basin would result in samples with either a uniform or random distribution of <110 Ma zircon fractions and not one that follows such a stark provenance distinction.

Sample LPB-145, deposited ~8 m.y. after our oldest Blackleaf samples, also bears a Group 1 (lower Paleozoic–Cambrian) signature as well as 103 and ca. 160 Ma age modes. Because lower Paleozoic input from central Idaho terminated ca. 105–100 Ma, we surmise that due to the high similarity between Group 1, Kootenai Type II, and Blackleaf-A Group 2, LPB-145 likely represents recycling of Cretaceous foreland strata into the basin from a more frontal thrust, possibly as far east as the Hawley Creek thrust system (Fig. 1). LPB-145 also has a dominant 92 Ma age mode in conjunction with a 92 Ma MDA, indicating a significant airfall contribution likely derived from Idaho batholith emergence in central Idaho. However, without higher resolution sampling in the basin during this time interval, conclusions about hinterland tectonics are limited.

Schwartz (1972) first recognized portions of the upper Blackleaf system as a prograding clastic wedge, where fluvial deposition of the lower Blackleaf-D sands in the west was coeval with Blackleaf-C shoreface deposition in the east. Therefore, sediment provenance cannot be evaluated by considering Blackleaf-C and Blackleaf-D as distinct temporal units, but rather they must be considered a singular depositional system separated into time slices. Because we infer that MDAs approximate true depositional age in our samples, they can be used to create a detailed tectonic and sediment transport model for the Blackleaf basin from ca. 106–100 Ma. MDAs and their respective provenance groups for each sample are presented in Figure 7, grouped by unit and plotted to demonstrate their geographic relationships. Samples are plotted by perpendicular distance from the Hawley Creek thrust front and demonstrate that MDAs within Blackleaf-C and -D young to the east. Blackleaf-A sample LPB-15 from Rosenblume et al. (2021a) is included in the sediment transport model because its MDA (105.8 ± 1.4 Ma) and zircon signature are very similar to some nearby samples from the upper Blackleaf Formation. The LPB-15 site (Fig. 1) is located slightly west of the inferred maximum transgression shoreline, and Blackleaf-C was not identified there, so we can consider it to have been deposited very shortly before the onset of Blackleaf-C in the western basin.

Stage 1 is based on the youngest Blackleaf-A sample (LPB-15; Rosenblume et al., 2021a) and the two oldest and westernmost Blackleaf-C samples (MM1-195, RRB-65) that define our Group 1 provenance, all with similar MDAs that are distinctly older than the other upper Blackleaf samples. In the southwestern corner of the basin, Rosenblume et al. (2021a) interpreted Blackleaf-C sample LPB-15 to be sourced from recycled upper Paleozoic–Triassic strata (akin to upper Blackleaf Group 2). Sparse Cretaceous grains (4%) in LPB-15 suggest that the sediment sources were located primarily away from the arc. To the north of the LPB locality at the onset of Blackleaf-C deposition (MM1-195, RRB-65), sediment was sourced primarily from lower Paleozoic strata and minor Cambrian strata in central Idaho, coupled with significant Cretaceous arc-derived detritus (15%–46%) from magmatism along the SRSZ (Fig. 8A). We interpret this pattern of sediment dispersal in the northern basin to represent continued provenance from units sourcing the Kootenai and Blackleaf-A, as well as the exhumation of younger arc sources, whereas the southern basin shows the first appearance of a younger, more easterly sediment source from upper Paleozoic–Triassic strata.

Stage 2 is based on Blackleaf-C samples LINCMB-68, NMB-175, CMB-85 and Blackleaf-D sample MM1-240 (Fig. 8B). Stage 2 was characterized by samples with Group 2 and Group 4 signatures (Fig. 6). Group 2 signatures typify all parts of the basin in the south and indicate recycling of upper Paleozoic–Triassic strata (LINCMB-68, CMB-85, and MM1-240); the Group 4 signature is present to the north and suggests a mixed Cambrian–upper Paleozoic–Triassic source (NMB-175). This shift in the south from primarily lower Paleozoic sediment sources in Stage 1 to upper Paleozoic–Triassic sources in Stage 2, as well as a reduction in <201 Ma arc grains (Fig. 4), suggests renewed, in-sequence fold-thrust belt propagation and uplift that resulted in unroofing of younger strata in east-central Idaho within the upper structural levels above the Lemhi arch (Parker and Pearson, 2021). In contrast, during this same interval to the north, <200 Ma grains and a mix of older and younger sedimentary sources indicate continued unroofing of the arc and older strata that sourced Group 1 as well as an influx from new, younger sedimentary sources. Clastic wedge progradation during this time resulted in the migration of the basal Blackleaf-C shoreface sandstone to the east into the medial and distal parts of the study area, whereas in the proximal region, nonmarine deposition of Blackleaf-D initiated.

Stage 3 is based on the easternmost Blackleaf-D samples LINCMB-85, NMB-250, and CMB-98. While MDAs are like the eastern Blackleaf-C from Stage 2, all samples in this stage were taken upsection from Stage 2 samples and are therefore stratigraphically younger and temporally distinct. In Stage 3, the nonmarine clastic wedge represented by Blackleaf-D continued to prograde into the basin (Fig. 8C) resulting in deposition of the medial and distal Blackleaf-D. At the LINCMB and CMB localities to the south, upsection changes show a shift from Group 2 to Group 4 resulting from increased mixing between older and younger sedimentary sources across the basin. Arc detritus is also found in the distal part of the basin in variable quantities (92% at LINCMB; 13% at CMB). At NMB in the north, a similar upsection trend of increasing flux from Cambrian strata resulted in a switch from Group 4 to Group 3, coupled with a slight increase in arc detritus. The exclusively Cambrian signature in Blackleaf-D at NMB indicates deeper levels of unroofing in central Idaho.

Overall, two distinct provenance shifts between 106 and 100 Ma occurred during deposition of the upper Blackleaf Formation. From Stage 1 to Stage 2, a pulse of recycled Triassic–Upper Paleozoic detritus (Group 2) in the southern part of the basin is consistent with thrust propagation and unroofing of younger material in the eastern Idaho fold-thrust belt, while increasing Cambrian components and consistent arc detritus in the north suggest further unroofing to deeper levels in central Idaho. These shifts occurred synchronously but have independent tectonic implications. Later, from Stage 2 to Stage 3, continued exhumation of the older strata west of the Lemhi Arch resulted in widespread mixing of older and younger sources throughout the basin.

Group 1 samples indicate a continuation of sediment provenance from central Idaho from the Kootenai through the Blackleaf, with an increasing Cambrian component representing deeper levels of unroofing. Conversely, the exhumation of younger strata that sourced Group 2 was likely occurring more inboard to the east. In east-central Idaho, the upper thrust system of the Idaho-Montana fold-thrust belt accommodated major folding and internal thrusting that mainly occurred prior to activation of the structurally deeper thrust system (Parker and Pearson, 2021). Within this upper thrust system, Mississippian–Pennsylvanian strata present in the Lost River Range accommodated significant folding and thrusting (Fisher and Anastasio, 1994; Hedlund et al., 1994; Anastasio et al., 1997), which would have resulted in uplift and erosion of the overlying Pennsylvanian and younger strata, providing a likely sediment source for Group 2 and limiting sourcing of arc grains from western Idaho seen in Group 1. Low-angle thrust faults in Mississippian strata crosscut by the Hawley Creek thrust system on the Lemhi Arch suggest that a Devonian/Mississippian-level shallow thrust system propagated across central Idaho prior to deformation along deeper detachment levels (Lund, 2018; Parker and Pearson, 2021), which is best constrained to ca. 83–68 Ma (Lund, 2018; Garber et al., 2020; Parker and Pearson, 2020, 2021).

Our Group 2 detrital zircon signature suggests derivation from the region above the Lemhi Arch rather than in central Idaho and indicates initial activity of the upper thrust system in east-central Idaho from ca. 105–101 Ma. West of the Lemhi Arch, Upper Mississippian strata are composed of sandy Antler foredeep deposits (Link et al., 1996) that were recycled from lower Paleozoic strata in central Idaho (Beranek et al., 2016); on the Lemhi Arch to the east, these sand bodies thin, and the Upper Mississippian strata consist of >1000 m of carbonate platform rocks (Skipp et al., 1979b). Therefore, unroofing of Upper Mississippian and younger strata west of the Lemhi Arch would produce a lower Paleozoic zircon signature (Group 1) in association with the Group 2 zircon signature. Sourcing from the Lemhi Arch region, however, would result in erosion of carbonate bank debris with limited zircon input from lower Paleozoic sources. Groups 2 and 4 have an upper Paleozoic zircon source; samples within those groups contain >15% limestone lithic fragments (Fig. 3), whereas Groups 1 and 3 with a lower Paleozoic source contain minimal limestone lithics (4%). Altogether, these observations strongly suggest that Blackleaf sandstones defined by Group 2 provenance signatures were eroded from the upper thrust system above the Lemhi Arch in east-central Idaho. Janecke et al. (2000) documented paleo valleys crosscutting the Lemhi and Beaverhead Ranges, sourced from central Idaho to southwest Montana; formation of these paleovalleys preceded deep-seated deformation and Lemhi Arch uplift. These paleovalleys were eroding older, deeper Mesoproterozoic material by the Late Cretaceous, and thus could have initiated in the mid-Cretaceous with the propagation of shallow thrusting across Idaho, providing a direct source-to-sink linkage for the Blackleaf. We interpret that deformation in this upper thrust system led to widespread unroofing throughout central and east-central Idaho (e.g., Parker and Pearson, 2021), explaining the notable present-day absence of Triassic strata.

Significant amounts (>15%) of arc grains in samples bearing a Cambrian or lower Paleozoic source component indicate catchments encompassing both Cordilleran arc-related rocks and Cambrian and lower Paleozoic strata in western and central Idaho respectively. In Stage 3 (102–100 Ma), the Lemhi Arch–derived Group 2 signature was overprinted by an increase in Group 3 (Cambrian) and arc detritus flooding the basin, suggesting an increase in sediment volume from Cambrian strata. This suggests renewed tectonic activity in the hinterland that was relatively stagnant between the Kootenai Formation and Stage 1 (ca. 110–106 Ma; Rosenblume et al., 2021a, 2021b). However, limited Cambrian exposures through central Idaho and overprinting of the Idaho batholith likely conceal a large area where Cambrian rocks once existed (Ma et al., 2016), allowing for the possibility of a nowobscured thrust system, possibly along-strike to the Pioneer thrust, exhuming Cambrian strata.

In nearby southeastern Idaho and Wyoming, Gentry et al. (2018) and Yonkee et al. (2019) documented enhanced fault slip rates, eastward fault propagation, and exhumation of Jurassic through lower Paleozoic strata on the Willard-Paris-Meade thrust system between 105 and 90 Ma, coeval with upper Blackleaf deposition in southwestern Montana. Slip rates on Sevier thrusts also increased margin-wide between 105 and 100 Ma, documented in the Alberta Rockies (Pană and van der Pluijm, 2015), southern Nevada and California (Wells, 2016), and Utah (Quick et al., 2020), likely due to increasing plate convergence rates (Yonkee and Weil, 2015). Overall, our data document an episode of exhumation within the Idaho-Montana fold-thrust belt, with renewed activity in central Idaho and eastward fold-thrust belt propagation above the Lemhi Arch that was synchronous with increased rates of deformation all along the Sevier margin.

The increase in sediment flux from renewed regional exhumation likely influenced basin sediment dynamics. Blackleaf-C and -D record progradation of a clastic wedge across the southwestern Montana foreland basin during a period of otherwise eustatic sea-level rise (Haq, 2014). Vuke (1984) notes that the clastic progradation of the Blackleaf in southwestern Montana did not align with deposition of timeequivalent strata elsewhere and was likely associated with increased sediment supply locally. Because our samples were collected at the first appearance of sandstone in each section (with the exception of LPB-145, excluded from this analysis), and our MDAs in both Blackleaf-C and -D young from west to east across the basin, we can approximate the rate of shoreline progradation during Blackleaf time, averaging MDAs in-section at sample localities to reduce error. The shoreline progradation rate is calculated based on MDAs (Table 1) and their perpendicular distance from the Hawley Creek thrust front to approximate relative positioning in the basin (Fig. 7), equating to ~24 km/m.y. for the upper Blackleaf system. This progradation rate is consistent with the numerical models of similar foreland systems. Zhang et al. (2020) shows a model of continental-scale sediment routing systems with coupled increased uplift and rising eustatic sea level that result in an ~20 km/m.y. increase in shoreline progradation rate from steady-state over ~5 m.y., roughly consistent with the rate increase from the regressive or static depositional systems associated with Blackleaf-A and -B at ca. 106 Ma to complete progradation of Blackleaf-C and -D in Stage 3 at ca. 101 Ma. However, this model is not accommodation-space limited and has a flat floor, which differs from the typical understanding of retroarc foreland basins.

This shift in sediment supply was driven by the renewed thrusting in both central and eastern Idaho. However, adjacent basins did not experience the same rates of shoreline progradation (Vuke, 1984). Different thrust geometries are likely to blame. Flexural load from shallow thrusting in the Idaho-Montana fold-thrust belt’s upper thrust system likely resulted in less accommodation space and more rapid filling of the foreland, aiding shoreline progradation; in contrast, along-strike thrusting in southeastern Idaho encompassing a thicker stratigraphic package (Yonkee et al., 2019) emplaced a greater load with increased foreland accommodation space, explaining the difference in shoreline position relative to the fold-thrust belt between Montana and Wyoming (Fig. 9).

A compilation of our inferred sediment provenance and equivalent foreland basin strata in surrounding areas is shown in Figure 9. In northern Montana, sediment from age-correlative upper Blackleaf sandstones was derived from the Omineca belt and recycling of the Mesoproterozoic Belt Supergroup and possibly Neoproterozoic through lower Paleozoic strata (Fuentes et al., 2011; Quinn et al., 2018). Fuentes et al. (2011) cite specific sources from the St. Mary and Hall Lake thrusts in northern Washington and southern British Columbia, whereas Quinn et al. (2018) do not infer specific thrust sheets as sources. Zircon signatures from both of these studies contain similar Jurassic–Cretaceous arc fractions; however, overall small n-values (n≈100) hinder the ability to make direct comparisons about the relative abundances of the different recycled populations present from the remaining grains (Saylor and Sundell, 2016).

Gentry et al. (2018) presents zircon spectra and mixture models for Blackleaf-equivalent wedge-top strata in southeastern Wyoming affiliated with the Willard-Paris-Meade thrust system. While n-values (n≈120) are similarly small, robust detrital zircon characterization of adjacent fold-thrust belt strata allow for refined provenance modeling, linking deposition to exhumation of lower and upper Paleozoic strata from southeastern Idaho, and lower Paleozoic through Jurassic strata from northeastern Utah.

In the Bighorn Basin to the east of the study area, zircon signatures in the coeval Muddy and Mowry formations are similar to Group 4 in southwestern Montana, with a high Cretaceous arc fraction, a ca. 1780 Ma age peak, and minor Neoproterozoic to Mesoproterozoic (900–1600 Ma) grains (May et al., 2013). Though May et al. (2013) inferred a Belt Supergroup source from Idaho and Montana, more recent work characterizing the detrital zircon signatures of Cordilleran strata and the timing of Belt Supergroup exhumation (e.g., Brennan et al., 2020; Garber et al., 2020), coupled with our study showing similar signatures in the Blackleaf, suggest that these age-equivalent strata in north-central Wyoming were instead sourced from Cambrian and minor Triassic–upper Paleozoic strata in central Idaho. Given the distal position of their study area, we infer that their data represent distal mixing of central Idaho Blackleaf detritus.

New large-n detrital zircon U-Pb geochronology and statistical modeling from upper Blackleaf Formation sandstones in southwestern Montana identify both primary arc-related and recycled sediment sourcing from the Idaho-Montana fold-thrust belt. Detrital zircon maximum depositional ages show that upper Blackleaf Formation deposition in southwestern Montana occurred between ca. 106–92 Ma, with most samples representing the earliest occurrence of a coupled shallow marine shoreface (Blackleaf-C) and fluvial (Blackleaf-D) system that prograded eastward across the basin into the Western Interior Seaway between ca. 106–100 Ma, consistent with previous sedimentological work.

Cretaceous zircon grains were sourced from the Cordilleran arc rocks in western Idaho along the Salmon River Suture Zone. Older grains were sourced via sedimentary recycling of lower Paleozoic–Cambrian strata in central Idaho and upper Paleozoic–Triassic strata above the Lemhi Arch in eastern Idaho.

The oldest zircon samples (constrained by MDAs) show a similar source as the underlying Kootenai Formation and Blackleaf-A, with arc and recycled lower Paleozoic detritus sourced from central Idaho. Between ca. 105–101 Ma, recycled upper Paleozoic–Triassic detritus from east-central Idaho was widespread across the basin, mixing with older units and arc grains in the north. Between 102 and 101 Ma, western Idaho sources were unroofed to exclusively Cambrian levels, bringing a renewed Cambrian signal across the basin mixing with the upper Paleozoic–Triassic detritus from eastern Idaho. Overall, sediment provenance was controlled by the tectonic evolution of these distinct source domains.

Sourcing of arc detritus was driven by arc volcanism near the Salmon River Suture Zone. The structural provenance of lower Paleozoic–Cambrian detritus is more elusive but inevitably tied to tectonic activity in central Idaho. Sourcing of upper Paleozoic–Triassic detritus was driven by deformation in an upper thrust system involving Mississippian and younger strata above the Lemhi Arch in east-central Idaho. Increased exhumation in the central Idaho fold-thrust belt and the upper thrust system in east-central Idaho initiated at ca. 105 Ma, consistent with coeval slip-rate increases observed along the Sevier margin. Increased sediment volume from two discrete provinces and decreased flexural loading drove basinward shoreline progradation compared to adjacent systems. This work emphasizes the power of large-n analyses to allow for more refined provenance modeling, drawing specific connections between source strata and individual basin samples, and refining the tectonic history of western North America.

We thank Chris Holm-Denoma, Benjamin Howard, and Jack Davis for assistance with fieldwork, sample preparation, and analysis. Funding for this research was provided by National Science Foundation (NSF) tectonics grants EAR-1727504 (Finzel) and EAR-1728563 (Pearson), NSF-EAR grant 1649254 (Arizona LaserChron Center), and graduate student awards to C. Gardner from SEPM (Society for Sedimentary Geology). Fieldwork was partly funded by grants from the Iowa Earth and Environmental Sciences Department and Tobacco Root Geological Society to J. Rosenblume.