Abstract

Characterizations of tsunami hazards along the Cascadia subduction zone hinge on uncertainties in megathrust rupture models used for simulating tsunami inundation. To explore these uncertainties, we constructed 15 megathrust earthquake scenarios using rupture models that supply the initial conditions for tsunami simulations at Bandon, Oregon. Tsunami inundation varies with the amount and distribution of fault slip assigned to rupture models, including models where slip is partitioned to a splay fault in the accretionary wedge and models that vary the updip limit of slip on a buried fault. Constraints on fault slip come from onshore and offshore paleoseismological evidence. We rank each rupture model using a logic tree that evaluates a model’s consistency with geological and geophysical data. The scenarios provide inputs to a hydrodynamic model, SELFE, used to simulate tsunami generation, propagation, and inundation on unstructured grids with <5–15 m resolution in coastal areas. Tsunami simulations delineate the likelihood that Cascadia tsunamis will exceed mapped inundation lines. Maximum wave elevations at the shoreline varied from ∼4 m to 25 m for earthquakes with 9–44 m slip and Mw 8.7–9.2. Simulated tsunami inundation agrees with sparse deposits left by the A.D. 1700 and older tsunamis. Tsunami simulations for large (22–30 m slip) and medium (14–19 m slip) splay fault scenarios encompass 80%–95% of all inundation scenarios and provide reasonable guidelines for land-use planning and coastal development. The maximum tsunami inundation simulated for the greatest splay fault scenario (36–44 m slip) can help to guide development of local tsunami evacuation zones.

INTRODUCTION

There were no direct observations of the seafloor deformation produced by the moment magnitude (Mw) ∼9 great Cascadia earthquake of A.D. 1700 that generated a Pacific-wide tsunami (Satake et al., 1996, 2003; Atwater et al., 2005). Although the impacts of the A.D. 1700 earthquake and tsunami and their predecessors left an enduring impression in the oral histories of native North Americans (Ludwin et al., 2005), great earthquakes and tsunamis escape mention in European written documents. Because of the lack of observations, simulations of Cascadia tsunami inundation require formulating hypothetical earthquake scenarios for the Cascadia megathrust. This paper describes earthquake scenarios developed using: (1) knowledge of the structure of the Cascadia megathrust (Goldfinger, 1994; McCrory et al., 2004); (2) coastal and offshore paleoseismic records, some extending back 10,000 yr (e.g., Goldfinger et al., 2012); (3) theoretical understanding of how megathrust ruptures deform the seafloor (Wang and Hu, 2006; Wang and He, 2008); and (4) observations after historical great earthquakes (e.g., Chlieh et al., 2007; Hsu et al., 2006; Subarya et al., 2006; Moreno et al., 2009; Ito et al., 2011).

Here, we characterize 15 earthquake rupture scenarios for the entire (∼1000-km-long) margin by considering a range of geological and geophysical information, including onshore paleoseismic evidence and deposits of turbidity currents triggered offshore by great earthquakes (Goldfinger et al., 2012). These scenarios build on the results of a companion study of Cascadia tsunami inundation hazard at Cannon Beach, Oregon (Priest et al., 2009, 2010), and they provide the foundation for tsunami inundation hazard maps for the entire ∼580-km-long coast of Oregon. To evaluate our approach and explore the range of tsunami hazards, we focused on Bandon, Oregon, where the continental shelf becomes narrower and the deformation front approaches the coast of southern Oregon and northern California (Fig. 1A). Southern Oregon (Fig. 1B) offers a rich geologic history of repeated coseismic subsidence (Kelsey et al., 1998, 2002; Witter et al., 2003) and tsunami inundation (Kelsey et al., 2005; Nelson et al., 2006) caused by past great Cascadia earthquakes.

This paper describes the approach and methods used to design geologically defensible Cascadia earthquake scenarios and tsunami simulations that span the range of variability of the tsunami-generation process. The findings highlighted here are condensed from a comprehensive technical report released by the Oregon Department of Geology and Mineral Industries (Witter et al., 2011) that applied our modeling approach by assessing tsunami hazards along the southern Oregon coast (Fig. 1B). The technical report includes supplementary modeling data, geographic information system layers, tsunami animations, and a tsunami evacuation map for the city of Bandon, Oregon.

BACKGROUND

Paleoseismic Evidence of Cascadia Earthquakes and Tsunamis

Onshore and offshore paleoseismic evidence in southern Oregon provides a long (10,000 yr) and complete history of great megathrust earthquakes that we used to construct hypothetical earthquake rupture models for the Cascadia megathrust (Fig. 1).

Evidence of Coseismic Subsidence in Estuaries

Stacks of buried tidal marsh and upland soils beneath wetlands surrounding the Coquille and Sixes River estuaries (Fig. 1B) chronicle 12 episodes in the past 6700 yr when great Cascadia earthquakes dropped the area by as much as 1.2–3 m (Kelsey et al., 1998, 2002; Witter et al., 2003). These buried soils offer a 6700-yr-long, on-land record of repeated coseismic subsidence during great Cascadia earthquakes and support speculation that some past Cascadia ruptures were limited to southern fault segments (Nelson et al., 2006). Sandy deposits overlying the buried soils imply inundation by Cascadia tsunamis generated by great earthquakes. At the Coquille River, sandy deposits left by the A.D. 1700 tsunami extend 2 km inland (Witter et al., 2003). At the Sixes River, 30 km south of Bandon, the A.D. 1700 tsunami blanketed the lower kilometer of the river’s floodplain with a sheet of sand as thick as 35 cm (Kelsey et al., 1998). At both sites, older tsunami deposits reach farther into the estuaries. Ten kilometers inland from the mouth of the Coquille River, the thickest tsunami deposit exceeds 60 cm. The spatial distributions of these deposits provide crude constraints on the minimum extent of tsunami inundation, but using their extents to validate tsunami simulations is complicated by probable changes in the coastal landscape over the last 6000–7000 yr.

Evidence of Tsunamis and Strong Shaking in Bradley Lake

Located above the Cascadia megathrust within 0.5 km of the Pacific Ocean, Bradley Lake (Fig. 1B) became a natural tsunami recorder ∼7500 yr ago when landward-advancing sand dunes dammed coastal streams. Distinctive beds of debris and sand, which interrupt faintly laminated lake mud, probably reflect seismic shaking and tsunami inundation in the lake caused by 17 great Cascadia earthquakes (Kelsey et al., 2005). In 13 cases, the sand beds are best explained by sustained ocean surges that transported beach and dune sand into the lake. Kelsey et al. (2005) argued that over the period when the lake was an optimal tsunami recorder, its elevation (∼6 m above sea level) was too high to be breached by Pacific Ocean tsunamis initiated by distant earthquakes. Other processes, like short-period waves (tens of seconds) driven by North Pacific storms, offered less likely explanations of the disturbances. For justification, Kelsey et al. (2005) pointed to short-lived peaks in lake salinity, evident from Holocene marine diatoms present in disturbed sediments that required marine water to flow into the lake for tens of minutes to hours—conditions consistent with inundation by Cascadia tsunamis, including the most recent waves in A.D. 1700.

Evidence of Earthquake-Triggered Turbidity Currents Offshore

Goldfinger et al. (2012) examined sediment cores from multiple submarine channels, slope basins, and abyssal fans along the Cascadia subduction zone to evaluate Holocene turbidites as evidence for great earthquakes on the Cascadia megathrust. They ruled out nonseismic processes by arguing that simultaneously triggered turbidity currents at sites separated by 50–150 km required a common event. In their view, the common event was a great earthquake, based on (1) agreement between onshore and offshore radiocarbon age chronologies, (2) consistent numbers of turbidites above and below channel confluences, and (3) stratigraphic correlation based on geophysical core logs. Goldfinger et al. (2012) correlated ∼19 sandy, relatively massive turbidites along 800 km of the Cascadia margin, from Barkley Canyon in the north to Eel Canyon in the south (Fig. 1A). Turbidite age estimates broadly overlap age ranges for Cascadia earthquakes estimated from shorter, coastal paleoseismic records (e.g., Atwater et al., 2004; Kelsey et al., 2002; Witter et al., 2003) (Fig. 2). The lateral continuity of turbidites along the Cascadia margin implied that 19 full-length or nearly full-length ruptures of the megathrust triggered widespread turbidity currents along most of the margin (Goldfinger et al., 2012).

In addition, sediment cores collected offshore southern Oregon contain ∼22 thinner, mud turbidites that appear to correlate between multiple sites along the southern margin, including Hydrate Ridge and Rogue Canyon (Fig. 1A). The along-strike continuity of the mud turbidites is supported by 29 new cores and high-resolution compressed high-intensity radar pulse (CHIRP) seismic-reflection profiles that image stratigraphy between core sites (Goldfinger et al., 2013b). Goldfinger et al. (2012) showed that each of the 41 turbidites in the 10,000 yr Cascadia record meets a variety of criteria that favor synchronous triggering by plate-boundary earthquakes. Mud turbidites deposited along the central and southern parts of the margin suggest that ∼22 additional earthquakes ruptured shorter segments of the southern megathrust.

APPROACH AND METHODS

Cascadia Earthquake Scenarios

The earthquake scenarios described herein depict full-length rupture of the Cascadia megathrust and the corresponding surface deformation used for tsunami simulations. The earthquake rupture models include slip partitioned to a hypothetical splay fault in the accretionary wedge and models that vary the updip limit of slip on the megathrust. Evidence supporting full-length rupture comes from offshore sandy turbidites that have been interpreted to record 19 Cascadia earthquakes (Goldfinger et al., 2012). We chose to develop full-length earthquake scenarios because they provide surface deformation models needed for regional tsunami hazard assessment. Tsunamis generated by smaller but more frequent ruptures on the northern (Atwater and Griggs, 2012) or southern parts of the Cascadia margin (Goldfinger et al., 2012) present additional hazards that we did not assess.

Coseismic Slip Estimates Determined from Turbidite Records and Plate Convergence Rates

Each Cascadia earthquake scenario incorporates the assumption that a slip deficit accrues at the plate convergence rate (i.e., coupling ratio = 1.0), and that this deficit is recovered by peak coseismic slip during long (>800 km) ruptures of the subduction zone megathrust. We estimate maximum coseismic slip for a range of scenarios by multiplying time intervals between earthquakes (inferred from the 10,000-yr-long Cascadia earthquake chronology of Goldfinger et al., 2012) by the plate convergence rate, which varies with latitude (Wang et al., 2003). Interevent time intervals that separate the 19 sandy turbidites range from as little as ∼110 yr to as long as ∼1150 yr (Table 1; Goldfinger et al., 2012). At the latitude of southern Oregon, calculations of maximum fault slip that use a convergence rate of 34 mm/yr imply slip of <5 m to over 40 m for the range of interevent times considered. The method to account for postulated slip during smaller southern Cascadia ruptures implied by the 22 muddy turbidites examined by Goldfinger et al. (2012) is explained later in this section.

Four time intervals were selected to represent four general earthquake size classes (Fig. 3; Table 2): small (SM, 300 yr), medium (M, 525 yr), large (L, 800 yr), and extra-large (XL, 1200 yr). A fifth extra-extra-large (XXL) scenario was used to simulate a maximum tsunami to guide evacuation planning. The time intervals represent mean values (rounded to the nearest quarter century) of binned interevent times estimated using radiocarbon ages with 2σ errors, many longer than a century, that underpin Cascadia earthquake chronologies (Atwater and Hemphill Haley, 1997; Nelson et al., 2006; Goldfinger et al., 2012). A complete uncertainty analysis of the turbidite radiocarbon chronology was described by Goldfinger et al. (2012).

The extra-large (XL) size class reflects a single event (T11) in the Holocene record of Cascadia turbidites (Table 1). T11 is the second largest turbidite by mass and has the longest postevent time interval (∼1150 yr). We rounded the interval for the extra-large scenario to 1200 yr for the purpose of modeling a hypothetical maximum earthquake with peak slip (>40 m) similar to the Mw 9.5 1960 southern Chile earthquake (Barrientos and Ward, 1990; Moreno et al., 2009). Large (L) events are defined by three turbidites with above-average masses and postevent times that ranged from ∼680 to 1000 yr (Table 1). We selected 800 yr, the approximate mean of the interevent times of turbidites with above-average mass, to estimate slip for large scenarios. Medium (M) scenarios reflect 10 Cascadia turbidites with typical masses and postevent intervals that range from ∼310 to 660 yr. We used 525 yr to define the medium earthquake scenarios because it approximates both the mean interevent time for all 19 turbidites and the mean time interval following turbidites of typical mass. Finally, slip estimates for small (S) earthquake scenarios use a 300 yr time interval, which reflects the average of postevent times ranging from ∼110 to 480 yr for the five turbidites with below-average masses (Table 1).

To account for potential slip released by hypothetical ruptures on southern segments of the megathrust postulated by Goldfinger et al. (2012) on the basis of the 22 mud turbidites, we constructed a slip budget that decreases slip from north to south along the length of the margin (Table 2; for more detailed description of the method, see Witter et al., 2011). If some Cascadia earthquakes involved ruptures of southern segments, then they might have accommodated some fraction of the total slip available from plate convergence between earthquakes. However, instead of using the extremely short interevent times for southern Cascadia turbidites to estimate this fraction of slip, estimates for southern fault segments used for rupture models and tsunami simulations adopted the minimum modeled fault slip consistent with tsunami inundation evidence at Bradley Lake (Witter et al., 2012). For ∼20 hypothetical ruptures of southern fault segments, we assigned ∼2 m of slip per event accumulated over 60–70 yr. Thus, our source model varies slip among the earthquake scenarios on segments of the megathrust so that the cumulative interevent time balances the total time span between T18 and T1, or ∼9600 yr (Table 1).

The XXL scenario is the exception to the bookkeeping described above, because for the maximum rupture model, we assumed no slip was accommodated by smaller ruptures of the southern Cascadia megathrust.

Updip Limit of Coseismic Slip

Varying the updip limit of coseismic slip in earthquake rupture models changes the amount of seafloor deformation in deep water, which controls the size of tsunami waves and how soon they reach the coast. If the shallowest portion of the megathrust is nearly flat, models that involve trench-breaking rupture predict less seafloor deformation than models that strengthen the updip edge of the megathrust (Wang and He, 2008). We propose two buried rupture models that use different updip limits of coseismic slip and assume that velocity-strengthening behavior of the fault’s shallowest segment will act as a barrier to rupture, preventing slip from reaching the seafloor (Wang and Hu, 2006). The first model, a shallow buried rupture model (e.g., Priest et al., 2010), simulates a buried rupture of the megathrust with slip tapering to zero at shallow depths (<2 km) near the deformation front (Fig. 4). The second model is a deep buried rupture model that tapers slip to zero on a deeper part of the megathrust. The updip limit of the second model is defined by the boundary between the inner (east) and outer (west) accretionary wedges offshore Washington and northern Oregon, which display distinctly different structural vergence (MacKay et al., 1992; Gutscher et al., 2001; Adam et al., 2004; Underwood, 2002). Because there are no direct observations of the behavior of the outer wedge in Cascadia subduction zone events, we treat the two models equally as possible sources for tsunami generation.

Slip Diverted to a Splay Fault

The third rupture model diverts slip to a hypothetical splay fault in the accretionary wedge that increases seafloor uplift and greatly amplifies the tsunami. Seismic lines and bathymetric maps of the accretionary complex offshore Washington and northern Oregon reveal physical evidence of a shallow splay fault system along the Cascadia forearc (Fig. 5). Seismic-reflection profiles across the continental slope (e.g., profile L-5-WO77–12; Mann and Snavely, 1984) show a fault zone separating the older accretionary complex on the east from younger accretionary wedge sediments to the west. The older accretionary wedge is deformed by seaward-vergent structures, whereas the younger wedge features widely spaced, landward-vergent folds (e.g., Adam et al., 2004). The spatial coincidence of seafloor scarps at the updip end of a structural boundary separating the older and younger accretionary wedges suggests deformation by a system of splay faults that may be activated by megathrust slip (Park et al., 2002; Wang and Hu, 2006).

We interpret this structural boundary and coincident seafloor scarps as a system of splay faults that deform Holocene sediments in the forearc and, in some areas, appear to break the seafloor (Goldfinger, 1994). The fault’s dip, ∼30° landward, comes from thrust faults in the triangle zone evident on a seismic profile offshore central Washington (lines SO 108–103; Adam et al., 2004). We generalize the splay fault system (Fig. 4) as a hypothetical fault that dips 30° landward and smoothly merges with the megathrust at depths less than ∼20 km. South of Cape Blanco, the surface trace of the hypothetical splay fault merges with the deformation front. Merging the splay fault with the updip edge of the megathrust has two consequences for the deformation models. First, seafloor uplift amplified by the splay fault gradually decreases southward as the splay fault approaches the deformation front. Second, there is no splay fault–related uplift in earthquake deformation models south of latitude 42.8°N.

Cascadia Earthquake Deformation Model

The amount and distribution of coseismic slip on the megathrust are the most important parameters controlling surface deformation produced by Cascadia subduction earthquakes (Wang, 2007). We assume these factors dominate the uncertainty in modeling coseismic surface deformation and ignore secondary effects related to material heterogeneity, inelastic behavior, and dynamic deformation. We use only the vertical component of seafloor displacement, since that is the principal control on tsunami generation for subduction zones like Cascadia that lack a steep slope at the oceanic trench (see review by Tanioka and Satake, 1996). Lacking direct observations of coseismic slip patterns produced by past great earthquakes, we use a regional slip distribution that simulates a narrow rupture along the entire length of the plate boundary (Fig. 6). Using a uniform rupture tends to exaggerate tsunami wave heights in areas where coseismic slip directly offshore is actually smaller than described by the uniform model. All models simulate surface deformation by numerically integrating the point source dislocation solution of Okada (1985) over a three-dimensionally (3-D) curved Cascadia megathrust fault in a uniform elastic half-space with a Poisson’s ratio of 0.25. Details of the modeling method, including triangulation of the fault for the integration and determination of convergence/rupture direction by correcting for a secular forearc motion, can be found in Wang et al. (2003).

Fault rupture models use a slip distribution that is symmetrical in the dip direction, following a modified version of the function of Freund and Barnett (1976) (Wang and He, 2008; Wang et al., 2013). Peak slip is located in the bottom part of the interseismic fully locked zone modeled by global positioning system (GPS) data (as described earlier; Wang et al., 2003), with slip tapering updip and downdip. For example, the downdip limit of rupture is defined by obtaining a best match to coseismic subsidence data summarized by Leonard et al. (2010) for the A.D. 1700 earthquake and more recent estimates by Hawkes et al. (2011), and to geological data that reflect the downdip limit of interplate coupling (Priest et al., 2009, 2010). For the buried rupture models, maximum slip at each location is obtained from the local plate convergence rate multiplied by the interevent time of the scenario being simulated (Table 2). In plan view, the splay fault slip model is simply the shallow buried rupture model truncated at the surface trace of the splay fault. The mesh used for the hypothetical splay fault was extremely dense to allow accurate modeling of its surface-breaching rupture. Seafloor uplift is amplified by the steeper dip of the hypothetical splay fault (Fig. 6).

Numerical Model of Tsunami Wave Propagation and Inundation

Hydrodynamic Model, SELFE

To simulate Cascadia tsunamis, we used the hydrodynamic model SELFE (Zhang and Baptista, 2008)—the Semi-implicit Eulerian-Lagrangian Finite Element code developed for modeling cross-scale ocean circulation, storm surges, and tsunamis. The algorithms used to solve the Navier-Stokes equations are computationally efficient and stable. All tsunami simulations were run assuming that prevailing tide was static (no flow) and equal to mean higher high water (MHHW) at the Port Orford tide gauge.

The computational grid was constructed by first compiling a digital elevation model (DEM) covering the project area and then retrieving from the DEM elevations at a series of points defining a triangular irregular network (TIN) utilized as the computational grid. Refer to Witter et al. (2011) for a complete description of the bathymetric and topographic data used to construct the grid. Grid spacing in the TIN differed from the detailed DEM in order to economize grid size while also accurately representing coastal features that control tsunami propagation and inundation, such as jetties, breakwaters, channels, and abrupt changes in slope. Grid spacing was adjusted so that at least five elements specified the width of jetties and breakwaters.

The final computational grid has ∼1.4 million nodes and 2.8 million triangles with finest resolution of ∼1 m (Fig. 7). Approximately 90% of all numerical elements (centers of triangles) have a resolution between 1 m and 20 m (Fig. 7B). High resolution (5–15 m) was used in the Coquille River estuary, from the shoreline to 20 m elevation near Cape Blanco and Port Orford, and around all major coastal lakes from Port Orford to Coquille. Resolution in the most densely populated areas was 5–7 m.

RESULTS

Logic Tree Evaluation of Earthquake Size and Rupture Geometry

Because information about the characteristics of past Cascadia earthquakes is incomplete, we used a logic tree to rank 15 rupture scenarios for the Cascadia subduction zone (Fig. 8). Logic trees have become useful tools for evaluating uncertainty captured by alternative model parameters in both probabilistic and deterministic seismic hazard analyses (Bommer et al., 2005). The branches of the logic tree presented here define two main rupture parameters: earthquake size and fault geometry. Each branch is weighted relative to other branches based on the degree to which the parameter is consistent with geological and geophysical data, theoretical models, and the judgment of the authors (Fig. 8; Table 3). The final weight for each scenario is the product of the weights of the two parameters, and it assigns higher rank, or likelihood, to rupture scenarios that are more consistent with existing data. We use the term likelihood to express relative and cumulative logic tree weights, which are not intended to represent the frequency or annual probability that a particular scenario will happen (Abrahamson and Bommer, 2005).

In the first node of the logic tree (Fig. 8), relative earthquake size is inferred from interevent time intervals from Cascadia turbidite paleoseismology (Goldfinger et al., 2012). Five size classes (SM, M, L, XL, and XXL) represent the relative amount of coseismic slip assigned to each scenario calculated from interevent times and plate convergence rates (Table 3). Higher slip assigned to an earthquake scenario corresponds to higher estimated moment magnitude (Table 4). We assign weights to branches of different earthquake size according to the number of events of that size inferred from the record of 19 turbidites correlated along most of the Cascadia margin, as described earlier (Fig. 3). Therefore, weights for the five earthquake size branches reflect the number of earthquakes in each size class over the last 10,000 yr and provide reasonable forecasts of the size of future great earthquakes.

From the analysis of turbidite mass and interevent times, we assign the weights for earthquake size as follows. Because five of the 19 turbidites have below-average masses, we infer they reflect small (SM) Cascadia earthquakes. So, we assign a weight of 0.26 (5 divided by 19) to the branch representing the small earthquake scenario. We assign a weight of 0.53 for medium (M) earthquakes based on the interpretation that 10 turbidites record earthquakes of moderate size. Large (L) scenarios are assigned a weight of 0.16 because three of the most massive turbidites had particularly long following times. Outsized events, including XL and XXL scenarios, were assigned equal weights that sum to 0.05 to account for a single turbidite in the Holocene record, T11, characterized by its large mass and unusually long postevent interval of ∼1150 yr along the northern margin.

Geodetic observations from leveling, tide gauge measurements, and GPS do not differentiate models that define the updip limit of coseismic rupture. To handle this uncertainty, we use the second node of the logic tree to evaluate three fault geometries that vary the distribution of slip on the updip edge of the rupture. The splay fault scenario is assigned a weight of 0.8 for the three largest size classes (L1, XL1, and XXL1), reflecting the authors’ consensus view that earthquakes with the greatest slip are more likely to involve simultaneous rupture of a splay fault (e.g., Plafker, 1972). Weighting was guided by the assumption that events with less slip (e.g., scenario SM1) are less likely to activate a splay fault (Table 3).

Simulated Cascadia Earthquake Deformation

Along-strike variations in the magnitude of coseismic slip and fault geometry produce corresponding variations in the 15 surface deformation models (Table 5). For example, variations in the vertical deformation profiles for the three earthquake scenarios weighted highest in the logic tree, models M1, M2, and M3 (Fig. 8; Table 2), reflect different fault geometries used in the rupture models (Fig. 9). Vertical deformation values for sister models of the same fault geometry scale proportionally with the simulated amount of coseismic slip (Table 5). In northern Oregon and Washington, the three rupture models produce very different deformation profiles (Olympics profile, Fig. 9A). In southern Oregon, at the latitude of Cape Blanco, the deformation profiles are similar because the model merges a splay fault into the deformation front (Fig. 9C), and there is no young outer accretionary wedge along the southern margin like there is off of Washington and northern Oregon. Along a profile near Bandon, the XXL1 scenario predicts maximum seafloor uplift of 10.2 m and maximum subsidence of 6.3 m (Table 5). In contrast, small earthquake scenarios (e.g., SM1) predict <2.6 m uplift and <1.6 m subsidence offshore Bandon.

The splay fault rupture model with the highest weight (scenario M1) simulates >7 m of peak seafloor uplift ∼70 km west of the Olympic Peninsula coast (Fig. 9A). Peak uplift, located offshore, gradually decreases southward to 5.2 m near Newport and 3.9 m off of Cape Blanco. The subsidence trough that parallels the Olympic Peninsula coastline is on land and exceeds 2 m, increasing to the north. Maximum subsidence decreases southward to ∼1.4 m in northern Oregon, where the trough swings offshore. The axis of subsidence remains offshore central Oregon until it crosses the coast at Cape Blanco, where subsidence reaches a maximum of ∼2.4 m.

The medium shallow buried rupture (scenario M2) simulates peak uplift of ∼2.5 m at ∼90 km off the coast of Washington. Seafloor uplift gradually increases southward as the width of the rupture zone narrows. Maximum seafloor uplift off of Newport is 3.4 m and slightly decreases to 3.2 m at Cape Blanco. The pattern of subsidence in the M2 model is very similar to the M1 model, except that maximum subsidence is slightly less (∼2.2 m) at Cape Blanco where the trough traverses on land.

The deep buried rupture with highest weight (scenario M3) simulates narrower seafloor uplift along the northern margin compared to the broad mound simulated in the shallow buried rupture model (M2), and peak uplift (>3 m) is shifted east to the boundary between the inner and outer accretionary wedges (Fig. 9A). Offshore of central and southern Oregon, where the outer wedge exhibits a steeper slope characterized by seaward-vergent structures, peak seafloor uplift is ∼3.6 m (Figs. 9B and 9C). Here, the peak of the bell-shaped mound is shifted east of the splay fault model (M1) and coincides with the break in slope that defines the boundary between the inner and outer wedges (Wang and Hu, 2006). Seafloor uplift shifted closer to land decreases tsunami arrival times.

Comparisons between Modeled Subsidence and Geological Observations

We compare longshore profiles of modeled coseismic subsidence for splay fault scenarios to estimates of coseismic subsidence during the A.D. 1700 Cascadia earthquake and earlier events compiled by Leonard et al. (2010) in Figure 10. Since the differences in coastal subsidence predicted by the buried rupture models generally fall within ∼0.5 m of the splay fault models (Fig. 9), comparisons in Figure 10B are relevant to all models. The subsidence profiles of all models are similar in shape to the profiles of subsidence from geological observations. The largest geologic estimates of subsidence reach a maximum of 1–2 m in southwestern Washington and >2 m in southern Oregon. In contrast, geological estimates of subsidence are lowest (<0.5–1 m) along the central Oregon coast.

The SM1 and M1 scenarios show the best agreement with paleoseismic observations compiled by Leonard et al. (2010), including the A.D. 1700 earthquake and earlier events. Although modeled subsidence of 1.5–2.5 m for scenario M1 slightly exceeds the range of observations compiled by Leonard et al. (2010) in southern Oregon, it falls below the maximum range of subsidence estimated at the Coquille estuary (Witter et al., 2003) and the Sixes River (Kelsey et al., 2002) (Fig. 10). New data for the A.D. 1700 earthquake by Hawkes et al. (2011) imply lower coseismic subsidence, which agrees well with the scenario SM1. Larger scenarios (L1, XL1, and XXL1; Fig. 10B) simulate greater amounts of subsidence in southern Oregon, which, in the case of the greatest scenarios, more than double geologic estimates from the last 6500 yr (Witter et al., 2003). Coastal paleoseismic records are too short at most localities to preserve evidence for the largest events in the early Holocene inferred from offshore turbidites.

The profile of coastal subsidence predicted by deformation models is relatively insensitive to the fault geometry used in our scenarios (Fig. 10C). Models M1, M2, and M3 all simulate 525–425 yr of slip with deformations that overlap the higher end of the subsidence range estimated from paleoseismology (Leonard et al., 2010).

Cascadia Tsunami Simulations

Variation of Hypothetical Tsunami Inundation and Runup

To aid tsunami hazard mitigation, we used the 15 Cascadia tsunami simulations to map tsunami inundation. Plate 1 shows simulated inundation for the five splay fault scenarios. Inundation exceedance lines delineate the likelihood (in percent) that Cascadia tsunamis will exceed the simulated inundation (Fig. 11). The exceedance lines are mapped by summing the cumulative logic tree weight at every grid node inundated by one or more tsunami simulations and then multiplying by 100. Inherent in the exceedance lines are weights for the earthquake size branches of the logic tree, which are based on the fraction of the 19 full-margin Cascadia earthquakes over the last 10,000 yr assigned to each size class. The exceedance lines allow decision makers to select specific inundation limits appropriate for different mitigation measures. For example, while residents may use the worse-case tsunami scenario for evacuation planning, policy makers may select a lower inundation limit for building codes. Inundation exceedance lines for the splay fault scenarios shown in Plate 1 are as follows: XXL1 <1%; XL1 ∼3%; L1 5%; M1 21%; and SM1 74%.

Comparisons of simulated tsunami inundation and runup at Bandon and New River show considerable variability (Fig. 12). Coseismic subsidence ranges from 1.5 to 6.3 m across both sites and controls whether or not tsunami waves breach coastal bluffs at Bandon. Runup elevations for the splay fault scenarios at Bandon vary between 5.0 and 16.7 m (Fig. 12A). Inundation stops at the coastal bluff for the smallest scenarios (SM1, M1, and L1), whereas the larger scenarios overtop the bluff and inundate 1–1.4 km inland. At New River, simulated runups vary from 0.4 to 19.5 m, and inundation reaches 1.2–3.4 km inland (Fig. 12B; Table 5). Waves produced by scenarios XL1 and XXL1 surge up topographic ramps at coastal foothills and reach runups higher than the maximum elevations at the shoreline, which vary between 5.6 and 14.8 m.

A longshore profile (Fig. 13) shows wide variations in maximum wave elevation at the shoreline and tsunami runup for scenario XXL1, revealing patterns related to nearshore wave shoaling and topographic effects. The 60 m isobath defines a broad concavity in the nearshore located below a 12-m-high mound reflected by maximum wave elevations (Fig. 13). The highest wave elevations at the MHHW shoreline exceed 20–25 m and occur predominantly along 45- to 60-m-high bluffs north and south of Cape Blanco and bluffs up to 60 m high at Fivemile Point (Fig. 13B). In these areas, and at 15- to 20-m-high bluffs in Bandon, the wave runup nearly matches the wave elevation at the shoreline (Fig. 13A).

The landward extent of inundation in lowlands is almost always greater than along bluff-backed shorelines. Simulated flooding in the Coquille River valley for the maximum considered scenario shows hydrologic effects of the tsunami that penetrate ∼40 km up river to near the town of Coquille (Fig. 13B). A review of model animations and time histories (Fig. 14) reveals that wave effects are negligible upriver of station 9 (Fig. 14A), ∼17 km from the mouth, and that flooding near the town of Coquille reflects coseismic subsidence predicted by the deformation model. In coastal lowlands and valleys near the Coquille River, Four-mile Creek, Floras Creek, and Elk River, runup descends to elevations 5–15 m lower than the maximum wave elevation at the shoreline. The inland extent of inundation stops short in areas protected by high coastal bluffs (e.g., Cape Blanco, Coquille Point, and Fivemile Point; Fig. 13B).

Time Series of Peak Tsunami Wave Elevation and Velocity

The arrival time of the first tsunami wave depression or peak at the shoreline depends on the shape and proximity of offshore uplift. Leading depression waves, which reflect subsidence of the seafloor, can amplify runup if the axis of maximum subsidence is significantly offshore (Tadepalli and Synolakis, 1994). No significant leading depression wave was generated for Cascadia tsunamis simulated for this project, because the axis of coseismic subsidence is near the coast in the Bandon area. The first tsunami peak is the largest wave in all scenarios at selected observation points. Peak waves arrive at the entrance of the Coquille estuary between 19 and 21 min after the simulations begin (Fig. 14). Arrival times of the first waves are nearly identical for scenarios that use the splay fault and shallow buried rupture models. Scenarios generated with the deep buried rupture model arrive 1 min earlier, because peak offshore uplift is located closer to shore. Maximum wave elevation and velocity time series are computed for selected sites along the Coquille River (Fig. 14).

DISCUSSION

Limitations

The hydrodynamic model SELFE solves 3-D nonlinear shallow-water wave equations using the finite-element method on unstructured grids. Like most hydrodynamic computer codes, the two-dimensional, depth-integrated configuration used to simulate tsunami waves simplifies tsunami flow dynamics. The model’s stable inundation algorithm uses zero bottom friction and viscosity, and because the model does not include realistic roughness parameters to account for buildings, vegetation, or other terrain variations, simulations of inundation may overestimate runup. All simulations were run using the MHHW tidal datum without accounting for ebb or flood tides, which may cause nonlinear amplification of tsunami waves (Myers and Baptista, 2001). Although negligible at the open coast, nonlinear tidal effects may influence wave dynamics in estuaries, especially for smaller tsunamis (Zhang et al., 2011). The model has performed satisfactorily when tested against three types of benchmarks (Zhang and Baptista, 2008) as recommended by Synolakis et al. (2007), including analytical solutions, laboratory tests, and field observations (National Tsunami Hazard Mitigation Program, 2012).

Earthquake rupture models, which supply information on the tsunami source, contribute the greatest uncertainties in near-source tsunami simulations. For the Cascadia rupture models developed here, we assume that Cascadia turbidite interevent times are proportional to the amount of slip released during past earthquakes. We selected four representative time intervals that encompass the uncertainty in turbidite interevent times to estimate different amounts of slip for a range of earthquake sizes (SM, M, L, and XL). Further uncertainty bears on whether or not the southern Cascadia plate boundary ruptures independently. However, the amount of coseismic slip released during earthquakes limited to southern segments of the megathrust is difficult to quantify. To account for this uncertainty, we assume that 15%–20% of the slip deficit is recovered during southern segment ruptures. This fraction of the southern Cascadia slip budget is loosely constrained by test simulations checked against tsunami deposits in Bradley Lake (Witter et al., 2012).

Our rupture models depict a simple regional rupture with a bell-shaped slip distribution along the full length of the margin. We assume 100% locking along a narrow patch of the plate interface that decays symmetrically updip and downdip consistent with earthquake deformation models at subduction zones (e.g., Wang, 2007; Wang and Hu, 2006). We do not model extremely high slip (>60 m) on the updip segment of the megathrust, like that observed during the 2011 Tohoku-oki earthquake (e.g., Yamazaki et al., 2011; Kodaira et al., 2012; Wei et al., 2012). Paleoseismic data along the Cascadia subduction zone neither support nor preclude earthquake scenarios involving such huge slips (Wang et al., 2013). Rather than the simple slip distributions modeled here, earthquakes on the Cascadia megathrust probably involve more complex rupture processes (Wang et al., 2013), like those of recent great earthquakes at other subduction zones in Sumatra (Chlieh et al., 2007), Chile (Moreno et al., 2009), and Japan (Iinuma et al., 2012; Shao et al., 2011; Simons et al., 2011).

Rupture models that partition slip to a shallow splay fault reflect geological evidence for a system of ∼30° landward-dipping thrust faults in the forearc observed in seismic profiles offshore Washington and northern Oregon (Fig. 5). The continuity of a splay fault system along an inferred length of 800 km is unclear, and its northern and southern ends are not well defined. The faults offset Holocene sediment and have clear surface expression in the seafloor (Goldfinger, 1994) that in some places defines the contact between the outer and inner accretionary wedge. However, we do not know whether one or more splay faults have ruptured with past megathrust earthquakes or whether any splay faults have broken independently.

Comparison of Regional Cascadia Tsunami Modeling Studies

Several studies have devised Cascadia earthquake scenarios to assess regional tsunami hazards in the Pacific Northwest (e.g., Hebenstreit and Murty, 1989; Ng et al., 1990; Whitmore, 1993; Priest et al., 1997, 2000; Geist, 2005). Early exercises proposed Mw 8.5–8.8 Cascadia earthquakes scenarios using rupture models with uniform slip on planar faults. Because near-source tsunamis are very sensitive to rupture details, including slip distribution and fault geometry, models using planar faults and uniform slip are hampered by oversimplified assumptions. Later models explored the effect of more detailed rupture parameters on the tsunami-generation process, including an ∼Mw 9 Cascadia earthquake scenario (Geist and Yoshioka, 1996; Priest et al., 1997, 2000; Geist and Dmowska, 1999; Geist, 2002, 2005), after discovery of the impacts of the 1700 Cascadia tsunami in Japan (Satake et al., 1996). These studies used more sophisticated earthquake rupture models with complex slip distributions on a curved fault plane, which avoid the inaccuracies of uniform slip models (Geist and Dmowska, 1999; Geist, 2002).

We ran regional tsunami simulations for the three largest Cascadia scenarios of this study (XXL 1–3) to estimate maximum nearshore wave amplitudes (wave elevation above tide at the 50 m isobath) along the Oregon and Washington coast. Figure 15 compares our results with peak nearshore tsunami amplitudes computed by Geist (2005), who used an earthquake source model similar to the “long-narrow” rupture geometry of Satake et al. (2003) and stochastic slip distributions intended to emulate rupture asperities (average slip 12.5–14.5 m with Mw = 8.8–8.9). Offshore southern Oregon, peak nearshore tsunami amplitudes calculated by Geist (2005) exceed 20–30 m and are roughly twice the amplitude of the highest nearshore tsunami waves simulated in this study (Fig. 15). However, the spatial positions of peaks and troughs in the longshore profiles show good agreement and likely reflect offshore bathymetric features on the continental shelf (e.g., banks and channels), as well as the configuration of the coastline. We find the differences in peak wave amplitudes surprising considering the average slip used in Geist’s sources is 57%–72% of the average slip used in our dislocation models (20–22 m for XXL scenarios; Table 4). One difference in modeling approach may explain the higher amplitudes estimated by Geist (2005): The crack model used by Geist (2005) concentrates slip at shallow depths beneath the outer wedge, leading to greater seafloor displacement in deeper water, which generates larger nearshore tsunami amplitudes compared to the lower nearshore tsunami amplitudes simulated by our models.

Comparisons of regional Cascadia tsunami simulations reveal two key rupture parameters with the greatest influence on tsunami size: (1) the amount of slip used in the rupture model and (2) whether or not the model diverts slip to a splay fault. Because these parameters are poorly known, we included end-member cases to explore the uncertainty of tsunami inundation. Rupture models incorporating a splay fault and rupture models with higher coseismic slip will produce greater tsunamis. The location of the updip limit of fault slip in buried rupture models, whether shallow near the trench or deeper downdip, produced less variation in the tsunami simulations.

Comparison of Simulated Tsunami Inundation to Paleotsunami Deposits

Sedimentary evidence of past Cascadia tsunamis near Bandon, Oregon (Kelsey et al., 1998, 2002, 2005; Witter et al., 2003), provides minimum estimates of the extent of inundation without considering changes in relative sea level, the configuration of late Holocene shorelines, the range of tide during past tsunamis, and gradual shoaling in estuaries. Symbols on Plate 1 differentiate sites that record sandy layers deposited by the A.D. 1700 Cascadia tsunami from sites where evidence of older tsunami deposits occurs. Most striking is the limited distribution of the A.D. 1700 tsunami deposit compared to the more inland distribution of older Cascadia tsunami deposits. Deposits from the most recent tsunami are only present near the mouths of the Coquille (Witter et al., 2003) and Sixes River estuaries (Kelsey et al., 1998). The A.D. 1700 tsunami also left a deposit in Bradley Lake that is much thinner than older deposits left by its predecessors (Kelsey et al., 2005). In contrast, older tsunami deposits in cores along Sevenmile Creek (Fig. 1B; Plate 1) extend ∼10 km from the mouth of the Coquille River estuary (Witter et al., 2003). The oldest tsunami deposits documented at the Sixes River are located ∼3.5 km from the river mouth (Kelsey et al., 2002). The more inland extent of older tsunami deposits likely reflects changes in the configurations of late Holocene estuarine shoreline sand barriers at the coast (Kelsey et al., 2002; Witter et al., 2003, 2012).

Comparisons among the distributions of Cascadia tsunami deposits and inundation predicted by tsunami simulations can identify scenarios that underestimate tsunami inundation, but they cannot preclude the largest scenarios. At the Coquille and Sixes Rivers, all Cascadia tsunami simulations flood beyond sparsely distributed deposits left by the A.D. 1700 tsunami (Plate 1). Tsunami deposits rarely mark the maximum extent of inundation, as shown by the 2011 Tohoku-oki tsunami, which spread sand over 57% to 90% of its inundation distance (Abe et al., 2012). Small and medium tsunami simulations (SM and M) envelop all but the most inland tsunami deposits laid down prior to the A.D. 1700 event. The largest scenarios in this study (L, XL, and XXL) encompass deposits found at sites farthest from the present shoreline, along Sevenmile Creek at the Coquille estuary and in an abandoned meander of the Sixes River (Plate 1). At Bradley Lake, the smallest Cascadia tsunami simulations that use 300 yr of accumulated slip (SM1–SM3) fail to reach the lake outlet. These results are consistent with simulations matched to tsunami deposits in Bradley Lake that required rupture models with minimum slip deficits of 360–400 yr to inundate the lake (Witter et al., 2012). The M2 and M3 tsunami simulations use buried rupture models with a slip deficit of 450 yr, yet these simulations do not inundate the lake because the model grid incorporates light detection and ranging (LiDAR)–derived topography that includes extensive modern foredunes and a shoreline farther west of the lake than at the time of the A.D. 1700 tsunami (Witter et al., 2012). The M1 scenario and all larger simulations inundate Bradley Lake (Plate 1).

Evidence Supporting a Range of Cascadia Earthquake Scenarios

The most likely scenario based on total scenario weights in the logic tree (Fig. 8; Table 3) is the medium rupture model that partitions slip onto a splay fault (M1). Coastal subsidence for the M1 model agrees well with paleoseismic estimates of coastal subsidence (Leonard et al., 2010; Hawkes et al., 2011; Fig. 10). Although important details differ, coseismic slip used for the M1 model (14–19 m) is similar to the amount of slip in the “most-likely” earthquake used by Satake et al. (2003) to model runup in Japan caused by the A.D. 1700 Cascadia tsunami. Their preferred “long-narrow” model simulated a Mw 9.0 earthquake with uniform slip on a 1100-km-long, full-slip zone with a linear downdip transition to zero. In contrast, our model features a bell-shaped slip distribution that decreases slip updip and downdip using the slip function of Wang and He (2008) and produces higher seafloor uplift and larger local tsunami waves.

Cascadia rupture models that invoke a splay fault system in the accretionary wedge along the northern margin may be analogous to ruptures that involved splay faulting along other convergent margins. Splay faults ruptured during the 1964 Alaska earthquake (Plafker, 1972; Johnson et al., 1996), the 1963 Kuril Islands and 1975 eastern Hokkaido earthquakes (Fukao, 1979), and possibly the 1944 Mw 8.1 Tonankai earthquake (Tanioka and Satake, 2001; Park et al., 2002; Moore et al., 2007; Baba et al., 2006; Strasser et al., 2009).

Our analyses include several outsized rupture models that may be very unlikely but provide reasonable maximum limits on potential rupture and resulting tsunami runup, information that is critical for tsunami hazard mitigation. The displacements predicted by the maximum rupture models are consistent with historical great earthquakes along other convergent margins, and the fault geometries accurately represent the structure of the Cascadia forearc.

The spectrum of Cascadia earthquake models presented here also may be applied as a suite of source scenarios for probabilistic tsunami hazard assessments. Our interpretation that the turbidite record reflects great earthquakes of variable size provides a means by which to estimate recurrence rates for small through extra-large scenarios. Of the 19 full-margin ruptures interpreted from the ∼10,000 yr turbidite paleoseismic record, five are classified as small (SM), 10 as medium (M), three as large (L), and one as extra-large (XL). This distribution of size classes corresponds to recurrence rates as follows: SM, 1/2000 yr; M, 1/1000 yr; L, 1/3333 yr; and XL, <1/10,000 yr. The recurrence rate for XXL earthquake scenarios is not known.

Slip distributions of historical megathrust earthquakes offer modern analogs that support the maximum slip estimates used in our models. Peak slip estimates for historical earthquakes include: 15–25 m offshore northern Sumatra during the 2004 Mw 9.1 Sumatra-Andaman Islands earthquake (e.g., Ammon et al., 2005; Lay et al., 2005; Chlieh et al., 2007; Fujii and Satake, 2008; Grilli et al., 2007); 20–25 m on two slip patches of the 1964 Mw 9.2 Prince William Sound earthquake (Johnson et al., 1996; Suito and Freymueller, 2009); and more than 40 m of peak slip during the 1960 Mw 9.5 southern Chile earthquake (Barrientos and Ward, 1990; Moreno et al., 2009). Estimates of peak coseismic slip for the 2011 Tohoku-oki earthquake exceed 60 m near the Japan Trench (e.g., Ozawa et al., 2011; Wei et al., 2012; Yamazaki et al., 2011). If subduction zones store energy over long intervals that span several earthquake cycles, higher slip deficits than inferred in this study might be achieved (i.e., supercycles; Sieh et al., 2008; Goldfinger et al., 2013a).

The 2011 Tohoku-oki earthquake highlights another uncertainty in Cascadia: whether or not rupture can break the seafloor, particularly along the southern margin. The huge coseismic slips at shallow depth near the Japan Trench (Iinuma et al., 2012) should motivate future seismic surveys to seek evidence, if present, of shallow ruptures along Cascadia’s deformation front.

SUMMARY AND CONCLUSIONS

We have presented a strategy for full-margin assessment of Cascadia tsunami inundation hazards by implementing a suite of megathrust earthquake scenarios that were ranked using a logic tree. We used the hydrodynamic model SELFE (Zhang and Baptista, 2008) to assess the potential variability of tsunami inundation exemplified by a series of 15 simulations along the southern Oregon coast near Bandon. We found that the two primary controls on tsunami inundation include (1) the amount of maximum slip used in the rupture model, and (2) whether or not slip is diverted to a shallow splay fault that amplifies seafloor uplift. We mapped the inundation simulated for the 15 tsunami scenarios and calculated inundation exceedance lines (in percent), which will allow decision makers to select appropriate scenarios for specific mitigation purposes.

Evaluating the likely variability of tsunami inundation has important implications for tsunami hazards mitigation. Although the maximum scenario (XXL1) of Mw 9.1 with 36–44 m maximum slip represents a very rare event, the state of Oregon has adopted it as the benchmark for tsunami evacuation maps. In contrast, the most likely Cascadia tsunami identified in our assessment (scenario M1) involves a Mw 8.9 earthquake with 14–19 m maximum slip. Because they encompass 80%–95% of simulated tsunami inundation scenarios, we suggest scenario M1, or the larger L1 scenario (22–30 m slip), as a credible alternative to consider for future revisions to coastal construction standards (e.g., Olmstead, 2003), land-use planning, and for engineering design of critical structures along the coast.

Finally, future studies should investigate the evidence for trench-breaking rupture, which led to the huge horizontal and vertical displacements of the seafloor near the Japan Trench (Wei et al., 2012) and continue to glean information from paleoseismology. More accurate estimates of coseismic deformation from paleogeodetic studies (Wang et al., 2013) and evidence for earthquake shaking in lake sediment (Morey et al., 2012) may help to refine estimates of rupture length for the next generation of coseismic slip models.

Supporting grants from the National Oceanic and Atmospheric Administration (NA08NWS4670028 and NA09NWS4670014) were awarded through the National Tsunami Hazard Mitigation Program. Critical and constructive reviews by H. Kelsey, P. Molnar, D. Brothers, and an anonymous reviewer strengthened the paper. This research is a contribution toward National Science Foundation award EAR-0842728 and International Geoscience Programme (IGCP) 588 coastal change.