Geodetic Coupling Models as Constraints on Stochastic Earthquake Ruptures: An Example Application to PTHA in Cascadia

Current stochastic rupture modeling techniques do not consider the potential influence of inter‐seismic coupling, a first‐order property of a megathrust, which can show correlation between areas of high coupling and areas of greater slip as seen in recent large ruptures globally. Therefore, it is reasonable to assume that it should be considered as prior information in rupture modeling. Here, we first present a mathematical formalism to introduce coupling models as prior information into stochastic rupture modeling. We then focus on how introducing slip deficit information into the stochastic rupture models influences slip distributions for the Cascadia subduction zone (CSZ). We compare rupture models created with two end member models of coupling, one with a shallow coupling and another with coupling deeper downdip. We also discuss the comparison to models created without assuming knowledge of the coupling distribution except for variation in the downdip limit of slip. Variations occur and correlate well with areas with the largest differences in slip deficit rates. The ruptures are then used for regional probabilistic tsunami hazard assessment. Overall, the tsunami amplitudes generated are much more hazardous in the northern extent of the CSZ where differences in coupling distribution are more prevalent. Models obtained from assuming a shallower downdip limit have tsunami amplitudes more similar to those from the geodetic coupling models. Although uncertainties are present in the accuracy of coupling, imposing either constraint created different hazard estimations when compared to those where no prior coupling information was used.

stochastic sources can even be used for analyzing the response of structures and critical infrastructure to potential future earthquakes (Bijelić et al., 2018). Likewise, stochastic sources and their resulting seismic and geodetic waveforms can be used to test the response of earthquake early warning systems (Goldberg & Melgar, 2020;McGuire et al., 2021;Ruhl et al., 2017) and tsunami early warning systems that rely on onshore data (Williamson, Melgar, Crowell, et al., 2020). This same approach is now becoming common practice in probabilistic tsunami hazard analysis (PTHA) where stochastic source modeling, when connected with advanced numerical models that simulate tsunami propagation, are also rapidly becoming a mainstay of that field (Grezio, Babeyko, et al., 2017).
Like the advances in stochastic source modeling, tectonic geodesy has progressed substantially in the last 30 years, spurred from the expansion of Global Navigation Satellite System (GNSS) capabilities. For example, it has become commonplace to use measurements of the long-term inter-seismic velocity field to identify which portions of faults are more or less coupled. A review of this can be found in Bock and Melgar (2016). Although our ability to predict future earthquakes is limited, we are able to relate these coupling models to long-term future earthquake potentials. Notably, in the California earthquake hazard model UCERF3 (Field et al., 2017) the inter-seismic surface velocity field is used to constrain long term "slip rates" of faults which are a primary constraint to determine the moment budget over broad areas along a fault zone. Overall, it is generally agreed that where coupling is higher and faults are accumulating a slip deficit at a faster rate, earthquakes are more likely to occur. Similarly, the areas with the greatest amount of slip during large ruptures often correlate with highly coupled patches of the fault zone. This relationship between coupling and slip is seen in previous large ruptures across a range of fault zones (e.g., Barnhart et al., 2016;Konca et al., 2008;S. Li & Freymueller, 2018;Moreno, Rosenau, & Oncken, 2010;Ozawa et al., 2011).
This correlation is not perfect, however, and the present-day pattern of heterogeneous coupling on a fault does not equate to the pattern of slip during the next earthquake. There is a myriad of other controlling variables which can be unique to each fault zone, such as the complex fault geometries, past and present stress regimes, and rupture dynamics to name a few. This creates large uncertainty in how a future earthquake will rupture. In spite of these, it remains true that highly coupled fault patches are accumulating a slip deficit at a faster rate and thus have a larger budget of available slip to use during the next event. Here we will show how to use that fault coupling to condition the likelihood of where slip should be expected in future earthquakes. Our proposed approach does not require, in a deterministic sense, that high slip occur in highly coupled patches, rather it increases the probability that it does. As a result, over many stochastic rupture simulations, on average, more slip will occur where inter-seismic coupling is higher and less where coupling is lower. However, for any one particular realization, slip can still be high in a low coupling region, and low in a highly coupled region.
To illustrate the impact of this we will show how assuming different coupling models impacts PTHA in the Cascadia subduction zone (CSZ) in particular. Subsequent to the 2004 M9.2 Sumatra earthquake and tsunami which led to 240,000 casualties, probabilistic tsunami hazard assessment (PTHA, Geist & Parsons, 2006;Grezio, Babeyko, et al., 2017) has become a rapidly expanding methodology used for assessing the hazard potential of future earthquakes and tsunamis. Unlike site-specific tsunami studies, PTHA is rooted in determination of the probability of exceeding some threshold of tsunami intensity measure (e.g., tsunami arrival height) for one or many target sites for a given return period (e.g., 100 years). A fundamental difficulty addressed by application of PTHA is that awareness of previous historical tsunamis past hundreds to thousands of years is limited and sometimes based on historical recounts rather than surficial expressions or direct measurements of their impacts. In recent advances to PTHA, it has become more prevalent to rely on modeling complex earthquake source rupture processes, such as the heterogeneous slip distributions, and use the resulting deformation as the initial condition for propagation modeling. This is then combined with some probabilistic scaling relations as well as earthquake recurrence rates for quantification of likely exceedance of tsunami hazard intensities (e.g., De Risis & Goda, 2017;Grezio, Cinti, et al., 2020;L. Li et al., 2016).
In this approach to PTHA, arguably, the largest source of uncertainty is the earthquake or rather the tsunami source, since propagation modeling is already highly advanced and bathymetry, both in the deep ocean and in the nearshore, is relatively well-known. Earlier tsunami forecasting methods assumed a uniform slip on a geometrically simple fault buried in a homogeneous elastic half-space. It is now broadly recognized SMALL AND MELGAR 10.1029/2020JB021149 2 of 20 that assuming homogeneous slip drastically underrepresents the tsunamigenic potentials when compared to heterogeneous slip models for the same magnitude (Melgar, Williamson, & Salazar-Monroy, 2019;Ruiz et al., 2016). The situation is different, however, for inundation modeling of overland flow. Many advances in modeling the fluid dynamics of propagation over an erodible substrate and through the built environment are still necessary and introduce equally, if not larger, sources of uncertainty than knowledge of the earthquake source. However, for the simpler problem of quantifying the expected tsunami at the shoreline, without considering inundation, better constraints on what sources can realistically be expected to occur is one of the improvements that can most significantly reduce uncertainties in the hazard estimate. Lastly, presence and contribution of splay faulting on tsunami hazards can also produce reasonable further uncertainties for PTHA studies (e.g., Gao et al., 2018). For the purpose of the paper, however, we do not focus on this aspect.
In this work we focus on the CSZ as it is extensively studied, has a suite of previously constrained coupling models, is well instrumented onshore, and is perceived as having high associated hazards. First, we present the mathematical formalism for introducing a geodetic coupling model as a prior constraint on the resulting stochastic rupture. We then model a total of 11,200 stochastic slip rupture scenarios for a magnitude range of M7.8-M9.1 on the CSZ in order sample a large enough range of tsunamigenic ruptures. We make four different assumptions on prior coupling which leads to for four distinct classes of ruptures. Two classes include different end member coupling models from Schmalzle et al. (2014) that vary in terms of their near-trench coupling estimates (noted as the "Gamma" and "Gaussian" coupling models). For the other two classes of ruptures, we produce stochastic slip in the traditional way, where, beyond assuming a down-dip limit of slip, there are no assumptions on where slip can occur, that is, with no coupling model. Following Frankel, Chen, et al. (2015), these last two classes have a deep ∼30 km depth limit associated with the top of the tremor zone, and a ∼15-25 km depth limit associated with the average 1 cm/yr slip deficit rate (25%) contour. For convenience, the class with a downdip limit associated with the top of the tremor zone will hereafter be called the "30 km depth" class. This last class follows what was proposed by Wirth and Frankel (2019) and although it is defined from regional coupling, no knowledge of the pattern of coupling is associated with the calculation of these ruptures.
For each of the rupture simulations we model tsunami propagation to the coast and analyze the resulting coastal tsunami amplitudes using hazard curves and hazard maps. Tsunami intensities vary considerably between classes, where these differences are best expressed off the coast of the Olympic Peninsula. Hazard curves for coastal points and hazard maps detail that with increasing distance north, the variations in hazards become much more distinct, where the coastal points experience greatest tsunami arrivals for the Gamma class. We stress that in this work we are not making an authoritative hazard assessment for the CSZ. Rather, we simply aim to demonstrate the impacts of conditioning the ruptures with a coupling model. Whether or not ruptures constrained by coupling are similar to past and future ruptures is a key question moving forward. It is our hope, however, that future authoritative hazards assessments more formally consider geodetic coupling and that our work be used as motivation to better constrain offshore coupling in particular through expanded seafloor measurements. The variability between models that we will show highlights the need for these improvements.

Geodetic Coupling Models
Fault coupling models (e.g., Bürgmann et al., 2005;Konca et al., 2008;Loveless & Meade, 2010;Moreno, Melnick, et al., 2011) are obtained from inversion of surface velocities measured by geodetic techniques such as global navigation satellite systems (GNSS) and tide gauges. Accurately constraining coupling models can be done through integrating paleoseismic data (e.g., McCaffrey et al., 2007), however, good paleoseismic records are typically sparse for a given region. In theory, coupling models provide an estimate of a fault's present-day ability to move. Coupling is typically quantified by the ratio of the on-fault slip rate versus the long term far-field plate rate (e.g., the convergence rate at a subduction zone). A ratio of 1 is fully coupled and represents no current motion along the plate interface during the inter-seismic period, whereas a ratio of 0 indicates aseismic stable sliding (creep). Conversely, coupling models also include by definition an estimate of the "slip deficit rate." The coupling fraction is therefore the ratio between the slip deficit rate and the local plate convergence rate. Areas that are fully coupled, are accumulating a slip deficit fastest, specifically at the plate convergence rate.
The reliability of the calculated coupling models is dependent on the abundance, distribution, and timescale of geodetic stations (primarily GNSS). The resolving power of a specific GNSS site falls rapidly with distance. Currently, GNSS coverage of subduction zones worldwide is quite good for the on-shore portion (e.g., Barrientos & Pérez-Campos, 2018) but is very limited offshore. Seafloor geodetic instrumentation, which was pioneered at the CSZ (Spiess et al., 1998) is currently only widely implemented in Japan (Yokota et al., 2016), however its use is slowly expanding. This lapse in GNSS coverage means that the uncertainty in the recovered interseismic coupling, which is largely in the offshore shallow portion of the subduction zone, can be quite high (e.g., Loveless & Meade, 2010;Schmalzle et al., 2014). Due to the sparsity of station coverage, a multitude of non-unique estimations of plate coupling can be determined for a region. As a result of this uncertainty, different modeling assumptions will lead to different results for the near-trench coupling. Although coupling models are non-unique, they currently provide one of the best approaches for understanding the influence of regional slip deficits on rupture heterogeneities in a given area.
In this paper, we consider two different coupling models ( Figure 1) for the CSZ from Schmalzle et al. (2014). Both represent the "decade scale" distribution of coupling where the effect of transient slip such as slow slip events is accounted for. The first model uses an a priori assumption that complete coupling occurs from the trench to some distance down dip and that at all points along strike coupling fraction decreases to free slip by a shape factor, gamma (K. Wang et al., 2003; further referred to as the "Gamma" coupling model). In contrast, the second approach uses a Gaussian distribution of coupling with depth, dependent on a model parameter, mean depth, and spread, or standard deviation, of the coupling (further referred to as the "Gaussian" coupling model). Due to the Gaussian nature of the coupling pattern, at shallow depths near the trench, the plate has low coupling and is mostly stably sliding. In this model coupling is centralized down dip ( Figure 1). These models were computed on an older megathrust geometry from McCrory et al. (2012). However, we translated these models to the updated slab geometry from Hayes et al. (2018) through nearest neighbor interpolation. It is important to note that both of these approaches fit the known on-shore GNSS velocities, tide gauge records, and geologically derived uplift rates with near identical confidence.

Stochastic Ruptures With a Geodetic Coupling Model Constraint
The first step in obtaining a rupture model is to determine which portions of the larger megathrust will contribute to a given earthquake. Since our target magnitudes span a range of M7.8 to M9.1 not all portions of the megathrust will participate in any given rupture. To select a subset of the megathrust we follow the approach detailed in Melgar, LeVeque, et al. (2016). Given a target magnitude, we determine the length, L, and width, W, of the rupture. We make a random draw from probabilistic scaling laws (Blaser et al., 2010) which state that length and width follow the magnitude dependent log-normal distribution with mean and standard deviation given by with standard deviations defined in the original work. By making these random draws, the objective is to obtain a length and width that is consistent with the behavior seen in earthquakes worldwide while retaining the observed variability as well. The probabilistic scaling law thus ensures that for a given magnitude we do not always employ the same fault dimensions. Once these dimensions are known we randomly select a portion of the larger megathrust models that is within these bounds. An example of this can be seen in Figures 2 and 3 where the fault dimensions for a potential M8.7 event have been obtained from the probabilistic scaling laws and used to define a subset of the megathrust to use for the stochastic rupture.
Next we define the rupture model itself. Slip can be conceptualized as a spatially random field whose heterogeneity can be described by statistical parameters. A number of unique slip realizations can then be determined as long as they are constrained by an underlying probability distribution. Mai and Beroza (2002) SMALL AND MELGAR 10.1029/2020JB021149 found that to best model the spatially random slip field the most suitable autocorrelation function (ACF) is the von Karman ACF. In their proposed approach, the von Karman ACF is enforced using a spectral representation, P(k), of slip in the Fourier domain defined as the ratio of the correlation lengths for along-strike where H is the hurst exponent describing the spectral decay at higher wavenumbers. The radial wavenumber, k, is then defined as The Gaussian decade-scale coupling model which imposes a Gaussian distribution of coupling with depth as well as penalties to constrain mean coupling above 30 km in depth. (c) Downdip limit of slip model defined by the Frankel, Chen, et al. (2015) 1 cm/yr coupling contour (25% coupling) limit of slip. Color defines either active regions that may participate in slip (above limit) or regions of the fault that will not (below limit). Contours are the slab depths from Hayes et al. (2018) at 5 km intervals. Inset at top right shows the location of the CSZ.
The correlation lengths described in the von Karman ACF determine the dominant asperity sizes for the model. Here, it is determined that as magnitudes increase, the correlation lengths increase as well following a log-linear scaling (e.g., Melgar & Hayes, 2019). The Hurst exponent in Equation 3 on the other hand, seems to be magnitude independent and typically is assumed to be between 0.4 and 0.7 (Mai & Beroza, 2002;Melgar & Hayes, 2019).
The approach from Mai and Beroza (2002) is well suited for an approximation of a rectangular fault geometry, however, some complexities such as a multi fault or 3D fault geometry (e.g., the large bend in northern CSZ) can be difficult to account for. Similarly, enforcing prior constraints on the rupture model, such as the geodetic coupling model, is not inherently straightforward. An alternate approach is to work directly in the spatial domain. LeVeque, Waagan, et al. (2016) use the Karhunen-Loeve (K-L) expansion with the von Karman ACF to the same effect as the spectral approach. The spatial representation of the VK-ACF is used, where C ij is the correlation between the i-th and j-th subfaults, K H is the modified Bessel function of the second kind and H is the Hurst exponent. r ij is the inter-subfault distance given by where r s is the along-strike distance and r d the along-dip distance. These are obtained using a spline interpolation of the 3D fault geometry as detailed by Melgar, LeVeque, et al. (2016). Once all the parameters of the correlation matrix are defined, the covariance matrix is obtained by where  is the standard deviation of slip, which we set here to 0.45, irrespective of magnitude, following Melgar and Hayes (2019). The K-L expansion then states that to obtain a random realization the slip vector, s, that contains each subfault's slip will be Here, μ, is the mean of s and the statistics of the VK-ACF are enforced by the eigenvalues,  k , and eigenvectors,  k, of C ij (Equation 8). k z are normally distributed random numbers with a mean of 0 and a standard deviation of 1 which introduce the desired stochastic variability. N is the number of eigenmodes which corresponds to the number of subfaults or elements of s. If all the eigenmodes are used, then the stochastic realization produces the same results as if the analysis was carried out in the wavenumber domain. For certain applications, such as PTHA, where long period features dominate the resulting tsunami signals, after the first few tens of modes the contributions to tsunamigenesis from the short length-scale modes becomes negligible (LeVeque, Waagan, et al., 2016). In these cases, the summation can be truncated. Here, however, we use all eigenmodes. In the absence of any external knowledge it is commonplace to assume a homogeneous mean slip model, . In this case, given the assumed fault dimension and rigidities from the reference Earth model, enough slip is distributed over all subfaults to match the desired target magnitude. In other words, slip is equally likely at all subfaults irrespective of location both along strike and down dip of the fault. Once more, if this is done then the results will be the same as if carrying out the stochastic slip SMALL AND MELGAR 10.1029/2020JB021149 7 of 20 realization in the wavenumber domain. However, here lies the critical advantage afforded by the K-L expansion, we can assume that the mean slip model, , is related to, or rather defined by, the geodetic coupling.
The process for defining the mean model, , based on the geodetic coupling is as follows and is shown in Figure 2. Once a megathrust segment is selected following the probabilistic scaling laws for the desired magnitude (M8.7 in the example in Figures 2 and 3) the geodetic coupling (Figures 2a and 3a) is re-scaled to slip (Figures 2b and 3b) to match the target moment. This mean slip model now has the same features as the coupling, with higher slip in high coupling areas and lower slip where coupling is low, and even 0 slip where coupling is 0. Figures 2c and 3c then show 8 realizations of stochastic rupture using the K-L expansion with the non-homogeneous mean model. Note that each realization does not look exactly like the coupling model, but on average slip is more likely where coupling is high, as desired.
To carry out these stochastic realizations constrained by a non-homogeneous mean model we modified the open access forward rupture modeling FakeQuakes module which is part of the forward modeling and inversion code MudPy (Melgar, 2020;Melgar & Bock, 2015;Melgar, LeVeque, et al., 2016). Two-hundred stochastic slip ruptures for 0.1 magnitude bins ranging between M7.8 and M9.1 are modeled for the Gamma, Gaussian, 30 km depth, and 1 cm/yr class scenarios ( Figure 1). Following the K-L expansion approach, slip on triangular subfaults are determined using all eigenmodes. Ruptures are fixed to the given desired magnitude and have a mean rake of 90 in order to account for pure thrust motion. Following Graves and Pitarka (2015), we introduce some stochastic variability around the rake value as well. The hypocentral location is also randomly assigned from the selected subfaults. The location of the slip area and the hypocentral location are unconstrained, allowing ruptures to be equally likely anywhere along-strike on the megathrust. Additionally, the degree of slip on a single subfault was given an upper bound of 100 m of slip. Given the two presented coupling models, accumulation of 100 m of slip deficit would take well over 300 years (the time since the last large rupture), however, we consider this to allow for ruptures within a given magnitude bin to vary quite drastically between scenarios. For the rupture classes that follow more traditional stochastic slip modeling, no tapering of slip downdip is applied. Figures 2 and 3 shows one example of the resulting rupture models calculated using the Gamma and Gaussian model implemented in the stochastic rupture modeling. Although the coupling model is applied as the desired mean slip, variations in rupture slip patterns and resultant displacements are present and are further discussed later.

Tsunami Modeling and Hazard Curves
For each rupture coseismic vertical displacements at the surface are determined using the analytical solution for angular dislocations for triangular subfaults in an elastic half space (Comninou & Dundurs, 1975). This method is an adaptation from the Okada equations (Okada, 1985), which focus on rectangular subfaults. The resulting vertical deformation patterns for 2 scenarios is depicted in Figure 4. This calculated deformation is then used as the initial condition for tsunami modeling. Here we use the finite volume 2D depth-averaged, non-linear tsunami modeling code GeoClaw (http://www.clawpack.org/geoclaw) (LeVeque, George, et al., 2011). Since rupture propagation velocities are much faster than tsunami wave velocities, we assume instantaneous rupture as the initial condition for the system of partial differential equations. This assumption has a negligible effect on near-source modeling as discussed in Williamson, Melgar, and Rim (2019). Topography and bathymetry from the SRTM15 relief model sampled at 15 arcseconds (Tozer et al., 2019) are used in the tsunami modeling. For each scenario we modeled wave propagation for 4 h after rupture initiation, as we are only interested in the variability of regional arrivals. One of the powerful functions of GeoClaw is the adaptive mesh refinement (AMR) which makes the simulations efficient so that more intricate and complex tsunami characteristics are represented by the finest bathymetric resolution and more simple waves are adaptively coarser. We used four levels of refinement with the coarsest level at 3 arcminutes and the finest level at the 15 s of the topography/bathymetry data. we use Green's law to re-scale the wave amplitudes to 1 m depth. Example tsunami models are shown in Figure 5 for the two rupture scenarios detailed prior.
Hazard curves and resulting hazard maps are calculated and reflect the probability of exceedance of tsunami arrival amplitudes for each coastal point over a given return interval. Inherently in the formulation of the hazard curves is the assumption of the time-dependency of earthquake occurrence. We assume a magnitude-time dependent relationship defined by the Gutenberg-Richter (G-R) distribution, where N is the number of yearly events for a given magnitude. The constants a and b are assumed to have values of 6.279 and 1, respectively. These values are used because they produce a return period of 526 years for M9 earthquakes. This aligns reasonably well with what is expected for the return period of M9 events for the CSZ from the paleoseismic record (e.g., Goldfinger, Galer, et al., 2017;Goldfinger, Nelson, et al., 2012). It should be noted that the applicability of the G-R distribution breaks down at lower magnitudes for the CSZ, however, since we are not focused on an authoritative hazard assessment for the region, we still employ it for the study. For results to be used for authoritative hazard assessment, other magnitude-frequency distributions should also be considered, for instance, the tapered GR distribution (Rong et al., 2014). Another potential approach is to use expert opinion to build logic trees with suitable rates for different magnitudes (e.g., Frankel, Chen, et al., 2015).  After assuming the rates at which earthquakes occur for every coastal location, we compute      c P , the probability of the tsunami amplitude at the coast, , exceeding a given threshold,  c. Following Geist and Parsons (2006), Here the product operates over K magnitude bins. In our case, there are fourteen bins for M i = [7.8,9.1] with 0.1 magnitude units between bins. t is the chosen return interval of interest and N i is the rate at which earthquakes in a given magnitude bin are assumed to occur from the G-R distribution.

Results
Here we present the results of the 11,200 rupture scenarios and tsunami models created. Our focus is on highlighting how different assumptions on coupling contribute to the estimated hazards, so only coastal tsunami amplitudes are considered. We initially present each model separately, highlighting key features so that comparison between model results is clearer.
Although the rupture area and distribution of slip vary in size for a given magnitude bin, both the rupture area and maximum amount of slip on a subfault does on average increase with magnitude for each model following known scaling laws (Figure 6). At larger magnitudes, events begin to saturate in both length and width due to the actual CSZ fault dimensions. As a result of the von Karman ACF, slip is focused around asperities of favorable length scales, and as stated prior, asperity size predominantly scales with increasing magnitude. The stochastic nature of the models partitions and redistributes the asperities throughout the rupture extent. Irrespective of the presence of coupling models, ruptures varied between one another. Inherent to the stochastic nature of the rupture modeling, presence and influence of the coupling models is not always immediately obvious for an individual rupture; however, when all ruptures were combined into an averaged slip model, the result subsequently resembles the given desired mean slip model  (Figures 2b  and 3b). Figure 7 shows summary violin plots of the coastal amplitudes for each of the coupling models assumed and for the no coupling model cases. Each violin represents the kernel density estimate of the distribution of tsunami heights that is possible for a magnitude range at that specific point. Unsurprisingly, within all models calculated, the coastal tsunami arrivals on average increase in response to increasing rupture magnitude irrespective of the coastal latitudinal location. At lower magnitude scenarios (e.g., M8.4) arrival height differences are not as obvious between models. Slight variations in tsunami arrival amplitudes are present at lower magnitudes. Once magnitudes become great enough (∼M8.6) and larger swaths of the SMALL AND MELGAR 10.1029/2020JB021149 11 of 20 megathrust participate, differences in the coastal amplitudes that correlate to the particularities of the assumed coupling model become more obvious. The southern extent (below 43° latitude) of the study area on average experiences smaller variations between all classes, even at larger magnitudes. However, the resulting tsunami amplitude variations become increasingly prominent in the northern extent (46°-48° latitude) of the study area off the coast of Washington in comparison to the southern extent.
The Gamma scenarios, which assume larger and shallower offshore coupling, experience the largest tsunami amplitudes in the northern extent, whereas the southern extent is more dominantly subdued in the observed arrivals. The fully coupled zone in the Gamma model may extend throughout the entirety of the trench, but the total area of full coupling is widest starting at about 46° latitude and continuing to the edge of the study area. This separation matches well with the divide between the two regions of southern and northern arrival heights. For instance, at M9.1, the arrival are typically ∼1-3 m in the southern region whereas they are ∼4-6 m on average in the northern region. Unlike the Gamma model, there is not as clear of a separation of zones of high and low estimated tsunami amplitudes for the Gaussian coupling model. Interestingly, the Gaussian decade-scale coupling model contains a prominent region of coupling apparent in both the southern and northern region of the subduction zone; yet unlike the Gamma model, coupling here occurs between 15 and 25 km downdip at its highest coupling ratio. There is a deficiency of slip along the trench in the northern region in comparison to the Gamma model (46°-48°). Although large slip occurs closer to the coast in these scenarios, it also occurs much deeper than the previous models, limiting the surface displacement response and therefore producing a smaller volume of displaced water leading to overall smaller amplitude tsunamis (Figure 7). Tsunami amplitudes between the Gaussian and Gamma models appear relatively similar throughout the model region except near the Olympic Peninsula where coupling model differences are more pronounced.
The 30 km depth class does not see a prominent distinction between regions for the tsunami arrivals. There is a slight increase in the amplitude of arrivals in the northern gauge points than the southern (Figure 7), similar to but less obvious than in the previous two models. Mean recorded gauge point amplitudes do not exceed 4 m even for M9.1 events. Since these models lack any form of constraints to favor slip in particular regions along the slab interface, there is no distinction in the probability of activation of slip at any given SMALL AND MELGAR 10.1029/2020JB021149 12 of 20 subfault. With that in mind, mean consequential tsunamis arrivals should only vary due to the bathymetry and shape of the coastlines and should therefore have smaller variability laterally along the coast as a result of the ruptures themselves. This is similarly true for the 1 cm/yr class, however, more heterogeneity in amplitudes is present due to the more variable nature of the downdip limit. In the southern region, this class is constrained by a much shallower limit (∼15 km) than the northern region (∼25 km), limiting the region of slip to occur almost completely offshore in the south (Figure 1). Tsunami arrivals are more similar between the 1 cm/yr class and the two coupling classes. The 1 cm/yr class produces tsunami amplitudes typically ∼1 m greater than the Gamma class for M9.1 ruptures in southern Cascadia. In the north, however, this distinction reverses, where the Gamma (and less abundantly Gaussian) class produces slightly larger arrival heights than the 1 cm/yr class.
Coastal subsidence patterns, which influences observed tsunami hazards, vary widely between classes. For all rupture classes except the deeper 30 km depth class, the coast predominantly experiences subsidence (Figure 7). The mean coastal subsidence measurements for the Gaussian and Gamma models are strikingly similar in the southern region. In the Olympic Peninsula (47°-49°) coastal subsidence values vary more appreciably between the two. The 1 cm/yr class produces on average stronger coastal subsidence throughout the coastal region. This is most notable in the southern extent (41°-43°).

Discussion
Regardless of rupture class, the ruptures varied considerably between model scenarios for any given magnitude bin. The two different assumed coupling models clearly show differences in the depth dependence of slip. Figure 8 depicts three depth-averaged slip cross-sections for each assumed coupling model for all M8.8-M9.1 ruptures. As expected, the ruptures with a homogeneous mean model (no assumed coupling) contain almost no depth dependency for slip other than the imposed depth limit. This then allows slip to propagate as deep as the downdip limits and therefore slip is as probable to occur in these deeper regions as it is close to the trench. Mean slip varies considerably between these depth constrained models. For instance, in the southern region, where the 1 cm/yr class has a depth limit near the 15 km contour (Figure 1c), mean slip is almost twice as high compared to the 30 km depth class. Meanwhile, for the cases where a coupling pattern is imposed, there is a clear connection between the amount of slip and depth of the slab. The depth dependency of both the ruptures with imposed Gamma and the Gaussian coupling models becomes synonymous to the distribution of coupling for the two regions. In the northern and southern regions of the CSZ where slip deficit is the highest for both models, a dominantly monotonic decrease in amount of slip is seen for the Gamma class and pseudo-Gaussian slip shape is seen for the Gaussian class. As a result of this, the Gamma class experiences far greater slip shallower in depth when compared to the other classes, directly affecting the tsunami generation as seen in the hazard curves and the hazard map (Figures 9 and 10). Areas associated with the smallest coupling ratio see very little slip, predominantly in the southern region of the CSZ. Along strike variations in slip concentration are also present between the models where coupling is included (Figures 2, 3 and 8).
Coastal subsidence is prevalent in three of the four classes to a varying degree, for the 30 km depth class, the coastal signal is dominated by uplift (Figure 7). In recent work by Wirth and Frankel (2019) where they model coseismic uplift for M9.0-M9.2 CSZ megathrust events and compare them with coastal subsidence records tied to the 1700 Cascadia event, a dominance of coastal subsidence is recorded for ruptures with down-dip rupture limits corresponding to the 1 cm/yr locking contour. Although at current, coseismic uplift data for the CSZ is restricted to a finite number of locales (∼20 in total) based on paleoseismic estimates (e.g., Kemp et al., 2018;P. L. Wang et al., 2013), we would expect that similar trends in our models would add confidence to constraining class applicability. Since the 30 km depth class is dominated by coastal uplift, it seems more likely that this reference model may produce potential ruptures that are less likely to occur in the CSZ. This then adds weight to the importance of coupling for constraining ruptures, since even the 1 cm/yr class is based in some aspect on local coupling patterns.
It is important to reiterate that the coupling fraction does not perfectly correlate to where the slip occurs for any given rupture. Rather, the coupling as implemented here (Equation 9) will define regions where SMALL AND MELGAR 10.1029/2020JB021149 slip is more likely or less likely to occur. Since there is a sparse record of the earthquake history of the CSZ, coupling should only be used as a first order constraint. Further research on more earthquake cycles, past and future, will then allow us to better constrain the accuracy and applicability of coupling as it relates to the rupture characteristics.
Additionally, in our implementation we have assumed that a coupling ratio of 0 (fully creeping) corresponds to parts of the megathrust that cannot participate in the coseismic process. However, recent modeling studies have argued that dynamic effects can push the rupture into transition regions of the subduction zone. Areas of 0 coupling may, under certain circumstances, participate seismically (Ramos & Huang, 2019). An example in recent is the July 2020 M7.8 earthquake off the Alaskan-Aleutian arc is believed to have ruptured in the eastern portion of the previously assumed, seismically uncoupled Shumagin Gap ( Melgar, 2020). Observational studies have also shown potential overlap between regions with slow slip and regions with coseismic slip (e.g., Lin et al., 2020). However, since the dynamic effects and their influences in rupture propagation into transition zones are not yet well understood in terms of what controls them and how to quantify their likelihood of occurrence, we do not include them in the rupture modeling process. This has little effect on PTHA since shallow slip contributes far more to the overall hazard. However, for other applications it is potentially important as it would allow for slip further inland and closer to large population centers on the CSZ. Following the rupture model patterns, the associated tsunami hazards tend to increase toward the north of the study area. Figure 9 shows hazard curves at six coastal locations (shown in Figure 10) for a 1 year return period. These curves give a good quantitative measure to the associated tsunami hazards for the rupture models. Although these are site specific, it becomes quite clear that there are pronounced effects in the associated hazards when coupling is included in modeling. Differences in probabilities of even a relatively small tsunami (<1 m) are apparent between the 30 km depth class and the other three classes, however, for several gauge points throughout the region, distinguishing between the other three are difficult. Because hazard curves are site specific, associated hazards can vary between points in larger regions. In the southern and central CSZ, the 1 cm/yr produces more probable large tsunamis and the two coupling classes produce similar probabilities (Figures 9 and 10). Recall that for the 1 cm/yr class, slip in this region is more shallowly constrained ( Figure 10. The four distinct classes are colored accordingly. Curves are determined using a magnitude time distribution defined by the Gutenberg-Richter scaling relation and the method is described previously. It should be noted that the GR distribution does not constrain CSZ seismicity well at lower magnitudes. For more authoritative hazard assessments, other distributions should be explored. interface between 44° and 46° latitude (Figure 1), the slip deficit is much smaller, and therefore the coupling ratio is much smaller.
Tsunami amplitude arrivals and associated hazards between the two coupling classes are variable throughout the coastal region, with the most pronounced variations in the southernmost region near California and the northern region of the Olympic Peninsula. This is highlighted by the 2% probability of exceedance in 50 years hazard maps for the four classes depicted in Figure 10. The stark difference in amplitudes in the north highlights the effects the coupling models have on the ruptures. Between 46° and 42° latitude the Gaussian and Gamma coupling models have similar patterns of coupling. At higher latitudes, where the slab bends and the slab dip shallows, coupling varies more noticeably, with the strongest coupling for the Gaussian model occurring between 15 and 25 km depth and the Gamma model is strongly coupled from 15 km all the way to the trench. With larger and shallower coupling regions in the Gamma model, more slip is observed closer to the trench ( Figure 8) and therefore the resulting tsunamis appear much larger.
The hazard variability along strike observed in the stochastic slip models without imposed coupling is associated with variations present in the downdip limit of slip, tsunami propagation path, and the site conditions (the shape of the coastline). The tsunami path and site conditions are dominantly controlled by the bathymetry and distance from the tsunami initiation loci to the coast. The influence of the surficial expressions in hazard generations are the only sources of variability for the downdip limit classes, whereas the models with coupling are more heterogeneous. Around 44° latitude, the distance of the trench to the coast begins to widen quite drastically and the bathymetry becomes more heterogeneous over a larger distance. The estimated tsunami amplitude for gauge points south of 44° are on average noticeably smaller than those north. For the 30 km depth class, tsunami hazards are lower than the other three classes for all coastal sites ( Figure 10). This is not the case for the 1 cm/yr class. Coastal tsunami hazards associated with the 1 cm/yr class are greater for the southern region. In the norther region (above 46° latitude), the Gamma class experiences overall greater hazards than the 1 cm/yr class. In general, the differences between the coupling and 1 cm/yr classes are no greater than 1-2 m for the 2% hazard maps, yet some coastal points notice differences SMALL AND MELGAR 10.1029/2020JB021149 16 of 20 of 4-5 m. This variability has large implications on potential hazards and required tsunami hazard preparedness. Tsunami hazards associated with either coupling model do not vary significantly between either the Gaussian or Gamma models, but two regions of variable hazards persist: the southern Cascadia region around 41°-42° latitude as well as northern Cascadia around 47°-48° latitude.
One important point for assessing tsunami hazards in this way is whether the 200 ruptures simulated for each magnitude bin are sufficient to capture all the possible variability in the resulting tsunamis. In order to assess the statistical stability of the models, we refer to Figure 11 where the conditional probabilities, , are plotted for four magnitude bins against the number of scenarios observed. The conditional probabilities are relatively stable even after only ∼100 scenarios are included, however, our results detailed prior are determined using all 200 scenarios. Although here we plot the Gamma class, the results are similar for the three other classes. A formal CSZ hazard assessment should consider whether the stability observed in these conditional probabilities is sufficient or whether more scenarios are needed.
We note as well that we have assumed a fairly simple geometry for the shallow most megathrust with no splay faults. These have been inferred to exist in the CSZ (Booth-Rea et al., 2008) and coseismic motion on them can have a significant effect on the resulting tsunami (e.g., Gao et al., 2018). Similarly, we have considered only simple elastic deformation. Distributed deformation of the heavily sedimented deformation front (Han et al., 2017) could affect the ensuing tsunami as has been noted in places like Indonesia (Hill et al., 2012). Likewise, plastic deformation of the wedge can increase the amount of vertical deformation and lead to larger than expected tsunami (Ma & Nie, 2019).
The results discussed above are a compelling argument to include the influence of coupling for both the possible slip and rupture geometry generated as well as the subsequent tsunami hazards. However, without stronger confidence in the specifics of the coupling for the region, large uncertainties are still present. The two coupling models included here differ significantly, especially in the near-trench region. There is currently no way to distinguish between them because they produce the same onshore inter-seismic velocity field (Schmalzle et al., 2014). Additionally, there are other potential coupling models for the region (e.g., D. Li (Figures 9 and 10). The other three rupture classes produce similar probable stabilities at around 100 scenarios. as long-term viscoelastic relaxation of the upper mantle. Nonetheless, the two models we have used here, can be thought of as depicting two end member cases of the possible coupling distribution throughout the CSZ. The other existing coupling models share similar features to either of these two. Although dominant similarities are present between the ruptures and hazards associated with the coupling models, larger-scale variations are present. This large uncertainty between models highlights the necessity of the community to further seafloor instrumentation and seafloor geodesy in order to constrain the shallow coupling. Constraining coupling using hybrid modeling of both on land and seafloor geodetic instrumentation has already been implemented for the Nankai Trough along the coast of Japan (Yokota et al., 2016). Use of the coupling model present by Yokota et al. (2016) has also been implemented in tsunami wave propagation modeling from the two Kii Peninsula earthquakes of 2004 (Watanabe et al., 2018).
Finally, we note once more that these results should not be used as an authoritative hazard assessment for the CSZ. Further work remains to be done in order to improve accuracy for future applications: we have used a number of simplifying assumptions, considered only one kind of tsunami source, and cannot at this point distinguish which of the two coupling models is more likely. Our results instead highlight the large variations in rupture patterns and the resultant tsunami hazards created by including inter-seismic coupling, a key property of the megathrust. This shows that disregarding important information, like inter-seismic coupling, in rupture modeling can lead to hazards estimates that are potentially underestimated.

Conclusion
Current stochastic rupture modeling techniques do not consider the influence of some first-order fault zone characteristics like inter-seismic coupling patterns. A correlation between high slip patches and areas of greatest coupling ratios has been noticed for recent large earthquakes on megathrusts globally and suggests that this information should be taken into account. Here we presented an updated mathematical formalism to introduce coupling models as prior information into stochastic rupture modeling. We compared two end member models of coupling, one with a fully coupled zone extending to the trench and another with coupling deeper downdip and determined the variable influence that this has on the resulting rupture models. These were compared with two classes of ruptures defined by variable downdip limits of slip: a depth limit defined by the top of the tremor zone and a shallower depth limit defined by the 1 cm/yr coupling contour (25% coupling). Large variations in slip distribution are present in the rupture models and correlate well with areas with the largest differences in slip deficit. To exemplify the importance of this the ruptures were then used in a probabilistic tsunami hazard assessment of the CSZ. We found that the tsunami amplitudes generated are more uncertain throughout the region when coupling is used to condition the rupture models when compared to the tremor zone depth condition. Hazards associated with the 1 cm/yr class were more similar to the coupling classes than the previous class, however, important variations in resultant hazards were present throughout. More specifically, in the northern extent of the CSZ, associated hazards are larger where differences in coupling distribution are more drastic. Although large uncertainties are present in the accuracy of modeled coupling, imposing either constraint created variable hazard estimations. To bridge the uncertainties in present inter-seismic coupling and therefore refine our knowledge of first order fault characteristics for future hazard assessments, expanded seafloor instrumentation is vital.