P-wave tomography of potential convective downwellings and their source regions, Sierra Nevada, California

The Sierra Nevada of California, stretching ∼500 km in length and rising to peaks above 4 km in elevation (Fig. 1), has proven to be a tectonic enigma. Underlain by Mesozoic granitoids and older metamorphic wall rocks (e.g., Fig. 2; Evernden and Kistler, 1970), the range’s elevation was initially thought to be in isostatic equilibrium with an Airy crustal root (Byerly, 1937; Lawson, 1936). Although later gravity measurements confirmed the presence of a compensating mass deficit (e.g., Oliver, 1977), the highest part of the range was found to have a crust only ∼35 km thick, only slightly thicker than that of the Death Valley region roughly 2 km lower (Fliedner et al., 1996; Jones and Phinney, 1998; Jones et al., 1994; Ruppert et al., 1998; Wernicke et al., 1996). Paradoxically, thicker crust is found under the low western foothills (Fliedner et al., 1996, 2000; Jones and Phinney, 1998; Ruppert et al., 1998). Support for at least the southern part of the range relative to the western foothills is therefore thought to at least partially originate in the mantle, an inference consistent with Quaternary xenoliths, magnetotellurics, and seismic wavespeeds (Ducea and Saleeby, 1996; Fliedner et al., 1996; Jones et al., 1994; Lee et al., 2000; Park, 2004) indicating the presence of asthenospheric material at or only a short distance below the modern Moho.

Compensation in the mantle indicated that material was removed from beneath the eastern, high part of the Sierra in the late Cenozoic. Initially, some workers had suggested that the mantle lithosphere had been extensionally thinned or removed from under the eastern Sierra (Jones, 1987; Jones et al., 1994; Mavko and Thompson, 1983), but most xenoliths from 8 to 12 Ma volcanics in the southern Sierra indicated that the crust had been perhaps 65 km thick at that time, with the lower half being a mafic rock containing considerable amounts of garnet (Ducea and Saleeby, 1996, 1998). Such xenoliths are absent from volcanics 3 Ma and younger as well as the 8.6 Ma site at Blue Knob, which do contain xenoliths from this depth range. Thus Ducea and Saleeby (1996, 1998) inferred that this dense material was removed since 8–12 Ma, a scenario that would cause the Sierra to gain elevation in the absence of any cryptic crustal thinning (e.g., Jones et al., 2004). Removal at ca. 3–4 Ma has been suggested because of a short-lived burst of potassic volcanism with unusual isotopic characteristics (Farmer et al., 2002; Lee et al., 2000; Manley et al., 2000). A late Cenozoic age for removal of this material is also consistent with heat-flow measurements in the Sierra that are far too low to be in steady state with such shallow asthenosphere (Crough and Thompson, 1977; Saltus and Lachenbruch, 1991).

With the suggestion that dense material had been removed from the southern Sierra, attention turned to where this material might be at present. A high-wavespeed body in the upper mantle under the southeastern Great Valley and southwestern Sierra had been identified previously (Benz and Zandt, 1993; Biasi and Humphreys, 1992; Humphreys et al., 1984; Jones et al., 1994; Raikes, 1980); variously termed the Southern Great Valley anomaly, the southern Central Valley anomaly, and the Bakersfield anomaly, we will call it the Isabella anomaly (reflecting its first identification by Raikes [1980] at station ISA). Data from a 1997 deployment led to new images (Boyd et al., 2004) and an interpretation that the Moho was distorted by this body (Zandt et al., 2004). The geometry of this body has varied between studies, from being near vertical (Benz and Zandt, 1993; Biasi and Humphreys, 1992) to perhaps plunging somewhat to the west (Jones et al., 1994) to plunging to the east or southeast (Boyd et al., 2004; Yang and Forsyth, 2006). Differences in the geometry are largely due to differences in coverage and resolution of the different studies.

The Isabella anomaly had previously been interpreted in many ways, including being a fragment of slab (Benz and Zandt, 1993) or a convective downwelling (Aki, 1982) of Great Valley lithosphere (Zandt and Carrigan, 1993), sub–Tehachapi Mountains lithosphere (Jones et al., 1994), sub-Sierran lithosphere (Biasi and Humphreys, 1992), or Coast Ranges and Mojave Desert lithosphere (Biasi, 2009). Subsequent to the recognition of late Tertiary removal of a mafic lower crust from the Sierra Nevada to the northeast, suggestions focused on the Isabella anomaly being that material, having been removed convectively (Molnar and Jones, 2004; Zandt, 2003) or through something closer to delamination (Le Pourhiet et al., 2006). These suggestions have been bolstered by the inference of Saleeby and Foster (2004) and Saleeby et al. (2013) that the Sierran foothills and San Joaquin Valley above and adjacent to the Isabella anomaly have subsided in Quaternary time. Even among these interpretations, there is disagreement about what exactly is in the Isabella anomaly: Zandt (2003) suggested that the upper mantle anomaly was lithospheric mantle entrained by removal of mafic crust that is now at much deeper levels in the mantle, but Saleeby et al. (2003) argue that the mafic crust must reside within the upper mantle anomaly. Biasi (2009) argues against a sub-Sierran source from calculated discrepancies in volumes between sources and sink; he instead prefers some downwelling induced by convergence in the mantle. A sub-Sierran origin has also been challenged by Wang et al. (2013), who resurrect the inference that this is a piece of Farallon slab (part of the mostly subducted Monterey subplate).

Because the geologic and geophysical constraints on lithospheric foundering were largely in the southern Sierra, focus on upper-mantle anomalies was directed to the adjacent Isabella anomaly, but a somewhat similar anomaly had been imaged in northern California by Benz and Zandt (1993). They interpreted this Redding anomaly as part of a stagnating slab associated with subduction of the adjacent Juan de Fuca plate, but Jones et al. (2004) suggested that the anomaly, which is unusually large in magnitude and volume when compared with images of the slab to the north, might also contain descending sub-Sierran lithosphere.

The possibility that the Sierra Nevada is the product of foundering of mantle lithosphere and/or mafic lower crust and that this foundering material might still remain distinct in the upper mantle led to a pair of major initiatives: the Sierra Nevada Earthscope Project (SNEP) focused on seismological imaging of the lower crust and upper mantle under the Sierran region, and the Sierra Nevada Drips Continental Dynamics Project focused on the geologic history, geodynamic constraints, and magnetotelluric imaging of the range. The overall goal of these projects is to determine the extent of lithospheric foundering, its impact on the geology of the region, and the mechanism of removal. We report here on the P-wave tomography facet of SNEP. A companion work growing out of the same field experiment has found that the Moho geometry of the southern Sierra appears to characterize much of the range, with a deep, poorly defined and probably gradational Moho under the western foothills contrasting with a sharp, shallower Moho under the eastern Sierra (Frassetto et al., 2011).

The main goals in the tomographic component of SNEP are (1) to determine the extent of the range currently lacking the Mesozoic mafic lower crust identified by Ducea and Saleeby (1996, 1998) and (2) to determine the characteristics of the Isabella and Redding anomalies and their possible relationship to the Sierra Nevada. To that end, 84 new sites were occupied for at least a year with broadband seismometers from the FlexArray pool from EarthScope (Supplemental Table 1¹; Figs. 1 and 2). Seismograms from these stations and USArray Transportable Array (TA) stations in the region were picked for P-wave arrival times, which were merged with arrival times to older, more localized deployments in the southern Sierra in order to develop a broader and more uniform picture of the wavespeed variations in the lower crust and upper mantle through the Sierra and adjacent areas. A companion S-wave study is still under way.

SNEP deployment from May 2005 to September 2007 was the first deployment of the EarthScope FlexArray instruments and overlapped with the USArray Transportable Array (TA) in this region. The seismometer network spanned the Sierra Nevada, extending 150 km in width (SW-NE) and 600 km in length (NW-SE). Southerly stations were occupied in phase I of the experiment (May 2005 to mid-2006) and northerly stations occupied in phase II (mid-2006 to September 2007). The SNEP array overlapped with the northern portion of two previous studies done in the area, the Sierran Paradox Experiment (SPE) and a 1988 deployment (Boyd et al., 2004; Jones et al., 1994) (Fig. 1). Within the combined SNEP and TA seismometer array in the Sierra, the majority of stations were ∼25 km apart (Supplemental Table 1 [see ^{footnote 1}] and Fig. 1).

Nearly all of the 84 locations occupied during the SNEP deployment had three-component, broadband CMG-3T seismometers recorded at 40 sps (samples per second) by Reftek R130, 24-bit, three-channel digitizer data acquisition systems. The exceptions were four 1-s CMG-40T seismometers used from September 2006 to January 2007 and a few stations set to 20 sps when access was expected to be too infrequent for flash disks to be replaced before filling. In general, one year’s worth of data was collected at each site; a few stations were occupied for two years to overcome equipment failures in an attempt to achieve this goal. Eleven additional broadband seismometers were deployed in the summer of 2007 along the Tioga Pass Road in Yosemite National Park at 5 km spacing. All SNEP data are archived with the Incorporated Research Institutions for Seismology Data Management Center (IRIS DMC) as data set XE in 2005–2007 (http://www.iris.edu/gmap/XE?timewindow=2005-2007).

The 2005–2007 SNEP deployment was supplemented with seismograms from permanent stations lying within a large polygon for the time period of May 2005–September 2007 (Figs. 1 and 2). The area of study measured ∼700 km east to west and 900 km north to south. Seismograms came from networks operated by USArray Transportable Network; Berkeley, California Institute of Technology; University of Nevada, Reno; University of Washington; University of Oregon; U.S. National Seismic Network; and the Leo Brady Network (Supplemental Table 2²). The sampling rate from these stations was usually 40 sps, but a few stations were sampled at 25 sps or 50 sps. Such seismograms were resampled to 40 sps using Gary Pavlis’s dbresample routine for Antelope 4.8 (xcor 1.2) (Antelope contributed software, http://antelopeusersgroup.org/repository).

Seismograms from 1 min prior to 2 min after the predicted P-wave arrival time of more than 2000 teleseismic events of magnitude 4.5 or greater were extracted from continuous records. All events at distances between 28° and 92° from N38.4°, W120° and with a magnitude 5.5 or greater were included. Events with a smaller magnitude from underrepresented backazimuths were also included to improve azimuthal coverage.

The older experiments providing P-wave arrival times were shorter in duration but occupied the same general area in the southern Sierra (Supplemental Table 3³ and Fig. 1). The SPE deployment used 22 three-component broadband seismometers recording continuously at 40 sps over the summer of 1997; these data are available from the IRIS DMC as data set XJ in 1997. The 1988 experiment recorded in a triggered mode a mix of 18 three- or single-component 1 s short period stations; additionally, 41 stations of the U.S. Geological Survey–California Institute of Technology Southern California Network were picked for the 1988 data set (Jones et al., 1994). The 1997 and 1988 experiments shared some sites, as did the SNEP and 1997 deployments.

New P-arrival picks were made using dbxcor, a cross-correlation, beam-forming program described by Pavlis and Vernon (2010). All waveforms were initially filtered with a Butterworth band-pass filter from 0.25 to 10 Hz (0.1–4 s), and additional filters were applied to some events in order to increase the signal-to-noise ratio. Seismograms were aligned to an optimal beam derived from all stations recording the event and then arrival times derived from the difference with the predicted arrival time from the IASPEI91 velocity model (Kennett and Engdahl, 1991).

A total of 81,210 P arrivals were picked using dbxcor. Because the number of P arrivals in the full data set was skewed toward a few highly active backazimuths (Table 1 and Fig. 3B), a subset of these picked P arrivals was chosen for use in the travel-time inversion in order to minimize azimuthal bias in locations that existed in the full data set. In total, we picked P arrivals for 653 events, 319 during phase I and 334 during phase II. We removed events from well-sampled backazimuths that were recorded by fewer stations than other events from the same backazimuth and distance. Events recorded at less than ∼50% of the stations were generally dropped. After this winnowing, the 206 events with 26,856 arrival times that remained sampled available backazimuths more evenly (Table 1 and Fig. 3; Supplemental Table 4⁴).

To increase the number of arrivals in the southern portion of the Sierra Nevada, we supplemented our picks with 1231 new dbxcor picks from 57 teleseismic events from the 1997 SPE deployment (Boyd et al., 2004) and 1099 picks from 41 events recorded in 1988 previously used by Jones et al. (1994) (Supplemental Tables 5⁵ and 6⁶; Fig. 3). All of the older events share some station locations in common with later deployments. The complete set of P-wave arrival times are included in Supplemental Table 7⁷.

Of the 29,186 picks used, about a third came from the core SNEP deployment area.

Patterns of P-wave residuals indicate the presence of large lateral variations in P wavespeeds at upper mantle depths (Fig. 4I). From backazimuths 50°–125°, the stations in northern California north of ∼39°N register residuals that are 0.5–1.5 s early. This is likely the effect that Benz and Zandt (1993) attributed to the subducting Gorda plate, as seismic waves travel relatively quickly up the dipping slab. In the southern portion of the study area, the early arrivals due to the Isabella anomaly can clearly be seen as a cluster migrating with changing backazimuth and incidence angle most clearly identified for southeastern and northwestern backazimuths. Greater changes in the position of the anomaly with backazimuth and incidence angle reflect a greater depth of the causative anomaly. A more diffuse and lower-amplitude set of early arrivals is seen in the western Sierra, dominantly for events from the west.

Delayed arrivals are less distinct, but a noticeable band of late arrivals is seen between the Sierran crest and the California-Nevada border. This anomaly is strongest from eastern backazimuths, with largest anomalies along the Sierran crest, while from western backazimuths, the amplitude of this feature decreases and shifts somewhat to the east to the state border.

We used an iterative, spherical, finite-difference (FD) tomography code, sphfd, created by S. Roecker (Rensselaer Polytechnic Institute) (Roecker et al., 2006), to invert the 29,186 P-arrival times for three-dimensional perturbations in P wavespeed (V_p). This code parameterizes a volume by nodes laid out on a rectangular grid at specified depths. Seismic slownesses vary linearly between nodes. Raypaths are found through a starting 3-D structure using the finite-difference formulation of Hole and Zelt (1995) based on that of Vidale (1988). Partial derivatives of travel times with slownesses are accumulated to produce a system of linear equations relating slowness changes to travel-time residuals. A damping term is added primarily to limit changes in any single iteration. Inversion of this system of equations uses the LSQR algorithm of Paige and Saunders (1982b). Adjustments to model slownesses are then added to the starting model, and the process is repeated starting with the calculation of new raypaths through the new model. We used the IASPEI91 earth model as an initial 1-D seismic velocity structure except where replaced by a starting 3-D structure (Kennett and Engdahl, 1991).

Final results depend on the starting model, node geometry chosen, the damping factor used, smoothing of the model, and exclusion of blunders from the data set. We discuss the effects of starting model more below; the following discussion applies to nearly all inversion runs investigated. Nodes were separated by 25 km horizontally, roughly the station spacing in the densest part of our network, and at depths of –30, 0, 20, 40, 70, 120, 170, 220, 270, 320, and 430 km below sea level. Separating the nodes by greater distances yielded poorer fits to data, while the 25 km spacing did not produce unacceptable artifacts.

The relative success in minimizing the two variances is controlled by a damping factor λ, which reduces deviations from the starting model as it increases and is implemented as a scalar in the LSQR matrix inversion routine used (Paige and Saunders, 1982a, 1982b). We explored a range of possible values, selecting a value of 250 as providing the best tradeoff of fit to the data to a smooth model (Fig. 5). In most inversions, λ was relaxed to 150 in the last few iterations as the model approached a stable solution. Relaxing λ below that value resulted in the inversion diverging from a stable solution.

We chose to further enforce a smooth model by applying an a posteriori smoothing function. The smoothing function is a simple moving-window average of the change in that iteration to a node’s slowness with that of the 26 adjacent nodes in the x, y, and z directions; this was applied twice after each iteration. Inversion runs with less smoothing yielded comparable reductions in the data variance but much higher model variances.

We also attempted to remove blunders in the form of cycle skips that yield arrival times generally more than 1 s different from a correct arrival time. Such outliers tend to distort the success of a least-squares algorithm and are best removed. We removed individual arrival times with residuals above a threshold τ_i prior to each iteration. The threshold τ_i declined from 3 s at the start to 1.5 s at the seventh iteration to 1.0 s at the 12th iteration. Toward the end of the inversion, this process removed ∼130 of the 29,186 arrival times (∼0.5% of the data).

We have made no attempt to remove any signal from variations in the depth of the Moho. This is mostly because of the challenge in identifying the Moho in the Sierran foothills (e.g., Frassetto et al., 2011; Zandt et al., 2004) and the strong indication that the velocity contrast across the Moho varies dramatically through this region. We do not think this produces significant difficulties in interpreting our results. A 10 km increase in the depth of a normal Moho contrast (6.8 km/s material over 8.0 km/s material producing a travel-time anomaly just over 0.2 s) represents a 10-km-thick layer with a 15% low-wavespeed anomaly; averaged over the 25 km that the 40 km layer of nodes controls, this would be reduced to a –6% wavespeed anomaly (other combinations are illustrated in Fig. 6). It is further reduced if vertical smearing reduces the amplitude, as discussed below. The generally poor correlation of anomalies near Moho depths with inferred Moho topography suggests that very little of our tomographic signal reflects such topography and that other variations are far more important.

As a general guide to the success of our inversion in resolving wavespeed anomalies, we present results of simple inversions of some simple alternating anomalies (checkerboards). Travel times are calculated through the trial model for event-station pairs used in the inversion of real data; this synthetic noise-free data set is then subjected to the same inversion procedure as the real data. The difference between the inverted result and the starting model provides an indication of the tendency of the inversion to smear compact anomalies and reveals areas where damping and inadequate ray coverage prevent a robust solution. A similar approach with specific hypothetical structures will be presented with our interpretations of the actual inversion.

For our checkerboard tests, we created a trial model with square seismic anomalies 75 km (four nodes) on a side with an amplitude of ±3% in a single layer of an otherwise homogeneous layered halfspace. One model had this layer at 170 km depth (Figs. 7 and 8; near the center of the Isabella anomaly), the other at 40 km depth (Figs. 9 and 10). The inversion was iterated 14 times; variance was reduced 92%–95% in each case. About 60% of the original anomaly amplitude was recovered in the 170 km depth layer; at 40 km, recovery was closer to 35%–45%. The vertical integral of the slowness anomalies, however, was between ∼50% and 140% of the original value, with the lowest integrals at the corner nodes of anomalies and the highest values in the center four nodes of each anomaly. When vertically integrated over the full depth range of the inversion and then averaged over the horizontal extent of the original checker, ∼75%–95% of the original slowness perturbation is recovered at 40 km depth, and ∼85%–97% of the original slowness perturbation is recovered at 170 km depth.

The general result of these and other checkerboard tests is that the horizontal position and centroid depth of anomalies are well recovered within the area of interest (generally the area where anomalies are recovered in Figs. 7 and 9), though the average amplitude of anomalies tends to be underestimated and some vertical smearing of perturbations into layers above and below the actual position of an anomaly occurs. The vertical smearing is not an effect of the smoothing operator used: inversions without smoothing yielded nearly identical vertical smearing. Instead, such smearing reflects a fundamental inability of teleseismic travel-time tomography to distinguish a vertically compact anomaly from one with greater vertical extent but somewhat smaller amplitude. The bias in a least-squares inversion constructed in a manner such as ours is toward minimizing the amplitude of anomalies; so the smeared, lower-amplitude anomaly is preferred. Some additional tests indicate that individual anomalies down to at least a 50-km diameter can be recovered.

Although this vertical smoothing might suggest that vertically alternating high- and low-wavespeed anomalies would smear together and be unresolved, a separate test found this not to be the case if anomalies are sufficiently separated. A checkerboard like those used above at depths from 20 to 70 km overlying an inverted version of the same checkerboard at depths from 170 to 270 km was used to generate synthetic data, which was then inverted. The resulting wavespeeds recover individual anomalies much better than the individual checkerboards shown here, with individual nodes in the horizontal center of anomalies recovering 95% of the original anomaly amplitude even in the 70 and 170 km layers. Based on this, we consider the individual layer checkerboard tests to be worst-case scenarios; we expect our results to be closer to the real earth than these tests are to the initial checkerboards.

We consider the effect of noise by adding to “perfect” synthetic data values drawn from a Gaussian distribution with a variance equal to that of the travel-time residuals at the close of our inversion of real data. The overall standard deviation of the difference between inverting perfect data and perfect data + noise is ∼0.3% of the mean wavespeed (i.e., ∼0.02 km/s in the crust and 0.03 km/s in the mantle). Thus 95% of the nodes have velocity perturbations within ∼0.6% of the perfect inversion results. The larger deviations tend to be at the periphery of the volume sampled by raypaths used in the inversion; from 40 to 320 km directly below the Sierra, errors from noise almost never reach 0.5%, and these largest errors are spatially limited to one or two adjacent nodes. We infer that the 95% uncertainty on anomalies under the Sierra and immediately adjacent terrain with dimensions of 50 km (three nodes across) is ∼0.3% (meaning that a velocity anomaly of 1% could be between 0.7% and 1.3%). As we show below, this is far below the magnitude of any anomalies we interpret; the greater uncertainty is caused by the smearing issues already discussed. The contour interval on our maps and sections of 0.5% is roughly equal to uncertainty from noise.

Thus the resolution tests indicate that the amplitude of anomalies recovered in this inversion is likely to be lower than in the real earth and can be more smoothed vertically than in the earth. The central depth of an anomaly should be recovered and its lateral extent correctly delineated. The integrated magnitude of an anomaly should be somewhat smaller than the actual anomaly, but only low by 5%–25%, with better recovery for anomalies below the crust. The magnitudes of smearing and decrease in anomaly amplitude depends on the exact geometry; we expect that in the real earth, both smearing and decrease in amplitude will be less than exhibited in inverting checkerboard patterns. Because of this, we will present some special tests in interpreting certain facets of the recovered wavespeed structure after first summarizing our findings.

A well-known characteristic of teleseismic body-wave tomography is that there is no vertical resolution near the surface (roughly from the surface to depths about equal to interstation distances); to explore the limitations of our tomography, we not only began our inversions from a 1-D model (with [GV 1-D] and without [1-D start] station corrections for sediments in the Central Valley), but we also conducted several other inversions starting from 3-D models derived elsewhere. Inversions spanning large regions often correct for an assumed crustal model and focus exclusively on deeper anomalies at some risk that errors in the crustal model will propagate down into the deeper structure. Long-wavelength crustal variations are a likely problem in this region, because the deep Great Valley sediments have much lower wavespeeds than the crystalline rock of the Sierra. With our unusually close station spacing for this aperture, we should have some resolution even at depths of 20 km, where rays from adjacent stations begin to cross, and our tectonic interpretations would be strengthened by constraining wavespeeds up into the crust. Surface-wave models, either generated from ambient noise (e.g., Moschetti et al., 2010; Yang et al., 2008) or earthquakes (e.g., Gilbert et al., 2012), generally provide a better vertical resolution than teleseismic body-wave tomography but suffer increasingly poor lateral resolution with depth. As such, they are a useful complement to our analysis.

We supplemented inversions starting from a 1-D structure with a series of inversions that started from three surface-wave models for this region (Table 2). Moschetti et al. (2010) used both ambient noise and surface waves from earthquakes recorded dominantly at EarthScope TA stations to construct an isotropic shear-wave model for the western United States. The version we use here was constructed without an assumed Moho (which is not pictured in the published paper). Gilbert et al. (2012) used teleseismic Rayleigh waves recorded at SPE, SNEP, and TA stations to construct an S_V shear velocity model of California and western Nevada. Stachnik et al. (2011) constructed a shear velocity model of the Sierra and immediately adjacent region from recordings of ambient noise at SNEP and TA stations; because the area was smaller than this tomography, it was also used embedded within the Moschetti et al. (2010) model. This tended to produce artifacts near the Stachnik model edges and yielded higher residuals and so discussion of this model will be minimal. These three structures have commonalities and differences that permit us to explore the effect of varying shallow structures on deeper interpretations (especially in the vicinity of the Moho).

A number of assumptions were required to adapt these models for use here. We converted these shear-wave models to compressional wavespeeds using the v_p on v_s relation of Brocher (2005), which captures the likely systematic variations of the v_p/v_s ratio between felsic and mafic rocks. The original models have a very different parameterization, with much thinner layers than in this study. We averaged the wavespeeds derived from the surface-wave models such that the vertical transit time of a P-wave would remain the same in the original and reparameterized structures. The errors bound to be introduced by these assumptions should be much smaller than the large (±10%) wavespeed variations common to the surface-wave models.

At the resolution used here, these three models are similar above 40 km depth but diverge below that depth (Figs. 11IA–11IC). All the a priori models have relatively high wavespeeds in the Sierra and very low wavespeeds in the Central Valley and adjacent Coast Ranges in the top layer (an average of 0–10 km depth), with about a 25% contrast between these extremes. In the 20 and 40 km depth layers, high wavespeeds underlie the Central Valley and Coast Ranges, while low wavespeeds are under the eastern Sierra and western Great Basin; the models differ in the position of the boundary between these regimes, varying from just west of the Sierran crest to near the boundary with the Central Valley. At and below our 70-km-depth layer, the a priori models diverge greatly in the Sierra, though they all tend to show slow upper mantle velocities to the east of the Cascades and northern Sierra and an east-west–oriented band of low-wavespeed material extending across the northern Sierra.

The 3-D starting models based largely on the Moschetti et al. (2010) model reduced the data variance relative to a 1-D starting model by ∼6%–7%, a small fraction of the ∼75% overall variance reduction achieved after nearly all the inversions (Table 2). They were better in the core SNEP area, where variance was reduced ∼20%. The Gilbert et al. (2012) model tended to have a higher variance than a 1-D starting model and a much higher model variance than the Moschetti et al. (2010) model. Presumably the reason for the relatively small reduction is a combination of the different sensitivities of the two techniques, imperfections in the conversion from S-wave to P-wave structures, and errors in the starting models. Although ideally the surface-wave and body-wave observations would be jointly fit, such an inversion is beyond the scope of this paper.

We generally conducted a pair of inversions for any particular starting model: one where the surface-wave–derived structure was used as a 3-D starting model, but otherwise the inversion proceeded as in the 1-D cases (termed a “Free” inversion), the other where some number of top layers were held fixed (usually nodes at –30, 0, 20, and 40 km depth) until the inversion converged (this result termed “fixed”), and then those layers were freed and the inversion resumed until it converged once more (termed “fixed-free”). Changes over this latter procedure after all layers are free tend to reveal travel-time variations that most strongly require differences between the surface-wave and body-wave tomographic models in the crust and the very top of the mantle, while inversions that allow all nodes to vary from the start are less apt to enforce a real difference between the P and S structures. The models collectively give us insight into the null space of the teleseismic P tomography as well as the likely tradeoffs of altering any particular part of the model. We will bring up these results below where they provide insight into the uncertainties of the teleseismic work.

Some general guidelines in examining these results might be helpful. Where the amplitude of anomalies in the top 40 km of the “fixed-free” model is much higher than in the “fixed” result, the P residuals require a higher-amplitude anomaly at shallow depths than was derived from the surface-wave models, meaning that the S-wave model underestimated the anomaly either because our v_s-v_p mapping used a poor v_p/v_s ratio or because of some limitation in the v_s model itself, such as a limited horizontal resolution. Larger-amplitude anomalies in the 70 km horizon in the “fixed” inversion than the “free” inversion similarly indicate a strong disconnect between the starting model and the P residuals. In some cases this results in a unimodal anomaly in an inversion starting from 1-D becoming bifurcated in a “fixed” inversion; frequently the amplitude of the bifurcation is reduced in the “fixed-free” inversion, indicating that the P residuals are in conflict with some anomaly in the S-wave model. In some cases, a unimodal anomaly in the 1-D inversions will either bifurcate or be paired with an anomaly of an opposite polarity in the “free” as well as the “fixed” inversions. In these cases the P residuals have little resolving power and are compatible with the more complex S-wave model. The reader can apply these guidelines to anomalies not explicitly considered below to gauge the reality of such anomalies.

We first discuss the properties of the inversions starting from a 1-D layered model, including the main anomalies defined in these inversions. These are features that are required by the travel-time observations and must be present in some form in all the inversions. We then consider the details of certain parts of the inversion in light of the inversions started from the 3-D models. In general we note that there are strong differences in the upper 90 km of the models (i.e., nodes down to 70 km layer differ) but few differences below that depth. This is to be expected as the surface-wave models have less resolution at greater depth and the P-wave tomography has much greater resolution below ∼40–70 km. Some differences do exist, though, that reflect the non-uniform distribution of both stations and sources.

When starting from a 1-D structure, the variance of the unweighted travel-time residuals of the 29,056 arrivals from 304 events that remained after 14 iterations of the inversion is 0.0404 s², a 77% reduction from the 0.178 s² variance of residuals relative to a radially symmetric Earth. This results in a standard deviation of the residuals of 0.20 s, well above the 0.05–0.10 s uncertainty we would anticipate solely from picking waveforms. At stations within the core area of the SNEP deployment, variance was reduced from 0.136 at the start to 0.0281 at the end, a 79% variance reduction and a standard deviation of 0.16 s. Thus a significant part of the overall variance remaining at the end of the inversion comes from the periphery of the model, where greater station spacing, poor azimuthal coverage, and a greater chance for miscorrelation of arrivals due to changes in the waveform shape with distance from the core all contribute uncertainty. Changes between the last few iterations suggest that some of the remaining residual might contribute to a larger amplitude of some of the anomalies clearly evident in the current results that remain unrealized because of the damping term chosen; we thus expect that we tend to underestimate the amplitude of anomalies by a small amount. Within the Sierra Nevada itself, we expect to have captured nearly all the wavespeed variations producing arrival time variations, with most of the uncaptured variance reduction reflecting stronger gradients and higher amplitudes than our damping and smoothing parameters allow, although some might reflect 3-D effects not captured in our approach (e.g., theoretical arrival times for near-zero energy paths). The residuals remaining after the inversion is complete can be viewed in Supplemental Figure 5¹¹.

The resulting velocity structure (Fig. 11Ie, “1-D start”) contains about four main, distinct positive anomalies (high wavespeeds) and a less distinct set of low-wavespeed anomalies. The southern deep anomaly is the Isabella anomaly; high wavespeeds of the Isabella anomaly are evident between 70 and 220 km in the vicinity of 36°N, 119°W. North of the Sierra Nevada lies another, larger high-velocity body visible from 70 km to 220 km centered around 40.5°N, 122°W, previously referred to as the Redding anomaly (Benz and Zandt, 1993). This anomaly lies against or within the southern edge of the subducting Gorda plate; it is mainly distinguished by its higher wavespeeds and shallower top (Fig. 12I, profile E). A shallower high-velocity anomaly extends between these two deeper anomalies along the western foothills of the Sierra Nevada above 70 km; we follow Biasi (2009) in referring to this as the Sierran Foothills anomaly.

In contrast to the well-defined, high-velocity bodies, the low-velocity anomalies are generally less cohesive and have smaller magnitudes. The large-amplitude, low-wavespeed body to the east of the Redding anomaly (41°N, 121°W), termed the Gorda backarc anomaly, extends from the surface to 120 km depth (especially 40–70 km, Fig. 11I and section S3 of Fig. 12I); this could be the western edge of a broader low-velocity body beyond the limits of the stations used. Another large low-velocity anomaly, the Central Great Valley anomaly, exists beneath Central California from 170 km to 270 km at ∼38°N, 120°W (Fig. 11I, 170, 220, and 270 km layers). Near the Isabella anomaly to the south and west, a partial ring of low-velocity material can be seen at 70 km and 120 km. A smaller-amplitude, low-wavespeed body generally trends along the eastern Sierra near 40 km depth.

Although the basic pattern of the main anomalies remains as described above, a number of differences emerge that affect our interpretations. Some general aspects of the surface-wave–based solutions deserve mention. First, vertical variations in wavespeed invisible to the body-wave tomography, such as the relatively low wavespeed sediments of the Central Valley over the relatively high wavespeed basement underneath, are far clearer in these solutions and fairly consistent between models. Second, differences between the surface-wave models below ∼40 km can cause the resulting solutions to have different mean values in the core part of the study area; this can produce differences in the color scale that may make some local wavespeed contrasts appear different for different models when the local contrast is the same. The statistics of all these models are summarized in Table 2.

The Isabella anomaly is a prominent body ∼4% higher in P wavespeed than average that is distinctly recognizable at depths from ∼100 km to 220 km below the southern Central Valley in all inversions. We find the Isabella anomaly extends ∼120 km east-west and ∼100 km north-south with a center at 120 km depth at 36°N, 119.3°W when using the +2.5% contour for reference (sections A [Figs. 12I, 14] and B [Figs. 12I, 13]). There appears to be a ∼60°–70° plunge to the east below ∼100 km. The size, depth extent, and magnitude of the anomaly agree well with most previous P-wave tomography results (Benz and Zandt, 1993; Biasi, 2009; Biasi and Humphreys, 1992; Jones et al., 1994; Zandt and Carrigan, 1993). Our synthetic tests (e.g., Fig. 13) suggest that the wavespeed anomaly could be somewhat more compact laterally with a wavespeed ∼5% higher than average; these also indicate that the plunge of the anomaly should be recoverable below 100 km depth. Some connection to the crust seems plausible with all but the 1-D inversion putting large (>∼+7%) anomalies between 40 and 70 km depth above the deeper, more compact anomaly.

Low-velocity anomalies border the Isabella anomaly. From 70 to 120 km depth, these are to the southeast and southwest and down to –4% in amplitude; this arcuate low decreases in amplitude below 120 km. From ∼120 km to 220 km, large (∼–3%) low-wavespeed anomalies are to the north-northwest and north-northeast, with a pronounced anomaly in central California having a geographic extent seemingly unrelated to patterns in the surface topography or geology (the central Great Valley anomaly mentioned above). The area due east of the Isabella anomaly (under the eastern Sierra) remains slightly (0 to –1%) to significantly (0 to –5%) lower in wavespeed on average from 50 to 150 km depth, with more pronounced anomalies in the inversions starting from surface waves. Tests with synthetic structures show that these low-wavespeed anomalies are not artifacts created by the presence of the Isabella anomaly.

The top of the Isabella anomaly is poorly constrained owing to the few raypaths through it and the exceptionally variable wavespeeds associated with the Great Valley sediments and underlying high-wavespeed basement (e.g., profile B, “Show rays” checked, Fig. 12I). All but the “1-D start” inversion (which is questionable as it smears low wavespeeds from sedimentary fill deep into the earth) bring the top of the Isabella anomaly above 70 km, merging with high-wavespeed crust to the west of the anomaly present in the various surface-wave tomographies. Although the geometry of the body below ∼100 km is consistent between inversions (e.g., bottom panels, Fig. 13), the overall geometric shape of the anomaly in an east-west section varies between a drip under the crust (e.g., Moschetti free inversion, thin dashed line, Fig. 13) to a near uniform slab of material bending steeply down from the lower crust (e.g., Moschetti fixed-free inversion, heavy dashed line, Fig. 13). These different images carry different connotations of the mechanism(s) producing the Isabella anomaly; because the different images are from inversions with similar overall fits to the teleseismic data, there is no way to identify the processes creating the Isabella anomaly from the appearance of the anomaly in cross section.

Beneath the western Sierran foothills, the Sierran Foothills anomaly extends discontinuously southward from ∼39.5°N to 36.5°N in the 1-D–based inversions. This anomaly is stronger and/or extends somewhat farther east south of 38°N. This parallels much, but not all, of the Cretaceous Sierran batholith to the east (Fig. 2). Resolution to the west is poor owing to the small number of stations in the Central Valley and the great thickness of sediment there. Surface-wave models all indicate that high wavespeeds at 20 and 40 km depth extend southwest from the teleseismically observed anomaly to the coast (Figs. 11If–11Im). Whether this reflects a single large anomaly or a combination of two or more anomalies requires careful examination of the data.

The inversions starting from a 1-D structure provide insight into the overall magnitude of this anomaly. In cross section, this shallow structure appears to be connected to the Isabella anomaly in the south and possibly the Redding anomaly in the north (Fig. 15, profiles F and V; Fig. 12I). This body has locally exceptionally high wavespeeds (+6%–8%) in the 1-D inversions at depths down to 40 km (Figs. 11Ie and 15). It continues with decreased magnitude (+2%–4%) to 70 km and is absent by 120 km depth. The Foothills anomaly also displays considerable variation in amplitude along its strike. A synthetic test suggests that a uniform anomaly might appear to vary by ∼40% in amplitude due to uneven coverage by the station-event pairs available to us (Fig. 15C). Even so, some variation seems likely to be real, especially the higher wavespeeds to the south and the lower wavespeeds near 38.5°N (∼250 km on Fig. 15).

Specific synthetic tests undertaken to learn the effects of vertical smearing and the fraction of the anomaly recovered suggest that the Foothills anomaly could be entirely within the crust and even higher in amplitude than shown in Figures 11Id–11Im. As illustrated in Figure 10D, a velocity anomaly at 40 km depth could have its amplitude reduced by about half but its depth extent expanded from 0 to ∼100 km depth. A more specific synthetic test starting with an anomaly of +5% in nodes at depths of 20 and 40 km yields somewhat similar results (Fig. 15). The peak values of the recovered anomaly are under 60% of the starting value, with most of the anomaly at or under 40% of the starting value. The anomaly smears preferentially downward into the 70 km layer. Vertically integrating the result from 0 to 95 km (i.e., nodes at 0, 20, 40, and 70 km) only recovers ∼70% of the total amplitude of the original anomaly; integration to 120 km improves this to 80%–85%. Comparing the anomalies from the inverted synthetic structure with the observed tomographic image (Fig. 15) suggests that the Foothills anomaly need not extend below 40 km depth, but concentrating the anomaly into these depths probably requires wavespeed anomalies 50%–100% higher than those present in our inversion (i.e., ∼10%–15% fast). Such laterally high wavespeed bodies are found in the surface-wave models (Figs. 11Ia–11Ic and 11If–11Im).

The surface-wave models and inversions derived from them are far more complex in the region of the Foothills anomaly. In general, these models compress much of the high-wavespeed material into the 0 and 20 km depth nodes, often leaving the lower crust (40 km nodes) with a near or below average wavespeed but then introducing a weaker high wavespeed body in the upper mantle, especially toward the south (e.g., profiles F and V, Figs. 12Ic–12Il). These upper mantle (∼70 km deep) anomalies are only present south of 38°N to the Isabella anomaly and north from ∼39° N to the Redding anomaly. Their amplitude varies from ∼+2 to +5%, depending on the inversion. Unlike high-wavespeed anomalies at 20 and 40 km, which extend to the coast, these deeper anomalies are more strongly focused in the foothills, particularly south of 38°N.

A major complication in interpreting the Foothills anomaly near 40 km depth is the variation in the Moho depth across the Sierra and into the Central Valley. The location of the Moho under the Sierran foothills has been uncertain for some time (e.g., Bolt and Gutdeutsch, 1982; Carder, 1973; Savage et al., 1994), and receiver function studies have suggested that the amplitude of a Ps conversion from the Moho is small or even nonexistent (e.g., Zandt et al., 2004). However, both refraction (Fliedner et al., 1996; Ruppert et al., 1998) and receiver function analyses (Frassetto et al., 2011) indicate that the Moho is deeper under the foothills than to the east or west. The Moho under the Central Valley is inferred to be only ∼25–30 km deep (Fliedner et al., 2000; Holbrook and Mooney, 1987; Miller and Mooney, 1994; Ruppert et al., 1998), while to the east, the Moho is inferred to be ∼35 km deep under the eastern Sierra (Fliedner et al., 1996; Frassetto et al., 2011; Jones and Phinney, 1998; Jones et al., 1994; Ruppert et al., 1998). The wavespeed variations from the travel-time inversion must first be corrected for these variations before lateral variations in wavespeed above and below the Moho can be discerned.

We extend our earlier calculations relating lateral changes in Moho depth (Fig. 6) to compute the variation in crustal and upper mantle wavespeeds from the overall variation at 40 km (reflecting wavespeed variations between 30 and 55 km) while accounting for variation in Moho depths by assuming that the time for a P wave to vertically traverse this layer must remain the same for both flat Moho and laterally varying Moho cases. We assume a 6.35 km/s lower crust under the eastern Sierra (Fliedner et al., 1996; Jones and Phinney, 1998; Ruppert et al., 1998) above mantle averaging 7.8 km/s (somewhat higher than the P_n wavespeeds reported in the area by Fliedner et al., 1996; Jones et al., 1994; Ruppert et al., 1998). If we assume that the variations (in percent) are the same for the crust and mantle, we can then estimate the wavespeeds between 30 and 55 km across the Sierra (Figs. 16A and 16B). The wavespeeds are exceptionally high in the mantle under the foothills, so an alternative shown in Fig. 16C is to limit the mantle to 8.4 km/s and adjust the crust accordingly. These calculations are carried into the Central Valley, where our computed estimate of 8.06 km/s is very close to the 8.1 km/s upper mantle wavespeed inferred by Holbrook and Mooney (1987).

Results of our calculations indicate that the lower crust and upper mantle under the Sierran foothills is distinct from the high-wavespeed crust and upper mantle of the Central Valley (Fig. 16). Even allowing for exceptionally high wavespeed mantle, the lower crust must have a higher wavespeed than areas to the east or west. Such high wavespeeds have been hinted at from seismic experiments within the foothills (Miller and Mooney, 1994; Spieth et al., 1981) and might bring the P-refraction Moho (Fliedner et al., 1996; Ruppert et al., 1998) more in line with the receiver function Moho (Frassetto et al., 2011). The distinct identity of the Sierran foothills anomaly is further reinforced by the coincidence of the higher wavespeed material at 70 km depth aligning with the deeper (>45 km) Moho inferred by Frassetto et al. (2011) (e.g., check Moho box in Fig. 11I and examine 70-km-deep slices).

Any extension of the Foothills anomaly north of ∼38.5°N is complicated by the uncertain Moho geometry and decrease in anomaly amplitude in that region. North from 38.5°N, the weak, deep Moho seen to the south weakens further, and a shallower Ps conversion increases in amplitude and depth to the degree that this shallower conversion is chosen as the Moho by Frassetto et al. (2011). If the Moho does shallow in this manner, the somewhat higher wavespeeds seen at the 40 km layer would reflect the shallower Moho and would not require such extreme high wavespeeds as inferred to the south. If the deeper conversion reflects a continued presence of crustal material, though, then this anomaly could continue north to ∼39°N, where it would seem to deviate to the west under the Central Valley.

Between ∼150 and 250 km below the Foothills anomaly is a broad zone more than 1.5% lower than surrounding mantle. This body is roughly circular in map view at 220 km and fills the area between the Gorda slab and the Isabella anomaly.

The Redding–Gorda Slab anomaly is the prominent high-velocity body at and below depths of 40 km extending from the northern end of the Great Valley (Figs. 11Id–11Im). The Redding–Gorda Slab anomaly as defined here exists between depths of 40 km and 220 km. It has a velocity perturbation exceeding +7% in the southern part of the anomaly at depths of 40–120 km. In most inversions, the magnitude decreases to ∼5% farther downdip. The southern edge of the downgoing Gorda slab is well defined at ∼39.5°N at a longitude of 121.5°W. The slab strikes nearly due north. Sections S1–S3 (Fig. 12I) are oriented roughly parallel to the slab edge and so mildly understate the true dip; section E (Fig. 12I) parallels the trench.

Although the bulk of the anomaly north of ∼39.5°N seems geometrically to be consistent with a subducting slab, there are some oddities to its geometry and amplitude that suggest that the southern part of the anomaly is more complex. First, the amplitude of the anomaly between ∼30 and 100 km depth in the south is 6%–8.5% fast, both shallower than equivalent parts of the slab to the north and higher in wavespeed (Figs. 11Id–11Im; section E, Fig. 12I). Synthetic tests indicate this is not an artifact of the data set or inversion. Second, the southern edge of the anomaly seems inconsistent with the geometry expected for this slab. The sections parallel to the perceived southern edge of the slab (sections S2–S3, Fig. 12I) show that at the southern edge (section S2), the anomaly is quite strong at shallow depths but decreases more rapidly downdip than areas to the north (section S3) despite the better ray coverage to the south. The next slice to the south, which is south of the perceived slab edge (section S1, Fig. 12I), shows a high-velocity body only at shallow depth. One might expect the opposite behavior given the northward migration of the Mendocino triple junction: the slab should be more intact farther downdip along the southern edge (e.g., fig. 13 of Severinghaus and Atwater, 1990). For these reasons, we separate the unusually high wavespeed southern part of the anomaly and term it the Redding anomaly, as this is essentially the same body identified by Benz and Zandt (1993).

When considered on the north-south profile (section E, Fig. 12I), the Redding anomaly at 122.2°W appears to be a distinctive feature of the upper part of the subducting slab. The bottom of the high-wavespeed slab between 120 and 170 km depth is uniform from south to north in amplitude and dips gently northward; above 120 km, however, the wavespeed structure is quite different. This indicates that the cause of this feature either arises within the slab or, more likely, is generated from the upper plate.

A smaller high-velocity body at depths from 0 to 40 km lies under the general area of Sutters Buttes near 39.25°N 121.8°W and seems to be connected to the northern edge of the Sierran Foothills anomaly above 70 km (section V, Fig. 12I). To the north and at depths of 40 and especially 70 km, this anomaly seems to be an extension of the Redding anomaly. Although it is possible another branch of high-wavespeed material could extend southwest from this body under the deep sediments masking the lower crust in that direction, the termination to the south and south-southeast is well controlled. The potential significance of these connections is discussed later.

To the east of the Gorda slab and Redding anomaly are pronounced low-wavespeed regions up to 5% slow at depths of 40–120 km. With a western edge roughly 10–30 km east of the arc volcanoes at Mount Shasta and Lassen Peak, this Gorda backarc anomaly is probably a zone of hot mantle and partial melt responsible for volcanism and high heat flow in this region. The southern edge of this anomaly coincides almost perfectly with the projected southern edge of the Gorda slab.

Although the density of seismometers from the SNEP deployment is greatest in the Sierra, surprisingly large ambiguities remain in the amplitude and exact depth of anomalies presumably related to the mechanism(s) supporting the range. All of the inversions tend to reveal higher wavespeed crust and uppermost mantle to the west than to the east within the Sierra, and they generally tend to show laterally relatively higher wavespeed material overlying laterally lower wavespeed material in the eastern Sierra. The surface-wave models uniformly contain a high-wavespeed uppermost crust in the Sierra (>∼3% above average; Figs. 11Ia–11Ic and 11If–11Im, 0 km) that reflects the crystalline batholith contrasting with sedimentary and more fractured bedrock elsewhere. These high wavespeeds are replaced with depth by low-wavespeed material such that by 40 km depth, the eastern Sierra overlies material ∼3%–4% slow with respect to average and ∼6% slower than the foothills to the west. This is more striking because the Moho is thought to be 5–15 km shallower under the eastern Sierra than the western Sierra (e.g., Frassetto et al., 2011), which should create a 3%–15% wavespeed high under the eastern Sierra in this inversion. A similar pattern is discernable in the inversions starting from a 1-D model (Figs. 11Ie and 12Ib, profiles B and Y) but isn’t as well defined. This relatively low wavespeed material seen in nearly every inversion at 40 km depth could reflect a combination of the low wavespeeds of the deep, silica-rich batholith relative to more mafic lower crust elsewhere along with low-wavespeed uppermost mantle (e.g., Ducea and Saleeby, 1996; Fliedner et al., 1996; Jones and Phinney, 1998; Ruppert et al., 1998).

Greater differences arise in the extent of the low-wavespeed material both vertically and horizontally. The Moschetti-based inversions carry low wavespeeds to the western edge of the Sierra, while the other inversions tend to limit low wavespeeds to areas under the eastern Sierra (Fig. 12I, profile Y; Fig. 17). The difference is compensated by the creation of a high-wavespeed body in the upper mantle beneath the western Sierran foothills in the Moschetti-based inversions, thus preserving the overall pattern of delay times observed.

North-south differences in the geometry of the low-wavespeed anomaly are also present. The 1-D starting models carry the low at 40 km southward through the southernmost Sierra, thus crossing from the high Sierra to much lower terrain (compare topography with 1-D start at 40 km in Fig. 11Ie); however, the surface-wave–based models uniformly increase wavespeeds in this depth range to the south as topography lowers, suggesting that this depth range might be providing the support for the higher elevations of the Sierra. From this and the generally poorer station coverage in the south, we prefer the surface-wave–based solutions.

Although the seismological structures in this region are now well defined, the interpretation of these anomalies in terms of the thermal, petrological, structural, and geochemical properties of the rocks within these anomalies is far less clear. We first examine certain geometrical aspects of parts of the anomalies to help clarify forms that might be related to the evolution of these bodies. We then consider broader geologic and other geophysical constraints to derive a plausible geologic and geodynamic interpretation of the modern seismological structure.

The clearest result of this work is that no high-wavespeed material remains in the uppermost mantle under the eastern Sierra Nevada (profile C, Fig. 12I). The details of the geometry of the low-wavespeed material that is present remain murkier. Under the Sierran crest, all of the inversions place relatively low wavespeed material below ∼10 km depth, but the surface-wave–based solutions put the lowest wavespeed material under some of the highest terrain (∼36.15°N to 38.25°N) in the 40-km-deep layer and contain an increase in wavespeed to the south of 36°N coincident with a decrease in elevation of the range (compare topography layer with 40 km layers in Figs. 11If–11Im). This property is present both in the original surface-wave models and in all of the inversions derived from them; consistency with the overall teleseismic delay times is maintained by lower wavespeeds in the crust of the southernmost Sierra, presumably reflecting a combination of sedimentary rock in some of the basins in this area and emplacement of oceanic-affinity Rand Schist under the Tehachapi Mountains (e.g., Malin et al., 1995). We thus suggest that topographic variations along the strike of the eastern Sierra are being supported by density variations compatible with the observed seismic wavespeed variations between ∼30 and 55 km depth. Geologically these low wavespeeds and densities probably reflect a combination of silica-rich lower crust and hot lower crust and uppermost mantle.

The relation of this low-wavespeed material to the higher wavespeed Foothills anomaly under the western Sierra is far more ambiguous both because of resolution limits at shallow depths and ambiguity introduced by uncertainty in the depth of the Moho. Most provocatively, Figure 17B shows that in the Moschetti fixed-free inversion, the west edge of the low-wavespeed material could be entering a wedge opening created as high-wavespeed material apparently delaminates. This geometry is not seen in the unconstrained inversions (e.g., Figs. 17A and 17C) where the low-wavespeed material only extends under areas above ∼1 km in elevation and seems to either abut or lie under the high-wavespeed material. In the fixed-free inversions, the wedge is smaller and extends less to the west than in the fixed inversion (compare e, i, and l in profile Y, Fig. 12I, with d, h, and k), suggesting that the wedge is probably an artifact reflecting either errors in the surface-wave structure (perhaps related to the profound gradient in the area) or in the conversion of S-wavespeeds to P-wavespeeds, where the P anomalies are being underestimated (v_p/v_s used was too low). Because the Brocher (2005) relationship tends toward somewhat lower v_p/v_s ratios at these wavespeeds than found by Christensen (1996) for mafic rocks likely present in the high-wavespeed body, we tend to favor the unconstrained inversion results and the more unimodal interpretation of the Foothills anomaly.

As noted above, the Foothills anomaly’s high-wavespeed material is most plausibly located above a Moho lying somewhere between ∼40 and ∼55 km depth, though some higher than average wavespeeds could extend into the mantle. If averaged over 50 km, the anomaly south of 38°N would be 5% high and higher, while averaged anomaly magnitudes are only ∼2.5% and 3.5% high from ∼38° to 39°N (Figs. 11Id–11Im); if the anomaly is thinner, the magnitudes would increase proportionally. Some of this difference could reflect the data coverage (e.g., Fig. 15 synthetics), but the bulk of this variation is probably real.

The anomaly generally underlies the more mafic western part of the Sierran batholith and the Foothills Metamorphic belt (Fig. 2), suggesting some relationship to the surface geology. But the northern and southern terminations of the anomaly do not coincide with the ends of geologic structures in the foothills, indicating that this is not simply a product of the observed surface geology. The western edge is, as noted above, tied into the poorly constrained Moho geometry under the eastern Central Valley.

Depths relevant to the geology of this anomaly were sampled by Miocene volcanics of the Kings River and San Joaquin volcanic fields (Fig. 11I, xenolith layer) that presently lie at the eastern edge of the Sierran foothills anomaly. The vertical section of lithosphere inferred from these xenoliths (Ducea and Saleeby, 1996, 1998) contains considerable garnet-rich eclogites from ∼40–100 km depth, with garnet peridotites dominating below. All of these rocks were unusually cold. Such eclogitic lithologies would have the high wavespeeds observed here, as would cold mantle peridotites.

Two observations help illuminate the origin and evolution of this body. First, the modern anomaly does not extend as deeply as would be expected were the vertical section sampled by 8–12 Ma volcanics still in place. Second, the late Cenozoic topographic changes above this anomaly are mostly upwards, with the chief exception being at the south end, adjacent to the Isabella anomaly. From ∼39°N to 37.25°N, Cenozoic cover rocks (classically called the Superjacent Series) are exposed along the western edge of the Sierra over a distance of ∼20 km west from the last exposures of pre-Cenozoic rocks (Figs. 2 and 13). This band of younger rock narrows dramatically south of 37.25°N and is absent from 37°N to 36°N, where these rocks reappear within the Kern arch. Saleeby and Foster (2004) interpreted the embayed topography and absence of Cenozoic rocks as indicating subsidence from 36° to 37°N. These observations suggest that the Foothills anomaly has been thinning north of ∼37°N.

The one area where the Foothills anomaly underlies subsiding crystalline rock is south of 37°N. The largest amplitude part of the anomaly is near 37°N, and the anomaly decreases in amplitude to the south (Figs. 14 and 15). Thus the center of the subsided western edge of the Sierra lies between the Isabella anomaly and the southernmost high of the Foothills anomaly. There is less in the way of high-wavespeed material under the center of subsidence near 36.5°N. If both Isabella and this southern high were equal loads attached to the crust above, we would expect the region of subsidence to include both anomalies instead of simply lying between them. We posit that the bulk of the loading is from the Isabella anomaly, but that its connection to the crust at its eastern side is dominantly through higher viscosity material within the high-wavespeed bridge seeming to connect the Isabella anomaly to the Sierran foothills anomaly (profiles V and F, Fig. 12I; Fig. 14). The development of this asymmetry may be responsible for the complex topographic evolution of the Sierra to the south of 36°N (Mahéo et al., 2009; Saleeby et al., 2013). Such a loading pattern is consistent with that proposed by Le Pourheit and Saleeby (2013).

Previous workers have suggested various origins for the Isabella anomaly, some of which depend on the present shape of the anomaly. In particular, an east dip to this body might suggest that it is either a remnant of the slab or comprised of pieces of lithosphere descending vertically from points of origin to the east. A delamination origin for the Isabella anomaly might also be expected to produce an east-dipping body, though this is less clearly required by this hypothesis. Various Rayleigh-Taylor instability models would seem most consistent with a vertical body, although the presence of vertical gradients in horizontal motion in the mantle (such as a mantle wind, e.g., Zandt, 2003) could produce a tilted body.

Most earlier work found the anomaly to be steep sided or fully vertical (Benz and Zandt, 1993; Biasi and Humphreys, 1992; Jones et al., 1994; Zandt, 2003), but Boyd et al. (2004) found a 50°–60° east dip to the anomaly from body-wave tomography, and a shallower 30° dip was inferred from surface-wave tomography by Yang and Forsyth (2006). As was clear from the discussion above, we find that there is a plunge of ∼60°–70° to the east (Fig. 13), a result that at least below 100 km depth and especially on the eastern side of the anomaly is robust. We do not find the body extending as far east as Yang and Forsyth (2006); the farthest east this anomaly gets when ≥2.5% fast is ∼118.4°W.

Another aspect of the Isabella anomaly is its geometry at the top (Fig. 13, profile B; Fig. 12I). A fossil slab or delamination would presumably look very similar to the Moschetti Fixed solution (Fig. 12Id), while a Rayleigh-Taylor drip might look more like the Moschetti Free or Gilbert Free solution (Figs. 12Ic and 12Ij). These differences are unresolvable with the data used here (as checking the “Show Rays” box on Fig. 12I demonstrates).

A third aspect of interest with the Isabella anomaly is its possible connection to the Foothills anomaly, which terminates to the south near the Isabella anomaly. A profile through the two anomalies shows wavespeeds in excess of 2.5% fast connecting the two anomalies in most inversions, though the Isabella anomaly appears separate in the Moschetti Free inversion (sections A, F, and V, Fig. 12I; Fig. 14). If these are connected, they could both represent a single body perhaps flexing or flowing downward to the south. If separate, either Isabella is detaching from the Foothills anomaly or these are unrelated bodies. Unfortunately, synthetic tests with structures having the two bodies connected or unconnected do not yield clearly different results (left and right columns, Fig. 14). In particular, the downward smearing of the crustal Foothills anomaly tends to generate the appearance of a connection even if one does not exist.

Although the Isabella anomaly has been connected to Quaternary subsidence creating the Tulare Lake basin and sedimentary onlap onto Mesozoic basement rocks of the southern Sierra (Saleeby and Foster, 2004), its position is not centered on the greatest total subsidence. The western part of the Isabella anomaly does directly underlie Tulare Lake and the deepest structural levels of the 0.615 Ma Corcoran clay (Saleeby and Foster, 2004), but the highest wavespeeds are farther east. This could simply reflect the interaction of a load from the Isabella anomaly with the general westward thickening of sediments in the Great Valley such that the subsidence unique to this latitude is being generated by the Isabella anomaly (Saleeby and Foster, 2004), and so this geometry appears consistent with the Isabella anomaly being an antibuoyant body that has grown over the past ca. 4 Ma. The eastern part of the Isabella anomaly underlies the transition from the embayed topography described by Saleeby and Foster (2004) to the north from uplifting sediments in the Kern arch to the south (Cecil et al., 2014; Mahéo et al., 2009; Saleeby et al., 2013). Such a position is unexpected, if the load of the Isabella anomaly is being directly applied to the crust above it, as was discussed above with respect to the Foothills anomaly.

To summarize, the overall picture of the Isabella anomaly is of a 60°–70° plunging body exceeding 4% high wavespeed relative to model average from ∼100 to ∼220 km depth, becoming less pronounced with a slight easterly dip below 220 km. It is probably denser than surrounding mantle. It lies below and is perhaps connected to a shallower crustal anomaly. Its northern edge is near the southern edge of the Foothills anomaly, with which it might or might not be connected. Low-wavespeed material mostly encircles the Isabella anomaly, with pronounced lows from Moho to ∼100 km depth on the southwest to southeast, and below that depth on the north to northeast side.

A fossil slab origin for the Isabella anomaly as advocated by Wang et al. (2013) seems implausible. For this slab to travel with the Pacific plate from its original position relative to North America, Pikser et al. (2012) required the slab to be neutrally buoyant. They did not explicitly determine the degree to which warming of the slab would lead to its erosion and deformation and overstated its likely thickness as 70 km, when the bulk of the material was only 3–4 million years old at subduction (Severinghaus and Atwater, 1990) and so only ∼30–40 km thick (e.g., Bohannon and Parsons, 1995). A thermal model considering the exposure of the slab to the asthenosphere (Bohannon and Parsons, 1995) found that the slab would lack coherence once it passed below the continental lithosphere. Neutral buoyancy would mean that the positioning of the Isabella anomaly under subsidence in the southern San Joaquin Valley would be purely coincidental and would require that the oceanic crust was removed as at the present depths, such crust would be converted to dense eclogite. However, removing the crust would remove nearly the entire thermal anomaly (e.g., Bohannon and Parsons, 1995), meaning that the large seismic anomaly is entirely warm but depleted oceanic lithosphere. Adding some antibuoyancy as proposed by Wang et al. (2013) opens the question of why the slab retains a 60°–70° dip and why it remains so far from the paleotrench. Connection with the Pacific plate also demands a 100–150-km-long décollement somewhere within the high-wavespeed material seen to the west of the Isabella anomaly in all the surface-wave inversions. Apparent connections to the Foothills anomaly would also be fortuitous. Additionally, coastal volcanism has long been considered to be a reflection of emplacement of hot asthenosphere into the slab gap created as the Monterey plate broke; these volcanics do not permit the slab to remain under the Coast Ranges, let alone farther inland, after the Monterey plate was captured by the Pacific plate (Dickinson, 1997; Wilson et al., 2005); recent papers advocating a deeper microplate fail to explain how volcanics were generated and erupted through a coherent Monterey microplate and cold forearc. While this has been an interesting rebirth of an old interpretation (Benz and Zandt, 1993), we find it a poorer explanation than some alternatives.

Delamination of continental lithosphere and in particular mafic batholithic residuum (“arclogite”) has been suggested by Le Pourhiet et al. (2006) from 2-D numerical modeling, although the numerical models exhibit a style of delamination that incorporates considerable deformation. Although the geometry of the body in profile A (Figs. 12I and 13) tempts acceptance of true peeling back of lower lithosphere (delamination senso stricto) (Bird, 1979; Göğüş and Pysklywec, 2008), there are difficulties. If it is a true delaminating body, then we should be able to roll it back into position under its source terrain to the east, which would place it across the Sierra and under Death Valley to the east. If we also restored crustal extension across this region, its paleoposition could have nearly crossed the Basin and Range in the middle Tertiary. While this might violate some suggestions that some mantle lithosphere remains in this footprint (e.g., Schulte-Pelkum et al., 2011), a greater issue is explaining why similar material does not remain to the north of ∼36.5°N (the northern edge of the Isabella anomaly). As noted above, the lithospheric column in the Foothills anomaly appears to be thinner than that represented in the Miocene xenolith section and so cannot represent similarly displaced material. A strictly delamination origin for the Isabella anomaly seems unlikely.

A greater possibility for conforming with observation is provided by allowing deformation within the Isabella anomaly, which in fact is suggested by some internal deformation within the “delaminating body” of Le Pourhiet et al. (2006) and Saleeby et al. (2012). However, nearly any planar 2-D reconstruction (including Rayleigh-Taylor instabilities) will produce the same contradiction: if this body is restored to some original lithospheric geometry in a 2-D section, the resulting edges to the body on the sides paralleling the section will be unexplained. (A radial 2-D solution, such as used by Elkins-Tanton [2007] encounters similar difficulties with the lack of axial symmetry to a lithospheric source.) In essence, there isn’t an obvious scar with the same geometry as this body. There are three solutions assuming this is from a lithospheric source: (1) the sides of the scar have plummeted to depths unobserved by us; (2) extensional deformation has altered the scar to make it unrecognizable; or (3) material has flowed in three dimensions. We cannot explicitly disprove foundering of adjacent lithosphere but view it as an ad hoc explanation. Extensional deformation (and associated magmatic modification of the lithosphere) requires mass balancing beyond the scope of this paper, but the fairly uniform lithospheric structure of the eastern Sierra seems difficult to recover from a combination of a limited foundering of lithosphere and extensional attenuation of the sides of a scar. We therefore suggest that some significant 3-D deformation is the most plausible means of creating the Isabella anomaly. We explore some constraints on such an interpretation below.

Interpretation of the Redding anomaly is seriously clouded by the presence of the Gorda slab and that this anomaly is largely outside the denser SNEP deployment. Despite this, two lines of evidence indicate that this body might be related to the evolution of the Sierra Nevada. First, the high wavespeeds seen near Redding are both thicker and higher amplitude than within the slab to the north as far as central Washington (Burdick et al., 2008; Roth et al., 2008; Schmandt and Humphreys, 2010; Xue and Allen, 2010; Yang et al., 2008); therefore, it seems curious that the largest anomalies would occur at the edge of the slab presumably exposed to asthenospheric upwelling as the slab moves north. These high wavespeeds are limited to the uppermost part of the slab edge (section S2, Fig. 12I). Second, high-wavespeed material extends well to the south of the edge of the Gorda plate and connects, in a fashion, to the Foothills anomaly (sections E and V, Figs. 12I; Fig. 15).

The high-wavespeed bodies between the north end of the anomaly under the Sierran foothills and the Redding anomaly are enigmatic and interpretable in multiple ways. The high-wavespeed body seen trending due south from the Redding anomaly could be a piece of imbricated slab against the base of North America. It could be some chilled piece of the upper plate now exposed to asthenosphere. Given the increase in amplitude under Sutter Buttes, an extinct Quaternary stratovolcano, it could reflect some depletion of the upper mantle associated with magmatism. Finally, it could reflect a path by which mafic crust or mantle lithosphere is removed from parts of the northern Sierra.

Although a definitive interpretation is not possible, several factors point toward some connection with the Sierra. The Foothills anomaly, which we have suggested is antibuoyant and connected to the Isabella anomaly, terminates its north-northwest trend within the Sierran foothills southeast of the Redding anomaly. At depths of 20 and 40 km, this anomaly seems to continue to the west to the vicinity of Sutter Buttes; at depths of 70 km, this appears to be continuous to the north with the Redding anomaly (profile V, Fig. 12I). Our synthetic tests of a direct connection along the strike of the Sierran foothills seem to indicate that there is probably not a connection on this line (Fig. 15). However, the connection from the Redding anomaly to the high near Sutter Buttes as shown on profile V (Fig. 12I) is unlikely to be an artifact like those seen in Fig. 15, which were generated by upward or downward smearing of anomalies at dissimilar depths. The connection from Sutter Buttes to the north end of the Sierran foothills anomaly, while having smaller amplitudes of +2%–2.5%, is also unlikely to be an artifact from smearing.

The northern end of the Foothills anomaly is itself quite curious; unlike the rest of the anomaly to the south, which occurs in areas with a deep (40–55 km), low-amplitude Moho as determined from receiver functions, the Moho inferred by Frassetto et al. (2011) in this area is quite shallow. One plausible interpretation of this part of the anomaly is that it represents a Moho ∼10–15 km shallower than areas to the east (e.g., Fig. 6). In fact there are two Ps conversions in this region, a deeper one that weakens to the north and a shallower one that strengthens and deepens to the north. Rather than this part of the Foothills anomaly being distinct from the anomaly to the south, we suggest it is indeed the continuation of the Foothills anomaly. The appearance of a shallower Moho conversion could reflect tectonism associated with the step to the west of this anomaly or, as suggested by Frassetto et al. (2011), a feature reflecting the changing architecture of the Foothills metamorphic belt as it widens, structural dips decrease, and the intensity of Cretaceous plutonism declines. The continued presence of anti-buoyant material through this region is supported by our analysis of the isostatic balance in this region (Levandowski et al., 2013).

We consider three different origins for the Redding anomaly: a complication within the slab, a superimposed crustal anomaly, and the entrainment of material from elsewhere. The clockwise rotation of the southernmost part of the Juan de Fuca plate (more informally called the Gorda plate) in the past few million years suggests that there has been considerable internal deformation within the plate (cf. Wilson, 1989). If the subducted slab also rotated in such a manner, some imbrication of the slab might have occurred, potentially near its southern edge. Wilson (1989) suggested that this produced a bent southern edge to the slab (slab edge check box, Fig. 11I) somewhat in line with earlier inferences from gravity (Jachens and Griscom, 1983). Section S1 on Fig. 12I should show the slab, if it extends as far south as shown by Wilson, but it apparently does not. Instead the slab edge is far closer to the projection inferred from local seismicity by Furlong and Schwartz (2004) and Verdonck and Zandt (1994). Either there is an error in the analysis by Wilson (1989) or the material that should have extended the plate farther south has remained farther north as the plate thickened internally, perhaps contributing to the Redding anomaly. In the latter case, though, the absence of equivalent thickening at greater depth seems peculiar; the high velocities distinctive of the Redding anomaly are absent by 220 km and quite possibly could be restricted to depths above 170 km.

A second possibility is that the Redding anomaly is simply the superposition of a crust high-wavespeed anomaly associated with the Great Valley onto the southern end of the Gorda slab; the small gap between these two anomalies might be undetectable in this more poorly sampled part of our model. Such an interpretation is tempting in view of the shallow upper edge of the Redding anomaly (sections D and E, Fig. 12I), but the southward termination of the large anomaly seems abrupt, if the crustal contribution is to continue south, as is usually inferred (e.g., Cady, 1975). There is no indication in the synthetic checkerboard (Fig. 9) that such a major anomaly would vanish as it continued into a more densely sampled part of the inversion. However, the great thickness of sediment on the west side of the valley (Fig. 2) could mask a lower-crustal, high-velocity body on that side of the valley. The Redding anomaly is to the northeast of the high-wavespeed material clearly imaged by the surface-wave analyses at 40 km depth between the Sierra and Pacific Ocean (Figs. 11Ie–11Im).

The third possibility is that some additional material has been entrained in the subduction zone. Certainly the flow field expected in the mantle wedge above the slab could entrain material, but the motion of material from south of the slab is less certain. If the slab is rising, material from within the mantle wedge might be expelled and entrainment above the slab would be difficult, but if the slab were steepening, material would flow into the mantle wedge from the side, making entrainment from the south plausible. Although this suggestion provides some explanation for the curious tendril of high-wavespeed material extending south from the Redding anomaly, it is difficult to demonstrate that such deformation would produce the seismic structure observed. Additional complications would be expected from the northward motion (relative to North America) of the slab edge.

The analysis above suggests that the Isabella anomaly is anti-buoyant and that high-wavespeed material within the Foothills anomaly is connected to it and probably also antibuoyant. In contrast, the eastern Sierra lacks any arclogitic root despite geological inference that it should have been there as recently as 8–12 Ma. We will first consider these anomalies before considering the implications for the Redding anomaly.

Accepting our earlier argument against a slab origin, the material in the Isabella anomaly probably originates in the Sierra Nevada. Removal of lithospheric material from the Great Valley to the north should produce uplift that is not seen. Areas to the west were almost directly on the slab prior to passage of the Mendocino triple junction ca. 20 Ma; the amount of lithosphere created by cooling since then would be small. To the southeast, lithosphere was removed early in the Laramide orogeny and replaced in part with schists with oceanic affinity such as the Rand schist (Malin et al., 1995; Saleeby, 2003). Although the exact lithospheric architecture after this tectonism is unclear, it seems unlikely to produce thick, high-wavespeed material such as that in the Isabella anomaly. This leaves the Sierra Nevada and, just possibly, the Pacific plate as potential sources for this material. Given the plausible connection between the Isabella anomaly and the Sierran Foothills anomaly, the Sierra seems the most likely source for this material.

It is difficult to fully compare the volume of the Isabella anomaly to the volume of plausible source areas. First, the volume depends on the degree of the anomaly chosen to define it (e.g., the 1% or 3% contour). Second, if the material is dominantly cold mantle, it will heat as it descends and it cools surrounding mantle. A third problem is the possibility that the regional mean seismic velocity structure changes with depth such that material with the same potential temperature at two different depths might have a different percent anomaly relative to average structure. Fourth, as noted in discussing the checkerboard tests, teleseismic tomography will smear a compact anomaly into a larger but lower-amplitude anomaly.

To get some measure of the volume of the Isabella anomaly, we have chosen to integrate the percent anomaly over the volume extent of the anomaly. In tests on synthetic structures, such integrals nearly perfectly recover the total anomaly at these depths. The area within the +1% contour at 120, 170, and 220 km depth is 18,000–20,000 km² (from GV 1-D inversion; 18,000–25,000 in the Moschetti Fixed-Free inversion). The volume within the 1% contour from these layers is 2.9–3.2 × 10⁶ km³. The mean anomaly within this volume is 2.1%–2.5% fast relative to mean. Note that we ignore the volume outside the +1% anomaly, and we have not attempted to integrate in some parts of the anomaly from 70 and 270 km node depths; so our estimate should be a minimum reflecting the anomaly between ∼95 and 245 km depth.

An area of ∼17,400 km² lies between 36° and 38°N, east of the high wavespeeds seen at 20 km depth and west of Owens Valley. Areas south of 35.7°N plausibly had most of their Mesozoic lithosphere removed during the Laramide (e.g., Malin et al., 1995; Nadin and Saleeby, 2008; Saleeby, 2003). Restoring the Isabella anomaly to this region, it would be a layer 350–450 km thick at 1% fast or, equivalently, a layer 10% fast only 35–45 km thick. Lherzolite with an average temperature of 600 °C would be ∼4.8% faster and ∼80 kg/m³ denser than the same material at 1300 °C, ignoring any possible melt and anelastic effects (Hacker and Abers, 2004). Converting the Isabella anomaly into a layer of this material would be ∼73–94 km thick. Removal of such a body should produce uplift of 1.8 km or more. Because garnet pyroxenite described by Ducea and Saleeby (1996, 1998) is far denser than the lherzolite for similar P wavespeeds, an even greater load could reside within the Isabella anomaly than calculated here. We defer a more thorough investigation of the possible lithologies in the Isabella anomaly and their implications to a separate paper.

Within these constraints, there appears to be no objection to the Isabella anomaly being the sink for lower lithosphere from a broad region. We continue along these lines by next considering the means by which this material could have arrived in the anomaly. If the Foothills anomaly is in place with respect to the upper crust, then the material in the Isabella anomaly had to flow from north to south. Delamination seems implausible with this geometry, as we would expect such a peeling body to be connected to the lithosphere to the south; instead, as noted previously, it appears that the Isabella anomaly is, on its east side, connected to lithosphere to the north. The shape of the Isabella anomaly is less inconsistent with a delamination peeling from east to west, but, as noted above, there are issues with a resulting lithospheric scar. A delamination progressing from south to north might be consistent with our observations, but this demands a large source area to the south, something we discounted above owing to emplacement of oceanic schists in this region. From this analysis we determine that there must have been internal deformation within the material as it traveled from its source area to the Isabella anomaly. Such internal deformation of presumably cold material has been shown to be physically plausible (Molnar and Jones, 2004) and also emerged in different form in numerical experiments (Le Pourhiet et al., 2006; Saleeby et al., 2012).

How material moved from the eastern Sierra to the Isabella anomaly depends on the history of material in the Foothills anomaly. If the anomaly represents material unmoved since creation in the Mesozoic, either eastern Sierra lower lithosphere went under it or around it (Figs. 18A and 18B). For material to move around the Foothills anomaly, some instability had to nucleate in the foothills near 36°N, and then flow was channeled by viscosity variations. It is hard to see why material from the eastern Sierra near 38° N would follow such a path without descending into the mantle somewhere sooner. Material moving under the Foothills anomaly might be more akin to a radial 2-D Rayleigh-Taylor instability, but the rapid descent of material moving from the eastern Sierra near 40 km depth to 70 km or more under the Foothills anomaly only to move horizontally again to the Isabella anomaly seems physically suspect.

A second path connecting the Isabella anomaly to an eastern Sierra source requires that the Sierran Foothills anomaly be, at least in part, allocthonous. In this construct, material under the eastern Sierra moves to the west, thickening dense material under the foothills as this material thins under the crest (Fig. 18C). This is very much like an asymmetrical Rayleigh-Taylor instability. The second part of the motion is from north to south: the Isabella anomaly initiates as an instability, possibly because of the free southern edge to dense material, and starts to draw in material from the north (Fig. 18D). These motions could be simultaneous such that material from the crest might move southwesterly more than west and then south. The remaining connection between the Isabella and Foothills anomaly would reflect this conduit, perhaps in the process of necking out. As the Isabella anomaly developed, it might have pulled in material from the other sides as well, a process that might have helped develop the low-wavespeed halo on those sides seen at 70 and 120 km depth (Figs. 11Id–11Im). Perhaps in this process some of the oceanic schists from the south were pulled into the southern top of the Isabella anomaly, effectively disconnecting the antibuoyant mass from the Kern arch to the south, permitting the uplift of this feature.

This path through the Foothills anomaly suggests that there have likely been changes in the buoyancy of this material as well. Presumably as this material is pulled into the Isabella anomaly, it starts to thin some distance away from the Isabella anomaly. At this position, the surface should rise up. The lower amplitude of the Foothills anomaly under the raised Tertiary section just north of 38°N (the central Foothills swell of Saleeby et al., 2013) might be such a feature.

Perhaps the Redding anomaly is, in part, material derived from the Sierra. The generally lower elevations of the northern Sierra suggest that less material needed to be removed, and geochemistry of the northern batholith seems to suggest less garnet was present in the residuum (Cecil et al., 2012). Because of the presence of the slab, a mass balance such as we attempted for the Isabella anomaly is difficult and will not be attempted here. The connection from the Sierra to the Redding anomaly is more problematic: there is no drowned topography, as there is to the south, though the width of the Tertiary exposures decreases significantly, and the connection is through a more convoluted path. Perhaps flow is influenced by the counterflows driven by slab descent and upwelling driven by northward motion of the slab edge. These effects are plausible, but we are unable to predict if they produce the geometry of high-wavespeed material observed. Although this is by no means proved, the information we have does not disprove such a connection.

Asthenosphere must replace any dense material removed in this manner. We suggest that this upwelling tends to occur near the Sierran crest or to the east, with the asthenospheric material then moving to the west in trailing the downwelling dense material. If our flowfield has been moving dense material as suggested, then there is a position somewhere between the Redding and Isabella anomalies where these flows diverge. Such a position might further enhance upwelling; this might be responsible for the Quaternary magmatism seen near Mono Lake, including the Long Valley caldera (Figs. 1 and 2). An area of deep crustal seismicity located to the west of the Long Valley caldera (∼37.4°N, 120°W, Ryan et al., 2008) might also reflect increased strain rates in this area of divergence.

P-wave tomography of the Sierra Nevada locates three large, high-wavespeed bodies: the Isabella anomaly, the Redding anomaly, and the Sierran Foothills anomaly. No high-wavespeed material was found under the eastern Sierra, indicating that any mafic residuum is now gone. We find that the Isabella anomaly plunges 60°–70° east from 100 to 220 km depth with peak amplitudes of ∼4%–5% higher wavespeed than average. The western part directly underlies historic Tulare Lake, previously inferred by Saleeby and Foster (2004) to be a locus of Quaternary subsidence. The Redding anomaly resembles a large bulge in the Gorda slab, with higher amplitudes and greater thickness than the slab farther north. The anomaly has peak wavespeeds more than 8% faster than average. The Foothills anomaly lies along the western Sierran foothills between the two deeper anomalies and appears to be distinct from the thin, high- wavespeed crust under the Central Valley just to the west. Unlike earlier studies, we place this body above or near the Moho and estimate its wavespeed anomaly to be more than 2%–5% faster than average, depending on the degree of vertical smearing that has occurred and the Moho depth.

We argue that these three anomalies are related to the convective removal of dense material from beneath the Sierra Nevada. The volume of material in the Isabella anomaly is consistent with the volume of material removed from under the eastern Sierra. We suggest that material originally under the Sierran crest traveled westward and toward the Isabella anomaly, if south of ∼38°N, and the Redding anomaly, if to the north. This interpretation is largely based on the geologic and seismological relations observed near the Isabella anomaly, where complications associated with the subduction of the Gorda plate are absent.

This inference indicates a potential inadequacy of numerical experiments exploring the removal of dense lower lithosphere made to date. The geometry of this inferred flow system lacks either radial or linear symmetry. The downwellings localized at two prominent breaks in structure: Isabella at the southern edge of preserved Mesozoic lithosphere, Redding at the northern edge of the “slab gap,” and the large left-step of pre-Cenozoic structure from the Sierra to the Klamath Mountains. This suggests that profound lithospheric discontinuities might be important in the development of downwellings.

We are indebted to the numerous public and private landowners who permitted us to place seismometers on their land. Conversations with colleagues in the Sierran Drips Continental Dynamics Project and others participating in our occasional workshops helped to shape our views and understanding of the region. We are indebted to G. Pavlis and S. Roecker for their help in getting codes to behave as desired. M. Moschetti shared his velocity model with us. Several figures were prepared with the aid of Generic Mapping Tools (Wessel and Smith, 1998). Two anonymous reviews and guest editor Saleeby helped to hone the text and improve the presentation. This work was supported by National Science Foundation EAR grants 0454535, 0454524, and 0454554 to University of Colorado, University of South Carolina, and University of Arizona, respectively.