Over the past 10 Ma, the high-relief landscapes of the Colorado Plateau–southern Rocky Mountains region have been shaped by erosional processes. Incision rates have increased in the southern Rocky Mountains, the Colorado River system has been superimposed across buried Laramide structures as it was integrated from the Rocky Mountains to the Gulf of California, the modern Grand Canyon formed, and there has been widespread denudation of the Canyonlands region of the Colorado Plateau. We examine the spatial and temporal distribution of erosion and its associated isostatic rebound since 10 Ma. Erosion estimates come from apatite fission track (AFT) and apatite (U-Th)/He (AHe) thermochronometric studies at 14 sites across the region, including recent AHe data with ages younger than 12 Ma, and from ca. 10 Ma 40Ar/39Ar dated basalt paleosurfaces at 55 locations on the perimeter of the Colorado Plateau and in the southern Rocky Mountains. Estimated eroded thickness is added to modern topography above numerous control points to reconstruct a 10 Ma paleosurface across the region (referenced to modern elevations); this also yields an eroded thickness volume. Erosion has been spatially variable since 10 Ma: we find widespread denudation with as much as 2 km of incision along rivers in the Canyonlands region of Utah, 1–1.5 km of incision along rivers exiting the Rocky Mountains onto the eastern piedmont since 6 Ma, ∼1 km removed across the high peaks of the southern Rocky Mountains since 10 Ma, and little net erosion in the Basin and Range.

Post–10 Ma flexural isostatic response to the eroded volume is calculated using known variable elastic thickness. This rebound caused much of the Colorado Plateau region to undergo more than 800 m of rock uplift, exceeding 1 km in local areas in the Canyonlands and southwestern Colorado. The Lees Ferry and Glen Canyon areas have been isostatically uplifted >500 m relative to the eastern Grand Canyon and the Tavaputs Plateau has been isostatically uplifted 400 m relative to Browns Park. This differential rock uplift driven by erosional isostasy has created or accentuated many of the features of the modern landscape. This component of rock uplift is “removed” by adding the eroded thickness onto modern topography, then subtracting the calculated rebound. The resulting (pre-erosion and pre-rebound) map provides a model of the 10 Ma landscape, neglecting any tectonic uplift contribution to regional elevations. This model suggests the presence of internal drainages on the Colorado Plateau, that the elevation of the Green River Basin and the Tavaputs Plateau were subequivalent, allowing the Green River to flow southward, and shows high topography in the Rocky Mountains that mimicked modern topography, but with potentially lower relief. Future refinements of both the timing and magnitude of differential erosion and rebound models provide an avenue for improved models for Cenozoic landscape evolution of the region.

This paper is an advance over previous studies that focused just on the Colorado Plateau. Here we evaluate isostatic response to erosion in an extended region that includes parts of the Basin and Range, Colorado Plateau, southern Rocky Mountains, and eastern piedmont of the Rocky Mountains. We find that erosion of the southern Rocky Mountains and eastern piedmont is comparable to that of the Colorado Plateau and that the flexural isostatic rebounds of all these regions are coupled and cannot be considered in isolation. Furthermore, we focus on the 10 Ma time frame, rather than the 30 or 70 Ma period of previous researchers, as the key time frame during which the modern landscape rapidly evolved. In addition, the use of AFT and AHe thermochronometric constraints on thicknesses and ages of now-eroded sediments has solved key problems that hampered previous erosion studies. Data and analyses of regional post–10 Ma differential erosion and its resulting differential isostatic rebound provide essential constraints for any viable models for landscape evolution in this classic region.


Uplift of orogenic plateaus and erosional topography are globally important topics related to how landscapes respond to mixed climatic, geomorphic, and tectonic influences. Debates about dynamic topography of broad plateau regions (Braun, 2010; Karlstrom et al., 2012), global climate forcing (Molnar, 2004), the influence of climate on uplift of mountain belts (Whipple, 2009), and resolving surface expressions of tectonism (Kirby and Whipple, 2012) are important and relevant to the study area discussed in this paper, and require careful geomorphic definitions and analysis of different types of uplift and erosion. The terminology used herein follows England and Molnar (1990) and is provided in Table 1. This paper highlights an extended region of the southwestern United States, including the Colorado Plateau and southern Rocky Mountains, that has undergone drainage integration of the Colorado River system and dramatic denudation and consequent erosion-driven isostatic rock uplift within the past 10 Ma. The timing and spatial distribution of these events are critical to understanding the forcing mechanisms behind this evolving landscape, and are the subject of this paper. We focus on the component of rock uplift driven by isostatic responses to erosional exhumation, a sometimes-overlooked component of the overall rock uplift history in plateau landscapes. All other rock uplift mechanisms are lumped into the category of tectonic uplift for the purpose of this paper. The terms isostatic rock uplift, isostatic rebound, isostatic deflection, and rebound are used interchangeably for flexural rock uplift and/or subsidence in response to surface unloading and/or loading. The term flexural refers to support of lithospheric loads by combined elastic restoring moments and buoyancy forces in the crust and mantle, and the contribution of each component is dependent upon the lateral scale of the surface load.

The Rocky Mountain–Colorado Plateau region (Fig. 1 and Supplemental File 11) is a spectacular erosional landscape with 3–4 km elevations of the southern Rocky Mountains in Colorado adjacent to the deeply eroded sandstone canyon lands in Utah, and the Grand Canyon in Arizona. This region has been uplifted since ca. 75 Ma, when it was at sea level and was part of the Cretaceous interior seaway. At that time, the stratigraphic section of the region included ∼3–4 km of Paleozoic and Mesozoic strata topped by thick Cretaceous marine shales and sandstones. There has been a vigorous century-long debate about how the Late Cretaceous low-elevation, low-relief region was transformed into the incised high-elevation plateau we see today. One of the major limitations common to all attempts to reconstruct the landscape evolution of this region is that erosion has removed much of the geological record of the intermediate, mid-Tertiary, stages in the evolution of the landscape. In particular, there are few preserved remnants of Tertiary aggradation surfaces across the Canyonlands region of Utah and in the southern Rocky Mountains that record the progressive erosional history.

The problem of the lack of preserved surfaces has been greatly aided by newly emerging data sets that provide better constraints on landscape evolution in this region. Thermochronologic data offer increasingly resolved images of now-eroded landscapes. Apatite fission track (AFT) and apatite (U-Th)/He (AHe) thermochronology, especially when modeled together (Lee et al., 2013), can constrain the cooling history of rocks within the interval of ∼110–40 °C as they approached the Earth’s surface. These data shed light on past landscapes that were several kilometers higher than present exposures, assuming geothermal gradients of 25 °C/km and modern surface temperatures of ∼10 °C (e.g., Fitzgerald et al., 2009; Kelley and Chapin, 2004; Flowers et al., 2008, Hoffman et al., 2011; Landman and Flowers, 2012). Improved databases on rates of river incision are emerging from dating of the fragmentary record preserved in elevated river terraces using combined geochronologic techniques (Ar-Ar, U series, cosmogenic surface exposure and burial dating, and optically stimulated luminescence [e.g., Pederson et al., 2002; Crow et al., 2008; Aslan et al., 2011; Darling et al., 2012]). Recent studies have also provided new insights on basin remnants such as the 34–26 Ma Chuska erg (Cather et al., 2008), 15–6 Ma Bidahochi Formation (Dallegge et al., 2001; Dickinson, 2012), Muddy Creek Formation, and the Imperial Formation that provide data on the sediment budget for the Colorado River (Dorsey, 2010; Dorsey and Lazear, 2013). An integration of these data sets provides the opportunity in this paper to quantify post–10 Ma erosional history and its flexural isostatic response in the Colorado Plateau–southern Rocky Mountains region.

Differential uplift and exhumation of the Colorado Plateau and southern Rocky Mountains regions are entwined and have occurred in three main stages: Laramide deformation (75–45 Ma), mid-Cenozoic (35–15 Ma), and Neogene (10 Ma–present) (Karlstrom et al., 2012; Cather et al., 2012). The first stage (75–45 Ma) is characterized by uplift and erosion of the southern Rocky Mountains and aggradation on the adjacent eastern piedmont and Colorado Plateau during Laramide deformation. Laramide contractional tectonism produced basement-cored uplifts of the Rocky Mountains as well as upwarps on the Colorado Plateau, including the Kaibab, Uncompahgre, Monument, Defiance, and Zuni uplifts (see Fig. 1). Deep basins, including the San Juan, Uinta, Denver, and Raton Basins, received synorogenic clastic sediment shed from progressively growing regional uplifts (Cather, 2004). A second uplift and/or exhumation stage (25–15 Ma) is characterized by erosion of the Kaibab Plateau (Flowers et al., 2008) and southwestern portion of the Chuska erg (Cather et al., 2008), and continued erosion of the southern Rocky Mountains (Feldman, 2010; Garcia, 2011). Uplift and/or exhumation stage three (10 Ma–present) is defined by widespread denudation of the Canyonlands region of Utah, continued erosion and incision of modern canyons in the southern Rocky Mountains, integration of the modern Colorado River system across the Colorado Plateau to the Gulf of California, and incision of the modern Grand Canyon.

A goal of this paper is to reconstruct a 10 Ma paleosurface that can be used as a datum relative to modern elevations to estimate eroded thickness and flexural isostatic rebound since 10 Ma across the Colorado Plateau–Rocky Mountains region. Our analysis extends east-west from the Basin and Range, across the southern Rocky Mountains, to the western Great Plains (long ∼102°W–115°W), and north-south from southern Wyoming to Mexico (lat 44°N–32°N). This paper advances similar previous studies (Pederson et al., 2002; Roy et al., 2009) by (1) focusing on the post–10 Ma time frame that encompasses increased incision in the region; (2) utilizing numerous ca. 10 Ma 40Ar/39Ar dated basalt paleosurfaces on the perimeter of the Colorado Plateau and in the southern Rocky Mountains to define our surface; (3) using new AHe thermochronology data that provide timing and eroded thickness for mapping differential erosion in regions without preserved aggradation surfaces; and (4) evaluating the isostatic response of the coupled southern Rocky Mountain–Colorado Plateau system, recognizing that drainages have been integrated between these two physiographic provinces since the Oligocene (Hansen, 1965; Larson et al., 1975; Aslan et al., 2010) such that flexural isostatic rebound of the Colorado Plateau, southern Rocky Mountains, and eastern piedmont are coupled and cannot be treated in isolation. Eroded thickness is treated as a surface load on the lithosphere, and flexural isostatic deflection is computed using a variable effective elastic thickness for the lithosphere (Te) to accommodate a factor of 5 change in Te (10–50 km) across the study area (Lowry et al., 2000). We present the methods and results of our study first, and provide a detailed comparison of our results with previous research later in this paper.


A variety of data are utilized herein to define a topographic model for the 10 Ma paleosurface. Details regarding data constraints in specific regions are provided in the following section. Here we discuss the characteristics and uncertainties associated with specific data types. Our data sources include (1) the present elevation of ca. 8–12 Ma basalt flows that armored and preserved erosional remnants of the paleosurface; (2) thickness and timing of erosion provided by ca. 10 Ma AFT and AHe thermochronometry dates; (3) the history of the Chuska erg (Cather et al., 2008) and Bidahochi Formation (Dickinson, 2012) on the southern Colorado Plateau; (4) incision of the Miocene–Pliocene Ogallala Group on the eastern piedmont of the Rocky Mountains (Leonard, 2002; McMillan et al., 2002); and (5) geologic histories of the Rio Grande Rift (Chapin and Cather, 1994), the Basin and Range (Faulds et al., 2001), and the Gilbert erosion surface (Hansen, 1986).

The primary source for basalt paleosurface constraints is the North American Volcanic and Intrusive Rock Database (NAVDAT; http://navdat.org), from which dated basalt flows were extracted within geographic coordinates 116°W–102°W and 31°N–44°N using an age range from 8 to 12 Ma. Additional basalt data points for the southern Rocky Mountains were supplied in Aslan et al. (2010, 2011) and Cole (2011), giving a total of 55 basalt control points documented in Supplemental File 22 and shown in Figure 2 (white circles). Basalt flows occupy local low points in the terrain at the time of eruption, and tend to armor the ground surface against further erosion and river incision, resulting in inverted topography. When river gravels underlie flows, the basalts represent points in a paleoriver system, and provide the most accurate estimates of incision when the modern counterpart river is still in close proximity to the flow-mantled paleoriver. The difference between the armored elevation and the surrounding river levels and erosion surface elevations provides an estimate of net eroded thickness since the time of the basalt. If differential rock uplift between the armored surface and modern topography has occurred due to faulting, salt collapse, or tilting, these areas are useful for differential incision and fault-dampened incision studies (Kunk et al., 2002; Crow et al., 2008), but are avoided for present purposes. Therefore, basalt data from a known salt collapse near Carbondale, Colorado, were excluded from our database (Kunk et al., 2002). When there is no association of the basalt with river gravels or a nearby modern river, the elevation of the basalt relative to the local topographic base level at the time of emplacement is unknown, and it is difficult to assess uncertainty in eroded thickness. For example, the 10 Ma Red Butte flow is close to the rim of modern Grand Canyon, but could have stopped flowing here even if a paleocanyon existed. In general, basalt flows are found on the perimeter of the Colorado Plateau, in the Basin and Range province to the south and west of the plateau, and in the southern Rocky Mountains. For ages different than 10 Ma, elevations are adjusted to 10 Ma using approximate erosion rates or are considered equivalent to 10 Ma elevations if they are in a region where 12–8 Ma denudation was slow or minimal.

The AHe and AFT thermochronologic data provide estimates of eroded thickness and timing of erosion by measuring the elapsed time since apatite crystals passed through specific temperatures as they were exhumed toward the surface. Estimates of eroded thickness derived from these data have a considerable uncertainty associated with the assignment of closure temperature (40–80 °C for AHe and ∼110 for AFT). The AHe closure temperature varies with grain size, uranium and thorium content, and cooling rate, and is complicated by zoning of uranium and thorium within apatite crystals (Shuster et al., 2006; Flowers et al., 2009). The AHe closure temperature is most strongly affected by U and Th content; grains with high concentrations of U and Th accumulate larger amounts of radiation damage, which slows helium diffusion and increases closure temperature (Shuster et al., 2006; Flowers and Kelley, 2007). The AFT closure temperature is also dependent on the chlorine content of apatite and cooling rate. Factors affecting the kinetics of helium diffusion and track annealing (and related closure temperature) have been incorporated into the radiation damage accumulation and annealing model (RDAAM) of Flowers et al. (2009) and in the kinetic model HeFTy (Ketcham, 2005). HeFTy and RDAAM can be used to generate thermal histories and constrain temperatures from several samples on an age-elevation traverse or in a drill hole at specific times, in this particular case at 10 Ma.

Converting thermochronologically constrained paleotemperatures to an inferred depth requires knowledge of past geothermal gradients, thermal conductivity, and interpretation of the cooling history in terms of thermal events and changes in heat flow. At the level of certainty needed for this paper there are two cases. The first is for Colorado Plateau control points, where we use a geothermal gradient of 20–25 °C/km, surface temperature of 10–25 °C, and a closure temperature of ∼50–70 °C for the minimum AHe ages. Thus, AHe ages of ca. 10 Ma in the Colorado Plateau are inferred to reflect ∼2 km ± 1 km of burial. In the second case, for the Rocky Mountains, we use a geothermal gradient of 40 °C/km, surface temperature of 4–7 °C, and closure temperature of 50 °C for minimum AHe ages. For these 10 Ma AHe ages, we infer that rocks were buried 1.0 ± 0.5 km. Systematic age-elevations transects suggest relatively stable geotherms during erosional exhumation. Thus, interpretation of thermochronologic data across the region provides important, although imprecise, constraints on both spatial and temporal erosional history that are not available from other data sets such as incision history along a river course, or point estimates of incision relative to basalt-armored surfaces. Thermochronologic data are especially useful in regions of our model where basalt flows are absent, such as in the Canyonlands and portions of the southern Rocky Mountains. In areas where thermochronologic data complement basalt surface data, as near Grand Mesa, the two systems agree well with landscape evolution from depths of several kilometers to the present surface (see following). The thermochronometric studies used in this paper are documented in Table 2 and their locations are shown by numbered red circles in Figure 1, and indicated by large red circles on the map in Figure 2.


We divide the study area into subregions that correspond to specific features in the reconstructed surface and provide detailed documentation of data constraints related to that subregion. The data sources described here constrain a generalized model for the 10 Ma paleosurface, and this model is encoded numerically in the form of elevation control points that define the surface. It is important to understand that control points do not represent specific data measurements, but instead define what is known regarding topography of the 10 Ma surface. In regions where constraining data indicate a low-relief surface, control points are sparse, while regions with higher relief paleotopography have denser control points that define transitions in elevation. All control points have equal status in defining the surface, but may have different uncertainty associated with their value. This set of control points can be modified and improved with further work, and as new data become available, but they encode our present state of knowledge of the 10 Ma paleosurface.

All control points are shown in Figure 2 and are color coded and numbered by group; 13 groups are represented. White-rimed points in Figure 2 indicate basalt flows and larger red circles indicate sites of thermochronologic studies. Documentation for all control points is provided in Supplemental File 2 (see footnote 2), and includes location, elevation, and uncertainty. Uncertainty in eroded thickness at each control point is given as a 1σ standard deviation, but the uncertainties are a statement of our state of knowledge rather than measurement error, and so are subjective estimates of the combined analytical and geologic uncertainties. Nevertheless, it is important to provide some assessment of relative accuracy among the diverse data constraints used to define our model. We assign much lower uncertainty where the 10 Ma surface is preserved, but greater uncertainty where eroded thickness is inferred from thermochronologic constraints alone. In the following, we discuss data constraints for each group.

Group 1—Basin and Range

Black control points (Fig. 2) constrain elevations of the 10 Ma surface in basins of the Basin and Range province to the west and south of the Colorado Plateau. The Grand Wash fault that forms the boundary between the Colorado Plateau and Basin and Range was initiated by extension ca. 16.5 Ma that peaked by 13 Ma, and mostly ceased by 11–8 Ma (Faulds et al., 2001). It is likely that elevations of the Basin and Range relative to the Colorado Plateau were near present values by 10 Ma. The province has been characterized by internal drainage basins with little net transport of material from the subregion. Numerous basalt flows in the Basin and Range along the western margin of the plateau are near present basin elevations, indicating little erosion since emplacement (see Fig. 2; Supplemental File 2 [see footnote 2]). Some basins were integrated by the Colorado River ca. 5.3 Ma (Dorsey, 2010), but incision rates of the lower Colorado River are ∼20 m/Ma over the past 5.5 Ma (Karlstrom et al., 2007). Therefore, to a good approximation, net eroded thickness since 10 Ma for this subregion is <100 m. Control point elevations are equal to elevations of 10 Ma basalt flows (white rimed dots in Fig. 2), and equal to modern elevations of basins throughout this region (small black dots). Control points are located close to the margin of the Colorado Plateau to spatially define the transition in elevation to the plateau, and are located along the boundary of the study area to assure that the 10 Ma surface is ultimately interpolated over the entire study area.

Isostatic uplift due to extensional unloading of the margin of the Colorado Plateau likely developed when the Grand Wash fault was active ca. 17 Ma (Faulds et al., 2001). Thermochronologic data also indicate that exhumation by normal faulting in the vicinity of the Virgin River northwest of the present Grand Canyon was initiated from 22 to 17 Ma, and that 3–5 km of overburden was removed from 17 to 15 Ma (Quigley et al., 2010). The Vening Meinesz model predicts a maximum of ∼1 km of rift flank uplift within ∼50 km of the margin for effective elastic thickness of the crust of Te = 10 km, crustal thickness of ∼30 km, and a 60° dip normal fault (Watts, 2001). The locus of extensional slip along the western margin of the Colorado Plateau has expanded and migrated eastward as post–6 Ma slip has been observed on the Toroweap, Hurricane, Grand Wash, and Wheeler fault zones (Bohannen and Howard, 2002; Karlstrom et al., 2008). As these blocks are lowered into the Basin and Range, the rift flank uplift should migrate eastward, but the relative elevation between the western margin of the Colorado Plateau and Basin and Range would not change. Our interpretation is that the western margin of the Colorado Plateau has likely been at similar elevations relative to the Basin and Range since 10 Ma, but the margin may have migrated eastward, lowering some 10 Ma basalt flows toward the Basin and Range while lifting those farther east as the margin approaches them. These tectonic changes to the 10 Ma surface may alter elevations along the margin, but do not change our estimated regional eroded thickness relative to present elevations.

Group 2—Kaibab Plateau

Light green control points (Fig. 2) define the elevation of the 10 Ma surface on the Kaibab Plateau south of the Grand Canyon. By 10 Ma the Mesozoic section had been removed south of the present Grand Canyon and the paleosurface was near the top of the Kaibab Limestone (Flowers et al., 2008). An ∼2-km-high escarpment just north of the present Grand Canyon marked the southern extent of the remaining Mesozoic section (Lee et al., 2011, 2013). Control points include the 10 Ma Red Butte basalt near south rim, which overlies the Shinarump Formation of the Triassic Chinle Group ∼100 m above the Kaibab surface, and similar relationships where 6–10 Ma basalts on the Shivwits Plateau overlie lower Moenkopi strata. Thermal histories derived from both AFT and AHe data in the Grand Canyon suggest that this area cooled to below ∼20 °C by 10 Ma (Lee et al., 2013), indicating an eroded thicknesses of less than a few hundred meters. Based on preserved basalts, we have estimated an erosion rate of 10–20 m/Ma in the past 10 Ma. Control points are located to define elevations along the present canyon rim and plateau perimeter, and basalt flows and modern elevations are adjusted to the 10 Ma surface using the 20 m/Ma incision rate, resulting in ∼200 m of eroded thickness with an uncertainty of ∼100 m.

Group 3—Southern Rim of Colorado Plateau

Yellow control points (Fig. 2) define elevations along the narrow southern rim of the Colorado Plateau. This area is dotted with volcanic fields of Oligocene–Pliocene age and elevations of the 10 Ma surface are constrained by two ca. 8 Ma basalt flows (Fig. 2; Supplemental File 2 [see footnote 2]). The southeast basalt at 1690 m elevation is likely lowered by faulting on the southern rim of the plateau. Cather et al. (2008) interpreted the Chuska erg to have been aggrading just north of here (group 8) at 30 Ma at a minimum surface elevation of ∼3000 m relative to present elevations of ∼2400 m (Cather et al., 2008, their figs. 13 and 14, cross-section point C). This rim of the plateau was likely eroded as the basin formed in the present Little Colorado River drainage from ca. 26 to 16 Ma prior to deposition of the Bidahochi Formation. This suggests a long-term incision rate of 20 m/Ma, which we have used to adjust present elevations and the 8 Ma basalt elevations to the 10 Ma surface, resulting in ∼200 m of eroded thickness since 10 Ma. The uncertainty in eroded thickness was chosen to be equal to the eroded thickness to reflect the lack of constraints in this region.

Group 4—Rio Grande Rift

Red control points (Fig. 2) define elevations in the Rio Grande Rift. The rift was initiated ca. 36 Ma (Chapin and Cather, 1994) and is likely still slowly spreading and filling with sediment. Four basalt flows ranging in age from 8 to 11 Ma are perched above the current valley floor and indicate incision rates ranging from just 5 m/Ma south of the San Luis Valley to 38 m/Ma for the southern points (Fig. 2; Supplemental File 2 [see footnote 2]). Basalt flow elevations were corrected to the 10 Ma surface using their respective incision rates and ages. Control points used to define the perimeter of the rift had their present elevations raised by as much as 350 m to approximate the 10 Ma surface. In the San Luis Valley we have assumed that infill has kept pace with extension and there has been no net change in surface load since 10 Ma (zero eroded thickness), but this is poorly constrained and we have assigned an uncertainty of 300 m.

Group 5—Southern Rocky Mountains

Dark blue dots define the 10 Ma elevations across the southern Rocky Mountains (Fig. 2). Basalt and thermochronologic measurements in the southern Rocky Mountains form a complimentary picture of the 10 Ma surface. Age-elevation transects in the southeast Gore Range (Landman and Flowers, 2012), at Mount Sopris and Mount Ragged (Garcia, 2011), and at Huron Peak in the Needle Mountains (Feldman, 2010) indicate that modern valleys were filled to elevations of ∼3000–3500 m, and that 1 km of additional overburden covered present topography of the high peaks at 10 Ma, raising them to present elevations of ∼5 km (Karlstrom et al., 2012). The entire landscape was denuding at average rates of 100–200 m/Ma since 12 Ma based on both age-elevation transects and average incision rates (140 m/Ma) of the Colorado and Gunnison Rivers (Aslan et al., 2010).

Outside of areas affected by salt collapse, ca. 10 Ma basalt flows in the southern Rocky Mountains range in elevation from 2900–3400 m (Kirkham et al., 2001; Aslan et al., 2010). At Grand Mesa, Colorado, a 10.76 Ma basalt (Kunk et al., 2002) flowed down a paleovalley and buried gravels of Colorado River provenance at a modern elevation of 2935 m (Czapla and Aslan, 2009; Aslan et al., 2010). Other broadly correlative basalt flows associated with river gravels are present in the Flat Tops and at Basalt Mountain in western Colorado (Aslan et al., 2010).

The AFT and basalt data from the Flat Tops show good agreement; AFT data indicate that the surface had cooled through 110 °C from 60 to 30 Ma, and both 25 Ma and 10 Ma basalt flows are present at similar elevations, which indicates little change from 25 to 10 Ma at this location (Larson et al., 1975). In the MWX well near Rifle, Colorado, an AFT thermochronology age-elevation transect shows that it was buried ∼3 ± 1 km until 8 ± 2 Ma (Kelley and Blackwell, 1990), and was then exhumed at apparent rates of ∼200 m/Ma (Karlstrom et al., 2012), in agreement with AFT ages and incision rates in the Elk Mountains (Garcia, 2011) and at Grand Mesa (Aslan et al., 2011; Cole, 2011).

Taken together these data suggest that the now-eroded 10 Ma landscape in the southern Rocky Mountains would be ∼1 km above the modern terrain, with peaks at 5 km elevation, and 10 Ma basalts occupying paleovalleys at 3–3.5 km elevation relative to present. Based on this evidence our control points define ∼1 km of overburden across the highest peaks of the southern Rocky Mountains ca. 10 Ma, forming a broad region at 4.0–4.5 km present elevation with higher peaks reaching ∼5 km in the vicinity of Mount Sopris and the Elk Range, the southern Gore Range, and the San Juan Mountains. This terrain transitions into the piedmont of the Rocky Mountains (group 13) to the east, and into a broad low-relief surface extending west from Grand Mesa across the Uncompahgre uplift (prior to incision of Unaweep Canyon) toward the Tavaputs Plateau and the Canyonlands region (group 12). An uncertainty of 500 m has been assumed for all thermochronologic constraints in this region, whereas 50 m uncertainty is assumed for points constrained by basalt flows.

Group 6—Browns Park

Brown dots (Fig. 2) show control points on the 10 Ma paleosurface in Browns Park and the Green River Basin. North of the Uinta Mountains the 10 Ma paleosurface is defined by erosional remnants in the Green River Basin and east to the piedmont. The Gilbert Peak erosion surface is mantled by Oligocene Bishop Conglomerate, which represents the minimum elevation of the 10 Ma paleosurface. This assumption is supported by scattered outcrops of probable Oligocene or Miocene deposits in the region that locally overlie the Bishop Conglomerate. In addition, Pliocene–Pleistocene (younger than 2.5 Ma) lava flows of the Leucite Hills are present at elevations higher than the distal remnants of the Oligocene Bishop Conglomerate, and older lava flows occur at higher elevations than the younger flows. Collectively this information indicates that the 10 Ma surface was at least locally higher and certainly no lower than the elevations of the present-day Bishop Conglomerate and Gilbert Peak erosion surface remnants. Thus the model for the Green River Basin (Fig. 2) represents a minimum 10 Ma elevation. The timing of dissection of this surface is poorly constrained, and while we use elevations of erosional remnants to approximate the 10 Ma paleosurface, we assign an uncertainty of 300 m to this region.

Browns Park, southeast of the Uinta Mountains, is a half-graben that formed along the northern boundary of the eastern Uinta Mountains, and began to fill with Browns Park sediments by 25 Ma (Izett, 1975). Episodic faulting and syntectonic filling continued until ca. 8 Ma. The 10 Ma surface in Browns Park is constrained by two ashes with ages of 11.3 Ma and 9.1 Ma that are part of the graben fill (Izett, 1975; Luft, 1985) and have likely been lowered by subsidence. Uplands to the north and south of the Browns Park graben fill are veneered by Oligocene Bishop Conglomerate or denuded bedrock surfaces representing the Late Eocene–Oligocene Gilbert Peak erosion surface (Hansen, 1986). These uplands are ∼500 m higher in elevation than the ca. 10 Ma ashes, and represent the minimum elevation of the 10 Ma surface in this region. We consider this area to be poorly constrained and have added 500 m to present elevations of the ca. 10 Ma ash beds and assign a 500 m uncertainty to this region.

Group 7—Western Margin of Colorado Plateau

Orange dots show control points for the 10 Ma paleosurface along the western margin of the Colorado Plateau. Basalt flows on high ridges (Fig. 2) indicate that the surface ranged from 2900 to 3000 m present elevations, and we assume that these elevations have not been altered by faulting since 10 Ma. Elevations transition into the Basin and Range to the west, and are continuous with the low-relief surface extending east over the present Tavaputs Plateau to the southern Rocky Mountains. We use basalt elevations where present, and have raised present elevations along the top of the ridge that were not armored by basalt by 200 m to approximate the 10 Ma surface, assuming 20 m/Ma erosion rate for the low-relief surface. Where basalt flows are in place we have assigned an uncertainty in eroded thickness of 50 m, and 100 m elsewhere.

Group 8—Little Colorado River Drainage

Light blue control points (Fig. 2) define elevations in the drainage basin of the present Little Colorado River eroded from the Chuska erg to the level of the Bidahochi Formation by ca. 16 Ma (Cather et al., 2008). Present elevations of the Bidahochi Formation range from ∼1750 m to ∼2000 m. We have defined 1770 m as the approximate present elevation of the 10 Ma surface, and control points delineate the extent of this basin. Present elevations are as low as ∼1270 m, so our 10 Ma surface gives a maximum eroded thickness of ∼500 m. We assigned an uncertainty in eroded thickness of 50 m for this group since the Hopi Lake level is fairly well constrained.

Group 9—Escarpment North of Present Grand Canyon

Magenta control points (Fig. 2) delineate a paleoescarpment in the Mesozoic section that was just north of the present Grand Canyon at 10 Ma, as inferred from thermochronologic data (Flowers et al., 2008; Lee et al., 2011, 2013). The Mesozoic section likely extended from the top of the paleoescarpment north and east, and was continuous with the low-relief plateau across the Canyonlands region. This paleoescarpment retreated to the north over the past 10 Ma to become the present Grand Staircase (Escalante National Monument, Utah). We use the present Kaibab Limestone north of the Grand Canyon at ∼2000 m as the elevation at the base of the paleoescarpment in the Mesozoic section with uncertainty in eroded thickness of 30 m, and use 3000 m elevation for the top of the escarpment, consistent with elevations along the western margin of the Colorado Plateau. Uncertainty in eroded thickness is 200 m for the top of the escarpment.

Group 10—Interior Southern Colorado Plateau

Dark gray control points define the 10 Ma paleosurface as the extension of the Mesozoic section from east of the Kaibab uplift to the Rio Grande Rift. This surface is constrained to be ∼2 km thick over Lee’s Ferry (Lee et al., 2011, 2013), and should remain below the elevation of the Oligocene Chuska erg throughout the region. Between the Chuska Mountains and the southern Rocky Mountains the elevation of the 10 Ma surface is poorly constrained, but must be high enough to prevent the paleo–San Juan River from flowing south. The control points along the San Juan River and parallel points to the south provide this constraint. One 8 Ma basalt flow southeast of the Chuska Mountains has an elevation of 2188 m. There are numerous erosional remnants that must have been as high as or higher than their present elevation and serve as a minimum elevation for the 10 Ma paleosurface. These remnants have elevations that are tightly clustered between ∼2400 and 2700 m over the region. Maximum elevation of the Chuska erg was ∼3000 m, so we assume that 300–600 m has eroded since 30 Ma in this region, giving an erosion rate of 10–20 m/Ma. Eroded thickness is assigned an uncertainty of 100 m for this group, based upon uncertainty in erosion rate.

Group 11—Front Range of Southern Rocky Mountains

Salmon colored control points (Fig. 2) define a high topographic ridge that separates the Rio Grande Rift on the west from the eastern piedmont. A single ca. 12 Ma basalt flow is at >3600 m elevation. This ridge has likely been subjected to rift flank uplift associated with extension of the Rio Grande Rift. For a Te of 15–20 km (Lowry et al., 2000), the Vening Meinesz model predicts a maximum of ∼0.6 km of rift flank uplift within ∼70 km of the margin (Watts, 2001). The rift initiated ca. 30 Ma, so the rift flank uplift was likely established prior to 10 Ma and we have assumed no additional contribution to elevations since 10 Ma. Present elevations on the narrow ridge have been corrected to the 10 Ma surface using a conservative erosion rate of 70 m/Ma (approximately half of the ∼150 m/Ma in the southern Rocky Mountains based upon AHe data; Kelley et al., 2010). The 10 Ma surface is poorly constrained, so we have assigned an uncertainty of 500 m of eroded thickness for this ridge.

Group 12—Canyonlands of Utah

Dark green control points define elevations across the Canyonlands region derived using eroded thickness estimates from thermochronologic data. Just north of the present Grand Canyon, and including the Lees Ferry region, a 2 km thick Mesozoic section was still preserved at 10 Ma that extended north across Canyonlands, with an escarpment that may have trended east to the Chuska Mountains (Lee et al., 2011, 2013). North of the paleoescarpment, Mesozoic rocks were buried by Cenozoic deposits that extended north to the Tavaputs Plateau, and east to the southern Rocky Mountains (Hoffman et al., 2011; Kelley et al., 2010; Karlstrom et al., 2012). The 2 km of Mesozoic overburden at Lee’s Ferry was removed after 6 Ma (Lee et al., 2011, 2013), 1–2 km of overburden at the Monument uplift was rapidly eroded after 6–7 Ma (Hoffman et al., 2011), 2–3 km was removed from the region of the present confluence of the Green and Colorado Rivers after 4–5 Ma, and 1.2–1.9 km was removed south of the Book Cliffs near Green River, Utah, after 5–10 Ma (Hoffman et al., 2011). Our control points define elevations that result from adding these eroded thicknesses to present elevations at the sites of these studies, and they are consistent with our other groups in defining a high-elevation (relative to present), low-relief surface that extended from the Wasatch front on the western margin of the plateau east to the southern Rocky Mountains. This surface was presumably undulating and bowed up above numerous Oligocene intrusions, including the LaSal, Abajo, and Henry Mountains. The extent of any volcanic component to these laccoliths is unknown.

Group 13—Eastern Piedmont of Southern Rocky Mountains

Purple control points in the Great Plains define the elevation of the 10 Ma surface across the eastern piedmont. Deposition of the Ogallala surface continued ca. 18–5 Ma (Izett, 1975) in the north and 13–5 Ma in the southern portion of the study area (McMillan et al., 2002; Gustavson and Winkler, 1988; Hawley, 1993), so the 10 Ma surface is below and subparallel to the present undissected Ogallala surface. Control points were established for the 10 Ma surface using elevations on mapped erosional remnants of the Ogallala Formation and cross sections through the piedmont and mountain front to project the remnant Ogallala slope over eroded regions. There are 10 Ma basalts in the Raton-Clayton and Ocate volcanic fields at elevations just above the Ogallala unconformity, indicating limited erosion outside of the river drainages (Nereson et al., 2013). The surface is well constrained by basalt remnants, so we assigned 30 m as the uncertainty in eroded thickness.


Control points on the 10 Ma surface in Figure 2 have an irregular spatial distribution and were interpolated using a triangular network with linear interpolation. This method produces linear ridges that we subsequently smoothed for aesthetic purposes using a two-dimensional Gaussian filter with 1σ radius of 15 grid points. The numerical grid is defined in Universal Transverse Mercator (UTM) coordinates (Zone 12, NAD83; grid origin 102534 m, 3476024 m; grid increment 1250 m, 1250 m; grid dimensions 1024, 1108). The surface that results from interpolation is displayed in Figure 3 and supplied in Supplemental File 33, and represents the present elevation of the 10 Ma paleosurface. This paleosurface differs from present topography by the eroded thickness implied by all geological constraints defined herein and is coincident with erosional remnants of the paleosurface, such as the Grand Mesa. Overall the present elevations of the 10 Ma paleosurface are 1 km above the southern Rocky Mountains, ∼2 km above river levels in the Canyonlands region, and 1 km above the south-central Colorado Plateau. We have not specifically included localized volcanic intrusions or laccoliths like the Henry Mountains, or LaSal Mountains of Utah in our 10 Ma surface, because erosion of these Oligocene features is poorly constrained and their individual size is insufficient to make a significant contribution to isostatic rebound of the region.


The spatial resolution of the reconstructed 10 Ma paleosurface in Figure 3 is determined from evidence in the constraining data, but is generally a smooth representation of relative paleotopography, while modern topography has much greater detail. Subtraction of present topography from the 10 Ma surface produces an estimate of eroded thickness that is far more detailed than is justified by the resolution of the paleosurface (Fig. 4). Some choose to filter the topography to the same resolution as the surface before subtraction (Roy et al., 2009). However, it is possible, perhaps likely, that the 10 Ma landscape had lower relief and was smoother than modern topography, and that some difference in detail is valid. Since the ultimate goal is to use eroded thickness as a lithospheric load in isostatic deflection calculations, we choose to leave the eroded thickness as detailed as possible and let the isostatic rebound filter determine the final smoothness of the isostatic response. This approach will lead to increasingly good eroded thickness estimates as future data improve the resolution of the regional 10 Ma paleosurface.

Our estimate of net eroded thickness since 10 Ma is displayed in Figure 4 and provided in Supplemental File 44, and it is evident that widespread denudation has occurred across the Canyonlands region, the southern Rocky Mountains, and the western Great Plains. A striking feature of Figure 4 is the insignificance of the Grand Canyon compared to denudation of the overall region. Integration of both eroded thickness and our assessment of uncertainty over the Colorado River catchment, shown by the white boundary in Figure 4, gives the total volume of eroded material that was transported through the Colorado River system since 10 Ma, i.e., 3.4 ± 1.2 × 105 km3. Because formation of the Grand Canyon and most of the denudation of the Canyonlands followed river integration through the Grand Wash Cliffs 5–6 Ma (Spencer et al., 2001; Dorsey et al., 2007, 2011) we would expect the majority of our estimated eroded volume to have been transported to the Gulf of California. Dorsey (2010) derived an estimate of 2.8 ± 0.6 × 105 km3 for the volume of sediment delivered to the Gulf of California by the Colorado River since integration ca. 5–6 Ma. This result provides a constraint for all erosion models, and agrees with our estimate within the uncertainty. (See Dorsey and Lazear, 2013, for a discussion of sediment budget of the Colorado River since 6 Ma.)


Flexural isostatic response to erosional unloading is strongly dependent on the effective Te of the lithosphere. Lowry et al. (2000) derived an estimate of the spatial variation of Te for the southwestern U.S. using spectral coherence between topography and gravity that shows variations of a factor of 5 in the magnitude of Te over our region of interest. It is, therefore, important to use a laterally varying filter to derive the flexural isostatic response. Figure 5 shows a map of the laterally variant Te supplied to the isostasy algorithm. When Te, and therefore flexural rigidity D = E Te3/[12(1 – ν2)] (E = Young’s modulus and ν = Poisson’s ratio), is not constant, the fourth-order differential equation for flexure of the lithosphere w due to a surface load distribution q(x, y) has the form (van Wees and Cloetingh, 1994): 
where ρm = density of the mantle and g = acceleration of gravity. The second and third terms are equivalent to: 

The terms in Equations 2 and 3 differ from the form in van Wees and Cloetingh (1994) because they assumed D was a function of E and ν, while in our equations the independent variable is Te.

The importance of using variable Te for calculating isostatic rebound is illustrated in Figure 6, which shows the isostatic rebound filter for three points indicated by white dots in Figure 5. The Te and diameter of the filter at each location is indicated. The amplitude of the filter diminishes and the radius increases sharply as Te increases, so the filters have been scaled to allow comparison, with the height of the center filter multiplied by 5, and the height of the right filter multiplied by 10 relative to the left filter. The flexural rebound filter for Te of 7 km has a diameter of only 100 km, while the filter for 48 km Te spans the Colorado Plateau from the eastern Grand Canyon to Grand Mesa, Colorado (450 km). The filter test (Fig. 6) verifies that the spatially variant rebound filter is working properly, and the software was also tested using a constant Te and compared to a Fourier transform solution algorithm to verify that results were equivalent.

The distribution of lithospheric load q(x, y) in Equation 1 is derived from the spatially variant eroded thickness using density of eroded sediment and acceleration of gravity. The density of shale can range from 2100 to 2700 kg/m3 and that of sandstone ranges from 1900 to 2500 kg/m3 (Turcotte and Schubert, 2002). Assuming equal thickness of each, the mean density would be 2300 kg/m3, which is the value used for isostatic calculations herein (see Table 3 for all solution parameters).

Equation 1 was solved using second-order finite-difference numerical approximations on a discrete UTM grid (Zone 12) with 10 km sample intervals in both dimensions. This grid spacing was chosen to reduce the number of equations in the solution and is sufficient to provide ∼10 samples across the discrete deflection filter corresponding to the smallest Te. Eroded thickness on the 1250 m grid was anti-aliasing filtered using a 17 × 17 grid point matrix (20 km × 20 km) running average that has a notch at the Nyquist wavenumber for a 10 km sample interval, and eroded thickness was resampled to the 10 km grid for this solution. Parameter values used in the solution are given in Table 3. The system of equations was solved by direct matrix solution algorithms in Wolfram Mathematica (http://reference.wolfram.com/mathematica/tutorial/TheAlgorithmsOfMathematica.html) to provide deflection of the lithosphere produced by the surface load. The solution was then interpolated back to the original 1250 m grid. The small white rectangle in the lower left corner of Figure 4 shows the size of the anti-aliasing filter. This filter has no effect on the isostatic rebound solution because the rebound filter is much wider than that of the anti-aliasing filter.

The computed flexural isostatic rebound solution is shown in Figure 7 and provided in Supplemental File 55. There are two maxima in the isostatic response, one equal to 1.1 km near the confluence of the Green and Colorado Rivers in the Canyonlands region of Utah, and the other of 1 km north of the San Juan Mountains of southwestern Colorado. Superposition of the isostatic response of the southern Rocky Mountains and the eastern piedmont has produced >800 m of rebound at the head of the Arkansas River drainage, with 600 m of differential uplift along the drainage over a distance of 150 km east of the mountain front. Isostatic rock uplift in the western Grand Canyon is only ∼200 m, the eastern Grand Canyon is subparallel to the 300–400 m rebound contour, and the Grand Mesa in western Colorado is subparallel to the 800 m rebound contour. The northern Colorado Plateau and southern Rocky Mountains as a whole are shown to have undergone 800–1000 m of isostatic rock uplift since 10 Ma, resulting in >500 m of differential uplift of the upper Colorado River system relative to the section below the Lee’s Ferry knickpoint. Between Browns Park and the lower end of the Desolation and Gray Canyons, the Green River flows into ∼800 m of differential isostatic rock uplift in the downstream direction. This differential isostatic uplift within the Colorado River system has ramifications for superposition of the rivers on Laramide structures, integration of the Green and Colorado Rivers, paleo–drainage systems, the distribution of incision rates, and the shape of longitudinal profiles of the rivers. We discuss some of these ramifications in detail herein.


Two previous studies focused on isostatic rock uplift due to erosion of the Colorado Plateau (Pederson et al., 2002; Roy et al., 2009) and two others examined erosion across a broader region of the southwestern United States (McMillan et al., 2006; Cather et al., 2012).

The innovative study of Pederson et al. (2002) used modern elevations of erosional remnants of Late Cretaceous marine facies to estimate total rock uplift (tectonic and isostatic) since ca. 75 Ma. Where the facies are missing due to erosion, Pederson et al. (2002) estimated the thickness of the missing section to produce the elevation of this paleosurface relative to present elevations. The structural relief on the Late Cretaceous (post-Laramide) marine facies is ∼13 km, or equal to one-third the thickness of the crust, with a range in elevation from –5 km in the Uinta Basin to +8 km over the Uncompahgre uplift, and a mean elevation of ∼2 km that is the mean net rock uplift of the Colorado Plateau since 75 Ma. In addition, Pederson et al. (2002) found net erosion on the Colorado Plateau since 75 Ma to be just 400 m, and the corresponding Airy isostatic rebound was 308 m; no assessment of uncertainty was provided with these figures.

Pederson et al. (2002) also used preserved remnants of Oligocene–Eocene deposits to construct a ca. 30 Ma paleosurface; as before, where the contact is missing due to erosion, they estimated the thickness of missing section to produce the elevation of these paleosurfaces relative to present elevations. The 30 Ma surface was defined by elevations at 69 control points, 57 of which involved estimation of eroded section ranging from hundreds of meters to >2 km, so their control point elevations have large, but unspecified uncertainty. In Figure 8, our present (post-rebound) 10 Ma surface from Figure 3 has been cropped to the perimeter of the Colorado Plateau (Fig. 8A) for comparison with the Pederson et al. (2002) 30 Ma surface (Fig. 8B). The Pederson et al. (2002) surface shows elevations on the Kaibab Plateau that are close to modern elevations and similar to our 10 Ma surface, yet recent thermochronometric data suggest that 1–2 km of Mesozoic section was removed from the plateau between 30 and 10 Ma (Flowers et al., 2008), so we would expect a greater difference between these surfaces in this region. Likewise, the area of the Oligocene Chuska erg is shown at ∼2500 m elevation in Figure 8B; more recent geologic evidence indicates that the surface of the erg was at a minimum elevation of ∼3000 m at 30 Ma (Cather et al., 2008). Figure 8A also shows high elevations along the eastern margin of the plateau where it adjoins the San Juan Mountains. The San Juan volcanoes were actively building 38–29 Ma (McIntosh and Chapin, 2004) and likely reached elevations exceeding 5 km by 30 Ma, so our higher elevations on this margin reflect the transition into the remnants of this volcanic terrain at 10 Ma. In summary, our 10 Ma paleosurface (Fig. 8A) has several important refinements over previous work because it matches geologic data for erosion of the drainage basin where the Bidahochi Formation was deposited ca. 16–6 Ma, has higher elevations in the San Juan Mountains based upon AHe data for the southern Rocky Mountains, and has near modern elevations on the Kaibab Plateau reflecting the stripping of 1–2 km of Mesozoic section from the Kaibab surface between 30 and 10 Ma.

In spite of the differences, the two paleosurfaces are probably still more similar than would be expected, given the 20 Ma difference in modeled time frames. Pederson et al. (2002) found mean net erosion since 30 Ma to be 843 m, while the mean of eroded thickness shown in Figure 4 for just the Colorado Plateau is 800 ± 208 m. These values are identical within the uncertainty of our estimate. An important implication of this result is that there may have been relatively little net erosion between 30 and 10 Ma on the northern portion of the Colorado Plateau. Instead, much of the eroded thickness was removed from the northern plateau in the past 10 Ma, and perhaps mostly in the past 5–6 Ma, once the Colorado River became integrated to the Gulf of California through the Grand Canyon (Spencer et al., 2001; Dorsey et al., 2007, 2011). A broader result of this paper is that erosion in the past 10 Ma has not been limited to just the Colorado Plateau. Comparable magnitudes of erosion also occurred in the adjacent southern Rocky Mountains, so analysis of flexural isostatic response to erosional unloading of the lithosphere must consider the region as a whole.

Pederson (2006) and Pederson et al. (2007) extended their work to include flexural isostasy that defined differential rock uplift across the Colorado Plateau and found a bullseye of isostatic rebound centered near the region of the modern confluence of the Green and Colorado Rivers in the Canyonlands (Fig. 9B). Our solution for flexural isostasy using the extended Colorado Plateau and southern Rocky Mountains region shown in Figure 7 does not exhibit a strong bullseye, but instead shows a broad region of rebound from central Utah to central Colorado. Erosion in the southern Rocky Mountains produces isostatic rebound that is coupled to and superimposed with that of the Colorado Plateau and eastern piedmont, so the region must be treated in a unified manner and the Colorado Plateau cannot be treated in isolation. Figure 9A shows an isostatic rebound model that would result from cropping our eroded thickness estimate in Figure 4 to the boundary of the Colorado Plateau. This, in effect, assumes that the eroded thickness outside the plateau is zero, and artificially closes contours parallel to the plateau boundary. The result is very similar to the post–30 Ma rebound bullseye (Fig. 9B) computed from the Pederson et al. (2002) eroded thickness using a constant Te of 21 km and crustal density of 2300 kg/m3, and compares with the results of Pederson et al. (2007). When compared to the full rebound solution (Fig. 7, reproduced as Fig. 9C for comparison) the main change caused by truncating eroded thickness at the plateau boundary (Pederson et al., 2002; Roy et al., 2009) is to ascribe maximum rebound to a local region of the Colorado Plateau rather than a much broader region that includes the southern Rocky Mountains.

Roy et al. (2009) attempted to refine the estimate for total rock uplift on the Colorado Plateau since ca. 75 Ma derived by Pederson et al. (2002) by removing contributions from shorter wavelength Laramide deformation using a 2° × 2° smoothing filter, and removing the combined contributions from flexural isostatic uplift due to erosion and extensional unloading of the western margin of the Colorado Plateau. Their eroded thickness values for much of the plateau (Roy et al., 2009, their fig. 1B) were <1 km, consistent with the mean erosion of 400 m, and their rebound result since 75 Ma (their fig. 1C) shows the west side of the Colorado Plateau uplifted by ∼1 km relative to the east side. Their residual rock uplift since 75 Ma (Roy et al., 2009, their fig. 1D) was centered in the Henry Mountains–Marysvale region; they modeled this residual (nonisostatic) uplift as resulting from conductive warming of the mantle, mainly in the mid-Cenozoic. This is in contrast to more dynamic uplift mechanisms proposed by some, including a young (younger than 10 Ma) uplift component (e.g., Crow et al., 2011; van Wijk et al., 2010; Levander et al., 2011; Karlstrom et al., 2012). Roy et al. (2009) also limited their analysis to the Colorado Plateau, thereby excluding contributions to the isostatic rebound of their study area that arise from surrounding eroded thickness.

McMillan et al. (2006) used the projected tops of basin fill deposits to estimate erosion following post-Laramide aggradation, and were the first to examine post–10 Ma erosion of the Colorado Plateau, Rocky Mountains, and eastern piedmont; they mapped the highest erosional remnants of these formations from Nebraska and Wyoming to Arizona and formed an envelope elevation surface. Partial erosion of the preserved remnants and complete removal of these deposits across much of the Colorado Plateau resulted in much lower estimates of magnitude and distribution of eroded thickness than those of Pederson et al. (2002), or our results. McMillan et al. (2006) showed an isolated region of 1.2–1.4 km of erosion in the vicinity of Grand Mesa in western Colorado, 0.4–0.6 km of erosion across the Canyonlands, and 1 km of erosion in the San Juan River drainage; they recognized that their reconstruction represented a minimum estimate of eroded thickness, and that data sets such as elevations of Oligocene–Miocene volcanics indicated greater eroded thickness. McMillan et al. (2006) made no attempt to estimate the isostatic response to their eroded thickness.

Cather et al. (2012) summarized the erosional history of the southwestern U.S. and northwestern Mexico (long 98°W–116°W, lat 22°N–44°N), and examined the entire time frame from the Late Cretaceous (75 Ma) to the present. Cather et al. (2012) identified four erosional episodes, and the last episode from 6 Ma to present overlaps our study; they utilized the same recent AFT and AHe thermochronometry data as in our study. Comparison of our eroded thickness map (Fig. 4) with Cather et al. (2012; Fig. 8) indicates very good agreement over most of the overlapping region. A notable difference is found in southeastern New Mexico along the Pecos River, where our group 11 control points (Fig. 2) defined a high rift flank uplift ridge on the east flank of the Rio Grande Rift that graded across the now-eroded Pecos River drainage in the vicinity of the Jemez volcanic lineament. A 12 Ma basalt is at an elevation of ∼3600 m on this ridge. The elevation of this ridge and graded surface gave ∼1.5 km of eroded thickness at the mountain front, compared with 0.25–0.4 km estimated by Cather et al. (2012). Our result is dependent upon an assumption that elevations of the high ridge have not been altered by normal faulting or additional rift flank uplift since 10 Ma.

Similarly, our elevations along the Front Range near Denver, Colorado, exceed 4 km on the 10 Ma (post-rebound) paleosurface of Figure 3, and grade onto the eastern piedmont, leading to eroded thicknesses ∼2 km at the mountain front. This is double the estimate of Cather et al. (2012), although our thicknesses are in agreement with theirs farther east on the Great Plains. Our high eroded thickness results from our estimates of ∼1 km of exhumation across the southern Rocky Mountains since 10 Ma based upon AHe age-elevation transects from recent theses and papers (study areas indicated by red circles in Fig. 1; Garcia, 2011; Feldman, 2010; Hoffman, 2009; Landman and Flowers, 2012; McKeon, 2009). Cather et al. (2012) did not show any erosion estimates within the southern Rocky Mountains, and indicated that erosion along the Front Range must have been modest, based upon preservation of the Eocene Rocky Mountains erosion surface.


The timing and spatial distribution of erosion and the associated isostatic response to widespread post–10 Ma denudation that we have documented herein have far reaching implications for the geomorphology of the region, and provide important constraints on the formation and integration of the Colorado River system. We examine some of the ramifications of our results in the following. It is important to emphasize that, in our methodology, all elevations in all our figures are referenced to modern elevation, with no attempt to address models for any tectonic surface uplift that may have contributed to changing paleoelevations through time. Recent summaries of the ambiguity of paleoelevation studies in the region were summarized in Cather et al. (2012); until consensus and new data emerge, absolute surface elevation change remains poorly constrained.


To restore an estimate of the pre-erosion 10 Ma surface, we reverse the process of erosion and rebound by adding eroded thickness back to present topography, and apply the isostatic subsidence produced by the restored surface load. This process restores the 10 Ma paleosurface to its original pre-erosion topography if no other processes, like tectonic uplift, have altered elevations. If known tectonic uplift has occurred since 10 Ma it could also be reversed to obtain the original absolute elevations at 10 Ma, but such information is generally lacking. Figure 10A shows the restored 10 Ma paleosurface, which is also provided in Supplemental File 66. The high plateau across the Canyonlands in the post-rebound surface of Figure 3 has been lowered to the same elevation as the Green River Basin, and the rebounded 5 km peaks in the southern Rocky Mountains are isostatically lowered back to ∼4 km elevation at 10 Ma. The earlier (16–10 Ma) erosion of the Chuska erg in the present drainage of the Little Colorado River, and stripping of Mesozoic section from the Kaibab Plateau, have resulted in lower elevations on the southwestern margin of the Colorado Plateau. This lowered topography may have played a critical role in the evolution of the Colorado River system and is discussed further in the following.


As shown in Figure 10B, we speculate that consequent streams drained east off the Rocky Mountains to the Great Plains and deposited the Ogallala Formation from 18 to 5 Ma (Izett, 1975; Gustavson and Winkler, 1988), and that these streams had not yet integrated to form the Arkansas and Platte Rivers ca. 10 Ma. West of the continental divide streams may have flowed into several internally drained basins, and deposition was occurring in Hopi Lake (Bidahochi Formation) from 16–6 Ma (Cather et al., 2008), Browns Park from 25–9 Ma (Hansen, 1986), and the Rio Grande Rift since ca. 30 Ma (Chapin and Cather, 1994). We know that at the present location of the Grand Mesa the paleo–Colorado River was flowing southwest at ∼2200 m (relative to present) and likely flowed over the Uncompahgre uplift ca. 10 Ma prior to incision of Unaweep Canyon. Internal drainage in the Basin and Range province included the Muddy Creek Formation and proto–Lake Bonneville. The paleo–Rio Grande and Chama Rivers exited the plateau to the south and flowed into playa lakes in the southern Albuquerque Basin (Connell et al., 2005), and were probably not integrated to the Gulf of Mexico until the last few million years.

Other speculative paleodrainage features include (1) large internally drained lake areas on the Colorado Plateau that are a possible solution for the lack of integrated 10 Ma drainage from the Rocky Mountains to the Pacific (Blackwelder, 1934; Hunt, 1969), and (2) a paleocanyon in the location of the paleo–Little Colorado River (Karlstrom et al., 2012, 2013) connected to an ancestral eastern Grand Canyon suggested by thermochronologic data (Lee et al., 2013) and possibly connected to the Virgin River drainage, as suggested by Lucchitta (1989), and the Crooked Ridge river system (Lucchitta et al., 2011).


Some previous scenarios for formation of the Grand Canyon and erosion of the Colorado Plateau involved relief production at the Grand Wash fault that either (1) produced a transient knickpoint when the Colorado River reached that location 5–6 Ma that propagated upstream at a rate of 100 km/Ma to arrive only recently at its present location near Lee’s Ferry; or (2) initiated headward erosion of a drainage into the Colorado Plateau that propagated at 15 km/Ma and captured the Colorado River east of the Shivwitz Plateau 5–6 Ma (Pelletier, 2010). Neither scenario would account for the widespread denudation of the Canyonlands region during this same ca. 6 Ma time frame prior to arrival of the change in base level to the Canyonlands region.

Our smooth reconstruction of the 10 Ma surface in Figure 10 suggests that earlier (pre–10 Ma) erosion of the Kaibab Plateau and southwest Chuska erg (Cather et al., 2008) created an elevation drop of ∼600 m (from 2200 to 1600 m) relative to Lee’s Ferry that lowered the base level of the region south of Lee’s Ferry relative to the internally drained regions of the Canyonlands to the north, and provided an impetus for top-down integration of the Colorado River from the north. If a paleocanyon existed in the Mesozoic section in this region, as inferred from thermochronometry data (Lee et al., 2011, 2013; Karlstrom et al., 2013), this would have provided additional relief for lowering the base level. In this scenario, built upon the similar ideas of Blackwelder (1934) and Meek and Douglass (2001), rivers on the western slope of the Rocky Mountains flowed into internal drainage basins 10 Ma to form shallow lakes analogous to Hopi Lake to the south (Fig. 10B). These local drainage basins ultimately linked from north to south through lake spillover or groundwater sapping, and eventually became integrated to the lower base level south of Lee’s Ferry through escarpments of Mesozoic rocks that formed the paleo–Vermillion Cliffs and may have rimmed parts of Hopi Lake basin. The ∼600 m lowering of base level along this paleo–Colorado River could have triggered a transient knickpoint that would propagate rapidly upstream and widen across the Canyonlands, even as the river continued to integrate downstream. This scenario not only provides for erosion of the Canyonlands region prior to complete river integration, but also generates contemporaneous differential isostatic rock uplift from widespread denudation of the Canyonlands and the southern Rocky Mountains region that could drive superposition of the river system on Laramide structures, and maintain river gradients through Marble Canyon and the Grand Canyon following river integration.


On the restored 10 Ma surface (Fig. 10), the Tavaputs Plateau and Canyonlands regions of Utah have nearly the same elevation as the Green River Basin (∼2250 m relative to present), resulting in a low-gradient surface with playa lakes separated by broad divides. Such a surface would allow developing rivers to meander over the top of buried Laramide structures and around relic Oligocene volcanic edifices and domes. The Green River must have crossed the buried Tavaputs Plateau prior to onset of rapid denudation and isostatic uplift of the Canyonlands region ca. 6 Ma. If it had not been integrated with the Colorado River prior to ca. 6 Ma, uplift of the Tavaputs Plateau erosional remnant would have blocked its path southward. Integration of the Green River with the Colorado River must have taken place after ca. 8 Ma, which is the youngest age for Browns Park deposition (Luft, 1985), and before 1.2 Ma, the oldest dated Green River terrace (Darling et al., 2012). However, we suggest that integration took place between 8 and 5 Ma, and, as denudation progressed, the region was isostatically uplifted under the incising Green River to carve the Desolation and Gray Canyons.


On the restored 10 Ma surface of Figure 10, modern elevations for the top of Laramide structures have been lowered to pre-rebound elevations and the structures buried by the restored section. What today are widely ranging elevations of the Green River Basin (∼2000 m), Tavaputs Plateau (∼2800 m), Canyonlands of Utah (∼1200 m), and Monument uplift (∼2600 m) that would form barriers to river integration were isostatically adjusted at 10 Ma to underlie a low-relief surface on which rivers could meander over these structures. As the 10 Ma surface was denuded, the Laramide structures were isostatically uplifted into incising rivers, allowing entrenchment of meanders in deep canyons that characterize the Colorado Plateau. Table 4 shows modern elevations of the rivers and the elevation of the top of exposed Laramide uplifts transected by the Colorado River and Green River systems. It also shows the elevations of the corresponding points at 10 Ma (from Fig. 10).

Undoing isostatic rebound and adding back the now-eroded section places the rims of modern canyons below the modeled river elevation at 10 Ma. Only the top of the Kaibab uplift exceeded the elevation of the restored 10 Ma paleosurface, and we speculate that this is the result of an earlier scenario of isostatic rebound associated with erosion of Mesozoic section from the Kaibab Plateau and erosion of the Chuska erg between 26 and 16 Ma, prior to formation of Hopi Lake (Cather et al., 2008). We estimate that this period of erosion produced ∼500 m of isostatic rebound of the Kaibab uplift. The Kaibab uplift was likely breached by a 25–15 Ma paleoriver as it was isostatically raised under the paleoriver’s path, which is consistent with a paleocanyon inferred from thermochronometry data (Lee et al., 2011, 2013; Karlstrom et al., 2013) that incised into the Mesozoic section and was coincident with the course of the modern eastern Grand Canyon.

These data provide an ironic resolution and support for both sides of the century-long debate about antecedence versus superposition to explain the observation that rivers cut sharply across fault-bounded uplifts (Douglass et al., 2009). The superposition model (Davis, 1899; Blackwelder, 1934) is that meandering rivers flowed on a low-relief surface and were lowered by erosion onto stratigraphically buried block uplifts. This model is supported by our analysis. The antecedence model (Powell, 1875) is that mountains were uplifted into the path of established rivers. This model also has validity if the uplift is envisioned to be, in part, isostatic rebound resulting from regional scale erosional exhumation.


Any change in surface elevation is related to uplift of the Earth’s surface minus the thickness eroded from the surface (see Fig. 11A). For our purpose we divide uplift into a tectonic and an erosional isostatic component. We define the erosionally driven isostatic component as the rock uplift produced by elastic bending moments and buoyancy forces in the lithosphere in response to a change in loading of the surface, and we group all other sources of uplift into the tectonic component. This general relationship can be written 
where Δe is change in surface elevation, UT is tectonic rock uplift, UI is isostatic rock uplift from erosional unloading, and I is incision or eroded thickness (negative in sign for erosion, positive for deposition). The equation holds for both magnitudes and rates of uplift and incision.

In this paper we have estimated two of the four parameters in Equation 4. Our derivation of eroded thickness since 10 Ma (I) is shown in Figure 4, and the associated isostatic rebound (UI) is shown in Figure 7, both of which can be measured or computed independent of any tectonic uplift and surface lowering (relative to sea level) that has occurred during this time frame. Since both UI and I are known from our analyses, and they are opposite in sign, we define UI + I to be a residual incision. Then, change in elevation of the surface equals tectonic uplift plus residual incision (see Fig. 11A).

In Figure 12 we have displayed cross sections of our solution surfaces along the river corridor for the Colorado River (Fig. 12A) and Green River (Fig. 12B) to allow comparison of the relationships between these surfaces. Below the graphs we have included tomographic cross sections of P-wave velocity perturbations in the lithospheric mantle and upper asthenosphere (60–250 km depth) (Schmandt and Humphreys, 2010). The tomographic images depict relative P-wave velocity in the lithospheric mantle beneath the river profiles. Lower P-wave velocity (red) is related to lower density and higher temperature, and a 3% perturbation is equivalent to ∼500 °C. The Colorado River has higher temperature mantle below the southern Rocky Mountains (left), and below the Grand Canyon (right), with colder mantle underlying the central Colorado Plateau. Hotter mantle translates into buoyancy that can produce uplift, but surface elevation is also dependent upon crustal thickness and density as well as thickness and density of lithospheric mantle, so any correlation between hot mantle and elevation is not simple.

The longitudinal profile of the Colorado River (purple curve in Fig. 12A) traverses the southern Rocky Mountains and the Colorado Plateau through the Grand Canyon. Channel slope is steep in the Rocky Mountains, gentle on the Colorado Plateau, and steep again through the Grand Canyon. The blue curve is the cumulative magnitude of incision over 10 Ma, although most of this incision likely occurred after 6 Ma. The red shaded curve is the computed cumulative, erosionally driven, isostatic rock uplift. Residual incision (black curve) results from adding isostatic uplift to the negative incision, and represents the net lowering of surface elevation in the absence of any tectonic uplift (see Fig. 11). Compared to incision, which jumps from 1500 m in the Grand Canyon to >2000 m in the Canyonlands, residual incision is less variable, is approximately linear (dashed line), and (apart from the Rocky Mountains) becomes more negative in the downstream direction. The final two curves are the original (pre-erosion and pre-rebound) 10 Ma surface (gray), and the present elevation of that surface following erosional isostatic rebound. This is where erosional remnants of the 10 Ma surface, like the Grand Mesa, are found today. At first glance this figure is paradoxical, because the greatest erosion has occurred where the modern river gradient is minimal. However, the modern river profile has been shaped by accumulation of 1 km of isostatic rock uplift that diminished the slope of the river channel in the central Colorado Plateau as the river adjusted to isostatically uplifting bedrock. The steeper channel slope downstream of Lee’s Ferry through the Marble and Grand Canyons has previously been attributed to harder bedrock of the Paleozoic section (Pelletier, 2010) or tectonic uplift of the western Colorado Plateau (Karlstrom et al., 2008), but differential isostatic uplift of the Canyonlands relative to the lower Grand Canyon may also have influenced channel slopes of the modern Grand Canyon.

The Green River (Fig. 12B) is analogous to the upper Colorado River because it also flows into the 1 km of cumulative isostatic rock uplift at the confluence with the Colorado River. In contrast to the Colorado River, however, its longitudinal profile (purple) has almost uniform slope with little concavity. The residual incision (black) exhibits an overall linear trend (dashed line) and becomes more negative with increasing river distance. The tomographic cross section of the Green River exhibits cold mantle underlying the entire profile, with a shallow region of warmer mantle under the eastern flank of the Uinta Mountains and Browns Park. The striking differences in these two river profiles and in the mantle temperature perturbations suggest that the southern Rocky Mountains underwent Neogene tectonic rock uplift relative to the Wyoming Rocky Mountains.

Because the surface paleoelevation history of the Colorado Plateau region is still debated, we discuss two specific scenarios based upon our residual incision analysis. The first is the “old uplift” scenario, in which there was no tectonic rock uplift after 10 Ma, and the river and landscape have lowered over time (Fig. 11C). Figure 12A for the Colorado River shows that erosion has removed as much as 2 km along the river channel since 10 Ma, but isostatic rebound restores ∼1 km, leaving a net residual incision of 1 km. To end up at modern elevations the river would have initially started ∼1 km higher in elevation, analogous to Figure 11C. The Green River likewise has lowered a net of 1 km at the confluence with the Colorado River, but the residual incision diminishes approximately linearly upstream. For this model with no tectonic uplift, the paleoelevation of the region would have been approximated by the elevations shown in Figure 10.

A second scenario is a “stationary” scenario, in which there is no change in elevation of the river, and incision is matched by differential total rock uplift, including both erosional isostatic rebound and tectonic uplift, for example, driven by mantle flow and buoyancy (Fig. 11B). For this model, the residual incision curve is the negative of the tectonic uplift required to maintain a stationary elevation along the river, and tectonic uplift would have to increase approximately linearly in the downstream direction for both rivers. In this model, the 10 Ma paleoelevations of bedrock would have been several hundred meters lower than modern elevations in the upper Green River, ∼800 m lower in the Rocky Mountains, 1 km lower at the Green River–Colorado River confluence, and 1200 m lower in the Grand Canyon. In this scenario combined isostatic and tectonic uplift raised bedrock under an incising river that maintained constant elevation. Obviously, these are just two possibilities, and other combinations of surface elevation change and tectonic uplift can be envisioned. However, these two scenarios, in which first tectonic uplift is zero, and then elevation change of the river is zero, provide examples for interpreting residual incision.


McMillan et al. (2002), Leonard (2002), and Duller et al. (2012) examined tilting of the Ogallala Formation on the eastern piedmont and found that 1 km of rock uplift at the mountain front was required to explain the tilt. Leonard (2002) resolved this uplift into contributions from tectonic uplift of the southern Rocky Mountains, and isostatic uplift from dissection of the piedmont by the integrated east flowing Arkansas and South Platte rivers since ca. 6 Ma. Leonard’s analysis of isostatic rebound involving only the piedmont was insufficient to explain warping of the Ogallala Formation, so tectonic uplift of the southern Rocky Mountains was invoked for the remaining uplift.

We reanalyzed the Leonard (2002, fig. 1 therein) north-south transect at long 104.5°W using our erosion and rebound results (presented in Figs. 3, 4, 7, and 10) and assuming that our 10 Ma paleosurface is coincident with the Ogallala Formation on the piedmont. The location of this section is shown in Figure 1. We reconstructed Figure 2C of Leonard (2002) using slices through our solution surfaces, and our result is displayed in Figure 13. The graph is a north-south cross section (north at right) and zero distance is on the southern Colorado state line at lat 37°N for the transect. The brown curve shows a modern topographic surface that is broadly arched with peak elevations over the Jemez lineament, which crosses the transect near the Colorado–New Mexico border. The dotted curve approximates Leonard’s (2002) polynomial fit to the projected sub-Ogallala surface that is used as an estimate for the elevation of the now-eroded sub-Ogallala surface and shows upwarping across the Arkansas River drainage, but not across the South Platte River drainage. The purple curve is a slice through our 10 Ma paleosurface (showing post-rebound modern elevations from Fig. 3); it is nearly everywhere above the older erosion surface (black dots), and shows upwarping of the Ogallala Formation over both river drainages. The light purple interval between the blue and brown curves is eroded thickness along the sections. The red curve is a slice through our isostatic rebound solution (contours in Fig. 7), and the gray curve is a slice through the original (pre-rebound) 10 Ma surface (Fig. 10) plotted relative to the original depositional surface elevation in the transect of 1470 m assumed by Leonard (2002).

Upwarping of the present 10 Ma surface (purple) across drainages is not evident once rebound is restored (gray curve), so isostatic rebound accounts for the two smaller wavelength upwarps, but does not account for the broad long-wavelength arch. Leonard (2002) made the assumption that the original depositional surface had a north-south strike and east dip at 1.7 m/km, and should be a horizontal line in the plane of this transect. If this assumption is correct, the configuration of the gray curve represents tectonic uplift since 10 Ma relative to the original depositional surface, but the result is very sensitive to this assumption of north-south strike. Another possible interpretation for this arching would be if the Ogallala was deposited across a pre–10 Ma Jemez uplift that caused the strike of the mountain front to change from north-south to northeast in the vicinity of the Jemez lineament. If so, the gray curve (Fig. 13) might represent the trace of the original depositional surface in the plane of the transect. However, in a parallel transect farther from the mountain front at 103.5°W, Cather et al. (2012) showed a broad upwarp of the base of the Ogallala as it beveled across broadly folded Mesozoic units. In this parallel transect, the sub-Ogallala surface is also prominently warped over the Jemez lineament, suggesting that the lineament is the cause of the arching rather than the strike of the mountain front.

Our solution for isostatic uplift at the mountain front west of the transect (Fig. 7) shows ∼850 m of rebound along the Arkansas River compared to the ∼750 m of Leonard (2002), ∼500 m of rebound at the mountain front along the South Platte–Arkansas River interfluve compared to the ∼300 m of Leonard (2002), and ∼250 m of rebound at the mountain front along the Cheyenne Tablelands between the North and South Platte Rivers compared to the ∼200 m of Leonard (2002) and ∼170 m (McMillan et al., 2002). While our rebound values are higher due to our inclusion of isostatic uplift of the Rocky Mountains, they still do not account for the broad arching that suggests that a component of tectonic uplift is required.

The mantle tomographic image (Schmandt and Humphreys, 2010) displayed below the transect in Figure 13 shows that the lowest velocity (hottest) upper mantle underlies the Jemez lineament in the plane of the transect, and correlates with the highest tectonic uplift. Our conclusion from this analysis on the north-south transect is that combined isostatic rebound from erosion in the southern Rocky Mountains and eastern piedmont can account for the observed short-wavelength upwarp of the sub-Ogallala surface over drainages on this transect. Rebound does not, however, account for the broad 500–700-km-scale arching correlated with the Jemez lineament. Thus, our preferred interpretation is that mantle-driven tectonic forcings related to the Jemez lineament starting ca. 3–5 Ma (Nereson et al., 2013) produced broad arching in the southern Rocky Mountains and may have driven the change from 13 to 5 Ma aggradation and surface stability to net erosion in the Great Plains.


In the past 10 Ma the landscape of the Colorado Plateau and southern Rocky Mountains region has undergone a dramatic evolution involving increased incision and canyon cutting in the southern Rocky Mountains, integration of the Green River across the eastern Uinta Mountains to the Colorado River, integration of the Colorado River system across the Colorado Plateau to the Gulf of California, formation of the modern Grand Canyon, widespread denudation of the Canyonlands region of Utah, superposition of the Colorado River system onto previously buried Laramide structures of the Colorado Plateau, and development of significant relief on the previously aggradational low-relief Great Plains surface.

A detailed understanding of the spatial distribution of erosion and its associated isostatic response has emerged from our analysis of the Colorado Plateau and southern Rocky Mountains region since 10 Ma, and our analysis provides critical evidence that must be honored by any viable landscape evolution model. We find that 800–1000 m of isostatic rock uplift of the entire northern Colorado Plateau–southern Rocky Mountains–Great Plains region since 10 Ma has played a major role in the evolution of the landscape. The reconstructed pre-erosion and pre-rebound 10 Ma surface (Fig. 10) provides insight into initial conditions at 10 Ma that guided the future evolution of the region. Improvement in resolution of the relative timing of events such as onset of denudation of the Canyonlands, integration of the Colorado River through the Grand Wash cliffs, and onset of incision of the Ogallala Formation on the piedmont is necessary to further resolve the sequence of events and their driving mechanisms.

Aspects of this paper that are different from previous erosion and rebound investigations include the following. (1) The use of AFT and AHe thermochronologic constraints to infer thickness and timing of now eroded sediments has helped resolve the problem of lack of remnants of Tertiary aggradation surfaces across the region that hampered previous erosion studies, and has been crucial to our analysis. (2) Use of the entire Colorado Plateau–Rocky Mountains region for analysis encompassed all significant contiguous eroded thickness and prevented introduction of artifacts in isostatic rebound computations that result from truncation of surface loads at artificial boundaries such as the margin of the Colorado Plateau. (3) The extended region of analysis required consideration of the spatial variation of effective Te that ranged from <10 km to >50 km within the study area. Comparison of the area of greatest eroded thickness in Figure 4 with variation of Te in Figure 5 shows that Te varies much less (Te = 15–25 km) within the region of greatest erosion than within the extended region as a whole, but the dependence of flexural rigidity on Te3 makes these variations significant. (4) Consideration of erosion and its rebound consequences over 10 Ma rather than 30 or 75 Ma time scales is more appropriate for comparison to modern mantle images, and encompasses a period of rapid evolution of the northern Colorado Plateau. Erosion and rebound estimates from 30 and 10 Ma are similar on the northern Colorado Plateau, and imply that there was little net denudation from 30 to 10 Ma in that region.

Several areas of our model for the reconstructed 10 Ma surface remain poorly constrained and will likely require further research to improve our understanding of their contribution to isostatic deflection. These areas include (1) the Uinta Basin and Browns Park, where the 10 Ma surface has been tectonically altered and modified by deposition in a subsiding basin; (2) the Green River Basin, where the relationship between the 10 Ma surface and the Gilbert erosion surface is not clear; and (3) the San Luis Valley and Rio Grande Rift, where spreading and infill have been ongoing since 30 Ma (Chapin and Cather, 1994), and where the location of the 10 Ma surface and the impact of spreading and infill on isostasy in this setting is not well understood. The model presented here is a work in progress that can be updated and refined as new geologic and geochronologic information becomes available.

This paper benefitted from collaboration with the Colorado Rockies Experiment and Seismic Transects (CREST), which was funded by the National Science Foundation Continental Dynamics Program (grant EAR-0607808); additional support was provided by grants EAR-0711546 (to Karlstrom) and EAR- 1119629 (to Karlstrom and Aslan). Discussions with numerous CREST colleagues are appreciated, including Eric Kirby, Dave Coblentz, Will Ouimet, Laurie Crossey, and Rex Cole. Brandon Schmandt provided the tomographic cross sections from a tomographic model that benefitted from the densified CREST array. We thank Ken Dueker, Rick Aster, Jon McCarthy, and Steve Hansen for work on the CREST geophysical aspect. Kelin Whipple, Jon Spencer, and Joel Pederson provided detailed reviews that helped improve the paper.

1Supplemental File 1. Xyz file of modern topographic elevation on 1250m grid. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00836.S1 or the full-text article on www.gsapubs.org to view Supplemental File 1.
2Supplemental File 2. Excel table of control points for 10 Ma surface. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00836.S2 or the full-text article on www.gsapubs.org to view Supplemental File 2.
3Supplemental File 3. Xyz file of modern topographic elevation on 1250m grid. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00836.S3 or the full-text article on www.gsapubs.org to view Supplemental File 3.
4Supplemental File 4. Xyz file of eroded thickness since 10 Ma. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00836.S4 or the full-text article on www.gsapubs.org to view Supplemental File 4.
5Supplemental File 5. Xyz file of isostatic deflection for variable Te and density 2300. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00836.S5 or the full-text article on www.gsapubs.org to view Supplemental File 5.
6Supplemental File 6. Xyz file of original preerosion pre-rebound 10 Ma surface from variable Te. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00836.S6 or the full-text article on www.gsapubs.org to view Supplemental File 6.