Mantle-driven dynamic uplift of the Rocky Mountains and Colorado Plateau and its surface response: Toward a unified hypothesis

Buoyancy associated with mantle dynamics and flow is increasingly recognized as an important driver of topographic change in continental settings (Braun, 2010; Faccenna and Becker, 2010). Western U.S. examples include plumes (Yuan and Dueker, 2005; Wegmann et al., 2007), lithospheric drips and delaminations (West et al., 2004; Zandt et al., 2004; Levander et al., 2011), and asthenospheric upwellings (West et al., 2004; Wilson et al., 2005; MacCarthy, 2010; Levander et al., 2011). The elevation response in the western United States to asthenospheric buoyancy (i.e., the topography deriving from buoyancy variations beneath the lithosphere) is discussed by Lowry et al. (2000, plate 3d), who concluded that dynamic topography in the western United States contributes a significant fraction to the total elevation, particularly in the northern Basin and Range Province. Their study also concluded that dynamic uplift associated with the Yellowstone hotspot is effectively limited to the northern Basin and Range and that melt buoyancy is likely inadequate to generate significant topography in other regions, such as the southern Colorado Plateau and southern Rocky Mountains. In contrast, other workers have suggested that mantle dynamical processes are important in shaping dynamic topography in the Colorado Plateau–Rocky Mountain region (Karlstrom et al., 2008). Mantle velocity patterns have been modeled in terms of whole mantle flow (Moucha et al., 2008, 2009), predicting long-wavelength dynamic uplift of the Colorado Plateau. However, shorter-wavelength, subregional features seem to require upper-mantle flow, such as the inferred edge-driven convection in the southwestern Colorado Plateau region (van Wijk et al., 2010), downwelling below the Escalante region (van Wijk et al., 2010; Levander et al., 2011), or wet diapirs from the 410 km discontinuity, perhaps driven by slab-edge–stimulated upflow (Richard and Bercovici, 2009; Faccenna and Becker, 2010).

In this paper, we use the term “mantle-driven dynamic topography” in a broad sense to mean topography that is actively being modified by young (younger than 10 Ma) tectonism as driven by mantle convection (differential pressures at the base of the lithosphere) and/or changes in mantle buoyancy. Attempts to test and quantify the timing, magnitudes, and rates of mantle-driven dynamic contributions to surface topography are often hindered by the difficulty in deconvolving tectonic, geomorphic, and climatic drivers for surface modification. Complications also arise from incomplete understanding of lithospheric heterogeneity and buoyancy structure, multiple uplift events, and the temporal evolution of mantle dynamical processes. The approach taken in this paper is to test models for mantle-driven surface uplift by integrating diverse geological and geophysical data sets in the context of the continental-scale Colorado River system. This river system provides a sensitive gauge of the complex interplay among tectonic, geomorphic, and climatic processes that have shaped the western U.S. landscape (e.g., Powell, 1875; Riihimaki et al., 2007; Cook et al., 2009).

The surface elevation of the Colorado Plateau and Rocky Mountain region (Fig. 1A) underwent uplift from sea level during the past 80 Ma, with the Colorado Plateau now at an average elevation of ∼2 km (e.g., Pederson et al., 2002a), and the Colorado Rocky Mountains at an average elevation of ∼3.2 km (Fig. 1A). The overall 2–3 km of total surface uplift likely reflects multiple episodes that each modified the lithosphere. Hypotheses regarding mechanisms for these episodes include: Laramide flat-slab subduction between 70 and 40 Ma that hydrated and potentially mechanically eroded the lithospheric base (Humphreys et al., 2003; Liu and Gurnis, 2010); fragmentation of the Farallon slab ca. 35–25 Ma that induced additional uplift (Humphreys, 1995; Spencer, 1996; Roy et al., 2009); and post–10 Ma convective sinking and detachment of the lower lithosphere and asthenospheric upwelling, as well as other dynamic mechanisms, that may be providing support for high elevations (Karlstrom et al., 2008; Moucha et al., 2008, 2009; van Wijk et al., 2010; Levander et al., 2011). Multidisciplinary regional studies are needed to understand the relative importance of each of these uplift episodes and the tectonic, as well as the geomorphic and climatic, processes that shape elevated topography (England and Molnar, 1990; Zhang et al., 2001; Whipple, 2009).

This paper focuses on the Neogene history of mantle buoyancy and its potential contributions to the present-day elevations of the western U.S. orogenic plateau. We integrate present-day channel morphology, measures of fluvial incision and regional denudation over the past 0.5–10 Ma, and geophysical constraints on the structure of the mantle. Our goal is to elucidate the timing, scales, and perhaps also the processes of lithosphere-asthenosphere interactions through their surface manifestations. The first part of this paper briefly summarizes debates about the uplift history of the Rocky Mountains and Colorado Plateau, emphasizing an increasing recognition of both the multistage uplift history and the youngest, potentially still-active, component. The second part of this paper presents regional mantle tomography and various other geophysical and surficial parameters to evaluate spatial associations between mantle domains and physiographic features that may be an indication of ongoing mantle activity. The third part of this paper presents new data on the Colorado River system, with analyses of its longitudinal profile, incision history, and thermochronology. This approach, like the Colorado River itself, provides a link between the uplift histories of the Rocky Mountains and the Colorado Plateau. The overall goal is to begin to integrate emerging geological and geophysical data sets from the CREST project (Colorado Rockies Experiment and Seismic Transects) that lead us toward testable hypotheses for mantle-driven dynamic uplift of the Colorado Plateau–Rocky Mountain region and the magnitude of its surface response in the past 10 Ma.

The Colorado River is the single river system that drains the western slope of the southern Rocky Mountains. Since the work of Powell (1875) and Dutton (1882), numerous workers have recognized the close linkages between understanding the evolution of the Colorado River system and understanding the timing and mechanisms of uplift of the Rocky Mountain–Colorado Plateau region (see summaries in McKee et al., 1967; Hunt, 1969; Young and Spamer, 2001; Riihimaki et al., 2007; Karlstrom et al., 2011). This regional approach is also pursued in this paper.

For the Rocky Mountains, as well as for Grand Canyon, debates about the timing and mechanisms of uplift focus on the relative importance of “old’ versus “young” uplift episodes. Various workers have emphasized: Laramide (Epis and Chapin, 1975; Wolfe et al., 1998), mid-Tertiary (Eaton, 2008; Roy et al., 2009), and Neogene (Lucchitta, 1979; Karlstrom et al., 2005; McMillan et al., 2002, 2006; Leonard, 2002; Aslan et al., 2010) uplift and the ways in which each may have contributed to present-day elevations in the Rockies. Precursor Oligocene rivers existed in the Rockies (Hansen, 1965; Hunt, 1969; Larson et al., 1975; Cather et al., 2008), but their geometries remain poorly defined. The oldest dated upper Colorado River system deposits are the river gravels beneath the 11 Ma Grand Mesa basalt (Aslan et al., 2010). Given that upper Colorado River drainage networks were not integrated through Grand Canyon to reach the Gulf of California until after 6 Ma (Dorsey et al., 2007), pre–6 Ma western-slope rivers probably flowed into an internally drained lake system (Hunt, 1969) and/or may have exited the Colorado Plateau region elsewhere than through Grand Canyon (e.g., Lucchitta, 1990; Ferguson, 2011; Lucchitta et al., 2011). The Green River tributary to the Colorado River is also quite “young,” becoming a south-flowing Colorado River tributary between 8 and 1.5 Ma, as bracketed by internally drained basin fill of the Brown's Park Formation at 8 Ma (Aslan et al., 2010); and by the oldest dated Green River gravels at 1.5 Ma, which are inset below higher, and probably significantly older, undated terraces (Darling et al., 2011).

Similarly, for Grand Canyon, views of an “old” (Laramide and mid-Tertiary) southwest-flowing Colorado River system that drained the Rockies through an “old” Grand Canyon (Powell, 1875; Elston and Young, 1991; Hunt, 1969) tend to be associated with views of the dominance of Laramide and mid-Tertiary surface uplift episodes where the Rocky Mountain–Colorado Plateau region achieved near-present elevations soon after the Laramide (Gregory and Chase, 1992, 1994; Wolfe et al., 1998; Pederson et al., 2002a; Huntington et al., 2010; Wernicke, 2011), and/or during mid-Tertiary uplift (Spencer, 1996; Flowers et al., 2008; Roy et al., 2009). In contrast, models for a younger than 6 Ma Grand Canyon cite evidence that the Colorado River did not exit western Grand Canyon until after 6 Ma (Blackwelder, 1934; Longwell, 1936; Lucchitta, 1990; see summary in Karlstrom et al., 2008). “Old” (Laramide to mid-Tertiary) paleocanyons and paleoriver systems drained northeast, onto the Colorado Plateau (Young, 2001; Potochnik, 2001), and some paleocanyon segments may later have been re-occupied by rivers that flowed southwest in response to landscape changes and drainage reversal (Young, 2001; Flowers et al., 2008; Wernicke, 2011). Support for the concept of a “young” rather than an “old” landscape and canyon morphology, comes from: cosmogenically estimated modern erosion rates of ∼100–200 m/Ma (Cleveland et al., 2006; Nichols et al., 2011), cliff retreat rates of 0.5–6 km/Ma (Schmidt, 1989), geomorphic modeling (Pelletier, 2010), differential incision along the river system (Pederson et al., 2002b; Karlstrom et al., 2007, 2008; Darling et al., 2011), petrologic evidence for asthenosphere replacing lithosphere around the margins of the plateau during the past several million years (Crow et al., 2011), and geophysical and modeling studies for uplift due to mantle flow (Moucha et al., 2008, 2009; van Wijk et al., 2010; Levander et al., 2011).

Climatic influences also have affected relief generation and landscape evolution in the Colorado Plateau–Rocky Mountain region. Potentially important climatic events and changes include: the mid-Miocene warm period (15–17 Ma—Zachos et al., 2001; or 12–17 Ma—Chapin, 2008) and cooling temperatures since then (Fox and Koch, 2004); an early Pliocene event at ca. 6 Ma related to opening of the Gulf of California and intensification of the North American monsoon (Chapin, 2008); and a mid-Pliocene period (ca. 3 Ma) that was about ∼3 °C warmer and more humid than today (Thompson and Fleming, 1996; Haywood and Valdes, 2004). Cooling temperatures since ca. 2.7 Ma are proposed to have led to an intensification of glacial climates (Haug et al., 2005) accompanied by increased seasonality of climate and higher erosivity of geomorphic systems, perhaps globally (Molnar, 2004; Molnar et al., 2006). At shorter time scales, oscillating Quaternary glacial-interglacial conditions that caused variations in the extent of mountain glaciers are typically cited as the cause of cyclic changes in downstream bedrock incision versus aggradation recorded by river terraces (cf. Bull, 1991; Hancock and Anderson, 2002; but cf. Finnegan and Dietrich, 2011). Similarly, spatial variations in long-term bedrock incision from place to place (over the same time interval) need to be examined in terms of geomorphic causes such as river integration versus differential dynamic uplift.

At present, paleoelevation studies do not agree about the surface uplift history of the Colorado Plateau and Rocky Mountain region (cf. Gregory and Chase, 1992, 1994; Wolfe et al., 1998; Dettman and Lohmann, 2000; Sahagian et al., 2002a, 2002b; Huntington et al., 2010). Thus, an integration of indirect data sets is needed. This paper examines: (1) the geophysical state of the western U.S. mantle to evaluate any potential “young” forcings, (2) the depositional record of early Colorado River gravels to reconstruct incision histories back to ca. 10 Ma, (3) thermochronologic evidence for differential rock uplift in the region since the Laramide, (4) isostatic models that may help quantify the expected rock uplift that can be derived solely from isostatic response to erosional denudation, and (5) geodynamic models of the expected magnitude of surface response to different types of mantle flow.

Figure 1 shows a set of regional maps of the southwestern United States. Most of the individual maps present published regional data, and some incorporate new results in Colorado and northern New Mexico from the CREST experiment; the distinction between previously published and new data is made clear in the subsequent discussions. These maps provide an opportunity to compare a wide range of geological and geophysical data sets that manifest measurable parameters from different depths and hence provide an opportunity to examine potential mantle-to-surface linkages at the regional scale for the Rocky Mountain–Colorado Plateau region.

Figure 1A shows an image of topography for wavelengths >250 km that was created using a low-pass filter with a cosine-taper between 200 and 300 km. In this image, the size of the chosen filtering window is at the wavelength at which elastic plate filtering effects become small—roughly 250 km in the Colorado Rockies–Colorado Plateau region (Coblentz et al., 2011). This treatment allows us to examine topographic variation at a scale that is commensurate with variations in geophysical parameters of the deep crust and mantle. Figure 1A reveals topographic subprovinces such as the >2.7-km-high Colorado Rocky Mountains dome, elevated topography (∼2 km) around the edges of the Colorado Plateau, lower topography (∼1.5 km) within the interior of the Colorado Plateau, Great Plains, and northern Great Basin, and lowest topography in the southern Basin and Range Province. The continental divide is shown to follow the broad wavelength swell that extends from Yellowstone to southwestern New Mexico.

Figure 1B shows the P-wave velocity anomaly of the Rocky Mountain region as derived from teleseismic P-wave residual measurements using the ∼70-km-spacing EarthScope Transportable Array (TA) seismic array and all available PASSCAL experiment data (Schmandt and Humphreys, 2010). The tomographic image in Colorado is refined based on data from the ∼25-km-spacing combined CREST/TA array (Aster et al., 2009; MacCarthy, 2010; MacCarthy et al., 2011). In this image, very sharp velocity gradients are associated with the Lees Ferry knickpoint on the Colorado River (LF of Figs. 1A and 1B), the high-velocity Escalante region of the western Colorado Plateau (E of Fig. 1B), and the low-velocity features in western Colorado associated with the San Juan and Aspen anomalies and the Jemez lineament (SJ, A, and JL of Fig. 1B; Aster et al., 2009; MacCarthy, 2010). Tomographically imaged large velocity contrasts at 100 km depth in Figure 1B are up to 6% (Vp) and have corresponding contrasts of 11% (Vs). These contrasts are robust features of these regularized tomographic inversions, which, being preferentially smoothed, may underestimate true velocity variation amplitudes. Using standard velocity-temperature scaling relations (e.g., Cammarano et al., 2003; Jackson and Faul, 2010), these velocity variations correspond to as much as 500–700 °C of thermal contrast over lateral distances as short as 100 km. If a melt phase is present in the low-velocity domains (e.g., Sine et al., 2008; Schmandt and Humphreys, 2010), the thermal contrasts would be less, but rheologic contrasts could be enhanced. Several geodynamic models that scale the observed velocity variations to density structure and use reasonable mantle flow laws suggest that these large velocity contrasts may reflect Neogene and ongoing upper-mantle convective flow (Schmandt and Humphreys, 2010) of perhaps 1–5 cm/yr, which, in turn, could cause surface uplift of several hundreds of meters (Moucha et al., 2009; van Wijk et al., 2010).

Other geophysical data sets further suggest deep crustal and upper-mantle influence on physiographic provinces. Studies of the geoid can be used to investigate the depth and mechanism of isostatic compensation for high topography (Coblentz et al., 2011). Figure 1C shows the result of applying a band-pass spherical harmonic filter to the geoid anomaly data in order to examine buoyancy contributions from density contrasts in the upper mantle at depths between ∼40 and 400 km. Buoyancy variations at these depths, which produce characteristic surface wavelengths (at 40°N) in the range of 85–1750 km, are examined by filtering the geoid signal from degree and order between 17 and 360 (using the filtering methodology described in Chase et al., 2002). This filtered geoid is referred to as the “lithospheric geoid” (Chase et al., 2002), and it isolates the geoid anomaly derived from upper-mantle density contrasts. The lithospheric geoid in the Colorado Rockies region is characterized by a positive 5–10 m anomaly relative to the Colorado Plateau and the Great Plains. Analyses of the geoid anomalies is preferred to gravity analysis for evaluating the depth of compensation because they are sensitive to the distribution of mass with depth and are proportional to 1/r (where r is the distance to the compensating mass), whereas gravity is proportional to 1/r² (Chase et al., 2002). Thus, geoid anomalies “see” deeper than gravity anomalies. Based on geodynamic models that scale the tomographically observed velocity variation to density and buoyancy variation, a 5–10 m geoid anomaly in the Colorado Rockies relative to the Colorado Plateau and Great Plains, and a 4–5 m geoid high around the southwestern edge of the Colorado Plateau relative to the central Colorado Plateau and southern Basin and Range are compatible with ∼400–800 m of Neogene (and ongoing) surface uplift of the Colorado Rockies and western Colorado Plateau (van Wijk et al., 2010). These geoid data support the possibility that upper-mantle processes such as small-scale convection add a significant dynamic component to the support of the regional topography.

The area of highest topography in the Colorado Rocky Mountains is also underlain by an area of high crustal seismic attenuation (Fig. 1D). The crustal (Lg) attenuation (or, Q) provides additional information about tectonic processes that have an influence on the surface topography. Lg is the predominant phase on most short-period regional distance seismograms along continental paths and results from the superposition of trapped crustal shear waves that sample the bulk properties of the crust (see review and references in Phillips and Stead, 2008). In general, high attenuation (low Q) values are associated with tectonically active regions that may contain partial melt in the crust (Fan and Lay, 2002, 2003; Xie, 2002a, 2002b). Lg waves sample the bulk properties of the crust, such that a plausible interpretation of areas that have both high attenuation (low Q) and high geoid anomaly values is that active upper-mantle processes (as indicated by the geoid high) correlate with bulk crustal properties (low Q), which are manifested in elevated surface topography. The correlation between high attenuation (low Q) values and positive geoid anomalies is strong throughout much of the western United States, particularly in Yellowstone, the Colorado Rocky Mountains, Rio Grande rift, and western edge of the Colorado Plateau. The presence of low Q in these regions is likely related to crustal partial melt, a conclusion supported by the presence of Quaternary basalts in these regions. Higher Q is associated with more stable crust, such as the interior of the Colorado Plateau and Great Plains.

Crustal thickness estimates are shown in Figure 1E based on values from the “Crust 2” model of Mooney et al. (1998), with superimposed CREST array Ps moho image results (Hansen et al., 2011). The Colorado portion of the crustal thickness map is derived from a P-to-S wave receiver function data set that includes over 200 teleseismic events (Hansen et al., 2011). These data were migrated to depth using a three-dimensional (3-D) surface wave shear velocity model and assuming Vp/Vs of 1.76 using a common conversion point (CCP) stacking methodology. Moho depths were picked by identifying the maximum stack amplitude between 35 and 55 km depth, and the resulting picks are not smoothed. Moho depths beneath western Colorado range from 38 to 56 km, with 40–44-km-thick crust underlying the highest topography in the San Juan and central Rockies. The thickest crust (52–56 km) occurs as two areas in SW and NW Colorado (Fig. 1E). The association of thinner crust with regions of high topography (Fig. 1A), low upper-mantle velocity (Fig. 1B), large geoid anomaly (Fig. 1C), and high crustal seismic attenuation (Fig. 1D) is inconsistent with crustal thickness support via Airy isostasy. This conclusion supports previous conclusions presented by Sheehan et al. (1995), Lee and Grand (1996), Lerner-Lam et al. (1998), and Gilbert and Sheehan (2004), and it is significantly strengthened here by the improved resolution of crustal thickness from the CREST experiment. However, low-velocity crust beneath high topography is consistent with an important component of Pratt isostatic support, as discussed later herein.

Negative Bouguer gravity anomalies beneath the Colorado Rockies of ∼300 mGal (Fig. 1F) are the largest in the continental United States (Cordell et al., 1991; Isaacson and Smithson, 1976). The region has a near-zero isostatic gravity anomaly (Woollard, 1972; Marsh and Marsh, 1976; Thompson and Zoback, 1979), which has been used, along with early seismic-refraction interpretations (Pakiser, 1963), to argue for Airy compensation of the surface topography (Parsons et al., 1994). However, the large negative Bouguer anomalies are spatially associated with large granitic batholithic complexes, implying Pratt-style isostatic compensation. The northern anomaly overlies the Laramide to mid-Tertiary Colorado Mineral belt (McCoy et al., 2005), and the southern one overlies the Oligocene San Juan volcanic field (Lipman, 2007). These anomalies are also approximately collocated with the Aspen and San Juan mantle low-velocity anomalies and the region's highest topography. As mentioned already, crustal thickness estimates of ∼42 km beneath much of the highest topography of the Colorado Rocky Mountains (Fig. 1E) indicate the absence of any crustal root that was overthickened with respect to the adjacent Colorado Plateau and Great Plains, such as would be needed to explain the high topography via Airy isostasy. Gravity modeling suggests, however, that granite batholiths alone are insufficient to account for the gravity field (McCoy et al., 2005). Hence, we use the term “rootless Rockies” for the interpretation that there is not an adequate thickness of crustal root and that high topography is instead supported by a combination of buoyant crust and upper mantle (Sheehan et al., 1995; Hansen et al., 2011).

In summary, we interpret the geophysical data sets in Figure 1 as support for the hypothesis that the high-elevation region of the Colorado Rocky Mountains, as well as the high average elevations around the edges of the Colorado Plateau are supported by a combination of low-density crust and upper-mantle buoyancy. In the Rocky Mountains, the spatial correspondence of highest topography (Fig. 1A), lowest velocity crust and mantle (Fig. 1B), a 5–10 m geoid high relative to the Colorado Plateau (Fig. 1C), and the thinnest crust (Fig. 1E) indicates that both low-density crust and upper mantle contribute significantly to the support of high topography.

Multiple hypotheses to explain the timing and processes of buoyancy generation in the mantle include the following. (1) Mantle buoyancy associated with many or all of the low-velocity domains in Figure 1B may be due to ongoing convection involving asthenospheric mantle (imaged as low-velocity domains) moving up relative to stable or downwelling lithosphere (imaged as high-velocity domains). (2) Mantle buoyancy may be derived from older compositional domains embedded in the lithosphere. These hypotheses are not mutually exclusive given that “young” and “old” buoyancy contributions are undoubtedly operating in combination to produce today's mantle structure and surface topography. Nevertheless, the relative importance of each can be tested through geologic data because each hypothesis predicts different timing and processes of uplift in the Rocky Mountain–Colorado Plateau region. In the first case, a component of modern topography may be very young and dynamically supported by mantle flow (e.g., Karlstrom et al., 2008; Moucha et al., 2009; van Wijk et al., 2010; Levander et al., 2011; Coblentz et al., 2011). In the second case, topography may be in isostatic balance at longer time scales and may be mainly due to lithospheric buoyancy variations; for example, heterogeneities in lithospheric buoyancy may be caused by the juxtaposition of Precambrian terranes (Karlstrom et al., 2005), domains of differential Laramide hydration (Humphreys et al., 2003), and/or effects of spatially variable mid-Tertiary conductive warming of the lithosphere (Roy et al., 2009).

To test the extent and importance of young dynamic support versus older compositional support for topography, we examined direct and indirect measures of differences in late Cenozoic erosion across the Colorado Plateau–Rocky Mountain region. We evaluated whether patterns in erosional proxies such as the shape of modern river profiles, rates of fluvial incision along the Colorado River and its tributaries, and rates of exhumation determined from low-temperature thermochronology are correlated to mantle velocity domains. Present geodynamic models using reasonable mantle viscosities suggest that, to achieve significant (hundreds of meters) dynamic topography, upwelling of asthenosphere would likely have been active over millions to tens of millions of years (Moucha et al., 2009; van Wijk et al., 2010). If dynamic surface uplift has been active over these time scales through to the present, modern-day spatial heterogeneity of the mantle velocity structure beneath the region (Fig. 1B) would predict young surface response where regional river systems are underlain by low-velocity, perhaps upwelling, mantle (e.g., the Rocky Mountains) and highly variable response where rivers cross different mantle velocity domains (e.g., near Lees Ferry).

The hypothesis of dynamic topography (hypothesis 1) predicts that: (1) high topography, steep rivers, and high topographic roughness should correlate with low mantle-velocity domains, and (2) patterns of long-term river incision should change across mantle-velocity gradients and be highest above low-velocity domains. Such an association and associated modeling would strongly bolster the contention that both topography and mantle velocity structure are undergoing young modification (e.g., Karlstrom et al., 2008; Crow et al., 2011; Levander et al., 2011). In contrast, the hypothesis involving compositional buoyancy (hypothesis 2) would be supported if: (1) topographic and mantle anomalies can be explained by older compositional variations in the lithosphere (Roy et al., 2009; Lowry and Gussinve, 2011); (2) differential incision can be explained by geomorphic processes such as river integration, knickpoint propagation, and variations in rock erodibility, perhaps in combination with increased climatic erosivity in the past 3–4 Ma (Molnar, 2004), all acting on a previously elevated plateau (e.g., Pederson et al., 2002a); and (3) high topography and low-mantle-velocity domains more closely correspond to older geologic and magmatic features known to be of Precambrian, Laramide, and mid-Tertiary age.

Figure 2 shows the longitudinal profile of the Colorado River system (in blue) and its main tributaries plotted above river-parallel cross sections of the mantle velocity domains beneath the Colorado and Green Rivers. The Colorado River profile (Fig. 2B) has two generally concave-up segments demarcated at the Lees Ferry knickpoint at the head of Grand Canyon. This “double-concave” profile contrasts with profiles of equilibrium or graded rivers (Mackin, 1948; Whipple and Tucker, 1999), and it is morphologically similar to profiles of rivers experiencing a transient adjustment to a change in relative base level (Whipple and Tucker, 1999) and/or to differential deformation (Kirby and Whipple, 2001). Above their confluence, the Green River profile is significantly less steep on average than the Colorado River. The steep Rocky Mountain and Grand Canyon reaches overlie low-velocity mantle (Fig. 2E), whereas much of the Green River (Fig. 2D) and the Colorado Plateau reaches of the Colorado River are underlain by higher-velocity mantle (Figs. 2A and 2E).

The knickpoint between the two concave-up segments near Lees Ferry (LF) coincides closely with the center of a sharp mantle velocity gradient (gray band of Figs. 2B and 2C), but also with an upstream-dipping bedrock ledge formed by the resistant Kaibab Limestone, which forms the substrate of the Kaibab Plateau. This knickpoint has been variably interpreted to be a transient incision wave propagating upstream from Grand Wash Cliffs (Cook et al., 2009; Pelletier, 2010), a remnant of lake spill-over integration of the Colorado River across the Kaibab Uplift ca. 6 Ma (Scarborough, 2001), and a bedrock ledge separating harder bedrock below from weaker bedrock above that affects river gradient. A hybrid model suggests that an incision pulse is diffusively bypassing the harder bedrock (Cook et al., 2009; Darling et al., 2011). The collocation of the largest knickpoint on the Colorado River system with the large mantle velocity gradient (Fig. 2B) suggests the additional possibility of dynamic differential uplift to explain the position of this regional knickpoint, for example, where the P-wave velocity gradient reflects differential mantle flow that is manifested by differential surface uplift.

Colorado River tributaries near the Lees Ferry knickpoint also show knickpoints and convexities at similar elevation (Cook et al., 2009). This pattern and their lack of consistent correlation with high-strength bedrock are consistent with a transient incision pulse moving through a river system (Wobus et al., 2006; Berlin and Anderson, 2007). The markedly different magnitudes of the convexities in the major tributaries downstream (e.g., Little Colorado River) versus upstream (e.g., San Juan River) of the Lees Ferry knickpoint may reflect the progressive upstream propagation of this pulse of recent incision (Fig. 2B).

To compare differences in river profiles more quantitatively throughout the Colorado watershed (Fig. 3), we utilized an index of channel gradient normalized for contributing drainage area (channel steepness index is denoted as k_sn; Wobus et al., 2006). Fluvial longitudinal profiles typically exhibit a scaling known as Flint's law where channel slope (S) and drainage area (A) are related through a power-law relationship: where k_s is typically referred to as the steepness index and θ as the concavity. The use of a reference concavity (θ_ref is typically 0.45 or 0.50) is widely accepted to derive a normalized channel steepness index (k_sn), allowing for the direct comparison of channel gradients despite widely varying drainage area (Wobus et al., 2006). Channel steepness is expected to vary with uplift rate, climate, and substrate erodibility (Duvall et al., 2004; Kirby et al., 2003; Lague and Davy, 2003; Snyder et al., 2000; Whipple and Tucker, 1999). Recent studies have documented systematic increases in k_sn as erosion rates increase, for example, in landscapes exhibiting a relatively uniform precipitation and lithology (DiBiase et al., 2010; Ouimet et al., 2009). In these studies, channels are argued to be locally adjusted to the differential trunk river incision (Ouimet et al., 2009) and differential rock uplift (DiBiase et al., 2010).

Figure 3 provides a summary of the k_sn results for the Colorado River system. Topographic data used in the analysis are from the U.S. Geological Survey (USGS) National Elevation Data set (NED). We began with ∼30 m resolution NED data and resampled to an ∼90 m resolution, which provides sufficient resolution while providing a more manageable data set. The river profile derived from the DEM is shown in Figure 3A; dashed sections of the profile indicate where we have removed artificial steps associated with Lakes Granby, Powell, Mead, Mohave, and Havasu. For the entire Colorado River system, we analyzed all channels with drainage area >150 km² and calculated k_sn values for river segments of >10 km using a reference concavity of 0.45 (Figs. 3B and 3C). The interpolated k_sn map (Fig. 3B) is created using a moving window of radius 20 km to average surrounding k_sn values. We analyzed individual slope-area arrays for the main trunk rivers within the Colorado River system. Figure 3D is the slope-area array associated with the Colorado River main stem, with regressions through portions of the data that demonstrate a lower average k_sn value of ∼73 for upstream reaches (e.g., above Lees Ferry) that step up to values of ∼220 in lower reaches (e.g., through Grand Canyon). We also examined variation in k_sn (Fig. 3A) by binning average k_sn segments (30 km streamwise distance). This analysis shows that k_sn is moderate in the Rockies, steps up briefly in Cataract Canyon, returns to lower values until Lees Ferry, and then steps up to values of 200–300 throughout Grand Canyon. Alternate possible causative explanations for the steep reaches include bedrock strength, decreased sediment supply, one or more transient incision waves, and/or rock or surface uplift.

Interpolated to the regional scale (Fig. 3B), this analysis shows that anomalously steep normalized profiles are present in Grand Canyon and the Colorado Rocky Mountains, with lower-gradient reaches in the center of the Colorado Plateau. In contrast, the Green River is characterized by lower normalized gradients than the Colorado River for most of its length, except for discrete knickpoints. Both the Green and Colorado Rivers have long low-gradient reaches from well above their confluence downstream to Lees Ferry. As mentioned already, the steep Rocky Mountain and Grand Canyon river segments overlie low-velocity mantle, whereas much of the Green River is underlain by higher-velocity mantle (Fig. 2). Although the effects of sediment supply, grain-size distribution, and bedrock substrate on channel steepness are well documented (e.g., Sklar and Dietrich, 2001), both the Colorado and Green Rivers traverse a similar broad spectrum of rock types from Precambrian crystalline rock (Fig. 2) to Cretaceous shale, and they have similar tributary and sediment supply systems. Thus, it is likely that differences in substrate erodibility and sediment supply are secondary controls on river gradient, at least at the regional scale considered here (e.g., Darling et al., 2011). Furthermore, discharge estimates using both drainage basin size and historical records (Darling et al., 2011) suggest higher discharge in the Colorado Rockies relative to the Green River such that the steeper normalized gradient for the Colorado River cannot be easily explained by differences in climate or precipitation.

At subregional scale within the Colorado Plateau, Figure 3B shows steep reaches that correspond to knick zones where the river crosses uplifts of the Uinta Mountains (UM), Tavaputs Plateau and Desolation Canyon (DC), and Cataract Canyon (CC). Although beyond the scope of this paper, these zones probably reflect a combination of mechanisms such as harder bedrock (e.g., Precambrian Uinta Mountain Quartzite), salt tectonics in Cataract Canyon (Huntoon, 1988), transient waves of incision that have been hung up in slightly harder bedrock (Cook et al., 2009), and/or surface uplift associated with upper-mantle low-velocity domains impinging into the Colorado Plateau from its western margins (as shown in Fig. 1B).

To extend the geomorphic analysis from channels to hillslopes, we integrated information about topographic roughness (Fig. 3E). This uses an eigenvalue analysis of surface-normal vectors for a given surface area with a moving window of radius 20 km, as calculated from a 90 m resolution DEM (Coblentz and Karlstrom, 2011). Rough areas are those that have statistically anticlustered surface-normal distributions. The similar patterns seen in Figures 3B and 3E are consistent with the interpretation that steep and rough regions are associated with rapidly incising and tectonically uplifting portions of the fluvial system (Ouimet et al., 2009).

To test whether differences in the shape of river profiles along the Colorado River and its major tributaries are associated with differences in the rates at which these reaches have been incising into bedrock, we compiled available data on incision rates over the past 10 Ma. These are summarized in Figure 4. Differential incision rate data above and below knickpoints offer another measure of the state of equilibrium or nonequilibrium in a river system. The location of a knickpoint potentially marks the boundary between the adjusted and unadjusted reaches of the channel (e.g., Wobus et al., 2006), with higher incision rates expected in the still-adjusting reaches below knickpoints. Figure 4A shows bedrock incision rates (m/Ma) and their locations (Fig. 4B) over time spans that range from 0.5 to 11 Ma The magnitude of incision (measured from a gravel-capped strath to the modern river level) and duration (Ma) are tabulated in the GSA Data Repository (Table DR1).¹ Data reported here were derived primarily from Ar-Ar dating of basalt flows that overlie river gravels in basalt-capped plateaus (inverted topography), U-Pb dating of speleothems from caves within Grand Canyon (Polyak et al., 2008), cosmogenic burial dating of gravels (Wolkowinsky and Granger, 2004; Darling et al., 2011), and occurrences of the 640 ka Lava Creek B ash (Andres et al., 2010). We use time spans of >0.6 Ma to minimize complexities in incision rate patterns due to glacial-interglacial cycles (Pederson et al., 2006; Aslan et al., 2010), as well as potential complexities associated with transient incision (Cook et al., 2009).

Working upstream from the Gulf of California, Figure 4 shows that incision rates downstream of Grand Canyon are very low and that the river has been within a few hundred meters of the same level since 5.5 Ma (House et al., 2008). In western Grand Canyon, between the Grand Wash Cliffs and the Hurricane fault, average incision rates of ∼78 m/Ma over the past 3.9 Ma have taken place over about the same time interval that the fault block has been lowered relative to the Colorado Plateau at a rate of 100 m/Ma, implying fault-dampened river “lowering” rates of 178 m/Ma. Eastern Grand Canyon incision rates are estimated at 233 m/Ma over the past 3.7 Ma (Karlstrom et al., 2008). Long-term incision rates above the Lees Ferry knickpoint are lower: 102–126 m/Ma over 1–2 Ma (Darling et al., 2011; Wolkowinsky and Granger, 2004). These are markedly slower than incision rates below Lees Ferry, especially if the published shorter-term (maximum) rates above Lees Ferry of >300 m/Ma in the past 0.5 Ma are realistic estimates of the more recent incision pulse in the Lees Ferry region (Cook et al., 2009; Darling et al., 2011).

Rates of differential incision between the Colorado and Green Rivers also vary strongly (Fig. 4). The Green River system has incision rates of 4–100 m/Ma over the past 12 Ma and relatively low incision magnitudes of <150 m over the past 10 Ma. In contrast, the Colorado River has rates of 133–164 m/Ma over the past 1–11 Ma, and the total depth of incision approaches ∼1500 m over this time period (Aslan et al., 2010). The observation of higher incision rates and magnitudes for the Colorado River relative to the Green River, coupled with the Colorado's steeper normalized channel gradients and higher discharge, is counter to what is expected for equilibrium fluvial systems. Rather, higher discharge (e.g., due to higher precipitation and climate patterns) would be expected to result in lower equilibrium channel gradients. Thus, our preferred interpretation is that the low-velocity mantle underlying the Colorado River in Colorado is driving differential rock uplift of the Colorado Rocky Mountains relative to the central Colorado Plateau, and that the observed increased steepness, roughness, discharge, and incision rates and magnitudes of the Colorado River relative to the Green River reflect adjustments to this mantle-driven uplift.

To evaluate patterns in longer-term exhumation, we plotted results of apatite fission-track analysis (AFT; Fig. 5A) and apatite (U-Th)/He analysis (AHe; Fig. 5B) from across the region (Fig. 5C). While not an exhaustive summary, the purpose of this figure is to show representative data sets that depict the Colorado Plateau–Rocky Mountain region's multistage regional exhumation history since the Laramide as recorded in transects that span present-day elevations from −1000 to 4000 m. These data extend denudation constraints back in time from the direct measurement of incision rates. For example, the heavy black lines on Figures 5A and 5B show the long-term average bedrock incision rate for the Colorado River based on the Grand Mesa basalt since 11 Ma (Aslan et al., 2010), and the incision rate needed to carve much of Grand Canyon in the past 6 Ma.

Figure 5A summarizes AFT data. Linear segments of the age-elevation transects can be used to estimate the “apparent exhumation rate” of rocks through the paleo–110 °C isotherm if one assumes there was a relatively simple thermal history in the transect region for the time period recorded by the ages. Although direct interpretation of exhumation rates from the slope of age-elevation data may be complicated by topographic and structural perturbations of near-surface isotherms (Stuwe et al., 1994; Ehlers et al., 2001) and by advection of heat during rapid exhumation (Moore and England, 2001), the relative simplicity of age-elevation data from the Rockies (Fig. 5A) and the consistency among different regions lead us to argue that cooling largely reflects regional denudation. Use of age-elevation traverse data to estimate exhumation rates has the advantage of being independent of assumed geothermal gradient, as long as the region has a semisteady isotherm structure for the region and time span under consideration. Also, comparing age-elevation data from adjacent regions may provide insight into the timing of faulting, as shown across the Colorado Mineral belt (Kelley and Chapin, 2004) and Laramide and mid-Tertiary faults in southern Wyoming (Steidtmann et al., 1989; Kelley, 2005).

The White River uplift data (WR of Fig. 5A) is representative of many areas in the Rocky Mountains that cooled through AFT closure temperatures of ∼110 °C during the Laramide, 70–40 Ma (see also Naeser et al., 2002; Kelley and Chapin, 2004). Apparent Laramide exhumation rates here were ∼60 m/Ma. In contrast, the western Gore Range cooled systematically from 35 to 15 Ma, with apparent exhumation rates of 40 m/Ma, and the eastern Gore Range cooled from 20 to 6 Ma, with apparent exhumation rates of 80 m/Ma (Naeser et al., 2002). These data suggest post–20 Ma faulting/tilting of the eastern side of the Gore Range and resulting differential exhumation due to Neogene extension in the region where the Colorado River crosses the Gore Valley (Naeser et al., 2002). Similarly, parts of the Grand Canyon region cooled through 110 °C in the Laramide, as shown by AFT age-elevation traverses in eastern Grand Canyon (Naeser et al., 2001; Kelley et al., 2001).

AFT data from the MWX well from along the Colorado River north of Grand Mesa (Fig. 5A) also show a kink in the AFT ages that marks the base of a 10 Ma partial annealing zone, indicating paleotemperatures (at 10 Ma) of ∼110 °C, corresponding to paleodepths of ∼3 km (Kelley and Blackwell, 1990). The base of this 10 Ma partial annealing zone is now at ∼33 m elevation with respect to mean sea level, ∼1.5 km below the present surface, indicating ∼1.5 km of exhumation since 6–10 Ma. In this well, shallow samples record landscape stability and the development of an AFT partial annealing zone during middle Cenozoic time, whereas deeper samples record a sharp increase in apparent exhumation rate to ∼200 m/Ma starting ca. 6–10 Ma (Kelley and Blackwell, 1990). Onset of rapid exhumation at ca. 6–10 Ma in some high-elevation regions is also recorded by an AFT age-elevation traverse from the northwestern San Juan Mountains (Kelley and McKeon, 2009).

Figure 5B shows emerging AHe age-elevation data sets that record progressive cooling through temperatures of 50–70 °C. Much of the Grand Canyon region was affected by a mid-Cenozoic cooling history (28–16 Ma) that followed a period of relative stability from 50 to 30 Ma (Flowers et al., 2008), and similar mid-Cenozoic cooling is observed in parts of the Rockies (e.g. Twin Lakes batholith of Feldman, 2010, and parts of the Elk Mountains of Garcia, 2011). A factor of direct interest to testing Neogene mantle forcings, late Cenozoic (10 Ma to present) AHe cooling through ∼50–70 °C is recorded in many areas of the Colorado Rockies and Colorado Plateau. Neogene exhumation seems to have occurred at a wide variety of elevations. Several areas show high post–10 Ma apparent exhumation rates: Mount Sopris and Ragged Mountain of the Elk Range (Garcia, 2011), Bear Peak of the northwestern San Juan Mountains (SJ; McKeon, 2009), Huron Peak of the Twin Lakes batholith area (Feldman, 2010), and Monument Uplift of the Colorado Plateau (Hoffman, 2009; Hoffman et al., 2011). These areas give high apparent exhumation rates of >150–200 m/Ma between about 12 and 5 Ma.

Thermochronology-derived exhumation rates of ∼150–300 m/Ma for the past 6–10 Ma from the MWX well and from several of the AHe transects, and erosion magnitudes of 1.5 km at MWX agree well with geologically derived long-term average river incision rates of 150 m/Ma and incision magnitudes of 1.5 km since eruption of the 11 Ma Grand Mesa basalt along the upper Colorado River (black line of Fig. 5B; Aslan et al., 2010). Additional data are needed to refine our knowledge of the timing of rapid onset of incision. Alternative mechanisms to explain increased denudation include: (1) integration of the Colorado River system ca. 6 Ma, (2) a climate change event ca. 6 Ma (Chapin, 2008), and (3) mantle-driven surface uplift due to mantle forcings starting 6–10 Ma. Based on existing data, the rapid onset of exhumation at 6–10 Ma seems to have predated the integration of the Colorado River through Grand Canyon that took place some 500 km downstream at 5–6 Ma (Lucchitta, 1990) such that downstream base-level fall probably does not explain the timing of onset of rapid exhumation in the Colorado Rockies. However, lake spillover and rapid downstream integration at 5–6 Ma remain possible, such that more precise timing constraints for the onset of rapid exhumation seen in the different data sets are needed to resolve the tectonic from geomorphic drivers that have contributed to post–10 Ma denudation.

It is important to distinguish different types and drivers of “uplift.” Simple isostatic considerations suggest that, in the absence of tectonically driven uplift and for broad loads that are locally compensated, a given amount of denudation (D) leads to surface lowering of ∼0.2D and 0.8D of rock uplift. When the scale of topographic loads is less than several hundred kilometers, however, the surface load is partially supported by elastic strength of the lithosphere, and the magnitude of rebound is diminished. Consequently, denudation on a lateral scale of hundreds of kilometers will produce a strong component of rock uplift due to isostatic rebound, whereas smaller-scale incision of narrow canyons results in little rock uplift due to isostatic response.

To evaluate the isostatic component of rock uplift in the Colorado Plateau–Rocky Mountain region resulting from widespread denudation since 10 Ma, we developed a set of flexural models based on a reconstructed 10 Ma paleosurface defined mainly by preserved basalt remnants, but also incorporating estimates of eroded thickness from thermochronology (Figs. 6A and 6B) and other constraints (Fig. DR2; Table DR2 [see ^{footnote 1}]). This modeling follows methodologies somewhat similar to other studies (Pederson et al., 2002a; McMillan et al., 2006; Roy et al., 2009), with the following differences: (1) We restrict the time period for estimating exhumation to the past 10 Ma (instead of post-Oligocene; Roy et al., 2009) based on the new thermochronology evidence for onset of rapid denudation across the region since 10 Ma. (2) We use a broader region, extending from the Great Plains to Nevada, and from Wyoming to Mexico, to reconstruct eroded thickness and the associated flexural isostatic rebound. This allows us to incorporate erosion of the southern Rocky Mountains and the piedmont of the Great Plains in our rebound models, and it provides a peripheral zone of data that removes edge effects of the ∼150-km-radius rebound filter from the regions of interest on the Colorado Plateau. In addition, it allows us to compare our results with evidence for tilting and deformation of Miocene strata on the eastern piedmont due to both recent uplift of the Rocky Mountains and isostatic rebound from erosion (McMillan et al., 2002; Leonard, 2002). (3) Our estimates of denudation magnitude since 10 Ma are constrained in regions around the margins of the Colorado Plateau by the abundance of ca. 10 Ma basalt flows that preserve remnants of the paleolandscape (Fig. DR2 [see ^{footnote 1}]). In areas without basalt paleosurfaces, we used the low-end values of eroded thickness estimates from the newly available thermochronology data discussed previously herein. In addition, we used geologic evidence for location and timing of events on the southern and western Colorado Plateau that define initial conditions for the 10 Ma paleosurface. These events include erosion of the Grand Canyon region to the Kaibab surface prior to 10 Ma (Flowers et al., 2008), development of the Basin and Range from 16 Ma to 10 Ma (Faulds et al., 2001), erosion of the Chuska Sandstone erg to near-modern elevations in the present Little Colorado River drainage by 16 Ma, and formation of Hopi Lake from 16 Ma to 6 Ma (Cather et al., 2008).

Elevations of the 10 Ma surface, measured relative to present elevations, were defined by control points (Fig. DR2 [see ^{footnote 1}]) that honor the eroded thickness constraints imposed by these data. Control points were interpolated to form a surface, and current topography was subtracted from this surface to produce the estimate of net eroded thickness shown in Figure 6A. The key features of this map are that ∼1500 m of denudation have occurred along the corridor of the Colorado and Green Rivers from the flanks of the Rocky Mountains to eastern Grand Canyon. On either end of the river system, incision is confined to narrow canyons that would elicit little isostatic response, while the Canyonlands region has experienced widespread denudation. Eroded thickness was treated as a surface load distribution, and we computed the flexural response to this load via convolution in the wave-number domain. Contours of isostatic rebound are displayed in Figure 6A overlying the map of eroded thickness. Our results are broadly similar to the isostatic rebound modeled by previous researchers, in spite of our shorter time frame, which suggests that denudation rates were very low prior to 10 Ma.

Figure 6B is a cross section through the 3-D isostatic solution, taken along the Colorado River channel. It shows that estimated rock thickness removed by erosion since 10 Ma (black line in Fig. 6B) is up to ∼1800 m near the Green River confluence, decreasing to lower values in the uppermost and lowermost parts of the river system. Since thermochronometry data suggest ∼2 km of exhumation near the Green River confluence since 4–5 Ma (Hoffman et al., 2011), our results are minimal estimates for eroded thickness. Effective elastic thickness also varies from 40 km in the plateau center to ∼10 km at each end of the river system (Lowry et al., 2000; Lowry and Gussinve, 2011), but our present calculations used a mean elastic thickness of 25 km. Calculated isostatic rebound (gray area, Fig. 6B) of up to 1 km in the center of the Colorado Plateau is ∼55% of total post–10 Ma incision (compared to 80% for Airy isostasy). This indicates that the eroded thickness was partially supported by elastic strength in the lithosphere. Likewise, isostatic rebound in Grand Canyon and southern Rocky Mountain regions amounts to 25% of post–10 Ma exhumation. Calculated isostatic rebound (gray line, Fig. 6B) of up to 1 km in the past 10 Ma in the center of the Colorado Plateau can account for ∼50% of total rock uplift since the Cretaceous. In contrast, most of the rock uplift in the Grand Canyon and southern Rocky Mountains regions must be explained by other mechanisms than isostatic rebound. The red curve in Figure 6B (exhumation–rebound) indicates that there is a greater discrepancy between total incision and isostatic response in Grand Canyon and Rocky Mountains than in the central Colorado Plateau. This would be consistent with the hypotheses that tectonic uplift may have been greatest at each end of the profile, for example, in areas overlying low P-wave velocity in the mantle. These results are similar to the “residual rock uplift” amounts of Roy et al. (2009), but we include additional data from the Rocky Mountains, and we envision a more dynamic causation than their conductive heating model.

Rebound models need to be interpreted with caution. Our model, as well as previous models (Pederson et al., 2002a; Roy et al., 2009), does not have firm constraints on absolute paleoelevation except for Cretaceous and present times (cf. Gregory and Chase, 1992; Proussevitch et al., 2007; Sahagian et al., 2002a, 2002b). Models that view isostatic response as a primary mechanism for rock uplift implicitly assume that mean elevation has decreased during rock uplift (Pederson et al., 2002a; Roy et al., 2009). Instead, our interpretation is that denudation may be a response to post–10 Ma tectonic surface uplift (of the river and the land surface). By this model, differential isostatic rebound is a response to tectonically driven differential surface uplift such that 0.5–1.0 km of post–10 Ma “residual incision” needs to be explained by mechanisms other than erosional isostasy (Fig. 6B). While nonstationary elevation of the river may account for a portion of the “residual incision,” it is also consistent with the possibility that 25%–50% of modern-day elevation of the southwestern Colorado Plateau and Colorado Rockies may be attributed to tectonically driven post–10 Ma surface uplift.

Petrologic evidence from Neogene basaltic magmatism suggests that asthenospheric upwelling is ongoing in some regions. In western Grand Canyon, basaltic volcanism and extensional faulting have migrated eastward onto the Colorado Plateau at rates of several cm/yr over the past 10 Ma (Best and Brimhall, 1974; Wenrich et al., 1995; Roy et al., 2009), and Nd isotopes show increasingly positive values with decreasing age (Wenrich et al., 1995), consistent with asthenospheric replacement of lithosphere (Crow et al., 2011). The 1 ka basanitic Little Springs flow (star in Fig. 2B) has an ε_Nd value of +7.7, which is indicative of a depleted mantle (asthenospheric) source (Crow et al., 2011). This suggests that asthenosphere has progressively replaced or infiltrated lithosphere at the 50–90 km melt-generating depth interval for these alkaline basalts. In the Rocky Mountains, 13 Ma to 3 ka basalt geochemistry suggests a dominantly lithospheric mantle source region (Beard and Johnson, 1993); however, Nd data show evidence for mixing of lithospheric and asthenospheric basalt sources (Leat et al., 1989), as would be consistent with early stages of replacement of lithosphere by upwelling or infiltrating asthenospheric melts (Crow et al., 2011).

Xenowhiffs (mantle gases entrained in groundwater; Crossey et al., 2009) also provide direct evidence for mantle-to-surface interconnections across the region (Newell et al., 2005). CO₂-rich hot and cool springs in Colorado show highest ³He/⁴He ratios (up to 25% of the ³He/⁴He ratio seen at mid-ocean ridges [MOR]) in areas above the lowest velocity mantle, compatible with translithosphere advection of primordial ³He from partial melting associated with asthenospheric upwelling.

Figure 7 shows new CREST seismic images for an E-W cross section through the San Juan Mountains at latitude 37.5°N. This section demonstrates the spatial (and temporal) distribution of highest elevation, high exhumation, low density, and low-seismic-wave velocity crust, and low-velocity upper mantle in Colorado. This association motivates the discussion of viable models for the geodynamic development of the Colorado Rockies. Surface, crust, and mantle observations from the Southern Rockies include: (1) regional tens- to hundreds-of-kilometer-wavelength differential rock uplift between the Colorado Plateau and the Basin and Range of ∼800 m in the past 2–4 Ma (Karlstrom et al., 2008), and similar-magnitude differential rock uplift of the Colorado Rockies relative to the Colorado Plateau in the past 6–10 Ma, both spatially associated with mantle velocity gradients and/or low-velocity domains in the asthenosphere; (2) corresponding changes in river gradient across the Lees Ferry knickpoint and between the Green and upper Colorado Rivers that seem to coincide spatially with sharp mantle velocity gradients and differential incision rates and magnitudes; (3) heat flow and gravity data (Reiter, 2008) that suggest highly radiogenic, low-density rocks in the crust, and elevated temperatures in the mantle lithosphere; (4) anomalously low seismic wave velocities in the crust and mantle below the Colorado Rockies (Schmandt and Humphreys, 2010; Dueker et al., 2011); (5) regional Nd (Crow et al., 2011) and xenowhiff (Crossey et al., 2009) data that support asthenospheric upwelling; and (6) evidence that Neogene tectonism was superimposed on Laramide and mid-Tertiary uplift and cooling events, as seen in thermochronologic data.

Many of these observations, and especially the thin, low-velocity crust overlying low-velocity mantle, might be explained by a model involving delamination or distributed foundering of dense lithosphere and its replacement by asthenosphere, similar to models by Zandt et al. (2004) and Levander et al. (2011). The timing of such a detachment could have been Oligocene and/or Neogene. Geodynamic models show that the rate of this process depends strongly on lithosphere rheology (Molnar and Jones, 2004; Göğüş and Pysklywec, 2008; van Wijk et al., 2010). A significant pre-Neogene lithospheric modification in Colorado involved high-volume Cenozoic magmatism in the Rockies that likely produced a chemically dense lower lithosphere during both 70–45 Ma Colorado Mineral belt and 35–25 Ma San Juan volcanic field magmatism (Farmer et al., 2008). Resulting magmatic differentiation and melt flux probably made the crust more radiogenic and felsic (Reiter, 2008), as reflected in low mid-/lower-crustal surface wave velocities of 3.5 km/s shown in Figure 7. Laramide to Oligocene creation of a more felsic crust modified the compositional and thermal structure of the lithosphere to set the stage for any Neogene mantle dynamic model. We speculate that high-density cumulates and restites in the mantle lithosphere that formed during the early magmatic phase may subsequently have foundered. Such removal of mantle lithosphere material, for example, to initiate the 35–25 Ma ignimbrite flare up, would have disturbed the Moho (van Wijk et al., 2010), and potentially thinned the crust (Fig. 7) via delamination, leaving the present continuous (healed), but shallow Moho. Removal of dense lower lithosphere in the Oligocene would also have resulted in surface uplift. Modeled amplitudes of uplift found in numerical models of delaminating/detaching lithosphere range between several hundreds of meters to more than 1 km (Göğüş and Pysklywec, 2008), which may be recorded in the Oligocene cooling histories in southern Colorado (Kelley et al., 1992).

If mantle lithosphere was removed recently, around 10 Ma, detachment-related uplift could have contributed to the 6–10 Ma phase of exhumation. Removal of large portions of mantle lithosphere material would result in decompression melting and related volcanism; perhaps reflected in scattered 5–12 Ma Cenozoic plutonism in SW Colorado.

An additional, or perhaps alternative, mechanism in Figure 7 shows low-velocity mantle that may be upwelling from 400 km depths below the San Juan Mountains (Dueker et al., 2011). These low-velocity “pipes” seem to emanate from the 410 km mantle discontinuity and resemble upper-mantle diapirs such as may be stimulated by Farallon slab segments sinking across the 410 km discontinuity (Faccenna and Becker, 2010). Such perturbations can trigger upward diapiric flow of hydrous melt that is present just above the 410 km discontinuity due to the different water contents of olivine (above) versus wadsleyite (below) (Bercovici and Karato, 2003; Richard et al., 2006; Leahy and Bercovici, 2007). Upward ascent of hydrous 410 km discontinuity upper-mantle diapirs might have contributed to renewed basaltic magmatism, surface uplift, and accelerated exhumation starting ca. 10 Ma across the region. Such a two-stage model involving mid-Tertiary delaminations below the Oligocene San Juan volcanic field followed by recent 410 km discontinuity–related diapirs that trigger additional uplift and magmatism would be consistent with the thermochronology evidence for multistage denudation pulses in the Rocky Mountains.

Figure 7 shows a complex sub-Moho structure zone beneath Colorado that we refer to as the lithosphere-asthenosphere mixing zone because it has characteristics of both lithosphere and asthenosphere, and perhaps reflects a process of heterogeneous alteration of Proterozoic chemical lithosphere by infiltrating asthenosphere. This zone is demarcated between negative velocity interfaces seen in S-wave receiver function images (Hansen et al., 2011). The upper interface, at 80–120 km, is within a transition from 4.4 km/s mantle above (lithosphere-like) to <4.2 km/s mantle below (asthenosphere-like) such that much of this lithosphere-asthenosphere mixing zone likely has the temperature (T >1000 °C) and rheological characteristics of mantle asthenosphere. However, the lithosphere-asthenosphere mixing zone is highly structured with negative S-wave receiver function interfaces (dashed lines) at 150–200 km that suggest relic compositional interfaces more similar to lithosphere than asthenosphere (Dueker et al., 2001; Karlstrom et al., 2005). The lithosphere-asthenosphere mixing zone concept is also supported by xenolith localities in the region, with examples of Devonian diamond-bearing xenoliths from 150 km depth near the Stateline, Colorado, district to the north (Fig. 7A; Farmer et al., 2005), lithospheric xenoliths at ∼900 °C in the Rio Puerco Necks volcanic field of the Jemez lineament in the eastern Colorado Plateau region (Porreca et al., 2006; Fig. 7B), and xenoliths at ∼140 km from Navajo volcanic field (Fig. 7C) that record cold sub-Moho temperatures at the time of Oligocene eruption (Smith et al., 2004).

To summarize, any viable mantle geodynamic model needs to not only explain correspondence of the highest elevations, thinnest crust, low-velocity crust, and mantle, but also to consider the apparent adjustment of erosional systems to these mantle domains in the past 10 Ma. Existing data suggest that patterns in channel steepness and incision rate are highest above high- versus low-mantle velocity domains, implying an ongoing surface response to young and ongoing mantle convection. The utility of Colorado River profile lies in the fact that it was graded to sea level starting ca. 6 Ma, such that differential long-term incision rates can be used to reconstruct paleoprofiles to distinguish surface uplift (relative to sea level) from only differential rock uplift for upstream reaches (Hamblin et al., 1981; Karlstrom et al., 2007, model 1). A continuing goal is to distinguish various scales of differential rock uplift from surface uplift. This study emphasizes that realistic solutions need to address and integrate multiple spatial and temporal scales to resolve the interacting mechanisms by which the mantle may be uplifting the western U.S. region while also causing differential surface response in individual subregions and channel segments.

Spatial and temporal correlations among surface parameters (physiographic provinces, elevation, roughness, channel steepness, and incision rates and magnitudes) and mantle velocity domains can be used to explore the ways in which mantle flow and buoyancy variations may be driving differential surface uplift of the western Colorado Plateau and Rocky Mountains relative to the central Colorado Plateau since 10 Ma. In the Rockies, the spatial correspondence of highest topography, lowest-velocity mantle, low-velocity crust, and thin crust, and a 5–10 m geoid high appears to necessitate that low-density crust and upper mantle contribute significantly to the support of the high topography. Thermochronology and river incision data indicate an onset of rapid regional exhumation regionally at 6–10 Ma, at a wide range of elevations and across the region, suggesting that present mantle velocity variations at least partly reflect Neogene mantle flow. Mantle flow and differential uplift of 0.5 to 1.5 km is also manifested by differential river incision during the past several million years. This amounts to 25%–50% of the post–Late Cretaceous uplift of the Colorado Plateau and Colorado Rocky Mountain regions that has occurred in the past 10 Ma. Nd data corroborate the hypothesis that convective replacement of lithosphere by upwelling asthenosphere within some low-velocity domains has taken place since 10 Ma. The importance of climatic influences needs to continue to be tested by identifying any temporal variations in incision rates across the entire region that may correspond to a given climate change event, versus more persistent forcings over the several- to ten-million-year time scales that might be more compatible with buoyancy associated with mantle flow. A goal of continuing work is to use surface evidence such as long-term incision histories of continental-scale rivers to elucidate geodynamic mechanisms and scales at which mantle flow may be shaping isostatic and dynamic topography at great distance from present-day plate boundaries.

The Colorado Rockies Experiment and Seismic Transects (CREST) was supported primarily by the National Science Foundation Continental Dynamics Program under award EAR-0607808. Additional support was from EAR-0711546 and EAR-0838575, and from the Institute of Geophysics and Planetary Physics at Los Alamos National Laboratory. We thank the Incorporated Research Institutions for Seismology (IRIS) PASSCAL Instrument Center at New Mexico Tech for facility support and field assistance. Data collected will be available through the IRIS Data Management Center. The facilities of the IRIS Consortium are supported by the National Science Foundation (NSF) under cooperative agreement EAR-0552316, the NSF Office of Polar Programs, and the Department of Energy (DOE) National Nuclear Security Administration. The CREST team also includes Lang Farmer, Colin Shaw, Eric Leonard, Clem Chase, Alex Nereson, Rex Cole, and others. We thank Rebecca Flowers and an anonymous reviewer for their reviews, which helped improve the paper.