Assessing Margin-Wide Rupture Behaviors Along the Cascadia Megathrust With 3-D Dynamic Rupture Simulations.

From California to British Columbia, the Pacific Northwest coast bears an omnipresent earthquake and tsunami hazard from the Cascadia subduction zone. Multiple lines of evidence suggests that magnitude eight and greater megathrust earthquakes have occurred - the most recent being 321 years ago (i.e., 1700 A.D.). Outstanding questions for the next great megathrust event include where it will initiate, what conditions are favorable for rupture to span the convergent margin, and how much slip may be expected. We develop the first 3-D fully dynamic rupture simulations for the Cascadia subduction zone that are driven by fault stress, strength and friction to address these questions. The initial dynamic stress drop distribution in our simulations is constrained by geodetic coupling models, with segment locations taken from geologic analyses. We document the sensitivity of nucleation location and stress drop to the final seismic moment and coseismic subsidence amplitudes. We find that the final earthquake size strongly depends on the amount of slip deficit in the central Cascadia region, which is inferred to be creeping interseismically, for a given initiation location in southern or northern Cascadia. Several simulations are also presented here that can closely approximate recorded coastal subsidence from the 1700 A.D. event without invoking localized high-stress asperities along the down-dip locked region of the megathrust. These results can be used to inform earthquake and tsunami hazards for not only Cascadia, but other subduction zones that have limited seismic observations but a wealth of geodetic inference.

Introduction

The Cascadia subduction zone megathrust dominates earthquake hazard in the United States Pacific Northwest. It is oft‐cited that the probability of a magnitude ∼9 (M9) event occurring in the coming decades is between 10%–14% (Petersen et al., 2014). The most recent megathrust rupture in Cascadia occurred in 1700 A.D. and generated a transoceanic tsunami (Heaton & Hartzell, 1987). Matching amplitudes of historical tsunami records from Japan requires a magnitude between M8.7–9.2 for this earthquake (Satake et al., 1996, 2003). While 321 years have elapsed since this last event, the Holocene (<12 kya) earthquake record onshore and offshore documents even older M > 8 megathrust events (e.g., Atwater & Griggs, 2012; Leonard et al., 2010; Kemp et al., 2018; Goldfinger et al., 2012, 2017). To recognize which geologic processes drive partial versus margin‐wide ruptures in Cascadia (e.g., segmentation), a synthesis of existing observations is needed to see how differences in subduction zone structure (asperities) or frictional behavior (aseismic deformation) may ultimately provide barriers to rupture propagation (Philibosian & Meltzner, 2020). Megathrust segmentation has played a role in other large ruptures (e.g., Alaskan and Chilean subduction zones; Ichinose et al., 2007; Moreno et al., 2011) and influences the nature of stress release.

Geological and Geophysical Inferences on the State of Megathrust Segmentation

Several geological and geophysical observations suggest the Cascadia megathrust exhibits along‐strike segmentation. For instance, there are systematic changes in the accretionary wedge backstop geometry, seismicity, and interseismic slip patterns (e.g., Stone et al., 2018; Bartlow, 2020; Watt & Brothers, 2020) that may indicate coseismic rupture patterns will also be variable along‐strike. The strongest observational constraints that may inform our understanding of future great earthquakes come from paleoseismic and geodetic observations. Underwater turbidite deposits, which can be generated from submarine landslides induced by strong ground‐shaking during megathrust earthquakes, have been extensively used to map along‐strike rupture extents (Goldfinger et al., 2003, 2012, 2017; Figure 1a). Analysis of the timing and spatial extents of turbidite deposits suggests that the recurrence interval (RI) between megathrust earthquakes could vary along the Cascadia margin. In particular, the RI estimated for northern Cascadia (>∼46latitude) exceeds 400 years whereas it is estimated to be less than 200 years for the southern portions of Cascadia (∼43latitude; Goldfinger et al., 2017). These differences in strain release and rupture extent suggested by offshore (e.g., Goldfinger et al., 2012) and onshore (e.g., Priest et al., 2017) hazard studies imply variable strain accumulation times along the Cascadia margin.

Figure 1

Cascadia subduction zone study area. (a) Megathrust segmentation (segments are separated by red lines) suggested from offshore turbidite deposits (Goldfinger et al., 2017) with corresponding estimated segment recurrence intervals in years. Primary morphotectonic regions identified by Watt and Brothers (2020) are superposed (blue dashed lines). (b) Gaussian and (c) Gamma coupling models from Schmalzle et al. (2014) projected onto the Slab2 megathrust geometry. The Gaussian coupling model assumes interseismic creep at shallow megathrust depths whereas the Gamma model assumes high strain accumulation. The inferred region of the creeping segment is denoted by a yellow box. Magenta stars denote rupture initiation locations in our dynamic rupture models. Thick white lines are megathrust depth contours (kilometers). JdF = Juan de Fuca plate.

Decadal scale interseismic velocities measured at the Earth's surface by Global Navigation Satellite System (GNSS) networks, tide gauge, and leveling data also find significant variations in coupling (slip deficit) distribution along the margin. Regions in northern and southern Cascadia have higher coupling suggesting they are accumulating strain that may be released in a future great earthquake (Schmalzle et al., 2014; Li et al., 2018; Yousefi et al., 2020; Figures 1b and 1c). However, the use of geodetic coupling inversions to place bounds on the future down‐dip or along‐strike rupture extent is complicated by heterogeneous frictional properties (Boulton et al., 2019) or the potential presence of stress shadows (Almeida et al., 2018; Hetland & Simons, 2010) and other factors (e.g., off‐fault deformation). Therefore, the down‐dip extent of coupling and coseismic rupture may differ even though these inversion results are the best available constraint on potential stress distributions for the Cascadia megathrust (Wang & Tréhu, 2016).

Another piece of evidence for segmentation comes from the behavior of episodic tremor and slip (ETS) events along the megathrust. In GNSS displacement records, ETS manifests as transient reversals in displacement indicative of slip on or near the megathrust at depths between 30 and 40 km. These slow slip episodes are often accompanied by a weak seismic signature known as nonvolcanic tremor (Rodgers & Dragert, 2003). The character of ETS events varies significantly along the margin. The northern (i.e., >47) and southern (<43) sections host more frequent slip episodes with average recurrence intervals of 10 and 14 months and have higher tremor density, whereas the central section of the megathrust hosts ETS approximately every 19 months (Brudzinski & Allen, 2007; Wech & Creager, 2008). Studies of the ETS source region find that the phenomenon occurs in regions of significantly elevated Vp/Vs ratios (e.g., Audet et al., 2009; Delph et al., 2021) and that tremor, and constituent low‐frequency earthquakes, are extremely sensitive to small magnitude stress changes such as those from the solid Earth tides (Royer et al., 2015). Collectively, these observations suggest that pore fluid pressures are nearly lithostatic in the ETS source region. The Cascadia megathrust also features a transition zone at depth that separates the ETS region from the region that is conventionally considered to be locked (20 km; Hyndman, 2013). Known as the gap, this spatial disconnect in slip behavior is also found in other subduction zones (Gao & Wang, 2017); the frictional behavior and shear stress accumulation levels in the gap may or may not allow for deeper rupture (Ramos & Huang, 2019).

Cascadia Earthquake Source Models

What are ways to anticipate how a future Cascadia megathrust earthquake may behave? One way to assess the hazard posed by large seismic events in Cascadia is to use kinematic rupture simulations. Kinematic rupture simulations are commonplace due to the straightforward relationship between fault slip and the recorded elastic displacement field once the Green's functions are known, allowing these types of models to be run at lower computational cost. Using the kinematic framework, potential locations of strong‐ground motion sources along the fault, sedimentary basin amplification or tsunami generation have been assessed (Delorey et al., 2014; Frankel et al., 2018; Melgar et al., 2016; Olsen et al., 2008; Roten et al., 2019; Wirth et al., 2018, 2019; Wirth & Frankel, 2019). Most of these kinematic rupture models calibrate first‐order rupture parameters (slip, slip‐rate, rise‐time or rupture speed) from the few large megathrust earthquakes observed in other subduction zones (e.g., M8.2 2003 Tokachi‐Oki, M9.2 2004 Sumatra‐Andaman, M8.8 2010 Maule, M9.1 2011 Tohoku‐Oki). Kinematic simulations provide important constraints on strong ground motions felt onshore, but because they must assume a slip distribution before computing the elastic wavefield they cannot answer what controls the final rupture size.

To account for source physics, fully dynamic rupture simulations can be used to investigate what controls the final rupture size and kinematic rupture properties like rupture speed. Dynamic rupture simulations are self‐consistent, physics‐based numerical models that describe the entire earthquake rupture process (nucleation, propagation, and arrest) that is coupled to a constitutive fault friction law (e.g., Madariaga & Olson, 2002). To date, 2‐D dynamic rupture models in Cascadia have focused on tsunami generation (e.g., Lotto et al., 2018) or how frictional and stress conditions in the transition zone may influence down‐dip rupture extent (Ramos & Huang, 2019).

Here we develop 3‐D dynamic rupture simulations to explore how variable strain accumulation rates, frictional behavior and hypocenter location influence megathrust rupture dynamics. We will use geodetic coupling results from Schmalzle et al. (2014), who utilized GNSS time series information spanning several decades for their coupling inversions. These coupling distributions are selected because they represent two possible end‐member scenarios for strain accumulation near the deformation front: either there is interseismic creep at shallow megathrust depths (hereafter referred to as the Gaussian coupling model; Figure 1b) or it is fully coupled (hereafter referred to as the Gamma coupling model; Figure 1c). We also use the Schmalzle et al. (2014) coupling models for consistency with the material rheology we adopt, which is a linearly elastic medium. The first order features of viscoelastic coupling models (e.g., Li et al., 2018; Yousefi et al., 2020) are not drastically different from the Gamma slip deficit rate distribution. Moreover, Li and Liu (2021) found that the Schmalzle et al. (2014) coupling models, when incorporated into quasi‐dynamic earthquake simulations, provide a stronger fit to the 1700 A.D. surface measurements. These coupling models will be used to constrain the dynamic stress drop, which is defined as the difference between the initial shear stress and dynamic fault strength. Dynamic stress drop is a key parameter determining how much energy is available for rupture propagation in dynamic rupture simulations (Kanamori & Rivera, 2004; Yang, Yao, He, Newman, et al., 2019, 2019b). Our dynamic stress drop levels are further constrained by strain accumulation times and segment locations adopted from paleoseismic studies (i.e., Goldfinger et al., 2017; Figure 1a). We compare the resulting coseismic uplift and subsidence patterns to available paleoseismic measurements and discuss which classes of models allow margin‐wide ruptures to develop. We find the final earthquake size is sensitive to earthquake nucleation location (e.g., northern vs. southern Cascadia) and the distribution of relative dynamic stress drop. The principal control on margin‐wide rupture, when using these particular end member geodetic coupling models, is the relative dynamic stress drop amplitude in the central Cascadia region (∼43–47latitude). The results also suggest that Gamma coupling models tend to produce larger earthquakes, even if shallow subducted sediment has a slip‐strengthening or velocity‐strengthening frictional behavior. Another intriguing question is if geodetic coupling models can inform our understanding about the 1700 A.D. earthquake when incorporated into a dynamic rupture simulation. To that end, we also present several rupture simulations that provide a close fit to the 1700 A.D. event.

Methodology

We solve for 3‐D elastodynamic earthquake rupture using SeisSol, a powerful open‐source software package that implements the Arbitrary high‐order DERivative‐Discontinuous Galerkin (ADER‐DG) approach to simulate wave propagation coupled to spontaneous dynamic rupture (De la Puente et al., 2009; Heinecke et al., 2014; Pelties et al., 2012; Uphoff et al., 2017). The capability of SeisSol to solve for complex source dynamics and incorporate realistic geometric features, such as bathymetry, topography and fault zone structure (e.g., Ulrich et al., 2019; Wollherr et al., 2019) nicely lends itself to our purposes of investigating how heterogeneous megathrust stresses influence rupture behavior.

We generate an unstructured 3‐D tetrahedral mesh for the Cascadia subduction zone that spans over 1,100 km along‐strike (39.0–51.0° latitude, −127.5 to −121.0° longitude) and we use static refinement to increase resolution locally. Note that we apply a coordinate transformation from longitude/latitude to easting/northing (X, Y) in kilometers that is centered at −128.00, 46.80 (longitude, latitude) for the entire model domain geometry. The average on‐fault element edge size (h) is 2.5 km, and the maximum depth of the fault mesh is 50 km (Hayes et al., 2018) and includes over 440,000 unstructured triangular elements. We account for the large‐scale variations in the free‐surface geometry by meshing the ETOPO1 topography and bathymetry data set to ∼1 km average element size near the coastline. In all of our simulations, we use ADER‐DG with fifth order accuracy (polynomial order p = 4) in time and space.

We ensure simulation results are sufficiently resolved by following the procedure established in Wollherr et al. (2018) to estimate the process zone, the region behind the rupture front where the fault strength drops from its static to dynamic level. For the 2.5‐km fault mesh, the median process zone width) is ∼1.1 km. The recommended number of elements needed to resolve in a purely elastic setup with depth‐dependent initial conditions is 2–3 (p = 4). The quadrature points approach utilized in SeisSol (Pelties et al., 2014) ensures each element edge length is sampled p + 2 times. Given our setup, is sampled by ∼2.7 elements which is within the recommended range. The expected relative percent error in the rupture arrival time, peak slip‐rate, and final slip are 0.9, 8.32, and 0.71, respectively (Wollherr et al., 2018). While the peak slip‐rate relative error is slightly larger than the 7% recommended by Day et al. (2005) for elastic rupture problems, we compare our model‐predicted slip and rupture size to higher resolution meshes with h = 1 km and h = 0.5 km and observe negligible changes, which gives us confidence that these first order rupture features are correctly resolved. The highest resolution mesh has more than 50 million elements and requires 22 h on 40 nodes of the supercomputer SuperMUC‐NG at the Leibniz Supercomputing Centre, Germany.

Constraining Dynamic Rupture With Geodetic Coupling Models

In our simulations, potential shear stress distributions are informed by geodetic coupling models. The Schmalzle et al. (2014) inversion for slip rate deficit was performed with respect to a Cascadia megathrust geometry predating Slab2 (McCrory et al., 2012) and as such, we first map the geodetic coupling models to our megathrust geometry through a bilinear interpolation using the cartesian horizontal plane coordinates. But the effect of this transformation does not distort the main features of the coupling models (Figures 1b and 1c).

We define the parameter T as the time needed for a certain level of slip deficit to accumulate on a section of the megathrust. The product of slip rate deficit (coupling) and T is slip deficit. T should not be interpreted as the RI, but rather as another way to quantify relative dynamic stress drop along the megathrust. Segments of the megathrust with higher T levels can be thought of as regions where interseismic strain is building. From these slip distributions, we estimate the static stress change using Poly3D, a three‐dimensional, polygonal element, displacement discontinuity boundary element method, which accounts for nonplanar megathrust geometry and the free‐surface effect due to buried slip (Thomas, 1993).

Initial shear stress is then estimated by adding the static stress change (the model output from Poly3D) to the dynamic fault strength. Calculating the initial shear stress in this manner is known as the complete stress drop assumption and assumes that slip deficit is accumulated linearly in the along‐dip fault dimension and will be entirely released during coseismic rupture (Yang, Yao, He, Newman, et al., 2019, b; Hok et al., 2011). This shear stress distribution is first resampled to an average grid spacing of ∼ 3 km and then linearly interpolated onto the fault mesh. We note that we initialize stress values and friction parameters in SeisSol with a high order sub‐element resolution, such that these parameters can be sampled by (p +2)2 Gaussian integration points within a tetrahedral element face (see section 4 in Pelties et al., 2014). This ensures that the gridded input data are accurately resolved without the need to adjust on‐fault element edge lengths. For all dynamic rupture simulations considered, we compare the results to the 1700 A.D. subsidence measurements along the coast where available (Wang et al., 2013), and to recorded subsidence amplitudes from other ∼ M9 earthquakes (e.g., 2011 Tohoku, 1960 Chile, 1964 Alaska).

Material Properties and Fault Strength

Wave propagation is simulated within a layered and linearly elastic medium where the elastic moduli (lame parameters) vary as a function of depth. The average 1‐D velocity structure is taken from the Cascadia 3‐D Community Velocity Model (3D‐CVM) for P and S waves (Figure 2a; Stephenson et al., 2017). Since the goal of this study is to calculate upper plate deformation and rupture extent (along‐dip and along‐strike) for a given dynamic stress drop distribution, we believe this is a satisfactory simplification to make. We estimate that we can resolve a cutoff seismic frequency up to ∼0.4 Hz in the near fault region. High frequency (>1 Hz) broadband ground motions can be calculated at a higher computational cost if an appropriate 3‐D velocity model is utilized. The current 3D‐CVM was developed with respect to an older Cascadia subduction zone geometry (i.e., McCrory et al., 2012) and thus, we leave direct extrapolation of this 3‐D velocity model to our model geometry for future work.

Figure 2

Material properties, strength, and frictional conditions for dynamic rupture simulations. (a) Smoothed 1‐D CVM velocity model for Cascadia (Stephenson et al., 2017). (b) Effective normal stress extended beyond 40 km depth from Ramos and Huang (2019). (c) Dynamic and static frictional coefficients with depth. (d) Frictional cohesion for the Gaussian and Gamma coupling models.

Effective normal stress accounts for pore pressure counteracting vertical lithostatic stress on the fault. We use the depth‐dependent effective normal stress distribution for Cascadia presented in Ramos and Huang (2019) that includes low strength levels (1 MPa) in the ETS region (Figure 2b). These incredibly low stress conditions in the ETS region are supported by observations on the sensitivity of tremor and low‐frequency earthquakes to small magnitude stress changes (e.g., Rubinstein et al., 2007; Royer et al., 2015), stress orientations in the ETS region (e.g., Newton & Thomas, 2020), and low stress drops of ETS events (e.g., Gao et al., 2012). For lack of in‐situ fault stress information, we assume a linear stress gradient above and below the locked region (10–20 km depth) that are consistent with other Cascadia megathrust simulations (Li & Liu, 2016; Liu & Rice, 2009). Such assumptions are simple but allow us to focus on how heterogeneous shear stresses on the megathrust contribute to first order rupture characteristics. Note that we do not explicitly calculate the pore fluid response from inertia during rupture simulation.

Fault Friction Law

The physics controlling the inelastic breakdown process in our dynamic simulations is given by a nonsingular linear slip‐weakening friction law (Ida, 1972; Palmer & Rice, 1973). This constitutive friction law allows us to idealize rupture as a propagating shear‐crack. It is described by the static () and dynamic () friction coefficients and a critical slip‐weakening distance (D ).

We set  = 0.6 and  = 0.1 within the locked region of the megathrust (5 km depth 20 km) (Figure 2c). Because Ramos and Huang (2019) showed that rupture can penetrate the gap or generate strong free‐surface reflections if its frictional behavior is slip‐weakening at depths >25 km and at depths <5 km (together with a highly negative stress drop), we set equal to or above in these regions (Figure 2c). D is set to a constant level of 1 m or 2 m. D c = 2 m is selected in the dynamic rupture model in which the stress and strength conditions of Ramos and Huang (2019) are extrapolated along strike, for consistency with the 2‐D dynamic rupture simulations. D c = 1 m is used for the dynamic rupture models based on the heterogeneous geodetic coupling prestress distributions. Our range of D values are consistent with those used in slip‐weakening simulations of the Tohoku‐Oki earthquake, which constrained D using the frequency range of back‐projection results (Huang et al., 2014). We make minimalistic assumptions for cohesion in the upper 5 km of the megathrust (Figure 2d). Due to the nearly zero dynamic stress drop amplitudes near the deformation front for Gaussian coupling models, the cohesion gradient can be low (Figure 2d). But in the case of the Gamma coupling models, relatively higher cohesion levels (∼5 MPa average) are locally needed at shallow fault depths to prevent fault failure at the start of the simulation (Figure 2d).

Rupture Initiation

Fault pre‐stress conditions influence the estimated critical nucleation size when using a linear slip‐weakening friction law. The theoretical critical nucleation radius that permits spontaneous dynamic rupture to initiate in a 3‐D linearly elastic and homogeneous media has been derived by Day (1982) and is given by, where G is the shear modulus, S is the relative fault strength defined as the ratio between strength excess (static fault strength minus initial shear stress) and dynamic stress drop ().

Expression (1) provides a sufficient means to initiate and sustain dynamic rupture propagation for the 3‐D dynamic rupture model that is adapted from 2‐D dynamic rupture simulations presented in Ramos and Huang (2019). For the prestress distributions derived from the heterogeneous coupling models, we determine the best numerical nucleation size through a trial‐and‐error approach. We find that critical nucleation radii are within ∼10% of the theoretically predicted value calculated from Equation 1. Rupture initiation is prescribed by a smooth space and time dependent function, leading to an imposed rupture velocity that decreases away from the hypocenter and allows a gradual transition from forced to spontaneous rupture (Harris et al., 2018). Rupture nucleation locations are chosen within the areas containing the highest dynamic stress drop distribution (see Figures 1b and 1c). Each dynamic rupture simulation is run for 420 s (7 min) to allow seismic waves to propagate to the edge of the model domain.

Results

Translating 2‐D Rupture Simulations to 3‐D

A 3‐D dynamic rupture model that assumes a simple dynamic stress drop profile extending along‐strike of the megathrust is presented in Figure 3. Previously developed 2‐D dynamic rupture simulations (Ramos & Huang, 2019) were relative to a specific location in northern Cascadia, which is where we initiate rupture (Figure 3a). Such a laterally uniform dynamic stress drop distribution is unlikely given observations of geophysical and geological megathrust segmentation (e.g., Watt & Brothers, 2020). However, we develop such a simulation to demonstrate 1) what a megathrust event would appear as if there was a strong gradient in shear stress‐rate from the locked to gap regions (20–30 km depth) across the margin and 2) how this scenario would influence coastal subsidence amplitudes.

Figure 3

Results of the dynamic rupture model in which the stress and strength conditions of Ramos and Huang (2019) are extrapolated along‐strike (RH). (a) Along‐strike dynamic stress drop and considered epicenter location. (b) Final megathrust slip‐distribution and moment‐magnitude. The black dashed lines indicate the 10 and 20 km depth whereas the solid black line denotes the coastline. (c) Coseismic uplift (red) and subsidence (blue) along the Cascadia margin. Squares signify the coastline. (d) Model predicted (black squares; same as panel C) and paleoseismic observations of estimated subsidence during the 1700 A.D. rupture (Wang et al., 2013).

In spite of the low stress drop (<−8 MPa) at depths shallower than 10 km and slip‐strengthening friction, coseismic slip is able to reach the deformation front with amplitudes exceeding 60 m in most locations along‐strike (Figure 3b). The along‐strike variation of slip at the deformation front exhibits two peaks north and south of the hypocenter ‐ even though the initial dynamic stress drop distribution is laterally invariant, the final coseismic slip pattern is not (Figure S1). This might be attributed to changes in the along‐strike megathrust dip angle. There are also small amounts of slip (<5 m) in the gap region. The coseismic hinge‐line, separating regions of subsidence from regions of uplift, is entirely offshore (Figure 3c). Subsidence levels exceeding 5 m are observed along most of the coastline (Figures 3c and 3d). This exceeds subsidence measurements from the 1700 A.D. event (Wang et al., 2013) by several times because the earthquake magnitude is 9.55 (Figure 3d). Such subsidence amplitudes are also much larger than the maximum levels observed for the 2011 M9.0 Tohoku‐Oki (∼1.1 m, Hashima et al., 2016), 1964 M 9.4 Alaska (∼2.4 m, Plafker, 1972) or the largest ever recorded event, the 1960 M9.5 Chile Earthquake (∼2.7 m, Plafker & Savage, 1970).

Interestingly, this model generates a down‐dip rupture front that can reach and “jump” across the gap region, despite the negative stress drop in the gap combined with a slip‐neutral frictional behavior (Figure S2). Due to dynamic stress perturbations carried by seismic waves (“dynamic unclamping”, Oglesby & Mai, 2012; Figure S2), this down‐dip rupture front is most likely triggered by temporal stress changes and is made possible by the incredibly low static fault strength here (i.e., 0.6 MPa).

Uniform Gaussian and Gamma Coupling Models

We now explore dynamic rupture scenarios based on the Gaussian and Gamma coupling models. We start with simulations that assume a uniform T level (Figures 4a and 4b). T is set to 320 years, the time elapsed since the most recent event (Goldfinger et al., 2017). In this parameterization, the highest dynamic stress drop amplitude is located in the northern Cascadia region for both Gaussian and Gamma distributions (Figures 4a and 4b), which is where spontaneous rupture is initiated. The location of highest dynamic stress drop is not coincident between the Gaussian and Gamma coupling models, and hence the hypocenter locations are slightly different. Uniform T for both coupling models generates margin‐wide rupture with coastal subsidence amplitudes that again exceed 1700 A.D (Figure 4c). The Gamma coupling model has higher dynamic stress drop than the Gaussian model near the deformation front, which leads to a 1 to 2‐m difference in subsidence amplitude for the northern (0 < Y < 200 km) region of the megathrust (Figure 4c). These subsidence amplitudes, while lower than the 2‐D extrapolated model, still surpass the estimated subsidence amplitudes of the largest recorded global megathrust earthquakes (i.e., 1960 Chile, 1964 Alaska). This result demonstrates that the uniform T coupling model overestimates the amount of slip deficit accumulated since 1700 A.D.

Figure 4

Comparison between Gaussian and Gamma dynamic stress drop distributions and the resultant subsidence patterns assuming the maximum strain accumulation time (T) of 320 years (i.e., time since the last great earthquake in 1700 A.D.) (a) Dynamic stress drop distribution for the Gaussian coupling model. (b) Dynamic stress drop distribution for the Gamma coupling model. Both ruptures are nucleated in northern Cascadia (magenta star). (c) Model predicted subsidence along the coastline compared to 1700 A.D. measurements.

When comparing the average along‐dip distribution of dynamic stress drop between the smooth (Figure 3) and Gaussian/Gamma 3‐D models (Figure 4), we note that the smoother model extends slightly deeper into the gap region. The amplitude of coseismic subsidence is probably more strongly controlled by the amount of dynamic stress drop closer to the coastline, though a larger overall stress drop along the seismogenic portion of the fault also generates a higher subsidence amplitude (Figure S3). A point to note is that the region of higher relative dynamic stress drop in the northern Cascadia region (0 Y 200 km) is also where there are limited paleoseismic measurements from 1700 A.D. Thus, while geodetic coupling models are well constrained here, the few along‐strike subsidence measurements limit rigorous comparison to physics‐based model predictions.

Segmented Gaussian and Gamma Coupling Models

We find that in order to produce coseismic uplift and subsidence amplitudes more consistent with the paleoseismic Cascadia measurements and data from other megathrust earthquakes (i.e., ±2 m), we must prescribe along‐strike variations of T, with T amplitudes lower than 320 years for a particular segment. This is especially needed for the northern and southern regions of the Cascadia megathrust, where both the Gaussian and Gamma coupling models predict higher subsidence amplitudes than observed if T is set to 320 years. We refer readers to the discussion section on the possible meaning of T values <320 years.

Our partitioning of the margin is informed by paleoseismic (Goldfinger et al., 2012, 2017), ETS (Brudzinski & Allen, 2007), and morphotectonic studies (Watt & Brothers, 2020) in Cascadia. The following dynamic rupture models are parameterized using at least three segments. This choice is conservative ‐ we found through trial‐and‐error that two segment models cannot match first‐order 1700 A.D. subsidence patterns as well as three‐segment models. We note that some geologic models may suggest up to five segments (e.g., Goldfinger et al., 2017; Wang et al., 2013) and thus there may be multiple ways to partition T levels along‐strike. We have compiled all rupture model results in Table 1.

Table 1

Summary of Dynamic Rupture Simulations for Segmented Models Generated From Gaussian or Gamma Coupling Distributions

Model no.	Nucleation location	Rupture length (km)	Mw	T Levels in northern, central, and southern segments (years)	Margin‐wide?	Central segment width (km)
1	Northern Cascadia, Gaussian	1,069	9.11	220, 250, 200	Yes	441
2		1,069	9.10		Yes	323
3		650	9.02		No	266
4		480	8.80		No	214
5	Gamma	1,069	8.81	94, 126, 82	Yes	266
6		1,069	8.80		Yes	214
7		1,069	8.79		Yes	164
8		1,069	8.78		Yes	133
9	Southern Cascadia, Gaussian	1,069	9.13		Yes	461
10		1,069	9.01	145, 242, 263	Yes	290
11		695	8.94		No	225
12		692	8.93		No	162
13	Gamma	1,069	9.11	100, 120, 150	Yes	342
14		1,069	9.10		Yes	229
15		1,069	9.10		Yes	119
16		173	8.25	60, 72, 90	No	342

Note. See Figure 5 for model number.

Figure 5

Gaussian and Gamma dynamic rupture simulations nucleated in northern or southern Cascadia (magenta star). In each plot, the colored rectangle corresponds to the along‐strike rupture extent. The solid or dashed lines correspond to the segment location, of which one is allowed to vary. (a) Gaussian ruptures where the width of the middle segment, containing the nucleation asperity, is varied until margin‐wide rupture no longer occurs. The T levels (relative dynamic stress drop) remains constant for each simulation. (b) Same idea as (a) but for Gamma rupture simulations. Moment‐magnitude is plotted along the x‐axis on all plots. (c) and (d) show Gaussian and Gamma ruptures nucleated in southern Cascadia, respectively. For each subfigure, the numbers refer to the model number, which are listed in Table 1.

We first study three segment rupture models that are nucleated in the northern Cascadia region to see how our choice of T and segment width affect final rupture length (Figure 5a). We find that placing segment limits near ∼46 and 43° latitude (Y ranges from 180 to −350 km; Figure 5a), together with T levels between 200–250 years (∼8–10‐m slip deficit), leads to margin‐wide rupture. In all segmented models, T levels are adjusted empirically based on comparisons of the model predicted subsidence amplitudes to those recorded from the 1700 A.D. event. We also note that the segment containing the hypocenter requires the highest T level for all segmented models. The position of these segment limits corresponds to changes in estimated RI level, tremor patterns and forearc morphology (Figure 1; Goldfinger et al., 2017; Watt & Brothers, 2020). The middle segment encompasses most of the creeping region offshore Oregon. Holding T levels constant, we move the location of the southern segment boundary northward by 0.5° latitude (∼56 km distance) until margin‐wide rupture is no longer observed. An average slip deficit of nearly 2 m over a width of ∼6 0 km (distance between dashed and dot‐dashed lines) is needed to drive rupture through the creeping section and into the southern end of Cascadia (Figure 5a). The higher coupling in the northernmost segment (Y > 200 km) allows for rupture to propagate north of the epicenter in all cases. In contrast, if we use the Gamma coupling model, margin‐wide rupture is much easier to attain even with lower relative stress drop (lower T values) (Figure 5b). Lower T levels are used in the Gamma rupture simulations as higher values are not required to achieve margin‐wide rupture with the Gamma distribution. This result demonstrates the sensitivity of margin‐wide rupture to the stress level in the shallow portions of the fault. As the length of the central segment becomes shorter, moment‐magnitude only weakly decreases (by ∼0.01) for Gamma ruptures. Gamma ruptures nucleated in northern Cascadia can feature shallow, narrow slip distributions and low rupture speeds ranging from ∼1 to 2 km/s in the central region of the megathrust (Figure S4).

For dynamic ruptures initiated in southern Cascadia, we repeat the T‐level experiment by fixing the segment limit near 43.5° latitude (Y = −312 km) and vary the northernmost segment limit. We found that slightly higher T levels (relative to ruptures nucleated in the north) are a necessary condition to sustain rupture propagation, particularly through the central Cascadia region (Figures 5c and 5d). Gaussian models that lead to a margin‐wide rupture required an additional slip deficit of 3 m over a length of ∼109 km in the central segment (i.e., Figure 5c, dashed line) compared to non‐margin‐wide rupture event (i.e., Figure 5c, dot‐dashed line). Similar to what was observed for ruptures initiated in northern Cascadia, Gamma coupling models tend to generate margin‐wide ruptures at much lower slip deficit (i.e., Figure 5d). We also present a Gamma model initiated in southern Cascadia that does not become margin‐wide (Figure 5d, gray dot‐dashed line) to illustrate that for the same segment locations but lower T, a rupture cannot penetrate through the central segment. The higher relative T levels in the southernmost segment allow rupture to initiate and propagate outside the region of spontaneous rupture initiation, given our slip‐weakening friction parameters (i.e., , , Dc) and effective normal stress that bound the fracture energy. The Cascadia megathrust dips more steeply below Oregon and this probably influences the initial stages of ruptures that propagate from south to north more than those that rupture north to south. In general, Gamma model results suggest that only a narrow region of concentrated, high dynamic stress drop is sufficient for promoting margin‐wide rupture, even if slip‐strengthening friction or higher cohesion levels are present. We will now discuss our assumptions about sediment friction in the shallow most portions of the megathrust, and its effect on rupture size.

Effect of Up‐Dip Frictional Behavior

In all simulations presented so far, we have assumed the influence of subducting sediments will lead to slip‐strengthening frictional behavior in the upper 5‐km of the megathrust along‐strike. We now relax this assumption and let the dynamic friction level vary from slip‐strengthening to slip‐weakening conditions (Figure 6) starting with the M8.80 Gaussian model shown in Figure 5a. In all simulations, we fix T levels and segment locations, while testing varying dynamic friction coefficients in the near‐margin region. Neither slip‐strengthening ( > 0.6) nor slip‐neutral (  =   = 0.6) friction leads to margin‐wide rupture for this particular parameterization (Figure 6a); only a slip‐weakening behavior at shallow depths allows rupture to spontaneously grow into a margin‐wide event. The effect of dynamic friction level on slip at the deformation front is shown in Figure 6b. We observe high slip amplitudes (>25 m) in northern and southern Cascadia and reduced slip in central Cascadia (Figure 6b). In the margin‐wide rupture case (i.e.,  = 0.1), this slip pattern is similar to other Gaussian coupling models.

Figure 6

Gaussian coupling models with variable sediment frictional behavior in the upper 5 km of the megathrust. (a) Along‐strike rupture lengths (colored lines) as function of dynamic friction coefficient. (b) Slip at the deformation‐front for each scenario shown in (a).

Effect of Down‐Dip Locking Depth

Estimating the seismogenic zone from the available geodetic data and paleoseismic measurements (Hyndman, 2013; Wang & Tréhu, 2016) is fraught with uncertainty because of their lack of offshore resolution. In both the Gaussian and Gamma coupling models, the down‐dip limit of coupling (positive stress drop) is near 20 km depth (Figure 1; Schmalzle et al., 2014), broadly consistent with thermal models proposed for this subduction zone (Cozzens & Spinelli, 2012; Hyndman, 2013; Wang et al., 1995). To assess how locking depth influences rupture width, length and subsidence amplitudes, we now relax this assumption and let locking depth vary. Note that the locking depth is meant as the maximum depth where slip‐weakening frictional behavior exists (i.e., < ). Again, we start with the three‐segment M8.80 Gaussian simulation (Figure 5a ), which does not break through the central Cascadia region (Figure 6a). Slip‐weakening behavior with  = 0.1 is initially set to end at 20 km depth and we systematically extend locking depth by two kilometers until 30 km (Figure 7a). A dynamic rupture simulation assuming a 15‐km locking depth is also shown for sake of comparison. We observe that moment magnitude increases (8.8 < Mw < 9.2) due to the propagation of rupture into the gap and slow‐slip regions of the fault. Ruptures progressively extend further south for greater locking depth, but do not become margin‐wide (Figure 7a).

Figure 7

Effect of down‐dip locking depth on coseismic slip distribution and uplift/subsidence patterns. (a) Final slip distributions for the Gaussian coupling model nucleated in northern Cascadia. The earthquake is nucleated at the red star and a profile of uplift/subsidence at the free surface is plotted in figure B (black line through red star). In each panel, the locking depth is systematically deepened by 2 km. (b) Model‐predicted coseismic subsidence and uplift for the range of locking depths studied. The coastline is plotted for reference. Each solid line represents the coseismic hinge‐line and is colored by its respective locking depth. We also show a shallower locking depth (15 km) in yellow for comparison.

We select a 2‐D profile near the hypocenter along‐strike to assess how the model predicts coseismic subsidence and amplitude patterns change in the margin‐perpendicular direction (Figure 7b). The maximum uplift and subsidence amplitudes increases by ∼1 m for every 2‐km increase in locking depth. For deeper locking depths, the coseismic hinge‐line moves closer to shore, although all hinge‐lines remain at least 100 km offshore for the profile selected in northern Cascadia (Figure 7b).

Fitting 1700 A.D. Subsidence Measurements

Previous elastic rupture models have shown that coastal subsidence measurements from 1700 A.D. can be well fit with high slip‐patches positioned along‐strike. Wang et al. (2013) used static models with four distinct asperities with T levels ranging from 450–550 years (18–22 m slip deficit) to reproduce the subsidence amplitudes. In these static models, the greatest locking depth was taken to coincide near the 350 C isotherm as this is where silica‐rich lithologies would be expected to transition from velocity‐weakening to velocity‐strengthening frictional behavior (Wang et al., 2003). 3‐D kinematic simulations used a range of locking depths (∼10–30 km) and determined that, in the presence of subevents, a locking depth near ∼15 km provided the strongest fit to the subsidence data (Wirth & Frankel, 2019).

We present four non‐unique 3‐D dynamic rupture scenarios derived from Gaussian and Gamma coupling distributions with a shallower locking depth at 15 km, but we also test a deeper (∼20 km) locking depth given in the Schmalzle et al. (2014) coupling models (see Discussion). The T levels and segment locations were selected through a trial‐and‐error approach. These dynamic source models show 1700 A.D. subsidence data can be reasonably fit without invoking high amplitude slip deficit or subevents (Figure 8). Ruptures initiated in northern Cascadia with modest T levels (250 years) can match subsidence data with three segments (Figure 8a) whereas we find that four segments are required for the Gaussian rupture model initiated in southern Cascadia (Figure 8b). The Gaussian‐type simulation initiated in southern Cascadia has a final rupture length ∼100 km shorter than the other ruptures. We note that the Gaussian simulation nucleated in northern Cascadia provides the best fit to the 1700 A.D. subsidence data (Figure 8b, blue line). There are probably several combinations of T levels and segment locations that can fit to the 1700 A.D. subsidence data, and moreover, the four exemplary models we show here are meant to show how shallower coupling distributions can affect the rupture dynamics.

Figure 8

Gaussian and Gamma coupling models with shallow locking depth (15 km) that match coastal 1700 A.D. subsidence measurements to first order. (a) Ruptures nucleated in northern Cascadia. (b) Ruptures nucleated in southern Cascadia. (c) Comparison of predicted coastal subsidence from simulations shown in (a) and (b). The white dashed line in (a) and (b) represents the 20‐km depth contour.

Discussion

What Allows Large Earthquakes to Develop Along the Cascadia Megathrust?

We observe margin‐wide ruptures under conditions of higher relative dynamic stress drop amplitudes in the inferred creeping region of the central Cascadia megathrust segment (i.e., higher T levels relative to the other segments). Alternatively, we also show that margin‐wide ruptures are promoted by slip weakening frictional behavior at shallow depth (Figure 6).

When the Gaussian or Gamma coupling models are used to generate heterogeneous shear stress distributions, there are two natural locations to initiate spontaneous rupture: in northern or southern Cascadia. Our results suggest that if rupture initiates in southern Cascadia, higher T levels are required to sustain rupture through the central creeping region for Gaussian stress distributions (Figure 5c). This is due to the combined effects from a lower slip‐rate deficit (inherent to both geodetic coupling models) and the generally narrower seismogenic rupture area offshore Oregon caused by an increasing megathrust dip angle in this region.

A notable feature of our dynamic rupture simulations is that large earthquakes (M8.8 and above) can be generated at much lower T levels than previously suggested from static models (i.e., Wang et al., 2013). An explanation for this comes from dynamic effects within the wedge. For instance, even though Gaussian simulations have little to no slip deficit extending to the deformation front, reflections within the accretionary wedge appear to drive rupture propagation along‐strike. While a more realistic rheology within the wedge would certainly affect wave propagation, our models suggest that wavefield interference at shallow depths could be a viable mechanism to sustain rupture (Huang et al., 2014). Velocity‐strengthening , or in our case, slip‐strengthening, friction is a common assumption in dynamic rupture simulations for subduction zones to represent the frictional behavior of sediments near the trench (Kozdon & Dunham, 2013). One may also explicitly incorporate a subducting sediment channel structure with depth‐varying rigidity using slip‐weakening friction (i.e., Ulrich et al., 2020). The presence of clays or fluids within the megathrust fault zone can complicate the frictional behavior, however (Saffer & Tobin, 2011). While Cascadia is well‐known to have significant sediment blanketing the trench along‐strike with variable states of consolidation (Han et al., 2017), there are no studies that directly sampled Cascadia megathrust fault gouge and subject them to high‐velocity friction experiments (Seyler et al., 2020). The assumption of slip‐strengthening friction in the upper 5 km in our dynamic rupture simulations is therefore modest and will greatly benefit from offshore drilling data. We do not repeat the exercise of lowering the dynamic friction level for ruptures in southern Cascadia, but we expect that a similar behavior would occur.

The bulk of this study explored T levels that are below 320 years and exhibit segmentation along‐strike. It is geologically reasonable to presume the convergence rates of the Juan de Fuca and Gorda plates have been stable for ∼0.78 Ma (DeMets et al., 2010). When megathrusts ruptures occur, they can relieve the accumulated slip deficit partially or completely (Hardebeck & Okada, 2018). Given this, our results could be interpreted to mean that the slip‐rate deficit is not accumulated the same everywhere along the Cascadia megathrust. This is consistent with suggested RI segmentation (e.g., Goldfinger et al., 2017), some of which may be expressed as along‐strike changes in fault vergence patterns and outer wedge geometry (Gulick et al., 1998; Watt & Brothers, 2020). There is also the possibility that the coupling models we used in our study overestimated the slip‐rate deficit due to inelastic off‐fault deformation (e.g., Baker et al., 2013). On the other hand, if the slip‐rate deficit from plate convergence is being accommodated along Cascadia at a constant rate, then in order to achieve reasonable coastal subsidence amplitudes (>−2 m) our models imply a partial stress drop during the 1700 A.D. event. Lower T values set a critical initial condition of available stress drop in our rupture simulations. In particular, the Gamma rupture models seem to require lower T amplitudes to generate margin‐wide rupture than RI values estimated from turbidite paleoseismology (i.e., Goldfinger et al., 2012; Figures 5b and 5d) whereas Gaussian rupture models initiated in southern Cascadia require more relative stress drop than those initiated in northern Cascadia (Figures 5a and 5c). But because there are no geophysical recordings of the 1700 A.D event, a detailed slip inversion will be remain out of reach. However, studies have inferred potential far‐field tsunami patterns of this earthquake using simple kinematic rupture models and Japanese written records (e.g., Satake et al., 2003, 2020). Seafloor deformation predicted from dynamic rupture models may be readily incorporated into kinematic tsunami propagation simulations ‐ differences in nucleation or peak slip locations (Figure 8) may alter our understanding of historic tsunamigenic earthquakes in Cascadia.

Explaining 1700 A.D. Subsidence Patterns With Dynamic Rupture Simulations

The geodetic coupling models we use show positive stress drop down to ∼20 km depth (Figure 1). On the other hand, 3‐D kinematic simulations were able to match 1700 A.D. subsidence data assuming positive stress drop does not extend deeper than a fixed coupling level closer to 15 km depth (i.e., 1 cm/yr contour from Burgette et al., 2009; McCaffrey et al., 2013; Wirth & Frankel, 2019). We show that dynamic rupture simulations that taper stress drop to 0 MPa below 16 km depth can agree well with the 1700 A.D subsidence data, particularly for ruptures initiated in northern Cascadia (Figure 8). While shallower locking depths generally provide a stronger fit to the paleoseismic data, we were also able to construct a dynamic rupture simulation with a 20‐km locking depth that fits the data just as well (Figure S5). This result suggests that if T levels are sufficiently low along the fault, the subsidence patterns can probably be fit by an even deeper locking depth (>20 km depth). The influence of a deeper locking depth is to move the coseismic hinge‐line landwards and increase the amplitude of subsidence and uplift (i.e., Figure 7b). As discussed in Kanda and Simons (2012), either the location of peak interseismic uplift‐rate or greatest coseismic subsidence can provide a stronger constraint on the extent of slip as opposed to the hinge‐line. Unfortunately, both the interseismic uplift (i.e., Krogstad et al., 2016) and paleoseismic subsidence data (i.e., Wang et al., 2013) are limited in the along‐dip direction for this subduction zone. We thus caution using only paleoseismic subsidence data to uniquely constrain the down‐dip rupture limit in Cascadia.

The Potential of Heterogeneous Down‐Dip Frictional Properties

The next Cascadia megathrust rupture may or may not include high‐frequency seismic energy radiated near the down‐dip limit of slip (e.g., Lay, 2015), but one way to accomplish this is to superimpose high stress drop subevents (>15 MPa) at several locations along‐strike (Figure 9). For these simulations, we assume a locking depth of 15 km. The influence of subevents, compared to a coupling model with no subevents, is to increase the subsidence amplitude and generate higher relative seismic frequencies.

Figure 9

Comparison between dynamic rupture simulations with and without high dynamic stress drop subevents positioned near the down‐dip edge of locking (∼15 km) along‐strike. (a) Gaussian dynamic stress drop distribution without subevents. White triangle denotes synthetic seismogram receiver location. (b) Gaussian dynamic stress drop distribution with superposed ∼15–20 MPa subevents (white boxes). (c) Coastal subsidence for both models with 1700 (a)D. observations for comparison. (d), (f), and (e) show the normalized spectral amplitudes of the x‐, y‐, and z‐component velocity seismograms, respectively. The influence of the subevents is to increase the high‐frequency amplitudes >0.1 Hz (bold colored lines).

To conceptually demonstrate that heterogeneous D can also generate relatively higher seismic frequencies in the specific case of the Cascadia margin, we also design a dynamic rupture simulation containing several 16 km2 asperities near the down‐dip edge of the locked megathrust that have lower D c = 1 m with D c = 2 m everywhere else (Figure S6). These D c levels are chosen to be consistent with already presented 3‐D rupture simulations that can resolve the cohesive zone widths. In this particular model, the effect of heterogeneous D is to slightly increase the high‐frequency energy content >0.1 Hz, with stations further away from the hypocenter showing this effect more clearly (Figure S6). The model with a heterogeneous D distribution is 0.28 magnitude units above the model with constant D .

What is unclear are properties most conducive to generating high frequency seismic radiation down‐dip, but dynamic rupture simulations for the 2011 Tohoku‐Oki earthquake showed that heterogeneous frictional properties or strength distributions (Galvez et al., 2014; Huang et al., 2014) might account for these observations. We note that Cascadia is remarkably different from the Japanese or Chilean subduction zones. In particular, the subducting interface of the Juan de Fuca plate is relatively smooth along most of the margin compared to in the aforementioned regions (van Rijsingen et al., 2018) and consequently, the interface topography may not provide an obvious explanation for future high seismic frequencies radiated down‐dip.

Limitations and Future Directions

Our study incorporates a physically consistent source model that emphasizes the importance of frictional and stress conditions necessary to generate M9‐type ruptures. For lack of detailed information on the velocity structure in the accretionary prism and the highly simplified 1‐D CVM used, our 3‐D dynamic rupture simulations do not capture accurate wave propagation effects along the Cascadia margin. Forecasting accurate ground motions during megathrust earthquakes is important, especially for subduction zones with limited or no seismic recordings (Frankel et al., 2018; Wirth et al., 2018). Developing dynamic rupture simulations that account for 3‐D source, site, and path effects is one future direction that would, for instance, lead to more physically informed hazard estimates (e.g., Wirth et al., 2020).

Another limitation of our study is the rheology assumed: a linearly elastic body. There is potential for off‐fault inelastic deformation in the wedge where there is a significant sediment volume (Ma, 2012). We note further that our choice of a linear slip‐weakening friction law allows us to assess first order along‐strike and along‐dip rupture limits, similar to Ramos and Huang (2019). Modifying the friction law (and adjusting the finite element mesh resolution accordingly) to account for strong rate‐weakening would permit us to test a wider range of rupture styles. Understanding what fault zone lithologies are present along the Cascadia megathrust would also be helpful in assigning realistic dynamic friction levels during coseismic rupture.

To improve the predictive capability of dynamic rupture simulations, offshore (e.g., near‐trench) geodetic and seismic measurements are needed. It would be particularly valuable if information about the interseismic uplift‐rate could be constrained offshore, to extend existing onshore leveling data (Krogstad et al., 2016). This would reduce the ambiguity in geodetic coupling models and improve our understanding of the spatial relationships between upper plate deformation and intra‐plate slip behaviors (Bruhat & Segall, 2017; Li & Liu, 2021; Malatesta et al., 2021; Watt & Brothers, 2020). Our study stresses the importance of the spatial variation in coupling, especially in the central Cascadia region where confirming the presence of lower coupling offshore Oregon is critical for both kinematic and dynamic rupture simulation results. Recent kinematic inversions for coupling suggest that high coupling may exist across Cascadia (i.e., Lindsey et al., 2021), which will have implications for future coseismic rupture propagation. A dense array of ocean‐bottom seismometers (for both active and passive seismic source studies) would provide further constraints on shallow megathrust geometry as well as better characterize interplate seismicity patterns offshore. A wider range of earthquake and tsunami hazard scenarios can be designed in the future as these data become available.

Other geophysical measurements that have not been incorporated in this suite of dynamic rupture simulations include inferences made about the seismogenic width from the arguably highest resolution geophysical data set available: the free‐air gravity field. Basset and Watts (2015) observed that trench‐parallel ridges in the free‐air gravity anomaly field correlate well to the top of slow‐slip and tremor across the Cascadia forearc. If such trench‐parallel features in the gravity field are a proxy for down‐dip rupture extent, then the transition from slip‐weakening to slip‐strengthening frictional behavior may extend to depths greater than 20 km in some regions of Cascadia. More work is needed to identify what geologic or geophysical features are most indicative of future coseismic rupture limits, especially in subduction zones like Cascadia that have not experienced megathrust events during the modern era of instrumentation.

Conclusions

Developing realistic seismic source models for the Cascadia megathrust is of paramount importance to assist with seismic and tsunami risk mitigation. We present 3‐D dynamic rupture simulations that incorporate different hypotheses for megathrust strain accumulation based on geodetic coupling models. We show that in order for margin‐wide, “M9” type ruptures to develop, there must be a sufficiently high relative dynamic stress drop in the central Cascadia region. Moreover, a slip weakening behavior or moderate slip deficit close to the deformation‐front can greatly facilitate margin‐wide ruptures. Along‐strike variations in the slip deficit pattern relative to the geodetic coupling models were required to match available paleoseismic data in our simulations. We note that strain accumulation times lower than those suggested from paleoseismic studies provide a better fit to the subsidence data, which might suggest coupling models are overpredicting the slip‐rate deficit or there was incomplete stress drop from the last megathrust rupture. A close fit to 1700 A.D. subsidence data can be achieved using Gaussian or Gamma coupling distributions with locking depths of 15 or 20‐km depth, obviating the need to call upon localized, high amplitude slip asperities along the down‐dip region of the seismogenic zone.

This work is a step forward in using fully dynamic rupture simulations for seismic hazard analysis where there have been no instrumentally recorded ruptures. Kinematic rupture properties (e.g., rise‐time, slip‐rate and rupture speed) and static seafloor displacement from our dynamic simulations can be readily incorporated into existing 3‐D kinematic rupture simulations or inform tsunami propagation and coastal inundation models.

Supporting Information S1

Click here for additional data file.