A Thermal Profile across the Idaho-Montana Fold-Thrust Belt Reveals a Low-Relief Orogenic Wedge That Developed atop a Pre-Orogenic Basement High

Continental fold-thrust belts accommodate horizontal shortening of the upper plate lithosphere at convergent plate boundaries. Due to the self-destructive nature of orogenesis, ancient orogenic belts are often deeply exhumed. In many cases, partial records of the eroded parts of these mountain belts are archived in the foreland basin stratigraphic record, which can be used to reconstruct the early deformational history and geometry of an ancient fold-thrust belt. However, early foreland basin sedimentary rocks may also be cannibalized during the progressive growth of the orogen (e.g., [1]). In some ancient orogenic belts such as the North American Cordillera, the depth of erosion and the structural configuration of the orogenic interior are preserved when regional unconformities developed during or shortly after orogenesis are overlain by sedimentary and/or volcanic rocks, enabling one to see which unit was at the surface when the overlying rocks were deposited. These unconformities can be used to create paleogeologic or “subcrop” maps, which can “see through” the effects of more recent faulting and erosion (e.g., [2–4]). In addition to providing information about the regional geologic relationships that existed during the development of the preserved unconformity, subcrop maps can also be used to constrain the cumulative exhumation magnitude achieved during orogenesis if an approximation of the preexisting stratigraphic thicknesses is known.

When unconformities are not preserved or the thickness of eroded overburden is unknown, exhumation and its associated cooling of rocks can be recorded as a thermal history that is accessed using thermochronometry or paleothermometry (e.g., [5]). Fold-thrust belts experience not only exhumation and cooling of uplifted rocks in hanging walls of thrusts but also burial and heating of rocks in thrust footwalls, resulting in a thermal profile that yields insight into the geometry of the eroded wedge geometry (Figure 1). For instance, erosion of a simple critically tapered wedge is expected to result in rocks that have experienced increasingly higher exhumation depths—and therefore higher maximum recorded temperatures—toward the orogenic interior (Figure 1(a)). More complicated wedge geometries may produce maximum temperature plateaus and/or localized increases in maximum temperatures that reflect where exhumation has occurred over a more steeply dipping section of the basal décollement, such as a regional footwall ramp (Figure 1(b)). Simply put, burial and exhumation produce a thermal profile of an orogenic wedge that approximates the geometry of the eroded surface slope.

Paleothermometers can independently constrain the maximum temperatures—and therefore maximum burial depths—of rocks, thereby providing an empirical way to constrain the geometries of ancient orogenic wedges, even with incomplete knowledge of preexisting stratigraphic thicknesses (e.g., [5–13]). In particular, Raman spectroscopy of carbonaceous material (RSCM) offers a simple way to investigate the thermal histories of carbon-bearing carbonate and fine-grained siliciclastic rocks that are not easily analyzed by conventional low-temperature thermochronometers (e.g., [11, 14, 15]). The temperature-dependent, one-way reaction of organic carbon to graphite [16] is utilized by RSCM to yield robust and reproducible maximum temperature constraints that capture thermal events on geological timescales [17] and at low temperatures <100°C [11, 18]. In addition, the acquisition of maximum temperature data by RSCM is much simpler than for conventional low-temperature thermochronometry (e.g., zircon and apatite (U-Th)/He)) because it is not affected by factors such as radiation damage, partial retention, and annealing (cf. [19]). This makes the RSCM method useful for independently constraining the geometries of orogenic wedges by offering data regarding maximum burial depths.

The Idaho-Montana fold-thrust belt is part of the Sevier and Laramide belts of the northern U.S. Rocky Mountains (e.g., [2]). Prior workers constructed a paleogeologic map of a large portion of the Idaho-Montana fold-thrust belt that approximates the erosion level in early Eocene time, shortly after the end of orogenesis [4]. In addition, limited peak temperature data are available across the fold-thrust belt (primarily conodont color alteration index (CAI) and vitrinite reflectance (R_o) data; [20, 21]). However, these data have not been evaluated in a regional context, and more data are needed in order to do so. Recent work on the Mesozoic record of orogenesis recorded by foreland basin rocks in southwestern Montana [22–24] has also bolstered our understanding of the timing and spatial pattern of deformation and erosion within the fold-thrust belt. Using RSCM to constrain maximum temperatures along a >350 km long transect and quantifying exhumation magnitude based on an updated paleogeologic map, we constrain burial and exhumation across the Idaho-Montana fold-thrust belt (Figures 2 and 3) and evaluate a recent model that hypothesized that early thin-skinned deformation of the orogenic wedge was overprinted by later thick-skinned phases [25]. With systematic sampling across the study area, we empirically constrain the pattern of maximum burial temperatures along the transect and use the data to determine stratigraphic burial depths and the approximate geometry of the orogenic wedge. The results highlight the thermal signature of the enigmatic early phase of deformation (~145-90 Ma) in the Idaho-Montana fold-thrust belt and document dismemberment and reorganization of the orogenic wedge during a transition to thick-skinned thrusting in Late Cretaceous to Eocene time (~90-53 Ma). Using RSCM helps to elucidate the geometry of the early orogenic wedge that was later overprinted during progressive deformation. The results suggest that the geometry of the pre-thrusting strata within the Idaho-Montana fold-thrust belt set limits on the deformed wedge geometry; once a threshold amount of shortening within the wedge was surpassed, a deeper midcrustal décollement was activated, resulting in widely distributed sedimentation and subsequent erosion. These results have implications for global orogenic wedges, particularly where fold-thrust belts deform a sedimentary succession with variations in thickness and mechanical stratigraphy.

The retroarc continental fold-thrust belt of the Late Jurassic to early Eocene North American Cordillera spans much of western North America (Figure 2; e.g., [2, 26–32]). Like elsewhere in the Cordillera, the Idaho-Montana fold-thrust belt developed mostly within rift and passive margin strata flanking the western edge of Laurentia (e.g., [33, 34]). However, weak Neoproterozoic to Cambrian strata that host the basal décollement in much of the Sevier fold-thrust belt (e.g., [2, 35–40]) are largely absent along much of the Idaho-Montana fold-thrust belt, where the lower (Neoproterozoic to Early Ordovician) part of the rift-passive margin stratigraphic package pinches out against a major basement high in east-central Idaho termed the Lemhi arch (Figure 4; [41]). The Lemhi arch consists of a thin succession of Middle Ordovician or Devonian strata that unconformably overlies tilted fault blocks of Mesoproterozoic quartzites of the Belt Supergroup in east-central Idaho and Paleoproterozoic and Archean crystalline basement and metasedimentary rocks within the Dillon block of the Wyoming craton in southwestern Montana [42–45]. The unconformities and condensed sections within the Neoproterozoic to Cambrian rift-to-drift strata in the vicinity of the Lemhi arch (e.g., [42, 45–49]) make the Idaho-Montana fold-thrust belt a suitable place to test how inherited rift architecture and mechanical stratigraphy may influence the geometry and kinematics of the fold-thrust belt developed within those rocks.

In east-central Idaho, an intra-Devonian unconformity (Figure 4(a)) marks the lowest regionally continuous passive margin strata on top of the Lemhi arch [46]. Above this unconformity, uppermost Devonian to Pennsylvanian mixed siliciclastic and carbonate rocks are regionally extensive, defining a west-facing shelf to slope transition that drapes the Lemhi arch (e.g., [46, 51, 52]). These middle Paleozoic strata host a basal décollement to many thrusts of the thin-skinned sector of the Idaho-Montana fold-thrust belt (e.g., [25, 53, 54]), with a hingeline [50, 55] that mimics both the geometry of the former rift margin and the later Sevier fold-thrust belt (e.g., [2, 26, 52, 56, 57]). In spite of a lack of preservation of Permian to Cretaceous strata throughout much of the central portion of the Idaho-Montana fold-thrust belt, it is clear that both the laterally discontinuous lower package of Neoproterozoic to Devonian rocks and the thin but laterally continuous upper stratigraphic package of Devonian and younger rocks surrounding the Lemhi arch severely limited the availability of weak detachment horizons in the later fold-thrust belt [25, 54].

The Idaho-Montana fold-thrust belt is a double-decker thrust system with a shallower, thin-skinned region of thrusting and folding that overlies a structurally deeper, thick-skinned system [25]. Structurally above the Devonian Lemhi arch unconformity, the upper, thin-skinned portion of the Idaho-Montana fold-thrust belt is characterized by kilometer-wavelength detachment folds that are rooted in Mississippian and Devonian shales, siltstones, and evaporites [46, 58–65]. Detachment folds are commonly upright in the Lost River Range (e.g., [63]), but overturned to recumbent in the Lemhi Range and Beaverhead Mountains (e.g., [59, 61, 64–66]), signifying increasing shear strain toward the foreland along a décollement just above the Lemhi arch [25]. Major regional thrusts such as the Thompson Gulch and Medicine Lodge (Figure 3) utilized this décollement (e.g., [54, 67]), with a frontal ramp ultimately continuing to the surface in a localized foreland basin (e.g., [22, 68]).

Throughout east-central Idaho and southwestern Montana, the Devonian and Mississippian décollement of the upper thin-skinned thrust system was tilted and locally crosscut by younger, more deeply detached structures that are characterized by a thick-skinned structural style [54, 59, 61, 65, 69, 70]. For example, the Thompson Gulch and correlative thrusts were tilted and crosscut by the thick-skinned thrusts of the Poison Creek-Baby Joe Gulch-Hawley Creek thrust system [25, 67, 71]. Associated with the structurally deeper thrust system, the regional-scale Patterson culmination exhumed a portion of the Lemhi arch and tilted the overlying detachment folds in the Lemhi Range [59, 61]. While there is a recognizable mechanical stratigraphy across the region that relates to the Lemhi arch, subsequent deformation of the Lemhi arch resulted in a complex and irregular exhumation pattern that eroded or overprinted evidence of older deformation at higher structural levels.

Though the older phase of thin-skinned deformation above the Lemhi arch is scantily preserved due to later exhumation associated with the deeper thrust system, the magnitude of exhumation related to shortening was documented by Rodgers and Janecke’s [4] subcrop map of the Eocene Challis Volcanic Group, which represents a regional, post-orogenic unconformity. The subcrop map highlights large wavelength (~75 km) folds that plunge toward the southeast, exhuming Mesoproterozoic to Permian rocks in and adjacent to the Lemhi arch. To better constrain the early Cenozoic geometry of the Idaho-Montana fold-thrust belt, we applied paleothermometry in the context of this subcrop map to constrain the maximum temperatures recorded by strata in the upper, thin-skinned thrust system.

Retroarc fold-thrust belts are often successfully modeled in the context of critical taper theory (e.g., [72]). This simple and predictive model views active mountain belts as growing wedges of deforming material, using the analogy of an advancing snowplow that piles a forward-tapering wedge of snow by scraping along a fixed base. Horizontal shortening and vertical thickening accommodated by folding and/or thrust faulting on the crustal scale result in a wedge-shaped geometry of the fold-thrust belt; the upper surface of the orogenic wedge is defined by the surface slope of the growing wedge, and a basal slope is governed by the dip of the basal fold-thrust belt detachment. The wedge taper angle is the sum of these two angles. As horizontal shortening and vertical thickening occur, the surface slope and taper angle increase until they surpass a critical value set by the compressive strength of the material within the wedge and the frictional resistance to sliding at the base. Once this threshold is surpassed, a new “in-sequence” thrust fault (or fold) propagates forward, thereby reducing the surface slope and decreasing the taper angle as deformation advances toward the undeformed crust. The cycle continues as “out-of-sequence” thrusts thicken the interior of the wedge, once again increasing the taper angle and promoting wedge failure. In sandbox models, this alternating process of internal thickening and forward advance of the wedge of deformed material is governed purely by changes in the surface slope because the décollement is usually fixed at the base of the model.

In nature, however, the deepest available décollement usually only defines the base of the orogenic wedge in the more interior part of the fold-thrust belt; closer to the foreland, thrust slip is instead fed to structurally shallower décollements toward the toe of the orogenic wedge. As the wedge grows, these structurally shallower décollements become inactive as the base of the wedge becomes defined by progressively deeper décollements. As such, the mechanical stratigraphy within and beneath the growing frontal part of the orogenic wedge dictates where the down-stepping décollement may activate and, in turn, may limit the ability of the wedge to grow. The freedom to activate a new décollement at relatively deeper stratigraphic levels means that the stress state of the wedge is balanced not only by changing the surface slope but also by changing the basal slope. It has long been observed that thrust faults cut across stratigraphic layers of varying strength in a stair step fashion as they follow the path of least resistance [73]. Therefore, pre-thrusting stratal geometries and the distribution of weak rocks like shales and silty carbonates may determine the basal slope, which in turn impacts the surface slope (e.g., [56, 74, 75]). In other words, the mechanical stratigraphy of the crust determines the path of least resistance, thereby limiting the angle and position of the basal slope to the growing wedge.

In order to employ critical taper theory in the study of ancient orogenic belts and to investigate the potential role of mechanical stratigraphy in regulating wedge geometry and growth, it is therefore crucial to constrain both the surface and basal slope through time. The basal slope may be inferred from surface mapping and through the creation of balanced and restored cross-sections. Inferring the surface slope requires additional constraints, such as the thermal history and depths of exhumation of sedimentary rocks within the fold-thrust belt.

Where magmatic activity is absent, the maximum temperature recorded in a sedimentary rock can yield a maximum burial depth, assuming a geothermal gradient and surface temperature. Whatever the burial mechanism, by calculating maximum burial depths along a transect, a slope can be generated relative to the sampled datum (Figure 1). The rate of change of the maximum temperature trend should scale with the time-averaged surface slope of the wedge. Thus, constructing a maximum temperature profile across the Idaho-Montana fold-thrust belt allows for estimations of burial depths and of the former wedge geometry.

In this study, we infer the wedge geometry for the Idaho-Montana fold-thrust belt by using RSCM to obtain maximum temperatures along a regional transect (Figures 1 and 3). A balanced and restored regional cross-section provides context and constrains the geometry of the pre-orogenic stratigraphic package and décollement [25]. Forty-two samples of dark, fine-grained carbonate and siliciclastic rocks were collected for RSCM along a ~350 km long transect across the study area (Figure 3). Samples within an intact stratigraphic column in the Lost River Range and Beaverhead Mountains constrain a reasonable maximum paleogeothermal gradient, which is presented in the Discussion section. With these geothermal gradient estimates, we use calculated maximum temperatures as a proxy for maximum burial depths. When possible, samples were collected directly above the Devonian Lemhi arch unconformity (Lemhi arch unconformity of Grader et al. [46]) and near the top of the regionally extensive Mississippian part of the stratigraphic succession. By targeting these two stratigraphic datums, we can use regional stratigraphic thicknesses to make first-order inferences about the original basin geometry and the dimensions of the critically tapered orogenic wedge.

Standard 30 μm thick, microprobe-polished thin sections were analyzed at the University of Colorado Boulder using a LabRAM HR Evolution Raman microscope. For each sample, 18 spots of carbonaceous material were analyzed using a 532 nm (green) laser. Spots beneath the surface of the slide were selected, in order to avoid possible effects from polishing. For each spot, two analyses with 5-second acquisition times each were collected at 25% power and averaged. A 10-second acquisition time at 50% power was used for samples that had significant noise in the spectra due to high fluorescence.

Iterative curve fitting was performed using IFORS software [76], downloaded with the author’s permission from http://www.sediment.uni-goettingen.de/download/. This method is preferred because its automated and iterative modeling approach removes user bias [77] and does not assume a set number or position of defect peaks (ex., D1-4; [78]), which themselves vary with maximum temperature [18]. The method quantifies temperature-dependent changes to the shapes of the peaks in the D and G bands (Figure 5(a)): D band peaks get narrower and relatively taller, and G band peaks get narrower and relatively smaller with increasing maximum temperature from ~50 to 400°C (e.g., [18, 78]). Rather than using the heights and widths of the graphite (G) and D1-4 peaks to calculate temperature (e.g., [15]), this method uses the scaled total area of the G and D bands to calculate temperature [78]. Simply put, this method sums all intensity values in the respective D and G bands, for Raman shift values between 1000 and 1800 cm-1, and divides them by the maximum intensity value in the respective band (Dmax or Gmax). Examples of the curve fitting process are illustrated in Figures 5(b) and 5(c), with shaded areas under the respective curves highlighting the areas used in the calculations. Using a calibration curve [76], scaled total area values were then converted to maximum temperatures (°C). Two spectra (samples 1 and 2) were analyzed using a more widely adopted method [15] but yielded the same result as the IFORS method used here, suggesting that the automated IFORS methodology yields accurate and reproducible results [78].

To calculate the scaled total area, a simple linear regression function was used as the baseline, within a narrow window of 1000-1800 cm^-1 by taking the difference between measured and modeled intensity values and interpolating between them (Figure 5). This simplified approach avoids erroneous overinterpretation of noisy spectra and yields comparable results to the default settings [76, 78] for spectra with a high signal-to-noise ratio. While selecting a baseline, the spectra were deconvolved into their principal G and D peaks (Figures 5(b) and 5(c)). This was done iteratively using pseudo-Voigt functions to model the G and D peaks and adding the baseline function, until the composite spectra could be recreated with residuals ≤1%. This iterative curve-fitting process was performed three times for each spot.

Using an empirically derived calibration curve from Lünsdorf and Lünsdorf [76], the averaged scaled total area values were used to calculate the maximum temperature. Uncertainties typically range from ~36 to 40°C for temperatures calculated for individual spots [78]. It is important to note that the particular calibration curve used relies on a linear projection for values below ~140°C. Spectra were accepted if (1) at least two peaks were visible, (2) all peaks occurred entirely within the window (1000-1800 cm^-1), and (3) no unimodal peaks overprinted the D band. Individual calculated temperatures were only accepted if they were >0°C and not an outlier by >100°C.

Two contour maps were created, showing exhumation magnitude (Figure 6) and maximum temperature (Figure 7). To create the exhumation map, the lithology below the sub-Eocene unconformity was identified at 298 locations, including all locations of Rodgers and Janecke [4] and additional locations on published maps. Minimum exhumation magnitude (km) was estimated at each point using stratigraphic thicknesses reported on geologic maps and regional stratigraphic columns of pre-Cretaceous strata. A complete list of all locations, exhumation values, and references can be found in the supplementary materials (available here).

The maximum temperature contour map uses the results of this study and 149 compiled maximum temperature values, inferred from CAI, Rock-Eval pyrolysis, vitrinite reflectance, and mineral deformation mechanisms that constrain deformation temperatures (Table 1). Compiled data were projected onto the cross-section line (Figure 8) for samples that have the same structural position as RSCM samples of this study (see Table 1). Reported CAI values were converted to maximum temperature and associated uncertainty following the key of Epstein et al. [79], (Figure 5), which is calibrated by field and laboratory data. Maximum temperatures were calculated from vitrinite reflectance (R_o) values using the calibration curve for burial heating of Barker and Pawlewicz [80], (Figure 1). Maximum temperatures were calculated from Rock-Eval pyrolysis values using the calibration curve of Zdanaviciute and Lazauskiene [81], (Figure 5 in Epstein et al. [79]). Both contour maps were created in ArcMap 10.6, using the inverse distance weighted toolbox, and trimmed to the area of interest.

Flexural modeling was carried out in FlexMC, following methods reported in Saylor et al. [87, 88]. In cross-section, a triangular-shaped load, composed of up to 10 blocks, and its predicted flexural profile were modeled to fit the observed flexural profile constrained by stratigraphic thicknesses of the Kootenai and Blackleaf formations projected onto cross-section line $B‐B′$ (Figure 2). A broken plate model was used to simulate a convergent margin setting, with the load pinned at the back (west) end of the orogen perpendicular line, sloping toward the basin. The western limit of the load was defined at the western Idaho shear zone (near the inferred plate boundary during Mesozoic time), and the eastern edge of the load was varied based on the extent of potential source areas constrained by provenance studies [23, 24]. Default values were used for the basin infill density (2400 kg/m³), load density (2500-2700 kg/m³), mantle density (3300 kg/m³), and Young’s modulus (70 GPa) and Poisson’s ratio (0.25) of the plate. Allowable effective elastic thicknesses ranged from 3 to 81 km, following values reported in Lowry and Pérez-Gussinyé [89], and were permitted to vary exponentially along the cross-section line. For each modeled scenario, 1000 randomly generated loads were produced and compared to the observed flexural profile.

Two-dimensional flexural modeling of the thrust load was constrained using measured stratigraphic thicknesses of the Kootenai (ca. 135-110 Ma) and Blackleaf (ca. 110-90 Ma) formations from eight localities in southwestern Montana grouped into one proximal and one distal region (reported in [23, 24]). These units correspond to roughly the lower half of the Cretaceous strata shown in Figures 3 and 8. In the proximal group, subsidence magnitude is 465 m with a horizontal spatial uncertainty of 30 km to account for the geographic spread of the included measured sections. In the distal group, subsidence magnitude is approximated at 235 m with a horizontal spatial uncertainty of 50 km. Randomly generated loads and their predicted flexural subsidence profiles were forward modeled, with acceptable models matching the observed subsidence magnitudes constrained by the measured thicknesses. Two scenarios were modeled using a minimum (105 km) and maximum (205 km) load width extending from the pinned western end of the cross-section profile. For the minimum width scenario, a maximum load height of 8 km was required in order to produce model results. For the minimum width scenario, a maximum load height of 4 km was permitted. These results test what range of surface slopes for the early Cretaceous (pre ~100 Ma) Idaho-Montana fold-thrust belt—assuming the minimum and maximum allowable widths and heights—recreates the observed flexural profile of the preserved foreland basin.

Figure 6 shows a map of cumulative exhumation magnitude, modified and updated from Rodgers and Janecke [4]. This exhumation map shows that exhumation along the transect is highly irregular, with large wavelength (~75 km) highs and lows. High exhumation magnitudes (>4-5 km) reflect erosion (inferred to be Cretaceous in age) of the Lemhi arch and its thick outboard rift succession that constitutes the largest wavelength folds of the Idaho-Montana fold-thrust belt. An important result is that the areas with the highest exhumation magnitudes do not correspond with mapped thrusts at the surface as expected for a simple, wedge-shaped orogenic belt.

Thirty-two of the forty-two RSCM samples yielded acceptable maximum temperatures (Table 2; Figure 3). Figure 8(a) shows the results projected on a line parallel to and beyond the margins of a regional balanced cross-section that was restored for Cenozoic normal faulting [25]. The inset in Figure 8 shows the results in their approximate stratigraphic positions. The maximum temperature reported for each sample is a weighted mean of all accepted spot values (⁠$n$⁠), not including their individual 36-40°C (mean of 36.8°C) uncertainties. Propagated $2σ$ uncertainty values (internal and external) are reported for maximum temperature values of individual samples. Calculated maximum temperatures and uncertainties as well as measured Raman shift (cm^-1) and associated intensity values and curve-fitting model results for individual spots are reported in the supplementary material (available here).

The temperature results (Figure 9) delineate three internally consistent groups, which we describe from the foreland (northeast) to the hinterland (southwest). The first group contains the lowest temperature samples, from the Ruby Range and Blacktail and Tendoy Mountains of Montana, which give a weighted mean maximum temperature of $52±18$°C. Though this value is below the sensitivity of our calibration curve (~140°C), it is within the range measurable by RSCM (e.g., [11, 15, 18) and is comparable to maximum temperatures of $~65±15$°C at similar stratigraphic levels independently constrained using available CAI data [21]. Samples 8-10 have weighted mean temperatures of $52±61$°C (⁠$n=6$⁠), $42±62$°C (⁠$n=3$⁠), and $60±29$°C (⁠$n=16$⁠). Sample 11 gives a higher result of $115±27$°C (⁠$n=18$⁠). The samples nearest to the Devonian datum (8, 9, 11) overlap within uncertainty with the sample nearest to the Mississippian datum (10) (Figure 9), demonstrating that rocks of group 1 were all heated to comparably low temperatures.

A sharp increase in maximum temperatures along the transect occurs in the Tendoy Mountains (Figure 9), defining the boundary between sample groups 1 and 2. Group 2 samples show a uniform maximum temperature across a >160 km region, giving a weighted mean of $243±9$°C for the Mississippian datum and $272±9$°C for the Devonian datum. Between samples 10/11 and 12 (Figure 3), maximum temperatures increase by at least ~60°C and perhaps as much as ~180°C over a distance of <15 km at the same stratigraphic horizon but across the Johnson thrust (Figure 8). Two Mississippian samples from the Tendoy Mountains (12 and 13) have weighted mean maximum temperatures of $231±26$°C (⁠$n=17$⁠) and $216±29$°C (⁠$n=18$⁠), respectively. Four samples from the central Beaverhead Mountains near the Idaho-Montana border were collected in the context of recently published 1 : 24,000 scale mapping [90]. These samples (7, 3, 5, and 6) have weighted mean maximum temperatures of $281±25$°C (⁠$n=17$⁠), $278±26$°C (⁠$n=18$⁠), $267±27$°C (n =18), and $244±25$°C (⁠$n=18$⁠), respectively. These samples span from precisely at the Devonian datum to hundreds of meters above the Mississippian datum, suggesting a general up section decrease in temperature, although the calculated temperatures overlap within uncertainty, giving a weighted mean of $267±13$°C. This interpretation is consistent with a local CAI value of 4 [21] that suggests a maximum temperature of $~245±55$°C nearest the Mississippian datum.

Five samples from the Lemhi Range and Donkey Hills of east-central Idaho (Figure 3) follow the Devonian $datum>30$ km across the regional-scale anticlinorium of the Patterson culmination (Figure 8). These samples (15, 14, 16, 17, and 18) have weighted mean maximum temperatures of $260±28$°C (⁠$n=17$⁠), $267±27$°C (⁠$n=15$⁠), $275±25$°C (⁠$n=18$⁠), $290±30$°C (⁠$n=16$⁠), and $244±30$°C (⁠$n=18$⁠), respectively. All of these values overlap within uncertainty, giving a weighted mean of $267±12$°C. The final eleven samples from group 2 come from the folded Devonian and Mississippian sections in the Lost River Range and White Knob Mountains of Idaho. These samples (19, 1, 2, 28, 20, 35, 34, 32, 33, 4, and 30) have weighted mean maximum temperatures of $229±31$°C (⁠$n=19$⁠), $259±26$°C (⁠$n=16$⁠), $280±26$°C (⁠$n=17$⁠), $273±26$°C (⁠$n=17$⁠), $263±24$°C (⁠$n=18$⁠), $263±28$°C (⁠$n=18$⁠), $239±29$°C (⁠$n=18$⁠), $258±27$°C (⁠$n=18$⁠), $241±25$°C (⁠$n=18$⁠), $218±26$°C (⁠$n=18$⁠), and $221±26$°C (⁠$n=18$⁠), respectively. Of the samples with tight stratigraphic control, two samples near the Devonian datum overlap and give a weighted mean of $277±18$°C, whereas seven samples near the Mississippian datum overlap and give a weighted mean of $244±10$°C. These values agree with a local CAI of 5 [20] that suggests maximum temperatures of $~300±50$°C. Temperatures calculated from illite crystallinity and deformation lamellae in calcite also suggest temperatures ranging from ~200 to 300°C [82].

The third group contains six samples from the Pioneer and Boulder Mountains of central Idaho (not shown in Figure 8) that give a weighted mean of $303±14$°C, when two outliers are excluded. This value is consistent with fifteen CAI samples with values of 5 that suggest maximum temperatures of $~300±50$°C [20]. The stratigraphic context is not as robust for most of these samples, which may explain some of the differences in groups 2 and 3. Four samples from group 3 (29, 23, 25, and 26) have weighted mean maximum temperatures of $292±25°C$ (⁠$n=18$⁠), $284±31$°C (⁠$n=14$⁠), $298±29$°C (⁠$n=18$⁠), and $339±31$°C (⁠$n=17$⁠), respectively. These samples come from Devonian and Pennsylvanian-Permian formations, respectively. Two additional outliers (22 and 27) are mostly distinct from this subpopulation with weighted mean maximum temperatures of $390±37$°C (⁠$n=18$⁠) and $392±52$°C (⁠$n=4$⁠). Both of these outliers are within ~500 m of mapped intrusive bodies.

Figure 7 shows a contour map of maximum temperatures of all available temperature data, from multiple independent paleothermometers (RSCM, CAI, and R_o; see Table 1 for compiled data sources). This maximum temperature map bolsters the trend that is apparent in our results (Figure 9). Sampled Paleozoic rocks delineate a relatively cold foreland (<100°C) and a broad, relatively warm (~250-300°C) hinterland, corresponding with groups 1 and 2 (Figure 9), respectively. These well-defined thermal domains are separated by a steep gradient near the Idaho/Montana border. A hotter (>300°C) portion of the hinterland is apparent but less well-defined, corresponding with group 3 (Figure 9).

Results of 2-D flexural modeling are consistent with the observed thicknesses of the Kootenai and Blackleaf formations. Our preferred model utilizes a broken (rather than infinite) plate because of the inferred proximity of the foreland basin to the Western Idaho shear zone and, more generally, the convergent plate boundary. However, we describe results from the infinite plate model to demonstrate that the choice of a broken vs. an infinite plate does not significantly affect the results. Results from our preferred broken plate models suggest orogenic wedge surface slopes ranging from ~1 to 4°. For the minimum width/maximum height load scenario using a broken plate model, two of 1000 model iterations (0.2%) yielded reasonable fits and a mean surface slope of $4.0±0.7°$ (Figure 10(a)); for an infinite plate, one of 1000 model iterations (0.1%) yielded a reasonable fit with a mean surface slope of 5.1°. However, we consider both of these minimum width/maximum height model results unrealistic because they require unreasonably thick loads with maximum load heights of $6.3±0.7$ km (broken plate) and 6.6 km (infinite plate) and mean topographic elevations across loads of $~2.9±0.5$ km (broken plate) and ~3.5 km (infinite plate) (e.g., [91]). In contrast, the maximum width/minimum height load scenario (Figure 10(b)) better fits the data; of 1000 models for a broken plate, seven (0.7%) were a reasonable fit, yielding a mean surface slope of $1.1±0.3°$⁠; for an infinite plate, 8 of 1000 models (0.8%) were a reasonable fit, yielding a mean surface slope of $1.3±0.4°$⁠. The maximum height of modeled loads for the broken plate scenario is $3.1±0.6$ km with a mean topographic elevation across the load of $1.4±0.3$ km; for an infinite plate, the maximum height of modeled loads is $3.5±0.5$ km with a mean topographic elevation across the load of $1.7±0.2$ km. These results suggest that the maximum width/minimum height load scenario is more plausible and that the surface slope at ca. 100 Ma was closer to 1° than 4°.

Using data from multiple independent paleothermometers (RSCM, CAI, and R_o), we document an abrupt increase in maximum temperature across the Johnson thrust (thrust D of DuBois [92]) and uniform maximum temperatures within much of the Idaho-Montana fold-thrust belt (Figures 7 and 9); this suggests that significant and approximately uniform burial occurred across much of the Idaho-Montana fold-thrust belt, with an abrupt front occurring across the Johnson thrust in the Tendoy Mountains. The RSCM results replicate results from CAI and R_o data, suggesting that the recorded temperatures are representative. Available low-temperature thermochronometric data (zircon (U-Th)/He, closure temperatures ~180°C; [93]) are generally reset for samples southwest of the Tendoy Mountains (e.g., [94–96]), but not for samples within the Montana foreland to the northeast (e.g., [97–99]). This regional trend and the internal consistency of the RSCM and CAI data suggest that the calculated maximum temperatures are not only accurate but also of tectonic origin.

There are several lines of evidence that suggest that the observed regional temperature trend is not an artifact of local temperature discrepancies and is instead a result of regional burial. Fluid flow and deformation likely had a negligible direct impact on our results because our sampling strategy avoided highly veined or deformed units and samples of different permeabilities yielded comparable maximum temperatures that systematically increase with stratigraphic depth (Figure 8). For instance, four samples (3, 5, 6, and 7) of limestone and dolostone mudstones and wackestones have a wide range of permeability, yet they yielded maximum temperature estimates that systematically increase with stratigraphic depth and overlap within uncertainty. Though gravity-driven fluid flow was likely active within the fold-thrust belt, this likely occurred through permeable units not sampled in this study [82, 100]. In the study area, fluid flow was apparently heterogeneous on the kilometer and smaller scale and closely associated with deformation [82], unlike our results which are remarkably uniform over the kilometer scale. Local magmatism likely affected only the two highest temperature samples in the distal part of the Idaho batholith (22 and 27), which are outliers at ~390°C. Magmatic heating has been shown to affect RSCM samples over fairly short distances of ~0.5 km (e.g., [101]). For these reasons, we argue that the regional trend of elevated temperatures between ~250°C and 300°C for carbonates along the transect is robust.

In order to infer burial depths from calculated maximum temperatures, the associated geothermal gradients must be known. Foreland basins typically have low geothermal gradients (~22°C/km; [102]). However, rather than assume a geothermal gradient to calculate burial depth, we can independently constrain maximum geothermal gradients by comparing maximum temperatures over a known vertical distance within well-constrained localities along our transect. Samples from the central Beaverhead Mountains and Lost River Range cover ≥1 km of section over short map distances (Figure 8), thereby offering some constraints regarding the maximum permissible geothermal gradient. In the central Beaverhead Mountains, temperatures range from 244°C (±25) to 281°C (±27) over a stratigraphic distance of ~1 km. Two samples from a vertical transect on Lost River Peak have temperatures of 259°C and 280°C (±26) over a stratigraphic distance of ~750 m. These observations are consistent with a maximum geothermal gradient of ~30-40°C/km. Although the uncertainty of the method precludes us from precisely constraining the maximum geothermal gradient, the fact that most Mississippian and Devonian samples overlap within uncertainty (Figure 9) suggests that a maximum geothermal gradient of ~30-40°C/km is reasonable. The upper bound of this estimate range overlaps with maximum geothermal gradients of ~40-60°C/km estimated for the Sevier hinterland in Nevada [103, 104] and Wyoming foreland [100] during Cretaceous time. For comparison, the region of active thrusting within the central Wyoming salient records a low Cretaceous geothermal gradient (~20°C/km; [100]).

Modern geothermal gradients in fold-thrust belts commonly vary from the hinterland to the foreland (Figure 1 of Husson and Moretti [105]), making it difficult to choose a unique geothermal gradient across our entire transect for calculating burial depth. However, in the Andes, modern geothermal gradients vary by less than 10°C/km between the hinterland and the foreland, over 250 km apart [105]. Geothermal gradients in fold-thrust belts may also show more localized along-strike variations. For example, in the Zagros fold-thrust belt, modern geothermal gradients show along-strike variation of 15-20°C/km over 50-200 km distances [106, 107]. Because independent estimates of paleogeothermal gradients across our transect are not available, we use a relatively high geothermal gradient to provide a conservative constraint on the minimum magnitude of burial and exhumation. We note that a variation of ~10°C/km across our transect would increase our burial estimates by as much as a third.

In order to make conservative maximum burial estimates (Figure 9(b)), we assume a high but reasonable geothermal gradient of ~40°C/km and a mean surface temperature of ~10°C (e.g., [108]). Although the calculated values for group 1 samples are limited by the sensitivity of our calibration curve and the maximum paleogeothermal gradient was likely <40°C/km, they could not have been buried significantly deeper than the ~1.5 km required by the stratigraphic thickness without being heated above $52°C±18$⁠. For group 2 samples, the mean maximum temperatures are regionally consistent for all the Mississippian (⁠$243°C±9$⁠) and Devonian (⁠$272°C±9$⁠) samples, respectively, giving average burial estimates of 5.9 km (±0.4) and 6.7 km (±0.3). Overall, burial estimates for group 2 samples range from 5.2 km (±0.7) to 7.0 km (±0.8), much greater than the ~2.5-4.5 km preserved stratigraphic thicknesses on top of the Lemhi arch (Figure 8; [25]). In order to constrain the discrepancy between estimated burial depths and predicted maximum temperatures based on stratigraphic thickness, we first grouped data points with overlapping temperature values from samples with similar stratigraphic thicknesses. For each cluster, the magnitude of the “burial deficit” was calculated by subtracting the pre-Cretaceous stratigraphic thickness (Figure 8(b)) from the mean maximum burial estimate. The results are shown in Figure 8(c). The burial deficit is absent east of the Tendoy Mountains and is most pronounced in the central Beaverhead Mountains (Figure 8), where burial estimates of at least 5.9 km (±0.6) for Pennsylvanian strata and 6.8 km (±0.6) for Devonian strata exceed the observed stratigraphic thickness by 4.9 km (±0.7). Even in the thicker section southwest of the Lemhi arch, burial estimates exceed stratigraphic thickness by 1.4 to 3.6 (±0.7) kilometers (Figure 9(c)). An inferred maximum depth of 8.2 km (±0.8) for Pennsylvanian/Permian strata (R26) in the Boulder Mountains, in the westernmost portion of our study area, overlaps with emplacement depths estimated for the Cretaceous Idaho batholith where it was intruded at a comparable stratigraphic level [109, 110], suggesting that our maximum depth estimates are both accurate and Cretaceous in age.

Stratigraphic uncertainties, major unconformities related to Paleozoic exhumation, and a lack of rocks younger than Permian make the magnitude of the burial deficit poorly constrained in the western part of the study area. Given this uncertainty, it is possible that samples in the Pioneer and Boulder Mountains have calculated temperatures that may be explained by stratigraphic burial alone, if basin development and exhumation during Antler and Ancestral Rocky Mountain orogenesis were localized [111–115]. Regardless of this uncertainty, the fold-thrust belt above the Lemhi arch was clearly buried ~1.5-5.0 km deeper than accounted for by pre-thrusting stratigraphic burial, meaning that the samples were likely buried by Cretaceous strata and/or thrusts that have since been eroded.

Along the transect, calculated maximum temperatures are consistent (Figure 9) and do not correlate with the heterogeneous pattern of cumulative pre-Eocene exhumation (Figure 6; [4]), suggesting that they reflect a regional trend, not local perturbations in heat flow. This suggests that maximum temperatures acquired during burial predate late-stage folding and faulting associated with the deeper, thick-skinned thrust system along the transect. For example, the presence of carbonate mylonite in the thin-skinned Thompson Gulch thrust, which is crosscut by the thick-skinned Baby Joe Gulch thrust, requires that the host rocks were hotter than the threshold for carbonate plasticity (>~250°C) during slip on the Thompson Gulch thrust and before slip along the Baby Joe Gulch thrust [25]. Like for the Baby Joe Gulch thrust, maximum temperatures are shifted across the Johnson thrust, demonstrating that maximum burial also predated this thick-skinned thrust (Figures 8 and 9). Similarly, Devonian samples above the Patterson culmination give overlapping temperature estimates with a mean of $267°C±15$⁠, despite over 3 km of structural relief (Figure 8), demonstrating they were achieved before deeply detached folding and thrusting in this region. The regionally consistent maximum temperatures, in all but the samples from the hanging wall of the McKenzie and Four Eyes Canyon thrusts, require consistent burial and subsequent exhumation of rocks in east-central Idaho (Figure 8). Therefore, rocks must have attained their maximum temperatures during the early phase of thin-skinned shortening and before late-stage folding, faulting, uplift, and exhumation that was spatially coincident with the Lemhi arch.

The absolute age of burial heating during the early phase of shortening can be constrained using available maximum depositional ages in adjacent foreland basin strata and low-temperature thermochronometric data from the hanging walls of regional thrusts. Crustal thickening in the arc and forearc from ca. 140 to 90 Ma [116–118] and detrital zircons in the foreland matching Cambrian to Mississippian sources from west of the Lemhi arch [23, 24, 119] demonstrate that the Idaho-Montana fold-thrust belt was already well-developed by mid-Early Cretaceous time. Specifically, maximum depositional ages and provenance studies for the Kootenai Formation, the oldest definitive foredeep deposit in the study area, demonstrate that shortening and exhumation down to lower Paleozoic levels occurred in the adjacent fold-thrust belt by ca. 145-125 Ma [23, 120]. Zircon (U-Th)/He thermochronometric data suggest that Belt basin quartzites beneath the Lemhi arch unconformity were cooled below ~180°C at ~87 Ma within the Patterson culmination [94] and ~68 Ma for the Poison Creek [95] and Hawley Creek thrusts [96], bracketing the youngest possible age of burial heating. This age range agrees with provenance studies that suggest that strata of the Lemhi subbasin, which make up the Lemhi arch, were not a sediment source to the foreland basin until during the deposition of the Frontier Formation at ~90-85 Ma [121]. Further, ~500 Ma zircons that are unique to the Beaverhead plutonic belt that is exposed in the hanging wall of the Baby Joe Gulch-Hawley Creek thrust system first appeared in the Beaverhead Group at ~67 Ma [22]. Therefore, maximum burial must have occurred prior to these exhumation ages. It is possible that burial was progressive, but the best available constraints limit the timing of maximum burial above the Lemhi arch to after ~145 and prior to ~90-70 Ma.

Restoring the late-stage faults and folds that disrupt the observed maximum temperature pattern allows us to describe the geometry of the early-stage burial envelope and make reasonable hypotheses regarding the geometry of the eroded orogenic wedge. Restoring the Johnson thrust (Thrust D of DuBois [92]) so that the Lemhi arch unconformity is continuous results in >4 km of burial difference between the samples, enough to account for the observed ~160°C temperature difference between groups 1 and 2, assuming a reasonable geothermal gradient of ~40°C/km. In other words, the sharp maximum temperature difference of 121°C (±57) over a distance of <15 km is likely a product of late-stage faulting on the Johnson thrust that juxtaposed different structural levels. The geometry during early deformation can be visualized by following the red asterisks in Figure 7, shown along the Medicine Lodge-McKenzie thrust system, which forms a continuous décollement when faults and folds of the Lemhi arch are restored. The initially continuous wedge of the thin-skinned Medicine Lodge-McKenzie thrust system was folded, faulted, uplifted, and exhumed by more deeply detached thick-skinned thrusts, such as the Hawley Creek and Johnson thrusts. This interpretation (Figure 11) is consistent with correlative hanging wall and footwall strata in the Medicine Lodge and McKenzie (and Four Eyes Canyon) thrusts (e.g., [53, 54]). Samples in the central Beaverhead Mountains come from the same structural position, within the footwall of the Thompson Gulch thrust, which links the flat décollement near the intra-Devonian unconformity in the Lemhi Range to the footwall ramp of the Medicine Lodge thrust [25]. Samples farther to the southwest, from the Lemhi and Lost River ranges, come from just above and below this flat regional décollement (White Knob fold belt of Beutner and Hait [59, 61]). Overlaying flat-lying isotherms and a 40°C/km geothermal gradient on this geometry recreates the observed temperatures (Figure 11(b)). In summary, restoration of late-stage thrust faults and folds highlights the décollement geometry of the Medicine Lodge-McKenzie thrust system, with a regional thrust flat near the intra-Devonian unconformity of the Lemhi arch at ~6.5 km depth and a frontal ramp between the Beaverhead and Tendoy Mountains.

Consistent burial of the Lemhi arch across the region during the early phase of deformation (prior to ~90-70 Ma) may be explained by several mechanisms. A now completely eroded (“phantom”; [122]) Early Cretaceous foredeep basin is one possible hypothesis, although strata of this age are not preserved and therefore had an unknown thickness. However, the observed temperatures and thicknesses require that the hypothesized basin would have increased in thickness northeastward, from ~1.4 km above the Lost River Range to ~4.9 km above the Beaverhead Mountains (Figure 9(c)). This apparent basin geometry is inconsistent with predictions of a phantom foredeep [122] because accommodation in foredeep basins is expected to reach a maximum value adjacent to the thrust load (e.g., [123]). Flexural modeling results for the maximum width load scenario show what a potential phantom foredeep may have looked like at ca. 100 Ma, prior to the inception of thick-skinned shortening and widespread exhumation along the Idaho-Montana fold-thrust belt (Figure 8). Whereas this model is consistent with preserved stratigraphic thicknesses in foreland basin strata of southwestern Montana, the predicted burial magnitudes in the more inboard part of the fold-thrust belt fall ~1-4 km short of explaining the observed temperatures. Because orogenic wedges grow progressively larger through time, the amplitude of the flexural deflection is expected to increase through time due to the growing load (e.g., [123]). However, major thrust faults between the Beaverhead and Pioneer Mountains are generally absent (Figure 3), and thus, there is nowhere to place the load adjacent to the hypothesized phantom foredeep. Further, because the maximum temperatures are interpreted to predate thick-skinned deformation in the study area and the modeled flexural wavelength is on the 100 km scale, samples adjacent to frontal thin-skinned thrusts 10s of km away (southwestern group 1 samples) should have also been deeply buried by the same foreland basin strata, yet they yield consistently low maximum temperatures (Figure 8). This suggests that a now-eroded foredeep basin cannot explain the calculated burial deficits.

Distributed folding and faulting at high structural levels (Pennsylvanian or younger) is an alternative mechanism for fold-thrust belt burial. Evidence of crustal thickening and erosion at high structural levels is suggested by a recent provenance analysis of the upper Blackleaf Formation, which suggests that upper Paleozoic and Triassic sources were eroded from central Idaho at ca. 110-100 Ma [119]. This would require doubling the thickness of the section by thrust burial in the footwalls of the Thompson Gulch and Medicine Lodge thrusts, which we hypothesize represented the décollement of the early fold-thrust belt [25]. However, exhumation during thrusting would make it unlikely that this entire thickness would have been preserved; therefore, this scenario alone is also unlikely to explain the burial deficits along the transect.

An alternative burial mechanism is a wedge-top basin. Wedge-top basins fill accommodation space next to the frontal regions of fold-thrust belts and may result in significant burial. For instance, the thickest part of the latest Cretaceous-early Cenozoic foreland basin succession in the study area occurs in Beaverhead Group strata that are interpreted to have been deposited in a wedge-top depozone [22]. Though these Beaverhead Group strata were deposited after the interval considered here, an older wedge-top succession was likely present within the study area during the time in question but is no longer preserved. Thus, a wedge-top (or hinterland) basin that buried a low-relief fold-thrust belt is a more permissible hypothesis that is consistent with both the wide and nearly uniform shape of the burial envelope (Figure 9(c)) and the necessary lack of exhumation on uplifted hanging walls. This preferred interpretation—of a low-relief orogenic wedge largely buried by wedge-top deposition—is schematically illustrated in Figure 11(b) (see also Fig. 1(b) of Simpson [124]). Similar, broad, low-relief orogenic wedges dominated by wedge-top sedimentation were produced by numerical models by Simpson [124] with very weak décollements (e.g., salt). Simpson [124] also documented an empirical relationship, where decreased surface slopes are correlated with proportionally more deposition in the wedge-top as compared to the foreland basin (Fig. 11(b) of [124]). This identifies the Jura, Apennines, Zagros, and Salt Range as potential analogs to the Idaho-Montana fold-thrust belt at ca. 100 Ma (Figure 11(b)), where exceptionally low surface slopes (<1.5°) may have promoted widespread wedge-top deposition. Flexural modeling results for the maximum width load scenario illustrate that a wedge with a surface slope of $1.0±0.3°$⁠, which overlaps with the surface slope constrained by our RSCM data, is compatible not only with preserved stratigraphic thicknesses but also with the dimensions of the required burial deficit.

Recently discovered remnants of an extensive wedge-top basin within the interior of the Canadian Cordillera [125] provide another potential analog. Cretaceous strata within this preserved wedge-top basin buried the underlying rocks, heating them by 115-135°C [125], an amount that could sufficiently explain the burial deficit for 75% of our samples. McMechan et al. [125] similarly concluded that an extensive wedge-top basin overlapped much of the active fold-thrust belt and that the corresponding orogenic wedge had a very low taper, indicative of a very weak décollement. Like the adjacent Canadian Cordillera, structural, stratigraphic, and temperature data from the Idaho-Montana fold-thrust belt are consistent with a “phantom” wedge-top basin that buried much of the low-relief orogenic wedge prior to ca. 90 Ma.

A combination of burial by folding, thrusting, and wedge-top deposition is a simple explanation for the observed burial deficits in the Idaho-Montana fold-thrust belt. As discussed previously, an in-sequence transition from a supercritical wedge and paired foredeep to a subcritical wedge and a “phantom” foredeep fails to explain our data for several reasons. Numerous crosscutting relationships between thin- and thick-skinned thrusts require a more sophisticated model with out-of-sequence deformation and a migrating décollement [25], which we propose created a “phantom” wedge-top basin. Relative timing constraints suggest that this burial predated Late Cretaceous thick-skinned thrusting. Thus, burial is thought to have occurred during earlier thin-skinned thrusting above the Lemhi arch unconformity, which had a shallow basal décollement; this, coupled with results from flexural modeling, demonstrates that the mid-Cretaceous orogenic wedge of the Idaho-Montana fold-thrust belt likely had a very low surface slope/wedge taper above the low-relief Lemhi arch basement high. Because the thickness of the wedge apparently mimicked the low-relief unconformity that draped the Lemhi arch and acted as the décollement (Figure 11), this suggests a causal relationship between the pre-thrusting mechanical stratigraphy and the resulting low-taper fold-thrust belt. In order to thicken the wedge and increase the taper angle toward its critical taper value required for further growth of the orogenic wedge, a step-down of the basal décollement into the Lemhi arch provided an alternative way to increase wedge taper (Figure 11(c)). In other words, deformation above the flat décollement following the top of the Lemhi arch apparently was not capable of increasing the taper angle for self-similar propagation of the wedge.

This model is appealing from a regional standpoint because it identifies the structure that accommodated the undocumented pre-90 Ma shortening in the Idaho-Montana fold-thrust belt: a décollement that followed the intra-Devonian Lemhi arch unconformity (i.e., Medicine Lodge-McKenzie thrust system). This model also suggests that thin strata overlying a basement high limited the dimensions of the wedge by not only determining the path of the basal décollement but also setting the threshold at which the stress state of the wedge could no longer be balanced without establishing a new décollement in the mechanical basement. In other words, thin-skinned thrusting in the thin passive margin strata may have accommodated its maximum amount of shortening before the décollement stepped down and began accommodating shortening in a more thick-skinned style.

An implication of this study is that the availability of potential décollements, in the form of weak sedimentary cover rocks, may fundamentally limit the ability of an orogenic wedge to grow in two related ways: (1) by determining the maximum allowable surface slope that can be built above a particular décollement and (2) by determining when and where the basal décollement must migrate in order to increase the taper angle by increasing the basal slope. Furthermore, low surface slopes may result in more distributed sedimentation that extends far into the fold-thrust belt, blurring the lines between wedge-top and foredeep depozones [124]. Interestingly, the Jura, Apennines, Zagros, and Salt Range (Figure 11 of [124]) are potential analogs to the Idaho-Montana fold-thrust belt not only in terms of low surface slopes and widespread wedge-top deposition but also in terms of structural style. Early thin-skinned and late thick-skinned thrusting has been proposed for all these potential analogs (e.g., [126–130]), remarkably similar to the double-decker model proposed by Parker and Pearson [25]. Unlike conventional fold-thrust belts, where shortening promotes exhumation above active thrusts and localized deposition in the foreland throughout orogenesis, low-relief wedges may promote widespread deposition, burial heating, and weakening of the crust; this leads to thick-skinned shortening of the underlying mechanical basement and widespread exhumation of the former thin-skinned fold-thrust belt.

The transition we describe (Figure 11) from a low-relief wedge detached in the upper crust that buries a basement high, to a higher-relief wedge detached in the middle crust that exhumes the basement high and former orogenic wedge, may initiate when thresholds are surpassed related to stratigraphic, rheological, or thermal properties [25, 126, 130–132]. This transition may occur once shortening can no longer be accommodated in the sedimentary cover rocks and strain propagates into the crystalline basement, as a viscous detachment in the middle crust is adopted (e.g., [25, 32, 133, 134]). Alternatively, the transition may occur once the foreland-advancing strain front first encounters viscous middle crust, where the continental crust is sufficiently thick (e.g., [130, 132]). The transition could also be thermally triggered, once increasing or decreasing temperatures alter the rheology of the lithosphere (e.g., [126]). These three end-member models, referred to as stratigraphic, rheological, and thermal inheritance by Horton and Folguera [131], may describe connected feedbacks within orogenic systems. Given the coincidence between the position of the Lemhi arch and the thermal gradient documented in this study, crustal rheology and stratigraphic thickness, inherited from the rifted continental margin, may predict where low-relief orogenic wedges form in the first place. However, the remarkable uniformity of the maximum temperatures along the low-relief Lemhi arch suggests a thermal trigger, with the transition occurring once the carbonate-hosted décollement exceeds ~270°C. Distributed deformation may lead to distributed sedimentation and widespread burial heating, facilitating the gradual transition from a thin-skinned fold-thrust belt detached in brittle sedimentary cover rocks to a thick-skinned belt detached in the viscous middle crust.

This study demonstrates not only that RSCM and other paleothermometers are useful tools for documenting and investigating switches from early burial to late exhumation phases in orogenic wedges but also that the thicknesses of sedimentary cover rocks, including “phantom” basins, may be a common attribute among many recently proposed tectonic inheritance models [25, 126, 130–132]. Thin preexisting strata may control the growth of an orogenic wedge by essentially dampening the ability of the surface slope to increase, resulting in a more diffuse flexural subsidence profile, widespread deposition across the wedge-top depozone, and burial heating throughout the fold-thrust belt. Both the inability of the surface slope to increase and widespread burial heating may work together to activate a deeper basal décollement in the ductile middle crust, in order to balance the stress. A balanced and restored cross-section of the Idaho-Montana fold-thrust belt will provide testable hypotheses regarding threshold values for transitional orogenic wedges. Continued investigations using paleothermometry in the shallow crust may be an effective way to investigate how the mere presence or absence of sedimentary cover rocks may set this threshold, setting the stage for the behavior of future growing orogenic wedges.

Raman spectroscopy of carbonaceous material (RSCM) of samples collected over a more than 350 km wide transect of the Idaho-Montana fold-thrust belt yields a low-relief maximum temperature profile (at ca. 145-70 Ma) averaging 272°C (± 9), over the ~160 km wide Lemhi arch basement high. Elevated temperatures are higher than can be explained by stratigraphic burial alone. Observed temperatures constrain a maximum geothermal gradient of ~40°C/km, suggesting that at least an additional 1.5-5.0 kilometers of burial must have occurred between ~145 and 90-70 Ma. We hypothesize that burial resulted from a combination of distributed folding and thrust faulting beneath a phantom (now eroded) wedge-top basin. The low-relief orogenic wedge had a surface slope that mimicked the geometry of the Lemhi arch basement high that hosted the décollement. After being buried to ~6.5 km, the décollement was abandoned for a deeper one, resulting in the observed cumulative exhumation patterns. Folding, faulting, and exhumation of the early low-relief wedge and the underlying Lemhi arch led to the disparity between the variable cumulative exhumation magnitudes inferred from the sub-Eocene unconformity and the consistent maximum temperatures recorded in the rocks above the Lemhi arch.

Our results suggest that thicknesses of pre-thrusting strata may limit the maximum possible thicknesses and surface slopes of orogenic wedges. Low-relief orogenic wedges formed in thin stratigraphic cover rocks may be subject to self-organized reorganization as the décollement is forced to migrate from weak horizons in the sedimentary cover rocks into the underlying mechanical basement in order to increase the taper angle of the system toward critical. Burial heating may facilitate this transition by weakening the middle crust. The Idaho-Montana fold-thrust belt illustrates how low-angle surface slopes result in more distributed flexural subsidence profiles, promoting widespread wedge-top deposition and burial heating. In this way, the thicknesses of weak sedimentary rocks in the upper crust may set thresholds that govern the growth of orogenic wedges, affecting not only their geometries but also their temperatures.

Data supporting the results of this study can be found in the manuscript text, the supplemental materials file, and publications referenced in the text.

This study represents a portion of SP’s PhD dissertation.

The authors declare that they have no conflicts of interest.

We thank Chase Porter for collecting some samples in the western part of the transect as well as Eric Elliston, Ryan Anderson, Alexis Ault, Kelly Bradbury, Tyler Paladino, and Keno Lünsdorf for their assistance during data collection and processing. We also thank Brian Horton for a constructive review. Funding for this research was provided by the Geological Society of America graduate student research grant, Tobacco Root Geological Society field scholarship, Idaho State University Center for Ecological Research and Education grant to SP, and NSF-Tectonics grants to DP [EAR-1728563] and EF [EAR-1727504].