The complex deformation history of the western United States since 36 Ma involves a dramatic transition from a subduction-dominated to a transform-dominated margin. This transition involved widespread extension, collapse of topography, and development of the interior Basin and Range region. The topographic collapse resulted in significant exhumation of deep crustal rocks exposed in metamorphic core complexes of the southwestern Cordillera. We use calculated position estimate changes for the western United States from previous work based on a comprehensive compilation of geological and structural information, and incorporation of constraints from Pacific plate motion history to determine lithospheric strain rates through time, and integrate these strain rates, to provide quantitative models of crustal thickness and surface elevation, along with formal standard errors, since 36 Ma. Our crustal thickness model at 36 Ma is consistent with a significant crustal welt, with an average thickness of ∼56.5 ± 2.5 km, in eastern portions of northern to southern Nevada as well as parts of eastern California and through Arizona. Our final integrated topography model shows a Nevadaplano of ∼3.95 ± 0.3 km average elevation in central, eastern, and southern Nevada, western Utah, and parts of easternmost California. A belt of high topography also trends through northwestern, central, and southeastern Arizona at 36 Ma (Mogollon Highlands). Our model shows little to no elevation change for the Colorado Plateau and the northern Sierra Nevada (north of 36°N) since at least 36 Ma, and that between 36 and 5 Ma, the Sierra Nevada was located at the Pacific Ocean margin, with a shoreline on the eastern edge of the present-day Great Valley.

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).

Figure 1.

Present-day geography of the western United States and the location of the Basin and Range province including the Great Basin and southern Basin and Range regions (gray area).

Figure 1.

Present-day geography of the western United States and the location of the Basin and Range province including the Great Basin and southern Basin and Range regions (gray area).

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.

Computation of a Velocity Gradient Tensor Field through Time within the Basin and Range

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.

Figure 2.

Western U.S. grid showing the position of fault displacement and velocity observation sites from McQuarrie and Wernicke (2005) (blue dots), and the location of points added in this study as constraints in the Colorado Plateau, eastern Rio Grande Rift, and Pacific Plate (gray dots). The red line represents the San Andreas fault system. The grid includes 2475 knotpoints and 1714 degrees of freedom associated with velocities.

Figure 2.

Western U.S. grid showing the position of fault displacement and velocity observation sites from McQuarrie and Wernicke (2005) (blue dots), and the location of points added in this study as constraints in the Colorado Plateau, eastern Rio Grande Rift, and Pacific Plate (gray dots). The red line represents the San Andreas fault system. The grid includes 2475 knotpoints and 1714 degrees of freedom associated with velocities.

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 forumla, and rotation rates about the vertical axis forumla (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 S11).

Figure 3.

Time-slice maps of western U.S. velocity field estimates relative to the North American frame. Contoured dilatation strain rates are obtained from inversion of position change estimates of McQuarrie and Wernicke (2005) using the inversion algorithm of Beavan and Haines (2001). Red vectors are model velocities (95% confidence error ellipse) at the coordinate sites where position change estimates have been obtained from the kinematic model of McQuarrie and Wernicke (2005). Green vectors are model velocities for the Pacific plate relative to the stable North American frame. The blue line represents the location of the continental margin through time, and the red line represents the location of the East Pacific Rise and, in later times, the location of the San Andreas fault system. Contour interval (white) is 25 × 10–9/yr.

Figure 3.

Time-slice maps of western U.S. velocity field estimates relative to the North American frame. Contoured dilatation strain rates are obtained from inversion of position change estimates of McQuarrie and Wernicke (2005) using the inversion algorithm of Beavan and Haines (2001). Red vectors are model velocities (95% confidence error ellipse) at the coordinate sites where position change estimates have been obtained from the kinematic model of McQuarrie and Wernicke (2005). Green vectors are model velocities for the Pacific plate relative to the stable North American frame. The blue line represents the location of the continental margin through time, and the red line represents the location of the East Pacific Rise and, in later times, the location of the San Andreas fault system. Contour interval (white) is 25 × 10–9/yr.

In the inversion procedure, the objective functional that is minimized is the following:
graphic
where σN and σE are the formal standard errors in velocity for north and east directions, respectively; here, VNobs is the observed velocity in north direction, VEobs is the observed velocity in east direction, and velocities superscripted with “mod” represent the model velocity for that component. Positive directions are east and north, respectively. In the second part of the functional, forumla are the model values of strain rates, and υ is a weighting factor. In practice, the second part of the functional, involving the second invariant of only model strain rates, is achieved by defining zero values for observed strain rates within all areas (Beavan and Haines, 2001). Note that the weighting factor υ consists of probabilistic constraints on these zero values of observed strain rates within areas, where:
graphic

Here, forumla) is the variance for zero-valued observed strain rates; Si is the area of the ith 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.

The shape of the continuous model strain rate field is controlled by how close we want to match the velocities. If the observed and modeled velocities are matched closely, then the field might be fairly irregular; however, if larger misfits are allowed with the observed velocities, the model field will be smoother. Hence, choosing a proper value for υ is of high importance. Beavan and Haines (2001) argued that an optimal value of υ was one that gave the sum of the squares misfit of the velocity components, SS(v), equal to the number of degrees of freedom (Ndof). Because there are two degrees of freedom associated with each horizontal velocity estimate, Ndof = 1714 for 857 velocity observation sites. However, we settle on an optimum value for υ that provides sum of the squares divided by the number of degrees of freedom approximately equal to 1.5. This allows for the fact that the velocity errors may be underestimated, and results in a smoother model solution. We quantify the misfit using the a posteriori formal standard error of unit weight (SEUW) (Beavan and Haines, 2001), which is:
graphic

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.

Figure 4.

Maps of the present-day western U.S. velocity field estimates relative to the North American frame with contoured dilatation strain rates. (A) Overdamped solution with standard error of unit weight (SEUW) = 1.7. (B) Intermediate damped solution with SEUW = 1.3. (C) Underdamped solution with SEUW = 1.0. The vectors are model velocities (95% confidence error ellipse) at the coordinate sites from the kinematic model of McQuarrie and Wernicke (2005). The green and red vectors are the same as in Figure 3.

Figure 4.

Maps of the present-day western U.S. velocity field estimates relative to the North American frame with contoured dilatation strain rates. (A) Overdamped solution with standard error of unit weight (SEUW) = 1.7. (B) Intermediate damped solution with SEUW = 1.3. (C) Underdamped solution with SEUW = 1.0. The vectors are model velocities (95% confidence error ellipse) at the coordinate sites from the kinematic model of McQuarrie and Wernicke (2005). The green and red vectors are the same as in Figure 3.

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 (forumla), the greater the degree of smoothing (with lower resolution), the smaller the formal errors (see Supplemental Animation S2 [footnote 1]).

Figure 5.

Maps of formal standard error for model contoured dilatation strain rates (see Fig. 4). (A) Overdamped solution with standard error of unit weight (SEUW) = 1.7. (B) Intermediate solution with SEUW = 1.3. (C) Underdamped solution with SEUW = 1.0. Note that each figure has its own specific range for the color scale. Dots are crustal velocity constraints obtained from McQuarrie and Wernicke (2005) as in Figure 2.

Figure 5.

Maps of formal standard error for model contoured dilatation strain rates (see Fig. 4). (A) Overdamped solution with standard error of unit weight (SEUW) = 1.7. (B) Intermediate solution with SEUW = 1.3. (C) Underdamped solution with SEUW = 1.0. Note that each figure has its own specific range for the color scale. Dots are crustal velocity constraints obtained from McQuarrie and Wernicke (2005) as in Figure 2.

Calculating Crustal Thickness and Position Changes through Time

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 forumla. 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.

Calculating Crustal Thickness

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).

In general, the relationship between a vector y in a deforming region at time t = 0, and time t, is:
graphic
where F(t) is the deformation gradient tensor (McKenzie and Jackson, 1983). For our specific case, y would be a vertical vector representing crustal thickness, H.
The time dependence of F is:
graphic
where L is the velocity gradient tensor (McKenzie and Jackson, 1983). We are interested in the component Fzz for a strain rate field that is considered vertically coherent. Using this vertical coherency for the vertical strain rates (forumla), we have from Equation 5:
graphic
Solving this first-order differential equation, we have forumla, where Δt is the incremental time step (0.5 m.y.). From Equation 5, the time-dependent crustal thickness, Ht, is
graphic
where H0 is the crustal thickness at the beginning of the time step, and the vertical strain rate (forumla) is appropriate for a specific time interval and coordinate position within the generally spatially varying (latitude and longitude only) strain rate field (Fig. 6).
Figure 6.

Schematic drawings of the finite strain estimate inferred from the distribution of vertical strain rates. (A) Instantaneous strain rate distribution. (B) Change in crustal thickness obtained from integration of vertical strain rates back in time over a specific time interval (Δt). Δx and Δy are the changes in length and width, respectively. forumla are the incremental changes in x and y directions, respectively. forumlazz is the vertical strain rate. H0 is the present-day crustal thickness. H1 is the newly calculated time-dependent crustal thickness for each 0.5 myr time interval. Green arrows depict extension direction in A, and reverse motion associated with reconstruction through time in B. See the text (Equations 7 and 8) for details.

Figure 6.

Schematic drawings of the finite strain estimate inferred from the distribution of vertical strain rates. (A) Instantaneous strain rate distribution. (B) Change in crustal thickness obtained from integration of vertical strain rates back in time over a specific time interval (Δt). Δx and Δy are the changes in length and width, respectively. forumla are the incremental changes in x and y directions, respectively. forumlazz is the vertical strain rate. H0 is the present-day crustal thickness. H1 is the newly calculated time-dependent crustal thickness for each 0.5 myr time interval. Green arrows depict extension direction in A, and reverse motion associated with reconstruction through time in B. See the text (Equations 7 and 8) for details.

We start with a dense set of coordinates of 0.1° spacing (longitude = ϕ, latitude = θ) and calculate a shift in position (Δϕ, Δθ) relative to the stable North American frame, using a time increment of Δt = −0.5 m.y. and the instantaneous velocity gradient tensor field solution. For the kth – 1 time interval, the position shift at the kth time interval is Δϕ = ϕk – ϕk–1, Δθ = θk – θk–1, where ϕk = ϕk–1 + Δϕ, and θk = θk–1 + Δθ. The crustal thickness (H) for the kth time interval is:
graphic
where forumla is the vertical strain rate for the jth velocity gradient tensor field solution, where k = 1–72 and j = 1–12. There are 72 time steps (k) because the interval over which we have strain rate constraints is 36 m.y. and the time increment is Δt = −0.5 m.y. There are 12 velocity gradient tensor field solutions (j), which cover the time periods of 0–2, 2–4, 4–6, 6–8, 8–10, 10–12, 12–14, 14–16, 16–18, 18–24, 24–30, and 30–36 Ma. Our starting crustal thickness data set, H00; θ0), is from Shen and Ritzwoller (2016), and all model solutions for crustal thickness and position at a given time in the past can only be obtained by starting at k = 1, j = 1 (present-day), and integrating backward in order. Because the Shen and Ritzwoller (2016) crustal thickness model covers only the entire U.S. and not Mexico (southern Basin and Range area) and offshore (oceanic plate), we merge the crustal thickness data set for Mexico and offshore from the CRUST 1.0 global crustal model (Laske et al., 2013) with the crustal thickness data set for western U.S. from Shen and Ritzwoller (2016) (Fig. 7).
Figure 7.

Time-slice maps of crustal thickness (in km) for the Basin and Range of the western United States, based on integrating strain rates backward to estimate crustal thickness over time. State boundaries are moved with all points for reference.

Figure 7.

Time-slice maps of crustal thickness (in km) for the Basin and Range of the western United States, based on integrating strain rates backward to estimate crustal thickness over time. State boundaries are moved with all points for reference.

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 kt = −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. S12, 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]).

Figure 8.

Time-slice maps of standard error estimates (in km) of crustal thickness for the Basin and Range of the western United States. State boundaries are moved for each time step as in Figure 7.

Figure 8.

Time-slice maps of standard error estimates (in km) of crustal thickness for the Basin and Range of the western United States. State boundaries are moved for each time step as in Figure 7.

Calculating Upper Mantle Density and Surface Elevation

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.

However, we adopt a simple model that assumes pressure equilibrium at an upper mantle depth of 100 km, consistent with the inference of negligible elevation variations owing to dynamic topography within the western U.S. (Levandowski et al., 2014). We start by investigating the present-day topography and crustal thickness estimates. We plot topography (ETOPO5 elevation data; https://www.ngdc.noaa.gov/mgg/global/etopo5.HTML) versus crustal thickness (Fig. 9A). If all topography points were compensated at 100 km depth relative to a reference column with a non-variable upper mantle density (e.g., 3300 kg/m3), then all points in a crustal thickness versus elevation plot would lie on a line defined by:
graphic
Figure 9.

(A) Scatterplot and a line of best fit based on present-day crustal thickness data (Shen and Ritzwoller, 2016) and ETOPO5, a global map of topographic elevation (https://www.ngdc.noaa.gov/mgg/global/etopo5.HTML). (B) Schematic illustration for our upper mantle density calculation and compensation of topography. H0 is the present-day crustal thickness, HA is the theoretical crustal thickness, HSL is the line intercept (crustal thickness at sea level), he is elevation, ρc is crustal density, ρm is the reference upper mantle density, ρm′ is upper mantle density, and D is the distance from base of crust to the compensation depth of 100 km below sea level.

Figure 9.

(A) Scatterplot and a line of best fit based on present-day crustal thickness data (Shen and Ritzwoller, 2016) and ETOPO5, a global map of topographic elevation (https://www.ngdc.noaa.gov/mgg/global/etopo5.HTML). (B) Schematic illustration for our upper mantle density calculation and compensation of topography. H0 is the present-day crustal thickness, HA is the theoretical crustal thickness, HSL is the line intercept (crustal thickness at sea level), he is elevation, ρc is crustal density, ρm is the reference upper mantle density, ρm′ is upper mantle density, and D is the distance from base of crust to the compensation depth of 100 km below sea level.

Here HA is the theoretical crustal thickness, HSL is the line intercept (crustal thickness at sea level), he 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.

If points do not lie on the line, then the crustal density, the mantle density, or both do not match with the reference model. While it is recognized that there are lateral variations in the average crustal density within the western U.S. (Levandowski et al., 2014), in our simple approach, we assume that all pressure variations are due to upper mantle densities relative to the reference model, and that the crustal density is constant. For example, if topography points lie above the best-fit line (Fig. 9), then they are supported by mantle densities that are greater than the reference density; points below the line are supported by mantle densities less than that of the reference model. Using an average crustal density of 2700 kg/m3 for the western U.S. (Mooney and Kaban, 2010; Levandowski et al., 2014), the slope of the best-fit line in Figure 9A yields a best-fit reference upper mantle density of 3150 kg/m3, and the intercept defines HSL of 27 km (thickness at sea level). As part of the assumption that topography is compensated at 100 km, we solve for the upper mantle density for each point. Allowing for a variable upper mantle density (ρm′), the revised relationship between actual crustal thickness (H0) (Shen and Ritzwoller, 2016), crustal density (ρc), elevation (he), reference upper mantle density (ρm), and upper mantle density (ρm′) is:
graphic
where D is the distance from base of crust to the compensation depth of 100 km below sea level (Fig. 9B). By subtracting Equation 9 from Equation 10, we solve for the variable upper mantle density (ρm′):
graphic
Solving for ρm′ yields a coherent pattern that shows lower densities in the Great Basin and southern Basin and Range provinces, normal densities in the Colorado Plateau, and higher densities in the Great Plains (Fig. 10) (see Supplemental Animation S6 [footnote 1]). The patterns in this upper mantle density model are consistent with patterns of seismic velocity structure for uppermost mantle at depths of 90–100 km in the western U.S. (Yang et al., 2008; Schmandt and Humphreys, 2010; Moschetti et al., 2010; Obrebski et al., 2011; Schmandt et al., 2015; Porter et al., 2016; Shen and Ritzwoller, 2016). That is, comparing with these seismic velocity models, the present-day density model (Fig. 10) shows values lower than the reference value in areas corresponding to slower-than-average upper mantle S- and P-wave velocities (Great Basin, southern Basin and Range, southern margin of the Colorado Plateau). Likewise, higher densities correspond to areas where upper mantle P- and S-wave velocities are generally higher (northern Colorado Plateau, western edge of the Great Plains). Our present-day density model is also in agreement with upper mantle temperature models (Lowry and Pérez-Gussinyé, 2011; Afonso et al., 2016). Having the crustal thickness and upper mantle density, using Equation 10, and because D = 100 km – H0 + he, we solve for the elevation (he) as:
graphic
Figure 10.

Time-slice maps of western United States upper mantle density (ρm′) calculated based on compensation at 100 km depth below sea level. Positions of density values are moved to new coordinates using the time-dependent velocity gradient tensor field solutions (see text). State boundaries are moved for each time step as in Figure 7.

Figure 10.

Time-slice maps of western United States upper mantle density (ρm′) calculated based on compensation at 100 km depth below sea level. Positions of density values are moved to new coordinates using the time-dependent velocity gradient tensor field solutions (see text). State boundaries are moved for each time step as in Figure 7.

As a first approach, to calculate the paleoelevation, we make the simplifying assumption that ρm′ does not change throughout time. We later discuss the effect of thermal perturbation on a new time-dependent upper mantle density model, with the same depth of compensation, and recalculate the paleoelevation based on the variation of ρm′ through time. We track the position of the ρm′(ϕ, θ) values (Fig. 10) in the same way that we track all crustal thickness points through time. Then, by substituting the time-dependent crustal thickness from Equation 8 into Equation 12, we arrive at the elevation (he) for the kth time interval (Fig. 11):
graphic
where again the integration must start at k = 1, j = 1, and proceed in order. Using the formal estimates of standard errors for crustal thickness enables us to calculate formal standard errors for surface elevation at all time steps (see the Appendix) (Fig. 12) (see Supplemental Animations S7, S8 [footnote 1]). To test the sensitivity to compensation depth, we solve for ρm′ and surface elevation assuming a depth of compensation of 150 km. This yields a mean value of 3170 kg/m3 for ρm′ and a slightly narrower range (3120–3220 kg/m3), and a nearly identical prediction for surface elevation through time. We also analyzed the CRUST 1.0 model (Laske et al., 2013) for comparison with results obtained using the model of Shen and Ritzwoller (2016) (see Supplemental Figs. S4, S5, S6, S7 [footnote 2]).
Figure 11.

Time-slice maps for the surface elevation (in m) of the western United States, achieved from compensation of crustal thickness. State boundaries are moved for each time step as in Figure 7.

Figure 11.

Time-slice maps for the surface elevation (in m) of the western United States, achieved from compensation of crustal thickness. State boundaries are moved for each time step as in Figure 7.

Figure 12.

Time-slice maps of surface elevation formal standard errors (in m) for the Basin and Range of the western United States. State boundaries are moved for each time step as in Figure 7.

Figure 12.

Time-slice maps of surface elevation formal standard errors (in m) for the Basin and Range of the western United States. State boundaries are moved for each time step as in Figure 7.

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]).

Figure 13.

Test of sensitivity for velocity field (top), final crustal thickness (middle), and topography estimates (bottom). Model velocity field (red vectors with 95% confidence) and dilatation strain rates for (A) overdamped solution with standard error of unit weight (SEUW) = 1.7; (B) optimal solution with SEUW = 1.3; and (C) underdamped solution with SEUW = 1.0. Note the presence of a high Nevadaplano for all solutions at 36 Ma. State boundaries are moved from their present-day positions as in Figure 7.

Figure 13.

Test of sensitivity for velocity field (top), final crustal thickness (middle), and topography estimates (bottom). Model velocity field (red vectors with 95% confidence) and dilatation strain rates for (A) overdamped solution with standard error of unit weight (SEUW) = 1.7; (B) optimal solution with SEUW = 1.3; and (C) underdamped solution with SEUW = 1.0. Note the presence of a high Nevadaplano for all solutions at 36 Ma. State boundaries are moved from their present-day positions as in Figure 7.

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]).

Influence of Thermal Perturbations on Western U.S. Upper Mantle Densities

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/m3 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 Files3).

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.

Crustal Thickness Model of the Western U.S.

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).

Figure 14.

(A) Map of crustal thickness at 36 Ma from integration of strain rates, showing a crustal welt in the western United States, similar to that proposed by Coney and Harms (1984). (B) Topography at 36 Ma from integration of strain rates, showing a high Nevadaplano with average elevation of ∼3.95 ± 0.3 km. Red dots are locations of metamorphic core complexes (Dickinson, 2002). State boundaries are moved from their present-day positions as in Figure 7.

Figure 14.

(A) Map of crustal thickness at 36 Ma from integration of strain rates, showing a crustal welt in the western United States, similar to that proposed by Coney and Harms (1984). (B) Topography at 36 Ma from integration of strain rates, showing a high Nevadaplano with average elevation of ∼3.95 ± 0.3 km. Red dots are locations of metamorphic core complexes (Dickinson, 2002). State boundaries are moved from their present-day positions as in Figure 7.

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.

Paleotopography Model of the Western U.S.

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/m3 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 δ18O 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 δ18O 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 δ18O 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/m3 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 δ18O 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).

PROPAGATION OF UNCERTAINTIES FOR CRUSTAL THICKNESS AND SURFACE ELEVATION MODELS OF WESTERN U.S.

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.

Computation of Formal Standard Error for Crustal Thickness through Time

As mentioned in the paper, strain rates and positions are tracked in temporal order, starting at the present day. The vertical strain rates for the kth position and time (ϕk = ϕk–1 + Δϕ; θk = θk–1 + Δθ, where ϕ is longitude, θ is latitude), and the jth strain solution, are denoted as forumla, where j = 1–12, k = 1–72, and time increment Δt = −0.5 m.y. When solutions are tracked in order, starting at k = 1, the crustal thickness at the kth time interval (Hk) is:
graphic
The formal variances in strain rates that are used in our error propagation calculation involve the spatially averaged quantities within each 1° × 1° area on the grid. As mentioned in the paper, this sized area is reasonable as it reflects our confidence in spatially averaged quantities and also represents the current resolution in crustal thickness determination inferred from seismology. Determining error propagation involves the need for estimating variances for each incremental change in crustal thickness (ΔH = HkHk–1) for the kth time step and the jth solution, where:
graphic
Using the general rule for calculating variance for a multivariable nonlinear function (Snedecor and Cochran, 1994), the variance for the change in crustal thickness from the kth – 1 to the kth time interval is:
graphic
The partial derivative of function F with respect toforumla:
graphic
and the partial derivative of function F with respect to Hk–1 is:
graphic
To calculate forumla, we use the formal variances forumla and covariances forumla for the spatially averaged strain rates within each 1° × 1° area of the grid for the jth solution, obtained in the least-squares inversion (Beavan and Haines, 2001). Considering the volumetric strain and volume conservation, the variance for vertical strain rate forumla is:
graphic
To calculate covariance between forumla and Hk, we sample the strain rate field and crustal thickness field on a dense set of points such that the covariance within each 1° × 1° area is defined as:
graphic
where N is the total number points in each 1° × 1° area in our western U.S. grid area, forumla is the average of vertical strains in each area on the grid, and forumla is the average of crustal thicknesses in each area at the kth time interval. The variance in crustal thickness at time t = ΔtN involves the sum of variances of incremental changes in crustal thicknesses for each time step, where N is the total number of time steps under consideration. By substituting Equations A4 and A5 into Equation A3, we arrive at the variance for crustal thickness at the kth time interval, var(Hk):
graphic
Finally, the formal standard error of crustal thickness at kth time interval [se(Hk)] is (Fig. 8):
graphic

Computation of Formal Standard Error for Surface Elevation through Time

To calculate and propagate the formal standard errors for surface elevations of the western U.S., incremental changes in surface elevation [Δhek, θk)] for each time step are needed. From Equation 13 in the main text, this can be written as:
graphic
where Δhek, θk) is surface elevation change, ΔH is crustal thickness changes from the kth – 1 to the kth time interval, ρm is the reference upper mantle density, ρc is the crustal density, ρm′(ϕk, θk) is the upper mantle density beneath the western North American plate, and HSL is the thickness of the crust at sea level.
The variance of function Δhek, θk) is
graphic
At present we do not consider errors in our density models for the calculation of error propagation for elevation. We only consider the influence of errors in crustal thickness. Hence, the partial derivative of function Δhek, θk) with respect to crustal thickness changes is:
graphic
By substituting Equation A12 into Equation A11, we have:
graphic
The variance in surface elevation {var[hek, θk)]} at time t = ΔtN involves the sum of variances of incremental changes in surface elevation for each time step, where N is the total number of time steps under consideration:
graphic
Thus the formal standard error of surface elevation {se[hek, θk)]} at the kth time interval is (Fig. 12):
graphic

Calculating Two-Dimensional Steady-State Conductive Thermal Perturbations of Western U.S. Upper Mantle Density

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]).

Using the Laplace equation and assuming constant thermal conductivity, the steady-state conductive heat distribution with no heat generation (Eckert and Drake, 1987) is:
graphic
Equation A16 yields the 2-D temperature (T) as a function of the two independent space coordinates (x = ϕk, y = θk). Using the Fourier equations (Eckert and Drake, 1987), then the heat flow (Q) in the ϕk and θk directions is calculated as:
graphic
graphic
where k is thermal conductivity (W/m·°C), A is the surface area (m2), and Ti is the temperature variation based on four different patterns of magmatism during the past 36 m.y. and the boundary condition. We use 3.5 W/m·°C as thermal conductivity of upper mantle (Robertson, 1988). The heat flow at any point on our 2-D grid area is the sum of the solutions for A17 and A18.

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.

In general, changing either the pressure or the temperature can change the density of the upper mantle. We ignore the pressure dependence for the density of the upper mantle lithosphere, and only consider density changes owing to temperature changes within the upper mantle lithosphere. Hence, we calculate a new time- and temperature-dependent upper mantle density model [ρTik, θk)] for each 0.5 m.y. time interval from 36 Ma to the present day (k = 72−1) based on thermal expansion of the upper mantle at constant pressure and differential temperatures as:
graphic
where ρTk, θk) is the upper mantle density at the kth time interval, α is the thermal expansion coefficient of the upper mantle, and ΔTk–1, θk–1) is the differential temperature between heat distribution solutions for the upper mantle (e.g., ΔT = Ti–1 – Ti, 4 < i < 1). Here we use 4 × 10−5 as the thermal expansion coefficient of the upper mantle (Suzuki, 1975; Robertson, 1988).

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 [ρTk, θk)] for the western U.S. from 36 Ma to the present day.

1Supplemental Animations 1–15. Time-slice animations of western U.S. contoured dilatation strain rates, crustal thickness, paleotopography, and formal standard errors for all calculations along with two-dimensional steady-state conductive heat distribution model of upper mantle and upper mantle density for every 0.5 m.y. from 36 Ma to present. Please visit http://doi.org/10.1130/GES01604.S1 or the full-text article on www.gsapubs.org to view the Supplemental Animations.
2 Supplemental Figures S1–S12. Time-slice figures of western U.S. crustal thickness, upper mantle density, paleotopography, and formal standard errors for all models using the CRUST 1.0 dataset (Laske et al., 2013) for 7 m.y. time intervals from 36 Ma to present. Please visit http://doi.org/10.1130/GES01604.S2 or the full-text article on www.gsapubs.org to view the Supplemental Figures.
3Supplemental Digital Files. Calculated velocity, strain rate, crustal thickness, and paleotopography fields for every 0.5 m.y. from 36 Ma to present. Please visit http://doi.org/10.1130/GES01604.S3 or the full-text article on www.gsapub.org to view the Supplemental Digital Files.
Science Editor: Raymond M. Russo
1.
Afonso
,
J.C.
,
Rawlinson
,
N.
,
Yang
,
Y.
,
Schutt
,
D.L.
,
Jones
,
A.G.
,
Fullea
,
J.
, and
Griffin
,
W.L.
,
2016
,
3-D multiobservable probabilistic inversion for the compositional and thermal structure of the lithosphere and upper mantle: III. Thermochemical tomography in the Western-Central US: Journal of Geophysical Research
:
Solid Earth
 , v.
121
, p.
7337
7370
, https://doi.org/10.1002/2016JB013049.
2.
Armstrong
,
R.L.
,
1982
,
Cordilleran metamorphic core complexes from Arizona to southern Canada
:
Annual Review of Earth and Planetary Sciences
 , v.
10
, p.
129
154
, https://doi.org/10.1146/annurev.ea.10.050182.001021.
3.
Armstrong
,
R.L.
, and
Ward
,
P.
,
1991
,
Evolving geographic patterns of Cenozoic magmatism in the North American Cordillera: The temporal and spatial association of magmatism and metamorphic core complexes
:
Journal of Geophysical Research
 , v.
96
, p.
13,201
13,224
, https://doi.org/10.1029/91JB00412.
4.
Atwater
,
T.
,
1970
,
Implications of plate tectonics for the Cenozoic tectonic evolution of western North America
:
Geological Society of America Bulletin
 , v.
81
, p.
3513
3536
, https://doi.org/10.1130/0016-7606(1970)81[3513:IOPTFT]2.0.CO;2.
5.
Atwater
,
T.
, and
Stock
,
J.
,
1998
,
Pacific–North America plate tectonics of the Neogene southwestern United States
:
An update: International Geology Review
 , v.
40
, p.
375
402
, https://doi.org/10.1080/00206819809465216.
6.
Badgley
,
C.
,
2010
,
Tectonics, topography, and mammalian diversity
:
Ecography
 , v.
33
, p.
220
231
, https://doi.org/10.1111/j.1600-0587.2010.06282.x.
7.
Badgley
,
C.
,
Smiley
,
T.M.
, and
Finarelli
,
J.A.
,
2014
,
Great Basin mammal diversity in relation to landscape history
:
Journal of Mammalogy
 , v.
95
, p.
1090
1106
, https://doi.org/10.1644/13-MAMM-S-088.
8.
Beavan
,
J.
, and
Haines
,
J.
,
2001
,
Contemporary horizontal velocity and strain rate fields of the Pacific-Australian plate boundary zone through New Zealand
:
Journal of Geophysical Research
 , v.
106
, p.
741
770
, https://doi.org/10.1029/2000JB900302.
9.
Becker
,
T.W.
,
Faccenna
,
C.
,
Humphreys
,
E.D.
,
Lowry
,
A.R.
, and
Miller
,
M.S.
,
2014
,
Static and dynamic support of western United States topography
:
Earth and Planetary Science Letters
 , v.
402
, p.
234
246
, https://doi.org/10.1016/j.epsl.2013.10.012.
10.
Best
,
M.G.
,
Christiansen
,
E.H.
, and
Gromme
,
S.
,
2013
,
Introduction: The 36–18 Ma southern Great Basin, USA, ignimbrite province and flareup: Swarms of subduction-related supervolcanoes
:
Geosphere
 , v.
9
, p.
260
274
, https://doi.org/10.1130/GES00870.1.
11.
Best
,
M.G.
,
Christiansen
,
E.H.
,
de Silva
,
S.
, and
Lipman
,
P.W.
,
2016
,
Slab-rollback ignimbrite flareups in the southern Great Basin and other Cenozoic American arcs: A distinct style of arc volcanism
:
Geosphere
 , v.
12
, p.
1097
1135
, https://doi.org/10.1130/GES01285.1.
12.
Bidgoli
,
T.S.
,
Amir
,
E.
,
Walker
,
J.D.
,
Stockli
,
D.F.
,
Andrew
,
J.E.
, and
Caskey
,
S.J.
,
2015
,
Low-temperature thermochronology of the Black and Panamint mountains, Death Valley, California: Implications for geodynamic controls on Cenozoic intraplate strain
:
Lithosphere
 , v.
7
, p.
473
480
, https://doi.org/10.1130/L406.1.
13.
Bird
,
P.
,
1988
,
Formation of the Rocky Mountains, western United States: A continuum computer model
:
Science
 , v.
239
, p.
1501
1507
, https://doi.org/10.1126/science.239.4847.1501.
14.
Bird
,
P.
,
1998
,
Kinematic history of the Laramide orogeny in latitudes 35°–49°N, western United States
:
Tectonics
 , v.
17
, p.
780
801
, https://doi.org/10.1029/98TC02698.
15.
Block
,
L.
, and
Royden
,
L.H.
,
1990
,
Core complex geometries and regional scale flow in the lower crust
:
Tectonics
 , v.
9
, p.
557
567
, https://doi.org/10.1029/TC009i004p00557.
16.
Burchfiel
,
B.C.
,
Cowan
,
D.S.
, and
Davis
,
G.A.
,
1992
,
Tectonic overview of the Cordilleran orogen in the western United States
, in
Burchfiel
,
B.C.
,
Lipman
,
P.W.
, and
Zoback
,
M.L.
, eds.,
The Cordilleran Orogen: Conterminous U.S.: Boulder, Colorado, Geological Society of America, The Geology of North America
 , v.
G-3
, p.
407
479
, https://doi.org/10.1130/DNAG-GNA-G3.407.
17.
Cassel
,
E.J.
,
Graham
,
S.A.
, and
Chamberlain
,
C.P.
,
2009
,
Cenozoic tectonic and topographic evolution of the northern Sierra Nevada, California, through stable isotope paleoaltimetry in volcanic glass
:
Geology
 , v.
37
, p.
547
550
, https://doi.org/10.1130/G25572A.1.
18.
Cassel
,
E.J.
,
Breecker
,
D.O.
,
Henry
,
C.D.
,
Larson
,
T.E.
, and
Stockli
,
D.F.
,
2014
,
Profile of a paleo-orogen: High topography across the present-day Basin and Range from 40 to 23 Ma
:
Geology
 , v.
42
, p.
1007
1010
, https://doi.org/10.1130/G35924.1.
19.
Cecil
,
M.R.
,
Ducea
,
M.N.
,
Reiners
,
P.
,
Gehrels
,
G.
,
Mulch
,
A.
,
Allen
,
C.
, and
Campbell
,
I.
,
2010
,
Provenance of Eocene river sediments from the central northern Sierra Nevada and implications for paleotopography
:
Tectonics
 , v.
29
,
TC6010
, https://doi.org/10.1029/2010TC002717.
20.
Chamberlain
,
C.P.
, and
Poage
,
M.A.
,
2000
,
Reconstructing the paleotopography of mountain belts from the isotopic composition of authigenic minerals
:
Geology
 , v.
28
, p.
115
118
, https://doi.org/10.1130/0091-7613(2000)28<115:RTPOMB>2.0.CO;2.
21.
Chapman
,
J.B.
,
Ducea
,
M.N.
,
DeCelles
,
P.G.
, and
Profeta
,
L.
,
2015
,
Tracking changes in crustal thickness during orogenic evolution with Sr/Y: An example from the North American Cordillera
:
Geology
 , v.
43
, p.
919
922
, https://doi.org/10.1130/G36996.1.
22.
Coney
,
P.J.
,
1980
,
Cordilleran metamorphic core complexes: An overview
, in
Crittenden
,
M.D.
, Jr
,
Coney
,
P.J.
, and
Davis
,
G.H.
, eds.,
Cordilleran Metamorphic Core Complexes: Geological Society of America Memoir 153
 , p.
7
31
, https://doi.org/10.1130/MEM153-p7.
23.
Coney
,
P.J.
, and
Harms
,
T.A.
,
1984
,
Cordilleran metamorphic core complexes: Cenozoic extensional relics of Mesozoic compression
:
Geology
 , v.
12
, p.
550
554
, https://doi.org/10.1130/0091-7613(1984)12<550:CMCCCE>2.0.CO;2.
24.
Cooley
,
M.E.
, and
Davidson
,
E.S.
,
1963
,
The Mogollon Highlands; their influence on Mesozoic and Cenozoic erosion and sedimentation
:
Arizona Geological Society Digest
 , v.
6
, p.
7
35
.
25.
Davis
,
G.A.
,
1988
,
Rapid upward transport of mid-crustal mylonitic gneisses in the footwall of a Miocene detachment fault, Whipple Mountains, southeastern California
:
Geologische Rundschau
 , v.
77
, p.
191
209
, https://doi.org/10.1007/BF01848684.
26.
Davis
,
G.H.
,
1980
,
Structural characteristics of metamorphic core complexes, southern Arizona
, in
Crittenden
,
M.D.
, Jr
,
Coney
,
P.J.
, and
Davis
,
G.H.
, eds.,
Cordilleran Metamorphic Core Complexes: Geological Society of America Memoir 153
 , p.
35
78
, https://doi.org/10.1130/MEM153-p35.
27.
DeCelles
,
P.G.
,
2004
,
Late Jurassic to Eocene evolution of the Cordilleran thrust belt and foreland basin system, western USA
:
American Journal of Science
 , v.
304
, p.
105
168
, https://doi.org/10.2475/ajs.304.2.105.
28.
DeCelles
,
P.G.
,
Lawton
,
T.F.
, and
Mitra
,
G.
,
1995
,
Thrust timing, growth of structural culminations, and synorogenic sedimentation in the type Sevier orogenic belt, western United States
:
Geology
 , v.
23
, p.
699
702
, https://doi.org/10.1130/0091-7613(1995)023<0699:TTGOSC>2.3.CO;2.
29.
DeMets
,
C.
, and
Dixon
,
T.H.
,
1999
,
New kinematic models for Pacific–North America motion from 3 Ma to present, I: Evidence for steady motion and biases in the NUVEL-1A model
:
Geophysical Research Letters
 , v.
26
, p.
1921
1924
, https://doi.org/10.1029/1999GL900405.
30.
Dickinson
,
W.R.
,
2002
,
The Basin and Range province as a composite extensional domain
:
International Geology Review
 , v.
44
, p.
1
38
, https://doi.org/10.2747/0020-6814.44.1.1.
31.
Dickinson
,
W.R.
,
2003
,
The Basin and Range Province as a composite extensional domain
, in
Klemperer
,
S.L.
, and
Ernst
W.G.
, eds.,
The Lithosphere of Western North America and Its Geophysical Characterization: The George A. Thompson Volume: Geological Society of America International Book Series
 , v.
7
, p.
213
250
.
32.
Dickinson
,
W.R.
,
2004
,
Evolution of the North American Cordillera
:
Annual Review of Earth and Planetary Sciences
 , v.
32
, p.
13
45
, https://doi.org/10.1146/annurev.earth.32.101802.120257.
33.
Dickinson
,
W.R.
,
2006
,
Geotectonic evolution of the Great Basin
:
Geosphere
 , v.
2
, p.
353
368
, https://doi.org/10.1130/GES00054.1.
34.
Eckert
,
E.R.G.
and
Drake
,
R.M.
, Jr
.,
1987
,
Analysis of Heat and Mass Transfer
:
New York
,
Hemisphere Publishing
,
806
p.
35.
Elston
,
D.P.
, and
Young
,
R.A.
,
1991
,
Cretaceous-Eocene (Laramide) landscape development and Oligocene-Pliocene drainage reorganization of transition zone and Colorado Plateau, Arizona: Journal of Geophysical Research
:
Solid Earth
 , v.
96
, p.
12389
12406
.
36.
Faulds
,
J.E.
,
Henry
,
C.D.
, and
Hinz
,
N.H.
,
2005
,
Kinematics of the northern Walker Lane: An incipient transform fault along the Pacific–North American plate boundary
:
Geology
 , v.
33
, p.
505
508
, https://doi.org/10.1130/G21274.1.
37.
Flesch
,
L.M.
,
Holt
,
W.E.
,
Haines
,
A.J.
, and
Shen-Tu
,
B.
,
2000
,
Dynamics of the Pacific–North American plate boundary in the western United States
:
Science
 , v.
287
, p.
834
836
, https://doi.org/10.1126/science.287.5454.834.
38.
Flesch
,
L.M.
,
Holt
,
W.E.
,
Haines
,
A.J.
,
Wen
,
L.
, and
Shen-Tu
,
B.
,
2007
,
The dynamics of western North America: Stress magnitudes and the relative role of gravitational potential energy, plate interaction at the boundary and basal tractions
:
Geophysical Journal International
 , v.
169
, p.
866
896
, https://doi.org/10.1111/j.1365-246X.2007.03274.x.
39.
Flowers
,
R.M.
,
Wernicke
,
B.P.
, and
Farley
,
K.A.
,
2008
,
Unroofing, incision, and uplift history of the southwestern Colorado Plateau from apatite (U-Th)/He thermochronometry
:
Geological Society of America Bulletin
 
120
, p.
571
587
, https://doi.org/10.1130/B26231.1.
40.
Garside
,
L.J.
,
Henry
,
C.D.
,
Faulds
,
J.E.
, and
Hinz
,
N.H.
,
2005
,
The upper reaches of the Sierra Nevada auriferous gold channels, California and Nevada
, in
Rhoden
,
H.N.
,
Steininger
,
R.C.
, and
Vikre
,
P.G.
, eds.,
Geological Society of Nevada Symposium 2005: Window to the World, Reno, Nevada
 , May 2005, p.
209
235
.
41.
Gébelin
,
A.
,
Mulch
,
A.
,
Teyssier
,
C.
,
Chamberlain
,
C.P.
, and
Heizler
,
M.
,
2012
,
Coupled basin-detachment systems as paleoaltimetry archives of the western North American Cordillera
:
Earth and Planetary Science Letters
 , v.
335
, p.
36
47
, https://doi.org/10.1016/j.epsl.2012.04.029.
42.
Ghosh
,
A.
, and
Holt
,
W.E.
,
2012
,
Plate motions and stresses from global dynamic models
:
Science
 , v.
335
, p.
838
843
, https://doi.org/10.1126/science.1214209.
43.
Ghosh
,
A.
,
Becker
,
T.W.
, and
Humphreys
,
E.D.
,
2013a
,
Dynamics of the North American continent
:
Geophysical Journal International
 , v.
194
, p.
651
669
, https://doi.org/10.1093/gji/ggt151.
44.
Ghosh
,
A.
,
Holt
,
W.E.
, and
Wen
,
L.
,
2013b
,
Predicting the lithospheric stress field and plate motions by joint modeling of lithosphere and mantle dynamics
:
Journal of Geophysical Research: Solid Earth
 , v.
118
, p.
346
368
, https://doi.org/10.1029/2012JB009516.
45.
Gregory-Wodzicki
,
K.M.
,
1997
,
The late Eocene house range flora, Sevier Desert, Utah: Paleoclimate and paleoelevation
:
Palaios
 
12
, p.
552
567
, https://doi.org/10.2307/3515411.
46.
Haines
,
A.J.
, and
Holt
,
W.E.
,
1993
,
A procedure for obtaining the complete horizontal motions within zones of distributed deformation from the inversion of strain rate data
:
Journal of Geophysical Research
 , v.
98
, p.
12,057
12,082
, https://doi.org/10.1029/93JB00892.
47.
Hamilton
,
W.
,
1969
,
Mesozoic California and the underflow of Pacific mantle
:
Geological Society of America Bulletin
 , v.
80
, p.
2409
2430
, https://doi.org/10.1130/0016-7606(1969)80[2409:MCATUO]2.0.CO;2.
48.
Hamilton
,
W.
, and
Myers
,
W.B.
,
1966
,
Cenozoic tectonics of the western United States
:
Reviews of Geophysics
 , v.
4
, p.
509
549
, https://doi.org/10.1029/RG004i004p00509.
49.
Henry
,
C.D.
,
2009
,
Uplift of the Sierra Nevada, California
:
Geology
 , v.
37
, p.
575
576
, https://doi.org/10.1130/focus062009.1.
50.
Holt
,
W.E.
,
Chamot-Rooke
,
N.
,
Le Pichon
,
X.
,
Haines
,
A.J.
,
Shen-Tu
,
B.
, and
Ren
,
J.
,
2000
,
Velocity field in Asia inferred from Quaternary fault slip rates and Global Positioning System observations
:
Journal of Geophysical Research
 , v.
105
, p.
19,185
19,209
, https://doi.org/10.1029/2000JB900045.
51.
Horton
,
T.W.
,
Sjostrom
,
D.J.
,
Abruzzese
,
M.J.
,
Poage
,
M.A.
,
Waldbauer
,
J.R.
,
Hren
,
M.
,
Wooden
,
J.
, and
Chamberlain
,
C.P.
,
2004
,
Spatial and temporal variation of Cenozoic surface elevation in the Great Basin and Sierra Nevada
:
American Journal of Science
 , v.
304
, p.
862
888
, https://doi.org/10.2475/ajs.304.10.862.
52.
Humphreys
,
E.D.
,
1995
,
Post-Laramide removal of the Farallon slab, western United States
:
Geology
 , v.
23
, p.
987
990
, https://doi.org/10.1130/0091-7613(1995)023<0987:PLROTF>2.3.CO;2.
53.
Humphreys
,
E.
,
2009
,
Relation of flat subduction to magmatism and deformation in the western United States
, in
Kay
,
S.M.
,
Ramos
,
V.A.
, and
Dickinson
,
W.R.
, eds.,
Backbone of the Americas: Shallow Subduction, Plateau Uplift, and Ridge and Terrane Collision: Geological Society of America Memoir 204
 , p.
85
98
, https://doi.org/10.1130/2009.1204(04).
54.
Humphreys
,
E.D.
, and
Coblentz
,
D.D.
,
2007
,
North American dynamics and western U.S. tectonics
:
Reviews of Geophysics
 , v.
45
, RG3001, https://doi.org/10.1029/2005RG000181.
55.
Humphreys
,
E.
,
Hessler
,
E.
,
Dueker
,
K.
,
Farmer
,
G.L.
,
Erslev
,
E.
, and
Atwater
,
T.
,
2003
,
How Laramide-age hydration of North American lithosphere by the Farallon slab controlled subsequent activity in the western United States
:
International Geology Review
 , v.
45
, p.
575
595
, https://doi.org/10.2747/0020-6814.45.7.575.
56.
Huntington
,
K.W.
,
Wernicke
,
B.P.
, and
Eiler
,
J.M.
,
2010
,
Influence of climate change and uplift on Colorado Plateau paleotemperatures from carbonate clumped isotope thermometry
:
Tectonics
 , v.
29
,
TC3005
, https://doi.org/10.1029/2009TC002449.
57.
Jones
,
C.H.
,
Unruh
,
J.R.
, and
Sonder
,
L.J.
,
1996
,
The role of gravitational potential energy in active deformation in the southwestern United States
:
Nature
 , v.
381
, p.
37
41
, https://doi.org/10.1038/381037a0.
58.
Jones
,
C.H.
,
Mahan
,
K.H.
,
Butcher
,
L.A.
,
Levandowski
,
W.B.
, and
Farmer
,
G.L.
,
2015
,
Continental uplift through crustal hydration
:
Geology
 , v.
43
, p.
355
358
, https://doi.org/10.1130/G36509.1.
59.
Kent-Corson
,
M.L.
,
Sherman
,
L.S.
,
Mulch
,
A.
, and
Chamberlain
,
C.P.
,
2006
,
Cenozoic topographic and climatic response to changing tectonic boundary conditions in western North America
:
Earth and Planetary Science Letters
 , v.
252
, p.
453
466
, https://doi.org/10.1016/j.epsl.2006.09.049.
60.
Laske
,
G.
,
Ma
,
Z.
,
Masters
,
G.
, and
Pasyanos
,
M.
,
2013
,
CRUST 1.0
: http://igppweb.ucsd.edu/∼gabi/crust1.html (last accessed August 2013).
61.
Lechler
,
A.R.
,
Niemi
,
N.A.
,
Hren
,
M.T.
, and
Lohmann
,
K.C.
,
2013
,
Paleoelevation estimates for the northern and central proto–Basin and Range from carbonate clumped isotope thermometry
:
Tectonics
 , v.
32
, p.
295
316
, https://doi.org/10.1002/tect.20016.
62.
Levandowski
,
W.
,
Jones
,
C.H.
,
Shen
,
W.
,
Ritzwoller
,
M.H.
, and
Schulte-Pelkum
,
V.
,
2014
,
Origins of topography in the western U.S.: Mapping crustal and upper mantle density variations using a uniform seismic velocity model: Journal of Geophysical Research
:
Solid Earth
 , v.
119
, p.
2375
2396
, https://doi.org/10.1002/2013JB010607.
63.
Liu
,
L.
, and
Gurnis
,
M.
,
2010
,
Dynamic subsidence and uplift of the Colorado Plateau
:
Geology
 , v.
38
, p.
663
666
, https://doi.org/10.1130/G30624.1.
64.
Liu
,
L.
,
Spasojevic´
,
S.
, and
Gurnis
,
M.
,
2008
,
Reconstructing Farallon plate subduction beneath North America back to the Late Cretaceous
:
Science
 , v.
322
, p.
934
938
, https://doi.org/10.1126/science.1162921.
65.
Liu
,
L.
,
Gurnis
,
M.
,
Seton
,
M.
,
Saleeby
,
J.
,
Müller
,
R.D.
, and
Jackson
,
J.M.
,
2010
,
The role of oceanic plateau subduction in the Laramide orogeny
:
Nature Geoscience
 , v.
3
, p.
353
357
, https://doi.org/10.1038/ngeo829.
66.
Liu
,
M.
,
2001
,
Cenozoic extension and magmatism in the North American Cordillera: The role of gravitational collapse
:
Tectonophysics
 , v.
342
, p.
407
433
, https://doi.org/10.1016/S0040-1951(01)00173-1.
67.
Liu
,
M.
, and
Shen
,
Y.
,
1998
,
Crustal collapse, mantle upwelling, and Cenozoic extension in the North American Cordillera
:
Tectonics
 , v.
17
, p.
311
321
, https://doi.org/10.1029/98TC00313.
68.
Lowry
,
A.R.
, and
Pérez-Gussinyé
,
M.
,
2011
,
The role of crustal quartz in controlling Cordilleran deformation
:
Nature
 , v.
471
, p.
353
357
, https://doi.org/10.1038/nature09912.
69.
McKenzie
,
D.
, and
Jackson
,
J.
,
1983
,
The relationship between strain rates, crustal thickening, palaeomagnetism, finite strain and fault movements within a deforming zone
:
Earth and Planetary Science Letters
 , v.
65
, p.
182
202
, https://doi.org/10.1016/0012-821X(83)90198-X.
70.
McQuarrie
,
N.
, and
Wernicke
,
B.P.
,
2005
,
An animated tectonic reconstruction of southwestern North America since 36 Ma
:
Geosphere
 , v.
1
, p.
147
172
, https://doi.org/10.1130/GES00016.1.
71.
Mix
,
H.T.
,
Mulch
,
A.
,
Kent-Corson
,
M.L.
, and
Chamberlain
,
C.P.
,
2011
,
Cenozoic migration of topography in the North American Cordillera
:
Geology
 , v.
39
, p.
87
90
, https://doi.org/10.1130/G31450.1.
72.
Mooney
,
W.D.
, and
Kaban
,
M.K.
,
2010
,
The North American upper mantle: Density, composition, and evolution
:
Journal of Geophysical Research
 , v.
115
,
B12424
, https://doi.org/10.1029/2010JB000866.
73.
Moschetti
,
M.P.
,
Ritzwoller
,
M.H.
,
Lin
,
F.-C.
, and
Yang
,
Y.
,
2010
,
Crustal shear wave velocity structure of the western United States inferred from ambient seismic noise and earthquake data
:
Journal of Geophysical Research
 , v.
115
,
B10306
, https://doi.org/10.1029/2010JB007448.
74.
Moucha
,
R.
,
Forte
,
A.M.
,
Rowley
,
D.B.
,
Mitrovica
,
J.X.
,
Simmons
,
N.A.
, and
Grand
,
S.P.
,
2008
,
Mantle convection and the recent evolution of the Colorado Plateau and the Rio Grande Rift valley
:
Geology
 , v.
36
, p.
439
442
, https://doi.org/10.1130/G24577A.1.
75.
Obrebski
,
M.
,
Allen
,
R.M.
,
Pollitz
,
F.
, and
Hung
,
S.-H.
,
2011
,
Lithosphere-asthenosphere interaction beneath the western United States from the joint inversion of body-wave traveltimes and surface-wave phase velocities
:
Geophysical Journal International
 , v.
185
, p.
1003
1021
, https://doi.org/10.1111/j.1365-246X.2011.04990.x.
76.
Parrish
,
R.R.
,
Carr
,
S.D.
, and
Parkinson
,
D.L.
,
1988
,
Eocene extensional tectonics and geochronology of the southern Omineca Belt, British Columbia and Washington
:
Tectonics
 , v.
7
, p.
181
212
, https://doi.org/10.1029/TC007i002p00181.
77.
Porter
,
R.
,
Liu
,
Y.Y.
, and
Holt
,
W.E.
,
2016
,
Lithospheric records of orogeny within the continental U.S
.:
Geophysical Research Letters
 , v.
43
, p.
144
153
, https://doi.org/10.1002/2015GL066950.
78.
Porter
,
R.
,
Hoisch
,
T.
, and
Holt
,
W.E.
,
2017
,
The role of lower-crustal hydration in the tectonic evolution of the Colorado Plateau
:
Tectonophysics
 , v.
712–713
, p.
221
231
, https://doi.org/10.1016/j.tecto.2017.05.025.
79.
Robertson
,
E.C.
,
1988
,
Thermal properties of rocks
:
U.S. Geological Survey Open-File Report
 
88-441
,
106
p.
80.
Roy
,
M.
,
Jordan
,
T.H.
,
Pederson
,
J.
,
2009
,
Colorado Plateau magmatism and uplift by warming of heterogeneous lithosphere
:
Nature
 , v.
459
, p.
978
982
, https://doi.org/10.1038/nature08052.
81.
Saleeby
,
J.B.
,
1983
,
Accretionary tectonics of the North American Cordillera
:
Annual Review of Earth and Planetary Sciences
 , v.
11
, p.
45
73
, https://doi.org/10.1146/annurev.ea.11.050183.000401.
82.
Saleeby
,
J.
,
Ducea
,
M.
, and
Clemens-Knott
,
D.
,
2003
,
Production and loss of high-density batholithic root, southern Sierra Nevada, California
:
Tectonics
 , v.
22
,
1064
, https://doi.org/10.1029/2002TC001374.
83.
Schmandt
,
B.
, and
Humphreys
,
E.
,
2010
,
Complex subduction and small-scale convection revealed by body-wave tomography of the western United States upper mantle
:
Earth and Planetary Science Letters
 , v.
297
, p.
435
445
, https://doi.org/10.1016/j.epsl.2010.06.047.
84.
Schmandt
,
B.
,
Lin
,
F.C.
, and
Karlstrom
,
K.E.
,
2015
,
Distinct crustal isostasy trends east and west of the Rocky Mountain Front
:
Geophysical Research Letters
 , v.
42
, p.
10,290
10,298
, https://doi.org/10.1002/2015GL066593.
85.
Schutt
,
D.
,
Lowry
,
A.R.
, and
Buehler
,
J.S.
,
2012
,
The temperature of the western U.S. lithosphere
:
Abstract V53A-2812 presented at 2012 Fall Meeting, American Geophysical Union
,
San Francisco, California
,
3
7 Dec
.
86.
Shen
,
W.S.
, and
Ritzwoller
,
M.H.
,
2016
,
Crustal and uppermost mantle structure beneath the United States: Journal of Geophysical Research
:
Solid Earth
 , v.
121
, p.
4306
4342
, https://doi.org/10.1002/2016JB012887.
87.
Sigloch
,
K.
, and
Mihalynuk
,
M.G.
,
2013
,
Intra-oceanic subduction shaped the assembly of Cordilleran North America
:
Nature
 , v.
496
, p.
50
56
, https://doi.org/10.1038/nature12019.
88.
Snedecor
,
G.W.
, and
Cochran
,
W.G.
,
1994
,
Statistical Methods (first East West Press edition)
:
New Delhi
,
Affiliated East West Private Ltd
.
89.
Snow
,
J.K.
, and
Wernicke
,
B.P.
,
2000
,
Cenozoic tectonism in the central Basin and Range: Magnitude, rate, and distribution of upper crustal strain
:
American Journal of Science
 , v.
300
, p.
659
719
, https://doi.org/10.2475/ajs.300.9.659.
90.
Sonder
,
L.J.
, and
Jones
,
C.H.
,
1999
,
Western United States extension: How the west was widened
:
Annual Review of Earth and Planetary Sciences
 , v.
27
, p.
417
462
, https://doi.org/10.1146/annurev.earth.27.1.417.
91.
Sonder
,
L.J.
,
England
,
P.C.
,
Wernicke
,
B.P.
, and
Christiansen
,
R.L.
,
1987
,
A physical model for Cenozoic extension of western North America
, in
Coward
,
M.P.
,
Dewey
,
J.F.
, and
Hancock
,
P.L.
, eds.,
Continental Extensional Tectonics: Geological Society of London Special Publication 28
 , p.
187
201
, https://doi.org/10.1144/GSL.SP.1987.028.01.14.
92.
Spencer
,
J.E.
,
1984
,
Role of tectonic denudation in warping and uplift of low-angle normal faults
:
Geology
 , v.
12
, p.
95
98
, https://doi.org/10.1130/0091-7613(1984)12<95:ROTDIW>2.0.CO;2.
93.
Spencer
,
J.E.
,
Smith
,
G.R.
, and
Dowling
,
T.E.
,
2008
,
Middle to late Cenozoic geology, hydrography, and fish evolution in the American Southwest
, in
Reheis
,
M.C.
,
Hershler
,
R.
, and
Miller
,
D.M.
, eds.,
Late Cenozoic Drainage History of the Southwestern Great Basin and Lower Colorado River Region: Geologic and Biotic Perspectives: Geological Society of America Special Paper 439
 , p.
279
299
, https://doi.org/10.1130/2008.2439(12).
94.
Suzuki
,
I.
,
1975
,
Thermal expansion of periclase and olivine, and their anharmonic properties
:
Journal of Physics of the Earth
 , v.
23
, p.
145
159
, https://doi.org/10.4294/jpe1952.23.145.
95.
Tipler
,
P.A.
, and
Mosca
,
G.
,
2007
,
Physics for Scientists and Engineers (sixth edition)
:
New York
,
Macmillan
,
1116
p.
96.
Wernicke
,
B.
,
1992
,
Cenozoic extensional tectonics of the U.S. Cordillera
, in
Burchfiel
,
B.C.
,
Lipman
,
P.W.
, and
Zoback
,
M.L.
, eds.,
The Cordilleran Orogen: Conterminous U.S
 .:
Boulder, Colorado
,
Geological Society of America, The Geology of North America
, v.
G-3
, p.
553
581
, https://doi.org/10.1130/DNAG-GNA-G3.553.
97.
Wernicke
,
B.
,
2011
,
The California River and its role in carving Grand Canyon
:
Geological Society of America Bulletin
 , v.
123
, p.
1288
1316
, https://doi.org/10.1130/B30274.1.
98.
Wernicke
,
B.
, and
Burchfiel
,
B.C.
,
1982
,
Modes of extensional tectonics
:
Journal of Structural Geology
 , v.
4
, p.
105
115
, https://doi.org/10.1016/0191-8141(82)90021-9.
99.
Wernicke
,
B.P.
,
England
,
P.C.
,
Sonder
,
L.J.
, and
Christiansen
,
R.L.
,
1987
,
Tectonomagmatic evolution of Cenozoic extension in the North American Cordillera
, in
Coward
,
M.P.
,
Dewey
,
J.F.
, and
Hancock
,
P.L.
, eds.,
Continental Extensional Tectonics: Geological Society of London Special Publication 28
 , p.
203
221
, https://doi.org/10.1144/GSL.SP.1987.028.01.15.
100.
Wernicke
,
B.
,
Axen
,
G.J.
, and
Snow
,
J.K.
,
1988
,
Basin and Range extensional tectonics at the latitude of Las Vegas, Nevada
:
Geological Society of America Bulletin
 , v.
100
, p.
1738
1757
, https://doi.org/10.1130/0016-7606(1988)100<1738:BARETA>2.3.CO;2.
101.
Wessel
,
P.
,
Smith
,
W.H.
,
Scharroo
,
R.
,
Luis
,
J.
, and
Wobbe
,
F.
,
2013
,
Generic mapping tools: Improved version released
:
Eos (Transactions, American Geophysical Union)
 , v.
94
, p.
409
410
, https://doi.org/10.1002/2013EO450001.
102.
Wolfe
,
J.A.
,
Schorn
,
H.E.
,
Forest
,
C.E.
, and
Molnar
,
P.
,
1997
,
Paleobotanical evidence for high altitudes in Nevada during the Miocene
:
Science
 , v.
276
, p.
1672
1675
, https://doi.org/10.1126/science.276.5319.1672.
103.
Wolfe
,
J.A.
,
Forest
,
C.E.
, and
Molnar
,
P.
,
1998
,
Paleobotanical evidence of Eocene and Oligocene paleoaltitudes in midlatitude western North America
:
Geological Society of America Bulletin
 , v.
110
, p.
664
678
, https://doi.org/10.1130/0016-7606(1998)110<0664:PEOEAO>2.3.CO;2.
104.
Yang
,
Y.
,
Ritzwoller
,
M.H.
,
Lin
,
F.-C.
,
Moschetti
,
M.P.
, and
Shapiro
,
N.M.
,
2008
,
Structure of the crust and uppermost mantle beneath the western United States revealed by ambient noise and earthquake tomography
:
Journal of Geophysical Research
 , v.
113
,
B12310
, https://doi.org/10.1029/2008JB005833.
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.