Probabilistic seismic hazard analysis provides the evaluation of annual frequencies of exceedance of ground motion intensity measures (IMs) [typically designated by peak ground acceleration (PGA) or by linear elastic damped response spectral ordinates at specific periods of vibration] at a site. The result of a PSHA is a seismic hazard curve [annual frequency of exceedance versus ground-motion intensity measure (IM) amplitude], a uniform hazard spectrum (spectral amplitude versus structural period, for a fixed annual frequency of exceedance), or conditional mean spectrum (Baker, 2011; Lin et al., 2013). First reports of PSHA date back to the 1960s. Since then, PSHA has become the basis for seismic assessment and design of new and existing engineered facilities ranging from civil structures, such as buildings and bridges, to critical facilities, such as nuclear power plants. In PSHA, all possible earthquake fault sources contributing to the hazard need to be characterized first. Second, ground-motion prediction equations (GMPEs) are used to relate ground-motion IMs to variables describing earthquake source, path, and site effects. Extensive research has been performed on GMPEs for use in PSHA. Douglas (2003, 2011, 2016) summarized over 400 GMPEs that were developed since 1964–2016 for estimation of PGA and over 250 GMPEs for estimation of spectral ordinates at a site. Douglas and Edwards (2016) provides a recent discussion of current and future trends in ground-motion prediction. Stewart et al. (2015) provides a discussion of the selection of GMPEs for hazard assessments for the three principal tectonic regimes: active crustal regions, SZs, and stable continental regions for a global earthquake model. Of interest to this paper, Stewart et al. (2015) recommended the use of three models for SZ ground-motion predictions, “BC Hydro” model of Abrahamson et al. (2016), the global earthquake model described in Atkinson and Boore (2003), and the model in the study by Zhao et al. (2006). The BC Hydro model was developed using different data sets of SZ strong-motion recordings (e.g., Crouse et al., 1988; Crouse, 1991; Youngs et al., 1997; Atkinson and Boore, 2003, 2008; Zhao et al., 2006; Lin and Lee, 2008). While GMPEs are most often used for PSHA, it is worth noting that other methods for generating the IMs that involve the generation of synthetic ground motions have also been recently proposed, however, all involving extremely computational intensive methods. These methods that could serve as alternatives to the GMPEs include kinematic earthquake models (e.g., Olsen et al., 2008; Frankel et al., 2014; Pulido et al., 2015; Iwaki et al., 2016), stochastic finite-fault ground-motion methods (e.g., Atkinson et al., 2009), or hybrid broadband ground-motion methods (e.g., Atkinson et al., 2011; Skarlatoudis et al., 2015).

Relative to PSHA, only recently has there been much work on PTHA. Mori et al. (submitted)¹ summarized 30 PTHA studies conducted worldwide, all but one of which were conducted after the 2004 Indian Ocean tsunami. Prior to the 2004 event, tsunami hazards were generally characterized by “worst credible” or “worst case” scenarios. Rikitake and Aida (1988) were the first to use historical records of past run-up events to characterize the probability of the tsunami run-up hazard in Japan. Subsequent to the 2004 event and with the more recent SZ tsunami events in Chile in 2010 and in Japan in 2011, there has been an increasing interest in developing PTHA. Studies have focused on regions throughout the Pacific Rim, including Japan (Burroughs and Tebbens, 2005; Annaka et al., 2007; Yanagisawa et al., 2007; Fukutani et al., 2015; and Goda and Song, 2016), the US Pacific Coast and Canada (Geist and Parsons, 2006; González et al., 2009; Thio and Somerville, 2009; Priest et al., 2010; Witter et al., 2013; Leonard et al., 2014; and Park and Cox, 2016), South China Sea (Liu et al., 2007; Li et al., 2016), New Zealand, and Australia (Power et al., 2007, 2013; Burbidge et al., 2008; and Mueller et al., 2015), as well as places in Europe (Tinti et al., 2005; Grezio et al., 2010; Anita et al., 2012) and the Northwestern Indian Ocean (Thio et al., 2007; Heidarzadeh and Kijko, 2011). As explained in the study by see text footnote 1, the PTHA generally uses one of three approaches for the tsunami generation: (1) historical record approach, (2) logic-tree approach, and (3) random phase approach. Generally, the historical record of tsunamis is not sufficient to build a credible probabilistic model. Therefore, the second two approaches are favored. The logic-tree approach (e.g., Park and Cox, 2016) is based on combinations of slip conditions (e.g., magnitude, peak slip location, and slip distribution), and these combinations are given a weighting based on expert opinion, historical record, or equal weighting in some cases. For the random phase approach (e.g., Goda et al., 2014), the slip distributions are created by an assumed slip wavenumber spectrum using random phases. For this paper, the logic-tree approach is used since it is straightforward to combine PTHA with the PSHA. In general, the output of the PTHA focuses on the flow depth at the shoreline, the maximum extent of inundation for planning, the evaluation of annual exceedance probabilities (AEPs), or the estimation of the extent of damage through a probabilistic tsunami damage analysis (e.g., Wiebe and Cox, 2014; Park et al., 2017).

The authors are aware of only one study by De Risi and Goda, 2016 (DG16 hereafter), which considers the combined PSTHA. In that paper, the authors present a consistent method to account for a common rupture process. They use GMPEs based on magnitude and rupture distance to quantify the shaking intensity at a particular location subjected to a SZ earthquake, and then solve the shallow water wave equations for tsunami propagation and inundation from the source to the site of interest. The main output of their work is seismic hazard curves and tsunami hazard curves that represents the mean annual rate (MAR) of exceedance of a given IM. Their generalized framework was applied for assessing combined earthquake and tsunami hazard at a single location on the coast line of Sendai City based on the subduction fault plane in the Tohoku region of Japan.

Similar to the work of DG16, the objective of this work is to develop a consistent framework for a multi-hazard analysis considering large magnitude SZ earthquakes and a subsequent near-field tsunami. However, the study performed herein is developed for a different site with a different methodology to define earthquake and tsunami sources when compared to DG16. The long-term objective is to be able to use the PSTHA as the basis for a probabilistic multi-hazard damage assessment to quantify the separate contributions of seismicity and tsunami hazards in the estimation of damage to the built environment at a community scale or for design of specific infrastructure elements such as a critical facility. A general methodology is first presented in Section “Methodology,” with a general review of the combined probabilistic hazard analysis (PHA) methodology proposed (see General PHA Methodology), the earthquake fault source models and their characteristics (see Earthquake Fault Source Models and Their Characteristics), the GMPEs used (see Earthquake Simulation), and the tsunami generation, propagation, and inundation (see Tsunami Generation, Propagation, and Inundation). Then, the multi-hazard logic-tree model is presented (see Multi-Hazard Logic-Tree Model), and the methods to estimate the AEP (see Estimates of the AEP) are presented next. In Section “Application: Full-Rapture Event of the CSZ Impacting the City of Seaside, OR,” the methodology is applied to an application example in which the combined seismic and tsunami hazard are quantified for the City of Seaside on the Oregon Coast of the US Pacific Northwest, conditional on a full rupture of the CSZ. In this application, the CSZ earthquake source model, and the geological, geographic, and morphological features of the City of Seaside, Oregon, are characterized first (see Characterization of the Source) and then details of the CSZ fault model and tsunami model (see CSZ Fault Modeling) are provided. In Section “Results,” results are presented in terms of the seismicity (see Seismicity) and tsunami intensity (see Tsunami Intensity) estimated for the area of interest. A spatial representation of both the seismic and tsunami intensity is presented, including the study of the granularity needed to characterize the hazards and joint distribution for different vector-valued IMs (see Spatial Representations Seismic and Tsunami Hazard). Finally, in Section “Summary, Conclusion, and Future Work,” discussion and summary of the results are provided, followed by the conclusion and recommendations for future research.

The general formulation presented here is aimed at developing the probabilistic earthquake and tsunami hazard at a site. In this formulation, the first topic that has to be addressed is the definition of an IM for the hazard. For example, with respect to earthquake ground shaking, typical measures of interest are spectral acceleration (S_a), while common IMs used for the tsunami hazard are the inundation flow depth (h), flow velocity (V), and specific momentum flux (M_F = hV²). Since there is great uncertainty in the quantification of the earthquake that induces ground shaking and possible tsunamis, be it in the focal mechanism including definition of the location (e.g., location of the epicenter, extension of fault rupture), size (magnitude), and resulting intensity of a future earthquake at a specific site of interest, PSHA (e.g., Cornell, 1968; McGuire, 1995; Kramer, 1996) was developed as an analytical tool to characterize the seismic hazard probabilistically. PSHA has become the most widely used method for assessing seismic hazard at a specific site. More recently, PTHA has also been developed (PTHA—e.g., Geist and Parsons, 2006; Annaka et al., 2007; Power et al., 2007) as reviewed by see text footnote 1 and summarized in Section “Background and Literature Review.”

The general PHA provides the MAR(λ) of IM exceeding an intensity measure value im that is computed using the Total Probability Theorem (Benjamin and Cornell, 1970), by integrating the contributions of all possible tsunami-seismogenic sources, and for each of the sources, all possible values of earthquake magnitude as (1)λIM>im(im)=∑i=1Nsourcesλi(M≥mmin)×∫Θ∫mminmmaxP(IM>imΘ,M)sΘi(ΘM)fMi(m)dm dΘ where N_sources denotes the total number of seismic sources contributing to the hazard at the site, λ_i(M ≥ m_min) is the MAR of occurrence of earthquakes with magnitude greater than a lower bound threshold value, m_min, of seismic source i, P (IM > im|Θ,M) represents the probability that intensity measure IM will exceed a given intensity value im at the site conditional on a given magnitude M and source parameters Θ, sΘi(ΘM) represents the characteristic probability density function (PDF) of earthquake source parameters Θ obtained from a hazard parameter prediction model conditional on the magnitude of the earthquake, and the function fMi(m) denotes the PDF of the magnitude M given the occurrence of an earthquake on seismic source i.

In Eq. 1, it is assumed that earthquake occurrences at different seismic sources are statistically independent (in terms of occurrence time, scaling relationship, M, etc.), implying that earthquake occurrences from all possible sources can be assumed to follow a Poisson process. It is also assumed that within each seismic source i, the magnitude M_i earthquake events are statistically independent. Thus, the summation in Eq. 1 considers the contributions from all seismic sources while the integrations over M_i and source parameters account for earthquakes of all possible magnitudes and source parameters conditional on the magnitude for each seismic source, respectively. In Eq. 1, the conditional probability P(IM > im|Θ,M), with im > 0, represents the complementary cumulative distribution function (CCDF) of the IM conditional on M and Θ. Figure 1 illustrates the steps that can be used to compute the probabilistic seismic and tsunami hazard relationships. In Step 1, the earthquake fault source models and characteristics of the earthquake source models are defined. In Step 2, the ground-motion IMs at a given site are obtained through the earthquake simulation, performed through either explicit source-to-site wave propagation models and synthetic ground-motion generation or, most commonly, with GMPEs. In Step 3, the tsunami is simulated, including tsunami generation, propagation and inundation modeling. In Step 4, the seismic-tsunami hazard curves and surfaces are determined. It is notable that Step 2 (seismic simulations) and Step 3 (tsunami simulations) are produced from a consistent earthquake event.

To describe the distribution of earthquake magnitudes in a given region of interest, the Gutenberg–Richter relationship (GR) is widely adopted. The GR relationship is given by (2)logλm=a−bm where λ_m is the MAR of exceedance of an earthquake of magnitude m, a represents the overall rate of earthquakes in a region of interest, and b represents the relative ratio of small and large magnitudes. The parameters a and b are estimated based on statistical analysis of the database of seismicity for the seismic source zone of interest. The GR relationship is developed from a regional dataset of seismicity accounting for many different source zones and has been found to be inadequate to represent the earthquake recurrence relationship for the tail-end of the magnitude–frequency distribution representing large magnitude earthquakes (Schwartz and Coppersmith, 1984; Wesnousky, 1994). Several models have been proposed to address the shortcomings of the GR recurrence law (Kagan, 1997, 2002a; Bird and Kagan, 2004). For SZs, such as the CSZ used as an example in this paper, a Tapered Gutenberg–Richter (TGR) distribution has been shown to be a robust model (Rong et al., 2014). The TGR is expressed as a function of the seismic moment M₀, instead of magnitude m and an exponential taper is applied to the number of events with very large seismic moment. The TGR CCDF is given by Kagan (2002a) (3)F(M0)=M0tM0βexpM0t−M0M0c,for M0t≤M0<∞ where β is the index parameter of the distribution, and β = (2/3)b, M_0c is the corner moment, and M_0t is the threshold moment above which the earthquake catalog is assumed to be complete. The conversion between seismic moment M₀ and moment magnitude m is given by (4)M0=101.5m+C where C = 9 − 9.1. The corner moment M_0c can be estimated using the seismic moment conversion principle (Kagan, 2002b) and given approximately by (5)M0c=χM˙T(1−β)αtM0tβΓ(2−β)1/(1−β) where χ is the seismic coupling coefficient, M˙T is the tectonic moment rate, α_t is the recurrence rate for earthquake with moment M_0t and greater, and Γ is the gamma function. Equation 5 can be used to derive the recurrence interval as (6)T(M0t)=11−βM0tβM0c1−βM˙TΓ2−βexpM0tM0c

The TGR CCDF can be rewritten in terms of magnitude m and given as (7)F(m)=101.5mt−mβexp101.5mt−mc−101.5m−mc

Rong et al. (2014) estimated the probable maximum earthquake that is likely to occur within a given time for circum-Pacific SZs using TGR distributions. For the CSZ, using maximum likelihood estimation method, the values of β = 0.59 and m_c = 9.02 were estimated considering the 10,000-year paleoseismic record based on the turbidite studies by Goldfinger et al. (2012), coupled with the limited number of instrumental earthquake data. In the implementation performed, it is convenient to convert the continuous distribution of magnitudes into a discrete set of possible magnitudes m_j, which are given by (8)P[M=mj]=G(mj+0.5Δm)−G(mj−0.5Δm) where the G(m) = 1 − F(m) is the cumulative density function and Δm is the adopted discretization interval.

Figure 2 shows the median earthquake recurrence relationship developed for the CSZ using a TGR distribution (Eqs 3–7) and with β = 0.59 and m_c = 9.02, as recommended by Rong et al. (2014) for the CSZ. According to the TGR distribution shown in Figure 2, m ≥ 8.8 earthquakes are expected with a return period of 500 years (λ_m = 0.002), while m ≥ 9.0 earthquakes are expected with a return period of 1,000 years (λ_m = 0.001). Goldfinger et al. (2012) reconstructed the large earthquake history of the CSZ for approximately10,000 years based on strong shaking-induced turbidite deposits in marine sediments and onshore paleoseismic records. The study suggested four types of earthquake rupture along the CSZ based on the interpretation of the turbidite data during the past 10,000 years: (1) 19–20 full-margin or nearly full-margin ruptures, (2) 3–4 ruptures along the 50–70% of the southern margins, (3) 10–12 southern ruptures from central Oregon southward, and (4) 7–8 southern Oregon/northern California ruptures. Though the turbidite data do not provide direct indication of the probable earthquake magnitudes, Goldfinger et al. (2012) estimated the earthquake magnitudes of different rupture events based on the relations observed among the rupture length (distance between offshore core sites containing turbidites from same events), turbidite thickness, and turbidite mass and estimated that full-rupture events constituted m = 8.7~9.3. Considering 20 full rupture over the past 10,000 years, the MAR of full-rupture events λfull−rupture≈2010,000=0.002, which is consistent with the λ_m≥8.8 = 0.002 estimated by Rong et al. (2014) using TGR distribution shown in Figure 2.

Near-field ground motions are strongly affected by the heterogeneity of earthquake rupture processes, such as slip distribution, rupture directivity, and the acceleration and deceleration of the rupture front. To estimate the ground motion quantitatively for a seismic hazard assessment, characterization of this heterogeneity is essential, and usually accounted for in ground-motion hybrid broadband simulation procedures (e.g., Somerville et al., 2012) or 3D simulations of earthquake event scenarios (e.g., Olsen et al., 2008; Delorey et al., 2014). While these explicit source-to-site wave propagation models and synthetic ground-motion generation models are extremely useful, especially in generating synthetic waveforms, these typically involve very large number of computations in large parallel computing centers, which is currently still not feasible when performing probabilistic seismic hazard analyses. Instead, when performing PSHA, GMPEs are still typically used today.

In this study, GMPEs, which are dependent on the local and regional site conditions, are used to simulate the ground-motion IMs. The recently developed Abrahamson et al. (2016) GMPE, which is based on the global datasets on SZ earthquakes described in Section “Background and Literature Review” is used. The functional form of this GMPE for interface SZ earthquakes is given by: (9)lnSainterface=θ1+θ4ΔC1+θ2+θ3M−7.8lnRrup+C4expθ9(M−6)+θ6Rrup+fMAG(M)+fFABA(Rrup)+fsite(PGA1,000,Vs30)+σϵ where θ_i are regression parameters, R_rup is the closest distance to the rupture area from site, Vs₃₀ is the shear wave velocity of the uppermost 30 m of soil, PGA_1,000 is the median PGA corresponding to Vs₃₀ = 1,000 m/s, σ is the total SD, which is obtained combining intra-event uncertainty ϕ and inter-event uncertainty τ, and ϵ is the standard normal error term. In this study, only the intra-event uncertainty ϕ is considered since only the CSZ is used for hazard analysis. The magnitude scaling term is given by (10)fMAG(M)=θ4M−C1+ΔC1+θ1310−M2,for M≤C1+ΔC1θ5M−C1+ΔC1+θ1310−M2,for M>C1+ΔC1 where C₁ = 7.8, and the ΔC₁ term represents the epistemic uncertainty around the distinct break in magnitude scaling between frequent smaller magnitudes events and rare large interface events. The forearc and backarc scaling term in Eq. 9 is given by (11)fFABA(M)=θ15+θ16lnmaxRrup,10040FFABA;FFABA=0,for forearc or unknow sites;1,for backarc sites.

The model for site response scaling is given by (12)fsitePGA1,000,Vs30=θ12lnVS∗Vlin−bln(PGA1,000+c)+b lnPGA1,000+cVS∗Vlinn,for Vs30<Vlinθ12 lnVS∗Vlin+b lnVS∗Vlin,for Vs30≥VlinVS∗=1,000,for Vs30>1,000;Vs30,for Vs30≤1,000.

All other model coefficients in Eqs 9–12 are listed in the study by Abrahamson et al. (2016).

In general terms, the tsunami hazard at a particular site requires the three steps of (1) tsunami generation, (2) propagation, and (3) inundation. For most modeling efforts, tsunami generation is given as an initial surface water displacement along the fault. The Okada (1985) model, which is based on the linear co-seismic dislocation of fault slips, is often used for simplicity, and the initial displacement is assumed to occur simultaneously along the fault.

Tsunami propagation is generally considered a solved problem in that the equations are well defined for long wave propagation in the open ocean. Of course, the propagation phase requires accurate knowledge of the underlying bathymetry which affects the wave through refraction, diffraction, and shoaling. The third phase, tsunami inundation, considers the flow of water over dry land. This is considered a difficult problem to solve because of the complex interaction of the flow with the built and natural environment that is also changing due to the destructive nature of the flow. However, most inundation models assume “bare earth” conditions, that is a digital elevation model (DEM) in which the natural and built environments are removed. Typically, the effects of the vegetation and structures are replaced by a suitable friction factor, although there are additional uncertainties in this step (e.g., Park et al., 2013; Bricker et al., 2015).

In this study, the logic-tree model by Park and Cox (2016) is applied for tsunami generation to characterize the fault slip. The slip model from the study by Park and Cox (2016) characterizes the randomness of fault slip distribution at the CSZ as a Gaussian shape, parameterized in terms of the moment magnitude, peak slip location, and a fault slip shape as (13)f(Y′/dL|α,β)=1β2πexp−Y′/dL−α22β2 where α and β are the slip distribution parameters along a rupture strike direction (Y ′), and dL is the unit length of sub-fault utilized in the slip model. Each parameter α and β controls the location of the peak slip and the shape of slips. A total of 72 scenarios from the three seismic moments, three slip shapes, and eight peak slip locations are proposed for the full-length rupture CSZ event. The occurrence rate of each seismic moment estimated from the paleoseismic data at CSZ (Goldfinger et al., 2012). Each of the fault slip distributions are determined from the slip model and applied to the ComMIT/MOST model (Titov et al., 2011) as an input to simulate the tsunami generation and propagation parts. Although MOST can also be used for the inundation phase, the software COULWAVE (Lynett et al., 2002) is used to model this final step (Park and Cox, 2016).

The generic multi-hazard logic tree is presented in Figure 3. The first step is to identify all tsunamigenic earthquakes that could potentially affect the site of interest. These earthquakes are then classified as near-field or far-field earthquakes depending on the source location with respect to the site. For each of the potential fault sources, the next step is to use an appropriate recurrence model. As mentioned in Section “Earthquake Fault Source Models and Their Characteristics,” the widely used modified-GR underpredicts the recurrence rate for large magnitude tsunamigenic earthquakes and hence TGR (Kagan, 2002a) and other characteristic moment–frequency distributions (e.g., Wesnousky, 1994) are preferred. Since most of the GMPEs use a distance metric of closest distance from the site to the rupture plane (R_rup), the depth of rupture which determines the position of the rupture edge close to the site of interest can have significant effect on the ground-motion intensity observed. Several geophysical models are available to determine the rupture depth and can be used for a specific fault of interest. As for example, for 2014 National Seismic Hazard Map of United States, three geophysical models were used to constrain the eastern edge of the CSZ rupture zone (Frankel et al., 2015).

Reliable magnitude scaling relationships for subduction earthquakes is a prerequisite for accurate estimation of earthquake and tsunami hazard intensities. There are several magnitude scaling relationships available in the literature that are derived from global observation of earthquakes on the plate interface of SZ (Papazachos et al., 2004; Strasser et al., 2010; Murotani et al., 2013; Goda et al., 2016; Skarlatoudis et al., 2016). Of these, Skarlatoudis et al. (2016) compiled an updated database of interface earthquakes that occurred worldwide in last decade in the major SZs (e.g., 2004 M 9.1 Sumatra, 2010 M 8.8 Chile, 2011M 9.0 Japan earthquake) and proposed new and improved source scaling laws with reduced uncertainty compared to other currently available source scaling laws for subduction earthquakes. Once the appropriate source magnitude and scaling laws are identified, the final step is to perform the earthquake and tsunami simulations for a multiple logic-tree branch. In the case of earthquakes, GMPEs are utilized, which typically make use of the M and R_rup as the input parameter to compute the ground-motion IMs at the site. On the other hand, for tsunamis, the details of the slip distribution (e.g., slip shape, peak slip location within rupture zones) and seismic moment for a given rupture are critical input parameters needed to compute the tsunami hazard intensity at a given site.

Four IMs were selected for the PSTHA, including the PGA, the 5% damped linear elastic spectral acceleration at a fundamental period of vibration of 0.3 s [S_a(T₁) = 0.3 s], the maximum flow depth (h_max), and the specific maximum momentum flux (M_F)_max because these IMs are often linked with structural damage at a structure scale (e.g., Faggella et al., 2013; Gidaris et al., 2016). Here, the flow depth is the net elevation of the free surface elevation above the local land elevation, and the specific momentum flux is given by the product of the flow depth with the square of the flow velocity (M_F = hV²).

The Poisson process (Cornell, 1968) is used to estimate the AEP of both earthquake and tsunami hazard curves at a specific location. The probability of IMs exceeding a certain level of that hazard conditional on a given time t is equal to the probability of at least one event occurring in time t, and is given by (14)PIM>im|t=1−e−λt where λ is the mean occurrence rate at which the IM will exceed a specific im at a given location during the time t. The MAR exceedance of each IM is computed using Eq. 1.

The AEP surfaces for the joint seismic-tsunami hazard can be computed using the formulation presented next, designated here as vector-valued probabilistic seismic and tsunami hazard analysis (VPSTHA). The presentation starts from the definition of the joint mean rate density (MRD) of the hazard (e.g., Bazzurro, 1998; Barbosa, 2011), here expanded to account for IMs related to ground-motion shaking intensity as well as tsunami IMs. For a vector of ground-motion and/or tsunami intensity parameters IM = {IM₁,IM₂}, the joint MRD of the hazard is given by (15)MRDIM1,IM2(im1,im2)=∑i=1Nsourcesλi(M≥mmin)×∫Θ∫MfIM1,IM2(im1,im2Θ,M)sΘi(ΘM)fMi(m)dm dΘ where fIM1,IM2im1,im2Θ,M is the joint PDF of IM₁ and IM₂ conditional on earthquake of magnitude M and source parameters Θ. Once the MRD is defined, the MAR of events at the site with IM₁ and IM₂ being between im_1,1 < IM₁ ≤ im_1,2 and im_2,1 < IM₂ ≤ im_2,2, respectively, is given by: (16)λIM1∈im1,1,im1,2,IM2∈im2,1,im2,2=∫im1,1im1,2∫im2,1im2,2MRDim1,im2dim2dim1

It is worth highlighting that the vector of intensity parameters IM = {IM₁,IM₂} may contain ground-motion scalar IMs and/or tsunami scalar IMs. Examples are vector-valued IMs such as IM = {S_a(T₁),h_max}, IM = {PGA,S_a(T₁)}, IM = {h_max,(M_F)_max} The formulation in Eqs 15 and 16 is generic, but its implementation requires the computation of joint MRDs, which can be computationally expensive.

The northern coast of the North American continent from Northern California in the United States to Vancouver Island in Canada is facing the threat of a megathrust earthquake event with near-field tsunami from the CSZ along the converging plate boundary between Juan de Fuca Plate and North American Plate. The Juan de Fuca Plate is sinking beneath the North American Plate with the mean rate of 0.04 m/year (Heaton and Hartzell, 1987) to the northeast direction (Figure 4). The accumulated potential energy between two plates is released in the megathrust earthquake events and can cause ground shaking and rapid displacement of the seafloor which generates the initial deformation of surface water. Each megathrust rupture of the converging plate boundary triggers the earthquake and tsunami event in both offshore and onshore directions. The last full-rupture event occurred on January 26, 1700, with moment magnitude estimated between 8.7 and 9.2 (Satake et al., 2003). It was also reported that there were smaller but more frequent partial rupture events on the north or south margins of the CSZ (Atwater and Griggs, 2012; Goldfinger et al., 2012). However, in this application example of the proposed multi-hazard assessment, the partial rupture events are not considered to reduce the number of computations needed for the tsunami inundation.

Seaside, Oregon, is the study area chosen to demonstrate the proposed PSTHA methodology. Figure 5A shows an aerial view of Seaside with approximately 4 km of shoreline facing the Pacific Ocean and two small rivers, the Necanicum River and Neawanna Creek, that run parallel to the shoreline. The city has a population of approximately 6,500 residents, and the number of people in the city can increase to over 20,000 during the summer tourist season. There are over 5,700 buildings in Seaside with the larger hotels constructed out of steel and reinforced concrete in the center of the city (highlighted by the red box in Figure 5A) and surrounded by mostly wooden residential structures to the north and south of the city center. Figure 5B shows the bathymetry and topography of the study area. The remnant coastal dune on which the city was built can be seen running parallel to the shoreline with a secondary rise between the Necanicum River and Neawanna Creek. To the east of the Neawanna Creek, there is a steep gradient leading to the foot hills and a large headland to the south west. Figure 5C shows the distribution of the soil classes for this region. The steeper mountain areas are considered as Class C, and a majority of the city area is Class D (ASCE, 2010). The soil class distribution is also aligned approximately in the shore parallel direction.

Seaside has been the subject of multiple tsunami studies in the past because it is considered to be one of the most vulnerable locations to a CSZ event. As reported by Wood et al. (2010), the low lying city has approximately 87% of its land within the inundation zone of a CSZ full-rupture event, and 89% of the employees work within this zone. Because of the high vulnerability of Seaside to a CSZ event, there have been several recent studies in this area as reviewed by Park and Cox (2016) and Park et al. (2017), including the work by González et al. (2009) that conducted a PTHA using 14 historic tsunami far-field scenarios and 12 scenarios of the CSZ event.

Figure 6 shows the logic-tree model utilized for both earthquake and tsunami for CSZ full-rupture event (Park and Cox, 2016). This logic-tree model characterizes the randomness of the fault slip in terms of the moment magnitude, peak slip location, and a fault slip shape. To simplify the scenarios, 27 sub-faults distributed along the CSZ were considered, which are the default sub-fault setup for the ComMIT model (Titov et al., 2011). The sub-faults are distributed from the northern Vancouver Island to Northern California. Each sub-fault is 100 km long by 50 km wide, as shown in Figure 4. The detailed geologic information of the location and slip information is summarized in Table 1. Each of the 27 sub-faults has a constant rake of 90° and has varied strike, dip, and depth conditions along the entire fault. The x and y coordinates indicate the middle point on the right edge boundary of each sub-faults in Figure 4.

The full-rupture event at the CSZ is discretized as three moment magnitude scenarios (M 8.8, 9.0, and 9.2), and the corresponding rupture areas (S) were estimated using the relationship between seismic moment M₀ and rupture area provided by Murotani et al. (2013) and given by (17)S=1.34×10−10M02/3

Park and Cox (2016) added upper and lower limits to the surface area moment magnitude relationships of Eq. 17 by multiplying Eq. 17 with 2.0 and 0.5, respectively, and these limits bounded most of historic recent tsunami events at SZ around the Pacific Ocean which are generated by earthquake scenarios with M > 8.5. Including those two limits, a total of three slip shapes along the strike direction were utilized as possible scenarios per moment magnitude condition. Additionally, 8 possible peak locations are considered along the rupture length, giving rise to 24 scenarios applied for each moment magnitude condition. A total of 72 scenarios (3 moment magnitudes × 3 slip shapes × 8 peak slip locations) are proposed here to characterize the full-rupture CSZ event. Park and Cox (2016) applied 72 scenarios as inputs to the ComMIT/MOST (Titov et al., 2011) model, which considers the elastic dislocation model (Okada, 1985) for tsunami generation and solves the non-linear shallow water equations implemented with a finite difference scheme. The ComMIT/MOST model consists of three nested grids (A, B, and C-Grid) shown in Figure 4. For the inundation modeling, the results of COULWAVE (Lynett et al., 2002) were used, which solves a set of Boussinesq equations with a high-order finite-volume method by using the output of ComMIT/MOST modeling at B-Grid as the input of COULWAVE at C-Grid. All bathymetry data were originated from NOAA’s National Geophysical Data Center and the DEM for the Seaside area whose resolution is 1/3 arc second was used for the C-Grid. Each A, B, and C-Grids has 1 min (400 × 400), 3 s (800 × 800), and 24 m (416 × 390) resolutions. As a tide condition, the mean high water level was fixed as conservative tsunami hazards estimation in this study. The model results provide surface elevation and velocity time series over the entire study area at the 24 m grid resolution. More details of fault slip distributions and tsunami simulations are available in the study by Park and Cox (2016). In case of earthquakes, the same fault model is used. However, the BC Hydro GMPE is utilized in this study, which is based on the closest rupture distance, R_rup, between the site of interest and rupture surface for the interface earthquakes considered.

Though larger magnitude earthquakes generally have higher damage potential compared to smaller magnitude earthquakes, smaller and frequent events can be dominant contributors to the structural damage risk measured in terms of mean annual frequency of exceeding a damage state such as the collapse damage state (e.g., Zareian and Krawinkler, 2007; Eads et al., 2013). Figure 7A shows the probability mass function (PMF) of the earthquakes in CSZ deemed capable of generating tsunamis. Tsunamigenic eartquakes with lowest magnitude of m_min = 7.3 is assumed and largest magnitude is assumed to be m_max = 9.3, consistent with the earthquake expected over a 10,000-year period in CSZ (Rong et al., 2014). A discretization interval of Δm = 0.2 is adopted to develop the PMF of Figure 7, resulting in 10 central magnitude values m = 7.4 to m = 9.2. The PMF for different magnitudes are computed based on the TGR distribution described in Section “Earthquake Fault Source Models and Their Characteristics” and using Eq. 8. Of these 10 magnitudes, m ≥ 8.8 is assumed to produce the full rupture along the CSZ following the recommendations by Goldfinger et al. (2012), which is also consistent with Park and Cox (2016). These full-rupture events (m = 8.8, 9.0, 9.2) are only considered for earthquake-tsunami hazard analysis for Seaside. According to Figure 7A, λ_m≥8.8 = 0.002 results in 500-year return period for m ≥ 8.8 representing full-rupture scenarios. To be consistent with 526-year return period used in Park and Cox (2016) for full-rupture CSZ events, λ_m≥8.8 = 0.0019 is used for the PSHA computations for the study area of City of Seaside.

One of the inputs required in the GMPEs for computing ground-motion intensities at a given site is the source-to-site distance. The BC Hydro GMPE used in this study is based on the closest rupture distance R_rup between the site of interest and the rupture surface for interface SZ earthquakes. To compute the R_rup for different locations in Seaside for different magnitude scenarios, the R_rup for all the C-grid locations are precomputed from the 27 sub-fault areas shown in Figure 4. The rupture surface for different magnitude scenarios of Figure 7A are then computed based on the magnitude scaling law by Murotani et al. (2013), which are randomly positioned along the 27 sub-faults corresponding to the approximately 1,000 km long by 150 m width fault plane. The randomly positioned rupture surfaces corresponding to different scenarios are subsequently used to interpolate the R_rup for each scenario from the precomputed closest distances. Figure 7B shows the PMF of R_rup computed for observation point 1 in Figure 5A. The closest distance computed for all the scenarios is 27.5 km, which is the distance computed from the fault area immediately below the C-grid. It can be observed in Figure 7B that smaller variability of R_rup is observed for larger magnitudes. For example, for m = 9.0 that is a full-rupture event results in a single source-to-site distance and the resulting R_rup can be considered as a deterministic event with P[R_rup = 27.5 km] = 1. For the partial rupture events, the rupture surface is randomly positioned within the fault plane, which results in variability in the closest rupture distance computed as shown in Figure 7B.

AEP for two earthquake IMs, PGA and S_a(T₁ = 0.3 s), are computed for three observation points (shown in Figure 5A) in Seaside. These points are aligned shore-normal along the urban center of Seaside with approximately 400 m distance apart from each other. The observation Points 1, 2, and 3 are located between shoreline and the Necanicum River, Necanicum River and the Neawanna Creek, and Neawanna Creek and the edge of inundation zone, respectively.

Figure 8 shows the AEP of PGA and S_a(T₁ = 0.3 s) at three observation points computed for the full-rupture scenarios (M = 8.8,M = 9.0,M = 9.2) using the BC Hydro GMPE. As per the soil site class map of Figure 5C, point 1 and point 2 fall in site class D where as point 3 is in site class C. The V_s30 values assigned to Point 1, Point 2, and Point 3 are 300, 360, and 420 m/s, respectively. A value of λ_m = 0.0019 corresponding to M ≥ 8.8 is used for the hazard curves and subsequent AEP computations. In AEP computations, the PMF values considered for the three magnitudes is consistent with those used in the study by Park and Cox (2016), which were obtained from the study by Goldfinger et al. (2012). Since the weights assigned to the M 8.8, M 9.0, and M 9.2 are 5/19, 13/19, and 1/19, respectively, the M 9.0 is the dominant contributor to the hazard curves presented in Figure 8 for both PGA and S_a(T₁ = 0.3 s), followed by M 8.8, and M 9.2, respectively. As shown in Figure 8, AEP of PGA for three points are almost identical whereas slight variation in AEP for S_a(T₁ = 0.3 s) is observed for those three points. The slight variation in AEP for S_a(T₁ = 0.3 s) is due to the different V_s30 values assigned to those points. Overall AEP of both IMs are insensitive to their geographical conditions at this 800 m distance between the three points, which are all approximately 27.5 km from the fault.

The joint MAR of exceedance of IMs contour plot of PGA and S_a(T₁ = 0.3 s) is shown in Figure 9 for the three points indicated in Figure 5 for the 1,000-year event. The joint MAR of exceedance of IMs is computed based on the formulation presented in Eqs 15 and 16, in which the joint probability function of the two IMs conditional on the magnitude at a given site is computed assuming that the IMs follow a joint lognormal distribution. For the jont MAR of exceedance computation of PGA and S_a(T₁ = 0.3 s), correlation between the two IMs are considered and computed based on the study by Baker and Jayaram (2008). It can be seen from Figure 9 that for all three observation points, the most likely joint intensities occur for values of PGA = 0.5 g and S_a(T₁ = 0.3 s) = 1.0 g. Moreover, there is no discernable difference in the joint hazard intensities among the observation points considered as can be observed in Figure 9, as expected.

Figure 10 shows the AEP of two tsunami IMs h_max and (M_F)_max at the three observation points similar to Figure 8. Generally, the tsunami IMs decrease as distance from the shoreline increases because of the dissipation of energy during the inundation process as shown for the AEP of h_max (Figure 10A) and (M_F)_max (Figure 10B). For example, the AEP of both h_max and (M_F)_max at Point 1 are higher than Points 2 and 3. However, for lower probability events in the range of 0 < h_max ≤ 3 m, Points 2 and 3 have nearly the same AEP. Point 3 has higher h_max at the higher probability events in the range of 3 < h_max ≤ 10 m and also has the highest h_max due to the local topographic effects, including effects of the river and the creek running parallel to the shoreline (Figure 5A), highlighting the sensitivity of the tsunami IM to site-specific conditions. On the other hand, (M_F)_max shown in Figure 10B has three distinct curves with a generally decreasing trend from the shoreline to inundation limits. Comparing Figures 8 and 10, it is noted that the variation of AEP of tsunami IMs among three observation points are qualitatively dissimilar to the results of earthquakes IMs in that there are strong spatial gradients for the tsunami IMs across the length scale of the city. This is somewhat expected because the AEP of both PGA and S_a(T₁ = 0.3 s) depend primarily on the soil types and R_rup that have relatively small variations over the study region at Seaside, especially for the full-rupture scenarios considered. Moreover, this difference underscores the differences in the fundamental physics of the propagation of the seismic energy through the subsurface and the propagation of the hydrodynamic tsunami energy.

Figures 11A–C shows the joint AEP of h_max and (M_F)_max at Point 1, 2, and 3, respectively, where the same number of bins for both IMs are utilized to calculate the vector-valued IMs of the joint h_max and (M_F)_max from Eqs 15 and 16. Figure 11 also includes the isolines for four Froude numbers (Fr = 0.2, 0.5, 1.0, and 2.0) to show the possible range of distributions of velocity fields at the given flow depth. In general, the results show the high correlation between h_max and (M_F)_max at all three observation points, and each joint AEP follows a particular Fr for that area. For example, at Point 1 when h_max is in the range 6 < h_max < 7 m, the corresponding range of (M_F)_max is 150 < (M_F)_max < 200 m³/s² and corresponds to a Fr range of approximately 0.6 < Fr < 0.7. The Fr generally decreases shoreward, and seems to be everywhere subcritical (Fr < 1.0) which is physically realistic based on observed inundations following the 2011 Tohoku tsunami. This is not to say that the flow is always subcritical, just that at the point of maximum flow depth or momentum flux, that the flow can be expected to be subcritical.

Figure 11 reveals that the different distributions of the joint surface of the AEP following the specific Fr at three observation points originate from the local bathymetric/topographic conditions and not from the fault slip distributions. These results are helpful to understand realistic input conditions of flow depth and velocity fields for different Fr regime that can be used to develop the tsunami fragility curves, which utilize the random combinations of h_max and flow velocity in the generation of fragility curves (Attary et al., 2017a; Attary et al., 2017b; Alam et al., submitted²). In addition, it is worth noting that the maximum flow depth and momentum flux typically do not occur at the same time (Park et al., 2013). Therefore, Fr at the maximum momentum flux and at the maximum flow depth can be expected to be slightly different. However, plots similar to Figure 11 (not shown in the interest of brevity) in which the flow depth h at the instant of the maximum momentum flux (M_F)_max or the instantaneous momentum flux is plotted against the corresponding maximum inundation flow depth lead to similar conclusions that the flow is generally subcritical for the maximum IMs.

Figure 12 shows the spatial distributions of the earthquake and tsunami IMs at Seaside, Oregon, for the AEP = 0.001 (often referred to as the “1,000-year event”) for the full-rupture scenario at CSZ. Figures 12A–D presents the spatial distribution of PGA, S_a (T₁ = 0.3 s), h_max, and (M_F)_max respectively, and the dotted contour lines in each panel show the maximum inundation limits (h_max = 0.3 m). Two distinct regions of earthquake IMs (PGA and S_a) are observed for the study area depending on the soil site class assumed in Figure 5C and within each region distributions of IMs are generally uniform. In the case of tsunami IMs (h_max and (M_F)_max), irregular distributions are observed depending on the bathymetry conditions. However, both the IMs generally decreases from the shore toward inland (i.e., in the positive x-direction).

The range of PGA and S_a(T₁ = 0.3 s) over the study region is 0.47–0.50 g and 1.03–1.08 g, respectively, while the h_max and (M_F)_max range from 0 to 12 m and 0 to 120 m³/s², respectively for the 1,000-year event. The values of PGA and S_a are uniformly distributed over the entire study area, while the area affected by h_max and (M_F)_max are limited to the maximum inundation limits. To understand the granularity of both IMs of the earthquake and tsunami at the same event, the spatial mean and deviation of IMs conditional on varying unit block size is computed. This analysis is performed on the rectangle region, shown in Figure 5A, which is a 1,990 m length and 3,000 m wide rectangle near the center of the City of Seaside. The smallest unit block size, designated as reference or unit block size here forth, is equal to the Cartesian mesh grid size (dx = dy = 24 m) considered for tsunami inundation modeling in the C-grid. The ratio of block size rB=Block sizereference mesh size is increased systematically, keeping the same block shape and the mean and SD of the IMs are estimated as a function of the changes in the block size. The mean of normalized IM is given by (18)IM′=∑j=1Nb∑i=1Ngimi,j′NbNg where N_b is the number of blocks depending on r_B conditions, and N_g is the number of grid points in a jth block (N_g = r_B²). The imij′ is the normalized IM of grid points of the jth block. The IMs in the jth block are normalized by the mean of the IMs for the jth blocks, and this calculation is performed for eight r_B conditions: r_B = 1, 2, 3, 5, 8, 13, 18, and 25. The corresponding block size and number of grid points per each block (N_g) are 24 m (1), 48 m (4), 72 m (9), 120 m (25), 192 m (64), 312 m (169), 432 m (324), and 600 m (625). For example, the total number of grid points in the study region is 10,625. For the case of r_B = 5, there are 425 block subsets over the study region, and each block is composed of N_g = 25 mesh grid points. The mean of IM at each of the 425 blocks is first computed, and then, each IM is normalized by the calculated mean.

Figure 13 shows the PDF of ln(h_max′) over the study region for AEP = 0.001 (1,000-year event). Figures 13A–D show the PDFs at r_B = 2, 5, 13, and 25, respectively. The red line in each panel shows the natural log-normal fitting curve. When ln(h_max′) is equal to zero, the mean of block subsets match exactly with the h_max of the unit grids, and each negative or positive value indicates an overestimation or underestimation of the block means. In the case of h_max, the PDF shape becomes wider and shows a larger deviation as the block size increases. This process was extended to the three other IMs: PGA and S_a and (M_F)_max.

Figure 14 summarizes the results of the four granularity tests for AEP = 0.001 (1,000-year event) by plotting the mean of the four normalized IMs with 90% confidence intervals for each unit block size. Both PGA and S_a (T₁ = 0.3 s) show almost insensitivity to the block sizes. All of the normalized means are essentially zero and the confidence intervals are small for all block sizes. However, both h_max and (M_F)_max results show significant variation with the 90% confidence intervals increasing as the block size increases. The largest deviation is found at (M_F)_max, and relatively smaller deviation if found h_max. The range of confidence interval of both h_max and (M_F)_max increases sharply with r_B = 13 (block size = 312 m). These results highlight the different sensitivity between two earthquake and tsunami IMs to the block size used to aggregate the IM results. In case of the PGA and S_a, the information of the site and location are less significant for determining each IM while h_max and (M_F)_max require detailed information of the site to minimize the uncertainty.

Lastly, the spatial distribution of joint IMs is analyzed next. Only h_max and (M_F)_max are considered herein since significant variations of tsunami IMs across the study area is observed in Figure 12. For the joint distribution computations uniform bins at 0.2 m interval are used for h_max, while a 4 m³/s² interval is used for (M_F)_max. The probability (%) of the joint distribution is computed by counting number of grids which involved in the joint h_max and (M_F)_max bin, for the whole study region, at the three AEP conditions, i.e., the 500, 1,000, and 2,500-year events, respectively. The results of spatial distributions of joint probability of h_max and (M_F)_max are shown in Figures 15A–C for the 500, 1,000, and 2,500-year events, respectively, with four Froude numbers plotted on each figure in a manner similar to Figure 11.

In the case of the 500-year event (Figure 15A), more than 45% of the joint h_max and (M_F)_max is distributed within 0.1 < Fr ≤ 0.5. The larger (M_F)_max is observed within 0.5 < Fr ≤ 1.0. The joint peak of h_max and (M_F)_max is located at 1.2 < h_max ≤ 1.4 m, and 0.0 < (M_F)_max ≤ 4.0 m³/s². Two distinct narrow banded regions of h_max are observed near h_max = 3.8 and h_max = 5.8 m, which correspond to different ranges of (M_F)_max even at the similar flow depth conditions. In the case of the 1,000-year event (Figure 15B), more than 60% of the joint h_max and (M_F)_max is distributed within 0.1 < Fr ≤ 0.5. Large (M_F)_max values are observed at both 0.1 < Fr ≤ 0.5 and 0.5 < Fr ≤ 1.0. The joint peak is located at 1.6 ≤ h_max < 1.8 m, and 0.0 < (M_F)_max ≤ 4.0 m³/s² which is a slightly higher h_max condition than the 500-year event. Several isolated islands of joint distributions are observed that have relatively higher (M_F)_max values. In the case of the 2,500-year event (Figure 15C), more than 68% of the joint h_max and (M_F)_max is located within 0.1 < Fr ≤ 0.5. Large (M_F)_max values are located primarily within 0.1 < Fr ≤ 0.5. The joint peak is located at 3.4 < h_max ≤ 3.6 m, and 16.0 < (M_F)_max ≤ 20.0 m³/s², which are significantly larger than the 500 and 1,000-year events. A few isolated points are observed for 2500-year event similar to that observed in 1,000-year event. Overall, the joint h_max and (M_F)_max is distributed within 0.1 ≤ Fr < 1.0 and are typically in the range 0.1 < Fr ≤ 0.5. It is also observed that each event has a different joint peak location. Further research is needed to understand the extent to which these are site-specific results, how this work can be generalized for other coastal archetypes (e.g., embayments), and how these results can be parameterized based on conditions offshore of the study area.