Colorado’s High Plains stand at anomalously high elevations (~1300–2100 m) for their continental interior setting, but when and why this region became elevated is poorly understood. The Cenozoic history of the High Plains is also likely linked with that of the Rocky Mountains, where the timing and cause(s) of uplift are similarly debated. We present apatite (U-Th)/He (AHe) data for 10 samples from Tertiary intrusives along a ~200 km west-to-east transect across the High Plains of southeastern Colorado to constrain the timing of exhumation and to gain insight into when and why regional elevation gain occurred. Mean sample AHe dates for the ~24–22 Ma East Spanish Peak pluton and associated radial dikes from the westernmost High Plains are 18.8 ± 1.4 to 14.1 ± 1.7 Ma, recording substantial postemplacement erosion. AHe results for the mafic to ultramafic Apishapa Dikes (oldest ~37 Ma, youngest ~14 Ma) located ~20–40 km farther north and east on the High Plains range from 12.0 ± 1.4 to 6.2 ± 1.9 Ma, documenting continued exhumation on the western High Plains during the ~12–5 Ma deposition of the Ogallala Formation farther east and suggesting that the western limit of Ogallala deposition was east of the Apishapa Dikes. In far southeastern Colorado, the Two Buttes lamprophyre was emplaced at 36.8 ± 0.4 Ma and yields a Late Oligocene AHe date of 27.1 ± 4 Ma. Here, the Ogallala Formation unconformably overlies Two Buttes, indicating that the regional ~12 Ma age for the base of the Ogallala is a minimum age for the exposure of the pluton at the surface. The AHe data presented here document that kilometer-scale erosion affected all of the southeastern Colorado High Plains in Oligo-Miocene time. While exhumation can have multiple possible causes, we favor contemporaneous surface uplift capable of elevating the region to modern heights.

The High Plains are a region of the Western United States that lies east of the Colorado and Wyoming Rocky Mountains and extends into Nebraska, Kansas, Oklahoma, and South Dakota (Figure 1). The High Plains of Colorado have surface elevations ranging from 1.3 to 2.1 km above mean sea level, even though the plains did not experience crustal shortening during the ca. 70–45 Ma Laramide Orogeny. When and how this vast, undeformed region was elevated is a major question in continental tectonics [1, 2]. Multiple, often conflicting, hypotheses have been proposed to explain the rise of the High Plains [2-6]. To deduce the mechanism behind how the High Plains reached their current height, we must first establish the time at which the plains reached these elevations [2]. Apatite (U-Th)/He (AHe) thermochronology provides a method to reconstruct the thermal histories of rocks, which can be used as a proxy for the erosion history of overlying material. Various mechanisms exist that could explain this erosion history [7], of which contemporaneous surface uplift is one. Different candidate surface uplift mechanisms should operate at different times, so that if AHe dates can be used to constrain when surface uplift occurred, some hypothesized mechanisms can be eliminated.

This thermochronologic approach has been used to deduce the surface uplift history of the Southern Rocky Mountains, which is similarly debated [6, 8, 9]. The Southern Rocky Mountains provide good targets for low-temperature thermochronologic studies due to abundant exposures of Precambrian crystalline basement and Laramide-age plutons, making such studies, using both AHe and apatite fission-track (AFT) thermochronology, possible there [8, 10-14]. Most thermochronologic dates from parts of the Southern Rocky Mountains that are not involved in the Oligo-Miocene Rio Grande Rift (RGR) are Laramide (~70–45 Ma) or pre-Laramide (pre-70 Ma). The abundant Laramide dates record rapid exhumation of the range in response to rock uplift. Clusters that give pre-Laramide dates reveal that some portions of the Southern Rocky Mountains experienced insufficient post-Laramide exhumation to expose rocks buried deeply enough to reset their dates [8, 9]. Exceptions to this data pattern are post-Laramide AFT and AHe dates from the San Juan and Elk ranges [6, 14] and post-Laramide AHe dates from a Late Cretaceous pluton in the Arkansas Hills [9].

In contrast, no low-temperature thermochronologic studies have been conducted on the Colorado High Plains with which to constrain the exhumation history and timing of uplift there. One likely explanation for the absence of such work is that the region’s exposures are dominated by Mesozoic sedimentary units in which high-quality apatite grains are scarce. Nevertheless, substantial post-Eocene exhumation on the High Plains is suggested by the pattern of the unconformity beneath sedimentary rocks of the ~12–5 Ma Ogallala Formation. This unconformity is underlain by northward-younging sedimentary rocks of Triassic age in northern and central New Mexico, Cretaceous age in southern Colorado, and Eocene age in northern Colorado [15, 16]. This pattern suggests a post-Eocene-pre-Ogallala (i.e., pre-12 Ma) exhumation event(s) that was most pronounced in southern Colorado and northern New Mexico. This exhumation has been attributed to a rock uplift event of diminishing northward magnitude [5, 17].

Here we use AHe dating of three clusters of apatite-bearing Cenozoic igneous intrusions in an ~200 km west-to-east transect across the Plains to investigate when exhumation across southeastern Colorado occurred (Figure 2). Apatite (U-Th)/He (AHe) thermochronology provides information about the thermal histories of rocks in the ~90–30°C temperature range, depending on accumulated radiation damage and grain size [18, 19], and can constrain the cooling history of rocks as they are exhumed through the upper ~1–3 km of the Earth’s crust. By targeting Cenozoic intrusions, we can more precisely resolve the timing of post-Laramide erosion of the sedimentary rocks that once covered the now exposed igneous bodies. We then use our thermochronologic results to decipher the broader spatial patterns of Cenozoic exhumation across southeastern Colorado, evaluate the implications for the original distribution of the Ogallala Formation in this region, and explore potential causes for this history.

Colorado’s Rocky Mountains and High Plains have a long and complex history of tectonic uplift, basin fill, and exhumation that influences the interpretation of our thermochronologic results. Here, we briefly summarize the key episodes in that history.

2.1. Ancestral Rocky Mountains and Synorogenic Fill

Proterozoic basement underpins the entire state. A thin (<500 m) passive-margin sedimentary sequence accumulated above the Great Unconformity during the early-middle Paleozoic. Uplift of the Ancestral Rockies during Pennsylvanian–Permian time stripped that cover sequence and removed some amount of Proterozoic crystalline basement from the ranges, whereas older Paleozoic rocks were preserved in the intervening basins and promptly buried by several kilometers of synorogenic sediment. In our study area, the Wet Mountains (Figures 1 and 2), which in places have Jurassic strata unconformably overlying basement, exemplify the location of an Ancestral Rocky Mountain range. The Sangre de Cristo Range, which consists mainly of a several-kilometer-thick Pennsylvanian–Permian sedimentary sequence belonging to the Sangre de Cristo Formation, illustrates the stratigraphy of an Ancestral Rockies basin. Several exposures of Sangre de Cristo Formation, or its lateral equivalents, exist in places where rivers of SE Colorado’s High Plains have incised deeply, indicating that this part of the plains was an Ancestral Rockies basin, not a range.

A thin Triassic–Jurassic cover sequence unconformably overlies both the Pennsylvanian–Permian synorogenic foreland basin deposits and the Precambrian basement that comprised those ranges, thus documenting that even the highest paleo-topography of the Ancestral Rocky Mountains had been eroded and buried by that time (Figure 3).

2.1.1. Western Interior Seaway Sediments

Southeastern Colorado’s High Plains, including our entire study area, preserve a thick accumulation of Cretaceous shales, sandstones, and limestones deposited in a marine setting during the transgression and subsequent regression of the Western Interior Seaway (WIS) [16, 20]. Rapid subsidence produced accumulations of between 1 and 2.8 km of shallow marine WIS sediment across Colorado, with the thickest accumulations in the northern part of the state and up to 2 km in the western portion of our study area, diminishing to less than 1 km in the eastern portion of the study area (Figure 3) [20, 21].

2.1.2. Laramide Orogeny

A second mountain-building event, the Laramide Orogeny, took place between ~70 and 45 Ma [22, 23] bringing much of the Precambrian material that is exposed in much of today’s Rocky Mountains to the surface. This includes the basement exposed in the Wet Mountains [8].

Colorado’s High Plains did not experience the same magnitude of uplift and deformation from Laramide tectonism as the Southern Rocky Mountains did. Instead, foreland basins developed on the western edge of what is now the High Plains in response to Laramide thrust loading. Cretaceous-Eocene synorogenic sediment accumulated in the Denver and Raton Basins (Figure 1). This synorogenic fill contains many Precambrian clasts, documenting that basement rock was exposed during the Laramide. Furthermore, detrital zircon ages from the Raton Basin formations have peaks at 1423–1430 Ma, and 1678–1687 Ma, all ages associated with crystalline basements from the nearby Wet Mountains [21, 24].

The westernmost samples in our AHe transect lie within the Raton Basin, where up to 3300 m of Late Cretaceous to Eocene arkosic sandstone and conglomerate units accumulated. These include the Poison Canyon, Cuchara, Huerfano, and Farasita Formations (Figure 3) [16, 21]. At the latitude of our transect, contours on the Dakota Formation show that it reaches a maximum elevation of ~1830 m at ~104° 40′ longitude. In the Raton Basin to the west, the Dakota Formation lies at significantly lower elevations, and to the east, out into the High Plains, it is only slightly lower [25]. This structural high defines the Laramide forebulge, thus delineating the eastern boundary of the Raton Basin. Transect samples east of 104° 40′ longitude (Figure 2) lie outside the Raton Basin, and thus are unlikely to have ever been covered by a substantial thickness of Laramide synorogenic fill (Figure 3).

2.1.3. Ignimbrite Flare-Up and Rocky Mountain Erosion Surface

Late Eocene erosion planed down much of the Laramide topography, creating a low-relief surface known as the Rocky Mountain Erosion Surface (RMES) [26].

Regional volcanism west of the study area, known as the ignimbrite flare-up, spread volcanics and volcaniclastics across the RMES and adjacent portions of the High Plains from the Late Eocene-Early Miocene [27]. The High Plains sedimentary units that compose this sequence are the White River Group [16] and their age equivalents, which are only exposed north of Colorado Springs (Figure 1) [17, 28]. They unconformably overlie both Denver Basin synorogenic and WIS sediments. Between Denver and Colorado Springs, the Wall Mountain Tuff, a 37-Ma ignimbrite flare-up ash flow [29], is interbedded between the Late Eocene Larkspur Conglomerate below and the Castle Rock Conglomerate above [9], both of which are WRE units on the High Plains [29]. The tuff also overlies the RMES in the Southern Rocky Mountains and has little or no displacement across the Southern Rocky Mountains-High Plains physiographic boundary [30]. This observation indicates that there was only gentle topography with local valleys cut into the RMES [31] across the Southern Rocky Mountains-High Plains boundary in Late Eocene.

No Eocene (WRE) material exists in our SE Colorado High Plains study area; the land surface consists everywhere of either older Raton Basin or WIS sediments or, from Two Buttes eastward, fluvial, Mio-Pliocene Ogallala Formation that lies unconformably atop the WIS strata (Figures 1 and 3), with an unconformity at the stratigraphic level of the White River. Are WRE strata missing because they were never deposited or because they were later eroded? Like previous authors [17], we strongly suspect they once existed and were later eroded away. The SE Colorado High Plains are closer to abundant ignimbrite flare-up volcanic sources than the White River Group of NE Colorado and adjacent states [28], so it is likely that WRE sediment once covered our study area. This supposition is supported by Smith et al. [32], who concluded that WRE strata exist in the Meade and Ashland basins of southwestern Kansas, where evaporite dissolution locally increased accommodation space [32].

2.1.4. Rio Grande Rifting and Emplacement of the Spanish Peaks

Around 28 Ma, extension began in Colorado and New Mexico. Normal faulting produced the RGR that was accompanied by alkaline magmatism [33, 34]. The RGR runs north–south through central New Mexico and Colorado, with the San Luis Valley, located immediately west of our study area, forming one segment of the RGR (Figure 1). Emplacement of the Spanish Peaks stocks and dikes between 27 and 22 Ma is associated with that magmatic episode [35, 36]. An elevation transect of Spanish Peaks intrusive rocks forms the western section of the AHe transect in this study (Figures 2 and 3).

The Jemez Lineament is a SW-NE trending series of volcanic fields that run through northern New Mexico and into eastern Colorado (Figure 1). Two Buttes, Colorado’s easternmost intrusive center and the easternmost sample in our transect, lies along the Jemez Lineament, but its emplacement age of ~37 Ma is considerably older than any other Jemez Lineament magmatic activity [37, 38].

2.1.5. Deposition of the Ogallala Formation

The Ogallala Formation in Kansas and NE Colorado began accumulating in the middle Miocene, by 12 Ma based on both faunal evidence and dated volcanic ash layers [39, 40]. An ash layer dated by K/Ar to 18 Ma was reported from Ellis County, Kansas, suggesting that Ogallala deposition in Kansas may have started as early as 18 Ma, but the abundant Upper Miocene fossils collected at the site indicate that this date is unreliable [40].

The Ogallala reaches a maximum thickness of 213 m [41] in Nebraska, and in Kansas, it ranges from 1 to 91 m thick [42]. The westernmost Ogallala exposure sits about 70 km east of the Front Range of the Southern Rocky Mountains. It is uncertain how far west of the current exposure the Ogallala used to extend, but the fact that a projection of the Ogallala outcrop base westward places it atop the RMES has led some workers to hypothesize that it once formed a continuous sheet that lapped onto the RMES at the range front [5, 33]. Others [32] argue that middle Cenozoic dynamic uplift and continental tilting caused the Ogallala to prograde eastward, resulting in a progressively younger onset of Ogallala deposition farther east.

Morgan et al. [43] conducted the only detrital zircon study to date on an Ogallala Formation outcrop from Colorado. The youngest single zircon grain in that sample is 23 ± 1 Ma and the youngest age population is 29 ± 1 Ma, providing the Ogallala Formation’s maximum depositional age; the actual age of deposition could be younger [32]. Smith et al. [32] point out that abundant sources of volcanically derived detrital zircons existed upwind of the High Plains throughout the Cenozoic, so they suggested that a sample’s youngest detrital zircon age population may closely resemble the actual depositional age in this region. Middle Miocene zircons are indeed found in Ogallala Formation outcrops from northern Kansas [39], so it is possible that the Ogallala in Colorado began accumulating in the Oligocene or Early Miocene rather than in the middle Miocene.

The Ogallala unconformably overlies progressively older units further south on the Colorado–New Mexico High Plains, from Eocene-Early Miocene in NE Colorado, to Cretaceous in SE Colorado, to Triassic in NE New Mexico [15, 16]. Roy et al. [17] concluded from this pattern that the High Plains of New Mexico and SE Colorado experienced a pre-Ogallala erosion episode of significantly higher magnitude than occurred in NE Colorado. This conclusion is consistent with Bogolub and Jones’ [2] finding that cryptic topography (which they define as topography that cannot be explained by observable processes and that must be supported in the subsurface) increases southwestward on Colorado’s High Plains. It is likely that whatever process created the cryptic topography is the same one that triggered the theorized exhumation episode.

3.1. Sampling Strategy

We targeted three clusters of Cenozoic igneous intrusions in an ~200 km-long east-to-west transect across the High Plains of southeastern Colorado for AHe thermochronology for four key reasons. (1) These lithologies contain high-quality euhedral apatite grains, as opposed to the apatite-poor sedimentary units that dominate this region. This means that grains have experienced similar thermal histories and can therefore be modeled together, unlike apatite grains that may be variably reset in a detrital dataset. (2) The fact that these Cenozoic intrusives were emplaced at some depth below the surface and are now exposed means that they underwent some magnitude of post-Eocene/post-Laramide cooling and erosion that we can potentially isolate using thermochronology. (3) These intrusives are likely to have undergone a cooling-only thermal history that is relatively simple to decipher with (U-Th)/He thermochronology, which contrasts with the more complex and variable thermal histories experienced by any detrital apatite grains in the exposed Cretaceous and Paleogene sedimentary rocks that may have been incompletely reset during burial [44]. (4) Finally, the west-to-east distribution of the samples across the southern Colorado High Plains provides the opportunity to evaluate spatial heterogeneity in the timing and magnitude of exhumation across this wide region.

In total we acquired AHe data for ten samples in three geographic groups: (1) the Spanish Peaks, (2) the Apishapa Dikes, and (3) the Two Buttes lamprophyre (Figures 1 and 2). Our westernmost study domain is the Spanish Peaks, consisting of Oligocene plutons that intrude Cretaceous and Paleogene sedimentary rocks. We collected five samples from these exposed granodiorite plutons and associated dikes to generate a ~1600-m elevation transect from 3867 m to 2272 m across a lateral distance of ~15 km (Figure 4). The top two samples of this transect are part of the East Spanish Peak pluton emplaced at 23.9 ± 0.08 Ma (Figure 5) (40Ar/39Ar dates from Penn and Lindsey [35]), while the three lower elevation samples are from radial dikes, which were emplaced between 21.9 ± 0.13 and 23.3 ± 0.04 Ma [35]. The central domain consists of four samples of the Apishapa Dikes (Cross, 1915) collected across a lateral distance of 60 km and encompassing an elevation range from 2054 to 1731 m (Figure 2). From west to east these samples are the 36.6 ± 0.12 Ma Little Black Hills), the 25.2 ± 0.25 Ma Huerfano Butte, the 14.3 ± 0.1 Ma Cuchara Reservoir Dike, and the undated Mica Butte (all emplacement ages are 40Ar/39Ar dates [35]). These isolated dikes are characterized by 10’s of meters of vertical relief and intrude into the Cretaceous WIS sedimentary sequence (Figure 3). The easternmost domain consists of a single sample of Two Buttes lamprophyre located ~200 km east of the Rocky Mountain front at an elevation of 1326 m. Two Buttes intrudes upper Paleozoic to Jurassic sedimentary rocks and is unconformably overlain by the Miocene Ogallala Formation (Figure 2). The lamprophyre has been K-Ar dated at 36.9 ± 0.4 Ma [37].

We additionally acquired zircon (U-Th)/He (ZHe) dates from the Two Buttes sample to confirm its emplacement age. The ZHe results should date lamprophyre emplacement because the stratigraphy (Figure 3) suggests that the lamprophyre was emplaced at relatively shallow depths, cooler than the ~200°C temperature sensitivity [45] of He diffusion in low-damage zircon.

3.2. Analytical Methods

Samples were prepared and analyzed at the University of Colorado Thermochronology Research and Instrumentation Laboratory (CU TRaIL). Whole rock samples were crushed to <500 µm size, then separated via water density, magnetic, and heavy liquid (lithium metatungstate) methods. For each sample, 4–9 apatite crystals were hand-selected to be euhedral and inclusion-free, with a minimum dimension larger than 60 µm. Crystal selection was done using a Leica M165 binocular microscope equipped with both reflected and transmitted polarized light. Two zircon grains were additionally selected from the Two Buttes sample for analysis. Dimensions of apatite and zircon grains were measured with a calibrated digital camera prior to packaging in Niobium tubes. Helium was degassed and analyzed on an ASI Alphachron He extraction and measurement line. U, Th, and Sm were measured using either a Thermo Element 2 ICP-MS or an Agilent 7900 ICP-MS. Detailed apatite and zircon analytical methods for He degassing and measurement, spiking and dissolution, and parent isotope analysis followed those described in Peak et al. [46]. Grain masses, nuclide concentrations, alpha-ejection corrections, and associated data were calculated using the methods described in Ketcham et al., [47]. AHe dates were calculated and uncertainties were propagated using the HeCalc python program [48]. Several aliquots were discarded from the dataset as they were deemed poor quality due to near-zero levels of parent U, Th, and Sm, and therefore do not produce usable (U-Th)/He dates.

3.3. Time–Temperature (tT) and Time–Depth (tZ) Modeling Methods

We carried out thermal history simulations to quantitatively constrain the timing and magnitude of postemplacement cooling and exhumation of the Spanish Peaks, Apishapa Dikes, and Two Buttes samples. Models were carried out using inverse modeling in HeFTy program v.2.0.9 [49] and employing the Radiation Damage Accumulation and Annealing Model (RDAAM) for apatite He diffusion [19]. For each sample, data were input using a standard approach in which individual analyses were binned into synthetic grains based on their effective uranium concentrations (eU), with the average equivalent spherical radius, eU, and AHe date of each synthetic grain input into HeFTy [50, 51]. AHe date uncertainties were either the standard deviation of the binned dates or 10%, whichever is larger. For each inverse model, either 10,000 or 100,000 potential tT paths were tested, the predicted dates compared against the input data, with “good” and “acceptable” paths that meet a statistical goodness of fit criterion being output and shown in the tT plots (Figures 6-8). We tested 10,000 paths for all single-sample models and 100,000 paths for the Spanish Peaks multi-sample models. Testing more paths was necessary for the multi-sample models because the additional data richness allows the model to reject a larger number of possible tT paths; these more stringent fitting criteria have the benefit of producing a tightly constrained model, but because the model accepts fewer good fits, more potential paths must be run to produce a representative suite of viable tT paths.

All samples were run in tZ mode, which uses a 1D thermal model to convert tZ to tT paths assuming a magnitude of topographic development and tilting, a geothermal gradient, initial sample depth, and other parameters [52]. We assumed full topographic development and no tilting over the course of the model, which adopts a flat surface above all samples at emplacement and no tilting during exhumation. We based the assumption of full topographic development on the widespread presence in the region of the Eocene RMES and the lack of displacement of the Wall Mountain Tuff across the Rocky Mountains-High Plains boundary (as discussed in the geologic background section), which suggests that little relief had developed prior to the ignimbrite flare-up [30]. The absence of tilting during exhumation is based on the subhorizontal layers of Paleogene-aged Cuchara and Huerfano Formations that form a contact aureole around the West Spanish Peak intrusion (Figure 3).

All tZ models assume surface temperatures based on the modern sample elevation assuming a sea level temperature of 20°C and an atmospheric lapse rate of 6°C/km. For the Spanish Peaks multisample models we carried out multiple trials with variations in (1) geothermal gradient (online Supplementary Material Figure S1), and (2) initial sample depth (online Supplementary Material Figure S3). For the easternmost Apishapa Dikes samples we performed multiple trials that tested both cooling-only and Ogallala reburial thermal histories (online Supplementary Material Figure S4). We discuss below the results of these multiple trials and in each case show our preferred model in the text and alternative models in the supplement.

All five samples from the East Spanish Peak elevation transect were modeled simultaneously in multi-sample tZ mode. The setup and outcomes of each model are discussed below, with additional model details in online Supplementary Material Table S1.

We report 51 new AHe analyses for ten samples across the High Plains of southeastern Colorado, as well as 2 new ZHe analyses for the Two Buttes sample. The data are reported in Table 1, with extended data details reported in online Supplementary Material Table S2. These tables follow the reporting protocols of Flowers et al., [50] [53]. Five apatite analyses that are reported in Table 1 are excluded from data plots because they appear to be outliers. Two of these grains are older than the crystallization age of the unit, two others are anomalously old compared with the rest of the sample’s data distribution, and one is anomalously young compared to the rest of the grains in the sample. A bias to older dates can occur because some difficult-to-detect factors such as micro-inclusions and He implantation yield excess He [50]. Uncertainties on single-grain dates are reported and plotted at 2-sigma and include the propagated analytical uncertainties on parent and daughter isotope measurements. All eU uncertainties are estimated at 15% of the eU value. Owing to the good reproducibility of dates within individual samples, we report mean sample dates for all samples with uncertainties reported as the standard deviation of single-grain aliquots. AHe date-eU plots for all samples are in Figures 6-8. Such plots are a common way to evaluate potential radiation damage effects on He retentivity and visually inspect data reproducibility. All samples in this study show uniform dates across the analyzed span of apatite eU.

In the west, the East Spanish Peak elevation transect (Figure 5) consists of 26 single-grain AHe dates from five samples across an elevation range of 1600 m. The AHe dates decrease from 17.3 ± 1.3 Ma at the mountain’s peak (3867 m elevation) to 14.1 ± 1.7 Ma from a radial dike at the base of the transect (2272 m elevation). All mean sample AHe dates are younger than 40Ar/39Ar crystallization age of their respective intrusion [35]. For example, the 17.3 ± 1.3 Ma AHe date of the summit sample and the 18.8 ± 1.4 Ma AHe date of sample CS13-7 are younger than their emplacement ages of 23.9 ± 0.08 Ma. The three lowest elevation samples have mean AHe dates from 16.9 ± 1.2 Ma to 14.1 ± 1.7 Ma, younger than the emplacement ages of these radial dikes at 21.9 ± 0.13 Ma to 23.3 ±0.04 Ma.

The Apishapa Dikes sample suite, which occupies the middle portion of the west-to-east transect, includes 19 AHe analyses for four samples across elevations of 2054 m to 1731 m over a lateral distance of 60 km. All Apishapa samples reside at lower elevations and yield dates younger than the lowest elevation sample of the East Spanish Peak transect. The Apishapa samples also all yield mean sample AHe dates younger than the emplacement age range of this dike suite, which is from 36.6 to 14.3 Ma [35]. The 36.6 ± 0.12 Ma Little Black Hills sample (CS13-1), ~30 km north of the Spanish Peaks, yields an AHe date of 11 ± 1 Ma. Farther east, the 25.2 ± 0.25 Ma Huerfano Butte lamprophyre (LA20-1) yields an AHe date of 7.7 ± 0.7 Ma. The 14.3 ± 0.1 Ma Cuchara Reservoir Dike (LA20-6) is AHe dated to 8.2 ± 1.3 Ma. Finally, the Mica Butte sample (SF21-5) has an AHe date of 6.2 ± 1.9 Ma. While the emplacement of the Mica Butte dike has not been directly dated, all dikes associated with the Spanish Peaks plutons have been dated in the 25–14 Ma age range [35], making it likely that Mica Butte was also intruded during this time period.

In the east, four apatite analyses for the Two Buttes (SK21-6) sample yield a mean AHe date of 27.1 ± 4.2 Ma, older than all samples from the Spanish Peaks and Apishapa Dikes. The AHe date post-dates the 36.9 ± 0.4 Ma [37] K-Ar (i.e., crystallization) age of this sample. The ZHe data yield a date of 36.4 ± 3.0 Ma, in agreement with the K-Ar data and consistent with the emplacement of the lamprophyre at temperatures cooler than the ~180–200°C closure temperature of undamaged zircon.

The fact that the AHe dates are younger than intrusion ages for all our samples across the entire 200 km transect (Figure 2) suggests that SE Colorado’s High Plains have experienced substantial exhumation during the Oligo-Miocene. Here we describe our inverse thermal history modeling of samples in each portion of the transect to quantitatively test this interpretation, summarize the constraints those models impose on the timing and magnitude of Oligo-Miocene exhumation, and discuss the geologic implications for the former extent and thicknesses of specific sedimentary units that have been denuded from the landscape since 24 Ma.

5.1. Western Portion of Transect—The Spanish Peaks Elevation Profile

5.1.1. Timing and Magnitude of Exhumation at Spanish Peaks

The AHe data for the five samples from the ~1600 m elevation profile from East Spanish Peaks are consistent with substantial and rapid cooling and exhumation in Early to Middle Miocene time. This is indicated by the consistency of the Miocene dates over the entire profile (18.8 ± 1.4 Ma to 14.1 ± 1.7 Ma, overall younger toward lower elevations; Figure 5) and the uniformity of AHe single-grain dates within individual samples across a broad eU range (Figure 6).

To use these data to quantitatively constrain the timing and magnitude of exhumation after ~24 Ma emplacement of the intrusives, we carried out inverse multi-sample tZ modeling in HeFTy. Multi-sample models take all five samples in the transect into account, testing pathways that honor the cooling histories of the adjacent samples and assuming all samples have remained at the same depth relative to one another through their history [52]. Our models used the highest elevation sample (LA19-9) from the summit of East Spanish Peak as the control sample [52]. We tested the role of both the geothermal gradient and the summit sample emplacement depth on the model outcomes. Besides the emplacement depth of the highest elevation sample, the only other constraint imposed on all samples was surface conditions at present day. Our preferred model is in Figure 6, with variants in online Supplementary Material Figures S1 and S2.

Geothermal gradients of 40, 45, 50, 60, and 70°C/km were tested in our multi-sample tZ models. These values were selected based on modern local Spanish Peaks geotherms of 40–50°C/km [54, 55] the interpretation that the Oligocene geotherm at Greenhorn Mountain in the Wet Mountains (Figure 2) was as high as 47°C/km [8], and the observation that low elevation portions of the Raton Basin southeast of the Spanish Peaks have modern geotherms >70°C/km [55]. Morgan [55] concluded that these elevated Raton Basin geotherms are due to west-to-east, down-dip groundwater flow in permeable units (i.e., away from the Spanish Peaks). We ultimately use 40°C/km in our preferred model (Figure 6) because the simulation using a geotherm of 50°C/km had very few good fit paths (online Supplementary Material Figure S1), the 60°C/km model had some acceptable fits but no good fits (online Supplementary Material Figure S2), and the 70°C/km geotherm produced no fits at all. Notably, all models with elevated geotherms (online Supplementary Material Figures S1 and S2) produce tT paths that do not differ appreciably from those of our preferred model (Figure 6). This sensitivity testing suggests that (1) the Spanish Peaks area did not experience a sustained geotherm of more than 40–50°C/km; (2) even if an elevated geotherm did exist, it does not change our fundamental conclusions.

For the emplacement depth estimates, the phaneritic texture of all East Spanish Peak samples points toward at least 1–2 km of overburden covering the modern-day East Spanish Peak summit at 24 Ma during pluton emplacement [34]. We therefore used a 0.5–2 km emplacement depth for our preferred multi-sample model as the most conservative exhumation-magnitude scenario (Figure 6). However, a deeper emplacement depth is not precluded by the data, so we additionally tested the effect of emplacement depth choice on the model results by running a comparable model with a greater emplacement depth of 2–3 km (online Supplementary Material Figure S3). Emplacement depths >3 km are unlikely based on consideration of the regional stratigraphy that may have been eroded (see discussion below). The model results using deeper emplacement depths (online Supplementary Material Figure S3) produce essentially the same range of good fits as our preferred model, so our results are insensitive to our choice of this parameter.

The good-fit paths for our preferred model (40°C/km geotherm; 1–2 km emplacement depth; Figure 6) show exhumation beginning sometime between pluton emplacement (24 Ma) and ~17 Ma, with the exhumation rate accelerating sometime between 21 and 17 Ma. The East Spanish Peak summit sample likely was not emplaced shallower than 1.7 km, because no good fit paths begin at depths shallower than that. That same sample was within 250 m of the surface by 16 Ma. These model results are consistent with a minimum of 1.45 km of exhumation at the summit of Spanish Peaks (i.e., 1.7 km - 0.25 km = 1.45 km) between 24 and 16 Ma, a rate exceeding 180 m/m.y. All samples show that the episode of rapid exhumation slowed considerably sometime between 18 and 16 Ma, with all of them subsequently cooling at a uniform rate. The modeled exhumation rates are consistent with the qualitative analysis of the age-elevation profile in Figure 5. The deepest sample (SK20-1) was at ~60°C at 18–16 Ma; given our assumed 40°C/km geotherm, it was located at ~1.5 km depth at this time and has moved toward the surface due to a steady ~100 m/m.y. exhumation rate ever since. Our model results suggest that the total amount of post-24 Ma exhumation at Spanish Peaks for the lowest elevation sample has been a minimum of 2.85 km (1.25 km exhumation of LA19-9 + 1.6 km vertical difference between samples).

5.1.2. Former Extent and Thicknesses of High Plains Sedimentary Units at the Spanish Peaks

We can use the regional stratigraphy (Figure 3), combined with the exhumation magnitude constraints from the AHe data and associated modeling (Figure 6), to infer the extent, nature, and thicknesses of the sedimentary units that previously buried the landscape at Spanish Peaks. The lower portions of the East Spanish Peak pluton are surrounded by contact metamorphosed Cuchara Formation [25], part of the Raton Basin syn-Laramide orogenic fill (Figure 3). Erosion has removed all the sedimentary rock that once surrounded the upper portion of that pluton, but the upper portions of the adjacent West Spanish Peak pluton (Figure 4) remain surrounded by contact metamorphosed, subhorizontal Huerfano Formation [25], demonstrating that the upper portions of both East and West Spanish Peak plutons were intruded at that stratigraphic level. The maximum exposed thickness in the Raton Basin of the synorogenic Huerfano and overlying Farasita formations is less than 1 km (Figure 3) [21]. That fact, combined with the minimum pluton emplacement depth constraint of 1.7 km from our HeFTy modeling requires that at least several hundred meters of material that is completely missing from today’s Raton Basin covered the Spanish Peaks at the time of pluton emplacement. Given the geologic history of the area discussed earlier in this paper, the only plausible candidates for that missing material are a once thicker Laramide synorogenic sequence, ignimbrite flare-up associated White River Equivalent (WRE) volcanics and volcaniclastics, and/or volcanic and debris flows from the Spanish Peaks (though few of the plutons display evidence that they ever vented to the surface [35].

The maximum total thickness of extant Raton Basin synorogenic deposits is ~3.3 km (i.e., the Poison Canyon through Farasita Formations in Figure 3) [21]. This thickness is comparable to the >2 km of synorogenic fill in the structurally similar Powder River Basin of Wyoming [56] and considerably thicker than the ~0.9 km thickness in the Denver Basin [57]. Therefore, it is possible but seemingly unlikely that the Raton Basin once contained a substantially thicker Laramide synorogenic fill than is currently exposed. We consider it much more likely that the now-eroded overburden that overlaid the summit of East Spanish Peak during pluton emplacement consisted mainly of WRE material. A 200-m thick 33.4 Ma ash flow [29] caps the summit of Greenhorn Mountain (Figure 2) [25], about 50 km north of the Spanish Peaks. As noted above, the Raton Basin lies close to the flare-up volcanic sources of the White River Group sediments, so it likely once hosted a thicker WRE stratigraphic section than the 300 m currently exposed in NE Colorado, South Dakota, and adjacent areas [58].

The total thickness of any preexisting WRE material that once overlaid the Spanish Peaks is unknown, so we cannot constrain the maximum burial depth of the pluton with confidence using just its stratigraphic position. However, based on the preceding arguments, a burial depth of 2 km would likely require a little over 1 km of now-eroded WRE material. Deeper emplacement depths, such as that modeled in online Supplementary Material Figure S3, would require abundantly thicker WRE sections. The High Plains generally lack accommodation space, making it unlikely that multiple kilometers of WRE ever accumulated there. When Cather et al. (2012) projected the gentle east dip of the Ogallala Formation from its westernmost outcrop on the Colorado High Plains near Limon (Figure 1) westward onto the RMES, they estimated that about 800 m of pre-Ogallala material has been eroded from beneath the basal Ogallala unconformity to the present land surface that consists of Cretaceous WIS sediments at the western limit of the High Plains [5]. Based on these considerations, we think it unlikely that more than 1–2 km of WRE sediments ever accumulated above the Spanish Peaks, hence our choice of a 1–2 km burial depth for the control sample LA19-9 in our favored HeFTy model (Figure 6). However, Cather et al.’s [5] analysis assumes that the Ogallala Formation once extended westward all the way to the Rockies, which our Apishapa Dikes data calls into question for SE Colorado, and in any case it applies only to post-Ogallala erosion; if a considerable thickness of strata was removed from the entire High Plains prior to deposition of the Ogallala Formation, that analysis would not serve to constrain the magnitude of WRE removal.

5.2. Middle Portion of Transect—The Apishapa Dikes

5.2.1. Timing and Magnitude of Exhumation at the Apishapa Dikes

The four Apishapa Dikes samples are at lower elevations than all samples at Spanish Peaks (~540 m lower for the lowest elevation sample) and yield younger dates than the AHe data from the Spanish Peaks. While the Apishapa Dikes are not part of an elevation transect because of the lateral distance from the Spanish Peaks and one another, samples continue the same trend of younger AHe dates with decreasing elevation as seen at the Spanish Peaks. This pattern provides evidence that the Early to Middle Miocene rapid exhumation event at Spanish Peaks also encompassed the Apishapa Dikes and thus was at least subregional in extent, with exhumation continuing until at least ~6 Ma (the AHe date of the youngest sample).

To better quantify this exhumation signal, we conducted tZ modeling of the individual Apishapa samples (Figure 7). We did not simulate the samples in a multi-sample model because the samples have different crystallization ages, which makes it challenging to appropriately set up such a model in HeFTy, and because their collection across a ~60-km distance suggests that a multi-sample model may not be well-justified. We assumed a geothermal gradient of 40°C/km and emplacement depths of 0.5–2 km for all sample models, using the same rationale as for our preferred model for Spanish Peaks, as well as additional evidence for the emplacement depth as outlined below. All models begin at the crystallization age of each sample, which varies for the different dikes, and ends with surface conditions at present-day.

We additionally carried out models for the two eastern Apishapa samples that explored whether the Ogallala was ever deposited in the area (online Supplementary Material Figure S4). This test was accomplished by allowing exhumation of the samples to the surface by 12–9 Ma, and then permitting (but not requiring) burial by up to 500 m of Ogallala. It is unlikely that the Ogallala was ever thicker than ~500 m because, as noted in the geologic background, accommodation space on the High Plains is limited and the maximum preserved thickness of the Ogallala is about 200 m [41]. These models yielded no good fits for the easternmost sample (Mica Butte, online Supplementary Material Figure S4(a)) and very few good fits for the other sample (Cucharas Reservoir Dike, online Supplementary Material Figure S4(b)), consistent with the notion that no substantial Ogallala deposition ever occurred here, as discussed further below.

The good-fit paths for the Apishapa models without Ogallala burial suggest the onset of rapid exhumation was recorded in the AHe data starting at ~13 Ma for the westernmost Apishapa sample (Little Black Hills), and as late as ~7 Ma for the easternmost and lowest elevation Apishapa sample (Mica Butte) (Figure 7). These model outcomes suggest that the exhumation event we documented at the Spanish Peaks (Figure 6) also affected an area of the High Plains >50 km farther east and continued until at least 7 Ma (Figure 7). Our models do not constrain whether the bulk of Apishapa Dike exhumation had ended by ~7 Ma or if it continued later, because good-fit tT paths exist for all samples that both place them on the surface by ~7 Ma and that include a component of significant post-7 Ma exhumation (Figure 7).

5.2.2. Former Extent and Thicknesses of High Plains Sedimentary Units at the Apishapa Dikes

We can use the stratigraphic level of each Apishapa Dike intrusion to analyze the identity of the sedimentary section that once overlaid each dike, in the same way we did for the Spanish Peaks. The Little Black Hills (CS13-1) and Huerfano Butte (LA20-1) samples were both intruded into the Pierre Shale unit of the WIS sequence (Figure 3), suggesting that between about 0.8 and 1.5 km of WIS sediment once overlaid the Little Black Hills and Huerfano Butte samples (using the maximum values for WIS unit thicknesses from [21]).

The Little Black Hills lie in the heart of the Raton Basin while Huerfano Butte lies on its periphery, so either sample could have once been covered by up to 3 km of synorogenic sediment. However, the Little Black Hills lie <30 km south of the 33.4 Ma Greenhorn Mountain ash flow (Figure 2; [29]) and Huerfano Butte lies <20 km east of it. The ash flow unconformably overlies Precambrian crystalline basement rock of the Wet Mountains at 3550 m elevation. That unconformity marks the 33 Ma land surface, which stood ~1.7 km above Huerfano Butte and ~1.5 km above the Little Black Hills. As noted previously, farther north in the Denver Basin the 37.6 Ma Wall Mountain Tuff unconformably overlies both the Precambrian basement and Denver Basin synorogenic fill at the Rocky Mountains–High Plains interface with little or no topographic offset [30]; by analogy, the High Plains topographic surface at 33 Ma likely stood ~1.7 km above Huerfano Butte and ~1.5 km above the Little Black Hills before both were buried by an unknown thickness of WRE volcanic and volcaniclastic material. By the time Huerfano Butte intruded at 24 Ma [35], at least 200 m of WRE had accumulated, based on its extant thickness at Greenhorn Mountain, placing the land surface >1.9 km above the sample height.

Electron microprobe analysis of equilibrium mineral assemblages in the Huerfano Butte lamprophyre suggests that the dike was emplaced at a pressure of ~500–700 bars (Shelby Litton and Aaron Bell, personal communication, 2020), which translates to an emplacement depth of ~2.0–2.9 km assuming an average overburden density of 2500 kg/m3. If this emplacement depth estimate is correct, that would indicate that no more than ~1.0–1.2 km of WRE material covered Huerfano Butte when the dike was emplaced. This thickness estimate is consistent with the independent estimate we made above for the Spanish Peaks.

The Cuchara Reservoir Dike (LA20-6) and Mica Butte (SF21-5) samples were both intruded at the stratigraphic level of the Benton Formation, implying that between 0.7 and 2 km of WIS sediment could have once covered both samples (Figure 3; using the minimum and maximum values for WIS unit thicknesses [21]). The HeFTy model for the Cuchara Reservoir Dike (LA20-6, Figure 7(c)) suggests that the dike was emplaced at least 1 km below the surface, so at least some and possibly all of the overlying WIS strata remained in the Apishapa Dikes area until at least 14 Ma. Both these samples lie east of the forebulge that marks the eastern boundary of the Raton Basin (based on structural contours drawn on the Dakota Formation as discussed previously in this paper [2, 25]), so it is unlikely that any substantial thickness of Raton Basin synorogenic sediment ever accumulated here. Any additional overburden that might have once overlaid these samples likely consisted of WRE material.

5.3. Eastern Portion of the Transect—Two Buttes

5.3.1. Timing and Magnitude of Exhumation at Two Buttes

The Two Buttes lamprophyre has an Oligocene AHe date of 27.1 ± 4.2 Ma (Table 1) which is younger than its Eocene emplacement age constrained both by previous K/Ar data (36.8 ± 0.4; [37]) and our new ZHe date (36.4 ± 3.0 Ma). This outcome indicates some magnitude of cooling and exhumation after ~37 Ma, even at this location 200 km east of the Rocky Mountain front.

To quantitatively evaluate the timing and magnitude of this exhumation signal we ran a tZ model for this sample (Figure 8). We assumed a geothermal gradient of 30°C/km based on the modern local geotherm in the area. This value is cooler than the geotherm assumed for our Spanish Peaks and Apishapa Dikes simulations, because today’s geotherm decreases systematically across southern Colorado’s High Plains, from 40°C/km in the west near Spanish Peaks to 20°C/km in the southeastern corner of the state. Our preferred model also assumes a 0.5–2 km emplacement depth, the same as assumed for the Spanish Peaks and the Apishapa Dikes. This depth estimate is calculated from the observation that the lamprophyre intruded into the Morrison Formation such that ~1 km of WIS sediment used to overlie it (Figure 3; [20]), although a greater emplacement depth due to burial by post-WIS strata cannot be precluded. The model begins at the Two Buttes crystallization age, assumes that surface conditions were reached by 12–9 Ma because Two Buttes is unconformably overlain by Ogallala sediments of this age [32], and allows (but does not require) up to 500 m of Ogallala deposition before returning to surface conditions by present-day.

For our preferred model (30°C/km; 0.5–2 km emplacement depth), no good-fit paths are yielded for emplacement depths shallower than ~1.7 km, indicating that is the minimum emplacement depth for Two Buttes and that 1. 7 km is the minimum amount of post-37 Ma exhumation. The good-fit paths indicate that exhumation began between ~29 and 24 Ma, distinctly after emplacement. Numerous good fit paths allow some magnitude of Ogallala burial up to 500 m depth, thus showing consistency with the presence of Ogallala on the landscape at this location.

As noted above, the youngest detrital zircon grain yet found in the Ogallala of Colorado is 23 ± 1 Ma [43], meaning it is possible that the Ogallala began accumulating in Colorado during the Early Miocene rather than at 12–9 Ma as our model assumes (Figure 8). Note that several good-fit tT paths in Figure 8 are consistent with Early Miocene onset of Ogallala deposition, as they place the Two Buttes lamprophyre at or near the surface by ~23 Ma. Nevertheless, we chose to assume as our model constraint the youngest reasonable age of onset for Ogallala deposition at Two Buttes (based on its youngest age farther east in Kansas). We did this because it produces the most conservative possible interpretation of when and how fast exhumation occurred on the High Plains of SE Colorado. If Ogallala deposition began at ~23 Ma, it does not change the requirement that kilometer-scale exhumation at Two Buttes began sometime between ~29 and 24 Ma; it does constrain more tightly the timing of that exhumation (to ~29–23 Ma compared with the ~29 Ma to as recent as 12 Ma episode shown in Figure 8) and it requires an exhumation rate in excess of 280 m/m.y. (>1.7 km in less than 6 m.y.), faster than that demanded by our conservative estimate of the maximum Ogallala depositional age.

5.3.2. Former Extent and Thicknesses of High Plains Sedimentary Units at Two Buttes

Given that Two Buttes must have been emplaced deeper than 1. 7 km and it is estimated that the WIS section in this part of Colorado was ~1 km thick [20], it is likely that post-WIS strata once overlay Two Buttes. It is equally likely that those strata consisted of WRE volcanic/volcaniclastic material. There are three reasons for this supposition: (1) Two Buttes lies far to the east of the Raton Basin, so it is unlikely that any substantial synorogenic sediment was deposited here (Figure 3); (2) Two Buttes was emplaced during the ignimbrite flare-up (it is by far the easternmost manifestation of that event in Colorado); it might have vented to the surface and constructed a local volcanic pile; (3) Two Buttes is unconformably overlain by the Ogallala Formation, which requires that the lamprophyre had reached the surface by the time Ogallala sediments began to accumulate, meaning the missing overburden must predate Ogallala deposition.

The Ogallala Formation was accumulating on portions of the High Plains east of Limon (Figure 1), including at Two Buttes, between ~12 Ma (or possibly earlier) and 5 Ma. The Apishapa Dikes have ~12–6 Ma AHe dates that coincide with this depositional interval, suggesting that erosion was taking place in the Apishapa Dikes area simultaneously with Ogallala deposition farther east. This observation suggests that the Ogallala may not have been deposited much farther west in Colorado than its current westernmost outcrop extent (Figure 1). That hypothesis conflicts with the conventional wisdom that the Ogallala Formation was deposited up to the edge of the Rocky Mountains and onlapped the RMES in Colorado [5, 33] in the same way it does today in southern Wyoming [59].

In order to test this hypothesis quantitatively, we ran HeFTy models for the Mica Butte and Cuchara Reservoir dikes that are identical to the models shown in Figure 7 except that they require the intrusions to reach the surface by 12–9 Ma and then allow up to 500 m of reburial by Ogallala (or younger) sediment (online Supplementary Material Figure S4). The Mica Butte model that assumed rapid exhumation to the surface and allowed reburial by the Ogallala was unable to find any good fits to the AHe data (online Supplementary Material Figure S4(a)). This is the easternmost Apishapa Dikes sample, closest to existing Ogallala outcrops, so it is the one that is most likely to have once been covered by Ogallala Formation sediment; the fact that we cannot obtain good model fits for this sample strongly suggests that Ogallala Formation was not deposited here. The simulation of the Cuchara Reservoir Dike that allowed Ogallala reburial obtained three good fits (online Supplementary Material Figure S4(b)), but in order to match the sample’s AHe dates, all good fit paths require the sample to exhume rapidly from >1.5 km depth between 10 and 9 Ma, be reburied by >350 m of Ogallala, and stay buried until it is re-exhumed sometime between 1 and 2 Ma. There are two implausible aspects to such a scenario. First, it prohibits significant exhumation of the dike until after 10 Ma but then requires a furious pulse of exhumation (>1.5 km/m.y.) to ensue between 10 and 9 Ma. Second, that rapid exhumation pulse must be followed immediately by rapid reburial; it is hard to see how the transition from extreme exhumation to deposition could occur so rapidly on the low-relief High Plains.

The AHe data combined with these model results suggest it is highly unlikely that any lasting and substantial accumulation of the Ogallala Formation was ever deposited as far west as the Apishapa Dikes in SE Colorado. It is possible that a thin veneer of Ogallala material was temporarily deposited here and then promptly re-eroded and transported farther east (as posited by Smith et al. [32]); our models would not be able to detect deposition of such a transient Ogallala veneer. However, we consider it more likely that Spanish Peaks magmatism generated a paleotopographic high relative to the surrounding landscape, such that no Ogallala was ever deposited in this region.

Our AHe data and modeling results document that a rapid pulse of exhumation began on SE Colorado’s High Plains, stretching at least 200 km from the Spanish Peaks eastward to Two Buttes, sometime between 27 and 19 Ma and removed >1.7 km of strata from the Spanish Peaks area by 16 Ma. The exhumation rate slowed at Spanish Peaks after ~17–16 Ma and it had ended by ~12–9 Ma (and possibly as early as ~23 Ma) at Two Buttes, as revealed by the Ogallala Formation unconformably overlying the Two Buttes intrusion. But exhumation continued in the Apishapa Dikes area, removing >1 km of additional strata from this area at ~9–7 Ma or thereafter. Figure 9 schematically summarizes the geologic events that transpired on SE Colorado’s High Plains from the latest Oligocene to the present.

Although previous workers have proposed middle Cenozoic exhumation of the SE Colorado High Plains based on stratigraphic and geologic observations [5, 34], no AFT data from this region or the nearby Wet Mountains are younger than Oligocene (i.e., none are younger than 23 Ma; [8, 17, 60]), and no AHe data have been previously published from SE Colorado’s High Plains, so the exact timing of this exhumation event has been unknown. Our AHe data are distinctly younger than all previously published AFT data from the High Plains and show that major exhumation occurred during the latest Oligocene to mid-Miocene.

Looking northward, no comparable low-temperature thermochronologic studies have yet been conducted on the High Plains of NE Colorado (north of Limon, Figure 1). However, the widespread presence in that region of Eocene-age sediments belonging to the White River Group and its age equivalents (i.e., Castle Rock Conglomerate) strongly suggests that the magnitude of Cenozoic exhumation there is much less than it has been in SE Colorado. Some geologic processess that triggered extensive Cenozoic erosion on SE Colorado’s High Plains were either absent or less active on NE Colorado’s High Plains. A variety of recent studies have offered models to explain why the High Plains are so anomalously high; all such models need to reconcile this fundamental discrepancy between the Cenozoic erosional histories of these two segments of the High Plains, which in other respects seem to be nearly identical.

7.1. Implications for Timing and Causes of Surface Uplift

We favor exhumation induced by Oligo-Miocene surface uplift as the primary explanation for the kilometer-scale Miocene exhumation that we document on SE Colorado’s High Plains. However, alternative explanations for a cooling signal in our AHe data include geotherm relaxation or exhumation induced by climate change or drainage reorganization. We first discuss these options, before outlining our arguments for Oligo-Miocene elevation gain.

A previous AFT study proposed that the southern Wet Mountains near our western Apishapa Dike samples experienced an Oligocene geothermal gradient of 47°C/km, much higher than areas to the north, and potentially marking the northern edge of a regional heat flow anomaly extending into northern New Mexico caused by Oligocene volcanism [8]. Relaxation of such an initially high geothermal gradient would cause shallow isotherms to migrate to a greater depth, locking in the AHe date when the cooling front passed downward rather than when the sample cooled as it was advected upward through the crust during erosional exhumation.

However, thermal relaxation cannot be the sole explanation for our data because some exhumation to expose our intrusive samples at the surface following emplacement is undeniable. Moreover, the hypothesized 47°C/km Oligocene geothermal gradient [8] is comparable to today’s 40–45°C/km geothermal gradient in the same area [54]. Although geothermal gradients are unlikely to remain perfectly static over 25 m.y., we adopt the most conservative interpretation of a constant geotherm over this time frame. If the geotherm did relax, the exhumation magnitudes required to bring the intrusive rocks to the surface would be even larger than those inferred from our models but would not otherwise affect our conclusions. We therefore consider major relaxation of the geotherm (>10°C/km) to not play a substantial role in the history recorded by our AHe data, and instead interpret our results in the context of exhumation.

By contrast, even higher geothermal gradients do occur in the modern Raton Basin due to west-east groundwater flow (40–90°C/km) [55]. However, this groundwater flow is directed southeastward away from the Spanish Peaks and the highest geothermal gradients are concentrated downflow ~30 km south of our samples, such that the extremely high geotherms of this groundwater system are unlikely to apply to our study area. Moreover, as shown in online Supplementary Material Figure S2, geothermal gradients >50°C/km do not produce “good” fits in models, thereby suggesting that such high gradients are unable to produce the AHe data found at the Spanish Peaks.

Climate change can trigger exhumation events even in the absence of contemporaneous rock uplift by increasing the efficiency of fluvial and hillslope processes [7]. The Middle Miocene Climate Transition (MMCT) is an important global cooling event that occurred ~15–13 Ma [61], with the potential to enhance exhumation due to the stormier climate produced. However, the onset of the MMCT at ~15–13 Ma postdates the ~29–19 Ma onset of rapid exhumation across southeastern Colorado, and thus cannot be its primary cause. Furthermore, one would expect climate change-induced exhumation to be regional in scope, which is not true across the larger High Plains. Whereas we argue that kilometer-scale Miocene exhumation removed late Eocene-early Oligocene White River Group strata from the High Plains of southeastern Colorado, abundant outcrops of White River Group and its age equivalents occur in northeastern Colorado (Figure 1) and in adjacent portions of Wyoming, Nebraska, and South Dakota [28]. The southward increase in exhumation magnitude on the High Plains is inconsistent with climate change as the cause.

Drainage reorganization is another hypothetical cause of exhumation. Our study area lies in the modern drainage basin of the Arkansas River; if a stream capture event in Oligo-Miocene time expanded the Arkansas drainage basin, the resulting increase in river discharge might explain the observed exhumation event. Several authors have suggested that the current Arkansas River captured its upper reaches, which flow south along the RGR (Figure 1) sometime during the Neogene [62, 63]. However, very little is currently known about whether or when a stream capture event took place. If the Arkansas drainage basin was affected by a Neogene capture event, its erosional signal should propagate upstream (i.e., from east to west) through time; the fact that the Apishapa Dikes AHe dates are younger than those farther west at the Spanish Peaks suggests that drainage reorganization is an unlikely explanation for our High Plains erosion signal.

The evidence that geotherm relaxation, climate change, and drainage reorganization cannot, by themselves, explain our data leads us to conclude that southeastern Colorado’s High Plains must have experienced rock uplift during the Oligo-Miocene (Table 1). If the amount of rock uplift during that episode exceeded the magnitude of erosion, the High Plains experienced surface uplift during that event [64]. Surface uplift can be associated with an orogeny or with broad epeirogenic doming. Most of the faults and monoclines mapped in SE Colorado are likely Laramide features that were inactive by the time the Ogallala Formation was deposited [65]. Given that the SE Colorado rock uplift we have documented began between ~29 and 19 Ma, long after the end of the Laramide Orogeny (ca. 70–45 Ma), it is likely epeirogenic in nature.

Previous studies have also concluded that Colorado’s High Plains experienced Cenozoic epeirogenic surface uplift [66]. Marder et al. [66] presented evidence for regional west-directed back-tilting of SE Colorado’s High Plains that they attributed to an eastward migrating wave of dynamic surface uplift that developed ~4 Ma. Ostenaa et al. [67] documented four surface rupture events since 19 ka on the Cheraw Fault (Figure 1), which lies within the footprint of that epeirogenic uplift. However, they noted that the fault has experienced no more than 24–30 m of total vertical displacement during a combination of early Cenozoic (Laramide) compression and late Cenozoic normal faulting. They conclude that recent activity is a response to erosional unloading by the nearby Arkansas River and its tributaries, not evidence of Cenozoic tectonism. This postulated episode of dynamic topography [66] started much too recently to explain epeirogenic rock uplift starting between ~29 and 19 Ma, but as noted above, our Spanish Peaks and Apishapa Dikes thermal models permit significant post-7 Ma exhumation, so it is possible that our data record superposed exhumation due to both Oligo-Miocene and post-4 Ma epeirogenic episodes.

Bogolub and Jones [2] used flexural analysis to quantify the magnitude of “cryptic topography” (i.e., topography whose cause is not clear from surface evidence and that must be supported by subsurface mechanisms) on Colorado’s High Plains between 39 and 40° latitude, in the Limon area (Figure 1) just north of our transect at ~37.5° [2]. They concluded that 750–1250 m of post-100 Ma surface uplift occurred at the western edge of Colorado’s High Plains and that the magnitude of uplift increases to the south and west at least as far as 38° latitude. They compared four possible causes of epeirogenic surface uplift with their contour map of cryptic topography: (1) removal of the Farallon slab; (2) tilting associated with the formation of the RGR; (3) doming in response to ignimbrite flare-up magmatism; (4) dynamic topography. They concluded that their contours of cryptic topography best fit the magmatic doming hypothesis but noted that the great eastward extent of the surface uplift probably requires at least a component of slab-derived or dynamic topographic surface uplift (hypotheses 1 and 4). They also noted that placing a tighter constraint on the timing of this surface uplift event would greatly facilitate discrimination of the mechanism involved. We submit that the rock uplift event we have documented that initiated sometime between ~29 and 19 Ma is likely responsible for most or all of that surface uplift. That time of exhumation initiation is contemporaneous with the initiation of Rio Grande rifting [10, 11] and with the proposed ca. 25 Ma detachment of the Farallon slab [11], so both of those mechanisms remain plausible candidates. Ignimbrite flare-up volcanism was ongoing at that time but considering that the volcanic episode had already started by 37 Ma [29], it seems a less suitable candidate.

Another important constraint on proposed models to explain High Plains surface uplift is the pattern of exhumation. As noted above, NE Colorado’s High Plains experienced considerably less Cenozoic exhumation than the SE High Plains, indicating diminished surface uplift to the north. Our Two Buttes data document less exhumation eastward, suggesting that surface uplift also diminished toward the east. An AHe date of 24.7 ± 4.0 Ma from ~3900 m elevation on Trinchera Peak in the Sangre de Cristo Range (Figure 1) [11] is older than our 17.3 ± 1.3 Ma AHe date at approximately the same elevation but ~14 km to the east, suggesting diminished erosion to the west. Removal of the Farallon slab [11] does not explain well the northward-decreasing exhumation magnitude, unless only the SE Colorado to New Mexico portion of the slab detached at that time. Ignimbrite flare-up magmatism was voluminous in southwestern Colorado and western New Mexico but not in northern Colorado, so the magmatically induced doming hypothesis [17] easily explains northward-diminishing exhumation. However, because the main locus of magmatism occurred far to the west of our study area, it has difficulty explaining the high-magnitude exhumation in far southeastern Colorado at Two Buttes, and it also cannot explain a westward decrease in Oligo-Miocene exhumation between Spanish Peaks and Trinchera Peak (although this difference could be explained if Trinchera Peak lies on the down-dropped hanging wall of a Rio-Grande Rift-associated normal fault).

In light of the difficulties existing hypotheses have in explaining the pattern of High Plains exhumation, we suggest an alternative candidate surface uplift mechanism—a mantle drip [9, 36]. A mantle drip was proposed to be responsible for alkaline magmatic activity in the Spanish Peaks area between ca. 25 and 21 Ma [36], which is consistent with the finding that such drips can trigger melting of foundering lithospheric rocks [68]. Our data reveal that exhumation began simultaneously across all of southeastern Colorado, even as far as 200 km east of Spanish Peaks, at that same time. Foundering of a lithospheric drip beneath southeastern Colorado would produce a dome of isostatically compensated surface uplift there. Uplift of this dome would likely trigger rapid erosion of the soft sedimentary rocks that comprise the High Plains. The maximum magnitude and duration of erosion should decrease away from the center of the dome, so the mantle drip hypothesis can comfortably explain both the increasing AHe dates and the diminishing exhumation magnitudes observed to the north, east, and west of Spanish Peaks.

Low-temperature thermochronologic data from this study indicate that an episode of rapid exhumation, with a minimum exhumation rate of 180 m/m.y., occurred on the High Plains of southeastern Colorado during the Oligo-Miocene, between ~29 and 16 Ma. AHe data at Two Buttes, located in far southeastern Colorado, record the onset of this event sometime between 29 and 24 Ma. At the Spanish Peaks, ~200 km west of Two Buttes, numerous alkalic plutons were being emplaced at the same time, between ~25 and 22 Ma. Multi-sample thermal models using AHe dates from five samples of the Spanish Peaks region show that over 1.7 km of exhumation had occurred there prior to ~16 Ma, with the exhumation rate slowing to ~100 m/m.y. thereafter.

Exhumation ended at Two Buttes by ~12–9 Ma at the latest, based on the youngest likely age of the Ogallala Formation that unconformably overlies the Two Buttes intrusion. However, exhumation continued in the Apishapa Dikes region, where an additional >1 km of strata was removed by ~9–7 Ma. Considering that the Apishapa Dikes region was being exhumed at the same time that the Ogallala Formation was being deposited farther east, it is unlikely that any substantial, durable thickness thickness of Ogallala Formation was ever deposited as far west as the Apishapa Dikes. The development of a paleotopographic high at Spanish Peaks due to the local plutonism can explain the absence of Ogallala deposition there, in contrast to Ogallala sedimentation occurring right up to the Rocky Mountains farther north on the High Plains in Wyoming [59] and as proposed elsewhere along the range front.

We favor a regional epeirogenic surface uplift event that began in the latest Oligocene or earliest Miocene as the cause for both the exhumation pattern we observe and the anomalous height of southeastern Colorado’s High Plains. Any model that proposes to explain the surface uplift of Colorado’s High Plains must explain the fact that over 3 km of material has been removed from southeastern Colorado since the latest Oligocene but considerably less erosion has occurred in the northeastern part of the state during that same time interval. Thermal effects associated with the ignimbrite flare-up and formation of the RGR can explain this dichotomy, but detachment of a mantle drip centered beneath southeastern Colorado better fits the pattern and spatial scale of the exhumation event.

All data to support this study is published in this paper or in the supplementary materials.

The authors declare that there is no conflict of interest regarding the publication of this paper.

We thank Rachel Landman for acquiring apatite (U-Th)/He data for the original four samples in this study as part of her CU-Boulder PhD thesis and for the idea in her thesis of using (U-Th)/He thermochronology of Cenozoic intrusions on the High Plains to constrain the High Plains exhumation history. Shelby Litton and Aaron Bell did emplacement depth work that kindled our interest in doing low-temperature thermochronology on Huerfano Butte.

We thank Rich Ketcham for his advice and assistance using the new multi-sample modeling capabilities of HeFTy. Matt Morgan, the Colorado State Geologist, and an anonymous reviewer provided inciteful comments that greatly improved this manuscript.

We thank the Undergraduate Research Opportunities Program at the University of Colorado for an award to SJK that supported this work. National Science Foundation awards EAR-1559306 and -1126991 to RMF funded the instrumentation used to acquire the (U-Th)/He data in this study. Funding for Skye Fernandez’s work on this project was provided by the UNAVCO Research Experiences for Solid Earth Science for Students (RESESS) program.

Publication of this article was funded by the University of Colorado Boulder Libraries Open Access Fund.

The supporting information presents extended single-grain apatite and zircon (U-Th)/He data, and thermal modeling settings and constraints used for all HeFTy models. Thermal models discussed as being less geologically viable are also included (Spanish Peaks transect at higher geothermal gradients (50 and 60°C/km and 2–3 km emplacement depth) and Apishapa Dikes samples exploring reburial by Ogallala deposition).

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

Supplementary data