Research ArticleGEOLOGY

Rainfall triggers more deep-seated landslides than Cascadia earthquakes in the Oregon Coast Range, USA

See allHide authors and affiliations

Science Advances  16 Sep 2020:
Vol. 6, no. 38, eaba6790
DOI: 10.1126/sciadv.aba6790

Abstract

The coastal Pacific Northwest USA hosts thousands of deep-seated landslides. Historic landslides have primarily been triggered by rainfall, but the region is also prone to large earthquakes on the 1100-km-long Cascadia Subduction Zone megathrust. Little is known about the number of landslides triggered by these earthquakes because the last magnitude 9 rupture occurred in 1700 CE. Here, we map 9938 deep-seated bedrock landslides in the Oregon Coast Range and use surface roughness dating to estimate that past earthquakes triggered fewer than half of the landslides in the past 1000 years. We find landslide frequency increases with mean annual precipitation but not with modeled peak ground acceleration or proximity to the megathrust. Our results agree with findings about other recent subduction zone earthquakes where relatively few deep-seated landslides were mapped and suggest that despite proximity to the megathrust, most deep-seated landslides in the Oregon Coast Range were triggered by rainfall.

INTRODUCTION

Landslides are one of the greatest secondary hazards of large earthquakes. Hillslope failures triggered by shaking can affect homes and infrastructure and dam rivers and cause flooding, sometimes in catastrophic outbursts (14). In the Pacific Northwest of the United States, magnitude 8 (M8) to M9 earthquakes occur every 300 to 500 years along the Cascadia Subduction Zone (CSZ) (5), a 1100-km-long megathrust fault that accommodates subduction of the Juan de Fuca Plate beneath the North American plate (Fig. 1). Empirical estimates for the number of triggered landslides for earthquakes of this magnitude range from 50,000 to 1,000,000+, with total predicted volumes of 10 to 110 km3 (6, 7), representing a risk to the millions of people living in the region. However, predictive models for earthquake-induced landslides from a potential CSZ event have high uncertainty and often fail to match the observations from similar subduction zone events. The few complete coseismic landslide inventories in subduction zones challenge assumptions about the number and volume of landslides triggered by these earthquakes.

Fig. 1 Tectonic setting and geography of study area.

Black outline depicts the central Oregon Coast Range (OCR) study area, which is underlain by the Eocene Tyee Formation and overlying Elkton Formation. The CSZ runs roughly north to south for over 1100 km just offshore. Shaded relief base map compiled from U.S. Geological Survey 30-m elevation data and National Oceanic and Atmospheric Administration bathymetry by the state of Oregon. Depth contours show the interface between the Juan de Fuca Plate and the overlying North American plate (27, 28).

Recent work has improved earthquake-landslide volume scaling relationships by accounting for rupture depth and topographic steepness; these models predict total coseismic landslide volumes of 0.6 to 2.0 km3 for moment magnitude (Mw) 9.0 CSZ earthquakes with rupture depths of between 20 and 10 km, respectively (8). Even these improved models overestimate the total coseismic landslide volume for recent subduction zone earthquakes by ~2 orders of magnitude, which may be an effect of the large distance between rupture hypocenters and subaerial landslide-prone topography. Only 1226 landslides were mapped following the 2010 Mw 8.8 Maule earthquake in Chile (9). The following year, the 2011 Mw 9.0 Tohoku earthquake struck Japan, causing peak ground accelerations (PGAs) over 1.0g at many sites but resulting in just 3477 mapped landslides (10). Most of the landslides triggered in both earthquakes were shallow, disrupted slides (defined here as landslides where the mobilized slide material mostly or entirely disintegrates during failure, resulting in flow or avalanche-style behavior) and lateral spreads. In Cascadia, the timing of the last great CSZ earthquake is well constrained to 1700 CE (11), yet no prehistoric landslides have been definitively dated to that time (12). Recent landslide dating efforts have precisely determined the age of several landslide-dammed lakes in western Oregon, but none date to 1700 CE (13). Complicating the search for past coseismic landslides in Cascadia is the characteristically wet climate of the Pacific Northwest, where heavy rainfall often triggers landslides (1416). Recent work in the Nepalese Himalaya, where landslides are triggered by earthquakes and monsoonal rain, suggests that rainfall is responsible for between 40 and 90% of the total landslide erosion on 104-year time scales (17). Because of the unique climate and tectonic setting of Cascadia, the challenge of parsing the relative importance of these triggers remains essential to understanding what controls landslide hazards and, on longer time scales, landscape denudation in the region.

One barrier to understanding the behavior and triggering mechanisms of past landslides is the difficulty of dating landslides on a regional scale. Landslide mapping inventories, however accurate, are of limited use in their ability to assess landslide triggering or event frequency because they often contain no age information for individual slope failure events. Recent work in Washington state has shown that for landslides in glacial sediment, surface roughness measured from lidar data can be used as a proxy for landslide age, enabling regional landslide chronologies to be constructed from a limited set of absolute dates (18, 19). This technique takes advantage of the topographic smoothing of landslide deposits over time as erosion and soil transport move material from convex areas to concave areas. Thus far, this technique has been applied to Quaternary glacial deposits (18, 19), and here, we demonstrate its effectiveness in dating bedrock landslides.

In this study, we apply surface roughness dating to an inventory of 9938 manually mapped, deep-seated bedrock landslides across 15,000 km2 in the sandstones and siltstones of the Tyee and Elkton Formations (20) of the central Oregon Coast Range (OCR) (Fig. 2). Landslide deposits were mapped from 0.9-m bare-earth lidar imagery based on the presence of an arcuate headscarp and adjacent hummocky topography (fig. S1). The OCR, a mountainous swath of western Cascadia, has a well-documented abundance of prehistoric and historic landslides and is prone to both deep-seated bedrock landslides and shallow translational slides and debris flows (2126). Thousands of landslides have previously been identified in the OCR (26), and automated landslide identification techniques suggest that up to 25% of the Tyee Formation is affected by deep-seated slides, underscoring the importance of this process to landscape evolution and morphology (24). Valley-spanning prehistoric landslides exist throughout the region, with remnant portions of deposits causing substantial upstream aggradation (22). Historic deep-seated landslides include a failure that occurred after days of heavy rain in 1975, damming Drift Creek and creating a lake that persists today (24).

Fig. 2 Landslide locations within the study area.

(A) Manually mapped deep-seated landslide deposits (n = 9938) from the study area (outline in Fig. 1) shown as black dots over a base map of elevation and shaded relief. Dated landslides shown as white circles with ID numbers (table S1). Inset map for (B) outlined in black. (B) Detailed map of landslide deposits and estimated age ranges. Deep-seated landslides shown as polygons draped over 0.9-m bare-earth lidar-derived hillshade, with estimated ages from deposit surface roughness identified by polygon color.

Numerical earthquake models suggest that the OCR, which overlies the CSZ and in some places is less than 100 km from the associated offshore trench, may experience PGAs >0.6g during M9 earthquakes (27, 28 Fig. 1). When these strong ground motions interact with steep slopes often exceeding 30° in our study area, many models would predict widespread coseismic landslide triggering [e.g., (29)]. However, the lack of historic records of coseismic landslides associated with the CSZ inhibits our understanding of where and how slopes fail during these large earthquakes. With the recent availability of <1-m-resolution bare-earth lidar coverage of the entire OCR, far more accurate and complete landslide mapping is now possible, as are surface roughness dating techniques that require high-resolution topographic data. Using 14 bedrock landslides with independent age constraints ranging from 7 ± 1 years to 41,800 ± 1000 years before present (ybp; here 2019) from this and other studies, we find an exponential relationship between landslide deposit surface roughness and age. We then use that calibrated roughness-age function to estimate the ages of all mapped landslides in the study area. With this landslide chronology, we first test whether predicted strong ground motions from CSZ earthquakes or mean annual rainfall exert control on the spatial distribution of deep-seated slope failures over the past 1000 years. We then examine temporal trends in landslide frequency over the same time span to determine the importance of past large earthquakes as landslide triggers.

RESULTS

Surface roughness dating of bedrock landslides

Previous work has tested the ability of numerous roughness metrics to both track known landslide ages and correctly identify relative ages among landslides with cross-cutting relationships (18, 19). For this study, we adopt the measure with the highest combined performance from these previous tests: the two-dimensional (2D) continuous wavelet transform (CWT) (19). We take the mean magnitude of the wavelet coefficient, which is equivalent to topographic curvature, over the entire mapped deposit for each of 14 landslides of known age to represent their roughness (table S1 and fig. S1). Similar to previous work in Washington state, the relationship between landslide deposit roughness and age is well fit by an exponential decay function (18, 19). We fit a robust (least absolute error) exponential regression using these roughness-age data, which yields the following function: t = 1.428 × 106 × e−412.8R, where t is the estimated landslide age, in years before 2019, and R is the average magnitude of the wavelet coefficient at a 20-m length scale within a mapped landslide polygon, in units of m−1 (Fig. 3). We use this function to estimate the ages and uncertainties of all 9938 landslides across the study area based on their measured roughness value (Fig. 2B). As a conservative estimate of uncertainty for any single estimated landslide age, we calculate prediction bounds at a 95% confidence level for all age values. Uncertainty in predicted ages increases with time. For landslides with historic estimated ages (less than 200 ybp), 95% prediction bounds yield an age range of 0 to 1300 years. This age range grows to 0 to 2400 years for landslides estimated to be 1000 years old. For landslides estimated to be thousands to tens of thousands of years old, uncertainty is plus or minus a few thousand years. While this technique is therefore not appropriate for pinpointing the triggering time for an individual slope failure, it offers a way to examine broadscale landsliding patterns in space and time across a landscape when applied to large landslide inventories.

Fig. 3 Plot of landslide age versus deposit surface roughness.

(A) Open circles represent ages of landslides with independent timing constraints from radiocarbon, dendrochronology, or repeat aerial imagery plotted against the surface roughness of those deposits, quantified by the average magnitude of the wavelet coefficient for each polygon using a CWT with a 20-m Mexican Hat wavelet. Exponential decay function fit to these data shown as red line, where t is the landslide age (years before 2019) and R is the landslide roughness (mean wavelet coefficient of mapped polygon). Functional bounds show 95% confidence of fitted regression, plotted as blue dashed lines. Prediction bounds show 95% confidence for landslide age prediction, plotted as solid black lines. (B) Same data points and regression plotted on a semilog axis to better illustrate variability in landslide ages for slides less than 1000 years old.

Spatial relationships between rainfall, CSZ earthquakes, and landsliding

To understand whether the spatial distribution of recent landslides is related to CSZ earthquakes or regional rainfall patterns, we identify all landslides with estimated ages less than 1000 ybp. For each slide, we extract values for modeled PGA from simulations of an M9 CSZ event (27, 28), 3D distance to the potential CSZ megathrust rupture plane (27, 28), and mean annual precipitation [MAP; (30)] at the centroid point of each slide (Fig. 4). We then plot landslide frequency per unit area against each of these metrics and test the statistical significance of each correlation by examining linear regressions and their associated confidence intervals. This analysis offers insight into the relative importance of earthquakes versus rainfall in controlling spatial patterns in the frequency of deep-seated landslides.

Fig. 4 Plots showing the relationship between landslide frequency during the past 1000 years and PGA, distance to CSZ rupture surface, and mean annual rainfall.

Area-normalized landslide frequency is plotted against (A) modeled peak ground acceleration (27, 28), (B) 3D distance to the modeled CSZ megathrust rupture surface (27, 28), and (C) MAP. Best-fit linear regression is shown as red line with shaded 95% functional confidence interval. Values for each of the three metrics considered were extracted at the centroid point of each landslide less than 1000 years old.

Although predicted PGA for a Mw 9.0 earthquake on the CSZ megathrust is substantial (ranging from 0.6g in the west to 0.2g in the east within the study area), our analysis shows no statistically significant linear correlation (Fig. 4A). A linear regression fit to the PGA-landslide frequency points yields the following equation with 95% confidence intervals: FLS = 1.07 × 10−3(±1.169 × 10−3) × PGA + 5.79 × 10−5(±12.2 × 10−5), where FLS is the yearly landslide frequency per square kilometer and PGA is unitless peak ground acceleration expressed as the proportion of g (R2 = 0.4). To account for the lack of any historic measurements of PGA from large-magnitude earthquakes along the CSZ, we also test these enigmatic results by examining a simpler and more objective measure: 3D distance from each landslide to the CSZ plate interface. We observe a similarly counterintuitive pattern of landslide frequency weakly increasing with distance from the fault before dropping off at distances greater than 30 km, once again resulting in a lack of statistically significant linear trend: FLS = 9.64 × 10−6(±17.3 × 10−6) × D + 1.48 × 10−5(±47.83 × 10−5), where D = 3D distance from landslide center to the rupture plane, in kilometers (R2 = 0.2; Fig. 4B). Using 3D distance from a potential rupture surface accounts for the relatively shallow dip angle of the CSZ megathrust (27, 28), although a similar overall pattern is observed when 1D map distance to the offshore CSZ fault trace is considered. When mean annual rainfall is plotted against landslide frequency, we do find a statistically significant positive correlation in the following form: FLS = 12.3 × 10−5(±6.15 × 10−5) × MAP + 5.79 × 10−5(±12.2 × 10−5), where MAP is mean annual precipitation, in meters (R2 = 0.7; Fig. 4C). The topography of OCR acts as an orographic barrier to offshore storms, creating a rain shadow that results in MAP that ranges from less than 1.5 m/year to over 2.5 m/year within the study area (30, 31). In wetter portions of the study area that receive ~2.5 m/year, landslide frequency is double that in the drier portions that receive less than 1.5 m/year (Fig. 4C). Collectively, these analyses give no indication that CSZ megathrust earthquakes control the location of landslides triggered during the past 1000 years but do show a strong correlation between rainfall and landslide frequency during the same period.

Estimating landslide frequency through time in the central OCR

The predominant pattern in landslide frequency within the study area is an apparent nonlinear decrease in frequency with age (Fig. 5A). This same pattern was observed in previous landslide chronologies (18, 19) and likely reflects both a real preservation bias of young slides and a bias effect associated with using an exponential function to estimate age (Supplementary Materials). The preservation bias should be expected as ongoing hillslope and fluvial processes tend to reshape relict landslide deposits, and as old slides are reactivated. Equilibrium hillslope adjustment time scales calibrated to the OCR suggest that the topographic signature of the deepest landslides in the study area may persist for tens of thousands of years after failure (32), while the legacy of smaller, shallower landslides is more short lived. This is supported by the distribution of landslide deposit area for young versus old landslides: Deposits estimated to be tens of thousands of years old have nearly double the median and mean area compared with landslides in the past 1000 years (table S2).

Fig. 5 Distribution of landslide ages and estimated landslide frequency.

(A) Histogram of landslide ages for the past 50,000 years across the study area shows a nonlinear decrease with time, which is likely the result of the preservation bias of younger-landslide deposits. (B) Mean yearly landslide frequency estimated for the past 1000 years. Frequency estimated from the best-fit roughness-age regression is shown as solid black line. A bootstrap analysis of 104 possible roughness-age regressions fit to our data is shown in blue, where the solid line is the mean and the dashed lines are the lower 5th and upper 95th percentiles of the bootstrap analysis. Timing of presumed full-rupture earthquakes shown as red solid lines with ±σ uncertainty (5).

Together, these processes lead to a preservation bias toward younger and larger landslides, similar to the well-documented apparent decline of sedimentation and erosion rates with geologic age that has been attributed to biased representation of the geologic record (33). As a result, we cannot directly compare the true frequency of landsliding from tens of thousands of years ago to present day, but we can statistically identify peaks in the landslide age-frequency data that differ significantly from the background trend and may represent periods of widespread landslide triggering due to extreme storms or large earthquakes. To test whether there is a signal from past CSZ earthquakes in the distribution of landslide ages, we focus on the age-frequency data for the past 1000 years, where the onshore and offshore geologic records suggest that three full-length (~M8.7 to 9.1) CSZ ruptures have occurred, including the most recent 1700 CE earthquake as well as earthquakes at 552 ± 83 ybp and 868 ± 58 ybp (5). After selecting a subset of 2676 landslides with estimated ages <1000 years, we bin landslide frequency into twenty 50-year intervals and then estimate slide frequency per year for each bin (Fig. 5B).

Landslide frequency over the past 1000 years reveals a nonlinear increase in landslide frequency toward the present (to the left on Fig. 5), consistent with the overall trend observed for the past ~50,000 years. Superimposed on this trend are two sharp peaks at 125 and 225 ybp, and a broad peak centered at 375 ybp (Fig. 5B). To test whether these peaks represent noise in the dataset or a real increase in landslide frequency at these times, we compare this best-fit frequency plot with a suite of frequency plots generated from a bootstrap analysis. This bootstrap analysis allows the roughness and age values for landslides of known age to vary within their individual uncertainties and then fits a unique regression for each of 104 bootstrap runs. Given the limited number of landslides with age constraints used to fit our roughness-age function, this analysis helps to address error in our landslide frequency plots caused by uncertainty in the roughness-age regression. When we plot the bootstrapped mean frequency as well as the lower 5th and upper 95th percentiles, a trend of decreasing frequency with age is observed, similar to the frequency plot generated by the best-fit regression; however, no peaks appear in the bootstrap mean distribution. When individual frequency plots from the bootstrap analysis are examined, peaks of similar amplitude as the peaks observed in the best-fit frequency plot appear at random times in the 1000-year history. Collectively, this analysis suggests that the sharp peaks observed in the best-fit frequency plot may not be statistically significant, and thus, deviations from a constant rate of landsliding during the past 1000 years should be viewed with caution. If the peaks in the observed landslide frequency plot are noise, as the bootstrap analysis suggests, then there is no justification for inferring widespread landslide-triggering events at those times.

Quantifying coseismic landslide triggering rates from CSZ earthquakes

Although we do not detect clustered landslide activity in our analysis, the overall shape of the curve is useful in determining the relative importance of earthquakes as a triggering mechanism in the OCR. Because of the large uncertainty in the roughness-age technique, pulses of coseismic landslides may not show up as discrete peaks in the resulting landsliding frequency data but could form broad undulations that control the overall shape of the dataset over the past 1000 years. Moreover, there has not been a large earthquake in Cascadia for 319 years, so we may expect that the frequency plot from a landslide inventory triggered primarily by earthquakes would peak around the time of the 1700 CE earthquake and decrease toward present day. To systematically explore these scenarios and to place bounds on the number of possible earthquake-triggered landslides in the past 1000 years, we generated simulated landslide chronologies and compared those with the observed landslide dataset. By superimposing pulses of coseismic landslides for each of the three large CSZ earthquakes in the past 1000 years (5) on a specified background rate of landsliding, we compared modeled frequency plots for three hypothetical landslide histories: (i) All landslides are seismically triggered; (ii) all landslides are triggered by background processes (we assume rainfall as the predominant background trigger); or (iii) an equal mix of seismic- and rainfall-triggered landslides. These simulated landslide chronologies allow us to discern the expected form of landslide frequency plots for the past 1000 years for different triggering scenarios. By accounting for uncertainty in the roughness-age dating method that arises from natural variability in roughness values for landslides of the same age, we present a statistically justified way of interpreting our observed frequency data.

The simulated landslide history composed entirely of “background” landslides, where no coseismic landslide pulses are added, most closely follows the observed frequency plot (Fig. 6A). This is expected as the background rate and preservation bias in the simulated histories are tuned to match the total number of slides in the observed frequency data. For this “purely background” simulated history, the observed data fall almost entirely within the 5th and 95th percentile bounds of the model runs. As we increase the number of slides that are triggered in each earthquake pulse in the simulated landslide histories, the resulting frequency plots deviate further from the observed data. In model runs where 50% of the total landslides are coseismic, corresponding to 600 landslides triggered per earthquake, the observed data fall outside the modeled bounds for two distinct intervals in the 1000-year history: from ~700 to 250 ybp and from ~150 ybp to present (Fig. 6B). The purely coseismic simulated history, for which 1300 landslides are triggered per earthquake, reveals the strongest deviation from the observed frequency data, with minimal overlap between simulated and observed landslide frequency (Fig. 6C). In the 100% coseismic simulation, the discrepancy between observed and simulated frequencies is caused by a shift in the overall shape of the curve, not by the formation of distinct frequency peaks at the times of earthquakes. The diffuse nature of these peaks causes them to overlap such that they are indiscernible. This effect is due to both natural variability in initial roughness conditions immediately after slope failure and variability in the quality of lidar data used to measure roughness, which together result in a distribution of roughness values for landslides of the same age. When we shorten the bandwidth of this roughness distribution by decreasing the SD of roughness values used in the model runs, frequency peaks do appear (fig. S2). However, this requires unrealistic SD values, which are substantially lower than the measured SD and much lower than the more conservative SD that we adopt for all model runs (see Materials and Methods and the Supplementary Materials).

Fig. 6 Comparison of observed versus simulated landslide frequency histories for the past 1000 years.

(A) Observed frequency calculated from age estimates for landslides mapped in the study area is plotted as a solid black line. Simulated landslide frequency for 104 modeled scenarios where no landslides are triggered during past full-rupture CSZ earthquakes is shown in red, where the solid line is the mean and the dashed lines are the lower 5th and upper 95th percentiles of all model runs. In the modeled scenarios, landslides are triggered at a steady rate. (B) Observed frequency shown as solid black line and simulated landslide frequency shown in red for modeled scenarios where 50% of the total modeled landslides are triggered by CSZ earthquakes, equivalent to 600 landslides per earthquake for three earthquakes during the past 1000 years. The other 50% of landslides are triggered at a steady rate. Individual coseismic pulses associated with each earthquake are indiscernible due to roughness variability for landslides of the same age. (C) Observed frequency shown as solid black line and simulated landslide frequency shown in red for modeled scenarios where 100% of the total modeled landslides are triggered by CSZ earthquakes, equivalent to 1300 landslides per earthquake.

DISCUSSION

In this study, we demonstrate that roughness dating of landslide deposits can be applied to large inventories of deep-seated bedrock landslides. In conjunction with previous studies of landslides in glacial sediments, this work suggests that roughness dating is a broadly applicable tool suitable for use in diverse landscapes underlain by either bedrock or Quaternary sedimentary deposits (18, 19). In the Tyee Formation of the OCR, we show roughness dating to be a powerful tool for rapidly estimating the ages of nearly 10,000 landslides across 15,000 km2, albeit with less accuracy than traditional landslide geochronometers such as radiocarbon dating. While this uncertainty may be too high to identify peaks in landslide frequency associated with individual earthquakes that occur every few hundred years in this region, the technique does offer an unprecedented means of parsing the relative importance of coseismic landslide triggering versus other triggering mechanisms.

Because there has not been a large earthquake in the region in over 300 years, the shape of the observed versus modeled landslide frequency plots differs substantially, with frequency dropping off rapidly toward the present in simulated landslide histories that incorporate a larger relative proportion of coseismic landslides. The observed landslide frequency over the past 1000 years shows the opposite: a slight increase in landslide frequency toward present day. Some of this increase is likely due to preservation bias, although we cannot discount the possibility of a real landslide frequency increase, potentially driven by widespread logging beginning in the late 19th century. The observed frequency is best reproduced with model simulations that include no pulses of coseismic landslides during the times of great CSZ earthquakes. In simulations where coseismic landslides account for half of the total landslides, the modeled frequency through time deviates substantially enough from the observed frequency that these simulations likely represent a conservative bound on the maximum proportion of coseismic triggering: equivalent to no more than 600 landslides per earthquake event. We cannot dismiss the possibility that CSZ earthquakes prime hillslopes for precipitation-induced landslides by damaging the uppermost few meters of regolith and bedrock without inducing catastrophic failure (17, 34). However, if this effect was substantial but short-lived, as is the case for earthquakes on shallow thrust faults (34), then we would still expect to see a peak in landslide frequency sometime after the last CSZ earthquake, but before present day.

The apparent absence of coseismic landsliding in the study area may be unexpected, given the proximity of the region to the CSZ megathrust and the intense ground shaking predicted by numerical earthquake simulation models (27, 28). However, these results are consistent with the conclusions drawn from recent coseismic landslide inventories mapped after subduction zone earthquakes in Maule, Chile (9) and Tohoku, Japan (10). These studies found that the total number of coseismic landslides fell one to two orders of magnitude short of expected values based on earthquake magnitude–landslide frequency scaling relationships developed for crustal faults (7). Moreover, coherent, deep-seated landslides made up the smallest number of landslides per failure style in the mapped inventories. Of the 1226 landslides that Serey et al. (9) mapped after the Maule earthquake, only 8 were identified as coherent, deep-seated landslides. Improved predictive models that account for depth to the rupture surface still overpredict landslide volume for recent subduction zone earthquakes (8), although the relatively shallow dip angle of the CSZ megathrust positions a greater area of landslide-prone terrain directly above the potential rupture interface compared with other subduction zones. Still, we see no evidence that predicted PGA or 3D distance to the rupture interface exerts any control on spatial patterns in landslide frequency across the study area over the past 1000 years (Fig. 4).

While catastrophic deep-seated landslides are large and may be very destructive, this study does not consider the potential for smaller shallow translational landslides or deep-seated slides that are only minimally displaced during a reactivation event. Topographic evidence for each of these landslide styles may not be perceptible in 1-m lidar data, especially after hundreds of years, so we should consider the potential for these types of landslides to be triggered on a widespread scale during large earthquakes. Specifically, shallow landslides smaller than our threshold of 5 × 103 m2 may be expected to occur in far greater number, on the basis of both landslide area-frequency scaling relationships (8, 35) and recent subduction zone earthquake–triggered landslide inventories (9, 10). Shallow slides and debris flows are particularly prevalent in steep catchments in the central OCR (36) and, importantly, could represent the dominant style of coseismic landslide in the study area. Colluvial hollows, a common debris flow source, refill with hillslope material, and shallow landslide deposits are less voluminous and more easily evacuated (37). This makes identifying and mapping these landforms difficult or impossible hundreds of years after the last large earthquake. In contrast, scars and deposits associated with deep-seated landslides may last tens of thousands of years, a time span that may encompass dozens or even hundreds of earthquakes in seismically active regions.

Given the well-established importance of rainfall as a landslide trigger (3840), the lack of evidence for seismic triggering of deep-seated landslides, and the strong positive correlation between MAP and landslide frequency, we propose that precipitation likely triggers most landslides considered in this study. While we do not attempt to quantify erosion or rates of sediment delivery by rainfall-triggered landslides in the OCR, a recent study of the Nepalese Himalaya found that, if the 2015 Gorkha earthquake is representative of landslide-triggering earthquakes over 104-year time scales, then earthquakes only trigger 10% of the total landslide volume in the region. If the entire suite of possible earthquake magnitudes and recurrence intervals are accounted for, then coseismic landslide volumes may account for up to 60% of the total. This attributes 40 to 90% of the total landslide erosional budget to rainfall-induced landslides, a similar conclusion to what we draw from the OCR, despite the substantial differences in tectonic setting, lithology, and climate. Our conclusions are bolstered by recent efforts to pinpoint the timing of individual landslide deposits in the OCR using dendrochronology, often yielding ages with uncertainties of less than 1 year (13). This ongoing work has revealed multiple deep-seated landslides that date to 1890 CE, when historical records indicate severe flooding in the region (41), but has not identified any 1700 CE coseismic landslides. Together, both methods converge on a similar dearth of 1700 CE deep-seated landslides and suggest that the primary hazard posed by these types of landslides in Cascadia is driven by rainfall rather than large-magnitude subduction zone earthquakes.

MATERIALS AND METHODS

Landslide mapping and dating

To construct a surface roughness–age relationship, we mapped 9938 deep-seated landslides in the Tyee Formation and overlying Elkton Formation. Because surface roughness dating requires the assumption that landslides smooth over time because of erosion and soil diffusion, landslide deposit polygons had to be mapped carefully to avoid headscarps, incised stream gullies, and river cutbanks. Mapping was completed manually using bare-earth lidar-derived hillshade and slope images. Only the presently intact surfaces of landslide deposits were mapped to capture areas dominated by diffusive smoothing processes. Major roads and anthropogenic structures were removed from polygons before roughness was calculated (fig. S1).

We mapped only coherent rotational and translational deep-seated landslide deposits with areas of at least 5 × 103 m2 and with at least one dimension over 100 m. Although we can discern smaller landslides in the lidar imagery, we use a 20-m wavelet to measure curvature for a roughness metric, so slides much smaller than this seemed inappropriate for consideration in this study. We define coherent landslides as slope failures where the mobilized material has not completely disintegrated to the point of allowing flow or avalanche-style runout. Coherent slides more consistently preserve roughness elements like intact blocks, minor scarps, hummocks, and compressional ridges needed for roughness dating. Wherever remobilized portions of older landslide deposits satisfying these criteria were observed, they were mapped as individual slides. To satisfy the assumption inherent to surface roughness dating that all landslides fall within a range of expected initial roughness values directly after they occur, certain failure modes were excluded from the analysis. Landslides with channelized travel paths, debris flows, and earthflows were not considered, as the mechanics of these landslide types do not universally produce classic deep-seated landslide morphologies such as hummocky topography, scarps, and displaced blocks characterized by high initial roughness, as is typical of coherent landslides.

A suite of landslides with independently constrained ages is necessary for roughness-age dating, as these landslides provide the data to regress against roughness to define the roughness-age function used to estimate the ages of other landslides in the study area. We used radiocarbon dating and historical repeat aerial imagery to ascertain the timing of failure for six landslides in the mapped chronology (table S1). In addition, we used dates from previous work, including recent dendrochronology dating, for 8 other landslides (21, 22, 25, 26, 4143) yielding a total dataset of 14 landslides of known age (table S1). Because we date multiple landslides younger than 1950 CE, we report dated landslide ages in this study as years before 2019, as opposed to the radiocarbon convention of years before 1950. In this manuscript, we use the phrase “years before present” and its abbreviation “ybp” to mean years before 2019.

Surface roughness analysis and landslide age estimation

For roughness dating to be applicable to a landslide-prone region, the initial roughness of landslide deposits is assumed to be similar, falling within a range of relatively high expected values. An additional assumption is that the material properties of the landslide and the climate of the study area are both similar enough that the rate of erosion and soil diffusion does not vary markedly among landslides within the inventory. Moreover, while this technique has been shown to be applicable to many landslides across large study areas (19), it has not yet been implemented on bedrock landslides. To test the utility of the approach on bedrock landslides, we focus our efforts on a swath of the central OCR, underlain entirely by Eocene marine sedimentary rock of the Tyee Formation and overlying Elkton Formation.

We measured the roughness of each mapped landslide using the 2D CWT with a Mexican Hat wavelet (fig. S1), which is equivalent to calculating topographic curvature (the Laplacian of elevation) at the specified length scale [(44) and the Supplementary Materials]. We first calculated the CWT at spatial scales of 4, 15, 20, and 45 m to find the scale that offered the best fit to the landslides of known age (table S3). For each CWT analysis, roughness values were computed for each cell in the 0.9-m bare-earth digital elevation model, and mean values of the magnitude for each landslide deposit of known age were extracted from mapped polygons and used to plot age against roughness. Two regression methods, robust (least absolute error or L1 norm) and least squares (L2 norm), were used to fit single exponential functions to these roughness-age datasets. We tested the performance of both methods in fitting the entire age range of dated landslides as well as a subset of young landslides (<1000 ybp). It was important to select a fit that did not result in a systematic overprediction or underprediction in age estimations for the youngest landslides, as we used this fit to examine frequency trends in the past 1000 years. When the entire range of ages in the dataset was considered, the SE of the least squares regression was slightly lower (~1.2 × 103 years) compared with the robust regression (~1.6 × 103 years). The least squares regression optimized the fit for the oldest (>40,000 ybp) landslides but systematically overestimated the ages of the younger-landslide data points. When just the 10 landslides less than 1000 ybp were considered, the least squares regression resulted in nearly double the SE as the robust regression: ~400 years versus ~200 years, and examination of residual plots showed less bias in the robust regression. Because the robust regression method is less sensitive to outliers than the least squares regression, it most effectively fits a function across the age range of the dataset, especially for the younger (<1000 ybp) landslides. This robust regression, which we deem the best fit to our observed data, yields the following function that we use to estimate landslide age across the study area: t = 1.428 × 106 × e−412.8R, where t is the estimated landslide age, in years before 2019, and R is the average magnitude of the wavelet coefficient at a 20-m length scale within a mapped landslide polygon, m−1 (Fig. 3).

Analysis of PGA, distance to CSZ megathrust, and MAP

To identify any correlation between the locations of recent landslides and patterns in modeled PGA, 3D distance to the CSZ, and MAP, we plot histograms of these data extracted from the center of each mapped landslide with estimated ages less than 1000 ybp. For modeled PGA, we adopt a mean of 30 simulated M9 CSZ earthquakes designed to capture a wide range of potential ground-shaking outcomes from variations in hypocenter location, down-dip rupture limit, slip distribution, and subevent locations (27, 28). This mean is downsampled to a 1-km grid using bilinear interpolation to smooth unrealistically large jumps between cells in the raw data. We also consider a simpler measure of CSZ seismic hazard: the minimum 3D distance from each landslide point to the CSZ megathrust rupture surface. We take the rupture surface to be the interface between the Juan de Fuca Plate and the overlying North American plate from depths of 5 to 30 km (27, 28). For rainfall, we use a 30-year record of MAP (30). Although MAP does not capture the hourly intensity of individual storms, the deep-seated slides that we consider in this study are more likely to be sensitive to total precipitation volume on time scales of days to months [e.g., (39)]. Because of the highly seasonal nature of precipitation in the OCR along with the spatial variability in winter storms, a 30-year record of MAP is a reasonable metric to capture the overall patterns in rainfall that are likely to control precipitation-induced landsliding.

For each of the three metrics, we plot histograms by first dividing the middle 95% of the total range of values across the study area into 10 equal intervals. Excluding the upper and lower 2.5% helps to avoid extreme values that affect little of the study area. The number of landslides for each value interval is then normalized by total area that interval represents within the study area and, lastly, converted into a yearly frequency. The resultant y axis shows the average frequency of landslides per kilometer for each equal interval of PGA, distance to CSZ, and MAP. To identify the correlation between each of the three metrics and the locations of recent landslides, we fit a line to the data in each plot and calculate the R2 value to test the goodness of fit (Fig. 4).

Simulated coseismic landslide chronologies

If large CSZ earthquakes do periodically trigger widespread landsliding in the study area, then we would expect the resultant peaks in the histogram to be diffuse rather than sharp, because of the variability of measured roughness values within each pulse of coseismic landslides. To account for this variability and to aid in a more robust interpretation of the observed landslide frequency data, we generated simulated landslide chronologies. These simulated chronologies offer a means to estimate the proportion of landslides in the study area triggered by large CSZ earthquakes. Specifically, they allow us to better determine what a coseismic signal should look like. By comparing the observed and synthetic datasets, we can also deduce whether the apparent steady increase in landslide frequency that we observe could be due to uncertainty in the roughness dating technique.

We assume that when a pulse of landslides occurs, such as from a large earthquake, landslides fail instantaneously, and their deposits have normally distributed roughness values. Over time, both the mean and the SD of the roughness values that compose this pulse decrease. The SD decreases because the roughest landslides have more regions of high curvature, and therefore smooth considerably faster than the landslides with lower initial roughness, effectively shrinking the distribution. This has been demonstrated numerically (19) and is supported by both geomorphic transport laws (45) and the exponential form of the roughness-age regression itself. To illustrate this effect in the bedrock setting of the central OCR, we model the topographic evolution of the roughest and smoothest modern landslides in our dataset that have failed in the past 20 years. This exercise illustrates how the initial difference in roughness values between the two slides rapidly shrinks after only a few hundred years (figs. S3 to S5).

To calculate the SD of roughness values for a specific age, a dataset of landslides of known age would ideally include many landslides of the same age. From our landslides with known ages, we calculate the SD for four landslides with ages that fall within 100 years of the 1700 CE earthquake, yielding a value of 8.4 × 10−4 m−1. However, to address the uncertainty in calculating the SD from a small sample, we chose to double the measured SD and instead use an SD of 1.7 × 10−3 for all modeled landslide distributions. Rather than attempt to scale down SD with age, we take a conservative approach to this uncertainty and use the same value for the entire 1000-year simulation. We acknowledge that the modeled roughness SD has a substantial effect on the resulting landslide frequency plots. Very high SD values yield estimated ages that are skewed toward present day, producing a steep nonlinear increase in apparent frequency toward present day. However, this does not occur until SD values exceed 6.8 × 10−3, a value eight times greater than that measured in the study area for landslides ~1700 CE (fig. S2).

Simulated landslide chronologies are run for the past 5000 years and consist of two components: a background rate of sliding and pulses of coseismic landslides. Although we only examine landslide frequency over the past 1000 years, we extend the simulated chronology to 5000 ybp to eliminate any effects from old landslides with erroneously young estimated ages. To build these chronologies, we first generate 10 coseismic landslide pulses at the times of past widespread turbidite events T1 (1700 CE) to T10. These widespread turbidite triggering events, identified by Goldfinger et al. (5), are thought to correspond to past, full-length CSZ ruptures with estimated magnitudes of 8.7 to 9.1. The number of landslides triggered during each earthquake can be adjusted, allowing for any proportion of coseismic landsliding to be examined. Next, we prescribe an adjusted background rate of landsliding, in slides/year. This represents a spatial and temporal average across the entire study area of the rate of nonseismic landslide occurrence per year. At each time step in the chronology, a randomly generated, normally distributed set of landslide roughness values corresponding to the time step age is added. The background rate is constant for each model run but is scaled depending on the number of coseismic landslides triggered in earthquake pulses so that, after adjusting for preservation bias, the total number of simulated landslides estimated to have ages of less than 1000 ybp is equal to the total number observed: ~2700 landslides.

Because of the error associated with converting a roughness value to an age, coseismic landslide pulses have tails that extend past 1000 ybp, so the prescribed total number of coseismic landslides triggered by the three large earthquakes in the past 1000 years is greater than the total number estimated to be less than 1000 ybp. To directly compare simulated landslide chronologies to our observed chronology, we normalize for the same number of landslides with estimated ages less than 1000 years, even though this resulted in a variable number of prescribed landslides for each simulation. For example, for the simulated 50% coseismic landsliding scenario, 600 coseismic landslides are prescribed per earthquake for a total of 1800 coseismic landslides in the past 1000 years, while only 1317 are estimated to have ages less than 1000 ybp (Fig. 6B). To account for the preservation bias in the observed frequency, the total number of landslides for each simulated chronology is scaled as a function of time. We adopt a normalized exponential fit to the observed frequency data as this scaling function, which allows for a direct comparison between simulated and observed frequency plots.

We show results from simulations of 0% coseismic triggering, 50% coseismic triggering, and 100% coseismic triggering, corresponding to 0, 600, and 1300 coseismic landslides triggered per CSZ earthquake, respectively. Each of these three sets of conditions is modeled 104 times, and then the ages are calculated for each run using the same best-fit roughness-age regression used to calculate observed ages. Last, the mean frequency as well as 5th and 95th bounding percentiles for all model runs are calculated and plotted (Fig. 6).

Statistical analysis of age estimates and sources of uncertainty

The total error in our estimated landslide ages stems from natural variability in initial roughness conditions, uncertainty in roughness calculations, uncertainty in interpreting radiocarbon-based landslide ages, and the inherently subjective nature of mapping landslides by eye. We attribute most of the uncertainty in our surface roughness calculations to variations in ground point density from the lidar point clouds used to generate the digital elevation models used in this study, an effect that has been noted in previous work (46). To quantify this effect and assign error bars to points on the roughness-age curve, we compare the roughness of two lidar datasets taken in 2011 and 2015 of the same region within the OCR (Supplementary Materials). We assign the ±2σ m−1 of the observed error between repeat lidar datasets (9.6 × 10−4) as the ±2σ m−1 for each point used in our roughness-age regression and plot this value as horizontal error bars on landslides of known age (Fig. 3).

Because our age estimations are determined by the parameters of the best-fit regression, we needed to account for uncertainty in the regression itself. Uncertainty in the regression would primarily arise from uncertainty in the roughness-age data points used in the fitting process. To address the effect of error in these, we use a bootstrap analysis. For 104 runs, the roughness and age for each of the 14 landslides of known age were randomly selected at a frequency dependent on the respective uncertainty distribution for each landslide, and a best-fit exponential regression was fit. Each bootstrapped regression is then used to estimate landslide ages and generate landslide frequency versus time plots for both the measured roughness values in the study area and the roughness values generated in our simulated coseismic landslide datasets (Fig. 5). This allows us to interpret trends in our observed frequency in a more statistically robust manner.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/38/eaba6790/DC1

https://creativecommons.org/licenses/by-nc/4.0/

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank S. Harbert and L. Wetherell for assistance in the field, V. Bright and K. Lowery for assistance mapping landslides, and P. Schoettle-Greene for thoughtful insight during numerous discussions. This manuscript benefitted tremendously from reviews by N. Hovius, C. Cerovski-Darriau, J. Perkins, A. Baltay-Sundstrom, and an anonymous reviewer. Funding: We acknowledge the National Science Foundation (grant EAR-1331412) and the Geological Society of America (graduate student grant to S.R.L.) for support of this research. Author contributions: S.R.L., A.R.D., A.M.B., A.G., W.S., and J.J.R. conducted fieldwork to date landslides used in the roughness-age regression. B.A.M. developed methods for data analysis to explore the sensitivity of frequency data to regression parameters. J.W. provided expertise regarding coseismic landslide initiation and helped refine related passages in the manuscript. S.R.L. manually mapped landslides, developed the simulated landslide history models, and wrote the bulk of the manuscript. All authors were instrumental in crafting the study and revising the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. For questions regarding access to the database of mapped deep-seated landslides considered in this manuscript, please email the corresponding author at seanlah{at}gmail.com.
View Abstract

Stay Connected to Science Advances

Navigate This Article