Abstract

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.

INTRODUCTION

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

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

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

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

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

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

METHODS AND RESULTS

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.

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

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.

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

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

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

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

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

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

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

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

DISCUSSION

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

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.

CONCLUSIONS

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.

ACKNOWLEDGMENTS

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
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.