Lithospheric density models reveal evidence for Cenozoic uplift of the Colorado Plateau and Great Plains by lower-crustal hydration

Paleozoic and Mesozoic marine sediments blanket the Colorado Plateau and Great Plains, requiring broadly uniform elevations near sea level through much of the Phanerozoic. Since the Late Cretaceous, however, these units have risen—generally undeformed—to average elevations of 1.9 km in the Colorado Plateau and 1.6 km in the western Great Plains (Fig. 1A). By contrast, on the eastern Great Plains, they remain only 300 m above modern sea level. Modern topographic relief therefore reflects differential uplift, so understanding what supports this relief will provide a window into the processes responsible for the widespread but cryptic Cenozoic modification of intraplate North America.

Although the Sevier and Laramide orogenies caused substantial horizontal contraction in the modern Basin and Range and Southern Rockies, respectively, the Colorado Plateau and Great Plains experienced <5% shortening (Davis, 1978; Tikoff and Maxson, 2001). Instead, the major impact of Farallon subduction was a period of subsidence and shallow marine sedimentation (Sloss, 1963). Sediments increased crustal buoyancy beyond presubsidence values, and Tertiary removal of the Farallon slab allowed the crust to rebound to correspondingly higher surface elevation (Mitrovica et al., 1989). Digitizing the isopach maps of Cook and Bally (1975), we find that an average of 900 m of post-Jurassic sedimentary rock is preserved on the Colorado Plateau, and across the Great Plains, their average thickness grades from 1.3 km near the Rocky Mountain front to minor net erosion in eastern Kansas, Oklahoma, and Nebraska.

Assuming a density of 2000–2200 kg/m³ for these sedimentary rocks (e.g., Spencer, 1996) and asthenospheric density of 3200 kg/m³ (e.g., Lachenbruch and Morgan, 1990), they account for ∼300–350 m of uplift of the Colorado Plateau relative to the eastern Great Plains and 400–500 m in the western Great Plains (Fig. 1B). In the Colorado Plateau, 1.3 km of Cenozoic uplift remain unexplained, and on the Great Plains, there is a roughly linear gradient in remaining uplift from 900 m near the Rocky Mountain front to 0 m on the eastern Great Plains.

Because the Colorado Plateau and Great Plains share similar Mesozoic histories and even overlie the same Proterozoic terranes (Whitmeyer and Karlstrom, 2007), systematic differences in modern density that support modern topographic relief most logically result from Cenozoic modification/deformation. The EarthScope Transportable Array (TA) has recently offered unprecedented seismic coverage of the central United States; here, we used TA-based seismic velocity models (Shen et al., 2013) along with heat-flow, gravity, and topographic data to develop three-dimensional (3-D) lithospheric density models. These estimates map and quantify the contributors to modern elevation differences—and by inference Cenozoic surface uplift—across the Great Plains and Colorado Plateau.

Surface uplift may occur in response to changes in asthenospheric flow (i.e., dynamic topography) or in response to decreases in the density of the lithosphere, which depends on both temperature and composition. Therefore, the possible causes of uplift relative to the Paleozoic–Mesozoic can be enumerated (following McGetchin et al., 1980; Morgan and Swanberg, 1985): crustal thickening, crustal heating, crustal phase changes, advective or conductive heating of the mantle, chemical or phase changes in the mantle lithosphere, or dynamic topography (here defined as a non-isostatic contribution from mantle flow, not simply mantle buoyancy).

Indeed, most of these mechanisms have been proposed for the Colorado Plateau and/or Great Plains. Nevertheless, crustal thickening (Bird, 1984; McQuarrie and Chase, 2000) seems unlikely because of the minimal shortening observed at the surface (Davis, 1978; Tikoff and Maxson, 2001). Moreover, crustal thickness (gray line in Fig. 1B) does not vary systematically across the Great Plains and is not correlated with surface elevation (r² = 0.19). While crustal heating decreases density and increases elevation, modern heat flow is a relatively uniform 50–60 mW/m² across the Great Plains (e.g., Blackwell and Richards, 2004). Lower-crustal and upper-mantle seismicity in the Colorado Plateau interior (Wong and Humphrey, 1989) also suggests low thermal gradients, and joint analysis of heat-flow data and P_n velocities estimates a Moho temperature of <700 °C (Schutt et al., 2018). Around the margins of the Colorado Plateau, heat flow is modestly higher, up to ∼75 mW/m² (e.g., Blackwell and Richards, 2004), and Moho temperatures average 900 °C (Schutt et al., 2018), meaning that the average temperature of the Colorado Plateau Moho is ∼800 °C. By contrast, temperatures of ∼550 °C typify the Great Plains (Schutt et al., 2018). Assuming a roughly linear geotherm from the Moho to the surface, our crude calculation suggests an average difference of 125 °C. For a coefficient of thermal expansion of 3.0 × 10⁻⁵/°C and a reference density of 2800 kg/m³, this difference accounts for a 10.5 kg/m³ density difference. By Equation 1, the 40-km-thick crust of the Colorado Plateau thus supports only 130 m of additional elevation.

Each of the remaining hypotheses—dynamic topography, mantle or crustal composition/phase changes, and mantle heating—makes testable predictions about modern lithospheric seismic velocity and/or density structure, but these predictions have yet to be investigated systematically.

To discriminate among these remaining possibilities, we derived 3-D density models of the crust and upper mantle following the approach of Levandowski et al. (2015), which jointly analyzed shear-wave velocity (Shen et al., 2013), gravity (Fig. 2B), topography (Fig. 1A), and heat flow. Velocity is scaled to density to create a 3-D starting density model (Fig. 3). Gravity and flexurally smoothed topography are then forward modeled (Figs. 2C–2D) and compared to observations. Finally, a random-walk Monte Carlo algorithm iteratively refines the density model until free-air gravity and flexurally smoothed topography are reproduced at all points in the study area to within 10 mGal and 100 m, respectively (Figs. 2G–2H). (L1/L2 [L1 is the mean absolute value; L2 is the root mean squared] norms are generally ∼2.5/4 mGal and 30/50 m.)

Shen et al. (2013) provided hundreds of velocity models at each of ∼1000 TA stations. These models comprise one-dimensional (1-D) velocity profiles and include a sedimentary thickness and crustal thickness, both of which also vary among the individual 1-D models at any given station. They were derived by joint Bayesian inversion of short-period Rayleigh wave dispersion from ambient noise, long-period dispersion from ballistic Rayleigh waves (a total range of ∼8–80 s), and receiver functions. Nevertheless, there is not a unique velocity model that is capable of reproducing these data (which themselves are subject to uncertainty). In an effort to minimize the bias imposed on our density models of the vagaries of a given S-velocity model, we created 1000 simulations of the 3-D density structure, and each simulation began by randomly selecting one of the acceptable Vs models at each station and scaling these to density.

This polynomial regression—or any reasonable regression—cannot faithfully reproduce the density of all lithologies because of the natural spread in density of rocks that have similar velocity; the range in density about a given seismic velocity is approximately ±150 kg/m³. In particular, this regression systematically underestimates the density of mafic rocks and overestimates the density of felsic units (Levandowski et al., 2015).

We explored the effect of different choices of υ₀ in the Supplemental Material¹, but this choice can be viewed as a hypothesis to be tested. If, for example, the predicted densities of regions such as the Basin and Range or Snake River Plain with lower velocities in the mantle are too high—as manifest in topography and/or gravity—we would conclude that a lower υ₀ would be indicated.

Equation 3 was derived for a composition of 30% Fo₉₀, 30% Fo₉₂, 25% orthopyroxene (opx), 10% clinopyroxene (cps), 2.5% garnet (Gt), 2.5% spinel (Sp), and 1 mm grains, where Fo indicates the proportion of forsterite in olivine. Because the anelasticity data (Jackson and Faul, 2010) are for single-crystal olivine only, we assumed that anelastic effects were similar for other minerals. The effects of grain size and solidus velocity are addressed in the Supplemental Material (see ^{footnote 1}), but we found that our results were robust with respect to changes in these two factors. For each mineral, we accounted for the pressure and temperature dependence of the shear and bulk moduli (Bouhifd et al., 1996; Hugh-Jones, 1997; Jackson et al., 2003; Afonso et al., 2005) and for the temperature dependence of thermal expansivity (Afonso et al., 2005). To account for anelasticity, we also accounted for the temperature, pressure, and seismic period dependence of the dynamic compliance, or the Laplace transform of the creep function (using scripts provided by Jackson and Faul, 2010). Finally, the representative period of surface waves increases with depth (from ∼30 s at 50 km to 80 s—the longest period used by Shen et al. [2013]—at 150 km), so these velocity-density relationships are depth dependent, both because of increasing pressure and increasing seismic period. In other words, the thermal velocity–density relation is depth dependent because of pressure- and period-dependent anelasticity.

The lithosphere was divided into 16 layers: surface to sea level, 12 layers of 5 km thickness from sea level to 60 km, and 30-km-thick layers from 60 to 150 km depth. The region was then divided into 20 × 20 km columns. Each cell (20 × 20 × 5 km above 60 km depth or 20 × 20 × 30 km below) was assigned a uniform density, interpolated from the initial estimate from the randomly selected seismic velocity model of Shen et al. (2013) and Equations 2–3.

We then forward modeled gravity and flexurally smoothed topography (using the same procedure detailed by Levandowski et al., 2015) and compared the predictions of this initial model to observed gravity and flexurally smoothed surface elevation. To account for uncertainty in elastic thickness, each of the 1000 simulations randomly chose one of three two-dimensional elastic thickness models (Kirby and Swain, 2009; Lowry and Pérez-Gussinyé, 2011; Watts, 2012). We note, however, that the primary features that we discuss are suitably long-wavelength (>100 km) that the specific elastic thickness model has limited impact.

The velocity-density scaling reproduces the broad patterns of gravity and topography across the western United States (L1/L2 norms of 12.2/22.3 mGal and 124/155 m). Nevertheless, some short-wavelength gravity anomalies are unexplained, and there are modest misfits to topography (Figs. 2E–2F). Short-wavelength residuals may simply reflect features below the ∼100 km resolution of the TA-derived models; considering gravity as well can sharpen images of lithospheric structure (e.g., Maceira and Ammon, 2009). Broader residuals plausibly reveal compositional variations in the mantle or crustal lithologies that do not conform to our velocity-density relation. For example, mantle melt depletion lowers density substantially and slightly increases S velocity (Lee, 2003; Schutt and Lesher, 2010): Assuming uniform composition, we would underestimate the elevation in areas with depleted mantle. Similarly, hydration (e.g., serpentinization) causes a greater loss in density than estimated from our velocity-density scaling (cf. Eq. 3 with figure 1 of Christensen, 2004). Finally, the crustal velocity-density regression (Eq. 2) systematically underestimates the density of mafic units and overestimates that of many felsic rocks, with misfits as large as ∼170 kg/m³ possible for some lithologies (Levandowski et al., 2015).

In order to reproduce gravity and topography, as well as to improve lateral resolution of lithospheric structure, we next allowed departures from the initial density estimate. These departures essentially relaxed the assumptions of homogeneous mantle composition and of a fixed crustal velocity-to-density scaling, and allowed imaging of shorter-wavelength features.

Density cannot be known from seismic velocity alone, for at least three reasons. (1) Our models seek finer-scale resolution than the ∼100 km horizontal resolution of the velocity models derived from TA surface wave data. (2) There is uncertainty in the velocity models, which can be quantified in terms of the range of velocity at any given depth across the hundreds of acceptable velocity models beneath any seismic station. (3) There simply is not a single-valued mapping of velocity to density that captures all lithologies or all chemical/compositional trends (e.g., melt depletion). Additional factors, such as Vp/Vs variations, departures from the Q model chosen by Shen et al. (2013), and anisotropy would also influence the conversion of Rayleigh wave phase velocities and receiver functions to velocity profiles, and the subsequent conversion of these models to density. Put differently, there are factors that affect seismic velocity but not density, factors that affect density but not velocity, and features of the density of the crust and upper mantle—with or without a seismic signature—that are simply finer than the velocity models can resolve. Patterns and biases that are of particular importance to the present study include melt depletion, serpentinization, and the bias of Equation 2 to systematically overestimate the density of felsic material and underestimate the density of mafic material, as discussed already.

The misfits between predictions of the initial density model and observed gravity and topography are generally small compared to the ∼3 km of relief across the study area and the large variations in free-air gravity (Fig. 2). Nevertheless, to produce more robust density models, we employed the random walk Monte Carlo algorithm of Levandowski et al. (2015) to refine the initial density structures until the density model reproduced gravity and flexural topography to within 10 mGal and 100 m at all points in the study area. The Monte Carlo proceeds by selecting one of the nodes from the 20 × 20 km grid at random. A cell beneath that node and a density perturbation (limited to ±150 kg/m³ in the crust and ±50 kg/m³ in the mantle) are chosen at random, and the attendant variance reductions (or increases) of gravity and topography residuals are calculated. The algorithm is offered a number of cell/density-perturbation choices, initially two, but increasing in number through the inversion to aid convergence, and it ultimately selects the best-fitting one (even if that increases residual variance) and applies that change to the density of that cell, and updates residuals. When the model reproduces gravity and flexurally modulated topography to within 10 mGal and 100 m at all points in the study area, the 3-D density model is accepted as a member of a posterior distribution of acceptable density structures. This process—randomly selecting one of Shen’s models at each TA station, converting velocity to density, selecting one of three elastic thickness models, forward modeling gravity and flexural topography, and iteratively refining the density model—is repeated 1000 times to embrace the non-uniqueness of gravity and topography. Uncertainties of density models are discussed in the Supplemental Material (see ^{footnote 1}), but a typical uncertainty of density in a region is ∼15 kg/m³.

The magnitude of the adjustments that are necessary to match gravity and topography (Fig. 4) average ∼14 kg/m³ in the crust and 5 kg/m³ in the mantle, and adjustments are typically less than ∼100 km in the lateral dimension. The crustal adjustments are small compared to the plausible range of ∼170 kg/m³ that a rock of known velocity may have. Moreover, these adjustments are well within the uncertainties of the velocity models. Because Shen et al. (2013) provided as many as thousands (generally hundreds) of velocity models at each TA station, the velocity-derived uncertainty in density at any point is readily quantified. The range on either side of the median that subsumes 95% of their models is ∼0.2 km/s in the midcrust and upper mantle, ∼0.3 km/s in the lower crust, and greater still in the uppermost crust. As such, the inherent range of density estimated from those velocity models exceeds ±100 kg/m³ throughout the crust and is approximately ±30 kg/m³ in the mantle.

The two areas with large-magnitude and laterally extensive adjustments are an ∼300 × 300 km portion of the upper mantle beneath the Wyoming craton that is some 25 kg/m³ less dense than estimated (Figs. 4E–4F) and an ∼100 × 200 km region near the Southern Oklahoma aulacogen in which the crust is 75–100 kg/m³ denser than estimated (Figs. 4E–4F). These departures from the seismically derived density estimates are within the typical uncertainties associated with the velocity models themselves, but they could also represent systematic deviations from Equations 2–3 such as the compositional anomalies discussed earlier. In the Wyoming craton, we speculate that the low-density but high-velocity Archean mantle is melt-depleted residuum from the initial, high-temperature extraction of crustal material >2 Ga. In southern Oklahoma, three lines of evidence suggest that very dense crust is likely partially eclogitized mafic intrusions. First, the anomaly lies in or near the Oklahoma aulacogen, which underwent extension and associated emplacement of mafic sills and dikes in the Eocambrian. The aulacogen was subsequently the focus of early-stage NE-SW Ouachita contraction, which would have thickened the suite of mafic units under horizontal compression. Second, it is the densest material in the study area but is above the receiver function–defined Moho. Third, it is denser than expected from Equation 2, which—as argued by Levandowski et al. (2015) for the northern portion of the Midcontinent Rift (Fig. 4A)—is a hallmark of mafic crustal lithologies.

Earlier herein, we enumerated the possible causes of Cenozoic uplift of the Colorado Plateau and western Great Plains relative to the eastern Great Plains. We now discuss the testable predictions that each hypothesis makes about seismic velocity and density structure. The simplest prediction that each hypothesis makes is where in the lithosphere modern topographic relief is supported. For example, arguments in favor of mantle heating would require that support for modern relief is derived from mantle depths. In addition to the depths of density variations, each hypothesis makes certain predictions about the relationships between velocity and density (i.e., the validity of Equations 2–3 and of our underlying assumptions, manifest as whether adjustments to our initial, seismically derived models are necessary).

Increases in average mantle temperature could either be conductive (e.g., Roy et al., 2009) or advective, such as lithospheric thinning and replacement by warmer asthenosphere. One version of the former hypothesis—conductive reheating of unthinned lithosphere following a period of insulation by the Farallon slab—can be dismissed because such reheating does not account for any uplift relative to the Paleozoic–Mesozoic (i.e., before cooling). If Cenozoic lithospheric thinning caused differential uplift, however, the density of the lower lithosphere should decrease from east to west across the Great Plains and/or be less beneath the Colorado Plateau than the eastern Great Plains. Additionally, since our initial assumption is that mantle velocities reflect temperature variations, this systematic trend in mantle density should exist in our initial, seismically derived model.

In a region of substantially melt-depleted lithosphere, a velocity-temperature-density scaling should also lead to underpredicted elevations. Although xenoliths from the Colorado Plateau record upper mantle that is ∼1% enriched in magnesium (Alibert, 1994; Lee et al., 2001), the predicted elevation in the Colorado Plateau and western Great Plains is generally well explained by our initial density models (Fig. 2), which consider only thermal variations. Moreover, melt depletion would have to have occurred during the Cenozoic, which is at odds with the paucity of volcanism in the Colorado Plateau and Great Plains.

Mantle hydration (e.g., serpentinization) is also recorded in some Colorado Plateau xenoliths (Usui et al., 2003; Smith and Griffin, 2005), but the effects of fluid flux on velocity and density depend on the specifics of its chemical and/or compositional impact. As noted above, serpentinization causes a proportionally greater decrease in density for a unit velocity decrease than predicted by Equation 3 (compare figure 1 of Christensen, 2004, with our Equations 3a and 3b); our initial density estimate would therefore be too great if substantial serpentinized material is present in the upper mantle. By contrast, incorporation of hydroxyl groups into nominally anhydrous olivine can have a strong impact on velocity, but unit cell volumes increase only slightly (Smyth and Jacobsen, 2006); in this case, our density estimate would be lower than the true density. Finally, fluid flux can deliver silica to the mantle lithosphere and possibly increase orthopyroxene at the expense of olivine, but this compositional trend has little effect on shear velocity or density (Schutt and Lesher, 2010). We argue against the first two possibilities because the mantle density beneath the Colorado Plateau and Great Plains accords so well with an estimate based on the assumption that velocity reflects temperature alone. Our models would be insensitive to the latter trend, but orthopyroxene enrichment is not correlated with density decrease anyway (Schutt and Lesher, 2010), so it would not explain uplift.

Although we argue against the three hydration-induced changes discussed here, many other effects of fluid flux are possible. It is crucial to note that any phenomenon that affects both velocity and density, and that does so in similar proportion to that derived in Equation 3 (i.e., with ∼6 kg/m³ density change corresponding to 1% shear velocity change) would be incorrectly ascribed to temperature variations. In other words, our null hypothesis is that all velocity variations reflect temperature variations. We test this hypothesis by scaling velocity to density using a thermal relationship, and then comparing predicted and observed gravity and topography. Because observations are reproduced rather well, we do not question the null hypothesis further, but we do acknowledge that if another factor controls velocity and density and produces a similar relationship to that described by Equation 3, our simple test would cause us to wrongly accept the hypothesis that mantle velocity and density variations are primarily functions of temperature variations.

Midcrustal xenoliths from the Navajo volcanic field on the Colorado Plateau record hydration-induced retrogression of garnet-bearing crustal assemblages to less dense mineralogies. Th-Pb dating of secondary monazite associated with these assemblages suggests that the majority of the retrogression occurred in the latest Cretaceous (Butcher, 2013; Butcher et al., 2017), contemporaneous with the arrival of the Farallon slab. A similar but undated retrograde reaction is documented in crustal xenoliths near the northern Great Plains: Southward from near the Canadian border to southern Wyoming, three xenolith localities increase in elevation from 1 to 2.5 km as xenolith density decreases by ∼500 kg/m³ (Barnhart et al., 2012; Farmer et al., 2012; Mahan et al., 2012). Noting a collocated southward decrease in seismic velocity in an ∼10-km-thick layer of lower crust, Jones et al. (2015) hypothesized (following Eq. 1) that hydration-induced retrogression produced much of the modern relief on the Great Plains. If so, then it follows that the modern density of the lower crust should decrease systematically from east to west. Additionally, the xenoliths discussed by Jones et al. (2015) and Butcher et al. (2017) have P-velocity–density trends that are broadly concordant with that of the lithologies used to develop Equation 2. Therefore, it is possible that such an east-west gradient in density would be manifest in velocities as well, though shear-velocity measurements are not available for those xenoliths. We discuss lower-crustal hydration specifically, but we note that any other mechanism of broad, dedensifying phase changes might produce similar patterns.

If surface elevation in a region were sustained by asthenospheric flow rather than lithospheric buoyancy (Moucha et al., 2008; Liu and Gurnis, 2010), seismic velocity would generally reflect the density of the crust and upper mantle: The predicted gravity would match observations, but such a region would stand at greater elevation than predicted. This pattern is not revealed by the residuals shown in Figure 2 in either the Colorado Plateau or the western Great Plains, where the initial model generally reproduces modern elevation. Therefore, we suggest that the impact of asthenospheric flow is minimal, perhaps <100 m. There is an important distinction between flow and density, however. We did not explicitly separate lithosphere and asthenosphere, and we included sublithospheric loads into our flexural-isostatic buoyant heights. We did indeed find modest variations in density at depths that likely correspond to asthenosphere, but we refer to any potential dynamic component as non-isostatic forces related to flow itself. Recently, Afonso et al. (2016) reached a similar conclusion—that dynamic topography is of secondary importance—even though their modeling explicitly calculated the vertical normal stress imparted on the base of the lithosphere by both the buoyancy/antibuoyancy of sublithospheric loads and by buoyancy-induced flow.

After accounting for the east-to-west increase in the thickness of Cretaceous sediments, a nearly linear increase in elevation across the Great Plains—more than 800 m at the Rocky Mountain front—remains unexplained. Our modeling reveals that the density of the lower crust—averaged longitudinally—decreases systematically from east to west. Little systematic difference is seen at other depths, casting doubt on the role of midcrustal alteration or mantle processes in supporting modern topographic slope. The average lower-crustal density (20–45 km) is 112 kg/m³ less in the western Great Plains than in the eastern Great Plains (Fig. 5), supporting ∼900 m of modern surface relief (Fig. 1B). Thus, density variations in the lower crust combine with Cenozoic sedimentation to account almost exactly for Cenozoic differential uplift across the Great Plains. Of the potential causes of uplift, crustal hydration is most consistent with this pattern, speculatively reflecting progressive dewatering of the Farallon slab with distance eastward under North America.

Finding, not surprisingly, a similar pattern in lower-crustal density, Levandowski et al. (2017) suggested that the lowest densities may coincide with Proterozoic continental sutures (the Cheyenne belt—the suture between the Wyoming craton and Yavapai terranes—in SE Wyoming and the Yavapai-Mazatzal suture in SE Colorado). If so, it is reasonable to expect that fluids may preferentially exploit these preexisting, lithospheric-scale fracture/suture zones, and lower-crustal hydration may consequently be greatest there. In fact, studies of hydrated mantle xenoliths in the Colorado Plateau indicate that fluid flow is controlled by the presence of fractures (Nielson et al., 1993). Although this line of reasoning presents an internally consistent explanation for the heterogeneous lower crust of the western Great Plains, there is little agreement on the location and nature of Proterozoic sutures on the Great Plains. The Cheyenne belt is the topic of comparatively less debate, but the Yavapai-Mazatzal suture as sketched by Levandowski et al. (2017)—based on Carlson (2007)—is near the northern edge of the suture zone depicted by Whitmeyer and Karlstrom (2007) and is substantially north of the location given by Magnani et al. (2004). Thus, the details of lower-crustal density structure provide—at best—speculative additional evidence in favor of lower-crustal hydration, specifically via the possibility that preexisting suture zones served as comparatively more efficient conduits for mantle-derived fluids than did intervening terranes.

Our purpose in this work was to investigate the broad-scale uplift of the Great Plains and Colorado Plateau. We argue that sedimentation and lower-crustal phase changes explain nearly all of the longitudinally averaged uplift. Smaller-scale topographic features may reflect other processes. For example, low-density mantle near the Jemez Lineament in northeastern New Mexico and southeastern Colorado supports ∼500 m of uplift relative to the eastern Great Plains (Fig. 3). It is compelling to note that this area also stands 200–600 m higher (up to 2200 m; Fig. 1A) than most of the rest of the western Great Plains (Fig. 1), perhaps providing evidence in favor of mantle-derived uplift (Nereson et al., 2013). We conclude that second-order topographic features may reflect other mechanisms of support, but the broad difference between the western and eastern Great Plains arises from sedimentation and lower-crustal density.

Between 20 and 40 km depth, the Colorado Plateau averages 106 kg/m³ less dense than the eastern Great Plains (Fig. 5), supporting 660 m of relief (Fig. 1B). As discussed earlier, there is likely some contribution from higher temperatures in the Colorado Plateau; following the same line of logic, we would estimate a 250 °C temperature anomaly at the Moho (∼21 kg/m³) and a ∼125 °C temperature anomaly at 20 km depth (∼10.5 kg/m³). Thus, the 20–40 km depth range may average 15 kg/m³ less dense beneath the Colorado Plateau than the eastern Great Plains because of its temperature, but the remaining 91 kg/m³ difference is better explained by compositional changes.

In addition, the 40–50 km depth range essentially subsumes the crust-mantle transition (Gilbert, 2012; Shen et al., 2013), and receiver functions generally image a gradient or low contrast in impedance across the Moho in the Colorado Plateau. The density of this transition zone is 55 kg/m³ lower than at the same depths beneath the eastern Great Plains, and this explains an additional 170 m of differential Cenozoic uplift. After accounting for sedimentation and the effects of lower-crustal density loss, only 400 m of Colorado Plateau uplift remain.

The Colorado Plateau mantle is similar in density to that of the southern Rockies and Basin and Range and accommodates ∼400 m of topographic relief relative to the eastern Great Plains (Fig. 1B). We infer that this density structure is thermal in origin (i.e., the Rockies, Colorado Plateau, and Basin and Range have warmer mantle than the Great Plains) because of the close correspondence between the elevation and gravity predicted by our thermal relation between velocity and density (Eq. 3) and observation. More obvious evidence comes from the elevated Moho temperatures in the Colorado Plateau relative to the Great Plains (Afonso et al., 2016; Schutt et al., 2018). We do not explicitly distinguish between lithospheric mantle and asthenosphere, so higher temperatures could reflect thinner lithosphere or warmer lithosphere; if the lithosphere and asthenosphere are in thermal equilibrium, the two are likely intertwined, because convectively thinned lithosphere will subsequently warm. A Cenozoic increase in mantle temperature is therefore a possible source of Colorado Plateau uplift, but we have already discounted purely conductive heating, so we now focus on advective heating.

Advective heating may occur by removal of mantle lithosphere and its replacement with asthenosphere or by intrusion. Limited volcanic activity and intrusive activity, as well as crustal thickness estimates, are at odds with large-magnitude igneous intrusion, so we suggest lithospheric thinning. A sinking Rayleigh-Taylor instability (Levander et al., 2011), delamination (Bird, 1979), or ablation by the Farallon slab (Bird, 1988) could thin the lithosphere; we do not explicitly discriminate among these. Additionally, we did find some heterogeneity in density of the Colorado Plateau mantle, with comparatively higher density in the north-central portion, lower density around the margins, and the lowest densities in the south (Figs. 5E–5F). Heating/thinning of the Colorado Plateau lithosphere is likely heterogeneous, as also evidenced by encroachment of volcanism from the edges (Roy et al., 2009), but we will illustrate the effects of lithospheric thinning with 1-D thermal models meant to be reflective of the average across the Colorado Plateau. Because buoyancy changes are convolved with the flexural response of the lithosphere (and the Colorado Plateau lithosphere is comparatively strong, with elastic thickness ∼30 km: Kirby and Swain, 2009; Lowry and Pérez-Gussinyé, 2011; Watts, 2012), it is indeed the average change in mantle buoyancy that is of interest to our present aim.

If the Mesozoic Colorado Plateau did indeed resemble the modern eastern Great Plains, the lithosphere would originally have been 150–200 km thick (van der Lee and Nolet, 1997; Yuan et al., 2014). In addition, if only the lower-mantle lithosphere was removed, then remaining melt-depleted material could later be sampled in xenolith suites (as argued by Spencer, 1996).

To estimate how much lithospheric thinning would be required in order to produce 400 m of modern uplift, we solved non-steady-state heat-transfer equations similar to those presented by Bird (1979). In that conception, a portion of lithosphere is removed at time 0 and replaced with uniform-temperature asthenosphere. The replacing material then remains fixed in space and cools conductively, essentially becoming thermal lithosphere. Here, we revisited this modeling and also solved similar non-steady-state heat-flow equations for the end-member situation that after lithospheric thinning, the replacing material does not cool off but rather is maintained convectively at a constant temperature, such that the lithosphere is permanently thinned. The former setup ignores convection, but the latter requires a long-term change in heat flux; neither is a perfect formulation. We view the two conceptions as bounding the reality that lies somewhere between these two idealizations.

If the asthenosphere cools as a semi-infinite conductive half-space (essentially becoming thermal lithosphere), the temperature profile returns to a linear geotherm from the surface (T = T_s) to z_old (T = T_a). If the lithosphere thins permanently, temperatures approach a steeper geotherm from the surface (T = T_s) to z_new (T = T_a).

The Fourier expansion of Equation 5 is (Boyce and DiPrima, 2003):

Here, L is the depth of the future conductive boundary (z_new or z_old), and w(z,0) = T(z,0) – v(z). The thermal diffusivity, κ, is 1 mm²/s.

The final unknowns are the initial lithospheric thickness—we used plausible values of 150 and 200 km, based on estimates of the modern eastern Great Plains (van der Lee and Nolet, 1997; Yuan et al., 2014)—and when thinning occurred. Figures 6A and 6B show a number of uplift versus time curves. (We used a constant coefficient of thermal expansion, α, of 3.2 × 10⁻⁵/°C.) If the lithosphere thins permanently, initial uplift is followed by a protracted period of asymptotic, ongoing surface uplift. If the replacing material cools, the elevation gain will decay away. Since there is no discernible heat flow anomaly in the modern Colorado Plateau, we also checked the predicted change in surface heat flow from each model (quantified as k[T_5km – T_s]/5, where k = 3 mW m^–2 °C^–1), shown in Figures 6C and 6D. The change for all models that produced 400 m of modern uplift relative to the initial state is a few milliwatts per square meter, which is allowed by observations. We next calculated the depth of the lithosphere-asthenosphere boundary immediately after thinning and the amount of lithospheric thinning needed to produce ∼400 m of modern uplift for the various scenarios (transient vs. permanent thinning and initial lithospheric thickness of 150 vs. 200 km) for thinning at 70 Ma (the arrival of the Farallon slab), 30 Ma (slab rollback/the ignimbrite flare-up), and 0 Ma (Figs. 7A–7B).

All acceptable scenarios left most of the mantle lithosphere intact to later be sampled as magnesian xenoliths (Alibert, 1994; Lee et al., 2001). Transient thinning required 45–70 km of lithosphere to be removed; permanent thinning required 25–50 km (Fig. 7B). Thus, we suggest that the Colorado Plateau and the western Great Plains experienced similar amounts of lower-crustal hydration by fluids exsolved from the Farallon slab, but the lithosphere beneath the Colorado Plateau also thinned by a few tens of kilometers during the Cenozoic, causing an additional 400 m of uplift.

Because the removal of mantle lithosphere by Rayleigh-Taylor instability, delamination, or similar mechanisms should create a number of small-scale convection cells, a more appropriate treatment of the problem might hold the lithosphere-asthenosphere boundary at a constant temperature for some time before allowing the replacing asthenosphere to begin to cool. The differential equations required to solve such a problem are the same as Equations 4–7. We conducted many trials with varying durations of what is effectively forced convection, and even a comparatively short duration of the constant-temperature-asthenosphere phase greatly protracted uplift and decreased the amount of thinning necessary. As such, the true amount of lithospheric thinning and the thickness of remaining lithosphere are likely closer to the permanent thinning case, <50 km.

An additional set of calculations was performed (not shown) to determine the effect of a period of refrigeration of the overlying lithosphere. Specifically (e.g., as posited by Roy et al., 2009), a shallowly subducting slab would isolate the base of the lithosphere from asthenospheric convection and cause the lithosphere to cool, densify, and subside. We find that even 20 m.y. of refrigeration (150-km-deep lithosphere-asthenosphere boundary cooled from 1350 °C to 900 °C and held fixed for 20 m.y.) would only produce ∼250 m of thermal subsidence. Since we also posit that the cooled lower portion of the lithosphere was removed, if thinning occurred after or as the slab was removed, the effect on surface topography would be inconsequential.

Cenozoic uplift of the largely undeformed Colorado Plateau and Great Plains was mainly achieved by hydration-induced retrogression of the lower crust by Farallon slab–derived fluids and accumulation of sediments during the Cretaceous and Early Tertiary. The 1.3 km of differential Cenozoic uplift across the Great Plains is due in part to an east-west gradient in the thickness of post-Jurassic sediments, which supports 0.4 km. A separate decrease in lower-crustal density, which can be explained by hydration due to mantle-derived fluids, caused the remaining 900 m. Of the 1.6 km of relative uplift of the Colorado Plateau, 250 and 950 m can be ascribed to sedimentation and lower-crustal hydration, respectively, and the final 400 m most plausibly resulted from Cenozoic removal of the lower 25–70 km of the mantle lithosphere.

We thank Derek Schutt, Juan Carlos Afonso, and an anonymous reviewer for their comments. Mahan received funding from National Science Foundation grants EAR-0746246 and EAR-1053291. The development of the density modeling algorithm was funded by a Mendenhall Postdoctoral Fellowship to Levandowski.