Rise of the Colorado Plateau: A Synthesis of Paleoelevation Constraints From the Region and a Path Forward Using Temperature-Based Elevation Proxies

The Colorado Plateau’s complex landscape has motivated over a century of debate, key to which is understanding the timing and processes of surface uplift of the greater Colorado Plateau region, and its interactions with erosion, drainage reorganization, and landscape evolution. Here, we evaluate what is known about the surface uplift history from prior paleoelevation estimates from the region by synthesizing and evaluating estimates 1) in context inferred from geologic, geomorphic, and thermochronologic constraints, and 2) in light of recent isotopic and paleobotanical proxy method advancements. Altogether, existing data and estimates suggest that half-modern surface elevations were attained by the end of the Laramide orogeny (∼40 Ma), and near-modern surface elevations by the mid-Miocene (∼16 Ma). However, our analysis of paleoelevation proxy methods highlights the need to improve proxy estimates from carbonate and floral archives including the ∼6–16 Ma Bidahochi and ∼34 Ma Florissant Formations and explore understudied (with respect to paleoelevation) Laramide basin deposits to fill knowledge gaps. We argue that there are opportunities to leverage recent advancements in temperature-based paleoaltimetry to refine the surface uplift history; for instance, via systematic comparison of clumped isotope and paleobotanical thermometry methods applied to lacustrine carbonates that span the region in both space and time, and by use of paleoclimate model mediated lapse rates in paleoelevation reconstruction.


INTRODUCTION
Scientists have long debated the uplift of the Colorado Plateau region (CP), including the modern Colorado River (CR) drainage from the Rocky Mountains to where it exits the modern plateau southwest of Grand Canyon (GC; Figure 1). This region was a tectonically stable sedimentary basin near sea-level until the Laramide orogeny initiated in the Late Cretaceous (Hunt, 1956;Nations et al., 1991) and now resides ∼2-3 km above sea level. Proposed mechanisms for plateau rise make testable predictions of uplift patterns due to processes including: lower-crustal flow (McQuarrie and Chase, 2000) or shear-removal of lithospheric mantle by the Farallon slab during the Laramide (Hernández-Uribe and Palin, 2019); lithospheric foundering (Bird, 1979) or hydration (Porter et al., 2017) in the mid-Cenozoic; and mantle upwelling (e.g., Moucha et al., 2009) or isostatic response to denudation (e.g., Lazear et al., 2013) in the late Cenozoic. Accurate paleoelevation constraints can evaluate hypothesized uplift mechanisms and provide context for CR evolution.
Here, we build on previous syntheses (e.g., Cather et al., 2012) to examine proposed CP paleoelevation estimates and highlight opportunities to refine the elevation history in light of recent proxy method advances. Decades of paleoelevation reconstructions combined with constraints from geology, geomorphology, and thermochronology support the interpretation that the CP was elevated through multi-stage, spatially differential uplift since the Laramide. Paleobotanical and geochemical proxy methods and the use of climate models to infer paleoelevation from proxy data have advanced significantly in the past 10-15 years and can be leveraged to refine and expand temperature-based paleoelevation reconstructions across the CP region.
Overview CP pre-Cenozoic stratigraphy is characterized by relatively flat-lying sedimentary sequences, which suggests that the region was a comparatively stable sedimentary basin near or below sea-level for most of the Paleozoic-Mesozoic (Hunt, 1956). Coincident with onset of the Laramide (≥90 Ma; Carrapa et al., 2019), CP uplift likely swept from southwest to northeast ( Figure 1B) suggested by the diachronous disappearance of Cretaceous Interior Seaway deposits (Nations et al., 1991;Raynolds and Johnson, 2003); barbed NE-flowing tributaries along sections of GC; and substantial (∼2 km) Paleozoic-Mesozoic to late Cretaceous-Paleogene denudation on the southwestern plateau indicated by erosional unconformities (e.g., Young, 2001;Hill et al., 2016). Imbricated gravels deposited in incised paleo-canyons (e.g., Milkweed-Hindu, Peach Springs, Salt River) on the southwestern CP margin in the late Cretaceous and Paleocene-Eocene support this view (e.g., Young, 2001;Potochnik et al., 2011), although the magnitude of Paleogene relief is debated (e.g., Karlstrom et al., 2014;Young and Crow, 2014).
The mid-Cenozoic ( Figure 1C) was a period of internal drainage from the Sevier highlands to the west (e.g., Davis et al., 2009;Ibarra et al., 2021) and uplifts to the northeast (e.g., Cather et al., 2012). Drainage reorganization and tilting are evidenced by widespread fluvial and lacustrine deposits across the CP (Hunt, 1956) and lacustrine carbonate oxygen isotopic records (Davis et al., 2009). Late Cenozoic ( Figure 1D) continental extension caused subsidence on the western and southern CP margins (McQuarrie and Wernicke, 2005). Drainage reversed to the southwest (e.g., Potochnik et.al., 2001), and a scarcity of late Cenozoic deposits suggests regional erosion (Hunt, 1956).
The modern CR was integrated from its headwaters in the Rocky Mountains through GC coincident with development of the Salton Trough by ∼6 Ma, recorded by CR deposits in the Muddy Creek Formation (e.g., Longwell, 1946;Lucchitta, 1972;Karlstrom et al., 2014). The CR reached its current base-level in the Gulf of California at 4.80-4.63 Ma (Crow et al., 2021).

Additional Uplift Constraints
Studies have investigated CP paleoelevation through reconstructions of erosion and river incision, which may occur in response to relief, drainage, and base-level changes that accompany surface uplift. Erosional exhumation may drive cooling of crustal rocks, which can be recorded by mineral thermal histories reconstructed from thermochronology (Reiners and Shuster, 2009), or constrained by other burial and thermal history information (Nuccio and Condon, 1996). For instance, thermochronometer cooling ages and peak burial temperature constraints from the southwestern CP surface suggests >1.5 km of denudation during the Laramide (e.g., Flowers et al., 2008), and ∼1.9-2.5 km of denudation prior to the Neogene (Ryb et al., 2021), supporting the interpretation that uplift of the southwestern CP margin initiated during the Laramide and may have continued into the late Cenozoic. Recent (<10 Ma) cooling of the southwestern CP surface rocks may reflect an erosional response to integration of the CR and a drop in base-level (e.g., Lee et al., 2013;Murray et al., 2016).
Cooling histories from canyon-rim and river-level samples track paleo-relief of GC (Flowers et al., 2008;Flowers and Farley, 2012;Flowers and Farley, 2013;Lee et al., 2013;Karlstrom et al., 2014;Winn et al., 2017) and constrain minimum plateau paleoelevation. A synthesis of thermochronologic datasets by Karlstrom et al. (2020) indicates segments of the GC were incised at different times before becoming integrated 6-5 Ma; the Western GC was likely a paleo-canyon with substantial relief incised by a northeast-flowing river by ∼65 Ma, whereas the Eastern GC was likely incised by a NW-flowing paleo-river by ∼20 Ma. These trends are compatible with multi-stage uplift prior to ∼20 Ma. Additional erosion constraints from volcanic deposits (e.g., Aslan et al., 2008;Donahue et al., 2013), cave deposits Polyak et al., 2008), and strath terraces show increased incision rates since ∼6 Ma, for example in eastern GC  and the Lower Colorado River , potentially reflecting small-scale (<400 m) uplift of the southwestern plateau related to Neogene volcanism (Crow et al., 2011;Crow et al., 2019).
River longitudinal profiles may also constrain CP uplift timing (Roberts et al., 2012). However, complications due to varying rock erodibility (Cook et al., 2009;Pederson and Tressler, 2012), normal faulting and paleocanyon integration , as well as drainage reorganization and resulting changes to stream power (Schwanghart and Scherler, 2020) limit straightforward correlation between river profiles and uplift. Nevertheless, an inverse correlation between channel steepness and upper mantle p-wave velocity in the Virgin River drainage suggests epeirogenic uplift affects modern CP river profiles (Walk et al., 2019). Projection of relict channels buried by Miocene and younger basalt (Hamblin et al., 1981), and knickpoint-bound upper reaches in CR tributaries (Darling and Whipple, 2015), indicate ∼1 km relief growth both prior to and since the Miocene. Such constraints do not distinguish between relief growth in response to CP uplift or subsidence at the plateau margin (Ott et al., 2018).

PREVIOUS PALEOALTIMETRY
Many workers have developed and applied direct paleoelevation proxies to resolve CP paleoelevation, though little agreement on methods or results has been reached ( Figure 2; Supplementary  Table S1). For instance, CP paleoelevation estimates based on basalt-vesicle paleobarometry (Sahagian et al., 2002a;Sahagian et al., 2002b;Sahagian et al., 2003), while elegant in theory, are problematic (Bondre, 2003;Libarkin and Chase, 2003), and the approach has not been replicated. Below, we evaluate CP paleoelevation estimates from oxygen-isotopic and paleothermometry-based methodologies.

Paleobotanical Thermometry-Based Paleoaltimetry
The steepest temperature gradients on Earth are vertical, making surface temperature proxies attractive paleoaltimeters. Proxies have been used to estimate depositional mean annual surface temperature (paleo-MAST), which may be combined with coeval near-sea-level temperature reference data and an estimated rate of temperature decrease with elevation ("lapse-rate") to infer paleoelevation.
Studies of key Eocene-Oligocene sites (Figure 1) highlight this variability and permit several CP uplift histories (Figure 2; Supplementary Table S1 and references therein). MAST estimates based on different methodologies applied to the same fossil assemblage span a ∼13°C range for the Florissant flora, and other key floral sites at Creede and Pitch-Pinnacle have similarly wide ranges of temperature estimates (Leopold and Zaborac-Reed, 2019). Existing paleoelevation reconstructions based on paleobotanical-proxy temperatures and estimated lapse-rates (either study-specific lapse-rates, or using Wolfe, 1992a;Meyer, 1986;Meyer, 1992; discussed in Challenges and Recent Advances section) contain significant uncertainty-in some cases estimates range from <0.5 to >4 km for the same flora (e.g., Florissant; Figure 2C). Additional estimates calculated here using the refined Meyer (2007) approach (i.e., sea-level correction, local lapse-rate), based on paleobotanical-proxy temperature estimates that lack published paleoelevation evaluations, show similarly broad possibilities (Supplementary Table S1). This highlights the need for an evaluation of paleobotanical methods and their application, as well as the independent verification of temperature estimates and paleoelevation approaches.

Clumped Isotope Thermometry-Based Paleoaltimetry
Carbonate clumped isotope (Δ 47 ) thermometry (Ghosh et al., 2006;Eiler, 2007;Eiler, 2011) is commonly used for reconstructing surface temperature and paleoelevation (Huntington and Lechler, 2015). The Δ 47 value refers to the ratio of "clumped" carbonate molecules containing heavy isotopes 13 C and 18 O, relative to the random distribution of isotopes among isotopologues. Heavy isotope clumping is directly temperature-dependent, and the carbonate temperature estimates (TΔ 47 ) are thermodynamically based; thus, the method can be applied in the geologic past without concern for plant evolutionary changes or variable paleobotanical-proxy methods, provided authigenic carbonates are available and well preserved. Clumped isotope data also constrain water δ 18 O values, and were used by Licht et al. (2017) to assess diagenesis and complement δ 18 O-based CP paleoaltimetry. TΔ 47 for the Miocene lacustrine Bidahochi Formation (16-6 Ma), combined with TΔ 47 of coeval near-sealevel deposits and a modern-lake-carbonate TΔ 47 lapse-rate, estimated a paleoelevation of ∼1.9 km, suggesting that nearmodern surface elevations of the southern CP were attained by 16 Ma (Figure 2; Huntington et al., 2010), consistent with thermochronology data .

Existing Estimates and Remaining Questions
Existing paleoelevation estimates are compatible with multistage, spatially differential uplift (Figure 2). Laramide-age deposits imply ≳1 km rise of the northwestern and southeastern CP by ∼50 and ∼38 Ma, respectively, consistent with paleo-canyon and paleo-relief studies. Geologic constraints from mid-Cenozoic (∼38 Ma) deposits suggest half-to nearmodern elevations in the northeastern CP (Steven and Ratté, 1965;Cather et al., 2012), consistent with some paleobotanical-Frontiers in Earth Science | www.frontiersin.org June 2021 | Volume 9 | Article 648605 based estimates. Late-Cenozoic (∼16-6 Ma) deposits suggest near-modern elevations in the southwestern CP. Altogether, this suggests half-modern surface elevations by ∼38 Ma, and near-modern surface elevations by ∼16-6 Ma. The timing supports hypothesized Laramide and mid-Cenozoic uplift mechanisms but does not preclude more recent modest surface uplift or subsidence, or rock uplift that is balanced by erosion (<6 Ma on the southern CP). Except for some paleobotanical-based estimates in the southern Rocky Mountains, for all regions and time periods paleoaltimetry data suggest that elevations were likely lower than or close to modern; this is consistent with geodynamic models and other constraints (e.g., Cather et al., 2012) and implies multi-stage uplift without a complex late-Cenozoic uplift and subsidence (c.f., Wolfe, 1994). The following remain to be resolved: the 1) amount and spatio-temporal pattern of Laramide-related surface uplift; 2) rate and timing of southern Rocky Mountains uplift, with implications for drainage reversal; 3) magnitude of possible ≲16 Ma CP elevation change; and 4) cause of variation in paleobotanical paleoelevation estimates.

Challenges and Recent Advances
Temperature-based paleoaltimetry methodologies are challenging, but have advanced substantially in the last decade and have potential to help address these remaining questions. Multiple factors independent of elevation influence MAST, including short-(e.g., El-Nino Southern Oscillation) and longterm climate changes (e.g., Cenozoic cooling). Additionally, climate variability on decadal to millennial time scales may be integrated differently depending on the resolution of the proxy, and seasonal bias may impact different proxy types and settings (e.g., Burgener et al., 2016). MAST varies spatially due to local vegetation cover, aspect, precipitation regime, or proximity to water bodies/the ocean (e.g., Kelson et al., 2020), which can change with secular or periodic changes in climate-potentially complicating comparison of known nearsea-level deposits with contemporaneous inland deposits of unknown elevation. Workers have addressed such complications in various ways; e.g., Huntington et al. (2010) analyzed the Bouse Formation as a near-sea-level anchor in two sub-basins, and opted to use estimates from the warmer sub-basin farther from the coast instead of the marine-biased colder subbasin proximal to the coast. Feng and Poulsen (2016) used paleoclimate simulations of Eocene North America to make improved predictions of contemporaneous sea-level temperatures for Cordilleran floral sites.
There have also been varied approaches to determining lapserates, such as applying the average modern global (5.5°C/km), regional (3°C/km; Wolfe, 1992b), or local (e.g., 5.9°C/km;Meyer, 1986) lapse rate for the CP (e.g., Meyer, 2007). Lapse-rate may change with aridity/humidity and latitude (e.g., Neumann, 1955;Li et al., 2015), leading to skepticism about whether modern or global lapse-rates are applicable to ancient reconstructions. Paleoclimate models have advanced significantly since initial rates were proposed and demonstrate that significant error can be generated by applying linear modern lapse-rates to ancient climates and paleogeographies; model-mediated lapse-rates should be a key component of future temperature-based paleoaltimetry studies (e.g., Feng and Poulsen, 2016;Farnsworth et al., 2021;Botsyun and Ehlers, 2021).
Confidence and accuracy in paleobotanical temperature estimates also have improved significantly. Advances include new standardizations in leaf physiognomic methods (e.g., Peppe et al., 2011), regional corrections and character reevaluations in multivariate methods (e.g., Spicer et al., 2021), larger comparative datasets and refined statistical methods in bioclimatic analyses (e.g., Chevalier, 2019), and improvements in paleofloral identifications and taxonomy (e.g., Manchester, 2014). Additionally, the number of available CP-region collections has increased; local and regional stratigraphy refinements have improved context and relationships between existing floras (e.g., Prothero, 2008), allowing for subsampling and clear division in collection units; and geochronology has refined ages and durations of distinct floras (e.g., Lipman and Bachmann, 2015), allowing for higher-resolution reconstructions and more precise application of corrections/ lapse-rates.
Challenges of TΔ 47 paleoaltimetry include: the need for up to tens of milligrams of carbonate to achieve precise temperature estimates, and relatively large elevation uncertainties (∼0.5 km; Huntington and Lechler, 2015); time integration during carbonate precipitation that can vary depending on climate (e.g., Kelson et al., 2020); and potential for carbonate temperature resetting by post-depositional diagenesis and heating (e.g., Huntington et al., 2011). Despite this, confidence and accuracy in Δ 47 measurements have improved significantly in the last decade-with capability for higher sample throughput, more robust calibrations, new carbonate standards, and improved instrument precision and data normalization (e.g., Bernasconi et al., 2018;Petersen et al., 2019;Saenger et al., 2021). Evaluation of carbonate diagenesis also has improved (e.g., Lacroix and Niemi, 2019), and recent work indicates that the burial/ exhumation history can be reconstructed from reset TΔ 47 values with implications for paleotopography (e.g., Ning et al., 2019).
While transfer functions have been proposed to translate carbonate TΔ 47 to MAST by assuming the seasonal timing of carbonate formation (Hren and Sheldon, 2012), using modern-carbonate-based lapse-rates may capture similar seasonal bias of ancient samples and reduce assumptions (e.g., Li et al., 2021). Quaternary lacustrine carbonates from southwest North America apparently record summer-biased temperature instead of MAST (Huntington et al., 2010). However, this bias may not hold for the paleo-record, particularly under Cretaceous and early Paleogene greenhouse climates. Novel isotopic approaches (tripleoxygen isotopes, 13 C-excess) can help correct for evaporative effects and estimate the unevaporated δ 18 O of catchment waters, which may aid in constraining paleohydrology and provide context for estimating carbonate temperature seasonality and lapse-rates.

Future Opportunities
These advances can be leveraged to refine temperature-based proxy reconstructions and expand their use in the CP region. The region is well suited to temperature-based paleoaltimetry because its paleo-latitude has changed little for ∼80 Ma; the mid-latitude setting allows for a steeper and better-constrained temperature-elevation gradient, thus minimizing elevation uncertainty; and target deposits are widespread in space and time and unlikely to have been deeply buried and diagenetically altered.
Essential to resolving paleobotanical estimates is independent comparison of carbonate TΔ 47 and wellstudied floral records from the same lacustrine formations (e.g., Florissant, Creede, Antero). Systematic comparison could expand the applicability and comparability of techniques across additional Laramide-age localities ( Figure 1A) that may contain either carbonates or fossil floras (Claron, Flagstaff, Uinta, Piceance, and San Juan Basins), and potentially provide context for interpreting 14-19°C early calcite cements from the ∼64 Ma southwestern CP Music Mountain Formation (Huntington et al., 2011). Additionally, there is opportunity to re-evaluate late-Cenozoic Bidahochi and Bouse Formation paleotemperatures (Huntington et al., 2010) with higher temporal/spatial resolution. In all cases, we recommend using paleoclimate model-mediated lapse-rates and sealevel corrections to minimize uncertainty and align methods for extrapolating temperature to elevation. Our analysis suggests that temperature-based paleoaltimetry methods are key to refining CP paleoelevation, with implications for understanding uplift mechanisms and CR drainage evolution.

AUTHOR CONTRIBUTIONS
EOH and KH (corresponding authors) conceived of the manuscript, with input from EGH. EOH did much of the literature review, first drafts of writing, and compilation of datasets. KH contributed to the writing and literature review on the Colorado Plateau, paleoaltimetry, and clumped isotopes. EGH contributed to the writing and literature review for paleobotanical methods. PSG drafted Figure 1 and contributed to the writing and literature review of river profile modelling. CB drafted Figure 2. All authors provided input on data analysis, interpretation, and writing.