Abstract

The magnitude of late Cenozoic rock uplift of the Sierra Nevada (California, USA) remains unresolved despite more than a century of investigation, with estimates ranging from essentially zero to ∼3 km of uplift at the range crest. Two sets of two-dimensional end-member mechanical models bracket how normal faulting along the eastern escarpment of the Sierra Nevada contributed to uplift of the range over a time span of millions of years. The short-term models are based on dislocations in an elastic half-space. The long-term models involve thin elastic beams resting on an inviscid fluid. Both sets of models predict that if the regional topography were entirely a response to faulting along the eastern escarpment, then the bedrock floors immediately east of the range should consistently lie thousands of meters below sea level, instead of thousands of meters above sea level as they generally do. Both sets of analyses indicate that although faulting would lift the range crest, it would drop the rock east of the range-front faults at least as much, and perhaps much more; model results suggest that ∼66%–85% of the current escarpment relief stems from subsidence of the grabens east of the Sierra Nevada, with only ∼15%–34% resulting from crestal uplift. Our results strongly indicate that range-front faulting in the last 3–10 m.y. uplifted rock at the Sierra Nevada crest by hundreds of meters to as much as 1 km, and that this uplift was superposed on high topography that predated the origin of the eastern escarpment. These conclusions are compatible with diverse geologic observations and measurements.

INTRODUCTION

The Sierra Nevada (California, USA) is one of the most impressive mountain ranges in the world, renowned for its scenic beauty, size, and height (e.g., Whitney, 1869; Muir, 1912; Matthes, 1930). It is the longest single range in the contiguous United States and contains the tallest mountain, Mount Whitney (4421 m). Its eastern escarpment, with locally as much as 3 km of topographic relief, is one of the most impressive on Earth (Fig. 1). The imposing nature of the range has motivated geologists for more than a century to understand its origin.

The magnitude and timing of Sierra Nevada uplift are debated even after more than a century of investigation (Henry, 2009; Putirka and Busby, 2011; Jones and Saleeby, 2013). Some advocate that the range underwent substantial uplift (∼500 m to more than 2 km) by westward tilting in the last 3–10 m.y. (e.g., Christensen, 1966; Huber, 1981; Unruh, 1991; Wakabayashi and Sawyer, 2001; Jones et al., 2004; Stock et al., 2004, 2005; McPhillips and Brandon, 2012; Wakabayashi, 2013). This view is based on a variety of stratigraphic and geomorphic data, including inclinations of bedding, analyses of rigid-body rotations (tilts), stream channel gradients, and incision rates, and exhumation rates derived from thermochronometric data. In contrast, others contend that the range has been high for perhaps as long as 40 m.y., and might have undergone little to no late Cenozoic uplift (e.g., Wernicke et al., 1996; DeCelles, 2004; Mulch et al., 2006; Cassel et al., 2009, 2012; Hren et al., 2010). That view is based on differing interpretations of the stratigraphic and geomorphic data, as well as paleobotanical data and isotopic analyses of meteoric water preserved in weathered volcanic ash and fluvial sediments. Although a clear picture of the late Cenozoic uplift history of the Sierra Nevada has yet to emerge, much of the evidence does not necessarily contradict either view; i.e., the range could have retained high elevations since the Late Cretaceous but have undergone additional uplift and westward tilting in the late Cenozoic (e.g., Wakabayashi and Sawyer, 2001; Stock et al., 2005; McPhillips and Brandon, 2012).

The contribution of faulting along the eastern escarpment to uplift of the Sierra Nevada is of fundamental importance to this debate. More than a century ago, Gilbert (1880) recognized that faulting along the eastern escarpment caused a tilting of the range to the west. Lindgren (1911, p. 41) subsequently argued that the relief along the eastern escarpment reflected a drop of the surface east of the Sierra Nevada, stating that “Nowhere can any evidence be found substantiating the theory that the fault scarp has been formed by an uplift of the western block.” Henry (2009, p. 575) noted that if the range was at least as high in late Miocene time as it is now, then “…late Miocene faulting on its eastern flank represents subsidence of the Basin and Range, rather than uplift of the Sierra Nevada.” Henry (2009, p. 576) also posed a key question, “Could uplift in part be relative to a subsiding Basin and Range?” We focus on that question here and attempt to answer it quantitatively.

Simple solutions from elasticity theory can quantify the absolute amounts of uplift and subsidence resulting from faulting along the eastern escarpment. This information is essential to assess the contribution that faulting there has made to the uplift of the Sierra Nevada. In contrast, geologic field measurements, which can provide valuable information on the relative displacement (i.e., slip) across faults, typically provide little or no direct information on the associated absolute displacements. Elasticity theory solutions provide insight into the absolute displacements that is currently difficult to obtain by any other means.

Previous mechanical analyses directed at uplift of the Sierra Nevada have focused primarily on isostatic effects and used thin plate theory. Chase and Wallace (1986) considered the uplift of the range as primarily, but not entirely, a late Cenozoic event driven by a buoyant crustal root emplaced with the Sierra Nevada batholith during Mesozoic time; they treated the Sierran crust as a 50-km-thick plate that was fractured at the San Andreas fault and at the Owens Valley fault zone east of the Sierra Nevada. An uplift mechanism involving a buoyant root has been called into question, however, because additional seismic data (e.g., Wernicke et al., 1996; Zandt et al., 2004) and petrologic data (e.g., Ducea and Saleeby, 1998; Saleeby et al., 2003) indicate a lack of a compensating root and seem to favor uplift being driven by the delamination of a dense crustal root beneath the Sierra Nevada. Small and Anderson (1995, p. 277) examined the mechanical contribution of erosion and sedimentation on the tilting and uplift of the Sierra Nevada and concluded that “…a large fraction of the primary evidence for uplift could be generated by the lithospheric response to coupled erosion of the Sierra Nevada and deposition in the adjacent Central Valley….” Thompson and Parsons (2009), after treating the Sierran crust as an elastic plate ∼15 km thick, concluded that footwall unloading cannot account for the entire elevation of the Sierran crest above sea level. They also noted that if range-front faulting started in an already elevated plateau, then a hybrid model of preexisting regional uplift and localized footwall unloading could account for the older and newer uplift phases suggested by the geologic record. These contributions still leave room for further illumination of how much faulting along its eastern escarpment has contributed to the uplift of the Sierra Nevada.

Our analyses here focus on how faulting along the eastern escarpment affects the Sierra Nevada and the bedrock immediately to the east. Our treatment builds upon the work cited above, as well as research on the mechanics of normal faulting by others begun decades ago (e.g., Heiskanen and Vening Meinesz, 1958; Stein and Barrientos, 1985; Chase and Wallace, 1986; King et al., 1988; Schultz and Lin, 2001). The simple mechanical analyses we employ in our theoretical treatment are based on elastic dislocation theory (e.g., Segall, 2010) and thin-beam theory (e.g., Hetenyi, 1961; Turcotte and Schubert, 2002). As we will show, the two sets of models yield some similar predictions, but the thin-beam models account better for the amount of slip along the eastern escarpment.

We begin by reviewing critical data on the generalized geology of the Sierra Nevada and the faults of the eastern escarpment. This introductory material leads to the crux of our contribution: a series of two-dimensional mechanical analyses that predict the absolute displacements associated with slip along surface-breaching normal faults. We focus in particular on how faulting would affect the relative and absolute amounts of uplift and subsidence across the faults, but also consider the effect of faulting on the profile of the range west of the crest as well as east of the range-front faults. We conclude with a discussion of how much range-front faulting probably has contributed to Sierra Nevada uplift over the past several million years.

GEOLOGIC SETTING

Certain aspects of the geology of the Sierra Nevada are particularly important for constraining our mechanical models. These are the geometry and general composition of the range, the topography, the thickness of the crust and effective elastic thickness of the lithosphere, the geometry of normal faults of the eastern Sierra, the amount of slip along the faults, and the age of the eastern escarpment.

General Geology of the Sierra Nevada

The Sierra Nevada is ∼640 km long and has an average width of ∼100–115 km. By width, we mean the distance from the range-front fault trace at the base of the eastern escarpment to locations in the western foothills where the elevation has dropped to 10% of the crest height (Figs. 2 and 3). The core of the range consists primarily of Mesozoic granitic rocks, which form one of the great batholiths in the world (Fig. 2). The granitic rocks intrude older sedimentary and metamorphic rocks. The Mesozoic and older basement rocks are overlain by Cenozoic volcanic and sedimentary rocks along the western flank and eastern escarpment (Bateman and Wahrhaftig, 1966). Although not shown in Figure 2, volcanic and volcaniclastic rocks mantle much of the granitic rock north of Yosemite National Park, and basement rock crops out amid volcanic rock and alluvium in many places near the base of the eastern escarpment. Some of the volcanic and volcaniclastic units provide markers for evaluating the amount of slip along the range-front faults.

West of the Sierra Nevada, the basement rocks project beneath the thick package of upper Mesozoic and Cenozoic sedimentary deposits of the Great Valley. In the San Joaquin Valley, the southern part of the Great Valley trough, this package thickens from ∼800 m in the north to more than 9000 m in the south; the total thickness of Neogene and Quaternary deposits locally exceeds 5 km (Bartow, 1991). Near the west end of profile D (in Fig. 3) the Neogene and Quaternary deposits are 1 km thick (Bartow, 1991).

The eastern boundary of the Sierra Nevada is structural and is marked by a system of normal faults that separates the Sierra Nevada from the westernmost grabens of the Basin and Range province. These grabens, which include the Lake Tahoe basin, Mono Lake basin, and Owens Valley, generally range in width from 10 to 30 km. Few faults are now considered active within the range.

Topography

The crest of the Sierra Nevada is one of the highest of any range in the United States. Between Lake Tahoe and Yosemite Valley, the crest generally reaches an elevation of ∼3 km, rising 1–2 km above the base of the escarpment (Fig. 3). Farther south, along the western margin of Owens Valley, the crest exceeds 4.2 km in elevation in several places, towering ∼3 km above the Owens Valley floor.

The range has long been known to have an asymmetric profile (e.g., Whitney, 1869; King, 1872; Lindgren, 1911). The average inclination of the western slope typically is 1°–2° (e.g., Fig. 3), whereas the average inclination of the eastern escarpment reaches 26° west of the Owens Valley near Bishop. The profiles of Figure 3, drawn across the central portion of the range between Lake Tahoe and Yosemite, show that elevation diminishes from ∼3 km near the crest to ∼10% of the crest height at a distance of ∼100–115 km west-southwest of the range-front fault traces. In contrast, the same profiles show that the range loses ∼1 km in elevation (∼33% of the elevation of the crest) over a distance of only ∼5–15 km to the east-northeast to the base of the eastern escarpment.

Thickness of the Crust and Effective Elastic Thickness of the Lithosphere

The most recent estimates of the thickness of the crust in the Sierra Nevada, as inferred from seismic wave velocity data, are generally in the range of 32–40 km, the maximum thickness being ∼44 km, with uncertainties of 1 km or less (Gilbert, 2012). In the western foothills the crustal thickness might reach 50 km (Frassetto et al., 2011), but these values have uncertainties of 5 km or more (Gilbert, 2012). Ostos and Park (2012) collected magnetotelluric data on a traverse across the central Sierra Nevada to gauge the crustal thickness. They inferred that the base of the batholith, and the base of the crust, is 30 km below sea level east of Yosemite National Park, and 45 km below sea level beneath the western Sierra Nevada. The magnetotelluric estimates thus are compatible with the seismic estimates of Gilbert (2012). We use 45 km as a maximum crustal thickness in our faulting models.

The effective elastic thickness of the lithosphere (Te) is used to describe the mechanical resistance of the lithosphere to bending and the wavelength at which bending effects die off (e.g., Turcotte and Schubert, 2002). This quantity cannot be measured directly; it arises in models for the bending of an elastic layer (Timoshenko and Woinowsky-Krieger, 1959), where it is coupled to the elastic parameters of Young’s modulus (E) and Poisson’s ratio (ν). Various studies have posited that, for the continents, the ratio of Te to the crustal thickness is likely to be close to unity (McKenzie and Fairhead, 1997), generally less than unity (Maggi et al., 2000), and generally greater than unity (Burov and Watts, 2006). In flexural analyses of the Sierra Nevada directed at various geologic issues, preferred values of Te range from 5 km (Granger and Stock, 2004), to ∼15 km (Thompson and Parsons, 2009), to 20–35 km (Small and Anderson, 1995), to as much as 50 km (Chase and Wallace, 1986). Values of E used range from 50 GPa (Chase and Wallace, 1986) to 100 GPa (Granger and Stock, 2004). Both Chase and Wallace (1986) and Granger and Stock (2004) used a value of 0.25 for ν. In a more recent study of the lithosphere in the western United States based on gravity measurements and topographic data, Lowry and Pérez-Gussinyé (2011) concluded that Te beneath the Sierra Nevada is between 5 and 15 km, but approaches 30 km in the Central Valley; they assumed E = 70 GPa (A. Lowry, 2013, personal commun). In light of the wide-ranging estimates cited herein, in our subsequent models involving beam flexure we consider values of Te between 5 km and 45 km, values of E between 50 GPa and 100 GPa, and ν = 0.25.

Geometry of Faults of the Eastern Escarpment

The normal faults of the eastern escarpment of the Sierra Nevada generally strike north-northwest, parallel to the range as a whole, but form a complex left-stepping echelon pattern (Lindgren, 1911; Knopf, 1918; Moore, 1963; Rinehart and Ross, 1964; Huber and Rinehart, 1965; Bateman, 1965; Kistler, 1966; Giusso, 1981; Armin and John, 1983; Stone et al., 2000). The range-front fault system typically is several kilometers wide and contains a series of normal faults, rather a single fault with a discontinuous trace (see also Fig. 2).

Typical dips of range-front faults are largely obscured owing to alluvial cover, but direct measurements as well as geophysical data consistently indicate that these faults dip steeply to the east in the bedrock. In Owens Valley, eastward dips for range-front faults have been recorded in several places: 60°–70° west of Bishop (Bateman, 1965); 78° west of the Big Pine volcanic field (Bateman, 1965); and 50°, 60°, and 80° west of Owens Lake (Moore, 1981; Stone et al., 2000). Based on geophysical surveys of Owen Valley, Pakiser et al. (1964) concluded that major faults bounding the valley had “very steep” dips, depicting the faults in cross sections as being vertical. Some normal faults near the range front dip to the west, but these tend to have trace lengths of a few kilometers or less. Most of the major range-front faults thus seem to dip between 50° to the east and 90° in bedrock. The reported range of fault dips in alluvial deposits is greater. Phillips and Majkowski (2011) calculated dips ranging from 26° to 87° using three-point solutions for points selected along fault scarps in alluvium at the base of the Sierra Nevada and the White Mountains. In addition, based on first-motion solutions of Waldhauser and Schaff (2008), Phillips and Majkowski (2011) postulated even shallower dips at depth. We have more confidence in the fault dip values that come directly from the bedrock measurements and geophysical data, but consider range-front fault dips as low as 20°.

Certain lavas, and basalts in particular, that erupted along or near range-front faults indicate that the faults, or their roots, might extend to great depth. For example, at several places in the Big Pine volcanic field, Quaternary basalts erupted along faults and/or are offset by faults (Moore, 1963; Bateman, 1965; Gillespie, 1982; Martel et al., 1987; Bierman et al., 1991). The basalt sources are generally considered to be at a depth of 45–50 km based on thermobarometric data and major and trace element analyses of the basalts (Wang et al., 2002; Mordick and Glazner, 2006). These depths correspond to the maximum thickness of the crust in the Sierra Nevada as inferred from seismic wave velocity data (30–40 km, Wernicke et al., 1996; 40 km, Das and Nolet, 1998; 35–42 km, Fliedner et al., 2000; 44 km, Gilbert, 2012). Near Sonora Pass, northwest of Mono Lake (Figs. 2 and 3), high-potassium latites and quartz latites erupted close to faults near the current crest of the range, and some of these lavas are offset by faults (Busby et al., 2008a, 2008b). Busby et al. (2008b) inferred that those faults penetrated deeply into the lithosphere. We interpret the association of deeply sourced lavas with the range-front faults to indicate that the faults could extend as deep as 45 km and perhaps penetrate the crust.

In their review of the 1872 Owens Valley earthquake, Hough and Hutton (2008) compiled data that help constrain the depth of faulting along the eastern escarpment. They expressed confidence in data indicating that most earthquake hypocenters at the southern end of the 1872 fault rupture were at depths less than 15 km, but with some at depths as great as 25 km. Based on a combination of geological field observations, macroseismic data, and instrumentally recorded seismicity, Hough and Hutton (2008) concluded the 1872 earthquake could have had a moment magnitude (Mw) as large as 7.9, which translates into a seismic moment of 7 × 1020 Nm (Hanks and Kanamori, 1979). The moment, in turn, is a product of the fault rupture area, the average slip, and the rock stiffness (Hanks and Kanamori, 1979). Field observations and measurements of Beanland and Clark (1994) allow the downdip extent of the 1872 fault rupture to be estimated. Using an average slip of 6 m, a rupture trace length of 100 km that equals the entire length of the fault zone as mapped by Beanland and Clark (1994), and assuming a rock shear modulus of 30 GPa, a vertical slip plane would need to extend to a depth of 40 km to yield a seismic moment of 7 × 1020 Nm. Hough and Hutton (2008) favored a rupture length of 130 km and a rupture depth of 20 km to yield their best estimate of a moment magnitude of 7.75. These considerations indicate that at least some faults along the eastern escarpment extend to depths of 20–40 km. We conclude that some faults penetrate much of the crust, and perhaps entirely through it.

Fault Slip

Slip along the normal faults of the eastern escarpment can be evaluated either by offset markers or from the topographic relief. The offset markers generally are late Cenozoic volcanic or volcaniclastic units in the northern half of the range. With two key exceptions, topographic relief has provided the primary means for estimating slip along the southern part of the range front.

At several locations along the northern range front, normal faults offset marker units vertically by 1 km or more. Near the headwaters of the Feather River (Fig. 2), vertical separations of as much as 1 km have been deduced for Mehrten Formation rocks with an age of ca. 5 Ma (Wakabayshi and Sawyer, 2000, 2001). Near Sonora Pass (Fig. 3), field relations indicate a vertical separation of at least 1.2 km on the Table Mountain Latite, with an age of 10.0 ± 0.2 Ma (Busby et al., 2008b). Near the headwaters of the San Joaquin River, west of Long Valley (Fig. 2), the vertical separation of andesitic and basaltic lavas with ages of 2.2–3.6 Ma is estimated as 1 km (Bailey, 1989; Wakabayashi and Sawyer, 2000). The current topographic relief from the crest of the range to the current basin floors to the east in these three areas is 1 km (Wakabayashi and Sawyer, 2000), 1.3 km, and 0.9 km (Bailey, 1989), respectively, within 10% of the respective vertical separation estimates. This coincidence should not be construed as indicating that faults have offset the modern topographic surface by an amount equal to the vertical separation of the Miocene units, but rather that the current topographic relief can be a good proxy for the net vertical component of post-Miocene slip along the range-front faults.

The topographic relief across the eastern escarpment provides a minimum value for the vertical component of slip across the range-front fault system if two conditions hold. First, the amount of vertical erosion of the rock dropped down east of the range-front faults must be no more than that along the range crest since the onset of faulting. This seems likely given the very low (≤∼10 m/m.y.) bedrock erosion rates determined both at the crest of the Sierra Nevada (Small et al., 1997) and at bedrock outcrops in the adjacent Owens Valley (Nichols et al., 2006), and in light of thermochronologic data indicating that the bedrock uplands of the southern Sierra Nevada eroded at rates of ∼40–60 m/m.y. over much of the Cenozoic (Clark et al., 2005). Second, the surface connecting rocks along the range crest to their offset counterparts east of the base of the escarpment must not have sloped appreciably to the east at the onset of faulting. This assumption also seems reasonable, so we infer that the topographic relief from the range crest to the base of the eastern escarpment provides a minimum estimate for the vertical component of slip. The topographic relief across the eastern escarpment is 1–2 km at the latitude of Echo Summit, a few kilometers south of Lake Tahoe, to Yosemite, increasing to 2.3–2.8 km along the western margin of Owens Valley. Where granitic rocks crop out amid alluvium and volcanic rocks near the base of the escarpment, as they do in many locations (e.g., Kistler, 1966; Bateman, 1965; Moore, 1981), the topographic relief across the escarpment might approximate the vertical component of fault slip.

Where bedrock is buried beneath alluvium, the vertical component of slip could be greater than the topographic relief. On the basis of gravity, magnetic, and seismic refraction data, Pakiser et al. (1964) concluded that the maximum relief on the pre-Tertiary bedrock surface now buried beneath Owens Valley is 5.8 ± 2.6 km, and they interpreted this as structural relief associated with slip on vertical or near-vertical faults. With a spacing of ∼15 km between the traces of range-bounding faults on the edges of the Owens Valley graben, this is one of the narrowest grabens bordering the Sierra Nevada. Boreholes near Mono Lake (Higgins et al., 1985) reveal complicated topography on the buried basement surface there, with as much as 2 km of relief over horizontal distances of several kilometers. Based on the borehole data and geophysical data, Higgins et al. (1985) concluded that some of the buried basement relief in the Mono Basin resulted from faulting. The available data thus point to substantial and complicated relief on the buried basement surface but neither permit the maximum amount of dip slip on the range-front faults to be well resolved nor indicate that basement rock generally is deeply buried along the base of the eastern escarpment.

Recent field research at two places along the southern part of the escarpment allows the topographic relief to be compared directly to the cumulative vertical component of slip across range-front faults. The first is at the western margin of Owens Valley near the Big Pine volcanic field (Fig. 2). A marker horizon there, the gently dipping (subhorizontal) roof of the leucogranitic Red Mountain pluton, is dropped down across several steep faults between the range crest and the head of the bajada at the base of the escarpment (Bartley et al., 2012). Although firm age control for the pluton is lacking, Bartley et al. (2012) considered it likely to be about as old as nearby leucogranitic plutons with ages of 161 ± 1 Ma and 164.5 ± 2.3 Ma. The field findings indicate that the cumulative vertical component of fault slip near the Big Pine volcanic field effectively matches the topographic relief across the eastern escarpment there. The second location is near Lone Pine. Based on tilt-corrected southern Sierra Nevada thermochronologic age-elevation relations, Ali et al. (2009) inferred that the granitic rock in the Alabama Hills underwent 2.6 km of vertical separation relative to the rock near Mount Whitney. This compares to 2.4 km of topographic relief between the summit of Mount Whitney and the head of the bajada at the base of the escarpment, and 3.0 km of topographic relief between the summit of Mount Whitney and the west side of the Alabama Hills. In both locations the current topographic relief from the crest of the range to the foot of the escarpment is within 15% of the vertical separation estimates. Accordingly, we infer that the vertical component of slip across the eastern escarpment is 1–2 km from near Lake Tahoe to Yosemite, and ∼3 km near Mount Whitney.

Escarpment Age and Fault Activity

Several factors indicate that faulting along the eastern escarpment began between ca. 3 Ma and ca. 10 Ma. In addition to the dated lavas mentioned here that are associated with faults, a series of basalts along the range front west of Bishop have ages of 11.78 ± 0.05 Ma and 3.40 ± 0.06 Ma; based on geomorphic data these basalts appear to have been offset by range-front faults beginning 2–5 Ma (Phillips et al., 2011). A basaltic lava southwest of Bishop, dated by Dalrymple (1963) as ca. 9.6 Ma, apparently flowed toward modern Owens Valley, indicating that a depression existed near Bishop at that time (Bateman, 1965). Also, sedimentary units from the Coso Formation at the south end of Owens Valley indicate that a shallow basin had formed there by ca. 6 Ma, with normal faulting occurring between 2.5 and 4 Ma (Bacon et al., 1982). Based on reconstructed topographic surfaces across Owens Valley, Jayko (2009) estimated the onset of faulting as ca.9 Ma. Geomorphic relationships in the Sierra Nevada indicate that accelerated stream incision and subsequent knickpoint creation occurred between 3 and 10 Ma (Wakabayashi and Sawyer, 2001; Stock et al., 2004, 2005; Clark et al., 2005; Wakabayashi, 2013). The geologic data collectively point to a period of enhanced faulting and volcanism along the eastern margin of the Sierra Nevada between ca. 3 and ca. 10 Ma. We interpret this as the interval when range-front faulting along the escarpment originated.

The range-front faults are considered active, offsetting the topographic surface and Quaternary features (e.g., moraines and stream channels) in many places. Offset markers have been carefully measured at several places along the base of the eastern escarpment (e.g., Le et al., 2007; Frankel et al., 2008; Rood et al., 2011), and Quaternary slip there is dominantly dip slip. In contrast, Quaternary slip appears to be dominantly right lateral along the Owens Valley fault zone (e.g., Lubetkin and Clark, 1988; Beanland and Clark, 1994; Le et al., 2007; Frankel et al., 2008), although locally slip appears to be pure dip slip (Martel et al., 1987). Faults bordering the Sierra Nevada on the east have generated several earthquakes larger than magnitude 6 in recorded history (Rinehart and Smith, 1982), including the 1872 Owens Valley earthquake (Mw = 7.8–7.9; Hough and Hutton, 2008). The modern interseismic uplift rate of the range is interpreted to be 1–2 mm/yr based on global positioning system and interferometric synthetic aperture radar measurements (Hammond et al., 2012).

MECHANICAL ANALYSES OF RANGE-FRONT FAULTING

The purpose of our mechanical analyses of range-front faulting is twofold. The first is to gain quantitative insight into the absolute values (and ratios) of maximum rock uplift and maximum subsidence across range-front normal faults. The second is to gauge the contribution that faulting is likely to have made to the elevation profile of the Sierra Nevada. We use the theoretical results to test a null hypothesis that the current regional topography of the Sierra Nevada and the bedrock to the east are entirely reflections of normal slip on the range-front faults.

We consider two sets of simple models (Fig. 4). Models of the first set are based on dislocations in an elastic half-space (Figs. 4A–4D), and those of the second set on thin elastic beams that overlie a static inviscid fluid (Figs. 4E, 4F). Within each set we treat single faults and symmetric grabens to help account for the structural variation along the eastern escarpment of the Sierra Nevada. The elastic half-space models might be considered as short term, where viscous relaxation (e.g., mantle flow) is not accounted for, and the thin-beam models might be considered as long term, where viscous relaxation is complete. By considering these two sets of models we strive to bracket the response of the Sierra Nevada to range-front faulting over a period of several million years. For the first set of models we consider two kinds of boundary conditions, one where fault slip is specified (Figs. 4A, 4C), and the other where stress boundary conditions are imposed and the slip is calculated in response to a loss of strength along a fault (Figs. 4B, 4D). We consider this range of boundary conditions to help draw useful conclusions from our short-term models. For the thin-beam models, which have been used for decades to gain insight into the long-term response of the lithosphere to faulting (Heiskanen and Vening Meinesz, 1958), we consider only slip boundary conditions.

To treat real fault systems with traces much longer than their downdip extents, we use quasi-static two-dimensional models. The reference frames for all the models have an origin at the surface, with the x-axis pointing to the right, the y-axis pointing up, and the z-axis pointing out of the page (Fig. 4). Displacements are confined to the x-y plane (displacements in the z-direction are zero), so the models are for plane strain. A model fault extends an infinite distance along strike, parallel to the z-axis, but a finite distance down from the surface. We model a real range-front fault system as a single feature of vanishing thickness and uniform orientation. To distinguish between a feature in a model and a real geologic structure, we use the terms model fault or dislocation, and fault, respectively. We assume that the elastic properties of the faulted media are homogeneous and isotropic for simplicity and befitting the nature of the granitic rocks in the Sierra Nevada. We do not account for thermal effects. The solutions herein describe displacements at the surface along profiles perpendicular to the strike of a model fault. Our focus is on the vertical component of displacement rather than the horizontal component given our interest in uplift and subsidence along the escarpment.

Elastic Half-Space Models

Single-Dislocation Models with Slip Boundary Conditions

Our first model contains a single dislocation of constant dip (δ) in an elastic half-space with a horizontal surface (Fig. 4A). The dip is less than 90° for a model fault that dips in the direction of the x-axis, and between 90° and 180° for one that dips the other way. The dislocation extends downdip from the surface a distance L. The boundary conditions along the dislocation are of specified uniform slip (s) along the dislocation, with the dislocation walls required to remain in contact and not interpenetrate. The relative displacement across the dislocation thus is purely parallel to dip. The amount of slip is constrained by the height of the Sierran escarpment and offset markers, as described above. Far from the model fault the horizontal and vertical displacements approach zero. Only the three model parameters δ, s, and L are required to evaluate displacements on the surface. The solutions for the displacements are independent of the elastic moduli of the material that contains the dislocation.

The solutions for the displacements associated with faulting are either taken from or adapted from solutions of Segall (2010); see Appendix 1 for details. The solutions are exact for a traction-free half-space with a horizontal surface. For an area with preexisting topography, however, the displacements calculated here for a horizontal surface would be a first-order perturbation to be superposed on the preexisting topography to obtain the approximate total elevation profile after faulting.

We begin by examining a vertical model fault that intersects the surface at the coordinate origin. The solution for the vertical displacement for this particularly simple case (Segall, 2010) is short, yet illuminates the character of the more general (and more complicated) solutions that follow. The vertical displacement (uy) and horizontal displacement (ux) are (Fig. 5): 
graphic
and 
graphic
where the sign function, sgn (x/L), returns the sign of x/L and yields a step in the vertical displacement across the model fault. These equations are functions of just the ratio x/L, not of both x and L independently. The displacements scale linearly with the slip (s), with positive values of s meaning that the positive wall of the fault (x/L > 0) moves down relative to the negative wall (x/L < 0). Both the horizontal and vertical displacements approach zero as x/L → ±∞. The vertical displacement profile (Fig. 5) exhibits (1) a discontinuity at the model fault that equals the slip; (2) an uplifted surface that is bent upward; and (3) a dropped surface that is bent downward. The amount of tilt along the surface varies. The faulted medium does not behave as two rigid blocks provided that L is finite. The vertical displacement is an odd function of x, so the ratio of the maximum uplift to the maximum subsidence is unity.

The parameter L provides the only length scale in this model. An inspection of Figure 5 shows that |uy| decays monotonically to 50% of its maximum at |x/L| ≈ 0.44. We call the distance x1/2 at which the uplift decays to half its predicted maximum the decay half-width. The distance x1/10 at which the uplift decays to one-tenth its predicted maximum is about 3 times greater: x1/10 /L ≈ 1.35.

In the more general case of a surface-breaking model fault with an arbitrary dip δ (δ ≠ 90°) and arbitrary surface trace location Ξ, the displacements (see Appendix 1) are: 
graphic
and 
graphic
where s is positive for normal faulting. The dimensionless term ζ is the equivalent of (x/L) in Equations 1 and 2; ζ is the scaled position of an observation point on the surface relative to the x-coordinate ξ of the lower tip of the model fault, with the depth (d) of the lower tip of the model fault (Fig. 4) being the scaling factor: 
graphic
Equations 3 and 4 address dips either in the direction of +x (0<δ<90°) or –x (90°<δ<180°). Both uy and ux scale linearly with the fault slip s but depend on L and δ in a nonlinear manner.

We now explore the sensitivity of the vertical displacement profiles to variation in model normal fault dip. Figure 6 shows uy/s at the surface for model fault dips of 20°–90°, with the slip s and model fault extent L being held constant. For a 90° dip, Equation 2 and Figure 6 show that the uplift on one wall is equal and opposite that of the downward displacement on the other wall. As the fault dip decreases, this symmetry dissolves. The vertical displacements on the uplifted footwall (x/L<0) decrease by more than 50% as the dip decreases from 90° to 50°, with a decrease of more than 90% as the dip decreases from 90° to 20°. The vertical displacements on the downdropped hanging wall, in contrast, change relatively little for x/L < 0.7 for dips of 50° to 90°. For dips of 20° and 30°, peak uplift on the hanging wall exceeds the peak uplift on the footwall. The vertical component of slip equals s sin δ and hence diminishes as the dip decreases, but for purely geometric reasons.

The change in the vertical displacement is depicted in a different way in Figure 7, where the vertical displacement values are scaled by the maximum uplift on the footwall (uy*). Four key points about the single-dislocation model stand out. First, fault slip results in both uplift of the footwall and subsidence of the hanging wall near the fault, but uplift occurs on the hanging wall for x/L > 1.3 for dips of 60° or less, and for x/L > 0.8 for dips of 30° or less. Second, the maximum drop of the hanging wall equals or exceeds the maximum uplift of the footwall by as much as a factor of 2.6 for plausible fault dips in the range of 50° to 90°, and by as much as a factor of 8 for dips of 20°. Third, as the fault is approached, the surface of the footwall bends upward, and the surface of the hanging wall bends downward for fault dips of 70° or greater. Fourth, the scaled uplift profiles for the footwall for the different dips collapse onto nearly the same curve. This means that the shape of the uplift profile on the footwall is not sensitive to the dip of the dislocation for dips of 20°–90°. Each scaled profile for footwall uplift decays monotonically to 50% of its maximum at |x/L| ≈ 0.44. Coseismic elevation changes consistent with these predictions accompanied normal faulting during the 1983 Borah Peak earthquake (Stein and Barrientos, 1985). This correspondence indicates that the model is geologically relevant. The footwall profiles for dips of 20° and 30° are quite unlike the topography east of Sierra Nevada, suggesting that faults of the eastern escarpment dip more steeply.

Single Frictionless-Fault Models with Stress Boundary Conditions

The frictionless-fault models (Fig. 4B) predict fault slip, whereas the single-dislocation models require the slip to be specified. The frictionless-fault models thus allow us to test whether the amount of slip in the single-dislocation models is physically reasonable, as well as to predict what the associated displacements of the half-space surface are.

In the frictionless-fault models, the model faults slip in response to an ambient regional stress field where the most compressive ambient stress is vertical. The ambient regional stress field prior to fault slip is taken as: 
graphic
 
graphic
and 
graphic
where ρ is the density of rock, g is gravitational acceleration, and k is the ratio of the ambient horizontal normal stress to the ambient vertical normal stress. The ambient vertical normal stress entirely reflects the weight of overburden. The value of y is negative in the subsurface (Fig. 4), and compressive stresses are negative. The walls of the fault sustain no shear traction and are required to remain in contact. A frictionless model fault will accommodate the maximum slip possible in a quasi-static model such as ours (Pollard and Segall, 1987).

The frictionless-fault models have seven parameters, the first six being ρ, E, ν, g, L, and δ, where E is Young’s modulus and ν is Poisson’s ratio. The choices of these six parameters are relatively straightforward. We assign values for the first three parameters based on the tabulations for granitic rocks by Hatheway and Kiersch (1990). We take ρ as 2.7 × 103 kg/m3, a standard value for granitic rocks. The tabulated values for E for granitic rocks range from 5 to 100 GPa. Over this range, the lower value of E = 5 GPa will maximize fault slip, because the slip in a frictionless-fault model is inversely proportional to E where an ambient regional stress field drives faulting (e.g., Pollard and Segall, 1987). The tabulated values of Poisson’s ratio (ν) for granitic rocks typically range from 0.2 to 0.3. We assume a standard value of ν = 0.25, but note that our results are relatively insensitive to small changes in ν. The value of g at the surface of the Earth is 9.8 m/s2. We consider faults that can extend as much as 45 km downdip, so this is the maximum value for L. We treat fault dips as shallow as 20° but pay greatest attention to faults with a dip of 80°–85°, because for some grabens considered herein, planar graben-bounding faults with dips of 80° or less would intersect if extended 45 km downdip, and vertical faults would be unable to slip under the ambient stress field of Equations 6–8.

The seventh model parameter is k in Equation 7. The value we consider is k = 1/3. This value is acceptable in terms of rock mechanics tests on laboratory samples and would yield no horizontal displacements far from a model fault (i.e., lateral confinement) for ν = 0.25 (e.g., Martel, 2000). The ambient shear stress parallel to any inclined fault plane in a two-dimensional analysis scales as the difference between the two principal stresses, which here are the ambient horizontal and vertical normal stresses. Thus, based on Equations 6–8, fault slip scales with 1–k. Over the interval 1/3≤k≤1, a choice of k = 1/3 maximizes fault slip, and k = 1 minimizes fault slip.

To model slip under the aforementioned conditions we employ a numerical method, the displacement discontinuity boundary element (BEM) method (e.g., Crouch and Starfield, 1983; Martel and Langley, 2006), because we know of no analytical solutions. In our BEM implementation, model faults are subdivided into a series of contiguous elements 100 m long. Interpreting our BEM solutions is more complicated than the single-dislocation solutions for two main reasons. First, like numerical solutions in general, BEM solutions do not provide a simple closed analytical solution such as Equations 1 and 2. Second, solving a boundary value problem for displacements, by any method, is inherently more complicated when some of the boundary conditions are in terms of stresses or tractions rather than just displacements. This greater complexity is reflected in the greater number of required model parameters.

Figure 8 shows vertical displacement profiles for dips of 20°–80° for frictionless model faults with L = 45 km, k = 1/3, E = 5 GPa, and ν = 0.25. The profiles show that frictionless single-fault models with dips of 40°–80° can account for a vertical component of slip (i.e., structural relief) of 1.9–3.1 km, commensurate with the observed relief of the eastern escarpment. The profiles for dips of 20° and 30°, in contrast, can only account for structural relief of 0.4 km and 1.1 km, respectively, and are inadequate to explain the escarpment relief. None of the profiles, however, can account for the elevation of the range crest: the peak model uplift is ∼1.2 km, well below the elevation of the Sierra Nevada crest. For a fault with a dip of 20°, even if L were 130 km and the fault extended to a depth of 45 km, the maximum uplift at the surface of the footwall would only be 0.4 km. Decreasing E would increase the slip and the peak uplift, but our choice of E was already at the low end of the experimental range, so the discrepancy between the model predictions and the elevation of the range is difficult to attribute to E being too high. Increasing k from 1/3 to 1 would decrease the slip and the faulting-induced displacements. Even frictionless model faults in an elastic half-space are unable to account for the height of the range for ambient stresses described by Equations 6–8.

Figure 9 shows scaled vertical displacement profiles, with the scaling factor being the maximum uplift along each unscaled profile, so the peak scaled values equal unity. A uniform change in E or k would not change the shapes of the curves in Figure 9, indicating that the scaled profiles are robust results for a frictionless-fault model. The scaled profiles reveal four key points, similar to those resulting from the single-dislocation models. First, fault slip results in both uplift of the footwall and drop of the hanging wall. Second, the drop of the hanging wall equals or exceeds the uplift of the footwall by a factor of at least 2.6 for fault dips in bedrock from 50° to 80°. Third, the surface of the footwall bends upward and the surface of the hanging wall bends downward as the fault is approached for fault dips of 70° or greater. Fourth, the scaled profiles on the uplifted footwall essentially collapse onto a single curve, indicating that the shapes of the uplift profiles are insensitive to the dip of model fault. Figure 9 shows the results for L = 45 km, but analogous results also arise if L were 25 km or 35 km. The decay half-width again is less than L. For L = 45 km, the dimensionless decay distance is |x/L| ≈ 0.7. Nearly the same result arises for frictionless faults if L were 25 km or 35 km: the dimensionless decay distance in those cases is |x/L| ≈ 0.6. For purposes of comparison, Figure 9 also shows the scaled profile from Figure 7 for a vertical fault modeled as a single dislocation of uniform slip. The uplift at the surface of the footwall decays more slowly for the frictionless faults than for the dislocation of uniform slip. This relatively small quantitative difference arises because the slip is uniform with depth for the dislocation models, whereas the slip is a bit greater at depth than at the surface for much of the downdip extent of frictionless model faults.

The last points we address regarding the frictionless-fault model are the maximum vertical component of slip and the maximum uplift at the surface. Table 1 shows model results for a sample of frictionless faults with steep dips, with a maximum value of L of 45 km and k = 1/3. Young’s modulus (E) is set to 5 GPa; this the smallest value we can defend based on laboratory measurements and so will maximize fault slip. The maximum vertical component of slip at the surface is 3.4 km (at a dip of 63°), with the maximum uplift of 1.4 km (at a dip of 65°). Models with steep planar single faults of plausible dimensions in a homogeneous elastic half-space thus conceivably could account for the escarpment relief, but not for the height of the range crest if the current regional topography of the Sierra Nevada were entirely a reflection of normal slip on the range-front faults.

Graben Models

In keeping with our use of simple models, we consider symmetric grabens bounded by with faults with equal amounts of slip and with a surface trace spacing W (Fig. 4). These models allow us to simulate not only the Sierra Nevada but also the grabens that border it on the east. These grabens typically are 10–30 km wide, and many are bounded on the east by ranges that approach or equal the height of the Sierra Nevada, such as the White and Inyo Mountains (Fig. 2). This provides partial justification for considering symmetric grabens. We consider two graben models that are counterparts to the single-fault models: a superposition model (Fig. 4C), and a mechanical interaction model (Fig. 4D).

Double-dislocation superposition models with slip boundary conditions. The first graben model simply superposes the solutions for two isolated dip-slip dislocations that dip toward each other (Fig. 4C). This superposition does not constitute a complete mechanical interaction between graben-bounding faults, but could be considered as a first step toward a mechanical interaction (Segall and Pollard, 1980). To help understand this model, consider Figure 10, which shows three solutions for vertical displacements, scaled by the slip s: one for an isolated vertical dislocation at x/L = −0.5 (long-dashed curve), a second for an isolated vertical dislocation of opposite polarity at x/L = +0.5 (short-dashed curve), and the third being a superposition solution for a graben (solid curve) that represents the sum of the displacements for the isolated dislocations. The superposition solution, obtained using Equation 1, yields uplifted graben flanks that are slightly lower in elevation than for a single dislocation. This is because the uplift on the flank of one fault is partially negated by the superposed drop of the other fault. The elevation of the graben floor is likewise lowered, the entire graben floor being below the base elevation of the flanks far from the graben. The graben floor is dropped everywhere more than the graben flanks are raised. The graben floor is much more uniform in elevation than the surface of the hanging wall in the single-dislocation models and thus accounts better for outcrops of basement rock over a narrow elevation range on the basin floors east of the Sierra Nevada. The graben model also has a convex crown that tilts toward the bounding faults. This tilt could explain why Mono Lake and Topaz Lake lap against the base of the Sierran escarpment rather than being sited more centrally in the basins they occupy.

Figure 11 shows the modeled vertical surface displacements for symmetric grabens bounded by vertical faults with 1 km of slip that extend to 45 km depth, with graben widths (W) of 10 km, 20 km, and 30 km. The W/L ratios of these model grabens are smaller than for the model graben of Figure 10 and presumably are more compatible with the grabens east of the Sierra Nevada. Figure 11 shows that for all three widths the maximum uplift is well below the current height of the range. In addition, as the model graben width decreases, the maximum elevation of the flanks and the elevation of the graben floor decrease as well, with the graben floor being many hundreds of meters below the base level far from the graben. This is in sharp contrast to profile D of Figure 3 through the Mono Lake basin, where the footwall surface approaches sea level far to the west of the escarpment but bedrock of the graben floor crops out at an elevation of ∼2000 m near Mono Lake. Similar situations occur elsewhere along the eastern escarpment. The discrepancy between the predicted and observed elevations of bedrock on valley floors east of the escarpment would only increase if the amount of slip were increased to yield a peak uplift that matched the current elevation of the range crest. These critical discrepancies imply that the null hypothesis is incorrect. The superposition solutions encourage us to conclude that the current elevations of basement floors immediately east of the Sierran crest, as well as the current elevation of the Sierra Nevada crest, reflect effects of range-front faulting superposed on preexisting high topography.

Although these graben models for an elastic half-space fail to explain the entire uplift of the Sierra Nevada, their predictions of how faulting contributes to rock uplift of the Sierra and subsidence to the east are worth exploring a bit more. The ratio of the maximum uplift to the maximum subsidence predicted by our superposition model can be explored easily because only three parameters are involved. Figure 12 shows how this ratio varies for values of L between 20 km and 45 km, and for dislocation dips of 90°, 80°, 70°, and 60°. The curves show the ratios for symmetric grabens with widths of 10 km, 20 km, and 30 km. The maximum uplift and maximum subsidence are both proportional to the slip s, so their ratio is independent of s. As graben width decreases and the dip of the bounding faults decreases, the maximum permitted value of L decreases in order to prevent the faults from intersecting; this is why Figure 12D has only the curve for a graben 30 km wide. By varying one of the three model parameters while holding the other two parameters fixed, we can see how each parameter affects the ratio of maximum uplift to maximum subsidence. First, each individual curve in Figure 12 shows that for a fixed graben width and fixed bounding-fault dip, the ratio decreases as L increases. Second, a comparison of curves for a fixed bounding-fault dip shows that the ratio decreases as graben width decreases. Third, a comparison of curves for different bounding fault dips but fixed graben width W (i.e., a comparison of curves of the same color) shows that the ratio decreases as the bounding-fault dip decreases. Vertical faults (Fig. 12A) thus give the greatest uplift to subsidence ratios. For all scenarios, however, the ratio is less than unity, so the maximum uplift is less than the maximum subsidence.

Graben models for mechanically interacting frictionless faults. We use BEM analyses of graben models with mechanically interacting frictionless faults to illuminate three central issues. First, the analyses can test whether graben-bounding normal faults that extend no more than 45 km downdip could reasonably accommodate enough slip at the surface to account for the height of the eastern escarpment given plausible stress conditions in the crust. Second, they can test whether graben-bounding normal faults could reasonably account for the entire height of the range. Third, they allow the mechanical interaction of normal faults to be evaluated; this could be important where the spacing of normal fault traces is less than the downdip extent of the faults (e.g., Pollard and Segall, 1987).

Figure 13 shows three vertical displacement profiles for model grabens with frictionless faults with a dip of 85° and L = 45 km in a homogeneous elastic half-space with E = 5 GPa. For grabens 10 km wide, dips of 80° or less do not permit faults to extend 45 km downdip without crossing. The profiles of Figure 13 show slip as great as ∼3000 m. This means that planar graben-bounding faults of plausible downdip extent theoretically could account for the relief of the eastern escarpment if E were 5 GPa. Even with a Young’s modulus that low, however, they cannot account for the elevation of the range crest: the peak model uplift is ∼400 m, well below the elevation of the Sierra Nevada crest. Most of the slip on the faults comes from subsidence of the graben floor rather than uplift of the graben flanks. As was the case for individual frictionless faults, decreasing E would increase the slip and the peak uplift, but our choice of E is already at the low end of the experimental range, so the discrepancy between the model predictions and the elevation of the range is difficult to attribute to E being too high. Increasing k from 1/3 to 1 would decrease the slip, so the faulting-induced displacements would be completely inadequate to account for the entire height of the range.

Figure 14 compares various model solutions for frictionless isolated single faults that dip at 85° to those for frictionless faults bounding grabens. The slip at the surface (Fig. 14A), the maximum uplift (Fig. 14B), and the maximum subsidence (14C) all increase as L increases from 20 km to 45 km. Interestingly, the slip and subsidence are greatest for the narrowest grabens (i.e., the most closely spaced faults), but the maximum uplift is greatest for isolated faults (i.e., faults with an infinite spacing). Steep graben-bounding faults slip thus would be more likely than steep isolated faults to account for the relief of the eastern escarpment. None of the graben models, however, can account for the height of the Sierra Nevada crest if its height were negligible prior to range-front faulting. Regardless of the height of the range prior to faulting, the maximum subsidence is greater for the grabens than for isolated faults (Fig. 14C), whereas the ratio of maximum uplift to maximum subsidence is less than unity (Fig. 14D).

Implications of Elastic Half-Space Models for Sierra Nevada Uplift

The elastic half-space models constitute a short-term end member for testing the null hypothesis that the current regional topography of the Sierra Nevada and the bedrock to the east is entirely a reflection of slip on the range-front faults. Three main results of the elastic half-space models stand out. First, if the regional topographic relief were negligible prior to range-front faulting and the crest were then lifted thousands of meters above sea level by faulting, then the graben floor should be thousands of meters below sea level at the base of the eastern escarpment. The maximum subsidence of the rock east of the escarpment would be at least as much as the maximum uplift of the rock of the Sierra Nevada. This is not consistent with the presence of basement rocks at elevations of 1300–2000 m above sea level in many places east of the escarpment (e.g., Kistler, 1966; Bateman, 1965; Moore, 1981). Second, although the frictionless fault models conceivably could account for the topographic relief along some of the escarpment based on the predicted vertical component of slip at the surface, none of them can account for elevations >1.5 km at the range crest. Third, each profile in Figure 3 shows that the elevation of the range decays to 50% of the elevation of the crest of the range ∼50–60 km west-southwest of the fault traces at the base of the eastern escarpment, whereas all our elastic half-space fault models with L≤45 km indicate that the elevations would decay to half the crest height at distances ≤45 km from the crest if the current regional topography of the range were entirely a response to uplift along the range-front faults. The short-term elastic half-space results indicate that the null hypothesis is incorrect, and therefore that the effects of range-front faulting are superposed on high topography that predated range-front faulting.

Thin-Beam Models

The actual response of the Sierra Nevada to faulting over millions of years could reflect mantle flow and/or viscous relaxation at the base of the crust. To help bracket this response, we now consider long-term models using thin faulted Euler-Bernoulli elastic beams resting on an inviscid fluid of uniform density ρm (Heiskanen and Vening Meinesz, 1958; Turcotte and Schubert, 2002). These simulate complete viscous relaxation below the lithosphere.

In the solutions here, the beams are initially along the x-axis and displacements are restricted to the vertical direction (y). We use the term beam to emphasize the two-dimensional nature of the model, but others use the term plate in their two-dimensional analyses (e.g., Turcotte and Schubert, 2002).

Single Fault Models

A fault in a thin beam of uniform effective elastic thickness Te can be simulated in two equivalent ways. The first is by imposing vertical displacements of opposite sign on the ends of a cut through the beam to produce a desired amount of slip. This is akin to the single dislocation models in an elastic half-space. The second way is to impose vertical shear forces of equal magnitude but opposite direction on the beam ends to cause the desired amount of slip, with the shear forces corresponding to the tractions associated with an average shear stress drop on the fault. We provide the results for both approaches.

We start by modeling a single fault as a vertical cut that completely divides an infinitely long beam into two adjacent semi-infinite beams (Fig. 4E). The left-side semi-infinite beam extends from x = 0 to x = –∞, and the right-side beam extends from x = 0+ to x = +∞. At the right end of the left-side beam, the vertical deflection is uy (x = 0) = uy* and the bending moment is zero. At the left end of the right-side beam, the vertical deflection is uy(x = 0+) = –uy* and the bending moment is zero. The slip s (i.e., vertical displacement discontinuity) thus equals 2uy*. We refer to the associated uniform vertical shearing forces (per unit depth) at the beam ends at x = 0± as ±V0/2 (e.g., Turcotte and Schubert, 2002). The magnitude of these forces balances the weight of all the fluid displaced by each semi-infinite beam, according to Archimedes principle. Thus, an isostatic force balance requires that the total uplift integrated along the left-side beam (the blue area in Fig. 4E) matches the total subsidence integrated along the right-side beam (the red area in Fig. 4E). The bending moment along a thin beam is proportional to the second derivative of the deflection with respect to x (Turcotte and Schubert, 2002), and because those bending moments are zero, the remaining two boundary conditions at the offset ends are 
graphic
Far from the model fault, the boundary conditions are that the deflection of the beam and its slope are zero, yielding the remaining four boundary conditions needed to solve the problem: uy(x→±∞) = 0 and 
graphic
For the boundary conditions described, the vertical deflection of the right-side beam (x>0) is (see Equations 3–140, 3–141, and 3–142 of Turcotte and Schubert, 2002): 
graphic
and the vertical deflection of the left-side beam (x < 0) is 
graphic
where V0/2 is positive, α is called the flexural parameter 
graphic
ρm is the density of the mantle (set to 3.3 × 103 kg/m3), and D is the flexural rigidity 
graphic
Equations 9 and 10 show that α provides a length scale dictating the spatial rate at which vertical displacements decay as well as a key wavelength. In the cases considered here, α generally is between 10 and 100 km.

Figure 15 shows scaled vertical displacement profiles for several models of thin faulted beams. The scaling factor is the maximum uplift along each unscaled profile, so the peak scaled values equal unity. The profile of the uplifted (left) beam shows a deflection of zero at a distance from the model fault of x0 = (π/2α) and the bottom of a trough at a distance xb = (3π/4)α, where the deflection is –0.067uy* (Turcotte and Schubert, 2002). The uplifted beam is concave up from the model fault to the bottom of the trough. On account of this concavity, the decay half-width, x1/2 = 0.54α = 0.34|x0|, is less than half the distance between the model fault (and the peak uplift) and the nearest point where the deflection is zero. Correspondingly, x1/10 = 1.22α = 0.78|x0| is <90% of the distance between the model fault and the nearest point where the deflection is zero. The similarity of the curves for E = 50 GPa and E = 100 GPa in Figure 15 shows that the vertical deflection is relatively insensitive to variation in E from 50 GPa to 100 GPa for the range of effective elastic thicknesses considered. The deflection profile of a downdropped right-side beam is identical in magnitude to that of the uplifted left-side beam, but opposite in sign.

The scaled profiles echo the first three of the four major points made about the elastic half-space models. First, fault slip results in both uplift of the footwall and drop of the hanging wall. Second, the drop of the hanging wall equals the uplift of the footwall. Third, the surface of the footwall bends upward and the surface of the hanging wall bends downward as the fault is approached.

The thin beam results differ from the elastic half-space models in some key ways. Whereas in an elastic half-space the decay half-width is less than the downdip extent of a model fault (L), the decay half-width for the thin beam models can be several kilometers greater than the depth to the base of the fault, which corresponds to the effective elastic thickness. As a result, the vertical deflections can decay more slowly with distance from the points of peak uplift for the thin beam models than for the elastic half-space models. In addition, the profiles of the thin beams (Fig. 15) exhibit the aforementioned troughs not seen in the single dislocation models (Fig. 6) or the frictionless fault models (Fig. 8).

Graben Models

A symmetric graben bounded by vertical faults can also be modeled using thin-beam theory. Consider a central beam of length W, which represents the graben floor, flanked by two semi-infinite beams that form a right-side and left-side pair (Fig. 4F). The left end of the right-side beam is at x = +W/2, and the right end of the left-side beam is at x = –W/2. Akin to the model for a single fault, the upward displacement of the ends of the semi-infinite beams at x = ±W/2 is uy*and no bending moments are applied there. The displacement and slope at the remote ends of the beams (x→±∞) are zero. The equation for the vertical displacement for the right-side beam is found by modifying Equation 9 by substituting (xW/2) for x. 
graphic
The equation for the vertical displacement for the left-side beam is found by modifying Equation 10 by substituting (x + W/2) for x. 
graphic
The central beam extends from x = –W/2 to x = +W/2. The symmetry of the problem requires equal vertical deflections of the two ends of the central beam, which are produced by downwardly directed forces (per unit depth) of magnitude V0/2 at each end of the central beam (recall that V0/2 also is the magnitude of the force that deflects the ends of the flanking semi-infinite beams upward). No bending moments are applied to the ends of the central beam. Hetenyi (1961) provided the static solution for the deflection of such a beam. His solution assumes that the beam rests on an elastic foundation that resists deflection by an amount proportional to the amount of deflection. That solution can be modified to solve for the case here where the weight of displaced inviscid fluid of uniform density provides the resistance to deflection by replacing the constant of proportionality, k, of Hetenyi (1964) with 4Dα–4. The solution for the deflection of the central beam (–W/2 < x < +W/2) on an inviscid fluid is: 
graphic
or 
graphic
The deflection for the central beam resembles that of a convex parabola, similar in shape to the graben floors of the elastic half-space models. This completes the deflection solution for a graben modeled by thin beams.
The thin-beam solution for a graben yields an isostatic force balance akin to the thin-beam solution for a single fault. Accordingly, in Figure 4F, the sum of the blue areas matches the red area. That relationship allows the mean vertical displacement (ūy) of the central beam to be found. The blue area beneath the uplifted beams can be equated to the area y of fluid displaced beneath the central beam (Fig. 4F), so: 
graphic
where the left-side and right-side integrals are solved by using the modified versions of Equations 10 and 9, respectively. The integrals in Equation 15 each equal uy*α/2 (Abramowitz and Stegun, 1972), so the mean vertical displacement of the central beam is 
graphic
This matches the mean deflection obtained from Equations 15.

Figure 16 shows vertical displacement profiles, scaled by the maximum uplift, for grabens modeled with thin beams; values of Te are 5 km and 45 km, and values of E are 50 GPa and 100 GPa. Poisson’s ratio is set to 0.25. We draw attention to two points. First, the graben floors subside well below the base level far from the grabens. Second, graben subsidence is sensitive to graben width, with the subsidence increasing markedly as the graben narrows. In addition, the models collectively show that as the flexural rigidity increases, either by an increase in E or Te , the graben flanks become less curved and broader, graben subsidence increases, and the graben floor flattens.

The ratio of peak uplift to maximum subsidence across the graben-bounding faults bears directly on the role of faulting on the formation of the eastern escarpment of the Sierra Nevada. This ratio is solved for by using Equation 15a, with xW/2: 
graphic
Figure 17 shows the modeled peak uplift/subsidence ratio as a function of beam thickness for graben widths of 10 km, 20 km, and 30 km for values of Te from 5 km to 45 km. All our models for Te < 15 km predict that a trough would develop within 100 km of the graben centers (e.g., Figs. 16A, 16B), yet no such trough is observed in the Sierra Nevada, so we are satisfied that regional values for Te are likely to exceed 15 km. For Te > 12 km, all the thin-beam models for grabens indicate that the maximum uplift of the graben flanks would be less than the maximum subsidence of the graben floor. The uplift/subsidence ratio decreases as effective elastic thickness increases and graben width decreases. The ratios in Figure 17 are generally lower than those for the grabens modeled by superposition for L = Te > 20 km (Fig. 12). For beams with an effective thickness of 35–45 km, the modeled uplift is about one-half to one-seventh the maximum subsidence of the graben floor.

Implications of Thin-Beam Models for Sierra Nevada Uplift

We now review the null hypothesis that the current regional topography of the Sierra Nevada and the bedrock to the east is entirely a reflection of slip on the range-front faults in light of the results from the thin-beam, long-term models. The three main results of the thin-beam models generally mirror those of the elastic half-space models, but paint a slightly more nuanced picture. First, if the regional topographic relief were negligible prior to range-front faulting and the crest were then lifted thousands of meters above sea level by faulting, then the graben floor should be thousands of meters below sea level at the base of the eastern escarpment. The maximum subsidence of the rock east of the escarpment would be at least as much as the maximum uplift of the rock of the Sierra Nevada for Te > 12 km (Figs. 15 and 17). These predictions are inconsistent with observations. Second, none of the thin-beam models for Te > 12 km that account for as much as 3 km of topographic relief along the escarpment yield elevations greater than 1.5 km at the range crest. Third, each profile in Figure 3 shows that the observed elevation of the range decays to 50% of the elevation of the crest of the range ∼50–60 km west-southwest of the fault traces at the base of the eastern escarpment, and to 10% of the crest height 100–115 km from the base of the escarpment. In contrast, the thin-beam models with Te ≤ 30 km predict that elevations should decay to half the crest height at distances no more than 40 km from the fault traces, and to 10% of the crest height no more than 90 km from the fault traces (Table 2). Thin-beam models with values of Te between 40 and 45 km and values of E between 50 and 100 GPa, however, could account for the two observed elevation profile decay characteristics of the western slope of the Sierra Nevada. In addition, the latter models also predict subsidence below sea level 120–170 km from the faults and thus could help explain the thickness of Neogene and Quaternary deposits in the Great Valley, which the elastic half-space models cannot. Still, the first two results from both the short-term elastic half-space models and the long-term thin-beam results indicate that the null hypothesis is incorrect. We therefore conclude that the effects of range-front faulting are superposed on high topography that predated range-front faulting.

DISCUSSION

Regarding the question posed by Henry (2009, p. 576), “Could uplift in part be relative to a subsiding Basin and Range?”, our mechanical analyses indicate that the answer is yes. Our model solutions indicate that slip along the range-front faults causes the east-side grabens to drop and the Sierra Nevada to rise. The graben models further indicate that the bedrock immediately east of the Sierra Nevada should have dropped more than the Sierra Nevada crest was uplifted, most likely at least twice as much and perhaps several times more. Given uncertainties in the downdip extent of the real faults, variations in the attitude and spacing of the faults, and variations in the slip along the faults, we expect that the ratio of uplift to subsidence is likely to vary along the Sierran escarpment. Nonetheless, our model results suggest that ∼66%–85% of the present escarpment height stems from dropping of the grabens east of the Sierra, with only ∼15%–34% of the present escarpment height resulting from uplift of the Sierra. Given that the topographic relief across the eastern escarpment in the central part of the range generally is 1–2 km, our results indicate that faulting along the eastern escarpment would likely yield a maximum uplift of 300–700 m at the range crest. The topographic relief across the eastern escarpment in the southern Sierra Nevada adjacent to the Owens Valley reaches 3 km, so crestal uplift arising from faulting there might be as much as 1 km.

The model results thus support a second conclusion that the range crest was within a few hundred meters to 1 km of its current height prior to the onset of slip along the range-front faults ca. 10 Ma. Apparently, the eastern escarpment developed in a region that was already relatively high. This conclusion is consistent with field-based studies suggesting that the Sierra Nevada and the areas immediately east of it retained high crestal elevations throughout the Cenozoic (e.g., Mulch et al., 2006; Cassel et al., 2009, 2012; Hren et al., 2010). It is also consistent with previous modeling work by Chase and Wallace (1986) and Thompson and Parsons (2009).

Our finding that the crest of the Sierra Nevada must have undergone at least some rock uplift associated with range-front faulting in the past ∼10 m.y. is also generally consistent with geomorphic evidence supporting uplift involving westward tilting (e.g., Huber, 1981; Unruh, 1991; Wakabayashi and Sawyer, 2001; Stock et al., 2004, 2005; Jones et al., 2004). Thus, our findings tend to support an emerging view that modern Sierra Nevada topography preserves elements of both Cretaceous uplift associated with arc volcanism and late Cenozoic uplift associated with faulting and westward tilting (e.g., Wakabayashi and Sawyer, 2001; Stock et al., 2005; McPhillips and Brandon, 2012). Our analyses suggest that the magnitude of late Cenozoic crestal uplift was not more than ∼1 km.

Our results have implications for how tilts of bedding bear on estimates of uplift of the range. Differences in dips between units of different age theoretically could give the tilting of the range if the units had the same original dip, the same dip direction, and occur the same distance from the traces of range-front fault (e.g., the units are in the same outcrop). These conditions, however, might be difficult to meet. If units of the same age and initial tilt crop out at decidedly different distances from the range-front faults, then the bending predicted by our models would result in different final tilts. Uplift analyses based on bedding tilts (e.g., Christensen, 1966; Unruh, 1991) typically assume that the range tilts uniformly as a rigid block, when in fact that is not likely to occur. Given the typical small slopes of the western Sierra Nevada, 1°–2°, tilt variations due to bending and variations in initial dips, as well as three-dimensional tilting effects not included in our two-dimensional model, could introduce large uncertainties in uplift estimates.

Our results have some interesting implications for graben development. Both the short-term elastic half-space models and most of the long-term thin-beam models predict that the maximum subsidence of rock in grabens would exceed the maximum uplift of the rock on the uplifted flanks. These theoretical predictions also are consistent with geodetic measurements of grabens that developed above dikes in rift zones of Iceland and Afar (Rubin and Pollard, 1988), suggesting that over a wide range of conditions the maximum subsidence in grabens would exceed the peak uplift on the flanks. In addition, the graben models consistently predict that the bedrock surface of the graben footwall would subside relatively uniformly, whereas the single-fault models predict a highly asymmetric footwall trench. The graben models appear to account much better for the presence of basement rocks that crop out several kilometers to the east of the base of the eastern escarpment of the Sierra Nevada.

CONCLUSIONS

Mechanical analyses indicate that the topographic relief of the eastern escarpment of the Sierra Nevada is due primarily to subsidence of the bedrock east of the range-front faults rather than to rock uplift of the Sierra Nevada. Uplift of the Sierra Nevada in response to slip along the range-front faults was probably no more than a few hundred meters to as much as 1 km at the crest, whereas the valleys to the east have probably been dropped down at least 1 km. Late Cenozoic faulting alone cannot account for the modern elevation and morphology of the range and the valleys and basins to the east. Our theoretical results strongly indicate that although faulting drove some crestal uplift, much of the present elevation of the Sierra Nevada predates formation of the eastern escarpment ca. 10 Ma.

We acknowledge Matt Barbee for his help with the initial handling of the topographic data. We thank Bob Burd for his help in identifying the peaks of Figure 1. We also thank Anthony Lowry for his comments on his research on effective elastic thickness of the lithosphere in the western United States, and Laura Wallace and Robert Anderson for discussions of flexure computer codes. Support from the National Science Foundation to Martel (grant EAR-05-38334) and the Yosemite Conservancy (Agreement CA#H80800090008) to Martel and Stock is gratefully acknowledged. This is SOEST (School of Ocean and Earth Science and Technology, University of Hawaii) contribution 8930.

APPENDIX 1. SOLUTIONS FOR A TWO-DIMENSIONAL PLANE-STRAIN DISPLACEMENT DISCONTINTUITY OF FINITE LENGTH

Segall (2010) derived the solutions for the displacement changes and strains changes at the surface (y = 0) that arise from slip along a displacement discontinuity that dips to the right (0° < δ < 90°) and intersect the surface at x = ξ (see text). These displacement and strain perturbations would need to be superposed on the displacement and strain fields prior to faulting to obtain total displacement and strain fields. The vertical displacement (uy) at the surface is (Segall, 2010, Equation 3.73; note that the term sgn(ζ) is the sign of ζ and corrects a typographical error in Equation 3.73): 
graphic
where the slip s is positive for normal faulting, and ζ is the normalized position of an observation point relative to the location of the fault trace, with the depth of the bottom of the displacement discontinuity being the normalization factor 
graphic
To obtain a more general solution in which faults can dip either to the left or right we invoke symmetry considerations. The vertical displacement at a point to the right (or left) of a fault that dips to the left is obtained by reflecting the solution for the vertical displacement at a point to the left (or right) of a fault that dips to the right: 
graphic
This constraint leads to the following solution for a fault that dips to the left: 
graphic
The solution for a fault that dips either to the right or the left is: 
graphic
The solution for the horizontal displacements for a normal fault that dips to the right is (Segall, 2010, Equation 3.73): 
graphic
with sgn(x) being replaced by sgn(ζ) to correct a typographical error. The horizontal displacement at a point to the right (or left) of a fault that dips to the left is obtained by reflecting the solution for the horizontal displacement at a point to the left (or right) of a fault that dips to the right: 
graphic
This constraint leads to the following solution for a fault that dips to the left: 
graphic
The corresponding equation for a fault that dips either to the left or the right, obtained using the procedure described here, is: 
graphic
The tilt is given by the derivative ∂uy/∂x. The solution follows from Equation A5: 
graphic

APPENDIX 2. BOUNDARY ELEMENT PROCEDURE

The general operation of the boundary element method we use was described lucidly by Crouch and Starfield (1983). The requirement we impose that the walls of the model fault remain in contact greatly simplifies how we implement the method. This is because only the shear tractions on the walls need be assigned, and only the shear (tangential) displacement discontinuities across the elements need to be solved for in the part of the method that involves an inverse solution.

If both the normal and shear tractions needed to be specified, the matrix equation that is solved to obtain the relative displacement distribution along a model fault would be of the following form: 
graphic
where the submatrices [A] contain influence coefficients that describe traction changes on an element due to a unit displacement discontinuity on another (or the same) element, the submatrices [Ds] and [Dn] contain the shear and normal displacement discontinuities (i.e., the relative displacements) across the boundary elements, and the submatrices [Bs] and [Bn] contain modified shear and normal traction boundary conditions on the boundary elements (Crouch and Starfield, 1983). The values of [A] depend on the boundary element geometry and the elastic moduli. Crouch and Starfield (1983) and Martel and Langley (2006) provided the formulas for calculating the influence coefficients. The modified shear and normal traction boundary conditions on the boundary elements are determined by the actual boundary conditions as corrected by the ambient stresses. The displacement discontinuities are solved for by an inverse method.
As we have formulated the problem, however, no normal traction boundary condition exists on the elements. This means that the second row of submatrices in Equation B1 can be eliminated. Also, because the walls of the elements remain in contact, all the values of [Dn] must equal zero, and hence the influence coefficient matrix [Asn] is immaterial. The only submatrices of Equation B1 that can contribute are: 
graphic
where [Ass] contains influence coefficients that describe shear traction changes on an element due to a unit shear displacement discontinuity on another (or the same) element. The values of [Ds] are solved for by an inverse method. The solution of Equation B2 takes much less time and computer storage than the solution of Equation B1.
Once the values of [Ds] are known, a new set of influence coefficient submatrices, [Axs] and [Ays], are calculated that describe the displacement changes at observation points due to unit shear displacement discontinuities at the elements. The following matrix is then solved for in a forward sense to yield the horizontal and vertical displacements at the observation points: 
graphic