Seismological estimates of means of isostatic support of the Sierra Nevada Open Access

The crust beneath the Sierra Nevada (California, USA; Fig. 1) is thinner (Fliedner et al., 2000; Frassetto et al., 2011; Jones et al., 1994; Ruppert et al., 1998) and more felsic than the worldwide average (e.g., Holbrook et al., 1992), and substantially thinner under the high parts of the range than under the low western foothills (Fig. 2B). Despite this thin crust, the mean surface elevation within the range is ∼2.5 km (Fig. 2A). Gravity measurements show the range to be in isostatic equilibrium or even somewhat overcompensated, i.e., having a greater gravity low than expected if the observed topography were compensated by an Airy root (Kennelly and Chase, 1989; Simpson et al., 1986). Such elevations in the absence of correlated variations in crustal thickness require lateral density variations within the crust and/or uppermost mantle (e.g., Jones et al., 1994).

The unusual crustal variations of the Sierra are matched by controversy over the origin of the topography. From the time of work by LeConte (1886) to more recent times, geologists inferred that the Sierra Nevada rose during the late Tertiary (e.g., Christensen, 1966; Huber, 1981; Unruh, 1991; Wakabayashi and Sawyer, 2001). This assertion was challenged first by U/Th-He dating of erosion (House et al., 1998, 2001) and later by isotopic measurements inferred to reflect isotopic fractionation of precipitation associated with topography (Cassel et al., 2009, 2012; Mulch et al., 2006, 2008; Poage and Chamberlain, 2002). Although this interpretation remains controversial (Galewsky, 2009a, 2009b; Jones et al., 2004; Molnar, 2010), inferences that the Sierran crest lost an antibuoyant root between 10 and 3 Ma (Ducea and Saleeby, 1996, 1998; Farmer et al., 2002; Manley et al., 2000) suggested a “Sierran paradox” of an old range with a young means of support (Wernicke et al., 1996). Quantifying the contributions of compositional and thermal variations of crust and mantle to the modern support of the range can help to resolve this controversy and highlight areas where a dynamic component (here defined as one that arises from the convective regime rather than directly from density variations) is necessary.

The 2005–2007 Sierra Nevada EarthScope Project (SNEP) deployment of FlexArray broadband seismometers across much of the Sierra Nevada has led to a much fuller three-dimensional characterization of the seismic structure of the range (Frassetto et al., 2011; Gilbert et al., 2012; Reeg, 2008). With such data we quantify the present-day distribution of density consistent with seismic and topographic constraints of the region. The simple assumptions that seismic wave speeds vary because of temperature variations in the mantle and because of a combination of compositional and temperature variations in the crust adequately recover topography. We thus do not find any need for dynamic topography in the region. Furthermore, the derived density variations can then be interpreted to limit the range of plausible means of supporting the topography and thus the tectonic history of the region.

H₀ is a correction term of 2.4 km to achieve isostatic equilibrium with an asthenospheric column (via mid-ocean ridges). The depth of the Moho as defined by receiver functions is z_c (Frassetto et al., 2011). In this study, we examine to a depth, z_a, of 195 km based on the loss of correlation between seismic velocity and topography below this depth. The density of the asthenosphere is ρ_a and is assumed to be 3250 kg/m³. Because the motivation of this study is to explore the source of heterogeneity in the region, the exact choice of reference parameters (including asthenospheric density and crustal velocity) is of second-order importance. The value ε is the elevation above sea level with elastic flexural effects removed (Fig. 2A) by smoothing by convolution with a zero-order Bessel function (Watts, 2001; discussed in the following).

In order to calculate the topography generated by the crust and mantle, we convert the P-wave speeds from the seismic tomography (Reeg, 2008; Jones et al., 2014) into a density profile for each column, with a division at the Moho. Separating the support for smoothed topography into crustal and mantle components is necessary because we use different approaches in the crust and mantle to derive densities from seismic wave speeds. Seismic data from SNEP are used to constrain crustal thickness (using receiver functions; Frassetto et al., 2011) and wave speed variations in the crust and upper mantle (from P-wave tomography; Jones et al., 2014).

Jones et al. (2014) present seismic velocity models derived from inversion of traveltime residuals subject to an array of starting models. Because teleseismic P-wave tomography has notoriously poor vertical resolution (especially in the crust), we choose a tomography model that was derived from an inversion where the top 55 km of the starting model was defined by ambient noise and earthquake surface wave tomography (Moschetti et al., 2010). The top of the starting model was held fixed for 14 of 21 iterations of the inversion, then allowed to vary as necessary to best fit teleseismic P-wave delay times. Because of the greater depth sensitivity of surface waves (at the cost of lateral resolution), this inversion scheme provides a more robust view of crustal velocity variations than would be available from teleseismic P-wave tomography alone. We have also completed our analysis with each of the other inversion schemes described by Jones et al.(2014) from a variety of starting models. Our results are generally robust with respect to seismic models, excepting the statistically significant difference arising in west-central Nevada, outside of the Sierra Nevada region with which we are primarily concerned.

Because we use different relationships between velocity and density (discussed in the following) in the crust and the mantle, we divide the tomography at the Moho depths determined in Frassetto et al. (2011). Nevertheless, the tomography of Jones et al. (2014) was calculated assuming a flat Moho. Thus, while we consider crustal thickness variations in estimated densities and crust- and/or mantle-derived topography, we must also adjust the teleseismic tomography for biases introduced by variable crustal thickness. This correction and its second-order effect on predicted topography are described in  Appendix 1.

We calculate the mantle-derived topography by assuming that density and wave speed variation is due to thermal heterogeneity. Isobaric heating will produce a decrease in both density and seismic velocity. Over a wide variety of lherzolite, harzburgite, and peridotite mineralogies the temperature derivatives of density are nearly equal, although the absolute densities vary considerably (Hacker and Abers, 2004). Therefore, we make no assumption of mantle mineralogy other than that it is laterally constant across the study area at any depth. Using a compositionally independent conversion of velocity anomalies to density anomalies, we can then constrain the mantle contributions to isostasy for a purely thermally varying mantle, interpreting P-wave speed variations reported by Jones et al. (2014) as temperature variations and calculating the resulting density structure.

Laboratory data (Jackson and Faul, 2010) show a nonlinear dependence of shear modulus on temperature, particularly within 150–250 °C of the solidus. For this reason, we employ a decay in the density response to decreasing velocity for negative velocity perturbations, assuming that mean velocities in the upper mantle reflect temperatures of ∼1050–1200 °C. This assumption is justified by estimates of minimum temperatures. At low temperatures velocity decreases <1% per 100 °C (Hacker and Abers, 2004). Thus, a region with perturbation 5.5% (the maximum observed at 120 km) has a temperature at least 600 °C colder than background. This would correspond to a temperature (T) of <600 °C, which we take as a reasonable lower limit at this depth.

Between 0% and –3% velocity anomaly, our estimate of δρ/δVp decreases from 10 to 6 kg/m³ per 1% velocity anomaly. This choice of parameters simulates the behavior of a variety of lherzolite, peridotite, and harzburgite mineralogies in the absence of anelasticity (i.e., at low temperature) (Hacker and Abers, 2004) and elastic response similar to that of millimeter-scale single crystal grains of olivine (Jackson and Faul, 2010). A schematic of the curvilinear relationship between velocity and density in the mantle is presented in Figure 3.

We also consider another end-member regression, a linear relationship between mantle velocity and density with no attempt to account for possible anelastic effects, but we find that our interpretations are relatively insensitive to the choice of regression (see discussion of results).

We note that P-wave speed variations caused by melt (–3.6% per 1% in situ melt fraction) (e.g., Hammond and Humphreys, 2000) produce disproportionately small changes in bulk density (between 0 and 4 kg/m³ per 1% in situ melt fraction) that will cause the magnitude of variation in H_m to be too great. Because melt production and ponding depends on many parameters, we do not explicitly incorporate this effect in our wave speed–density relationship, but consider this effect later.

With this mantle topography (Fig. 4A) we examine the thermal and chemical variations in the crust.

Seismic anomalies in the crust likely reflect variations in composition or temperature. We first convert P-wave speeds to density within the crust using regressions of density onto P-wave speed. We use velocities and densities reported by Christensen and Mooney (1995), correct the room-temperature densities reported for the temperature at which velocities are reported (using the coefficients of thermal expansion from the same compilation and using their average, ∼15 °C/km geotherm), and regress density onto velocity for all nonvolcanic and polymineralic rocks and for dunite and pyroxenite. From 0 to 10 km, ρ = 827.5 + 316.6 × Vp. From 10 to 30 km, ρ = 643 + 337 × Vp; from 30 km to the Moho, ρ = 552.0 + 350.5 × Vp. Because the velocity structure reported by Jones et al. (2014) is a perturbation from mean wave speeds, we also choose a reference velocity structure of 5.8 km/s from the surface to 10 km, 6.5 km/s from 10 to 30 km depth, and 6.8 km/s from 30 km to the Moho based on a tomographic model determined using regional seismicity (Thurber et al., 2009).

The initial assumption of an isothermal crust maximizes crustal density variations. For example, a 3% increase in velocity from our reference at 20 km (Vp = 6.5 km/s, ρ = 2879 kg/m³) due to compositional variations predicts a ∼65 kg/m³ higher density.

We employ the thermal-compositional basis in velocity-density space to generate a thermal structure that reproduces both velocities and topography.

For example, a region with 3% velocity perturbation might need a maximum 47 kg/m³ average density anomaly in the crust relative to the background structure to be in isostatic equilibrium within uncertainty. As discussed herein, a purely compositional source of heterogeneity would predict a 65 kg/m³ density anomaly. To decrease this density anomaly to 47 kg/m³, at least 50% of the velocity perturbation (1.5%, or 0.098 km/s for a 6.5 km/s background) is assigned to temperature. Following Equation 3, this crustal section is then at least 195 °C cooler on average than background.

We limit the thermal perturbation to a geologically plausible range of ±250 °C. Miocene Moho paleotemperatures of ∼350 °C (Molnar and Jones, 2004) are probably a lower bound on plausible Moho temperatures and are ∼1000 °C below asthenospheric potential temperatures (1350 °C). Thus, the Miocene crustal column would have an average temperature anomaly of 500 °C below the hottest plausible crust if geotherms are approximately linear. Given these limits, we calculate crust-derived topography (Fig. 4B) that is a combination of thermal and chemical factors. This quantity represents an attempt to satisfy Equation 5 subject to the limits on allowable thermal variations in the crust discussed.

The result of these calculations compares well with thermal constraints from surface heat flow (Figs. 4C, 4D). The regions in which colder crustal geotherms are employed match almost exactly with areas of low reduced heat flow (Saltus and Lachenbruch, 1991) and low surface heat flow (Erkan and Blackwell, 2009), <25 mW/m² and 40 mW/m², respectively (Fig. 4C).

We recognize that loads in the sublithospheric mantle are only partially coupled to the overlying crust and surface. Deeper loads (from velocities at 120 km and 170 km depth nodes) are reduced in amplitude to account for viscous response in the upper mantle (50% and 25% coupling, respectively) (Parsons and Daly, 1983). This reduction is an approximation, as the degree of reduction in topographic expression is a function of wavelength and viscosity structure. We have chosen these factors to be approximately correct for loads with wavelengths of ∼100 km in an isoviscous mantle. Between end-member cases of full coupling and complete decoupling of the lower 100 km of the model from the surface, the mean change (absolute value) in predicted elevation at any given node is 140 m. Our approximation cannot produce errors of greater mean absolute value than 100 m, the difference from the fully coupled end member.

By smoothing the elevation field, we remove the effect of short-wavelength topographic variations generated by erosion and/or deposition, fault displacement, and volcanism (i.e., removing the surface loads of Lowry et al., 2000). The elastic filter is based on an estimate of flexural rigidity (3.5 × 10²², elastic thickness 18 km) from basin geometry of the San Joaquin Valley, adjacent to the range (in accord with Saleeby and Foster, 2004). Although elastic thickness may differ across the region, varying the elastic thickness from 5 to 30 km changes our predicted topography by <250 m and not in a manner affecting the overall pattern of our results. Because flexural loads adjoining our study area will influence topography, we examined these edge effects by extending our study outward 100 km into the region with tomographic coverage but no reliable receiver functions (and assuming a crustal thickness similar to the nearest receiver function). These effects are <200 m for models with 15 km elastic thickness and only affect the margins of our model.

The crustal and mantle topography are smoothed by the same flexural filter.

Uncertainties in seismologically based estimates of topography are calculated from the errors of the parameters used. The most important parameters in our error analysis are the velocity-density relations (Christensen and Mooney, 1995; Hacker and Abers, 2004; Jackson and Faul, 2010), wave speed perturbations, crustal thickness variations, and background asthenospheric density. Uncertainties (2σ) are estimated as 60 kg/m³ for the velocity to density conversion (Christensen and Mooney, 1995), 0.25% for velocity perturbation in the mantle, 3 kg/m³ per 1% Vp perturbation, and 50 kg/m³ for asthenospheric density. Crustal thickness uncertainties were given in Frassetto et al. (2011) and range between 1 and 5 km. These values yield typical 2σ uncertainties of the resulting synthetic topography of ∼600 m (Fig. 5A).

Systematic biases in our analysis may also be present to the degree that the seismic data (tomography and receiver functions) are systematically biased. While traveltime delay along a ray path is a robust observation in traveltime tomography, the depth extent of velocity anomalies is less certain. We therefore consider the effects of the limited vertical resolution of teleseismic tomography, in particular the placement of anomalies from the mantle into the crust or vice versa. We have therefore conducted the same analysis under the physically unreasonable end-member cases that all traveltime anomalies produced in the top 95 km are due to wave speed variations solely in the crust or solely in the mantle. Because the velocity to density scaling relationships are different in the two domains but of the same polarity, the effect is one of magnification. Heterogeneities in predicted elevation (and thus residual topography) are maximized if the mantle is homogeneous and the crust is responsible for all traveltime delays. Because forcing all wave speed variations into the crust generates crustal velocity heterogeneity of several tens of percent, values well beyond plausible, this calculation bounds our estimated topography’s sensitivity to systematic biases in the tomography. The mean magnitude of change is <100 m, and ∼95% of the elevation predictions are affected by <400 m.

A similar end-member investigation explores the robustness of our observations in the face of receiver function errors. In our analysis, we have treated the crust-mantle boundary as laterally variable. Nevertheless, identification of the Moho in many areas, especially those with a Moho significantly different from ∼35–40 km depth, is challenging (Frassetto et al., 2011). Therefore, we explored similar topography calculations with a constant crustal thickness. To a great degree, the results of calculations with this architecture are consistent with those honoring receiver functions, as the practical effect of this change is to change the domain (crust or mantle) in which velocity anomalies are located, and consequently predicted elevations differ by as much as ∼100 m.

Furthermore, anelastic effects are not considered in the crust. The signature of low crustal velocities (thus high H_c) and high estimated crustal temperatures might be indicative of supersolidus temperatures. If melt is present in the crust, density estimates will be too low (as described herein), and positive gravity residuals will be produced.

As seen in Figure 5B, elevations in much of the study area are nearly matched by a combination of compositional and (geologically reasonable but not directly observed) thermal variations in the crust and thermal variations in the mantle. Furthermore, nearly all of the misfit is within the uncertainty of topography estimates (Fig. 5C).

To explore possible errors from our scaling of mantle velocity into density, we also calculated results assuming a linear velocity-density relationship (i.e., one lacking correction for possible anelastic effects). The results are nearly identical (Fig. 5D), but with an additional –400 m of residual topography in the southern Cascades and southwesternmost Great Basin, near Death Valley. This residual is consistent with small amounts of partial melt or with anelastic effects, similar to the interpretations made here.

A logical test of our three-dimensional density model is a comparison with the observed Bouguer gravity anomaly, although we note that the limited ability to resolve the vertical position of density anomalies within the crust will lead to large uncertainties, even at the >100 km wavelengths of greatest interest here, and the limited width of the study area means that significant gravity anomalies could be generated from outside the model. Plausible ranges in density structures fitting the gravity field alone were illustrated in Kennelly and Chase (1989) and Jones et al. (1994). We have calculated the Bouguer gravity field predicted by our seismically derived density model (Fig. 6B), which largely recovers the observed Bouguer gravity variations (Fig. 6A). The most striking feature of a Bouguer gravity map of the region is the ∼300 mGal decrease from the foothills to the range crest. Our model recovers >200 mGal decrease into the range. Our model also recovers Bouguer gravity along range strike both north-northwest and south-southeast of the range crest. Remaining differences reflect a combination of masses outside the model bounds and the uncertainty in the vertical placement of density anomalies in the crust.

We note here that our gravity model seems to preclude the presence of variations in dynamic topography. To illustrate, consider a region at sea level with isopycnic mantle lithosphere (i.e., H_m = 0) and 40-km-thick crust of uniform density. If in isostatic equilibrium (i.e., H_c = 2.4 km), the crust must be 3008 kg/m³. If a ∼1 km of dynamic topography (basal normal force of ∼30 MPa) is being generated by asthenospheric convection (i.e., H_c = 1.4 km), then crustal density is 3088 kg/m³. The difference in the gravity signal from these two crustal columns is ∼135 mGal in the infinite slab limit and 118 mGal if active over a 300 km wavelength. Thus, given the absence of large-magnitude gravity residuals, we argue that the density structure that we estimate, and not dynamic topography, is responsible for the modern topographic relief of the Sierra Nevada.

In our model, the mean crustal density is estimated to be between 2900 and 3000 kg/m³ beneath the foothills and ∼2750–2800 kg/m³ beneath the range crest (Fig. 7B). The ∼50–55-km-thick foothills crust provides ∼5 km of topography while the thinner (40–45 km thick) but lighter crust under the range crest provides ∼6 km of topography (Fig. 4B).

The remaining ∼1 km topographic difference between the two areas is generated within the upper mantle (Fig. 7D). There is a 10% velocity increase (Jones et al., 2014) over a lateral distance of ∼150 km from beneath the range crest southwestward into the foothills (Fig. 2C). With the anelastic scaling described here (following Hacker and Abers, 2004; Jackson and Faul, 2010), this velocity anomaly (if due as assumed to thermal variation alone) represents a density contrast of 84 kg/m³. Thus 1.0 km of relief is caused by this upper mantle anomaly, a magnitude comparable to the relief caused by crustal buoyancy variations.

As an alternative to decomposing the topography field spatially, into crustal and mantle components (Figs. 4A, 4B), one can separate topography mechanistically into thermal, compositional, and dynamic components. We present such decomposition in Figure 7, where all values are demeaned. Because we succeed in explaining topography (within uncertainty) in nearly all of the study area and closely recovering gravity, with purely thermal variations in the mantle, we limit the need for chemical variations in the mantle or a spatially varying dynamic component to topography to less than our uncertainty of ∼600 m. Thus, compositional topography (Fig. 7A) is simply crustal topography (Fig. 4B) with the crustal thermal topography (Fig. 7C) removed. Comparing the foothills to the range crest, we find little (∼100 m) difference in crustal thermal topography. Conversely, the mantle underlying the foothills provides ∼1 km less thermal support than under the range crest (Fig. 7D). The remaining 1 km of discrepancy is due to crustal composition and thickness (Fig. 7A).

These findings are consistent with the conceptual model proposed previously (e.g., Ducea and Saleeby, 1996; Jones et al., 2004; Le Pourhiet et al., 2006; Molnar and Jones, 2004; Saleeby et al., 2012; Zandt et al., 2004), that negatively buoyant material has been removed from the mantle beneath the southern Sierra Nevada and is now foundering to the west of the range. If the mantle beneath the range crest had resembled that of the foothills until the Miocene, elevations would have been ∼1 km lower. If either a thicker thermal boundary layer or compositionally antibuoyant mantle lithosphere had existed beneath the batholith, more uplift is admissible.

The ∼35-km-thick crust in the Sierra Nevada has relatively low velocities that are mainly due to lithology, rather than to high temperatures. The crust is thinner by 5–10 km but slower by ∼0.25 km/s beneath the Basin and Range and Cascades than beneath the Sierra Nevada. We propose that the lower wave speeds in the Death Valley region are due not to composition but to higher crustal temperatures and perhaps small average amounts of partial melt. The crust in the Death Valley region may be chemically denser (and thinner) than beneath the Sierra Nevada, and thus supports 1 km less topography than crust under the Sierra.

Beginning from seismic velocities, we have created a density structure in the Sierra Nevada region that reproduces modern topography and is nearly concordant with the observed gravity field. The model assumes thermal and chemical heterogeneity in the crust and temperature variations in the mantle, but we find no need for topography derived from mantle chemistry or asthenospheric convection (i.e., dynamic topography) in order to recover observed elevations. Uppermost mantle densities account for ∼1 km of relief between the Sierran crest and foothills. If such material was present beneath the range in the Miocene, we estimate ∼1 km of surface uplift upon its removal and regional isostatic equilibration. Our work is consistent with the hypothesis that cold upper mantle has been removed from the southern Sierra Nevada but not the western foothills. Chemically buoyant crust of moderate thickness beneath the high southern Sierra Nevada supports >1 km of topography relative to the thinner, warmer, chemically denser crust of the Basin and Range and southern Cascades.

The Sierra Nevada Earthscope Project was funded by NSF grant EAR-060783. We particularly thank landowners for permission to place seismometers on their property. The SNEP deployment was made possible by the hard work of a large team of students and researchers. We thank Tony Lowry and Gene Humphreys for thoughtful reviews.

We here detail the means of reconciling the variations in Moho depth with the P-wave structure derived assuming a flat Moho. Consider two locations with identical teleseismic P-wave speed anomalies but different Moho depths. To honor the teleseismic delay times, the integrated traveltime along teleseismic ray paths must be the same as that calculated with no Moho variation. Where the Moho is deeper, crustal and/or uppermost mantle wave speed anomalies must be higher than directly computed from the P-wave tomography, because rays must traverse a thicker section of crustal material that is initially assumed to have a lower Vp than the mantle. For example, a 10-km-thick section of lower crust (Vp = 6.8 km/s) is expected to result in a 0.23 s delay relative to 10 km of asthenosphere (Vp = 8.04 km/s).

We assign this deviation uniformly through the seismic model from the surface to 95 km. We use such a broad area for three reasons. First, this avoids generating large, probably unrealistic velocity anomalies adjacent to the Moho. Second, this approach is consistent with inversions of synthetic data sets (Reeg, 2008; Jones et al., 2014) that show that vertically compact velocity anomalies at these depths will tend to be smeared vertically over several adjacent layers in the tomographic model. Third, because the pertinent quantity in P-wave tomography is integrated delay time along the raypath and the pertinent quantity is prediction of topography is the integrated density through the lithospheric column, the predicted topography is virtually unaffected by this smearing and so is a robust result of the tomographic inversion. Loss of vertical resolution is accompanied by a decreased magnitude of the anomaly such that the predicted topography integral is within 2% in the case of a checkerboard test (without any discrepancy between crust and mantle); the low pass filter of lithospheric flexure further smooths the surface response.

Modulating absolute Vp would have a small effect on receiver function interpretations conducted without such changes. The greatest Moho depth (55 km) results in a dVp_Moho of 3.0%. If the mean crustal Vp is 6.5 and Vp/Vs is 1.8, then the difference between Ps and direct P travel time ts-tp=6.77s. Honoring this delay time and perturbing Vp and Vs by 3% (to 6.695 and 3.72 km/s, respectively) will lead to a Moho depth estimate of 56.65 km. In turn, a 1.65-km-thick section of lower crust with dVp = 3% provides of 43.5 m of topography while a 1.65-km-thick section of upper mantle with dVp = 3% supports –15.5 m of elevation. The most extreme changes necessary would therefore be 1.65 km crustal thickness and 59 m in predicted elevation.