We present quantitative estimates of near-surface rock strength relevant to landscape evolution and landslide hazard assessment for 15 geologic map units of the Longmen Shan, China. Strength estimates are derived from a novel method that inverts earthquake peak ground acceleration models and coseismic landslide inventories to obtain material properties and landslide thickness. Aggregate rock strength is determined by prescribing a friction angle of 30° and solving for effective cohesion. Effective cohesion ranges are from 70 kPa to 107 kPa for 15 geologic map units, and are approximately an order of magnitude less than typical laboratory measurements, probably because laboratory tests on hand-sized specimens do not incorporate the effects of heterogeneity and fracturing that likely control near-surface strength at the hillslope scale. We find that strength among the geologic map units studied varies by less than a factor of two. However, increased weakening of units with proximity to the range front, where precipitation and active fault density are the greatest, suggests that climatic and tectonic factors overwhelm lithologic differences in rock strength in this high-relief tectonically active setting.

Large earthquakes in mountainous regions commonly result in catastrophic and widespread landsliding, which is a significant hazard during and in the immediate aftermath of ground shaking (Keefer, 1994). In recent years, unprecedented geomorphic data sets of large coseismic landsliding events have helped refine our understanding of the relationship between strong ground motion and mass wasting (Meunier et al., 2007, 2013; Hovius et al., 2011; Parker et al., 2011). With these developments, the opportunity exists to exploit large earthquakes and the consequent landscape response to quantify near-surface rock strength under dynamic conditions at the hillslope scale (103–106 m2). From this perspective, the earthquake provides a regional shear-strength test where the earthquake imparts a measurable forcing (strong ground motion) and an observable response (slope failure) that is dependent on material strength.

The material strength of rock sets limits on topographic relief and modulates feedbacks between surface processes, tectonics, and climate (Schmidt and Montgomery, 1995; Burbank et al., 1996; Montgomery and Brandon, 2002). The strength of upper crustal earth materials generally is described using Mohr-Coulomb failure criteria, where shear strength is defined as the shear stress at failure, which is equal to the cohesion plus the product of the normal stress and the tangent of the angle of internal friction. Cohesion arises from the physical and chemical bonds in materials and friction describes the forces resisting motion between particles. The aggregate shear strength of soil or rock is described by the additive contributions of cohesive and frictional strength. Fractures, joints, and other discontinuities set the upper limit of rock strength, where their integrated effect acts to weaken rock at increasing spatial scales (Schmidt and Montgomery, 1995; Hoek and Brown, 1980). Thus, natural hillslope strength arises from the interplay of lithology, tectonic history, and climate.

Standard laboratory samples are too small to incorporate the scale and aggregate effect of strength-limiting discontinuities, and laboratory testing of large specimens that capture the true spatial heterogeneity of rock is not practical. These limitations have forced researchers to rely on qualitative measures to extrapolate laboratory measurements or use ranked classification schemes based on field observations to approximate rock strength at large spatial scales (Selby, 1980; Hoek and Brown, 1997). Although the Culmann (1875) criterion has been proposed as a way to estimate bedrock strength from hillslope height, this might apply only to deep-seated landsliding, and might not represent the integrated strength of the soil, regolith, and bedrock that characterizes many near-surface environments (Schmidt and Montgomery, 1995). Despite long recognition of these limitations, it remains notoriously difficult to quantify near-surface rock strength at geomorphically meaningful scales (Hoek and Brown, 1980, 1997), yet doing so is required to further our understanding of the role of rock strength in setting the pace of landscape evolution and to improve the accuracy of landslide hazard assessments.

Here we show how earthquake-triggered landslides and modeled peak ground acceleration (PGA) can be inverted to quantitatively document spatial trends in hillslope-scale rock strength using an example from the A.D. 2008 Mw 7.9 Wenchuan earthquake in Sichuan Province, China. The availability of high-resolution geologic and geomorphic data sets and geophysical models of seismic ground motion makes the Wenchuan earthquake the ideal natural laboratory to demonstrate the viability of using earthquake-triggered landslides to determine rock strength (Fig. 1). The earthquake occurred along the range front of the Longmen Shan and triggered >56,000 coseismic landslides over an area of >40,000 km2 (Dai et al., 2011). Coseismic landslides occurred in an array of different rock types and topographic conditions affected by a range of ground accelerations (Fig. 1; Ministry of Geology and Mineral Resources, 1991; Burchfiel et al., 1995; Worden et al., 2010). These slope failures also contributed to a large number of fatalities, and have had long-lasting societal and ecological impacts as rivers slowly transport material from affected areas (Huang and Fan, 2013).

We develop a mechanistic coseismic landslide model adapted from the Newmark (1965) sliding-block model commonly used for assessment of earthquake-induced landslide potential given a particular ground-shaking scenario (Wieczorek et al., 1985). The determination of strength parameters is relatively straightforward (Fig. DR1 in the GSA Data Repository1). Initial points of failure are identified in the digital elevation model (DEM) from the Newmark analysis, and depend on topographic slope and material properties; a synthetic landslide inventory is generated from these failure points using prescribed failure geometry. We perform a brute-force search to solve for the best-fit strength parameters that reproduce the frequency-area statistics of the mapped landslide inventory. In addition, we compare model results to total landslide volume and spatial distributions of landslide area and frequency as further model validation (for a full model description, see the Data Repository).

Prediction of hillslope failure is modeled using an infinite-slope stability analysis in which shear strength can be described using the angle of internal friction and cohesion. Coseismic surface displacements are calculated based on a critical (or yield) acceleration necessary to initiate displacement, from which material strength and topographic slope determine landslide susceptibility. Newmark displacements consider only rigid-block movements (i.e., no internal deformation) and plastic deformation at the base of the displaced mass. While other approaches consider internal deformation and elevated dynamic stresses that arise during movement, the Newmark analysis has been argued to be appropriate for regional-scale studies of coseismic landslides (Jibson et al., 2000; Godt et al., 2008).

Synthetic landslides are generated by first identifying target grid cells that exceed a threshold surface displacement (Newmark displacement of ≥5 cm) considering the critical acceleration necessary to initiate slope failure (Wieczorek et al., 1985; Jibson et al., 2000; Godt et al., 2008). Target cell displacements depend on local topographic slope (Shuttle Radar Topography mission DEM; ∼90 m resolution) and material strength properties, which are calculated for each geologic map unit (age and lithology). We consider a range of reasonable material strength values for each geologic map unit based on the literature (Schmidt and Montgomery, 1995; Hoek and Brown, 1980, 1997; Jibson et al., 2000).

Comparison of grid-based simulations with an observed landslide inventory requires some scheme to aggregate the one-dimensional (1-D) slope-stability solution to determine landslide area and volume. We apply geometric criteria to predict landslide area and volume for each target cell; this effectively expands the 1-D approximation into a 3-D synthetic landslide inventory. We considered a suite of landslide geometries constrained by different topographic boundary conditions, such as slope of the target cell and upslope distance to the nearest divide (Fig. 2A). All landslide geometries are modeled to extend a failure plane upslope from the target cell, which excludes consideration of transport of material and run out; in this way failure of the target cell acts to destabilize upslope cells beneath a failure plane. This is analogous to undercutting a slope and is an approach similar to that used to model deep-seated bedrock failures (Densmore et al., 1998).

The power-law exponent of the modeled frequency-area statistics is sensitive to the choice of geometric criteria to define the 3-D geometry of the landslide (Figs. 2A and 2B). We find that the best-fit model considers a concave-upward failure plane (i.e., bowl shaped) that extends from the target cell upslope until it reaches the surface, matching not only the power-law exponent but also the total landslide volume (Figs. 2A and 2B). The advantage of this geometry is that we obtain landslide depths and volumes independent of empirical scaling relationships (Hovius et al., 1997; Larsen et al., 2010); this enables comparison of total model volume to that of the measured inventory volume.

The frequency-area power-law scaling constant is determined by material strength and landslide thickness (Fig. 2C). This finding provides a means to quantify hillslope-scale rock strength by using a brute-force search to find the strength parameters for each geologic map unit to minimize χ2 (goodness-of-fit) between the observed and modeled landslide frequency-area distributions. Cohesion and thickness are treated as a ratio to avoid making a priori assumptions about landslide thickness.

Because the cohesion:thickness ratio and friction angle affect the power-law scaling constant similarly, no unique combination of these parameters minimizes χ2 (Fig. DR2; see the Data Repository). We use a constant friction angle of 30° (approximately the average modal hillslope angle; Fig. DR3) and solve for an effective cohesion. While the actual partitioning of friction and cohesion may differ from our calculations, the aggregate rock strength determined from the combination of these terms is the same as natural rock strength. Thus, effective cohesion normalized to a fixed fiction angle reflects the true relative differences in total rock strength.

We find statistically meaningful fits for the 15 geologic map units that we use to model the regional distribution of coseismic landslides (Figs. 2D and 3). Mean effective cohesion varies between 70 and 107 kPa, within the bounds of values estimated for hillslope-scale strength (Schmidt and Montgomery, 1995), and is roughly an order of magnitude lower than typical laboratory measurements (Fig. 4A) (Hoek and Brown, 1980, 1997). When we set cohesion to zero and solve for a best-fit effective friction angle, the same relative pattern in rock strength is found, and mean (±1σ) values are 50° ± 4° (Fig. DR4). These friction angles would be considered high even for intact rock strength (Hoek and Brown, 1997), suggesting that cohesion is an important component of hillslope-scale strength. This finding challenges the notion that rock strength at large scales is almost entirely frictional (cf. Burbank et al., 1996; Korup, 2008).

Due to the resolution of the digital topography (∼90 m) and PGA maps used in our study, we do not attempt to reproduce the extent and location of individual landslides; rather, our goal is to mimic drainage basin–scale to regional-scale statistics and spatial trends. The spatial pattern of our best-fit model matches the observed landslide inventory with mean misfits (±1σ) of 1.0% ± 10.2% and −3.6% ± 10.7% for landslide point and area densities, respectively (Fig. 3; Figs. DR5 and DR6). Principle component analysis shows that the model misfits exhibit no systematic correlation to the model inputs (Table DR1). Because we limit the observed inventory to the resolution of our model space (15,000–3,000,000 m2), it is reasonable to assume that errors associated with poorly resolved, small landslides within the study area are minimized (e.g., Li et al., 2014). Model misfits are likely a result of preforming an inversion at a coarse resolution (individual geologic map units cover areas between ∼400 and 22,000 km2) that likely averages local variations in rock strength. In addition, the PGA map used in our inversion might underestimate ground motions in the northeastern potion of the fault rupture where displacements were largely strike slip (Shen et al., 2009), and it does not take into account topographic site effects that can amplify seismic waves and cause local peaks in PGA and landslide concentrations (Meunier et al., 2008).

Model mean landslide thickness ranges from 6.8 m to 8.6 m and the total landslide volume is ∼4–7 km3, on par with empirically derived estimates of 6–13 m depth and 3–7 km3 volume estimated from the mapped inventory over the restricted model space (Fig. DR7) (Parker et al., 2011; Li et al., 2014; Larsen et al., 2010). The range of depths for observed and modeled landslides (one to tens of meters) indicates that failure occurred in a range of materials, from soil to weathered and fresh bedrock. We interpret our calculations as reflecting the strength of regolith and weathered or fractured bedrock that characterizes many near-surface environments (Anderson et al., 2013).

No obvious trend is observed between lithology and shear strength, implying that lithology plays little role in dictating the spatial variability of near-surface rock strength in the Longmen Shan (Fig. 4A). However, a positive correlation is noted between effective cohesion and distance from the range front (Figs. 4B and 4C). Both active fault density and the mean annual precipitation increase approaching the range front (Fig. 4D) (Burchfiel et al., 1995; Hijmans et al., 2005), so we speculate that in this humid, high-relief, tectonically active setting, tectonic and climatic history overwhelm intrinsic material properties and set limits on near-surface rock strength. If this inference is correct, reduced rock strength approaching the range front might be common in active orogenic belts owing to increased deformation toward the mountain front and the development of orographic precipitation, which would promote more rapid weathering and increase saturated conditions nearest the fault zone.

The low shear strength of rocks determined here underscores the difficulty of extrapolating laboratory strength tests to geomorphic problems. Because a single geotechnical measure or proxy is unlikely to fully capture the integrated strength at hillslope scales, we suggest that using extreme events that exert a large and measurable forcing on the landscape is a way forward to understand rock strength in the near-surface environment. We anticipate that rock strength, influenced by the effects of tectonic fracturing and climate, will quantitatively relate to erodibility because reduced rock strength suggests that the substrate is more easily detached and mobilized into the overlying soil column and available to transport by gravity, rivers, and glaciers (Molnar et al., 2007; Anderson et al., 2013). Advances in our ability to quantify hillslope-scale rock strength will also improve the implementation of landslide hazard models, thereby supporting efforts to reduce fatalities associated with seismically induced landslides.

This work is supported by the University of Michigan Department of Earth and Environmental Sciences, Turner Postdoctoral Fellowship to Gallen; a visiting faculty CIRES (Cooperative Institute for Research in Environmental Sciences) fellowship from the University of Colorado to Clark; and in part through computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor. We thank F.C. Dai for providing access to the landslide inventory used in this study, and Isaac Larsen, Randall Jibson, and two anonymous reviewers for comments that improved the manuscript.

1GSA Data Repository item 2015016, Table DR1 and Figures DR1–DR7, is available online at www.geosociety.org/pubs/ft2015.htm, or on request from editing@geosociety.org or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.