Abstract

The putative Pliocene–Quaternary removal of mantle lithosphere from beneath the southern Sierra Nevada region (California, USA) is investigated by the iteration of thermal-mechanical models that incorporate and are tested against a range of data that are geologically observable, including rock uplift and basin subsidence data, structural and compositional data on crustal architecture, and a synthesis of seismic data that image lower crust–upper mantle structure of the region. The primary focus is testing model results with rock uplift and basin subsidence data. The initial state of our models recognizes that (1) the sub–Sierra Nevada batholith mantle lithosphere, including a substantial thickness of eclogitic cumulates that were produced during high magma flux arc activity, termed arclogite, was cooled to a conductive geotherm by amagmatic flat slab subduction at the end of the Cretaceous; and (2) the gravitationally metastable mantle lithosphere was thermally mobilized from beneath in the Neogene by the opening of an underlying slab window. Based on a detailed synthesis of appropriate rheologies of the multiphase system, a preferred class of models correctly predicts (1) the ca. 10 Ma inception of the Sierra Nevada microplate due to a lithospheric separation event along the eastern Sierra Nevada region as a result of the mobilization of the mantle lithosphere as a Rayleigh-Taylor instability; and (2) the subsequent delamination of the arclogite root of the Sierra Nevada batholith that appears to be in progress. Our preferred model also predicts focused rock uplift and basin subsidence resulting from delamination, both of which are anomalous to uplift and subsidence patterns of all other regions of the microplate. The rheology of the Great Valley crust is found to control rock uplift patterns across the Sierra Nevada, and tectonic subsidence in the Tulare Basin of the Great Valley. The Tulare Basin is uniquely situated over the region where the principal residual arclogite root remains attached to batholithic crust. The anomalous rock uplift and tectonic subsidence data are best satisfied by modeling a bulk rheology for the Great Valley crust that is similar to that of the Sierra Nevada batholith. These results are consistent with a recent synthesis of basement core and geophysical data showing that much of the Great Valley basement consists of the western Early Cretaceous zones of the Sierra Nevada batholith. The existence of this batholithic domain within the Great Valley subsurface is also in agreement with recent seismic data that resolve additional residual arclogite root materials along the base of the crust of this region.

INTRODUCTION

The removal of mantle lithosphere from beneath continental crust has been called on to explain a number of tectonic features, including volcanism and anomalous heat flow, upper mantle seismic velocity and gravity anomalies, extensional tectonism, and both positive and negative topographic transients (cf. Bird, 1979; Kay and Kay, 1993; Ducea and Saleeby, 1998a; Jones et al., 2004; Morency and Doin, 2004; Pysklywec and Cruden, 2004; Le Pourhiet et al., 2006). These and a number of other geologic features, when considered together, strongly suggest the post–10 Ma removal of the mantle lithosphere that formed beneath the southern Sierra Nevada batholith (California, USA) in conjunction with Cretaceous high-flux arc magmatism (Ducea and Saleeby, 1998a; Saleeby et al., 2003; Zandt et al., 2004). Paralleling the geophysical and petrologic data in support of this putative event are geomorphic and basin subsidence data that reveal distinct vertical displacement pulses in the Earth’s surface that correlate in time and space with mantle lithosphere removal in the region (Wakabayashi and Sawyer, 2001; Saleeby and Foster, 2004; Stock et al., 2004; Clark et al., 2005; Maheo et al., 2009; Figueroa and Knott, 2010; McPhillips and Brandon, 2010; Hammond et al., 2102). Placing such observations in a physical context, these findings are augmented by a number of thermomechanical modeling studies of mantle lithosphere removal that have both generic and southern Sierra Nevada case-specific pertinence (Jull and Kelemen, 2001; Molnar and Jones, 2004; Morency and Doin, 2004; Pysklywec and Cruden, 2004; Le Pourhiet et al., 2006; Elkins-Tanton, 2007; Göğüş and Pysklywec, 2008). These offer the opportunity to compare surface geologic observations, remote geophysical imaging, and physical modeling in order to gain a fuller understanding of the geometry, physical processes, and geologic consequences of the southern Sierra Nevada mantle lithosphere removal event.

The focus of this paper is the testing of thermomechanical model results against geologic observations of rock uplift and tectonic subsidence. Thermomechanical modeling employed here was presented in an initial form in Le Pourhiet et al. (2006). We pursue these models, with new directions posed by sensitivity tests and derivative questions that were raised in our initial work regarding the influence of crustal rheology on surface displacement patterns. Observations and model predictions reveal transients in such vertical displacements; these transients are referred to as epeirogenic because they reflect relatively large wavelength to amplitude ratio surface displacements resulting primarily from radial force components that evolved during the course of mantle lithosphere removal. The geographic focus is on the southern Sierra Nevada and adjacent Great Valley (Fig. 1). In that the Sierra Nevada and Great Valley form a north-northwest–south-southeast elongate mountainous uplift and adjacent basinal system, modeling of both the principal seismic velocity anomaly and the dynamics of the removal process have focused on transverse (east-northeast–west-southwest) profiles. We diverge from this practice, and strive to consider both the seismic anomaly and the vertical displacements in three dimensions. The need for such a three-dimensional analysis will become evident upon considering both the three-dimensional geometry of the mobilized mantle lithosphere and the map distribution of the resolved transients (cf. Figs. 1 and 2). Experimentation with three-dimensional models that, like our two-dimensional models, entail visco-elasto-plastic rheologies and topographic changes is technically challenging and in progress. The experiments we have performed with our two-dimensional procedures, however, offer insights that can be applied to the observed three-dimensional patterns of uplift, subsidence, crust-mantle structure, and volcanism. An in-depth analysis of the three-dimensional pattern of uplift and subsidence across the southern Sierra Nevada region will be given in Saleeby et al., (2013, this themed issue), hereafter referenced as Part II. The most pertinent results of this analysis are compared here to our model results.

MANTLE LITHOSPHERE REMOVAL IN THE SOUTHERN SIERRA NEVADA REGION

Discovery and Resolution

The profound heterogeneity of the upper mantle beneath the U.S. Cordillera has been documented by seismic studies of the USArray (http://www.usarray.org/; Yang and Forsyth, 2006; Sun and Helmberger, 2010; Schmandt and Humphreys, 2010a, 2010b). A distinct high-velocity anomaly was discovered beneath the southern Sierra Nevada–Great Valley region (Biasi and Humphreys, 1992; Zandt and Carrigan, 1993; Jones et al., 1994). Referred to as the Isabella anomaly, this body is ∼100 × 200 km in diameter, extends from the base of the crust to ∼250 km depth, and has a P-wave speed anomaly (dVp) locally in excess of 5% fast. With the higher resolution inversions based on USArray data (Reeg, 2008; Gilbert et al., 2012; C.H. Jones, 2011, written commun.), the Isabella anomaly is one the more distinctive features beneath the U.S. Cordillera. In contrast to this high-velocity anomaly, the high-elevation axial and eastern Sierra Nevada is underlain at shallow Moho depths (∼35 km) by anomalously low seismic velocity upper mantle (Carder et al., 1970; Carder, 1973; Jones and Phinney, 1998; Ruppert et al., 1998; Fliedner et al., 2000; Zandt et al., 2004; Savage et al., 2003; Frassetto et al., 2011; Gilbert et al., 2012), indicating a lack of mantle lithosphere in this region. Petrologic study of mantle xenoliths derived from beneath the southern Sierra Nevada batholith indicate at least partial replacement of mantle lithosphere by modern asthenosphere within the time interval of 10–3.5 Ma (Ducea and Saleeby, 1996, 1998a, 1998b; Lee et al., 2001; Saleeby et al., 2003). Xenoliths of the subbatholithic mantle lithosphere suite are characterized by abundant eclogitic rocks that formed as cumulates within the batholithic source regime, and coexisting mantle wedge peridotites, the deepest of which equilibrated at ∼125 km. Recognizing the unique geodynamic and petrogenetic significance of the eclogitic rocks, Anderson (2005) suggested the term “arclogite” (adopted herein). Thermobarometry and relative abundance of rock types in the xenolith suites indicate that the arclogites were concentrated between ∼40 km and ∼75 km deep within the mantle wedge section, interlayered with subordinate spinel peridotite, and at greater depth spinel and garnet peridotites dominate the section with subordinate arclogite layers. The modern asthenosphere suite is characterized by a dominance of spinel and subordinate plagioclase peridotites, commonly with melt inclusions, and thermobarometric constraints that yield high potential temperatures, similar to mid-ocean ridges. Rocks of the asthenospheric suite are immediately beneath the base of the crust in the axial to eastern Sierra Nevada region, corroborating the seismic results referenced here.

The geophysical and petrologic data that suggest the removal of mantle lithosphere from beneath the southern Sierra Nevada region and its gathering into the Isabella anomaly are corroborated by a number of surface geologic observations. First was the outpouring of late Neogene basalts in the region (Fig. 1), and in particular a concentration of late Pliocene potassic basalts in the area where mantle lithosphere has been clearly removed (Van Kooten, 1980, 1981; Manley et al., 2000; Farmer et al., 2002). A phase of Pliocene rock uplift and elevation increase along the axial and eastern southern Sierra Nevada was coincident (Unruh, 1991; Wakabayashi and Sawyer, 2001; Stock et al., 2004; Clark et al., 2005). This rock uplift and elevation increase are paralleled to the west, above the principal seismic anomaly, by Pliocene–Quaternary anomalous subsidence in the Tulare Basin that occurs as a subbasin in the southern Great Valley (Saleeby and Foster, 2004; Saleeby et al., 2009; Part II). Coseismic strain patterns and late Cenozoic faulting of the southern Sierra Nevada region reveal extensional tectonism concentrated across the area that appears to have recently lost its mantle lithosphere (Unruh and Hauksson, 2009; Saleeby et al., 2009; Amos et al., 2010; Nadin and Saleeby, 2010; J.R. Unruh, 2012, written commun.). All of these features are readily explained in the context of the models that we present in the following.

Structure of the Isabella Anomaly

Some of the important geometric features of the Isabella anomaly and its vertically projected map relationships with surface geology are shown in Figures 1–3. We focus first on correlation between the geometry of the seismic anomaly, and what we define as two distinct epeirogenic domains of the region. The first domain consists of the Tulare Basin (Saleeby and Foster, 2004), an ∼100-km-long zone of anomalous Pliocene–Quaternary tectonic subsidence located between ∼36°N and ∼37°N within the San Joaquin Basin of the southern Great Valley (Part II, Fig. 1). Along the northeast to southwest perimeter of the Tulare Basin are areas of the Sierra Nevada that exhibit evidence of significant Pliocene–Quaternary rock uplift more clearly than any other areas of the range (Part II). We develop the rationale for designating this region of anomalous rock uplift as the delamination bulge. Comparison of Figures 1 and 2 shows that the principal Pliocene area of the delamination bulge coincides with an area of distinct 4 Ma or younger volcanic rocks, and that the principal Quaternary area of the bulge coincides with an area characterized by a modern anomalous thermal transient. The young volcanic and geothermal areas are interpreted as mature and juvenile thermal responses, respectively, to mantle lithosphere removal.

Lithospheric-scale structure sections of Figure 3 show that the Isabella anomaly is attached to the base of the crust beneath the Tulare Basin, and that it extends eastward and southward deeper into the mantle beneath the region of the delamination bulge. The structure sections are constructed from a synthesis of P-wave tomography, receiver function, and refraction data (Mereu, 1987; Jones and Phinney, 1998; Ruppert et al., 1998; Fliedner et al., 2000; Zandt et al., 2004; Yan and Clayton, 2007; Reeg, 2008; Frassetto et al., 2011; C.H. Jones, 2011, written commun.). Section A is along a near transverse trace, relative to Sierra Nevada–Great Valley convergent margin structure, and shows a high-velocity body attached to the base of the crust and dipping steeply eastward into lower velocity mantle. Along this section dVp is as much as ∼5% fast to ∼125 km depth, and ∼2%–3% fast to ∼250 km depth. Sections B and C rotate progressively southward from the trace of section A through diagonal and longitudinal traces, the latter oriented along the eastern Great Valley. Section B shows the steeply east-dipping anomaly with a core area of dVp >5% fast extending down to ∼75 km, and dVp ∼3%–4% fast extending down to ∼190 km. The deep 3%–4% fast appendage is enveloped to the north by the ∼2%–3% fast appendage of section A. Section C shows a zone of much more extensive lower crustal attachment for the high-velocity mantle lithosphere rocks than sections A or B, with the northern part of the anomaly along C restricted to relatively shallow levels (<100 km). These laterally extensive high-velocity rocks are interpreted as the residual arclogite root of Early Cretaceous high-volume batholithic rocks that are widespread along the western Sierra Nevada foothills and subsurface of the eastern Great Valley (May and Hewitt, 1948; Saleeby, 2007, 2012; Lackey et al., 2012). The high-velocity rocks extend northward as far as ∼38°N; peridotite-rich remnants of the Cretaceous mantle wedge extend farther north to ∼39.5°N (Gilbert et al., 2012; Fig. 4B). North of ∼36.5°N the Isabella anomaly is restricted to relatively shallow levels, the steeply east-dipping appendage fading out northward over a velocity gradient similar to that of the section B to A transition. The depth to Moho in the Figure 3 sections is well constrained by refraction and receiver function data across the Sierra Nevada, and reasonably constrained by refraction data across the Great Valley. All sections show a westward thickening of the crust across the zone of lithosphere detachment, resulting in the greatest crustal thicknesses for the Sierra Nevada along its lowest elevation western margin. All sections also show low Pn with underlying low Vp mantle to the east and southeast of the zone of crustal attachment. These same areas show negative conversions in receiver functions, interpreted to be zones of partial melt in the mid-crust and upper mantle.

The apparent locus of separation of the Isabella anomaly from the base of the crust is designated as the delamination hinge, as denoted by black dots in the Figure 3 sections. An approximate surface trace for the delamination hinge is projected onto the Figures 1 and 2 map surfaces from the Figure 3 relations, and from additional tomographic sections (Reeg, 2008; C.H. Jones, 2011, written commun.). In the analysis that follows, we use these syntheses as the best available approximations for the geometry of the Isabella anomaly, and its surface manifestations.

Dynamic Analyses

In general two mechanisms are considered for the removal of mantle lithosphere from the base of the crust. The most commonly pursued mechanism is convective removal as a Rayleigh-Taylor (RT) instability (cf. Houseman and Molnar, 1997), characterized by viscous drip-like behavior. The RT mechanism is widely cited for the Isabella anomaly (Biasi and Humphreys, 1992; Zandt and Carrigan, 1993; Jones et al., 1994; Saleeby et al., 2003; Zandt et al., 2004; Elkins-Tanton, 2007). The second mechanism is delamination, or the peeling away of the mantle lithosphere from the Moho with the infill of asthenosphere (Le Pourhiet et al., 2006). Several features have been suggested as being diagnostic of which of these two mechanisms has operated, or dominated in a given case (cf. Göğüş and Pysklywec, 2008). For example, delamination should peel away the mantle lithosphere from the lower crust, whereas RT behavior in most cases should leave an attenuated residual lens of mantle lithosphere beneath the crust. At surface levels epeirogenic transients should be symmetrically situated over an RT instability, whereas for delamination such transients evolve asymmetrically across the margin of the instability.

The geometry of the Isabella anomaly is most consistent with that of a mantle lithosphere slab that peeled off of the Moho in both east to west and south to north directions (Fig. 3). Such asymmetry is atypical of RT behavior. Furthermore, Pn values immediately to the east and south of the anomaly reveal low-velocity asthenospheric mantle immediately beneath the Moho, with no sign of an attenuated residual layer of mantle lithosphere. Similarly, mantle xenolith data indicate that where asthenosphere has ascended to shallow levels beneath the crust, no residual lens of mantle lithosphere remains (Ducea and Saleeby, 1996, 1998a). In the analysis of epeirogenic transients (presented elsewhere; Part II), highly asymmetric behavior is resolved as well. These features suggest a strong role for delamination-type behavior in the southern Sierra Nevada case. Nevertheless, our modeling suggests behavioral patterns entailing essential features of both the RT and delamination mechanisms, as discussed in the following.

In the following analysis, we adopt and expand upon the models of Le Pourhiet et al. (2006) that predict interplay between mantle lithosphere RT behavior and arclogite root delamination. Upon the imposing of a basal thermal perturbation and modest horizontal extensional strain, the load of the arclogite root induces a perturbation that instigates convective mobilization of the lower mantle lithosphere. Such RT behavior leads to the delamination of the arclogite root, which is facilitated by a low-viscosity channel that develops along the lower felsic crust. Strength drops of structural character within the lower crust and upper mantle are critical predictions of our model, and facilitate geologically rapid delamination (<10 m.y.). The mode of delamination that is predicted by our modeling may be considered a form of RT behavior, but complicated by multiple mantle lithosphere phases and structural strength drops, which are rarely incorporated into numerical models of RT behavior.

MODEL SETUP AND RESULTS

Setup

The numerical models in Le Pourhiet et al. (2006), developed further here, use a thermomechanical approach whereby well-constrained geological parameters such as initial geotherm, structural geometry, and elastic and brittle rheologies are fixed, and the most poorly constrained creep parameters are varied and evaluated in terms of their influence on mantle instability behavior and the geological resultants. The code employed models topography, visco-elasto-plastic temperature-dependent rheologies, and temperature- and composition-dependent density (Cundall and Board, 1988; Poliakov et al., 1993; Le Pourhiet et al., 2004). The model space is along a transverse section (Fig. 1) extending normal to the Cretaceous convergent margin of central California. The regional geologic structure that is used to help constrain the initial conditions is summarized in Figure 4, and the initial and boundary conditions are summarized in Figure 5. Figure 5 also shows model notations, and Figure 6 shows the results of our preferred model (model 4).

In our models the gravitationally unstable sub–Sierra Nevada batholith mantle lithosphere is assumed to have remained intact and attached to the base of the crust throughout early Cenozoic time, supported beneath by flat-slab subduction of the Farallon plate (Fig. 4). We use the term “arclogite root” to designate that part of the subbatholithic mantle lithosphere that was between ∼40 km and ∼75 km, and was dominated by garnet-pyroxene cumulates that developed during high magma flux phases of batholith development (Ducea, 2001; Saleeby et al., 2003). Figure 4A also recognizes the likely existence of underplated Franciscan eclogite adjacent to the arclogite root beneath the Great Valley (Saleeby, 2012). The initial thermal structure of the mantle lithosphere is a critical input to the model. Arc magmatism of the Sierra Nevada batholith was terminated in the Late Cretaceous, and its mantle wedge began cooling to a conductive geotherm in conjunction with flat slab subduction (Saleeby, 2003; Liu et al., 2010). Suprasubduction magmatism of latest Cretaceous–early Cenozoic time was accordingly displaced far inland to the foreland of the batholith. The crust of the Sierra Nevada batholith is well constrained to have been felsic and hydrous to ∼40 km depths (Ruppert et al., 1998; Fliedner et al., 2000; Saleeby et al., 2003). We combine three different lithospheric peridotite domains within the mantle lithosphere phase of our model. Mantle xenolith and late Cenozoic volcanic rock geochemical data indicate the likely existence of three such peridotite domains in the model space (Fig. 4): (1) wedge peridotites interlayered with and directly beneath the arclogite root; (2) Proterozoic subcontinental mantle to the east; and (3) underplated Farallon plate mantle nappes to the west (Ducea and Saleeby, 1996, 1998a, 1998b; Van Kooten, 1980, 1981; Beard and Glazner, 1995, 1998; Lee et al., 2001; Farmer et al., 2002; Luffi et al., 2009). Other than a possible lower thermal gradient for the mantle wedge suite and underplated Farallon plate nappes, relative to the subcontinental mantle to the east, the existing xenolith compositional and thermobarometric data indicate that in the late Cenozoic potential physical contrasts between these peridotitic domains are substantially less than contrasts between each of these domains and the arclogite root.

Middle to late Cenozoic perturbation of the subbatholithic mantle lithosphere is modeled as the result of the opening of an underlying slab window as the Mendocino Triple Junction migrated along the California margin (Atwater and Stock, 1998; Wilson et al., 2005). The approximate northern edge of the slab window through time is shown in Figure 1. Given the prolonged history of flat slab subduction inferred for the U.S. Cordillera, and the south to north opening pattern of the window, normal to the model trace, we assume a regionally flat initial thermal structure (Fig. 5). Erkan and Blackwell (2008) constructed thermal models that they claimed refute the slab window hypothesis, in favor of a stalled slab beneath the region following triple junction migration. We reject the principal interpretation of their models, considering that they ignored the preservation of the subbatholithic mantle lithosphere to ∼125 km depths based on the xenolith data (Ducea and Saleeby, 1996, 1998a, 1998b; Lee et al., 2001; Saleeby et al., 2003). Substitution of the xenolith-based subbatholithic lithospheric section with the thermobarometrically determined geotherm into the Erkan and Blackwell (2008) model setup results in an initial state indistinguishable from that of their preferred stalled slab model, which negates their discrimination between the slab window and stalled slab hypotheses.

The mantle xenolith data are used as a constraint for the initial mantle lithosphere geotherm. Utilizing the slab window hypothesis, a sharp inflection in the geotherm is imposed near the base of the mantle lithosphere, denoting the interface between previously cooled lithosphere and the underlying slab window (Fig. 5). We adopt the timing of slab window migration modeled by Atwater and Stock (1998) that places the window beneath the model space starting ca. 20 Ma (Fig. 1). We note that there are other models published for slab window opening, but the Atwater and Stock (1998) model best satisfies the timing of distinct volcanism and concurrent extensional tectonism along the southern Sierra Nevada and adjacent Great Valley (Fig. 1). We also show that the 20 Ma initiation of the basal thermal perturbation yields model results that are close in timing and kinematics to a number of distinct geological observations.

The model space is initially defined as 400 km in depth, and 600 km in cross section. The model base as well its eastern and western margins are set as free slip boundaries (Fig. 5). A 5 mm/yr extensional boundary condition is imposed along the west margin, intended to simulate an extensional component of strain along the early Neogene central California margin (Atwater and Stock, 1998). The top of the model space is a free surface with an initial topographic step imposed along the western Sierra Nevada. The topographic step is intended to simulate the position of the Sierra Nevada batholith along the western shoulder of the Nevadaplano (Fig. 4A), a regional Late Cretaceous–early Cenozoic orogenic plateau that characterized much of the U.S. Cordillera (DeCelles, 2004). The model initiates in isostatic equilibrium using a Pratt model by imposing density variations in the felsic Sierra Nevada batholith and Basin and Range crust in order to balance the Great Valley crust (Table 1). Strength envelopes used for the different model phases are summarized in Figure 7, and were discussed in detail in Le Pourhiet et al. (2006). In our current models we keep the arclogite root, lithospheric peridotite, and asthenospheric peridotite envelopes constant, and vary the most poorly constrained crustal creep parameters as represented by the upper family of envelopes. This family of crustal envelopes simulates a range of quartz and water contents, as well as the effects of structural weakening. A viscosity increase of ∼102 is imposed in the mantle at 300–400 km depth (Fig. 8A), following the common reasoning used in mantle convection models (Richards and Hager, 1984; Mitrovica and Forte, 2004; Forte et al., 2010; Cammarano et al., 2011). Other viscosity contrasts shown in Figure 8 reflect phase geometries and thermal structure.

Overview of Model Results

The initial modeling of Le Pourhiet al. (2006) focused considerably on defining suitable initial and boundary conditions (summarized in Fig. 5), and on selecting suitable rheological parameters for the various phases (summarized in Fig. 7). A series of tests were also performed to explore the potential controls of coupling versus decoupling across the Moho in that direct observations of lower crustal batholithic rocks at the southern end of the range (Fig. 4B) reveal that it is hydrated and quartz rich (Saleeby et al., 2003), rendering the lower batholithic crust weak and susceptible to flow (cf. Lowry and Perez Gussinye, 2011). These tests revealed the importance of structural strength drops in the crust and upper mantle, thereby predicting the importance of discrete geologic structures in governing mantle lithosphere mobilization and separation from the lower crust. Some of these predictions can be tested against current crustal structure (Fig. 9).

The preferred model detailed in Le Pourhiet al. (2006) runs for ∼24 m.y., and in terms of epeirogenic deformation focuses on the vertical displacement in time of a point along the eastern Sierra Nevada crest. This paper broadens its focus to include epeirogenic deformation patterns across the map area of the southern Sierra Nevada region as well as subsidence in the Tulare Basin. Subsidence data of the Tulare Basin are taken as critical geologic observations, and as discussed in Le Pourhiet et al. (2006), the crustal rheology assigned to the Great Valley has a strong influence on a number of model results, including Tulare Basin subsidence, eastern Sierra Nevada rock uplift, and delamination kinematics. Accordingly, we experimented with the rheologic parameters assigned to the Great Valley, and in doing so we have identified model parameters that influence strongly crustal deformation patterns, and have constrained those parameters to within some useful uncertainty level. We iterated further from the preferred model of Le Pourhiet et al. (2006) with focus on subsidence data of the Tulare Basin (synthesized elsewhere; Part II). Table 1 shows the rheologic parameters used in four distinct model runs that vary Great Valley crustal rheology, and Table 2 shows the simulated geologic materials for these contrasting rheologies and model results for vertical displacements of points in the center of the Tulare Basin and along the eastern Sierra Nevada crest. Model 4 best satisfies data on eastern Sierra Nevada crest rock uplift and Tulare Basin tectonic subsidence, as well as the timing of a number of distinct geologic events. This (preferred) model is similar to model A of Le Pourhiet et al. (2006), differing primarily in having a single crustal phase for the Great Valley versus in model A the addition of a second transitional phase to the felsic batholith.

Discussion of our preferred model focuses on a series of time steps in the Figure 6 phase solutions, initiating at 0 m.y. as defined in Figure 5. As the model progresses the load of the arclogite root depresses the base of the mantle lithosphere, inducing an RT instability. Convectively induced shear within the mantle lithosphere over the first ∼5 m.y. mechanically erodes a vertical arclogite pipe, which in subsequent time passively deforms within the descending instability. The effects of initial crustal structure and topography further induce asymmetry to the RT updraft pattern (5–10 m.y. transition). The eastern updraft thermally weakens and thins the mantle lithosphere beneath the upward-flexing Death Valley region. A combination of root loading and greater gravitational potential to the east of the batholith induces west-directed channel flow along the base of the felsic batholith.

Within 5 m.y. of model initiation the lower crust of the eastern Sierra Nevada–Basin and Range transition intrudes plastically westward along the felsic batholith–arclogite root interface (Fig. 6). Shear strain within the channel accelerates over the first 14 m.y. of model time (Fig. 6, insets), with the initially vertical strain markers (Fig. 5) displaying the magnitude of deformation within the channel (Fig. 6). Sensitivity analyses (Le Pourhiet et al., 2006, Fig. 6) suggest a minimum ∼10 km thickness for the lower crustal channel. Above the channel, strain localization in the elasto-plastic crust leads to the formation of a major low-angle detachment that surfaces in the Death Valley region, and that roots westward beneath the Sierra Nevada (Jones and Phinney, 1998; Park and Wernicke, 2003; Zandt et al., 2004; Fig. 9). This is shown in Figure 6 insets (10 m.y. and 14 m.y.) by the ascent and surface breaching of the high shear strain color tone in the Death Valley region. By ∼10 m.y. the mantle lithosphere beneath the Death Valley region is horizontally necked off, placing ascended asthenosphere in contact with the lower crustal channel. Shear strain accelerates both in the RT eastern upwelling and along the lower crustal channel, driving further extension above the west-dipping detachment surface. Extension above the lower crustal channel and the eastern RT upwelling produce the Death Valley–eastern Sierra Nevada extensional province at upper crustal levels (Fig 1). The shear strain insets of Figure 6 exhibit the westward jumps in the surface breaching of normal faults by high-strain-rate color tones, those of the 17 m.y. and 20 m.y. steps corresponding to the eastern Sierra Nevada escarpment system along Owens Valley (Figs. 1 and 9). At the time of the lithosphere separation event the coupled Sierra Nevada and Great Valley are calved off the Cordilleran interior as the Sierra Nevada microplate (cf. Argus and Gordon, 1991). The predicted ∼10 m.y. timing of this event is in accord with the ca. 10 ± 2 Ma timing of microplate inception based on broad geological constraints (Busby and Putirka, 2009; Saleeby et al., 2009; Part II). The westward propagation of the lower crustal channel that preceded this event, and the juxtaposition of asthenosphere with the eastern margin of the arclogite root resulting from lithosphere separation, set the stage for the rapid westward delamination of the root. We note here that there have been speculations of east-directed channel flow of weak felsic lower crust originating beneath the Sierra Nevada batholith, as a means of maintaining a relatively flat Moho between the eastern Sierra Nevada and Death Valley regions (Wernicke et al., 1996). This is not supported by our modeling, and violates strong observational constraints posed by mantle xenolith and lower crustal exposure data, both indicating that the primary batholithic crust transitioned rapidly at ∼40 km depths to drier garnet-rich mafic assemblages (Ducea and Saleeby, 1996, 1998b; Saleeby et al., 2003). The Cenozoic erosion of 7–10 km of overburden off the eastern margin of the batholith (Nadin and Saleeby, 2008) and the current crustal thickness of ∼30 km (Figs. 3B, 3C, and 9) account for the entire primary felsic crustal layer of the batholith, leaving nothing to distribute eastward.

Progressive westward delamination of the arclogite root is well expressed for the 14 m.y. to 24 m.y. time steps (Fig. 6). As the root peels away from the lower crustal channel it hinges down to the east and undergoes dip-down stretch. The model predicts that sheaths of lithospheric peridotite and lower crust flow into the area of root separation and armor the east margin of the arclogite appendage adjacent to ascending and infilling asthenosphere. The model also predicts that a substantial return flow current of previously descended lithospheric peridotite follows the ascending asthenosphere, and that these together advect relatively high temperatures to shallow upper mantle levels. In the 17–24 m.y. time steps the actively delaminating root is shown to undergo necking in conjunction with down-dip stretching. This is considered an important prediction of the model. First, this complicates attempts to retrodeform the partially delaminated root, based on the geometry of the Isabella anomaly, to its premobilized state beneath the batholith. Furthermore, fragments of the arclogite root are predicted to have completely necked off, and to have potentially reascended with lithospheric peridotites along RT return flow currents coupled to ascending asthenosphere. Such necked off and reascended lithosphere could provide the source of the apparent ∼100-km-deep high-velocity floor to the shallow asthenosphere beneath the Owens Valley and eastern Sierra Nevada region, as suggested by the eastern deep-level appendage of the Isabella anomaly (Fig. 3A), and as imaged by a number of other seismic studies (Savage et al., 2003; Boyd et al., 2004; Yang and Forsyth, 2006; Schmandt and Humphreys, 2010a; Sun and Helmberger, 2010; Gilbert et al., 2012).

Another noteworthy prediction of model 4 is the initial rapid descent of the RT over the first 5 m.y. of model time and its subsequent strong deceleration (Fig. 6); these are not artifacts of the depth of the model space. In Figure 8 the lower levels of the RT are tracked by the 1300 °C isotherm, and it is shown that the downward deceleration is related to descent through the lower level viscosity gradient. This may be an important finding of our modeling, in that USArray data for the U.S. Cordillera reveal many dangling high-velocity anomalies extending downward into the mantle from the crust (Schmandt and Humphreys, 2010b). Such a regional pattern for the U.S. Cordillera anomalies supports the hypothesized viscosity increase invoked by Richards and Hager (1984), Mitrovica and Forte (2004), Forte et al. (2010), and Cammarano et al. (2011).

Alternative Models

We found that one of the most critical, yet poorly constrained, parameters in controlling model results is crustal rheology. In the initial models of Le Pourhiet et al. (2006) and in our subsequent experiments it was found that all else being equal, for strain rates of 10−12 to 10−15 s−1, the application of the model 2–4 strength envelopes (Fig. 7) to the Basin and Range and Sierra Nevada felsic crust leads to similar results. Sensitivity analysis (Le Pourhiet et al., 2006) and model experimentation presented here show that the rheology of the Great Valley crust has a strong influence on delamination kinematics as well as the resulting vertical displacements. We have addressed the significance of the bulk crustal rheology for the Great Valley by experimenting with the four crustal strength envelopes (Fig. 7), which are intended to simulate controls by lithologic variation, degree of hydration, and structural weakening (Tables 1 and 2).

Our experiments with Great Valley crustal rheology began with a test of the common assumption (cf. Godfrey and Klemperer, 1998) that the Great Valley province is a cool rigid ophiolitic forearc. Accordingly we tested first (model 1) a dry diabase rheology for the Great Valley crust (after Tsenn and Carter, 1987). The principal result was a lack of arclogite root delamination, and epeirogenic signals of wrong sign including ∼250 m net subsidence of the eastern Sierra Nevada crest and ∼150 m net rock uplift in the Tulare Basin (Table 2). Nevertheless, the deeper mantle lithosphere is mobilized as an RT instability (Fig. 10A). These results are considered significant for two reasons: (1) in conjunction with the other model results discussed in the following, they show that the correctly polarized and significantly greater observed vertical displacements are arclogite root delamination (versus RT-driven lower lithosphere displacements); (2) they are consistent with basement core petrologic and crustal structure data showing that a weak rheology dominated by quartzose and/or granitic and/or hydrated and structurally weakened material is most appropriate for the Great Valley crust (May and Hewitt, 1948; Saleeby, 2007, 2012).

Additional experimentation entailed a progressively weaker rheology for the Great Valley crust relative to model 1 (models 2–4; Tables 1 and 2). In each case root delamination occurs within the geologically constrained time interval, much like that shown for model 4, with delamination linked to a lower crustal channel that extends beneath virtually the entire felsic batholith (Fig. 10). However, in the more viscous models (2 and 3) the load of the root is more efficiently coupled to the crust, leading to inordinately high subsidence results in the Tulare Basin (Fig. 11; Table 2). Model 2 and 3 predictions for eastern Sierra Nevada crest rock uplift are similar to the model 4 predictions, but notably less (Fig. 11; Table 2). The model 4 predictions for amount and timing of eastern crest uplift are more in line with observation (Part II). Suppression of root delamination by a highly rigid Great Valley crust (model 1) is explored further in Figure 10. In Figure 10A the 20 m.y. step of models 1–4 are compared in phase solution. Model 1 is the only model of the four for which the lower crustal channel is limited to the eastern half of the batholith. The shear strain rates at the 5 m.y. time steps for models 1–4 are compared in Figure 10B, which shows that a highly rigid Great Valley crust (model 1) seals off the lower crustal channel along the western margin of the batholith; this further suppresses flow for a significant distance eastward beneath the batholith. The model 1 flow pattern fails to promote root delamination, which is in contrast to models 2–4, which are characterized by west-directed channel flow across the entire cross-sectional extent of the batholith, in continuity with channel flow beneath the western Basin and Range.

In the following discussions on the regional patterns of late Neogene–Quaternary rock uplift and tectonic subsidence, we show that model 4 produces results that are closest to observations. This model incorporates a crustal rheology for the Great Valley that is similar to that of the Sierra Nevada batholith, justified by Great Valley basement core and crustal structure data (May and Hewitt, 1948; Saleeby, 2007, 2012; Figs. 4A and 9). We consider the correspondence of the model 4 vertical displacement predictions with observations, and the correspondence of our chosen crustal rheology for the Great Valley with observational constraints for that of the Great Valley as reasonable validations of model 4, and proceed with more detailed focus on the predictions of this model. Inasmuch as slab window migration beneath the model region progressed over a period of ∼10 m.y. we ran our models forward for 24 m.y. as a means of assuring ample time for instability nucleation and growth. The first 20–23 m.y. of model time produce results that most closely resemble rock uplift and tectonic subsidence data and the current geometry of the partially delaminated arclogite root. The slab window migration pattern of Atwater and Stock (1998) predicted the opening of the window beneath the model space ca. 20 Ma, which is in accord with the results of our modeling. This migration pattern is also in agreement with the ca. 22 Ma initiation of a distinct belt of slab window–related volcanic rocks (Fig. 1), and associated extensional structures that extend across the southernmost Sierra Nevada region (Evernden et al., 1964; Bartow and McDougall, 1984; Coles et al., 1997; Sharma et al., 1991; Monastero et al., 1997; Maheo et al., 2009; Part II). Therefore, we correlate here ca. 20 Ma in geologic time with the initiation of our model time (0 m.y.), and interpret 22 +1/–2 m.y. model time with the geologic present.

Model Predictions for Volcanism, Heat Flow, and Seismicity

Late Cenozoic volcanism is widespread in the southern Sierra Nevada region (Van Kooten, 1980, 1981; Moore and Dodge, 1980; Manley et al., 2000; Farmer et al., 2002). A first order test for our preferred model is a reasonable prediction of the timing and source characteristics of such volcanism. Late Cenozoic volcanism of the region is fundamentally basaltic, with scattered late-stage silicic centers. Our preferred model can be evaluated with respect to when and where the upper mantle is predicted to have ascended vigorously enough to have undergone decompression partial melting, which mantle phases were likely involved in a given melting regime, and whether, and when, the production of significant silicic magma is predicted. There are three distinct pulses of late Cenozoic volcanism in the region (Moore and Dodge, 1980; Manley et al., 2000; Farmer et al., 2002; Phillips et al., 2011): principally late Miocene (8–12 Ma), late Pliocene (ca. 3.5 Ma), and late Quaternary (after 1 Ma). Each age group shows evidence for hybridization of mantle lithosphere and depleted mantle source components, with local increase of depleted mantle components eastward with time (Van Kooten, 1981; Ducea and Saleeby, 1998a; Farmer et al., 2002; DePaolo and Daley, 2000; Blondes et al., 2008; Putirka et al., 2012). The depth of melt extraction decreases with time from ∼100–125 km to ∼40–75 km, and finally into the lower felsic crust for silicic members (Van Kooten, 1980; Ducea and Saleeby, 1996, 1998a; Putirka et al., 2012). Figure 1 shows the areas over which the ca. 3.5 Ma basaltic and Quaternary bimodal volcanic rocks occur, grouped together as 4 Ma and younger. The 8–12 Ma volcanic rocks are dispersed over a slightly greater region, but appear to have been less voluminous than the products of the younger eruptions.

The Figure 6 model predictions are in line with the timing and source characteristics of late Cenozoic volcanism of the region. Figure 12 shows a plot of geotherms predicted for selected time steps of model 4 beneath the eastern Sierra Nevada crest, along with the peridotite melting field and dry granite solidus. Partial melting of peridotite with bulk H2O as low as ∼0.05% at depths of ∼100 km is an initial state of our model, dictated by the initial geotherm. Geologically this corresponds to the emplacement of asthenosphere into the slab window. Over the first 10 m.y. of model time the geotherm relaxes slowly, and then it is perturbed by RT-driven convective overturn, resulting in a temperature inversion that focuses mantle melt production at a depth range of ∼100–125 km, in line with observations. By 17 m.y. model time the inversion has grown and ascended, and melt production in the upper mantle is focused to the observed depth range of ∼75–40 km. The 17 m.y. geotherm also crosses the dry granite solidus at basal crustal levels, while the 20 m.y. geotherm becomes steeper over lower crustal depths and approaches the dry granite solidus at basal crustal depths. This is in line with the appearance of silicic volcanism late in the youngest cycle of basaltic volcanism. Because the modeled geotherms of Figure 12 are for a single point along the model trace, it is also useful to study potential melt production sites, along with potential source components, by direct inspection of the Figure 6 thermal and phase structures. The 10 and 14 m.y. time steps show rapid ascent of asthenosphere and related attenuation of mantle lithosphere leading to lithosphere separation. The 1300 °C and 900 °C isotherms become highly compressed across the mantle lithosphere-asthenosphere transition in the axial to eastern Sierra Nevada; this promotes the interaction and hybridization of partial melts derived from the two mantle phases at ∼100 km depths. To the east of the batholith, depleted mantle becomes a greater source component, as observed, while the eastern RT updraft ascends beneath the region and necks off the mantle lithosphere. Detailed comparison of the 14 and 17 m.y. time steps offers a further rationale for the interpretation of the ca. 3.5 Ma volcanic pulse. As the arclogite root rapidly peels back and hinges down to the east, sheaths of mantle lithosphere and entrained lower crust armor its upper detached surface as asthenosphere intrudes in along the zone of progressive delamination. Both a temperature inversion and steep horizontal gradient bound the upper and eastern borders of the partially delaminated root and its labile sheaths. This offers an ideal setting for the generation of variably hybridized to highly enriched partial melt products at depths as shallow as ∼40 km.

Study of the 17–24 m.y. time steps reveals that fronts of the 900 °C isotherm migrate well into the lower crust beneath the eastern Sierra Nevada–western Basin and Range region in what we broadly consider modern time. Heating of the lower crust to these temperatures should result in substantial silicic melt production, as a function of H2O content. This provides a mechanism for Quaternary silicic volcanism of Long Valley, the Coso Range, and scattered silicic plugs and domes of the Owens Valley and adjacent southeastern Sierra Nevada region (cf. Moore and Dodge, 1980). The 20 m.y. step also shows lingering high temperatures in interlayered asthenospheric and lithospheric peridotites in the uppermost mantle beneath the eastern Sierra Nevada region (closed 1300 °C isotherm). The persistence of high upper mantle temperatures in such materials immediately beneath the base of the crust offers an ideal circumstance for the generation of modest amounts of variably enriched basaltic melt, given a suitable perturbation such as rapid normal fault displacements, or crustal unloading from rapid glacial retreat or catastrophic draining of large glacially fed lakes (e.g., Glazner et al., 1999).

The results of our modeling address the apparent conflict between seismic data that indicate the complete removal of mantle lithosphere from beneath the eastern Sierra Nevada region (Figs. 3A, 3B), and geochemical data indicating that basaltic rocks erupted during and following lithosphere removal contain substantial continental lithosphere components (cf. Putirka et al., 2012). Our preferred model predicts complex regional interlayering of ascended asthenosphere, mobilized mantle lithosphere, delaminated root, and entrained lower felsic crust, as well as compressed isotherms that migrate across the layering as it develops. Such a dynamic structure will render basaltic melts enriched in lithosphere components from asthenosphere that is defined solely by mechanical criteria. In contrast to model 4, models 1–3 do not produce the extensive interlayering of mantle lithosphere and asthenosphere of model 4 (Fig. 10A); this is explained in sensitivity tests (Le Pourhiet et al., 2006) that explore the influence of coupling versus decoupling across the Moho. Lower degrees of coupling promote more vigorous convective flow in the upper mantle, as well as zones of structural weakening in the upper mantle that become shear zones. Both of these effects increase tectonic disruption of the primary upper mantle structure, and promote more extensive interlayering of the major phases. The progressively weaker whole-crust rheologies encompassed in models 1–4 promote progressively lower degrees of coupling across the Moho, making model 4 the best solution for the generation of hybridized basaltic melts during the course of mantle lithosphere removal.

Figure 1 also shows a region across the southern Sierra Nevada that is characterized by a modern anomalous thermal transient. This region is uniquely characterized by the coincidence of low heat flow measurements in batholithic exposures, typical of much of the Sierra Nevada microplate interior, with the widespread occurrence of warm and hot springs and warm-water wells within the exposed batholith, and with anomalously high downhole temperature measurements in Tertiary strata of the Kern arch (Part II). The area of the anomaly also, at least in part, overlies negative conversions in seismic receiver function data at shallow mantle and mid-crustal depths that could correspond to zones of partial melt (Yan and Clayton, 2007; Frassetto et al., 2011). We note here the continuous peripheral distribution of the 4 Ma and younger volcanic rocks and modern thermal transient zone relative to the delamination hinge trace.

A high state of modern deviatoric stress is predicted in our preferred model for Moho levels proximal to the delamination hinge, and for a large domain of the crust down to ∼30 km depth across the Great Valley–western Sierra foothills transition (Fig. 6; 20 and 24 m.y. steps). Focal depths beneath the western foothills and eastern Great Valley, between 36.5°N and 37.5°N, cluster between 12 and 38 km, defining some of the deepest events in the U.S. Cordillera (Wong and Savage, 1983; Gilbert et al., 2007; Frassetto et al., 2011). Focal mechanisms are varied, including numerous oblique normal and oblique reverse events, but overall define approximately north-south σ1 and approximately east-west σ3. Calculation of bending stresses arising from flexure of the crust above the delamination hinge and residual root-loaded Tulare Basin indicate a high state of approximately east-west tensile stress below ∼10 km, and modest approximately east-west compressive stress in the upper crust above the neutral fiber (Le Pourhiet and Saleeby, 2013). Microseismicity strain inversions for the region (Unruh and Hauksson, 2009; J.R. Unruh, 2012, written commun.) also resolve approximately east-west extensional coseismic strain across much of the region.

In summary, model results for rapid vertical displacement components in the upper mantle related to mantle lithosphere removal, along with the modeled deformational geometry of the principal mantle phases and resulting thermal structure, reasonably predict the timing and principal source components for late Cenozoic volcanism of the region. Our model predictions also yield a relatively high state of deviatoric stress for distinct areas of the crust beneath the western Sierra foothills and adjacent Great Valley. These areas correspond to distinct patches of modern seismicity. We consider these to be reasonable validation tests for our preferred model; further validation tests by observations of surface vertical displacements through time and space are provided in the following.

COMPARISON OF MODEL RESULTS WITH OBSERVED EPEIROGENIC DISPLACEMENTS

Epeirogenic Displacements

Our modeling results show that vertical displacements throughout the course of mantle lithosphere mobilization and arclogite root delamination are driven by both flexure of the elasto-plastic crust and by buoyancy changes in the upper mantle and lower crustal channel. Flexure is driven primarily by the load of the root, with a small component related to footwall uplift along the eastern Sierra Nevada escarpment system. Vertical forces related to the changing distribution of the root load during the course of delamination and buoyancy changes as the mantle lithosphere separates and the lower crustal channel evolves drive vertical displacements across the Sierra Nevada uplands that are not necessarily unidirectional in time at any given point (Fig. 13). We focus first on the model 4 predictions for uplift of the eastern Sierra Nevada crest through time (Fig. 11). The modeled rock uplift is a minimum because it does not account for exhumation during rock uplift, nor does it account for sediment loading in the Tulare Basin, which induces an increment of eastern crest rock uplift via flexure. In that the principal phase of rock uplift has occurred over the past ∼6 m.y., and that both long- and short-term erosion rates of Sierra Nevada interfluve areas are ∼0.01–0.05 mm/yr (Small et al., 1997; Riebe et al., 2000; Stock et al., 2004; Clark et al., 2005; Phillips et al., 2011), the effects of interfluve area exhumation are considered to be trivial. In contrast, the flexural-isostatic responses to major drainage basin widening and deepening across the Sierra Nevada, and sediment loading in the Tulare Basin, are likely to be important (Small and Anderson, 1995; Pelletier, 2007; Part II). These will be discussed elsewhere (Part II) along with regional U.S. Cordillera epeirogenic patterns that were likely additive to more focused southern Sierra Nevada rock uplift, on which we focus here.

The rock uplift curve for model 4 (Fig. 11) shows ∼300 m of subsidence at the modern Sierra Nevada crest during the first 10 m.y. of model time, as a result of extension and loading of the arclogite root. Subsidence immediately east of the Sierra Nevada is more intense during this time interval, reflecting the underlying progression of mantle lithosphere separation (Fig. 6). As the Sierra Nevada microplate separates at ∼10 m.y. the area of the modern Sierra crest undergoes little vertical displacement, but by ∼15 m.y. the modern crest region begins to slowly rise as the eastern Sierra uplift initiates along the proximal footwall of the eastern Sierra escarpment at ∼17 m.y. (Fig. 13). From ∼17–20 m.y. maximum uplift rolls westward to the modern crest area, and culminates at ∼22 m.y., after which the crest region subsides. The eastern Sierra Nevada crest rock uplift curve during the 15–22 m.y. interval of model time has a small inflection reflecting changing flexural and buoyancy components. A total of ∼800 m of crest rock uplift occurs over this time interval, progressing at an average rate of 1.14 mm/yr. Elsewhere (Part II) we develop a rationale for deriving ∼400 m of additional crest uplift as a result of drainage basin exhumation and Tulare Basin sediment loading, stimulated by the ∼800 m of uplift driven by delamination. We consider the resulting ∼1200 m of ∼14–22 m.y. modeled rock uplift to be correlative to the mainly Pliocene–Quaternary uplift pulse that is hypothesized for along the Sierra Nevada crest region (Huber, 1981; Wakabayashi and Sawyer, 2001; Stock et al., 2004).

Direct quantitative constraints of rock uplift for the eastern Sierra Nevada crest are lacking. Incision data for Kings River, based on cosmogenic dating of cave sediments, at a location roughly half-way between the eastern crest and the edge of the Great Valley along the Figure 3A section trace (Fig. 1), indicate ∼400 m of rapid incision between 2.7 and 1.5 Ma (Stock et al., 2004). Stock et al. (2004) used a stream power–based numerical model initiating with 1500 m of eastern crest uplift at 9 Ma to explain their data. The model produces the amount and timing of incision observed by headward erosion and knickpoint migration from the western range front to the interior of the range. A dramatic post–1.5 Ma drop in incision rates resolved along the Kings River is interpreted as having resulted from the protective mantling of the main river channel by coarse glacial debris (Stock et al., 2004). The model 4 results offer a potential alternative, or additional mechanism for this change in incision rates, with the predicted postdelamination eastern crest subsidence phase that initiates at ∼22 m.y. (Figs. 11 and 13).

Anomalous subsidence in the southern Great Valley is predicted by our modeling to have developed by a combination of regional flexure through the course of root delamination, and loading by the western margin of the arclogite root over its zone of residual attachment. The anomalous domain currently corresponds to the Tulare Basin (Fig. 2), although we assert that in Pliocene–early Quaternary time the anomalous subsidence zone also extended across the area of the current Kern arch (Part II). The tectonic subsidence curve predicted by model 4 for the Tulare Basin (Fig. 11) shows as much as ∼200 m rock uplift over the first 5 m.y. of model time, followed by 5 m.y. of little change. Some portion of this early uplift of the basin could be an artifact from the smoothing out of the initial topographic step along the western Sierra Nevada (Fig. 13), although Loomis and Glazner (1986) pointed out that the southern San Joaquin Basin (Fig. 2) appears to have undergone 200–400 m regional uplift at ca. 15 Ma. Loomis and Glazner (1986) attributed this to the passage of the Mendocino Triple Junction. Our modeling offers the potential additional mechanism of broad flexural uplift of the basin sympathetic to the initial loading of the lithosphere by the arclogite root, prior to delamination. Model 4 subsidence in the Tulare Basin begins at ∼10 m.y. during lithosphere separation, and proceeds at an average rate of 0.09 mm/yr until ∼22 m.y. (Fig. 11).

The resolution of tectonic subsidence components in the Tulare Basin that may be related to mantle lithosphere removal is complicated because in the area of the subbasin, and farther south along the San Joaquin Basin, a Neogene age southward slope culminating in a deep marine basin was superposed over the regional west-southwest slope of the Cretaceous Great Valley forearc basin (Part II). This is further confounded by the late Cenozoic regional west tilt of the microplate and the building out of a westward-thickening sedimentary prism across the Great Valley. We address this problem by resolving tectonic subsidence residuals unique to the Tulare Basin, relative to Great Valley–wide subsidence patterns for late Neogene–Quaternary time (Part II). Our best estimation of this residual is 625 m since 7 Ma; our best estimation of equivalent model time is ∼14–22 m.y., for which our modeling predicts ∼680 m of tectonic subsidence (Fig. 11 and Part II–Table 2).

An additional useful comparison between observed and modeled tectonic subsidence is by graphic overlay of the model 4 subsidence curve onto published tectonic subsidence curves for the Tulare Basin (Moxom and Graham, 1987; Part II). In Figure 14 we compare subsidence curves from wells drilled at various locations along the Tulare Basin, and a well from the Kern arch. Wells used for Figures 14A and 14B are from the northern and eastern margins of the basin, and the well used for Figure 14C is from the central part of the basin. The record of the well from the Kern arch (Fig. 14D) indicates a similar subsidence history with the 14C well, until abruptly terminated by Quaternary rock uplift. The modeled curve is shown in magenta on each plot; the curves in Figures 14C and 14D show similar amounts of tectonic subsidence over similar time intervals for the model and observed curves. Curves in Figures 14A and 14B show observed subsidence to be less than that modeled, which is consistent with the two wells being located on the margin of the anomalous subsidence domain (Part II–Figs. 5 and 8). Each observed curve shows a modest inflection between 10 and 20 Ma, which could be interpreted as basin uplift related either to passage of the Mendocino Triple Junction (Loomis and Glazner, 1986) or to upward flexure related to the initial loading of the arclogite root. Robust paleobathymetric data on which Loomis and Glazner (1986) based their analysis are restricted to the deeper southwest portions of the San Joaquin Basin, and thus the uplift signal that they recognize is best recorded for the region of the Figure 14D well (Fig. 2). The model predictions for this early uplift phase of the basin may be partly contaminated with a modeling artifact.

Neogene–Quaternary Topographic Evolution

The model results for southern Sierra Nevada rock uplift pertain to uplift resulting primarily from arclogite root delamination. The model does not account for other possible far field induced epeirogenic effects that may have operated through the Cenozoic on the U.S. Cordillera as a whole (cf. Suppe et al., 1975; Unruh, 1991; Humphreys, 1995, 2008; Lowry et al., 2000; Pierce et al., 2002; Mix et al., 2011). The implications of such far field driven epeirogeny for Sierra Nevada microplate topographic evolution are pursued elsewhere (Part II), where a rationale is developed to isolate distant and more local delamination uplift components. Here we focus on topographic evolution arising primarily from root delamination, which is considered distinct from the regional elevation increase of the U.S. Cordillera.

Topographic evolution profiles derived from model 4 are shown in Figure 13. The area of the Sierra Nevada crest and regions to the east first undergo a period of subsidence driven by distributed crustal extension and flexure off the margin of the internally loaded arclogite root. The spikes along the eastern segments of the profiles represent the surface breaching of normal faults that root to the west into the lower crustal channel (Fig. 6 strain insets; Fig. 9). The region is primarily under subsidence for the first ∼15 m.y., which in geologic time culminated with the production of widespread Late Miocene–Pliocene lake basins of the eastern Sierra Nevada–Owens Valley region (Bachman, 1978; Bacon et al., 1982). By 16 m.y. subsidence proximal to the eastern Sierra Nevada was replaced by vigorous rock uplift, which we equate to the widely hypothesized phase of Pliocene eastern Sierra crest rock uplift (Huber, 1981; Unruh, 1991; Wakabayashi and Sawyer, 2001; Stock et al., 2004).

Eastern Sierra Nevada crest rock uplift accelerates and migrates westward between 17 and 22 m.y., while broad subsidence and distributed normal faulting continue to the east (Figs. 11 and 13). Eastern crest elevation increase, as controlled by root delamination, peaks at 22 m.y., and then the eastern crest region begins its postdelamination subsidence phase. The topographic culmination of the crest continues to migrate westward through the peak elevation phase. The 14–22 m.y. modeled growth of Sierra Nevada crest topography is equated to the growth and migration of the delamination bulge, and is not considered to be present throughout the entire Sierra Nevada (Fig. 2; Part II).

The topographic evolution curves of Figure 13 translate into tectonic subsidence and superposed rock uplift evolution in the Tulare Basin. The eastward and subsequent westward migration pattern of the maximum subsidence is consistent with the timing and geometry of the embayment of the eastern Tulare Basin into the western Sierra foothills, and the subsequent ongoing incision pattern of the aggraded fluvial sediments that filled the accommodation space created by the eastward migration pattern (Saleeby and Foster, 2004; Part II). Distinct provenance and facies relations for Pliocene–early Quaternary strata along the eastern San Joaquin Basin that appear to record a similar transgressive-regressive cycle are discussed elsewhere (Part II).

The Figure 13 curves suggest that as much as ∼500 m of rock uplift occurred in the Coast Ranges as a result of flexure related to root loading and delamination. We refer to this component of Coast Range uplift as the sympathetic (flexural) bulge to the delamination bulge. Additional topographic growth of the Coast Ranges is driven by east-directed blind thrusts that root into the San Andreas fault (cf. Wentworth et al., 1984; Fig. 9). Such east-directed thrusting is consistent with the force balance of our model, which suggests transverse shortening of the Coast Ranges driven by excess gravitational potential of the Sierra Nevada uplift (cf. Jones et al., 1994, 2004). Note that the Quaternary rise of the central foothills swell (Fig. 2) is consistent with the Quaternary south to north pattern in root delamination beneath the Kern arch region (Fig. 3C) (discussed further elsewhere; Part II). In this configuration the central foothills swell represents the sympathetic bulge to the Kern arch delamination bulge (cf. Figs. 2, 3C, and 13). In Figure 2, we merge this sympathetic bulge with the outer shoulder of the principal Pliocene–Quaternary delamination bulge rendering the composite outer shoulder. This is under further investigation using our three-dimensional modeling procedures.

Remnants of the arclogite root are still intact beneath the western foothills and adjacent Great Valley as far north as ∼38°N (Schmandt and Humphreys, 2010a; Reeg, 2008; C.H. Jones, 2011, written commun.; Figs. 1 and 3C). Additional tomographic inversions based on Rayleigh wave shear velocities (Gilbert et al., 2012) indicate that the peridotitic lithosphere continues northward from the area of residual arclogite root attachment along the eastern Great Valley and western foothills. The northward pinching out of the arclogite-rich mantle lithosphere layer is shown in Figure 4B as a primary feature of the Sierra Nevada Cretaceous mantle wedge. Such pinching out of the arclogite, in conjunction with its tectonic truncation along the Coast Range–Rand megathrust lateral ramp (Fig. 4B and Part II), produced the early Cenozoic ∼300-km-long, ∼100-km-wide, ∼35-km-deep tabular shape of the arclogite root (Fig. 4). As shown in Figure 4, this tabular body had gradational (primary) northern and western borders with adjacent peridotitic lithosphere, and sharper tectonic boundaries along its eastern and southern borders. From this we posit that surface manifestations of root removal would be more pronounced along its eastern and southern margins, and more subdued along its western and northern margins. Pronounced features along its eastern and southern margins, such as the Death Valley–eastern Sierra Nevada extensional province, Kern Canyon fault system, and Kern arch are quite evident.

Surface geological features that are readily attributed to root delamination, such as distinct basaltic volcanism across the southern Sierra Nevada and anomalous subsidence in the Great Valley, diminish northward from ∼37°N–37.5°N, and we accordingly fade out the surface trace of the delamination hinge in this region (Figs. 1 and 2). We likewise fade out the principal delamination bulge northward across this area, which is consistent with regional relief and geodetic patterns (Part II; see Hammond et al., 2012). Approaching the region of the Stanislaus River drainage (Fig. 2) from the south, and continuing northward, the geomorphology of the Sierra Nevada changes to a much more regular profile with a clear regional west tilt pattern that is common to Tertiary strata dipping off the west flank of the range, and to well-defined planar interfluve surfaces that can be traced up to the eastern Sierra crest (Wakabayashi and Sawyer, 2001; Saleeby et al., 2009; Part II). Southward from the Stanislaus drainage the topography becomes progressively more complex until the drainages of the Kings and Kaweah Rivers (Fig. 2); the west-tilt pattern is not readily visible. The principal geomorphic difference in the Sierra Nevada uplands between the region of the San Joaquin to Kaweah drainages, and the region of the Stanislaus to Mokelumne drainages, is the superposing of the delamination bulge across a northern Sierra–like regional topography (Part II). This is analogous to the difference between the Pliocene–Quaternary subsidence in the Tulare Basin and the rest of the Great Valley to the north being the superposing of delamination-related subsidence across regional Great Valley subsidence patterns.

DELAMINATION IN THREE DIMENSIONS

General

In Figure 2 three epeirogenic domains are differentiated for the southern Sierra Nevada region that may be considered anomalous relative to the rest of the range. Anomalous is defined in this context as rock uplift components that are not clearly expressed in all other regions of the Sierra Nevada (Part II). These three rock uplift domains are analogous to the anomalous tectonic subsidence domain of the Tulare Basin, but are more aerially expansive due to broad flexure off the margin of the residual arclogite root. The two principal anomalous rock uplift domains are differentiated as areas over which landscape features primarily reflect erosional responses to Late Miocene(?)–Pliocene rock uplift (Pliocene maxima) versus late Pliocene(?)–Quaternary rock uplift (Quaternary maxima). We also define a third outer shoulder domain, which we consider a composite feature (Part II). Comparison of Figures 1 and 2 shows that the area where anomalous rock uplift was dominant in the Pliocene coincides with the principal area of 4 Ma or younger delamination volcanism, and the area dominated by Quaternary rock uplift coincides with the area of the modern anomalous thermal transient. This thermal transient is interpreted as the initial thermal pulse resulting from ongoing delamination at depth. The thermal signal is carried primarily in hot aqueous fluids that are circulating within southern Sierra Nevada basement rocks and the deeper portions of Tertiary strata of the Kern arch (Part II).

Considering the coincidence of the principal area of 4 Ma or younger volcanism with the area of apparent waning rock uplift, and the coincidence of the anomalous thermal transient with the area of most recent rock uplift together, in relation to the surface trace of the delamination hinge, leads to the following hypothesize for regional delamination kinematics. In Miocene–Pliocene time east to west–directed delamination initiated along the entire eastern edge of the arclogite root and progressed into early Quaternary time. By mid-Quaternary time delamination along the southern end of the root began to transition into a south to north separation pattern (Fig. 3C), as east to west separation along the principal hinge slowed. In this context the anomalous thermal transient represents a thermal head wave that is perhaps the precursor to near future volcanism in the southern Sierra Nevada–Kern arch region. The brittle response of the upper crust follows a similar temporal-spatial pattern (Fig. 2). Normal fault scarps of the eastern Sierra Nevada fault system become more subdued and are mostly obscure southward from ∼36°N, which is the approximate latitude that active west-side-up normal faulting along the Kern Canyon system starts and continues southward to the end of the range (Amos et al., 2010; Nadin and Saleeby, 2010). This step-over pattern in active normal faulting also controls the west tilt pattern of the range and the attitude of Tertiary strata that are along the western foothills and extend into the Great Valley subsurface (Saleeby et al., 2009; Part II). The southwest terminus of the delamination hinge trace coincides with the active Pond-Poso fault system, which circuits southeastward into normal faults along the Kern front and Kern range front system (Fig. 2). This complex zone of normal and transfer motion terminates southeastward near the southern terminus of active scarps of the Kern Canyon system, and also bounds the southwest margin of the area affected by the anomalous thermal transient. Active normal faulting of the Kern Canyon and Pond-Poso-Kern front system is thus partitioning active epeirogenic uplift related to ongoing delamination along the curved southwest segment of the delamination hinge.

Considering the preceding discussion in conjunction with the initial geometry of the arclogite root and the subsidence history of the Kern arch lends further insight into regional delamination kinematics. The initial root had a tabular geometry, its southern margin tectonically truncated at ∼35.3°N (Fig. 4). Rock uplift and sediment dispersal data for the southeasternmost Sierra Nevada record latest Miocene–Pliocene relief generation (Part II) that is analogous to that documented and modeled for the Kings River drainage (Stock et al., 2004). Furthermore, subsidence data for the Kern arch (Fig. 14D) indicate that the area of the arch shared subsidence history similar to that of the Tulare Basin, prior to Quaternary uplift of the arch. Thus the subsidence and uplift data support a regional pattern of east to west delamination of the root along virtually the entire southern Sierra Nevada. When this transitioned in Quaternary time to the south to north delamination of the southern end of the residual root, the Tulare Basin was partitioned off from the rest of the San Joaquin Basin by the rise of the Kern arch.

Implications of Thermomechanical Modeling for Three-Dimensional Delamination

Neogene volcanism across the southernmost Sierra Nevada and adjacent Mojave plateau region initiated ca. 22 Ma and progressed until ca. 16 Ma (Evernden et al., 1964; Bartow and McDougall, 1984; Coles et al., 1997; Sharma et al., 1991; Monastero et al., 1997; our research). Geochemical data suggest that the lavas were produced by asthenosphere decompression partial melting and lower crustal assimilation as the slab window opened beneath the region (Sharma et al., 1991). This is in accord with the 23–24 Ma arrival of the slab window edge beneath the region (Atwater and Stock, 1998). Thereafter the window opened progressively northward, reaching the Figure 6 model trace at ca. 20 Ma (Fig. 1). Inasmuch as the arclogite root existed along most of the southern Sierra Nevada batholith (Fig. 4B), why didn’t delamination initiate in the south, and progress northward along the long axis of the root as the slab window opened? Our modeling results together with initial crustal structure offer a viable answer. The crust along the eastern edge of the root from ∼35.5°N northward was under significantly greater gravitational potential of the Nevadaplano than either the western margin of the root beneath the Great Valley (Fig. 4A), or the southern margin of the root adjacent to the highly extended terrane of the southernmost Sierra Nevada–Mojave region (Fig. 4B). In all of our model runs (and in those of Le Pourhiet et al., 2006), the higher topography of the Sierra Nevada batholith and its foreland, relative to the Great Valley forearc, served to drive westward flow in the lower crustal channel, and to effectively weaken the crust such that extension was concentrated to the east of the batholith (cf. McKenzie et al., 2000). The low regional elevations of the highly extended terrane to the south thus inhibited the development of a north-directed lower crustal channel along the felsic batholith-arclogite root interface, thereby inhibiting an initial south to north delamination pattern.

We pose the following scenario for the progressive delamination of the arclogite root to its current state. Based on the predictions of the Figure 6 model, the entire eastern margin of the root began mobilization ca. 10 Ma (∼10 m.y. model time), in response to lithosphere separation and initiation of the Sierra Nevada microplate. Starting at this stage extension to the east of the Sierra Nevada batholith operated in conjunction with left-slip transfer motion along the Garlock fault, with the development of the Death Valley–eastern Sierra Nevada extensional province, as outlined in Figure 1 (after Monastero et al., 1997; Snow and Wernicke, 2000). According to our model, upper crustal extension in this province was linked to a regional underlying lower crustal channel along which lower crust of the western Basin and Range province was drawn into the migrating zone of the delamination hinge. The magnitude of extension predicted by our model does not equate to the large magnitudes suggested by Snow and Wernicke (2000), but our model does not discount the possibility of additional components of extension arising from far field driven plate motions. Furthermore, the magnitude of extension in the eastern Sierra Nevada to Death Valley region is debated (cf. Applegate and Hodges, 1995; Renik et al., 2008).

Between 10 and 6 Ma the entire along-strike extent of the arclogite root began to peel away from the base of the crust, its separation progressing in a westward direction. By ca. 5 Ma a regional-scale limb of the root had rotated down to the east while undergoing downdip stretch (Fig. 15A). As this regional fold developed its hinge area increased the flexural rigidity of the root in a longitudinal direction, further suppressing south to north components of delamination. In Figure 15B we hypothesize that at ca. 4 Ma a large fragment of the partially delaminated root in the region between ∼36°N and ∼38°N necked off and foundered. We envisage the necking off of a megaboudin that was two to three times larger in downdip extent than that shown necking off between the 17 and 20 m.y. time steps of Figure 6, perhaps similar to the megaboudin produced in the 20 m.y. time step of model 3 (Fig. 10A). A combination of upside-down melting of the arclogite megaboudin (cf. Elkins-Tanton, 2007), inclusive of entrained felsic lower crust and peridotitic lithosphere layers, and the hybridization of these melts with decompression partial melts of the complementary asthenospheric upwelling produced the enriched ca. 3.5 Ma volcanic pulse. This necking and foundering event briefly enhanced the shallow upper mantle melting conditions predicted for Pliocene time (17 m.y. model time) in Figure 12.

Once the arclogite megaboudin detached, the residual limb of partially delaminated root to the south began peeling away from the base of the felsic crust with progressively greater south to north components (Fig. 15C). According to the analysis given in conjunction with Figure 10, the lower crustal channel had already been established along the southern reaches of the root during the east to west phase of root delamination. By ca. 1 Ma the south to north delamination pattern was in motion (Fig. 15E), and the Kern arch began its assent, disrupting the southern zone of anomalous subsidence and leaving the residual Tulare Basin over the principal residual root. The south to north pattern of delamination that is suggested in the longitudinal structure section of Figure 3C should produce a flexural bulge sympathetic to the delamination bulge, the most active part of which propagated to the Kern arch area by ca. 1 Ma. Geologic features suggesting that the central foothills swell (Fig. 2) is the expression of this sympathetic bulge are discussed elsewhere (Part II).

Both the thermal and epeirogenic evolutions of the southern Sierra Nevada region suggest the progression of events outlined in Figure 15. The 1 Ma or younger emplacement of hot asthenosphere beneath the crust above the northward-delaminating root segment is actively driving the thermal transient beneath the southern Sierra Nevada–Kern arch region as well as contributing to rapid rock uplift of the region (Figs. 1 and 2). In contrast, along the northern segment of the delamination zone the ca. 3.5 Ma volcanic pulse records the upper mantle climax in that region’s delamination-related thermal transient; subsequent Quaternary, mainly silicic, volcanism reflects lower crustal responses to the lingering high thermal gradient. The apparent waning of rock uplift in the northern region (Fig. 2) signals the initial stages of decay of the ca. 3.5 Ma thermal transient, and the onset of the postdelamination dynamic regime in that region.

In Figures 15C and 15D we further posit that the arclogite root megaboudin that detached at 4–3 Ma was partially entrained in the eastern RT upwelling, suspending it at ∼200 km beneath the eastern Sierra Nevada–Owens Valley region. Such entrainment would in theory also carry with it a selvage of mobilized peridotitic lithosphere, as displayed in the 20 and 24 m.y. time steps (Fig. 6). The entrainment and hosting return flow pattern could account for the eastern appendage of the Isabella anomaly that is beneath ascended asthenosphere of the eastern Sierra Nevada–Owens Valley region (Figs. 3A, 3B). We further note that the scale of the modeled thermal anomaly resulting from the RT mobilization of the lower lithosphere, as approximated by the depressed 1300 °C isotherm (Figs. 6 and 8), is similar in scale to the (greater) Isabella anomaly (dVp ≥ +1%) in transverse and diagonal sections (Figs. 3A, 3B). Comparison of the anomaly velocity structure and RT delamination structure of the 20 and 24 m.y. time steps (Fig. 6) suggests that the deeper portions of the anomaly are primarily thermal in origin, and the strongest compositional components are concentrated in the shallower core areas of the anomaly.

SUMMARY AND CONCLUSIONS

In this paper, we have developed further the thermomechanical modeling (initiated in Le Pourhiet et al., 2006) of mantle lithosphere removal in the southern Sierra Nevada region. We refined the initial models primarily by applying compositional and crustal structure constraints for the parameterization of the bulk rheology of the Great Valley crust (after Saleeby, 2007, 2012). As noted in Le Pourhiet et al. (2006), the Great Valley crustal rheology is of primary importance in the dynamics of delamination and in the generation of related vertical displacements. In another paper (Part II), we develop a database and derivative graphic formats of rock uplift and tectonic subsidence for the region for comparative purposes with model results of the same (presented here). We have shown that the observational data constraints and model results agree reasonably well, after adopting a relatively weak rheology for the Great Valley crust.

The model results and observational data suggest that a complex array of positive and negative epeirogenic transients migrated across the southern Sierra Nevada region in late Neogene–Quaternary time, in response to initial mantle lithosphere mobilization, and then there was geologically rapid (<5 m.y.) delamination of the arclogite root that developed in the Cretaceous beneath the Sierra Nevada batholith. The vertical displacement patterns are complex (Figs. 11 and 13) because crustal flexure and upper mantle–lower crustal buoyancy changes through the course of mobilization and removal are out of phase. The clearest example is exhibited in the 5 and 10 m.y. model steps (Fig. 6), where the thickening of the lower crustal channel beneath the batholith bows the elasto-plastic crust upward while the underlying arclogite root sags downward, and the deeper mantle lithosphere undergoes RT-type sinking. More subtle interplays between crustal and mantle buoyancy forces and crustal flexure occur in later time steps, resulting in the inflections in the vertical displacement curves of Figure 11. The migration of positive and negative epeirogenic transients across the southern Sierra Nevada is at odds with common debates regarding Sierra Nevada relief generation, which typically focus on unidirectional changes in regional relief through time, either positive or negative. In that analogous transients affected the southern Great Valley as well, an added complexity is posed to basin analysis, which is typically based primarily on eustatic and plate tectonic forcing.

Epeirogenic transients of the southern Sierra Nevada region are not axisymmetric about the Isabella anomaly. The anomaly is asymmetric, sympathetic with the asymmetry of the epeirogenic transients. These relations, among other factors, suggest viscous slab delamination over RT instability behavior as the principal driving mechanism for the observed transients. Nevertheless, RT instability behavior of the lithosphere as a whole is found to be a critical initiating factor in the delamination of the arclogite root, with RT mobilization driven by the high thermal gradient imparted below from the slab window, and the initial negative buoyancy perturbation of the arclogite root. Other factors that are found to be critical for this specific case are the development of a low-viscosity lower crustal channel along the locus of root delamination, and structural weakening of the upper mantle lithosphere and the crust with the development of regional ductile shear zones.

The modeling presented here and in Le Pourhiet et al. (2006) pertain specifically to initial conditions and kinematic and/or dynamic relations along transverse (east-northeast–west-southwest) sections across the Sierra Nevada batholith and its underlying mantle lithosphere. Seismic imaging of the Isabella anomaly and surface mapping of mantle lithosphere removal–related geological features show, however, that the southern Sierra Nevada case is complex in three dimensions. Nevertheless, the range of conditions encompassed by our various model iterations guided by the light of surface geological observations offers qualitative insights into how mantle lithosphere removal could, in this case, proceed in three dimensions.

This research was supported by National Science Foundation grant EAR-0606903, a grant from the George and Betty Moore Foundation, and Faculty Research Funds from the Université Pierre et Marie Curie (France). This is Caltech Tectonics Observatory Contribution 141. Communications with Craig Jones, Jeff Unruh, Mihai Ducea, George Zandt, Don Helmberger, Cyn-Ty Lee, Peter Molnar, Gene Humphreys, and Don Anderson helped stimulate our work.