Early Oligocene Surface Uplift in Southwestern Montana During the North American Cordilleran Extension

The topographic history of the North American Cordillera holds the key to understanding the tectonic and geodynamic processes of orogenesis and post‐orogenic collapse. In southwestern Montana and its adjacent regions, previous studies have suggested that the early middle Eocene rollback of the Farallon flat‐slab caused the latest topographic growth, after which regional mean elevations were progressively lowered due to extensional collapse. In this study, we report 60 new Eocene‐Miocene stable hydrogen isotopic values (δD) of hydrated volcanic glass from southwestern Montana to examine regional topographic history during the extension. The δD values show a two‐stage negative shift of ∼35‰ during the late Paleogene, the timing of which is further refined by published ages and new zircon U‐Pb ages. About 19‰ of the negative shift occurred during the late Eocene‐earliest Oligocene (∼40–32 Ma) when global and regional climate experienced cooling, and the shift can be mostly accounted for by climate change. The remaining 16‰ negative shift occurred during the early Oligocene (∼32–30 Ma), which can be best explained by regional surface uplift of 0.9 ± 0.5 km. This inferred surface uplift occurred more than 17 Myr after the initiation of the Cordilleran extension in this region, and thus challenges the traditional assumption of progressive topographic lowering during the post‐orogenic collapse. This surface uplift, if occurred, may reflect lower mantle lithosphere removal in southwestern Montana.

. Later, the regional mean elevations were inferred to reduce gradually because of subsequent extension, although local relief may have been amplified by lowering the floors of normal fault-bounded basins (Schwartz et al., 2019). However, no paleoelevation data have been reported to support this interpretation. In this study, we use the stable hydrogen isotopic values (δD) of hydrated volcanic glass to study the Eocene-Miocene paleoelevation history of southwestern Montana. The new δD data show a two-stage negative δD shift of ∼35‰ during the late Eocene-early Oligocene. While the first stage of this negative shift in the late Eocene-earliest Oligocene can be explained by climate change, the second stage in the early Oligocene most likely indicates a pulse of regional surface uplift. This finding contradicts the previous assumption of gradual elevation loss during the post-orogenic extension.

Tectonic Setting
The North American Cordilleran orogen includes the Sevier hinterland, the thin-skinned Sevier fold-thrust belt, and the thick-skinned Laramide province from the west to the east (Figure 1). In central Idaho and southwestern Montana, the fold-thrust belt contains a series of northto northwest-trending thrusts (Skipp, 1988). It has been inferred that these thrusts propagated generally eastward since the Late Cretaceous (∼90 Ma) and continued into the Paleocene (Skipp, 1988). The Laramide province was once the foreland of the Sevier fold-thrust belt and was partitioned by basement-cored uplifts into discrete basins when the Farallon slab changed from normal to low-angle subduction (DeCelles, 2004;Lawton, 2019;Schwartz & DeCelles, 1988;Yonkee & Weil, 2015). Laramide deformation in southwestern Montana was active no later than the latest Cretaceous (∼80 Ma), earlier than the Laramide deformation in Wyoming, suggesting that the Laramide deformation in southwestern Montana may not be related to flatslab subduction (e.g., Carrapa et al., 2019;Garber et al., 2020). The Idaho, Pioneer, and Boulder batholiths intruded the region during the Late Cretaceous (100-65 Ma) (Gaschnig et al., 2011;Lund et al., 2002), concurrent with the contractional deformation (Kalakay et al., 2001).
During the early Eocene, southwestern Montana transitioned from contraction to extension, characterized first by the exhumation of metamorphic core complexes as early as ∼53 Ma (e.g., Foster et al., 2007Foster et al., , 2010, and then the formation of grabens and half grabens since ∼50 Ma, which reactivated pre-existing reverse and thrust faults (Constenius, 1996). In adjacent areas, intensive eruptions of the Challis, Absaroka, and Lowland Creek volcanic fields started at the same time and ended during the middle Eocene (55-43 Ma) (Fritz et al., 2007). These magmatic events have been suggested to result from the Farallon flab-slab removal (Constenius, 1996;Smith et al., 2014). Local magmatism in southwestern Montana, such as the Dillion volcanics, lasted into the Neogene (∼17 Ma) (Fritz et al., 2007). Regional extension in southwestern Montana also lasted into the Neogene, concurrent with the Basin and Range extension in Nevada and Utah (Constenius, 1996;Sears et al., 2009).
The Renova Formation and its equivalents consist mainly of alluvial-fan, fluvial, and lacustrine sandstone and mudstone, with minor conglomerate (Janecke et al., 1999;Kuenzi & Fields, 1971;Monroe, 1981;Portner et al., 2011;Schwartz & Schwartz, 2013); while the Sixmile Creek Formation is dominated by coarsegrained alluvial-fan and fluvial sandstone and conglomerate (Fritz & Sears, 1993;Kuenzi & Fields, 1971;Sears et al., 2009). Both formations outcrop abundant volcaniclastic ash-fall and ash-flow materials, containing glass shards and zircons. Detailed descriptions of sedimentology and lithofacies of each studied basin are summarized in Text S1. The geological map and stratigraphic columns of the Sage Creek Basin and field photos of most sampled outcrops are also provided in Figures S1-S4.
The subdivisions and nomenclatures of the Renova Formation and its equivalents are variable among basins in southwestern Montana ( Figure 2). The Renova Formation in the Sage Creek Basin has the longest and most continuous sedimentary records, and the formation is divided into the Sage Creek, Dell, Cook Ranch, White Hills, and Blacktail Creek members, from old to young (e.g., Schwartz & Graham, 2017). A ∼5 Myrs unconformity exists between the Dell and Cook Ranch members in the Sage Creek Basin (Figure 2). This local unconformity belongs to the major regional unconformity between the Cretaceous-lower Eocene contraction-related deposits and the upper Paleogene extension-related deposits (Schwartz & Graham, 2017). Deposition of the Renova Formation in the other basins started later and has different nomenclatures (Rasmussen, 2003). Details of the nomenclatures of each studied basin are summarized in Figure 2 and discussed in Text S1.
(Figures S1 and S2). Except for the Blacktail Creek Member, samples from the other members were collected from sections that were previously measured and studied for carbonate stable oxygen isotopes (Kent-Corson et al., 2010Methner et al., 2016;Schwartz et al., 2019).
Five samples from the key intervals of the δD shift were studied for zircon U-Pb geochronology. The samples include three tuffaceous sandstone samples from the Sage Creek Basin, one tuffaceous sandstone sample from the Smith River Basin, and one tuff sample from the Gravelly Range (Table 1).

Analysis of Volcanic Glass Hydrogen Isotopes
The tuff and tuffaceous sandstone samples were gently crushed and soaked in weak hydrochloric acid (HCl) to remove carbonate. After wet sieving, fractions between grain size 63-125 μm were treated with heavy liquid and magnetic methods to separate and concentrate glass shards to >95% purity. The separated shards were soaked in 5% hydrofluoric acid (HF) for ∼30 s to remove surfaces that may have been altered into clays. The δD values of purified glass shards were analyzed on a Finnegan Delta V Mass Spectrometer with a Finnegan TC/EA peripheral in the Environmental Isotope Laboratory at the University of Arizona. See Text S2 for detailed sample pretreatment and analytical methods, and Table S1 for stable hydrogen isotopic data, and Table S2 for all raw data.

Analysis of Zircon U-Pb Ages
For U-Pb analysis, zircon grains were extracted from sandstone and tuff samples using the standard gravity, heavy liquid, and magnetic separation methods. Mounted and polished zircons were analyzed using a New Wave 213 nm solid-state (Nd: YAG) laser ablation system coupled with a Thermo Finnigan Element 2 Single-collector Inductively Coupled Plasma Mass Spectrometry at the Washington State University. See Text S3 for detailed analytical method and data processing, and Table S3 for zircon U-Pb age data.

Analysis of Major Elements of Volcanic Glass
We also analyzed the major elements of glass shards of nine randomly selected tuff samples using a Shimadzu XRF-1800 Sequential Wavelength-dispersive X-ray Fluorescence Spectrometer at the University of Texas at Arlington. See Text S4 for detailed analytical methods, and  (Vermeesch, 2018). c The weighted mean age of the youngest distinct age group including at least two grains. d TuffZirc age algorithm from Ludwig and Mundil (2002).

Age Constraints of Volcanic Glass Samples
Detailed age constraints of each basin and sample are described in Text S1. In brief, for samples from the Renova Formation in the Sage Creek Basin (Figures S1 and S2), the depositional age of each sample was determined based on linear extrapolation between the ages of the top and bottom of each member by assuming a constant sediment accumulation rate. This assumption may not be accurate given sediment accumulation rate is unlikely constant through time. However, because the glass samples were collected from several different basins, this approach offers the best way to display the stable isotopic data in figures to examine the temporal patterns. On the other hand, the age range of each member of the Sage Creek Basin is relatively small (i.e., <5 Myr; Figure 2), thus the age errors associated with this linear extrapolation approach are relatively small.
Depositional ages in the other basins are constrained by regional geological maps and published literature ( Figure 2 and references therein). As in most cases, only one glass sample was analyzed from each member in each basin, the sample age and uncertainty are assigned as the age mean and range of each unit. If a sample was collected adjacent to stratigraphic levels with published dates, we use these dates to refine the age of the sample. These approaches of age determination do not influence the interpretation of the timings of isotopic shifts, because the timings are further clarified by ages of the key stratigraphic levels published before and newly dated in this study (Table 1).
Four of the five zircon U-Pb samples studied here yield depositional ages between 23 and 32 Ma (Figure 3 and Table 1). The ash-fall tuff sample from the Gravelly Range (17GR01) yields a weighted mean age of 31.1 ± 0.1 Ma (n = 39, MSWD = 1.0). The tuffaceous sandstone samples from the Smith River Basin (17SR04) and the Cook Ranch Member of the Sage Creek Basin (17SC52) yield abundant young grains. We adopt the youngest deconvoluted peak age constrained by the unmixing method (Vermeesch, 2018) as the maximum depositional ages (MDAs) of these two samples. The MDAs are 30.4 ± 0.7 Ma and 31.9 ± 0.2 Ma, respectively. Several previous studies have shown that detrital zircon U-Pb MDAs of late Paleogene strata in the western USA nearly represent the true depositional ages because of the presence of ash-fall zircons from synchronous volcanic eruptions (e.g., Fan et al., 2015). These newly determined MDAs are equivalent to the ages constrained using the stratigraphic approaches described above, suggesting that the age assignments of the other samples are appropriate.
The tuffaceous sandstone sample (17SC60) from the Blacktail Creek Member of the Sage Creek Basin yields only two young grains, and we use the weighted mean of the two ages as its MDA: 23.8 ± 0.3 Ma (n = 2, LI ET AL. 10.1029/2020TC006671 5 of 17  (Schwartz et al., 2019). Vertical axis represents the number of zircons; and n = number of zircons <100 Ma/total number of zircons analyzed. The plots are prepared using online software IsoplotR (Vermeesch, 2018). MDA, maximum depositional age; YSG, youngest single grain (see Table 1 for details). MSWD = 0.1). This MDA is not as robust as the other two sandstone samples because of the small number of young grains (Dickinson & Gehrels, 2009); but it is reasonable given it agrees with a previously published tuff age of 23.0 ± 0.7 Ma in the lower Blacktail Creek Member at a different site in the Sage Creek Basin (Schwartz et al., 2019). The other sample (17SC64) from the Blacktail Creek member and the same location as 17SC60 has one single youngest zircon age of ∼27 Ma. This age is older than the MDA of 17SC60, which is stratigraphically below, suggesting that the youngest zircon grain of 17SC64 was recycled from older strata, for example, White Hills Member, but not a synchronous ash-fall volcanic grain.

Stable Hydrogen Isotopic Results
The glass shards separated from all the samples exhibit no apparent difference in color or texture under a binocular microscope. The nine randomly selected glass samples are all felsic based on that their SiO 2 contents varying between 64% and 74% (Table S4), suggesting that the other volcanic glass shards are also likely felsic. Most of the glass samples have water contents between 2 and 8 wt.% (Table S1). Six samples with <2 wt.% water contents were excluded from further discussion because they were not hydrated completely and the isotopic values were possibly influenced by magmatic water (Friedman et al., 1993;Seligman et al., 2016), although their δD vg values are not higher than those of other samples at adjacent stratigraphy levels.
The volcanic glass δD vg values range between −197‰ and −73‰ (Table S1 and Figure S5), and the corresponding δD mw values of paleo-meteoric waters range between −169‰ and −41‰ (Figure 4). The magnitude of δD mw change is nearly the same as the magnitude of δD vg change; as a result, in this manuscript, if no subscript of δD is denoted, it means that the description is valid for both the volcanic glass and meteoric water.
The first-order depositional environment (e.g., alluvial-fan, lacustrine, and fluvial) of each volcanic glass sample was summarized from previous sedimentological studies (Janecke et al., 1999;Kuenzi & Fields, 1971;Monroe, 1981;Portner et al., 2011;Schwartz & Schwartz, 2013). The δD values and depositional environments of the samples are not correlated (Figure 4). The lacustrine samples show similar δD values as those of the alluvial-fan and fluvial samples, suggesting that the lakes were hydrologically open and experienced little evaporative enrichment of hydrogen isotopes. Nine samples from alluvial-fan and fluvial environments have significantly higher (>20‰) δD values than other samples from adjacent stratigraphic levels ( Figure 4). The high values are scattered, thus unlikely caused by climate change or tectonics. These values are likely resulted from the evaporation of soil water or short-lived (e.g., thousands of years) pond water in which the glass shards were hydrated. Evaporated water bodies (e.g., swamp, vadose zone) are common on river floodplains and alluvial-fan plains. This speculation can be tested by detailed lithofacies analysis and micro-depositional environment interpretation of each sample in future studies. The nine samples were excluded when calculating the mean to display the general trend of the isotopic data.
The major observation made from the data set is that the δD values show a ∼35 ± 6‰ negative shift during the late Eocene-early Oligocene (Figure 4). The value is calculated as the difference of the means of samples between 42-40 Ma and 31-30 Ma, and the uncertainty of the difference is propagated from the uncertainties of the two means. This negative shift occurs in two stages, including a ∼19 ± 7‰ decrease during the late Eocene, and a ∼16 ± 6‰ decrease during the early Oligocene. These two stages show up both in the LI ET AL.
10.1029/2020TC006671 6 of 17 consistent with the trend of carbonate δ 18 O c values that was interpreted to reflect gradual aridification (Schwartz et al., 2019).

Timings of the Isotopic Shifts
In the Sage Creek Basin, the late Eocene negative δD shift (first stage) initiated at least by the initial deposition of the Cook Ranch Member. The shift likely initiated earlier because the Cook Ranch Member overlies a regional unconformity (Figure 4). The Cook Ranch Member initiated just before 35.4 ± 0.5 Ma based on one tuff zircon U-Pb age (Kent-Corson et al., 2006). The top of the Dell Member underlying the regional unconformity is 40.1 ± 0.8 Ma (Methner et al., 2016). These ages constrain the initiation of the negative δD shift to 40-35.4 Ma, during or after the formation of the major unconformity between the Dell and Cook Ranch members. Samples from the other basins also show that the shift began at the same time as in the Sage Creek Basin (Figure 4). This late Eocene stage of negative shift probably lasted across the Eocene-Oligocene transition and into the earliest Oligocene because of the uncertainty of age constraints (Figure 4).
In the Sage Creek Basin, the peak of the negative δD shift (second stage) occurs at ∼24 Ma in the lower Blacktail Creek Member, which is determined by the new MDA of 23.8 ± 0.3 Ma. In contrast, in the other basins, the peak of the negative δD shift is at ∼31 Ma based on one tuff zircon age and one tuffaceous sandstone MDA presented in this study ( Figure 4). The two most negative δD values at ∼24 Ma in the Blacktail Creek Member of the Sage Creek Basin are from tuffaceous sandstones, while one of the most negative δD values at ∼31 Ma in the other basins is from a tuff (Figure 4). The tuff samples show no evidence of reworking, and the volcanic glass in tuffs must be hydrated in-situ soon after deposition, and thus the δD vg values reflect synchronous local surface water δD mw values. Unlike those in tuffs, the volcanic glass shards in tuffaceous sandstones could be recycled shards that were deposited and hydrated before erosion and transported as the detritus of sandstones.
The different timings of peak δD shift between the Sage Creek Basin (determined by tuffaceous sandstones) and the other basins (determined by both tuffaceous sandstones and ash-fall tuffs) may be either a true signal or an artifact of glass shards recycling. Previous studies have shown that some Neogene distal tuffs in the western USA contain abundant glass shards but scarce zircon grains (e.g., Fan et al., 2015). Our zircon U-Pb data show that the two tuffaceous sandstone samples from the Blacktail Creek Member (17SC60 and 17SC64, Figure 3) have very few to none syndepositional ash-fall zircon grains, but abundant recycled zircon grains older than 60 Ma, suggesting that the glass shards in these samples are either syndepositional ash-fall shards or recycled from contemporaneous ash falls in highlands or slightly older strata of the Renova Formation, such as the White Hills Member. Comparison of the detrital zircon U-Pb age spectra between the two Blacktail Creek Member samples of this study and one White Hills Member sample of Schwartz et al. (2019) show great similarity (Figure 3), supporting the inference that the Blacktail Creek Member contains recycled detritus, including both zircons and volcanic glass shards, from the White Hills Member.
In the case of recycling, the δD vg values of the Blacktail Creek Member reflect the δD mw values of meteoric water during the deposition of the White Hills Member between 30 and 25 Ma, rather than the δD mw values of surface water during the deposition of the Blacktail Creek Member between 25 and 20 Ma. As a result, the peak of the negative shift was estimated to be as early as ∼31 Ma and may be extended to ∼25 Ma.

Volcanic Glass Hydrogen Isotopes Record Primary Earth Surface Signals
When volcanic ashes fall on Earth's surface and enter the depositional environment, the glass shards begin to get hydrated through the diffusion of surrounding waters. The hydration can be completed within a few tens of thousands of years, and the hydration rate decreases significantly after the first few thousand years (Friedman et al., 1993). Considering the common sediment accumulation rates of 0.01-1 mm/yr, the completion of glass hydration would finish when the shards are only shallowly buried. The sediment accumulation rate of the Renovo Formation is generally ∼0.1 mm/yr (Schwartz & Graham, 2017), thus the hydrogen isotopes of volcanic glass in the Sage Creek Basin should have recorded the hydrogen isotopes of surface water, and the temporal trend after excluding the nine samples with scattered high δD values reflect local surface water δD mw trend with minimal influence of evaporation.
A complication is that laboratory experiments have shown that after the initial hydration, the hydrogen isotopes of glass shards can exchange with those of water when the shards are in contact with water of different isotopic composition (Nolan & Bindeman, 2013). A recent study of both Eocene and Holocene ashes shows that the effect of hydrogen isotopic exchange can be removed by treating glass shards with HF (Cassel & Breecker, 2017). This is because the hydration of volcanic glass forms a dense gel layer around glass shards (Cailleteau et al., 2008) that limits weathering and secondary isotopic exchange to the outer rim and the HF treatment can effectively remove the outer rim (Cassel & Breecker, 2017). However, the lack of independently constrained water isotopic values makes it difficult to understand how well the HF treatment works for glass samples that are tens of millions of years old. The negative shift of volcanic glass hydrogen isotopic values and the corresponding negative shift in the carbonate stable oxygen isotopic values around 47 Ma ( Figure 5) seems to show that hydrogen isotope exchange is less a concern in our data set.

Controlling Factors of Carbonate and Volcanic Glass Stable Isotopic Values
Soil carbonate δ 18 Ο c values have been extensively reported in the Sage Creek Basin and several other basins in southwestern Montana and western-central Idaho (see summary in Schwartz et al., 2019). The stable oxygen isotopes of soil, lacustrine, and shallow groundwater carbonate are generally in equilibrium fractionation with those of ambient water; and the δ 18 Ο c value is determined by both the temperature of carbonate formation and the δ 18 Ο mw value of the paleo-meteoric water in which the carbonate is precipitated (Kim & O'Neil, 1997). Evaporation of surface waters can increase the δ 18 Ο mw values of soil and lake water, and thus the δ 18 Ο c values of carbonates. The δ 18 Ο c value of soil carbonate is specifically determined by the δ 18 Ο mw value of soil water, which typically represents local precipitation with the potential influence of evaporation, as well as the soil carbonate formation temperature that may be seasonally biased (Breecker et al., 2009;Hough et al., 2014;Quade et al., 2013).
Because the δ 18 Ο mw and δD mw values of meteoric water are highly related to each other (Craig, 1961), the glass δD vg record can be compared with the previously reported carbonate δ 18 Ο c record to understand the controls on the isotopic shifts ( Figure 5). Different from oxygen isotopes, the fractionation of stable hydrogen isotopes between volcanic glass and ambient water is driven by diffusion; and thus the δD vg value is controlled only by the δD mw value at surface temperatures (Cassel & Breecker, 2017;Friedman et al., 1993). For ash-fall glass shards, their δD vg values record the δD mw values of water in the depositional environment. For tuffaceous sandstone samples, as the glass shards may be recycled, the environmental water that their δD vg values record can be various. For example, depending on the rates of erosion, transport, and burial, the glass shards can be hydrated in the sediment source region to record highland δD mw values before transport and deposition, or during transport to record river water δD mw values, or after deposition to record the δD mw values of waters in the depositional environment. Despite these complexities, in this study the temporal δD vg pattern of the in situ ash-fall tuffs is similar to that of the tuffaceous sandstone samples during the late Eocene-Oligocene (Figure 4), suggesting that the reconstructed δD mw patterns from the two groups of samples were controlled most likely by the same mechanism, for example, decrease of regional precipitation δD mw values.

Similar Trends During the Early Middle Eocene (48-46 Ma)
Both the published soil carbonate δ 18 Ο c and the new volcanic glass δD vg records show an abrupt, negative shift at ∼47 Ma, between the lower and upper Sage Creek Member ( Figure 5). The negative shift is accompanied by a soil carbonate δ 13 C c increase of at least 2‰ (Kent-Corson et al., 2010;Schwartz et al., 2019). Although the negative δ 18 Ο c shift was not associated with any sediment provenance change, Kent-Corson et al. (2010) and Schwartz et al. (2019) suggested that it reflects a significant increase of groundwater or river water into soil water in association with changing fluvial catchment from a localized low-elevation catchment to a distal high-elevation hinterland catchment in central Idaho. This drainage reorganization is nearly contemporaneous with the catchment change of the Lake Gosiute from a local catchment in southwestern Wyoming to a distal catchment in central Idaho at ca. 49 Ma (Carroll et al., 2008). The inferred drainage shifts have been suggested to be temporally associated with the surface uplift of central Idaho, due to the magmatism of the Challis volcanic field, which is part of the southwestward magmatism sweep generated during the Farallon flat-slab rollback (e.g., Kent-Corson et al., 2010;Schwartz et al., 2019;Smith et al., 2014).

If this drainage shift interpretation is correct, it suggests that the analyzed carbonates in the upper Sage
Creek Member should be, at least partially, groundwater cement carbonates, to record the isotopic values of river water. This inference is consistent with the stable carbon isotopic data of these samples. The high δ 13 C c values (−3 to 1‰) in the upper Sage Creek Member (Kent-Corson et al., 2010;Schwartz et al., 2019) could not be solely attributed to lowered soil respiration rate in a C 4 plant-free ecosystem during the Eocene. More likely, dissolution of marine carbonate rocks, which have higher δ 13 C c values, in river water or shallow groundwater also contributed to the high δ 13 C c values of the upper Sage Creek carbonates. This interpretation suggests that the carbonates, at least partially, are shallow groundwater cement carbonates.
The uncertainties of carbonate type (soil vs. cement carbonate) and water that hydrated the tuffaceous volcanic glass shards (water in source terrane vs. water in depositional environment) permit alternative explanations for the negative δ 18 Ο c and δD vg shifts between the lower and upper Sage Creek Member in the Sage Creek Basin. In addition to the previously proposed drainage shift to reflect highland precipitation isotopes in the hinterland (Kent-Corson et al., 2010;Schwartz et al., 2019), these isotopic shifts could also reflect the decrease in local precipitation oxygen and hydrogen isotopic values, associated with regional surface uplift in the basins. These two possible explanations are not mutually exclusive, it is likely that removal of the Farallon flat slab from southwestern Montana (Chamberlain et al., 2012;Mix et al., 2011) would have caused regional surface uplift that not only influenced the hinterland, but also the basins. The negative δ 18 Ο c and δD vg shifts cannot have resulted from reducing evaporation from the lower to the upper Sage Creek Member because less evaporation cannot explain the accompanying positive δ 13 C c shift.

Different Trends During the Late Middle Eocene-Early Oligocene (45-30 Ma)
A major difference exists between the two datasets for the middle Eocene-early Oligocene ( Figure 5). The carbonate δ 18 Ο c values show a positive shift of ∼1‰ from the late middle Eocene to the earliest Oligocene and then remained nearly the same during the early Oligocene (Schwartz et al., 2019). In contrast, the δD vg values show a two-stage negative shift of ∼35‰, including a ∼19‰ shift during the late Eocene-earliest Oligocene (including age errors) and another ∼16‰ shift during the early Oligocene (32-30 Ma).
Global temperature decreased gradually during the late Eocene and abruptly across the Eocene-Oligocene transition (EOT) (Zachos et al., 2008). Regional climate cooling has been reported in Wyoming based on carbonate clumped isotope temperatures, including a ∼5°C decrease during the middle-late Eocene , and an additional ∼7°C decrease across the EOT (Fan et al., 2017. The cooling should lead to a gradual decrease of meteoric water isotopic values (Dansgaard, 1964), which has been documented in shallow groundwater δ 18 Ο mw values in Wyoming (Fan et al., 2017. Recent modeling studies also suggest that precipitation δ 18 Ο mw values are more sensitive to temperature changes in continental interiors than in coastal areas (Kukla et al., 2019;Winnick et al., 2014). However, cooling also leads to aridification and addition of recycled moisture in continental interiors, both tend to increase the precipitation isotopic values (Winnick et al., 2014). The balance between the amount of moisture recycling and temperature drop may determine the net change of isotopic values in a given region.
Our gradual decrease of δD mw values follows the same trend of the reconstructed δ 18 Ο mw values in Wyoming , suggesting that the negative δD mw shift during the late Eocene-earliest Oligocene was most likely a result of global cooling. If the magnitude of temperature decrease in southwestern Montana followed that in Wyoming, a ∼12°C temperature decrease during the late Eocene and across the EOT can cause a ∼67‰ decrease in δD mw record if the modern temperature-δD mw relationship of −5.6‰/°C (Dansgaard, 1964) is applicable to the Eocene. The estimated value of shift is much greater than the observed ∼19‰ decrease. The difference may be attributed to (a) uncertainty of temperature decrease in southwestern Montana; (b) moisture recycling in a dry climate; and (c) increase of oceanic vapor source isotopic values associated with Antarctic glaciation. However, the effects of the three factors cannot be quantitatively determined.
In contrast to our late Eocene-earliest Oligocene δD vg record, soil carbonate δ 18 Ο c values show at most a ∼1‰ increase during this period ( Figure 5), which was interpreted to be mainly a result of increasing moisture recycling when atmospheric moisture content decreases because of climate cooling (Schwartz et al., 2019). While we agree progressive addition of recycled moisture associated with drying is likely a cause for the slight increase of δ 18 Ο c values, the influence of temperature on the δ 18 Ο values of both precipitation and carbonate cannot be ruled out in the continental interior setting (Kukla et al., 2019). A decrease in carbonate formation temperature increases δ 18 Ο c values by ∼0.24‰/°C if the δ 18 Ο mw value remains stable (Kim & O'Neil, 1997). If the carbonate formation temperature in southwestern Montana followed that in Wyoming, a ∼12°C temperature decrease can cause a ∼2.9‰ increase in the δ 18 Ο c record. On the other hand, assuming the modern global mean slope of 8 for the δ 18 Ο mw and δD mw relationship (Craig, 1961) is applicable to the Eocene, the observed ∼19‰ decrease of δD mw values would indicate a 2.4‰ decrease of δ 18 Ο mw values, nearly canceling out the temperature effect on the isotope fractionation in carbonate crystallization. Despite the fact that the δ 18 Ο mw −δD mw relationship and temperature change in southwestern Montana during the late Eocene-earliest Oligocene were not well constrained, the reasoning here shows that the different trends of δD vg and δ 18 Ο c can be reconciled.
During the early Oligocene (∼32-30 Ma), global temperature did not experience major changes ( Figure 4; Zachos et al., 2008); the δD vg values decreased ∼16‰ while the δ 18 Ο c values remained nearly stable or slightly decreased ( Figure 5). Below we attribute the δD vg decrease to ∼0.9 km of regional surface uplift. This magnitude of uplift could lead to 3°C-5°C cooling following the modern elevation and air temperature relationship (Wolfe, 1992). This cooling can cause a 0.7‰-1.2‰ increase in δ 18 Ο c values, and the increase can be canceled out by the decrease in δ 18 Ο mw values associated with the inferred ∼0.9 km surface uplift. These two contrasting effects could lead to the observed δ 18 Ο c trend during the early Oligocene.

Early Oligocene Surface Uplift in Southwestern Montana
The reconstructed δD mw values decreased ∼16 ± 6‰ in the upper Cook Ranch Member in the Sage Creek Basin as well as its age-equivalents in the other basins. This stage of shift occurred after the EOT and before 30 Ma (∼32-30 Ma; Figure 4), which is well constrained by the new zircon U-Pb ages.
This early Oligocene δD mw shift cannot be explained by glass composition change or water content difference. A recent study shows that the isotopic separation factor (ΔD vg−mw = δD vg −δD mw ) of felsic and mafic glasses are −33‰ and −23‰, respectively (Seligman et al., 2016). The major element data of the nine randomly selected tuff samples indicate all the glass shards are felsic (Table S4). The pristine volcanic glass contains 0.1-0.3 wt.% magmatic water (Friedman et al., 1993). Incomplete hydration of volcanic glass leads to low meteoric water content and an increase of δD vg values by the influence of magmatic water (e.g., Seligman et al., 2016). The data do not show a correlation between δD vg values and water content (Figures 6a  and 6c), suggesting that the shift was not related to changing degree of hydration.
The early Oligocene δD mw shift is less likely explained by glass hydration with waters at different locations.
Although it is difficult to differentiate if the glass shards in the tuffaceous sandstone samples in the Sage Creek Basin were hydrated in highland before transport or within basin during deposition, the negative isotopic shift is also present in the ash-fall tuff samples. This indicates that the δD mw shift reflects the isotopic change of local surface water.
The early Oligocene negative isotopic shift is also less likely explained by glass hydration in different types of water. River, lake, and shallow ground waters may have different δD mw values because the water source and influence of evaporation are variable. The upper Eocene-Oligocene Renova Formation in southwestern Montana was deposited in a range of depositional environments, including alluvial-fans along basin margins and rivers and ephemeral lakes in basin centers (Janecke et al., 1999;Monroe, 1981;Portner et al., 2011;Schwartz et al., 2019). In the Sage Creek Basin, new observations and a previous sedimentology study (Schwartz & Graham, 2017) indicate that the Dell and Cook Ranch members were both deposited in alluvial-fan environments. In the other basins, in addition to alluvial-fan and fluvial environments, shallow lacustrine environments were occasionally present (see Text S1 for details). However, we do not observe a clear relationship between the δD vg values and their depositional environments (Figure 4), suggesting that this stage of δD vg shift is not related to changing depositional environments.
The early Oligocene negative isotopic shift is less likely explained by moisture source change. It has been suggested that southwestern Montana was in the rain shadow of vapor from the Pacific because high Cordilleran hinterland blocked moisture from the west, and the North American summer monsoon brought Gulf of Mexico moisture and recycled continental moisture to the Cordilleran foreland as far north as southwestern Montana during the Eocene (Feng et al., 2013). The high Cordillera hinterland in the western USA lasted into the late Oligocene (Cassel et al., 2018), suggesting that the dominant moisture source in southwestern Montana unlikely changed.
After excluding these potential causes, we suggest that the early Oligocene negative δD mw shift most likely reflects regional surface uplift. This uplift was not localized to a few mountain ranges because the glass samples were collected from a broad range of elevations ( Figure 6b) and basins in southwestern Montana. Although it is likely that an unknown portion of glass shards in the tuffaceous sandstone samples was hydrated in highland areas before being transported to basin floors by rivers, the ash-fall tuff samples without reworking also show the negative shift initiated at the same time ( Figure 4). Modern stable hydrogen isotopic lapse rate in the Cordilleran orogen varies from ∼15 to ∼30‰/km in the semi-arid central Rockies (Zhu et al., 2018) and Canadian Rockies (Yonge et al., 1989), and the relatively wet Sierra Nevada in southern California (Friedman et al., 1992). When a rate of 22.5 ± 7.5‰/km is adopted for the early Oligocene, the ∼16 ± 6‰ δD mw decrease reflects surface uplift of ∼0.9 ± 0.5 km. This estimate is likely conservative given that the peak shift may occur between 30 and 25 Ma in the Sage Creek Basin, which probably has even more negative δD mw values (Figure 4).

Implications for Regional Geodynamics
Various geodynamic mechanisms have been proposed to explain the topographic growth of mountain belts and plateaus. Surface elevations of an orogenic belt can be increased by isostatic compensation of a deep crustal root due to crustal thickening as a result of distributed horizontal shortening in a fold-thrust belt (e.g., DeCelles & Coogan, 2006), underthrusting of the whole or lower portion of crust (Powell & Conaghan, 1973), injection or flow of lower crustal materials (Clark & Royden, 2000), and magmatic intrusions (Perkins et al., 2016). Surface elevations can also be increased by increasing lithosphere buoyancy through hydration, or removal of dense materials, or by thermal expansion (Axen et al., 2018;Bird, 1979 Previous studies have indicated that southwestern Montana has obtained >2 km mean elevation mostly by crustal shortening during the Late Cretaceous-Paleocene Sevier-Laramide deformation (Figure 7a; Schwartz et al., 2019). The better and newer constraint of this study documents an additional ∼0.9 km of surface uplift in the region during the early Oligocene (∼32-30 Ma; Figure 7c). This uplift is at least 17 Myr after the initiation of extension in southwestern Montana since ∼53-50 Ma (Constenius, 1996;Foster et al., 2010). Crustal shortening and underthrusting can be excluded as drivers of this stage of surface uplift because contraction had ended by the early Eocene (Schwartz & DeCelles, 1988;Skipp, 1988). Lower crustal materials can flow from areas with thick crust to areas with thin crust and thus increase the elevations of the latter due to isostasy (e.g., Clark & Royden, 2000). However, southwestern Montana should have achieved thickened crust by the end of contraction (Gaschnig et al., 2011;Skipp, 1988), and thus lower crustal flow was unlikely during the early Oligocene. Magmatic intrusion cannot be a major driver because of the scarcity of late Paleogene magmatism in southwestern Montana (Fritz et al., 2007). Thermal and isostatic effects related to flat-slab subduction and rollback can also be excluded because the Farallon flat slab, if it had impacted southwestern Montana, should have passed the region by the middle Eocene (Figure 7b; Cassel et al., 2018;Smith et al., 2014). This same reason can exclude hydration of the lower lithosphere as a driver of the surface uplift (Jones et al., 2015) because the water of hydration is typically derived from the subducted oceanic plate.
We postulate the early Oligocene surface uplift in southwestern Montana as a result of removing the lower mantle lithosphere and its ensuing upwelling of the asthenosphere (Figure 7c). Lithosphere removal is common in orogens (e.g., England & Houseman, 1989;Kay & Kay, 1993;Krystopowicz & Currie, 2013) and can cause rapid surface uplift (e.g., Garzione et al., 2008;Sundell et al., 2019). We infer that the prolonged Cretaceous-early Paleogene contraction and lithosphere thickening may have set the stage for lithosphere removal. The subsequent heating associated with flat-slab rollback may have facilitated the removal by weakening the lithosphere. However, Saylor et al. (2020) suggest that the lithosphere in southwestern Montana and Wyoming has weakened since the Late Cretaceous. Regional magmatism during the Oligocene (Fritz et al., 2007) is scarce in southwestern Montana, which seems not to support the inference of lithosphere removal. However, the scale of the inferred lithosphere removal may be small and thus the ensuing asthenosphere upwelling was not extensive for magmatism (Drew et al., 2009).

Conclusion
This study reports 60 new volcanic glass δD values and four zircon U-Pb ages to constrain the middle Eocene-Oligocene topographic history of southwestern Montana. The δD data show a negative ∼18‰ shift at ∼47 Ma in the Sage Creek Basin, consistent with the carbonate δ 18 O c decrease recorded in previous studies. This shift may be due to drainage reorganization in the Sage Creek Basin caused by surface uplift in central Idaho or regional surface uplift that raised both the basin floor and adjacent highlands, associated with the Farallon flat-slab rollback. The major observation made from the new data set is the late Eocene-Oligocene negative δD shift up to 35‰. The shift initiated between 40 and 36 Ma and lasted into the early Oligocene (∼30 Ma) based on published and new ages presented in this study. About 19‰ negative shift of the first stage can be explained by climate change during the late Eocene and across the Eocene-Oligocene LI ET AL.

10.1029/2020TC006671
13 of 17 Figure 7. Tectonic evolution model of southwestern Montana and adjacent areas from the latest Cretaceous to early Oligocene. The dashed gray line represents the surface topography of the previous stage; white upward arrows indicate surface uplift. FTB, fold-thrust belt; VF, volcanic field. (a) During the latest Cretaceous-early Eocene, the Cordilleran orogen experienced shortening and crustal thickening due to subduction of the Farallon oceanic plate. Due to isostasy, the belt has attained elevations >2-3 km. (b) During the middle Eocene, in response to the rollback of the Farallon slab, asthenosphere upwelling caused eruptions of large volcanic fields and minor surface uplift, accompanied by regional drainage reorganizations. (c) During the early Oligocene (>17 Ma after the initiation of extension in southwestern Montana), the lower mantle lithosphere was removed after long-term crustal thickening and ensuing weakening, which caused a further increase of regional mean elevations.
transition. The remaining ∼16‰ negative shift occurred during the early Oligocene (∼32-30 Ma) and was most likely caused by regional surface uplift. This interpretation, if correct, contradicts the previous assumption that after the small amount of surface uplift associated with the rollback of the Farallon flat slab, the Cordillera orogen gradually lost elevation during the ensuing extension. Instead, our results suggest that localized surface uplift may have occurred during the post-orogenic collapse, which is >17 Myr after the initiation of regional extension in southwestern Montana. Although we favor lower lithosphere removal as the driver of the inferred early Oligocene surface uplift, more studies should be conducted to verify the surface uplift signal and test the geodynamic driver(s).