Classically held mechanisms for removing mountain topography (e.g., erosion and gravitational collapse) require 10-100 Myr or more to completely remove tectonically generated relief. Here, we propose that mountain ranges can be completely and rapidly (<2 Myr) removed by a migrating hotspot. In western North America, multiple mountain ranges, including the Teton Range, terminate at the boundary with the relatively low relief track of the Yellowstone hotspot. This abrupt transition leads to a previously untested hypothesis that preexisting mountainous topography along the track has been erased. We integrate thermochronologic data collected from the footwall of the Teton fault with flexural-kinematic modeling and length-displacement scaling to show that the paleo-Teton fault and associated Teton Range was much longer (min. original length 190-210 km) than the present topographic expression of the range front (~65 km) and extended across the modern-day Yellowstone hotspot track. These analyses also indicate that the majority of fault displacement (min. 11.4-12.6 km) and the associated footwall mountain range growth had accumulated prior to Yellowstone encroachment at ~2 Ma, leading us to interpret that eastward migration of the Yellowstone hotspot relative to stable North America led to removal of the paleo-Teton mountain topography via posteruptive collapse of the range following multiple supercaldera (VEI 8) eruptions from 2.0 Ma to 600 ka and/or an isostatic collapse response, similar to ranges north of the Snake River plain. While this extremely rapid removal of mountain ranges and adjoining basins is probably relatively infrequent in the geologic record, it has important implications for continental physiography and topography over very short time spans.

Studies of mountain ranges commonly invoke erosion and extensional collapse to explain the reduction of topographic relief [15] and rarely do such studies consider less common, but geologically significant mechanisms such as supercalderas. In regions that have experienced a prolonged history of explosive volcanism and particularly regions impacted by supercalderas eruptions, it is useful to consider whether or not these cataclysmic mechanisms could denude or even completely diminish mountain topography where it intersects the caldera boundaries. In the Basin and Range province in the western US, the ranges generally represent the uplifted footwalls of crustal-scale normal faults [6]. In the northern Basin and Range of Idaho and Wyoming, the NNW striking crustal-scale normal faults and their associated uplifted footwall mountain blocks and adjoining half-grabens terminate where they intersect the anomalous low relief of the Snake River Plain (SRP; Figure 1). The SRP marks the ENE migration of the Yellowstone hotspot (YHS), and the truncation of Basin and Range structures on either side of the SRP leads to the hypothesis that the YHS may have somehow removed the associated mountain topography [710]. The Teton, Gallatin, Madison, and Centennial ranges all abut and appear to be truncated by the 2.0-0.6 Ma Yellowstone caldera and two studies [9, 10] previously speculated that the Teton and Gallatin ranges may have been initially continuous across the caldera based on their similar topographic grains and footwall strata. However, this hypothesis has remained untested and, if supported, the potential implications for fault growth, regional landscape evolution, and topographically controlled continental drainage remain unclear but could be significant.

Although it is expected that the faults and the footwall mountain ranges transected by the hotspot track would be removed and thus cannot be directly observed, it is possible to determine if there is missing fault length based on determination of maximum displacement for each fault, as there is a robust empirical relationship between the length of faults and the maximum accumulated displacement (e.g., [11] and references therein). In this scenario, if the maximum displacement accumulated on a fault can be independently determined, the displacement can be compared with empirically determined length-displacement scaling relationships to identify any potential missing fault length that may have existed prior to migration of the hotspot along the track. Although such an exercise could be carried out for any of the ranges abutting the hotspot track, the Teton normal fault is an ideal candidate for such a test, due to the age control and delineated caldera extents of the three supercaldera events that may have affected Teton Range paleotopography. However, estimates of displacement magnitude on the Teton fault and, to a lesser extent, the long-term fault slip history, have remained varied and elusive. Additionally, the position of the Teton fault at the confluence of four tectonically unique provinces, including the Basin and Range extensional province, the Yellowstone caldera and associated Snake River Plain, the Sevier fold-and-thrust-belt, and the Laramide Gros-Ventre-Wind River uplift, has further complicated a genetic understanding of Teton fault motion.

Previous displacement estimates across the Teton fault range from 2 to 11 km ([12] and references therein). This large degree of uncertainty reflects the lack of subsurface control, the challenges of imaging through hanging wall glacial till and coarse fluvial deposits, and the various criteria used to interpret offset. Despite this, most geophysical and stratigraphic studies yield maximum displacement estimates of 7-11 km [1317], with lower estimates mostly confined to studies that examine stratigraphic tilting of Pliocene-Quaternary tuffs [12, 18, 19]. Because the latter studies focus on relatively young features, they may only record more recent increments of motion and thus yield lower displacement estimates.

The idea that the Teton fault could continue much farther north than traditionally interpreted on the basis of extent topography also has potential implications for seismic hazard evaluation. The Teton-Yellowstone region represents one of the most seismically active areas in the Intermountain US [20]. However, the Teton fault remains enigmatically quiescent and has not produced an earthquake >Mw 3.0 in recorded history [2123]. The two most recent fault ruptures along the Teton fault were interpreted as Mw 6.8 and 7.1 events resulting from 1.3 m to 2.8 m of slip at ~4,800 and 7,900 ka, respectively [23], based on trenching studies [2426]. If the slip rate (0.16 mm yr-1) determined from those events is extrapolated to present-day, the predicted slip deficit of ~2 m suggests the Teton fault is capable of generating an Mw~7.0 earthquake on known segments [12, 2428]. However, if it can be demonstrated that the active Teton fault is much longer than previously recognized, this could have major implications for understanding fault length-related seismic hazard potential in this region (e.g., [29]).

Here, we synthesize the motion history of the Teton fault and explore the Teton Range removal hypothesis with a combination of new analyses and integration of previous datasets, including (1) new apatite (U-Th)/He (AHe) data for two new transects (Static Peak and Eagles Rest) and the upper part of the Mount Moran transect (Figure 2), which are integrated with previously reported AHe and apatite fission track (AFT) data [10]; (2) 1-D numerical models that constrain the likely range of geothermal gradients that can influence the closure depth of AHe and AFT data; (3) development of inverse thermal history models for all transects based on the integration of new thermochronology data and geothermal gradient models; (4) integration of previously reported [30] flexural-kinematic models of crustal-scale normal fault evolution to constrain the various contributions of footwall uplift and hanging wall drop to total displacement; and (5) displacement analyses for transects based on the integration of inverse thermal history results, calculated geothermal gradients, flexural-kinematic modeling, and stratigraphic/structural relationships. Combined, these analyses allow for predictions of northern (and southern) Teton fault projections based on fault growth models, length-displacement scaling relationships, and displacement gradient analyses. Finally, we discuss possible mechanisms for removing the northern paleo-Teton Range topography and the implications that a northern fault extension may have for future seismic hazard assessments.

2.1. AHe and AFT Thermochronology

AHe data from two new transects collected at Static Peak and Eagles Rest Peak are here added to previously reported AHe and AFT data [10, 31, 32], which includes both multisample transects and isolated lower elevation samples collected along strike along the range front (Tables 16, Figure 2). Additionally, because a previous study [10] originally included multi-grain aliquots in their study (i.e., multiple grains run as single aliquots) for multiple samples along the higher elevation parts of the Moran transect, new grains were separated for single-grain analysis of those samples. These new single-grain analyses are here combined with previous single-grain analyses from that same study to provide analytical consistency for input into inverse thermal history and displacement models. To make the dataset internally consistent, ages for the previously reported data [10, 31, 32] were recalculated using the same procedure in QTQt [33] that was used for calculating ages for the new data reported here. These recalculated ages rarely varied more than 1-2% from the originally reported ages. For both the new and previously reported AHe data, apatite was separated using standard gravity and magnetic techniques. For data derived from the [10] study, He was measured at Virginia Tech and U-Th were measured at the University of Arizona for single-grain aliquots. For the new analyses presented here, He and U-Th-Sm were measured for single-grain aliquots at the University of Illinois Helium Analysis Laboratory using standard techniques [34]. Each grain was outgassed via heating to 950-1050°C for three minutes using a diode laser. U, Th, and Sm abundances were analyzed using a PrismaPlus QMG 220 quadropole mass spectrometer. A minimum of six aliquots were analyzed for each new sample, although due to persistent challenges with apatite yield and quality for the Teton samples, AHe grain ages >±1σ from the mean were removed from the final mean calculation (~14% of the total aliquots run). AFT data for the Moran transect was derived directly from [10].

2.2. Inverse Thermal History Modeling

All inverse thermal history models were run using the software package QTQt v.64R5.6.0 [33], which employs a Markov chain Monte Carlo (MCMC) inversion approach. All data and parameters are reported here in accordance with the recommendations of previous methodological studies [35, 36]. AHe data inputs and grain measurements are shown in Tables 14. AFT data that were considered in the Moran transect are shown in Tables 5 and 6. All models include a present-day surface temperature of 4°C (average annual temperature for the Jackson Valley) and an atmospheric lapse rate of 6.5°C km-1. The geothermal gradient prescribed in the QTQt models is fixed through time (fixed temperature offset) and is derived from the modeled geothermal gradient values of 25, 27, and 29°C km-1, for the southern (Rendezvous), central (Static Peak and Grand), and northern (Moran and Eagles Rest) transects, respectively. Model priors were 60±60 Myr for time and 70±70°C for temperature, which represents the entire range of ages and closure temperatures for the AHe and AFT systems. All models included consideration of the radiation damage model [37]. AFT annealing is modeled using the multicompositional algorithms of [38, 39]. Each model was run for 10,000 burn-in iterations and 20,000 postburn-in iterations, with a uniform birth proposal distribution. Models presented are the expected T-t histories, which represent the weighted average of accepted models. The acceptance rates for each model are included in Figures 3 and 4. For a complete treatment of the MCMC inversion approach employed by QTQT, the reader is referred to [33]. The comparison between modeled and observed ages, which are essential for evaluating the validity of inverse model results, are included for each model (Figures 3 and 4).

2.3. Geothermal Gradient Modeling

Because the models presented here must include a necessary dependency on the observed or assumed thermal structure of the upper crust (~15 km), it is useful to consider any available constraints on these modeling variables, particularly because the thermal structure of the lithosphere in the Teton-Yellowstone region is complex, owing to migration of the Yellowstone hotspot. In the continental crust, the geothermal gradient is controlled by the mantle heat flow (Qm) into the base of the crust, the internal radiogenic heat production (A), and thermal conductivity (K) of the crust, the surface heat flow (Qo), and any advectionary heat transport processes.

Depth filtered P-wave seismic inversion shows that, despite being immediately adjacent to the Yellowstone hotspot, the Teton Range is underlain by a relatively fast P-wave velocity zone from 2 to 14 km depth, likely indicating a zone of relatively cold crust adjacent to the low-velocity partial melt zones beneath the present-day caldera [40]. Because the modern Teton Range appears to be underlain by relatively cold crust despite being as close to the Yellowstone hotspot as at any point in its history, we assume in our models that the present-day geothermal gradient represents the maximum geothermal gradient experienced by these rocks in the timeframe being considered (>50 Ma-present). These P-wave tomographic models agree closely with previous work [41], which showed a present-day geothermal gradient of 18-27°C for Teton region. Here, the measured surface heat flow values (Figure 5) [42] are used as constraints, and mantle heat flow and radiogenic heat production are varied to derive a range of possible geothermal gradients for the Teton Range. The range of steady-state geotherms is calculated by solving the heat conduction equation with a contribution of radiogenic heat production:
(1)Ty=To+QoKy+A2Ky2,
where T is temperature for a given depth (y), To is surface temperature, and Qo is the surface heat flow. In the geothermal gradient models presented here, heat advection is not considered, as thermal-kinematic models of Teton uplift scenarios by [30] showed minimal advectionary heat transport and yield geothermal gradients of 20-27°C km-1 for the upper 10 km of the crust. Here, the geothermal gradient is constrained using surface heat flow values of 0.075, 0.081, and 0.087 W m-2 for the southern, central, and northern Tetons, respectively [42]. By using fixed Qo values for geothermal gradient calculations, the range of possible Qm and A values can be explored, and those values for each model are shown in Table 7. Although there are multiple methods for characterizing the distribution of radiogenic heat production (A) in the crust, including layered and exponential models, it is here modeled as a constant average value throughout the entire 30 km crustal thickness. Two-dimensional thermal-kinematic models of [30] that include exponential models of heat production in the crust still yield near surface (<10 km) geothermal gradients of 20-27°C km-1, similar to modern geothermal gradients calculated by [41].

2.4. Flexural-Kinematic Modeling

Flexural-kinematic models [30] are included here to evaluate how a range of subsurface geometries for the Teton fault and the effective elastic thickness of the crust in the region could influence the footwall isostatic response during fault displacement. To constrain the range of viable geometries and crustal parameters, the topographic profiles of each model result were compared to east-west swath topographic profiles from the Gros Ventre uplift across the Teton Range near Mount Moran to the Teton Basin west of the range. Swath profiles were produced using the Swath Profiler ArcGIS add-in [43] and a 10 m DEM. The swath was aligned through Mount Moran, which was identified as the locus of maximum displacement [10]. The SwathProfiler calculates a maximum, minimum, and mean elevation profile based on topographic profiles of 50 equally spaced transects within the 6 km swath width, and these calculated profiles can then be compared with topographic profiles produced by various iterations of the flexural-kinematic model. These comparisons were then used to select a reference case for modeling the various contributions of footwall uplift and hanging wall drop to the total displacement. The details of this methodology and the results of all models are described in [30]. Initial flexural-kinematic models included near-surface fault dips ranging from 45-70°, listric detachment depths of 15-20 km, and effective elastic thickness (Te) values of 5, 10, and 15 km, based on interpretations from multiple studies of Basin and Range normal faults [4446]. Footwall erosion was not included in the simple model presented here, but [30] considered incremental footwall erosion following each displacement step that yields the range of total eroded thickness interpreted in this system. Additionally, the models of [30] included ongoing sedimentation in the modeled hanging wall to match that observed in the Jackson Hole basin. Those models of [30] yield similar results as our simple model included here. Flexural-kinematic modeling was completed using Move from Petroleum Experts. The flexural-isostatic response in Move is described by the equation for a continuous 2D beam [47]:
(2)q=Dd4wdx4+ρmρcgw,
where applied vertical load (q) required to produce a deflection (wx) is a function of the flexural rigidity (D). In this equation, ρm is mantle density (3300 kg m-3), ρc is the density of the removed crustal material (2750 kg m-3), and g is the acceleration of gravity (9.81 m s-2). Flexural rigidity (D) of the lithosphere is
(3)D=ETe3121ν2,
where E is Young’s modulus, Te is effective elastic thickness, and v is Poisson’s ratio. Only the simple reference case model (near-surface fault dip=70°, listric detachment depth=15 km, Te=5 km) determined by [30] is considered here.

3.1. Thermochronology and Inverse Thermal History Modeling

All analytical results are shown in Tables 16. As expected, AHe ages generally increase with increasing elevation along each transect; however, the age elevation gradients are highly varied between the five transects (Figure 6). The Moran transect yields the steepest age-elevation gradient. In the southern and central Teton Range, inverse thermal history models of the Rendezvous, Static Peak, and Grand Teton transects yield periods of relatively rapid Miocene to present cooling with onset ages of ~8 Ma, ~15 Ma, and ~9 Ma, respectively (Figure 3). Inverse thermal history models of the Moran and Eagles Rest transects in the northern Teton Range yield a similar footwall cooling event with onset ages of 10-8 Ma (Figure 4). With the exception of the Eagles Rest transect, inverse thermal history models of the other four transects also yield an earlier Eocene-Miocene cooling event; however, the expected T-t histories show a substantial variation in age for the earlier cooling event(s). The expected model for the Rendezvous, Static Peak, Grand, and Moran transects yield Eocene-Miocene cooling events of 48-41 Ma, 31-23 Ma, 50-26 Ma, and 25-19 Ma, respectively, and the expected model for the Moran transect also yields a separate cooling event at 61-40 Ma. However, the credible intervals increase substantially for the earlier cooling histories predicted in the Rendezvous (48-41 Ma) and Grand (31-23 Ma) transects. The Static Peak and Moran (50-26 Ma, 61-40 Ma) inverse thermal history models show narrower credible intervals for the earlier cooling event at 31-23 Ma and 25-19 Ma, respectively, and this is addressed in the discussion.

3.2. Geothermal Gradient Modeling

Twelve separate models encompassing the likely range of present-day Teton region geothermal gradients were produced (Figure 7). Models are separated into three groups representing surface heat flow values (Qo) of 0.075, 0.081, and 0.087 W m-2 that match the observed conditions in the southern, central, and northern Teton Range, respectively [42] (Figure 5). Despite the incorporation of a range of mantle heat flow (Qm) and radiogenic heat production (A) values, models yield a relatively small range of geothermal gradient values for the upper 5 km of the crust. Gradients ranged from 22 to 25°C km-1, 24 to 27°C km-1, and 26 to 29°C km-1 for the southern, central, and northern Teton Range, respectively. In order to produce the most conservative values for footwall exhumation to be interpreted from the cooling history, the maximum modeled geothermal gradient for each region was used for inverse thermal history models and footwall cooling calculations.

3.3. Flexural-Kinematic Modeling

The flexural-kinematic response of the modeled Teton Range to the entire range of parameters is reported in [30]. Here, we only report the results from displacement studies of the base case model (near-surface fault dip=70°, Te=5 km, Zd=15 km) that produced the closest match between the modeled and observed flexural-topographic profile (Figure 8). In the flexural-kinematic model, footwall uplift increases with increasing displacement (Figure 9), with the footwall uplift contribution to total displacement decreasing from ~35% to ~23% as total displacement increases from 2 to 20 km (Figure 10). Based on the flexural-kinematic modeling, footwall uplift (Ufw) can be related to total displacement (x) by a simple linear relationship:
(4)Ufw=0.2457x.

This equation is used to calculate total displacement (x) from estimated footwall uplift values determined from inverse thermal history and geothermal gradient modeling.

4.1. Onset and Slip History for the Teton Fault

Two periods of relatively rapid cooling were identified in the inverse thermal history models, with the exception of the Eagles Rest transect, which currently yields only one pronounced cooling event (Figures 3 and 4). In all transects, the more recent event is interpreted to represent the final phase of Teton fault footwall exhumation that initiated in response to fault slip from ~10-8 Ma to present. Given that our sample locations lie in the immediate footwall of the Teton normal fault, we argue that this assumption is valid. These recent cooling trends, which are recognized in all five transects, can mostly be distinguished from the earlier cooling event which is discussed below. If the more recent event represents final Teton fault motion, it yields a Teton fault slip onset of 15-8 Ma, with the majority of transects yielding onset ages of 10-8 Ma.

The older cooling events predicted by some of the inverse thermal history models (Figures 3 and 4) are more difficult to reconcile, particularly because the credible intervals increase substantially during the earlier cooling history. Also, because the inverse thermal history models include the implicit assumption that these transects are vertical, the ratio of the transect relief to the lateral distance of the transect from the fault plane will control how much of the history the transect records (greater relief records more of the cooling history) and the degree to which the AHe ages reflect fault motion (samples further from the projected fault plane will have a cooling history that is less influenced by fault slip-related cooling). In the Teton example, the most recent 15-8 Ma uplift event is predicted from the relatively young low-elevation AHe ages, but because the higher elevation samples along each transect were collected at substantially different distances from the fault plane, the response of the ages to fault motion (and the subsequent cooling and landscape response) could be impacted by such transect artifacts.

Because of these possible sampling-related artifacts, we chose to focus on the two transects (Static Peak and Moran) that include the highest total relief and are laterally closest to the interpreted Teton fault plane when interpreting this earlier cooling event (Figures 3(b) and 4(a)). Inverse thermal history models of the Static Peak and Moran transects both yield a pronounced Oligocene-Miocene cooling event (30-20 Ma) that precedes the most recent 15-8 Ma cooling event. This cooling pattern, which we interpret to be related to slip on the Teton fault, can be clearly distinguished in along-strike plots of AHe ages that highlight the variable magnitude of cooling along strike (Figure 11). Because this earlier cooling event significantly postdates the most recent phase of Laramide orogenesis in the region (Laramide, Cretaceous-Paleocene) [48], we provisionally interpret this earlier cooling event to result from the earliest normal motion on the Teton fault. If correct, this would make the Teton fault one of the oldest Basin and Range structures [49]; however, we emphasize that this interpretation remains provisional until further data and analysis can be completed in a future study.

4.2. Footwall Exhumation from Inverse Thermal History Modeling

For the most recent cooling event predicted by the inverse thermal history models, which is interpreted to reflect Teton fault motion from 15-8 Ma to present, footwall cooling estimates can be integrated with the modeled geothermal gradients to calculate footwall exhumation magnitudes. For clarity, footwall cooling is the total cooling of the highest elevation footwall sample that yields a reset AHe age (i.e., an age younger than the age of fault slip onset) from each transect and footwall exhumation is calculated by integrating the maximum modeled geothermal gradient for each transect with total cooling. These calculations yield base case footwall exhumation magnitudes of 2.6-3.0 km (25° C km-1, 64-76°C) for Rendezvous Peak, 2.7-3.0 km (27° C km-1; 73-81°C) for Static Peak, 2.4-2.6 km (27° C km-1; 66-70°C) for Grand Teton, 2.8-3.1 km (29° C km-1; 82-91°C) for Mount Moran, and 2.1-2.8 km (29° C km-1; 61-81°C) for Eagles Rest Peak. These estimates only consider the total cooling that occurred during the most recent rapid cooling event commencing at 10-8 Ma (Figures 3 and 4). If the earlier Oligocene-Miocene cooling event predicted in the Moran and Static Peak models is considered a part of the total Teton fault slip history, then these two transects yield footwall exhumation estimates of 4.1 km (29° C km-1; ~120°C) and 3.5 km (27° C km-1; ~95°C), respectively. As the Moran transect yields the greatest footwall cooling and consequently the greatest calculated footwall exhumation of all the transects analyzed, the Moran region is here interpreted to represent the approximate center of the Teton fault, following [10]. Thus, maximum displacement (Dmax) for the Teton fault is calculated below using the Moran footwall exhumation values.

The cooling magnitudes derived from the inverse thermal history models and the calculated footwall exhumation magnitudes can also be compared with similar estimates derived from AHe age-elevation relationships (Figure 6), which can be useful for simple determinations of fault slip rate variations and total exhumation magnitudes. In this relatively simple analysis, it is assumed that the closure T of the AHe system is ~70°C [50, 51], the geothermal gradient is the same as that used for the inverse thermal history model of each transect, and any sample with an AHe younger than the age of fault onset as predicted by the inverse thermal history was buried at a depth greater than the AHe Tc prior to the onset of the latest Teton fault motion (~15-8 Ma). Thus, the difference between the elevation of the highest sample with an AHe age younger than the age of fault slip onset and the predicted depth to the AHe Tc yields the predicted total exhumation for each transect. In the southern and central Teton Range, this exercise yields footwall exhumation magnitudes of ~3.1, ~2.6, and~2.7 km for the Rendezvous, Static Peak, and Grand transects, respectively. In the northern Teton Range, these calculations yield footwall exhumation estimates of ~2.9 km for both the Moran and Eagle Rest transects. If the Oligocene-Miocene cooling event is included for the Static Peak (31-23 Ma) and Moran (25-19 Ma) transects, a similar analysis yields total exhumation magnitudes of 3.3 and 4.1 km, respectively. All of these estimates are very similar to the footwall exhumation estimates calculated from the total cooling.

4.3. Total Displacement Analyses

The calculated footwall exhumation derived from the inverse thermal history models to determine the range of possible Dmax values for the Teton fault. To do this requires the implicit assumption that the magnitude of footwall uplift due to fault displacement is approximately equal to the calculated magnitude of exhumation of the lower elevation samples along the transect that yield reset AHe ages. This assumption should generally hold true here, as the samples being exhumed in the lower part of the transect were hotter than the AHe Tc when the final phase of Teton fault slip begins. As discussed above, [30] evaluated how parameters such as surface fault dip, Te, Zd, and displacement will influence the varying contributions of footwall uplift and hanging wall drop contribute to Dmax. Using the linear relationship between footwall uplift and total displacement derived from the base case flexural-kinematic model of [30] that provided the best match to the modern day flexural profile (Figure 8; fault dip=70°, Te=5 km, Zd=15 km) yields minimum calculated Dmax values of 11.4-12.6 km for Mount Moran footwall uplift of 2.8-3.1 km for the most recent phase of footwall cooling. If the exhumation estimate included the Oligocene-Miocene cooling event identified in the inverse thermal history models, this would yield total displacement values of at least 16.7 km. Because that earlier event remains poorly constrained and will be the focus of a future study, we currently consider the displacement estimates of 11.4-12.6 km to represent a minimum Dmax for the most recent phase of Teton fault motion. Displacement estimates for the other transects are included in Figure 10.

4.4. Length of the Paleo-Teton Fault

If the calculated displacement for the Moran transect represents Dmax for the entire Teton fault, that value can be used to estimate the original length of the paleo-Teton fault. Studies of fault dimensions spanning more than eight orders of magnitude ([11] and references therein, [52]) demonstrate a near linear scaling between fault length (L) and displacement and show that Dmax is always ~1 to 6% of L for large faults (Dmax>104 m; Figure 12). Using the conservative assumption that Dmax is 6% of L in this case, displacement estimates from the Mount Moran transect would yield a minimum estimated fault length of 190-210 km, which is consistent with length-displacement scaling relationships of other major Basin and Range normal faults (e.g., [53]), including the Wasatch fault [54, 55]. Provided that Mount Moran represents the approximate paleo-fault center, the Teton fault would have extended a minimum of ~85-105 km north and south of Moran prior to the onset of Huckleberry Ridge volcanism in the Yellowstone hotspot at ~2 Ma [56, 57]. These projections assume that fault propagation was not limited by interaction with other major structures. South of Mount Moran, the Teton fault continues for ~40 km prior to intersecting and offsetting the Laramide Cache Creek thrust [58]. South of that point, detailed mapping by multiple studies indicates that the Teton fault continues at least 30 km further south [5961], yielding a minimum mappable southern extent at least 70 km south of Mount Moran.

To the north, the modern topographic expression of the Teton fault footwall only extends ~20 km and modern fault scarps that are possibly associated with modern Teton fault motion only extend another ~5 km [28], which is substantially less than fault lengths predicted by LDmax relationships. As such, we interpret that the Teton fault and its associated footwall topography originally continued at least 85-105 km north of Mount Moran, 55-65 km farther than currently recognized. Because inverse thermal history models indicate that most of the Teton fault motion and related footwall topography development had accrued prior to multiple VEI 8 (volcanic explosivity index) supercaldera eruptions, including the Huckleberry Ridge (~2.0 Ma), Mesa Falls (~1.3 Ma), and Lava Creek (~0.6 Ma), we interpret that this northern paleo-Teton Range was removed/diminished following these events. In this scenario, if the Teton fault were to follow the same approximate trend, it would have extended at least as far north as Yellowstone Lake. At Mount Moran, the inverse thermal history models indicate that the lowest elevation sample had cooled to ~20°C by the onset of Huckleberry Ridge volcanism at ~2.0 Ma, and thus, all but <1.0 km of total footwall uplift was in place at that time.

4.5. Removal of the Paleo-Teton Mountain Range

Although limited direct evidence currently exists to determine what happened to the northern extent of footwall topography, it is possible that some of the range was removed by collapse into the caldera either following the initial Huckleberry Ridge and later eruptive events or as a result of posteruptive isostatic collapse (e.g., [8]). Detailed mapping has shown that both the Huckleberry Ridge and Lava Creek [62] calderas would have transected the projected northern extension of footwall topography (Figure 2), but unfortunately, any structural or deformational evidence of the paleotopography removal north of this point would be buried by the subsequent rhyolite flows [28].

Flexural modeling was used to show that arching and subsidence of Paleozoic strata and Mesozoic fold hinges in the ranges north of the eastern Snake River Plain could be accommodated by isostatic adjustment to a relatively dense midcrustal sill [28]. Models from that study showed downward flexure of 2.8-4.2 km of the upper crustal strata was possible, depending on the depth of compensation. P-wave seismic inversion was used to interpret a lower-crustal basaltic body beneath the present-day Yellowstone hotspot [40] that may produce a similar isostatic response to that modeled in [8]. Traditionally, the arching of the north dipping base Paleozoic unconformity observed on many of the Teton peaks was interpreted to act as a key pinning horizon that tracked normal motion on the Teton fault. However, contouring of the compiled AHe data (Figure 11) shows that the highest magnitude of recent cooling does not correspond to the arched shape of this unconformity, and thus, it may be possible that the northern dip on the unconformity, like the Paleozoic units north of the Snake River plain, may be the result of such subsidence. Additionally, the contoured AHe data shows a pronounced and abrupt break in AHe ages between the Mount Moran and Eagles Rest transects, with much older AHe ages observed along the Eagles Rest transect at equivalent elevations to samples on the Moran transect. Provisionally, we interpret this to indicate that Eagles Rest Peak and the associated topography to the north may have experienced postfault motion subsidence (down-to-the-north) along a structure normal to the Teton fault. Coincidentally, this break, which is approximately located between Mount Moran and Eagles Rest Peak, also corresponds to (1) a pronounced eastern shift in the modern Teton fault scarp on the south side of Eagles Rest Peak at Snowshoe Canyon (Figure 13) and (2) an abrupt drop of the northern Teton average elevation by ~0.5 km between Mount Moran and Eagles Rest Peak, as indicated from along-strike topographic swath profiling (Figure 11) [63]. Although the nature, genesis, and kinematics of this topographic break will form part of a future study, it is possible that this structure may act in a similar fashion to the hotspot track parallel normal faults recognized along the Snake River plain [61].

New inverse thermal history models based on a compilation of new and previously existing AHe+AFT data is combined with new models of variable geothermal gradients along-strike of the Teton Range to define the timing of Teton fault slip and the total magnitude of exhumation, as a function of cooling. These results are combined with flexural-kinematic modeling to determine the total displacement magnitude accrued on the Teton fault at multiple transects, and to identify the maximum accumulated displacement (Dmax) that occurred in the vicinity of Mount Moran, which leads us to interpret that area as the approximate center of the paleo-Teton fault. These displacement estimates are then combined with fault length-displacement scaling relationships of global normal fault datasets to demonstrate that the Teton fault and its associated footwall topography comprising the Teton Range may have originally been much longer than the present-day topographic expression. In this scenario, we interpret that the northern paleo-Teton Range extended 85-105 km north of Mount Moran, well into the footprint of the modern Yellowstone hotspot track. Additionally, because our inverse thermal history models indicate that Teton fault slip onset started at least by 15-8 Ma (and possibly even earlier in the Oligocene), we interpret that most of the displacement had accumulated on the Teton fault prior to migration of the Yellowstone hotspot into its current position. In this scenario, the northernmost paleo-Teton Range would have been removed following the Huckleberry Ridge (~2.0 Ma) and later eruptions and/or subsided as a part of an isostatic response to a dense basalt body in the lower crust beneath Yellowstone, similar to the interpreted evolution of multiple mountain ranges along the northern edge of the Snake River plain. Some of this topographic reduction may have been accommodated by down-to-the-north structures provisionally indicated by multiple lines of evidence. This relatively rapid removal of mountain topography (<2 Myr) represents a potentially significant control on continental physiography, biogeography, continental-scale fluvial drainage, and climate over relatively short time spans.

All data for this study are included in the manuscript tables.

Portions of the abstract were presented at the GSA Annual Meeting in Phoenix, Arizona, USA-2019.

The authors declare that they have no conflicts of interest.

This work was supported by a UW-NPS seed grant and NSF-EAR 1932808; GSA Student Research Grants to RMH, MLS, and ALH; and an AAPG Student Grant-In-Aid to MLS. The Overcash Field Fund at UK supported undergraduate student participation in fieldwork. WRG acknowledges NSF-EAR 1735788, which supports analytical facilities used for this work.

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).