Reconstruction modeling of crustal thickness and paleotopography of western North America since 36 Ma

Crustal thickness and topography are a direct result of interactions among mantle convection, continental dynamics, and climatic and erosional processes. Hence, the topographic evolution of mountain belts and continental interiors reflects directly upon the coupling between mantle and surface processes. The Basin and Range province of the western United States, located between the Sierra Nevada and the Colorado Plateau, is a unique modern intracontinental extensional province (Dickinson, 2002, 2006) (Fig. 1). This extension follows a protracted phase (155–60 Ma) of crustal shortening and mountain building associated with a history of subduction and terrane accretion (Hamilton, 1969; Saleeby, 1983; Burchfiel et al., 1992; Saleeby et al., 2003; Liu et al., 2008; Sigloch and Mihalynuk 2013). Since ca. 40 Ma, the western North American plate tectonic history has involved a complex transition from early shallow- to flat-subduction of the east-dipping Farallon slab along the western margin to its present transtensional environment associated with the evolution of the San Andreas fault system (Bird, 1988, 1998; Liu, 2001; Dickinson, 2002, 2003, 2006; Humphreys and Coblentz, 2007; Humphreys, 2009). Through this evolution of changing boundary forces, the lithosphere underwent a profound transformation (Atwater, 1970; Coney and Harms, 1984), leading up to the present-day setting involving both shear and extension. Over time, as the boundary conditions evolved, topography, along with crustal thicknesses, were dramatically altered from probable high elevations of orogenic plateaus, with corresponding thick crustal welts (Coney and Harms, 1984; Wernicke, 2011), to the current Basin and Range–style topography and thin crustal structure (Wernicke et al., 1987; Atwater and Stock, 1998; McQuarrie and Wernicke, 2005; Roy et al., 2009).

While the Basin and Range landscape has increased in area by a factor of two over the past ∼20 m.y. (Hamilton and Myers, 1966; Wernicke et al., 1988; Burchfiel et al., 1992), the southern part of the Great Basin has undergone 200% west-east extension and the northern portion has undergone 50% since the middle Miocene (Wernicke et al., 1988; McQuarrie and Wernicke, 2005). The average present-day continental crust thickness within the Basin and Range area of the western U.S. is ∼30 km (Shen and Ritzwoller, 2016). Although it is recognized that the widespread extension across the Basin and Range impacted the regional climatic, faunal, floral, and mammal evolution across the North American southwest (Spencer et al., 2008; Badgley, 2010; Badgley et al., 2014), considerable controversy surrounds the mechanisms responsible for such extension (Sonder and Jones, 1999; Flesch et al., 2000, 2007; Humphreys and Coblentz, 2007; Ghosh and Holt, 2012; Ghosh et al., 2013a, 2013b).

The eastern and southern margins of the Great Basin contain belts of metamorphic core complexes. These unique features are thought to be a result of the extreme thinning of crust within the Basin and Range province (Davis, 1980). In the crustal thinning process, deep crustal levels are pulled from underneath an area of extensional detachment and exposed at the surface (Wernicke and Burchfiel, 1982; Davis, 1988; Dickinson, 2004). Coney and Harms (1984) argued that the extreme extension occurred above zones of crustal welts, which resulted from pre-extensional crustal horizontal shortening.

We reconstruct the crustal thickness and paleotopographic evolution of the Great Basin and southern Basin and Range regions using time-dependent crustal strain rate estimates with formal uncertainties obtained from position estimates of geologic reconstruction by McQuarrie and Wernicke (2005). Quantitative models of crustal thickness evolution between 36 Ma and the present are achieved by integrating the strain rates back through time. We then adopt a compensation model in order to calculate a paleoelevation. The compensation model involves a variable upper mantle density with a depth of compensation of 100 km below sea level for the present-day western U.S. To produce the quantitative models of surface elevation evolution from 36 Ma to present, we track the positions of this upper mantle density field through time in the same way that we track the positions of crustal thickness. We discuss the methodology and the propagation of errors for the strain rates, crustal thicknesses, and surface elevations. We also discuss the methodology for testing the effect of thermal perturbations as a result of slab rollback and the migration of volcanism on our upper mantle density and surface elevation models beneath the Basin and Range area.

The Basin and Range province is an ideal locality to test the validity of our assumptions for calculating topography and lithosphere thickness through time due to: (1) the widespread and detailed geologic mapping that has been completed, resulting in a fairly complete palinspastic reconstruction of Tertiary extension (Coney and Harms, 1984; McQuarrie and Wernicke, 2005), as well as some tests of paleoelevation from paleothermometers like stable isotopes, and clumped isotopes in lacustrine carbonates (e.g., Huntington et al., 2010; Lechler et al., 2013); (2) Earthscope USArray seismic network experiments that provide valuable new constraints on crustal and mantle structures (Schmandt et al., 2015; Porter et al., 2016; Shen and Ritzwoller, 2016); and (3) studies directed at lithosphere evolution through the Cenozoic using thermochronology, geochronology, and geochemistry (e.g., Flowers et al., 2008; Bidgoli et al., 2015).

Our crustal thickness and topography models for the western U.S. will provide constraints for a high-standing topography prior to Basin and Range extension. The estimates of paleotopography and crustal thickness evolution, with formal uncertainties through time, will provide constraints for quantifying the magnitude and distribution of lithospheric body force distributions and their influence on collapse and development of the Basin and Range within the Pacific–North American plate boundary zone. Our topography model also may help by providing constraints for studies that explore linkages between topographic evolution, landscape evolution, and floral and faunal changes. Furthermore, our time-dependent crustal thickness model provides a regional context for understanding the spatial distribution of exhumation histories, reflected in the vast and growing number of thermochronological studies.

McQuarrie and Wernicke (2005) provided position estimate changes with uncertainties based on a comprehensive compilation of structural information, synthesis of east-west profiles, and incorporation of constraints from Pacific plate motion history. For a number of points, they provided the direction and magnitude of displacements, which occurred between specified time intervals. Although this deformation is accommodated by a large number of fault-bounded blocks (resulting in the Basin and Range province), we model the field through time as a continuum. Our continuum approach is justified due to a number of practical limitations in our knowledge of the true field of finite fault slip history. For example, whereas the slip rate of an individual fault may, at a specific point in time, possess high uncertainties, the integrated offset across several faults over a finite time interval will be more precisely known. Thus, the estimate of the average strain rate within a finite volume is going to be more precise than a point estimate of strain rate at a specific fault location. We argue that the strain rates within finite volumes of the field provided by McQuarrie and Wernicke (2005) can be reliably quantified, and a continuum treatment is the simplest approach for modeling the strain history in the western U.S.

To reconstruct the paleotopography of the western U.S. since 36 Ma, it is necessary to obtain estimates of the kinematics of the lithosphere through time. To this end, we use position estimates from the fault displacement compilation of McQuarrie and Wernicke (2005) for 857 sites in the western U.S. to calculate displacement rates (Fig. 2). McQuarrie and Wernicke (2005) provided values in 2 m.y. bins for periods between the present time to 18 Ma, and in 6 m.y. bins for periods between 18 and 36 Ma. Using these position estimates, we obtain crustal velocities, with errors, at 857 sites at 1, 3, 5, 7, 9, 11, 13, 15, 17, 21, 27, and 33 Ma. We take the midpoint between the latitude and longitude of beginning and end points of each site for a specified time interval, along with the total magnitude and direction of displacement, to obtain a velocity vector for that location. This vector represents the average velocity of the deforming lithosphere for that location during the specified interval of time. We define a grid where we can use all velocity estimates to determine a continuous model of a velocity gradient tensor field for the time periods above (Fig. 2). Because locations of the stable blocks like the Colorado Plateau and Great Plains are the firmest constraints for delineating important boundaries between relatively undeformed and deformed lithosphere during the time interval of interest (0–36 Ma) (McQuarrie and Wernicke, 2005), we put a grid of points (0.5° × 0.5°) with zero velocities and 1 mm/yr standard errors as additional constraints for the strain rate model within these stable blocks (Fig. 2). During the time period that the Rio Grande rift was opening, based on McQuarrie and Wernicke’s (2005) model, we add the constraint of velocities on the Colorado Plateau that favor little to no internal deformation and that are consistent with the direction of opening of the rift. We also add a dense set of points (1.0° × 1.0°) with 2 mm/yr standard errors on the Pacific plate as constraints based on the direction of motion of the Pacific plate relative to stable North America between 16 Ma and the present (Atwater and Stock, 1998; DeMets and Dixon, 1999; McQuarrie and Wernicke, 2005) (Fig. 2). Having constraints for the Pacific plate through time provides the added benefit of resulting in a correct prediction of the timing of the coastal fragment of California being transferred to the Pacific plate. Additionally, these constraints help inform strain rate and crustal thickness models during the time interval that Baja California is transferred to the Pacific plate and the Gulf of California opens up.

The continuous approach we adopt enables us to quantify formal errors in strain rates and rotation rates averaged within the areas in Figure 2, given a set of velocity values with uncertainties. The smoothed velocity gradient tensor field solution is obtained using a damped least-squares inversion method (Beavan and Haines, 2001; Holt et al., 2000). In this method, the horizontal velocity field within a specified frame of reference is expressed as u(r) = W(r) × r, where W(r) is a continuous rotation vector function, and r is the radial position vector, passing through the Earth’s center and defining a point on the Earth’s surface. In the least-squares inversion, the values of W(r) are determined at the knotpoints of the grid, the grid line intersections (2475 points on the grid in Fig. 2). A continuous field of velocities throughout the domain is provided by bicubic spline interpolation of W(r) parameters (Beavan and Haines, 2001). Spatial derivatives of these continuous functions, W(r), provide strain and rotation rates (Haines and Holt, 1993). The parameterization thus gives us a continuous model field of horizontal velocities, strain rates ⁠, and rotation rates about the vertical axis (Beavan and Haines, 2001; Holt et al., 2000) throughout the domain for each given time step on the surface of the sphere (Fig. 3) (see Supplemental Animation S1¹).

Here, ⁠) is the variance for zero-valued observed strain rates; S_i is the area of the i^th cell on the grid. Thus, 1/υ is an adjustable strain rate variance parameter controlling the degree to which the continuous parameters can fit the input velocities (point values) or the zero strain rate constraints (area averages) (Beavan and Haines, 2001). The value of υ is a single adjustable parameter, assigned as a constant for all areas.

Large variances in strain rates impose small weighting on the second part of the functional in Equation 1. Small variances in strain rates impose larger weighting on the second part of Equation 1. Minimization of Equation 1 provides the minimum weighted sum of squares misfit to the velocities, under the constraint that there is a global minimum for the second invariant of the model strain rates, weighted by υ. Minimization of the functional of Equation 1 with small weighting factor υ implies that large model strain rates could be produced in order to minimize the weighted sum of squares misfit to the velocities.

We explore three different levels of smoothing for solutions within all time bins, with SEUW equal to 1.0 (underdamped), 1.3 (intermediate), and 1.7 (overdamped) (Fig. 4). If the uncertainties in velocities account for all unknowns in position and timing for accumulated fault motions over a given time interval, then what we call the underdamped solution would be optimal. Our solution of choice (SEUW = 1.3) provides results that are slightly more smoothed than the one where SEUW = 1.0, and accounts for the possibility that average uncertainties for all velocities are ∼30% higher than the estimates used.

The formal uncertainties in model strain rates are controlled by the spatial density of velocity observations, the velocity errors, and the degree of formal smoothing in the velocity field interpolation (Fig. 5). Note from Figure 5, which shows the formal standard errors for vertical strain rates (⁠⁠), the greater the degree of smoothing (with lower resolution), the smaller the formal errors (see Supplemental Animation S2 [^{footnote 1}]).

Given the solutions for the different time bins, we time-integrate the instantaneous flow field estimates to obtain coordinate changes and crustal thickness changes through time. In our analysis, we assume that the lithospheric deformation is vertically coherent, or that vertical variations in horizontal velocity within the lithosphere are small in comparison with horizontal gradients of horizontal velocity. We also approximate zero volume change, ∇ u = 0, and thus the vertical strain rates ⁠. Furthermore, we also ignore effects of erosion and igneous crustal addition. Given these assumptions, the vertical strain rates are used to calculate crustal thickness evolution throughout the southwestern U.S.

Later we will discuss the possibility of lower crustal flow within regions that have experienced extreme crustal extension and mid-crustal exhumation (metamorphic core complexes). We argue here that the assumption of vertical coherence is a reasonable first approximation even in the presence of lower crustal flow as long as (1) this lower crustal flow is not pervasive, but rather is isolated beneath zones of high topography and thick crust that underwent extreme extension, and (2) horizontal strain rates, averaged over 200 km length scales, do not differ by more than a factor of two between upper and lower crust. Not accounting for a factor of two difference between upper and lower crust, for extensional strain rate magnitudes found in this paper, results in only a ∼10% error in final crustal thickness estimates. Yet such a difference in strain rate, when integrated over a geological time scale (10 m.y.), can result in substantial differential offsets, much like that observed in the vicinity of metamorphic core complexes.

The calculation of crustal thickness evolution involves (1) tracking coordinate changes through time, based on the time-dependent velocity field, and (2) tracking the crustal thickness changes of those corresponding coordinates.

The parameterization of a continuous velocity field enables us to track position changes using relatively small time steps (0.5 m.y.). We use a fixed grid and velocities in North American frame, and we assume that the strain rates do not change temporally within each time interval. There is no advection of fault source structures during each time interval (2 m.y. intervals between 0 and 18 Ma, and 6 m.y. intervals between 18 and 36 Ma), but the advection of source structures from one time interval to the next time interval is accounted for. Using the horizontal velocity gradient tensor field for each specific time interval, we integrate points back in time to obtain past position coordinates for a dense grid of points (47,616 points, 0.1° × 0.1° spacing).

We make the simplifying assumption that the strain rate field is stationary during the time interval covering the particular velocity gradient tensor field solution (j). Because strain rates are spatially averaged and smoothed, the change in coordinates (Δϕ, Δθ) for any single time increment k (Δt = −0.5 m.y.) is small in comparison to distances over which there are substantial changes in strain rate values. Thus, the inaccuracies are small relative to the exact solution, which would involve tracking strain rate values in an advecting field containing fault source terms. With a change to a new velocity gradient tensor field solution (j = j + 1), the spatial shift in the strain rate field properly accounts for the advection of source faults (McQuarrie and Wernicke, 2005) relative to the North American frame.

Crustal extension and exhumation of the middle crust, as exemplified in core complex exposures, likely involved some flow of weak lower crust (e.g., Block and Royden, 1990). Such a flow of a weak and hot lower crust, enhanced by magmatic activity (Armstrong and Ward, 1991), can be a result of tectonic denudation and isostatic adjustment to crustal stretching, topographic forces, and sedimentary loading by erosion (Wernicke, 1992). Yet today, the Moho is relatively flat below regions containing core complexes (Snow and Wernicke, 2000). A present-day flat Moho suggests that the paleo–crustal root must have flattened as middle crustal rocks were exhumed. Moreover, a weak lower crust would have likely facilitated a broadening of the crustal roots prior to and during the extension, and consequently a broadening of the topography. As discussed above, our treatment using vertical coherence cannot account for lower crustal flow. We have also argued that not accounting for lower crustal flow is unlikely to lead to errors of more than 10% for paleo–crustal thickness. Nevertheless, we approximate the influence of the lower crustal flow in the following way and present crustal thickness models that incorporate this approximation in supplementary for comparison. We approximate the effect of lower crustal flow by applying additional spatial smoothing of the model strain rates. This smoothing has the influence of spreading the strain rates out spatially, but preserving the total integral of those strain rates (total extension). We will show that in comparison with the unsmoothed strain rates, the smoothed solution produces broader crustal welts and slightly lower elevations.

Integration back to 36 Ma shows a substantial crustal welt with an average crustal thickness of ∼56.5 ± 2.5 km within eastern Nevada, eastern California, and northwestern, central, and southeastern Arizona (Fig. 7). The formal estimates of standard errors from the model strain rate fields enable us to calculate error propagation for crustal thickness H (ϕ, θ, t) at all time steps (see the Appendix) (Fig. 8) (see Supplemental Animations S3, S4 [^{footnote 1}]). The smoothed solution that approximates the effect of lower crustal flow shows a slightly wider zone of crustal welt that with an average of ∼54.3 ± 2.5 km, which, as expected, is wider and not as thick as the unsmoothed solution (see Supplemental Fig. S1², and Supplemental Animation S5 [^{footnote 1}]). We also analyzed the CRUST 1.0 model (Laske et al., 2013) for the western U.S. for comparison with results obtained using the model of Shen and Ritzwoller (2016) (see Supplemental Figs. S2, S3 [^{footnote 2}]).

To create topographic maps of the western U.S. since 36 Ma, it is necessary to have a compensation model. Becker at al. (2014) argued for significant contributions from dynamic topography in the presence of a thinned lithosphere-asthenosphere boundary below the Great Basin and southern Basin and Range. Time-dependent topographic variations owing to possible temporal variations in the convecting regime are beyond the scope of the study, and such a treatment would require input from time-dependent models (e.g., Liu et al., 2008; Moucha et al., 2008). Because predictions for dynamic topography in the region do not exceed ∼1 km since 36 Ma (Moucha et al., 2008; Liu et al., 2010), our overall findings are not compromised because such a signal from dynamic topography would constitute only 25% of the total paleoelevation.

Here H_A is the theoretical crustal thickness, H_SL is the line intercept (crustal thickness at sea level), h_e is elevation, and [ρ_m / (ρ_m – ρ_c)] is the slope of the line, where ρ_c is crustal density and ρ_m is the reference upper mantle density. For this perfect compensation model (Airy compensation), both ρ_c and ρ_m do not vary laterally.

The final integrated result leads to a high Nevadaplano with an average elevation of ∼3.95 ± 0.3 km at 36 Ma for the optimal solution (SEUW = 1.3). This result of a high Nevadaplano at 36 Ma appears to be robust, as it is also obtained for both the underdamped (4.67 ± 0.4 km) and overdamped (3.75 ± 0.2 km) strain solutions (Fig. 13) (see Supplemental Fig. S8 [^{footnote 2}]).

By taking into account the fact that lower crustal flow was likely a factor, the smoothed surface elevation models yield a broader plateau and at slightly lower elevation of ∼3.57 ± 0.3 km (see Supplemental Fig. S9 [^{footnote 2}] and Supplemental Animation S9 [^{footnote 1}]).

Major thermal perturbations occurred in the upper mantle during the Oligocene and Miocene under the Basin and Range, as indicated by the mid-Tertiary ignimbrite flareup and the initiation of Basin and Range extension (Armstrong and Ward, 1991; Humphreys, 1995; Best et al., 2013, 2016). The thermal perturbation associated with this magmatic history likely impacted the density field of the upper mantle through time, as the present-day mantle densities beneath the Great Basin reflect warmer upper mantle that both postdates and likely results from the Farallon slab rollback during the Cenozoic. Therefore, to address the effect of thermal perturbation on the upper mantle density variation and compensated topography, we derive a thermal model for the upper mantle density of the western U.S. using the North American Volcanic and Intrusive Rock Database (NAVDAT; http://www.navdat.org) for magmatism that occurred in the western U.S. from 36 Ma to present. To this end, we use a finite element package to solve the steady-state conductive heat flow equation in COMSOL Multiphysics software (version 5.2) with imposed sources for temperature perturbation linked to the temporal and spatial history of magmatism from the NAVDAT data set (see the Appendix).

Our model for the thermal perturbation history for the upper mantle is an approximation. It is a superposition of four steady-state conductive heat distribution models. The four models divide the 36 m.y. time histories into four roughly equal parts (see Supplemental Fig. S10 [^{footnote 2}], and Supplemental Animation S10 [^{footnote 1}]). Our result can only be an approximation because we are obtaining solutions to the steady-state heat flow problem and not the time-dependent heat flow equations. Moreover, our model history between 36 Ma and 0 Ma is not constrained to match present-day heat flow. Instead, our approach is designed to test the sensitivity of topography models to the lithosphere’s upper mantle density changes arising from the estimates of thermal perturbations associated with the temporal and spatial history of magmatic activity (NAVDAT data set). The areas with densities of ∼3100 kg/m³ for this new upper mantle model at 36 Ma would represent density values following slab rollback, considering that the bulk of the north-to-south slab rollback was complete across southern Nevada (north of 36°N) by 34 Ma (Humphreys, 1995; Dickinson, 2003, 2006; Liu et al., 2010) (see Supplemental Fig. S11 [^{footnote 2}], and Supplemental Animations S11, S12 [^{footnote 1}]).

This new variable upper mantle density model yields a new surface elevation model showing slight differences from the original upper mantle density model, which lacked density changes associated with temperature variations. The main differences are slight increases in Nevadaplano elevation (by ∼1 km) as temperatures increase there (by 150 °C) owing to magmatic activity between 36 and 30 Ma. Following 30 Ma, the effect of topographic collapse dominates over temperature increase, and elevations for the Nevadaplano gradually decrease. There are also slight increases in elevation of the Rocky Mountains and even the Colorado Plateau (by <300 m), owing to volcanic activity in the Rockies and some conductive heat increase in the mantle lithosphere below the Colorado Plateau. One important outcome of our investigation using the NAVDAT data set is that the collapse of topography correlates in space and time with the distribution of active volcanism during the Oligocene and Miocene within the Basin and Range of the western U.S. (see Supplemental Fig. S12 [^{footnote 2}], and Supplemental Animation S13 [^{footnote 1}]). This collapse of topography was presumably facilitated, in part, by the thermal perturbation in the upper mantle associated with the magmatism (see Supplemental Digital Files³).

The timing and extent of the strain history embedded in our models are linked directly with McQuarrie and Wernicke’s (2005) model, who carefully compiled the existing structural data and estimates of net movements through time, along with estimates of uncertainties. Our contribution is to use this information to provide quantitative models of crustal thickness and surface elevation with associated formal uncertainties. Our formal uncertainties do not account for erosion and volume addition from igneous activity. Erosion reduces crustal thickness if the sediments leave the system. Therefore, starting with the present-day crustal thicknesses and integrating them back to estimate original thicknesses provides minimum estimates of pre-erosion crustal thicknesses. However, the basins of the Basin and Range province have significant thicknesses of locally derived sediments, which are part of the present-day crustal thickness distribution estimated through seismology and used in this study (Shen and Ritzwoller, 2016). Therefore, if most of the sediments are locally deposited, we do not expect to grossly underestimate original crustal thicknesses. On the other hand, igneous addition adds material to crustal thicknesses over time. Thus, not accounting for igneous addition causes an overestimation for crustal thickness through time. Because these two unknown factors have opposite signs, they may tend to cancel one another if the mass of sediments that have left the system through erosion is roughly equivalent to the mass of crustal igneous addition.

Our time-dependent crustal thickness model at 36 Ma is consistent with a significant crustal welt at that time, looking remarkably similar to the welt proposed by Coney and Harms (1984) from northern to southern Nevada (Nevadaplano), eastern California, and Arizona. However, our model does not predict the 50–55 km welt north of 38°N within eastern California predicted by Coney and Harms (1984) (Fig. 14).

A roughly 55-km-thick crust in Nevada is also consistent with high Sr/Y ratios of Jurassic–Cretaceous intermediate continental calc-alkaline magmatic rocks (Chapman et al., 2015). Decreasing Sr/Y ratios suggest mid-Eocene to Oligocene extension and decreased crustal thickness to 30–40 km by the Miocene as a result of extension (Chapman et al., 2015).

Such a thick crustal welt at around 36 Ma could have driven extension, crustal thinning, and collapse of topography, resulting ultimately in the present-day geometry of the Great Basin and Basin and Range provinces. Because this crustal welt would have been present long before the initiation of collapse (Sonder et al., 1987; DeCelles, 2004; Liu and Gurnis, 2010), such a collapse was probably initiated by reduction of viscosity by a mantle-derived heating event (Liu and Shen, 1998; Liu, 2001), thermal relaxation of the overthickened crust (Sonder et al., 1987), or collapse and steepening of a previously shallowly dipping Laramide Benioff zone, which may have reduced the regional stress and possibly started extension (Coney, 1980; Armstrong, 1982; Spencer, 1984; Humphreys, 1995; Jones et al., 1996; Dickinson, 2002, 2003).

Armstrong (1982), Coney and Harms (1984), and Parrish et al. (1988) suggested that metamorphic core complexes were extensional in origin and mainly Tertiary in age (65–2 Ma). Coney and Harms (1984) argued that the metamorphic core complexes formed above sites of extreme crustal thickening. Tracking the positions of core complexes back to 36 Ma locates them above the thickest crust and the Nevadaplano (Fig. 14) (see Supplemental Animations S14, S15 [^{footnote 1}]), confirming the findings of Coney and Harms (1984). However, Coney and Harms (1984) predicted a large crustal welt in eastern California (north of 38°N) where no core complexes are observed. By contrast, our model predicts only a modest crustal welt there.

Our final integrated topography model shows a highland with an average elevation of ∼3.95 ± 0.3 km in central, eastern, and southern Nevada, western Utah, parts of easternmost California, and for northwestern Arizona (Nevada plano). The Mogollon Highlands (Cooley and Davidson, 1963; Elston and Young, 1991) are also present within central and southeastern Arizona at 36 Ma (Fig. 11). This topographic belt stretching from northern Nevada to southeastern Arizona and northern Mexico results from a significant crustal welt that was likely a consequence of Sevier-Laramide convergence (Dickinson, 2002; DeCelles, 2004; Sigloch and Mihalynuk, 2013).

The calculated surface elevation estimates at ca. 36 Ma represent elevations following slab rollback, where the crustal root supporting much of that elevation was inherited prior to 36 Ma from Sevier-Laramide convergence (DeCelles et al., 1995; Dickinson, 2002; DeCelles, 2004). Our surface elevation models show that the crustal topography changes of the Nevadaplano, owing to the transition in upper mantle density associated with slab rollback and thermal perturbation, constitutes only ∼25% of the total elevation of the Nevadaplano. This heating of the upper mantle lithosphere accompanied the migration of volcanism following slab rollback and is represented as a heating event in our models (temperature increase of 150 °C in upper mantle lithosphere) between 36 Ma and 30 Ma. This temperature increase results in a reduction of 60 kg/m³ in upper mantle density, which yields an elevation increase of ∼1 km for Nevadaplano. Whereas 1 km is a substantial uplift, the bulk of the elevation of the Nevadaplano (average 3.95 ± 0.3 km) results from the crustal welt as opposed to the upper mantle density changes associated with slab rollback and thermal perturbation (e.g., Mix el al., 2011).

Our topography model also shows that between 36 and 5 Ma, the Sierra Nevada was located adjacent to paleo–sea level, with a shoreline on the eastern edge of the present-day Great Valley (Fig. 11). The model for the northern Sierra Nevada (north of 36°N) shows little to no elevation change throughout the deformation period since at least 36 Ma (Fig. 11), consistent with idea that the majority of uplift occurred prior to 36 Ma, in the Late Cretaceous to early Cenozoic (Cassel et al., 2009; Cecil et al., 2010). However, by contrast, our model for the southern Basin and Range at 14 Ma involves a steep gradient in crustal thickness from ∼30 km in eastern California to ∼50–55 km just east of the California–southern Nevada border region (Fig. 7). The topography model for this time period shows variable elevation in eastern California ranging from 0.5 to 1.5 km and reaching high elevations of ∼3.5–4 km in the area of the 50–55-km-thick crust in southern Nevada (Fig. 11).

Our paleotopography estimates can be compared to a variety of proxies that have been employed by many studies to estimate paleotopography and paleoaltitude across the western U.S. Because the isotopic composition of precipitation somewhat directly scales with elevation, it can be used to reconstruct topographic histories of mountain belts (Mix et al., 2011), though the basins in which the carbonates that can be analyzed are found clearly represent a minimum elevation. Maps of the δ¹⁸O of precipitation for different time bins during the Cenozoic highlight the major isotopic shifts observed within individual basins along the Cordillera (Mix et al., 2011; Horton et al., 2004; Kent-Corson et al., 2006). Using spatial and temporal patterns of δ¹⁸O of precipitation, Mix et al. (2011) determined paleoelevations for Elko Basin in northeastern Nevada (36–28 Ma, ∼3.4 km), Copper Basin in northeastern Nevada (37 Ma, ∼3.4 km), and the Sage Creek Basin in southwestern Montana (38.8–32.0 Ma, ∼3.7 km). They also show an Eocene–Oligocene highland of ∼3.4 km in eastern Nevada and western Utah. The paleobotanical study by Wolfe et al. (1998) also suggests the presence of an Eocene–Oligocene highland of ∼4 km or higher in the Copper Basin region (northeastern Nevada).

Paleobotanical results for the late Eocene House Range flora in the Sevier Desert (west-central Utah) yield a paleoelevation of ∼3–4 km at ca. 31 Ma (Gregory-Wodzicki, 1997). Based on δD measurements in hydrated glass from ignimbrites, Cassel et al. (2014) indicated the presence of a high, broad orogen that stretched across northern to southern Nevada during the Eocene to Oligocene, with the highest elevations of 3.5 km in the late Oligocene. Gébelin et al. (2012) results based on δ¹⁸O and δD values calculated from muscovite from Snake Range Mountain (northwestern Nevada) yield an elevation of ∼3850 ± 650 m between 27 and 20 Ma. These are in agreement with our paleotopography models for 36–20 Ma, showing a Nevadaplano with an average elevation of ∼3.95 ± 0.3 km in central, eastern, and southern Nevada, western Utah, and parts of easternmost California (Fig. 11).

Paleobotanical analyses obtained by Wolfe et al. (1997) for several mid-Miocene floras in eastern Nevada suggest a paleoaltitude of ∼3 km at ca. 15–16 Ma. Their assumption that the highland had collapsed by ca. 13 Ma is in agreement with our topography model showing the altitudes similar to the present-day in eastern Nevada at 13 Ma and afterward (Fig. 11).

Mix et al. (2011) suggested that slab rollback may be the primary source of a wave of uplift of the Nevadaplano that swept from north to south during the Oligocene. They argued that the southernmost measurement achieved an elevation change of ∼2.5 km late in the time between 36 and 28 Ma, possibly when slab rollback was occurring beneath the region. As we have explained above, we estimate a 60 kg/m³ change in upper mantle density beneath a ∼55-km-thick crust associated with thermal heating, and this produces only ∼1 km of uplift. Although Mix et al. (2011) may have indeed resolved a component of uplift associated with slab rollback, we show that elevations are expected to have still been high (in excess of 3 km) above the substantial crustal welt that was present just prior to the slab rollback. Another potential source of uplift prior to 36 Ma may have resulted from hydration of the lithosphere from the Farallon slab during the flat slab subduction phase (Humphreys et al., 2003; Jones et al., 2015; Porter et al., 2017).

Cassel et al. (2009) also showed that northern Sierra Nevada elevation (100 km east of the paleo-coastline) in the Oligocene was similar to that of the present day (∼2800 km). This is consistent with results of Chamberlain and Poage (2000) who measured little change in δ¹⁸O of smectites from the northern Sierra Nevada since 16 Ma, suggesting that the Sierra Nevada has been approximately at the same elevation since 16 Ma. This is in agreement with an elevation of ∼2.5 km shown in our model from 36 Ma to the present day. Our model also shows a Sierra Nevada adjacent to the paleo–sea level shoreline before 5 Ma. This is consistent with the idea that volcanic materials originating in northern and central Nevada were deposited during the Paleogene into the Pacific Ocean within what is now the Great Valley (Faulds et al., 2005; Garside et al., 2005; Henry, 2009).

Clumped isotope (Δ47) thermometry of lacustrine carbonates (∼17–24 °C) from the central Basin and Range and the southern Sierra Nevada Bena Basin indicate that middle Miocene paleoelevations in the Death Valley region were <1.5 km (Lechler et al., 2013), consistent with our estimates. However, our results do not agree with the 1.5 km elevations in the middle Miocene for the Spring Mountains region obtained by Lechler et al. (2013). They acknowledged that the 1.5 km elevation was exceptionally low for the expected pre-extensional crustal thicknesses there (>52 km), and that final post-extensional elevations can only be explained by igneous volume addition. Supporting evidence for the topography in our model lies in the agreement we see for the exhumation of core complexes along the belt near the border region between eastern California and southernmost Nevada. Possible discrepancies may be explicable by the sampling of lacustrine carbonates within low-lying basins that are not resolvable with our modeling technique, which addresses the smoothed components of the deformation field and thus smoothed topography. Huntington et al.’s (2010) paleoelevation estimates from apatite (U-Th)/He data of the Grand Canyon reveal that most of the Colorado Plateau’s lithospheric buoyancy and uplift occurred at ca. 80–60 Ma. This is in agreement with our surface elevation model showing almost no elevation change for Colorado Plateau at least during the past 36 m.y.

Our paleotopography reconstruction models provide a framework for considering the timing of the topographic development of the Sierra Nevada and Basin and Range province. Future models can be refined through the incorporation of results from isotopic and thermochronological studies.

We use the compilation of McQuarrie and Wernicke (2005) to determine a horizontal velocity gradient tensor field solution for the lithosphere through time. Using present-day crustal thickness estimates for initial conditions, and assuming volume conservation, we integrate these solutions to produce models of crustal thickness and surface elevation, along with formal uncertainties, through the past 36 m.y. Our crustal thickness model at 36 Ma is consistent with a significant crustal welt, with an average crustal thickness of 56.5 ± 2.5 km, looking similar to the welt proposed by Coney and Harms (1984) from northern to southern Nevada and through central and southeast Arizona, which was likely a consequence of Sevier-Laramide convergence. Such a thick crustal welt at ca. 36 Ma could have played a major role in driving extension, crustal thinning, and collapse of topography, resulting ultimately in the present-day geometry of the Great Basin and Basin and Range provinces.

Our final integrated topography model shows a Nevadaplano of ∼3.95 ± 0.3 km average elevation in central, eastern, and southern Nevada, western Utah, parts of easternmost California, and northwestern Arizona. Highlands of significant elevation (∼3–3.5 km) are also present along a belt in central and southeastern Arizona at 36 Ma (Mongollon Highlands). Our model, between 36–5 Ma, shows a Sierra Nevada adjacent to paleo–sea level, with a shoreline on the eastern edge of the present-day Great Valley. Moreover, based on our model, the Colorado Plateau and the northern Sierra Nevada (north of 36°N) show little to no elevation change since at least 36 Ma.

We would like to express our gratitude to reviewers Brian Wernicke, an anonymous reviewer, and Science Editor Raymond Russo for their careful reading of our manuscript and their many insightful comments and valuable suggestions, which greatly improved the quality of the manuscript. We thank Yuanyuan Liu and Rubin Smith for assistance with early calculations during the initial stages of this work. We are very grateful to Weisen Shen for kindly providing his most recent database for crustal structure of the United States. Grateful thanks are also expressed to Nadine McQuarrie, Daniel Davis, Catherine Badgley, Tara Smiley, Ryan Porter, and Gavin Piccione for their general interest, and fruitful discussions, which greatly assisted in the development and refinement of this study. This research was supported by the National Science Foundation under grant numbers EAR-1246971 and EAR-1052989, and NASA ESI under grant number NNX16AL18G, as well as Southern California Earthquake Center under grant numbers 16291 and 14226. Some data products used in this study were made possible through EarthScope (www.earthscope.org; EAR-0323309), supported by the National Science Foundation. All figures provided in this manuscript are created using GMT (Wessel et al., 2013).

To calculate the formal standard error for crustal thickness and surface elevation through time, we start with present-day crustal thickness of the western U.S. (Shen and Ritzwoller, 2016; Laske et al., 2013), and then we calculate the crustal thickness and surface elevation variations back in time to 36 Ma for every 0.5 m.y. time increment.

We generate time-dependent two-dimensional (2-D) steady-state conductive solutions of heat distribution for each 0.5 m.y. (k = 72−1) from 36 Ma to present day based on four different patterns (i = 1–4) of magmatism during the past 36 m.y.: 72 < k < 61, i = 1; 60 < k < 37, i = 2; 36 < k < 17, i = 3; and 16 < k < 0, i = 4 (see Supplemental Fig. S10 [^{footnote 2}]).

We use 150 °C as the differential temperature for the areas affected by magmatism and ignimbrite flareups. This choice of a thermal perturbation linked to active magmatism comes from constraints for present-day upper mantle lithosphere obtained by Schutt et al. (2012). Their study, inferred from present-day seismic and heat flow constraints, shows an average of ∼150 °C temperature variation of upper mantle lithosphere between undeformed blocks like the Great Plains and Colorado Plateau and magmatically and tectonically active regions such as Yellowstone, the margin of the Colorado Plateau, and the western margin of the Great Basin. To compute density variations owing to an increase in lithosphere upper mantle temperature of this amount, we only need to consider solutions for a relative temperature increase. Thus, we define boundary conditions of ΔT = 0 °C around our grid area, and also for the continental margin of the western U.S., as well as ΔT = 150 °C for the heat sources, based on the NAVDAT data set of magmatic source locations (see Supplemental Fig. S10 [^{footnote 2}]). The temperature distribution models for upper mantle are produced in COMSOL. We then advect the heat distribution model for each solution time interval (i = 1–4), and produce time-dependent temperature models for the upper mantle density.

Omitting details, we constrained the starting mantle density values prior to 36 Ma to be such that the thermal perturbation history gives back the present-day density values for the upper mantle, obtained through compensation of present-day topography. That is, the starting density values prior to 36 Ma are higher and the thermal perturbations through time cause density reductions over time, particularly in the vicinity of magmatic activities (Supplemental Fig. S11 [^{footnote 2}], and Supplemental Animations S11, S12 [^{footnote 1}]). Based on the differential temperature solutions for each 0.5 m.y. [ΔT(ϕ_k, θ_k)] and the thermal coefficient of expansion for the upper mantle, and using Equation A19, we calculate a new time-dependent upper mantle density [ρ_T(ϕ_k, θ_k)] for the western U.S. from 36 Ma to the present day.