A regionalized partially nonergodic ground-motion model for subduction earthquakes using the NGA-Sub database

In this study, we derived a regionalized partially nonergodic empirical ground-motion model (GMM) for subduction interface and intraslab earthquakes using an extensive global database compiled as part of the NGA-Subduction project. The model can be used to estimate peak ground acceleration (PGA), peak ground velocity (PGV), and ordinates of 5%-damped pseudo-spectral acceleration (PSA) at periods ranging from 0.01 to 10 s for M ≥ 5.0, M≤ 8.5 for intraslab events, M≤ 9.5 for interface events, ZTOR≤ 50 km for interface events, ZTOR≤ 200 km for intraslab events, 10 ≤RRUP≤ 800 km, and 100 ≤VS30≤ 1000 m/s. Besides a global version of the model, the GMM accounts for regional differences in the overall amplitude (constant), anelastic attenuation, linear site response, and basin response for seven subduction-zone regions: Alaska (AK), Central America and Mexico (CAM), Cascadia (CASC), Japan (JP), New Zealand (NZ), South America (SA), and Taiwan (TW). The functional form of the model is structured such that the breakpoint magnitude, the magnitude at which the magnitude-scaling rate (MSR) transitions from a steeper to a shallower slope, is an adjustable parameter in the model. This makes it possible to take epistemic uncertainty in this parameter into account or adjust it based on other empirical or physical information, such as when the model is applied to a subduction zone not considered in the GMM. Besides the traditional mixed-effects aleatory between-event standard deviations and within-event standard deviations, within-model epistemic standard deviations in the median prediction for each region is quantified from a posterior distribution of model coefficients, standard deviations, and coefficient correlations using a Bayesian regression approach. Our full 800-sample posterior distribution can be used to account for epistemic uncertainty in the model coefficients, standard deviations, and predicted values. We also provide a simplified epistemic model using magnitude- and distance-dependent within-model standard deviations that can be used to facilitate the inclusion of within-model epistemic uncertainty directly in a probabilistic seismic hazard analysis. The within-model standard deviations can also be used to scale the GMM using a backbone modeling approach.


Introduction
Subduction zones can produce very large megathrust earthquakes and, thus, pose a significant seismic risk in many regions of the world. Since ground-motion models (GMMs) are an integral part of probabilistic seismic hazard analysis (PSHA), the development of GMMs for subduction earthquakes is of major importance to accurately quantify the seismic hazard from such events. This has led to the development of several empirical and simulation-based GMMs over the years for both interface and intraslab subduction earthquakes. The models based on stochastic or physics-based ground-motion simulations are typically developed for a specific application to a localized region, such as the CASC in the Pacific Northwest of the United States (Atkinson and Macias, 2009;Gregor et al., 2002), whereas many empirical models are developed for more global applications (Abrahamson et al., 2016;Abrahamson and Gulerce, 2022;Atkinson andBoore, 2003, 2008;Parker et al., 2022) or for regions where there are an abundance of recordings, such as Chile (Fayaz et al., 2023;Idini et al., 2017;Macedo and Liu, 2022;Montalva et al., 2017Montalva et al., , 2022, Ecuador (Arteta et al., 2021), Greece (Kkallas et al., 2018), JP Ghofrani and Atkinson, 2014;Hassani and Atkinson, 2021;Morikawa and Fujiwara, 2013;Si et al., 2022;Zhao et al., 2006Zhao et al., , 2016aZhao et al., , 2016b, and TW Phung et al., 2020).
In this article, we summarize a regionalized partially nonergodic empirical GMM that was developed by Kuehn et al. (2020a) as part of the Next Generation Attenuation (NGA) research project on subduction earthquakes , hereafter referred to as NGA-Sub. We refer to this model as KBCG20 throughout the article. KBCG20 can be used to predict the RotD50 horizontal component (Boore, 2010) of peak ground acceleration (PGA), peak ground velocity (PGV), and ordinates of 5%-damped pseudo-spectral acceleration (PSA) at 21 periods ranging from 0.01 to 10 s. Details about the model and modeling approach is available in the study by Kuehn et al. (2020a), although we note that some of the coefficients from that study have since been modified. We regionalized the GMM for seven subduction-zone regions with special attention to CASC, because of the limited number of earthquakes from this latter region and its importance to seismic hazards in the Pacific Northwest of the United States (Petersen et al., 2020). In particular, KBCG20 was developed as a partially nonergodic model similar to the approach described in Stafford (2014), which means that some terms in the model vary by region while other terms are shared among the regions. The model parameters are estimated using Bayesian inference (Gelman et al., 2013;Spiegelhalter and Rice, 2009), which allows the inclusion of prior information in a probabilistic way. We also place strong emphasis on the estimation of within-model epistemic uncertainty, which is estimated from a posterior distribution of model parameters derived from the Bayesian regression.

Database
We used a subset of the NGA-Sub database based on a set of selection criteria designed to include only those data that are deemed to be accurate, reliable, and usable for GMM development (Kuehn et al., 2020a). The predictor variables included in KBCG20 are summarized in Table 1. Our data selection and exclusion criteria are listed below.
Inclusion criteria: 1. Interface (F S = 0); intraslab (F S = 1), or the lower part of a double seismicity intraslab zone (F S = 5). The latter two are combined and referred to simply as intraslab or slab events (F S = 1). Contreras et al. (2022) provides details on event classification in the NGA-Sub database. 2. Magnitudes of M . 4.0 which, although may not be of engineering interest, includes smaller magnitudes and an increased number of recordings that helps to constrain many of the model parameters in the regression analyses. 3. Ratios of the largest to smallest rupture distance larger than 2, which ensures that each event has a reasonable distance range that allows a reliable estimate of event terms and attenuation parameters. 4. Recordings with non-missing values of moment magnitude (M), rupture distance (R RUP ), time-averaged shear-wave velocity in the upper 30 m of the site (V S30 ), and depth to the top of the rupture surface (Z TOR ). 5. Recordings with a minimum rupture distance of R RUP \ 800 km or R RUP \ R MAX if smaller, where R MAX is the minimum distance to non-triggered recordings for each event, which arises when an instrument triggers only above a certain amplitude threshold (see Contreras et al., 2022).
Exclusion criteria: 6. Predictor variables other than basin depths with missing values. 7. Stations with an instrument depth . 2 m to avoid depletion of high-frequencies due to embedment effects. 8. Recordings with PGA . 10 g, which likely represent incorrect instrument gains or other errors. 9. Recordings with a multiple event flag equal to 1 representing recordings with more than one recorded earthquake. F S Flag indicating interface event (F S = 0) or intraslab event (F S = 1) F X Flag indicating no arc-crossing path (F X = 0) or arc-crossing path (F X = 1) M Moment magnitude PGA 1100 Median predicted value of PGA on rock with V S30 = 1100 m/s R 1 Distance traveled within Subregion 1 (Backarc) (km) R 2 Distance traveled within Subregion 2 (Forearc; Global, Japan Trench) (km) R 3 Distance traveled within Subregion 3 (Forearc; Nankai Trough) (km) R RUP Closest distance to rupture plane (km) V S30 Time-averaged shear-wave velocity in upper 30 m of a site (m/s) Z 1.0 Depth to 1.0 km/s shear-wave velocity horizon below the site (km) Z 2.5 Depth to 2.5 km/s shear-wave velocity horizon below the site (km) Z xx Generic representation of either Z 1.0 or Z 2.5 Z TOR Depth to top of rupture surface (km) 10. Recordings with a visual quality flag equal to 2 representing a late S-trigger or 9 representing non-useable data. 11. Recordings with a GMX site classification first letter code of N, Z, or F representing non-free-field stations.
We do not discard recordings from stations with no estimate of Z 1.0 or Z 2.5 . These socalled basin or sediment depths are only available for a limited number of stations and regions. Discarding these stations would severely reduce the number of usable data. Furthermore, where they are available, including them in the regression was found to significantly reduce the within-event aleatory variability at long periods.
For CASC, two events were classified as interface earthquakes in the database. Located off the central Oregon coast, they have magnitudes of M = 4.7 and M = 4.9. After discussions with the NGA-Sub project researchers and the USGS National Seismic Hazard Model Program, and considering the questionable classification of these events as interface events, these earthquakes were excluded from our database.
We also excluded the 4 October 1994 M 8.28 Kuril event because of its location on the outer rise of the Kurile Subduction Zone and its thrust-oblique focal mechanism (Tanioka et al., 1995), which is inconsistent with other Benioff-zone normal-faulting intraslab events included in the database. Based on an initial regression, we also removed all recordings with absolute residuals that deviated more than four times the within-event standard deviation from the median prediction at eight selected periods of engineering interest. This latter criterion was a NGA-Sub project team decision to remove those events and recordings with potential errors in instrument properties or other metadata. After all of the aforementioned criteria were applied, we selected those events that had at least five recordings to allow the calculation of reliable event terms. No such criterion was applied to the number of recordings per station.
The selected database has a maximum of 16,045 three-component recordings from 238 events and 3769 recording stations. There are 6864 recordings from 113 interface events and 9181 recordings from 125 intraslab events. At long periods, the number of recordings is reduced due to a smaller useable bandwidth. Table 2 lists the number of recordings, number of events, number of stations, and region abbreviation for each of the subduction regions for the subset of NGA-Sub data selected for analysis. JP dominates the data set with over 50% of the recordings. TW and SA also have a significant number of events but fewer recordings per event than JP. The number of available recordings for CASC, CAM, and NZ is relatively sparse. The magnitude-distance distribution of the data is shown in Figure 1.

Ground-motion model
The GMM is developed as a Bayesian multi-level/hierarchical model (Gelman and Hill, 2006). This means that we recognize that there are different levels in the data that can be exploited during the model-building process. Each level consists of data that can be grouped according to some criterion, such as by event or geographic region. The levels form the following hierarchy from broadest to narrowest: global ! regional ! event ! recording. Within each level, some parameters are shared within groups, such as the event term for all recordings from the same event, while others are taken from higher levels. The inclusion of the regional, event, and recording levels allows us to relax the ergodic assumption that all recordings come from the same population by associating them with groups or hierarchical levels (Anderson and Brune, 1999), leading to what is referred to as a partially nonergodic model (Stafford, 2014).

Base model
At the recording level, the model assumes that the target ground-motion variable ln Y is distributed according to a normal or Gaussian distribution given by its mean m and standard deviation f: where m is the sum of the natural logarithmic median base function f base (ũ,x) and an event term dB. The event terms have the following normal distribution with mean zero and standard deviation t: The mean prediction f base (ũ,x) depends on the model coefficientsũ (see Table 3) and the predictor variablesx (see Table 1) defined by the relationship: Anelastic attenuation within backarc region 1 after crossing a volcanic arc Yes u 6,x2 Anelastic attenuation within forearc region 2 after crossing a volcanic arc Yes u 6,x3 Anelastic attenuation within forearc region 3 after crossing a volcanic arc Yes u 6,1 Anelastic attenuation within backarc region 1 without crossing a volcanic arc Yes u 6,2 Anelastic attenuation within forearc region 2 without crossing a volcanic arc Yes u 6,3 Anelastic attenuation within forearc region 3 without crossing a volcanic arc Yes  Table 3 also identifies which of the coefficients are regionalized. In addition to the estimated coefficients, there are a few coefficients that are fixed (see Table 4). To ensure a physically meaningful spectrum, the predicted median PSA at short periods is not allowed to be smaller than PGA as given by the equation: where u ! is replaced with u ! (T ) to indicate that the coefficients are a function of period. The parenthetic T is dropped in the remainder of the article for simplicity.

Logistic hinge function
A bilinear logistic hinge function is used in several of the terms listed in Equation 4. It provides a smooth transition from one slope to another. It is defined by the equation: where x 0 is the slope breakpoint, b 0 is the slope for x \ x 0 , b 1 is the slope for x . x 0 , and d controls the smoothness of the transition between slope b 0 and b 1 with smaller values of d leading to a sharper transition. The bilinear logistic hinge function is used to model the magnitude-scaling term f mag (M, F S ) and the depth-scaling term f depth (Z TOR , F S ). Compared to a bilinear function that includes a non-continuous logical statement to model different slopes, the logistic hinge function is differentiable everywhere, making it more tractable in a nonlinear regression analysis.

Constant term
The constant term is given by the equation: which allows a different constant depending on whether the earthquake is an interface or an intraslab event.

Magnitude term
The magnitude-scaling term is modeled by the following bilinear logistic hinge function: The magnitude-scaling functional form is the same for interface and intraslab events, but the model coefficients and fixed parameters are different (see Tables 3 to 5). In particular, the breakpoint magnitude (M B ) and the MSR below the breakpoint magnitude (u 4,if and u 4,slab ) are different for interface and intraslab events. The MSR above the breakpoint magnitude (u 5 ) is assumed to be the same for both interface and intraslab events lacking information to the contrary. The smoothness of the transition between the two linear segments of the scaling relation is fixed at d M = 0.1. Preliminary trials to estimate this parameter demonstrated that it was not well constrained and that its posterior distribution was almost unchanged compared to its prior distribution. The selected value leads to a reasonably smooth transition between the two linear segments. The magnitude term has a value of zero at M ref = 6.0.
During an initial regression, it became apparent that at long periods (T . 1 s), the event terms for interface events at large magnitudes were biased low compared to the observations. As a result, an adjustment to M B,if for interface events at T . 1 s was modeled similar to that proposed by Abrahamson et al. (2016Abrahamson et al. ( , 2018. This adjustment is zero at T < 1 s and 20.4 at T ø 3 s with a log-linear interpolation in between (i.e. logarithmic in period and linear in the adjustment) as given by the equation: This adjustment is only applied to JP and SA interface events in the regression, since these regions are the only ones with large-enough magnitudes to be effected by this adjustment. However, we believe that there is no reason why it should not be applied to every subduction region in a forward prediction of large-magnitude ground motions.

Geometric attenuation term
The geometric attenuation term (i.e. geometrical spreading) is modeled by the equation: where the finite-fault or ''pseudo-depth'' term is defined by the equation: The regression coefficients u nft,1 and u nft,2 are given a relatively strong prior distribution because they are not well-constrained by the data.

Source-depth term
The source-depth term is given by a bilinear logistic hinge function with different slopes and depth breakpoints for interface and intraslab events as follows: where Z if,ref and Z slab,ref are reference depths (see Table 4) and d Z is the smoothness transition parameter between the two linear slopes. The bilinear DSR has breakpoints at Z B,if + dZ B,if and Z B,slab + dZ B,slab determined from the regression. The breakpoints were set to Z B,if = 30 km and Z B,slab = 80 km based on initial regressions. The perioddependent adjustment coefficients dZ B,if and dZ B,slab were determined by regression because the DSR breakpoints have no theoretical basis. The coefficients for the DSR up to the depth breakpoints are u 9,if and u 9,slab . Both the coefficients for the DSR above the break points were fixed at u 10 = 0, consistent with the study by Abrahamson et al. (2016Abrahamson et al. ( , 2018, because there are not many events with depths larger than these breakpoints and those depths that do exist do not show any systematic trend with ground-motion amplitude.

Site-response term
The site-response (site-amplification) term is parameterized as in the study by Campbell and Bozorgnia (2014): The linear site-scaling coefficient u 7 is a regression coefficient that varies by geographic region. The period-dependent parameters k 1 and k 2 , as well as the period-independent constants c and n, are fixed to the values of Campbell and Bozorgnia (2014). Hence, we make the implicit assumption that the nonlinear site amplification for subduction events is similar to that of crustal events, consistent with the assumptions of Abrahamson et al. (2016) and Montalva et al. (2017), although we note that Zhao et al.

Anelastic attenuation term
The anelastic attenuation term is defined by the equation: where F X = 1 when the travel path passes through a volcanic arc (i.e. the transition between the forearc and backarc attenuation subregions) and F X = 0 otherwise. The regression coefficient u 6,i represents the anelastic attenuation coefficient for attenuation subregion i before it crosses a volcanic arc. The regression coefficient u 6,xi represents the anelastic attenuation coefficient for attenuation subregion i after it crosses a volcanic arc. The distance R i is the distance traveled within each subregion, the sum of which equals R RUP . The coefficient u xc is the offset as the travel path crosses the volcanic arc and u 6,2 is the attenuation rate for those subduction regions that do not have different attenuation subregions.
The difference in how the subduction regions are modeled is based on inspection of residual plots and ground-motion scaling of trial regressions for each subregion that did not include forearc and backarc attenuation coefficients. Subregion index 1 indicates a travel path within the backarc of JP, CAM, or SA; Subregion 2 indicates a travel path within the forearc of CAM, SA, or northeastern JP; and Subregion 3 indicates a travel path within the forearc of southeastern JP. The northeastern JP forearc sits atop the JP Trench megathrust and the southeastern JP forearc sits atop the Nankai Trough megathrust. Figure 2 shows the attenuation subregions for JP, SA, and CAM.

Basin-response term
The effect of the deep crustal structure beneath the site is modeled with the following basin-response term: Seattle basin min u 11 + u 12 d ln Z , u 11, SEA ð Þ All other basins where Z xx is either Z 1.0 or Z 2.5 depending on the region, d ln Z = ln Z obs À ln Z ref (V S30 ), Z obs is the observed value of basin depth, and Z ref (V S30 ) is the reference value of basin depth for a given subduction region and value of V S30 . Z ref (V S30 ) is defined by the equation: where the regional coefficients u z1 through u z4 are different for CASC, JP, NZ, and TW. Basin depths are not available for all subduction regions or for all recording sites within a region. Although both Z 1.0 and Z 2.5 are available for CASC, the values of Z 1.0 are typically associated with the depth to the relatively shallow high shear-wave velocity horizon at the quaternary-tertiary boundary (Stephenson et al., 2019) and is not a good representation of the deeper basin structure that amplifies intermediate-to-long period ground motion (Chang et al., 2014;Frankel et al., 2009). Both basin depths are also available for JP, but we selected Z 2.5 to be consistent with the depth parameter used by Si et al. (2022) in their basin-response term and the relatively deep basin depths that Morikawa and Fujiwara (2013) found gave the smallest standard deviations at T . 0.3 s when a basinresponse term was included in their GMM. Only Z 1.0 is available for NZ and TW and no basin depths are available for AK, CAM, and SA. For regions where basin depths are available, the basin-response term is defined as a function of the difference between the observed and a reference basin depth defined by the regional relationship between basin depth and site velocity given by Equation 16. This is done to minimize the impact of the correlation between basin depth and site velocity in the regression.  Figure 4 shows a plot of V S30 versus Z 2.5 for CASC. For sites in CASC, we distinguish between non-basin sites, for which the basin-response term is zero, and sites located within the Everett, Georgia, North Portland, Seattle, Tacoma, and Tualatin basins (Kuehn et al., 2020a). In this latter figure, the sites inside the Seattle basin are found to have basin depths of Z 2.5 ' 7000 m independent of the value of V S30 . If these sites were included in the regression for Equation 16, it would bias the slope of the ln Z ref and d ln Z terms. To avoid this bias, we use a constant basin-amplification term (u 11,SEA ) for sites located inside the Seattle basin.

Geographic regionalization
As indicated in previous sections of this article, we account for regional differences in the constants u 1,if and u 1,slab ; the linear site-response coefficient u 7 ; the anelastic attenuation coefficients u 6,x1 , u 6,x2 , u 6,x3 , u 6,1 , u 6,2 , and u 6,3 ; and the basin-response coefficients u z1 , u z2 , u z3 , and u z4 (see Table 2). Regional differences in ground-motion scaling are expected as a  result of differences in tectonic and geological conditions. In particular, regional differences in the anelastic attenuation coefficients are expected to be related to differences in the seismological quality factor Q, and regional differences in the linear site-response and basin-response terms are expected to be related to differences in the average shear-wave velocity profile. In addition, differences in the average stress drop can translate to differences in the constants at short periods. Although basin response is also regionalized, because this term is fit to the within-event residuals after the GMM is developed, it is not part of the regionalization included in the regression. There are only enough data to estimate regional basin-response terms and their epistemic standard deviations for CASC, JP, NZ, and TW and not for a global model plus these regional adjustments.
In the Bayesian regression, the regionalized coefficients are modeled as regional random effects (Stafford, 2014) according to the following equation: where m u is the global coefficient and d u ! reg is the regional adjustment factor. The subscript reg has been dropped from the regional coefficients throughout this article for simplicity in nomenclature. One can refer to Table 3 for those coefficients that are regionalized. The regional adjustments are joint normally distributed with mean zero and marginal standard deviation c u . For further details, see Chapter 4 of the study by Kuehn et al. (2020a) and the Supplemental material to the article.
Past studies have found that the small-magnitude short-period CASC intraslab ground motions are significantly lower compared to intraslab events from other regions (Abrahamson et al., 2016(Abrahamson et al., , 2018Atkinson, 1997;Atkinson and Boore, 2003). Although we have confirmed this result, we also found that the ground motions from the two largest intraslab events in CASC are comparable to, although somewhat lower than, our global model. We did not want to bias the regional CASC constant u 1,slab toward low values based on ground motion data from these smaller earthquakes yet there are too few of them to define a CASC-specific MSR. Therefore, we used only the two largest events, the M 6.8 Nisqually and the M 6.55 Ferndale earthquakes, in the random effects regression to derive the regional CASC constant. For the other 10 CASC intraslab events, we centered the ground motions from each event using an event-specific constant so they could be used to estimate regional linear site-response and anelastic-attenuation coefficients without impacting the regional intraslab constant.
We also found that the within-event standard deviation f was smaller for those regions with a basin-response term. This was not the case for the between-event standard deviation t that was not regionalized. We also evaluated whether the standard deviations were dependent on magnitude, distance, and subduction type (i.e. interface and intraslab), but no significant biases or trends were found due to limited data and a significant overlap of their posterior distributions. However, we did find and included a strong regional dependence of epistemic uncertainty. Additional details are available in the study by Kuehn et al. (2020a).

Regionalization of breakpoint magnitude
It is clear from modern empirical subduction interface GMMs (Abrahamson et al., 2016(Abrahamson et al., , 2018Campbell et al., 2022;Morikawa and Fujiwara, 2013;Zhao et al., 2016b) and other literature reviews (Campbell, 2020;Stewart et al., 2013) that there is a magnitude at which the MSR of ground-motions from subduction interface earthquakes must become smaller. In KBCG20, this magnitude is referred to as breakpoint magnitude (M B ) after the study by Campbell (2020). The values of M B depend on the geometry, seismogenic characteristics, and age of the subduction zone. We base the values of the breakpoint magnitude used in KBCG20 on the studies of Ji and Archuleta (2018) for intraslab events (M B,slab ) and Campbell (2020) for interface events (M B,if ). For more information, we refer the reader to these studies and Section 8.2 of the study by Kuehn et al. (2020a). The suggested values of M B for the seven subduction zones and related subregions within these subduction zones are listed in Table 5.

Model estimation
The model parameters estimated in the Bayesian regression are the global coefficients u ! , the aleatory standard deviations f and t (corresponding to within-event and betweenevent variability, respectively), the event terms dB, and the regional adjustment terms d u ! together with their standard deviations and correlation coefficients. The model parameters are estimated using Bayesian inference. In Bayesian inference, the posterior distribution of a model parameter is proportional to its prior distribution times its likelihood. Kruschke (2015) and Spiegelhalter and Rice (2009)

Posterior distributions
The posterior distributions of the model parameters are estimated from Markov Chain Monte Carlo (MCMC) sampling using the Bayesian computer platform Stan (Carpenter et al., 2017). These posterior distributions are then used to calculate summary statistics such as the mean, median, fractiles and correlation coefficients of the coefficients, predicted values, and standard deviations. In total, we generate 800 samples from the posterior distribution. In Equation 1, the recording level of the GMM is modeled as a normal distribution with mean m and within-event standard deviation f, which means that the ground-motion parameter of interest Y is distributed according to a lognormal distribution. In a forward application of the model, this is how the model should be used even though the fit is done assuming a Student's t-distribution as discussed next.
During exploratory regressions, we observed several apparent outlier recordings. Including such outliers can lead to a biased estimate of the median prediction and an increased estimate of the standard deviations. To mitigate this issue, we performed a Bayesian robust regression. Instead of modeling the within-event variability with a normal distribution as is traditionally the case, a Student's t-distribution is used (Gelman and Hill, 2006;Kruschke, 2015). Compared to the normal distribution, the Student's t-distribution has wider tails, which makes it less susceptible to outliers. For more details, see Chapter 5 in the study by Kuehn et al. (2020a) and the Supplemental material to this article.
To ensure a smooth predicted response spectrum, we smoothed the means of the model coefficients and aleatory standard deviations (i.e. the mean of the 800 posterior samples for each model parameter) using a Gaussian process (GP) regression (Rasmussen and Williams, 2006) with a squared exponential covariance function. For more details see Section 4.3 in the study by Kuehn et al. (2020a) and the Supplemental material to the article. The regional standard deviations are not smoothed because they do not affect the predicted median spectrum. For the same reason, the parameters of the correlation matrix of the regional coefficients are not smoothed. We re-centered the posterior distribution of each coefficient by subtracting its smoothed value from its distribution mean and adding this difference to each posterior value. Not only does this retain the range of values in the original posterior distribution, it also retains the correlation between the posterior distributions of the coefficients.

Prior distributions
We classify the prior distributions of the regression coefficients and the other model parameters as being a mix of weakly informative and informative depending on the available data. For those regression coefficients that are not well constrained by data, it is important to define informative prior distributions using other types of information (e.g. physical constraints or physics-based ground-motion simulations) for the results to be meaningful. For other model parameters, a weakly informative prior distribution is used to downweight ranges of the parameter space that can be ruled out a priori. Although the term weakly informative is not well defined, it is usually interpreted to mean a probability distribution that is not very wide. Stafford (2019) suggests defining the prior distributions based on previous GMMs, which can help to stabilize regression results and lead to more reasonable predicted values (Kowsari et al., 2020;Kuehn and Scherbaum, 2016). However, because most GMMs do not report the full distributions of their coefficients, this information is generally not available. Furthermore, there is considerable overlap between the NGA-Sub database and data sets used in previous GMMs, which means that prior distributions based on published models might double-count some data.
We used informative prior distributions for the MSR coefficients u 4,if and u 4,slab for small events and u 5 for large events, as well as for the coefficients u nft,1 and u nft,2 of the near-fault saturation term. The prior distributions for the other regression coefficients are described in Section 4.1.8 of the study by Kuehn et al. (2020a) and in the Supplemental material to the article. The prior distribution of u 5 needs to be informative because there are not many events to empirically constrain it. Therefore, we base its prior distribution on results of the simulation-based GMMs of Atkinson and Macias (2009) and Gregor et al. (2002). The effective large-magnitude MSRs of these two models are shown in Figure 5. For comparison, we also plot the MSR estimated by Ghofrani and Atkinson (2014) from M . 7.0 empirical interface earthquakes in JP, although we note these authors found that their MSR at long periods was not statistically significant. To enforce magnitude saturation, but not allow oversaturation, of ground motions at large magnitudes and short distances, we constrain u 5 to be positive and smaller than the value of u 4,if for small-event MSRs. The standard deviation of the prior distribution is set to 0.2, corresponding to the range of MSRs across different periods, resulting in the following prior distribution for u 5 : where T(a,b) indicates a truncated normal distribution with lower limit a and upper limit b.
The prior distributions of the MSR at magnitudes below the breakpoint magnitude are defined by the equations: These prior distributions are informed by the physics-based simulations for intraslab events by Ji and Archuleta (2018). We performed a simple regression on the simulated ground motions and found the estimated MSR to be approximately one across all periods. Since the simulated data range is rather limited, we used the regression results as guidance to impart a wider standard deviation on the prior distribution based on the regression. The GMMs of Gregor et al. (2002) and Atkinson and Macias (2009) are valid for M . 7.5 and cannot be used to determine a prior estimate of the MSR at small magnitudes. However, we note that Campbell et al. (2022) empirically found an MSR of approximately 1 at magnitudes below the breakpoint magnitude for PGA using data from Japanese interface events.
Another parameter that needs to be constrained with an informative prior distribution is the near-fault or ''pseudo-depth'' term h(M) that controls the near-fault distance scaling of the GMM. This term is difficult to constrain empirically because of the limited number of recordings at short distances. We set the prior distributions of the two coefficients u nft,1 and u nft,2 that define this term based on the finite-fault terms in the studies by Abrahamson et al. (2016) and Parker et al. (2022). The latter study uses ground-motion simulations from the finite-fault stochastic simulation program EXSIM (Motazedian and Atkinson, 2005) to constrain the magnitude-scaling of h(M). To set the prior distribution for the finite-fault term in KBCG20, we require that the average of the prior distributions at M 6.0 and M 9.  The prior distributions for the regional standard deviations are exponentially distributed. The rate parameters of the exponential distributions are chosen such that large deviations from the global model are penalized, consistent with that proposed by Simpson et al. (2017). This approach is akin to Occam's razor, in which the simpler model is the global model and deviations from the global model are only modeled if there is a strong regional signal in the data. Since the number of regions is small, the estimate of the regional standard deviations can be quite noisy, which means it is important to set reasonable prior distributions for these parameters (Kuehn et al., 2020a).

Regional coefficients
All of the mean model coefficients and their posterior distributions are provided in the Supplemental material to the article. Herein, we briefly summarize the results of the Bayesian regression analyses with respect to regional differences, since this is one of the unique features of KBCG20 compared to previous subduction GMMs. A more comprehensive evaluation of the other coefficients is available in Section 4.2.3 of the study by Kuehn et al. (2020a) and the Supplemental material to the article. Figure 7 shows the values of the regional adjustment values of the interface constant (du 1,if ), the intraslab constant (du 1,slab ), the linear site-response coefficient (du 7 ), and the anelastic-attenuation coefficient (du 6,2 ). We only show du 6,2 because this coefficient is used in the anelasticattenuation terms of all regions as indicated in Equation 14.
Figure 7 indicates that regional differences in the coefficients can be relatively large, especially at short periods. Differences at long periods are mainly restricted to linear siteresponse scaling. As seen in Figure 7a, the regional adjustment for the CASC interface constant is slightly negative even though there are no CASC interface events in the database. The reason for this is that the regional adjustment coefficients are modeled as correlated between interface and intraslab events and the regional adjustment of the intraslab constant is negative. This leads to a negative interface constant adjustment for CASC.

Adjustment to AK regional coefficients
The database for the NGA-Sub project was finalized on 22 April 2019 (Bozorgnia and Stewart, 2020). The GMM presented in the study by Kuehn et al. (2020a) was developed using this database. Some changes were made to the database between this date and its final publication that we reviewed to determine which aspects of our model required modification. The modifications included values of V S30 for some of the recording stations in JP and SA and values of the distance metrics in AK . Preliminary analyses showed that changes in the predicted values due to differences in V S30 were small and could be neglected, but that changes in the predicted values due to differences in the distance metrics in AK were non-negligible. Therefore, we decided to refit the Kuehn et al. (2020a) model based on the updated AK distances. To avoid computational complexities, we only refit the AK regional adjustment coefficients to avoid having to refit the entire GMM. In the refitting process, we fixed the global coefficients to their estimated mean values and estimated the new regionalized AK adjustment constant, site-response coefficient, and anelastic-attenuation coefficient using the revised database. The outcome of this process was a set of new posterior samples of these regional coefficients. We then calculated the mean of the posterior samples and re-centered the posterior distribution obtained from the first fit to be centered around the mean value obtained from the refitting process. In this way, we keep the original posterior range and correlations. The resulting value of the kth posterior sample for coefficient i was obtained from the equation: where i is either the adjustment constant for interface events, the adjustment constant for intraslab events, the anelastic-attenuation coefficient, or the site-response coefficient, k is an index representing the 800 posterior samples, m i,origfit is the mean regional coefficient from the original fit, and m i,refit is the mean regional coefficient from the updated fit. This refitting process did not change the coefficients for the other subduction regions or the global model.

Aleatory variability
The total aleatory variance is calculated from the equation s 2 = f 2 + t 2 . We do not account for any systematic difference in site terms because the limited azimuths of many of the recordings means that anelastic attenuation could be inadvertently mapped into site effects which would bias the standard deviations. We developed the model under the assumption that the values of t and f are the same for all subduction regions, except for the smaller value of f for long-period PSA for those regions with basin-response terms. Due to a very different number of events and recordings in the different regions (see Table 2), this is a reasonable assumption as it would be difficult to obtain a reliable estimate of betweenevent and within-event variability for regions such as CASC or CAM where data are sparse. Figure 8 shows a plot versus period of both the original and smoothed mean values from the posterior distribution of the within-event standard deviations f, the betweenevent standard deviations t, and the total standard deviations s. The uncertainty in the standard deviations is shown as the 5% and 95% fractiles of the posterior distribution.

Model evaluation
As part of the evaluation process, we plot residuals for a limited number of groundmotion parameters and predictor variables as an example of the model validation we performed. Other residual plots are available in the study by Kuehn et al. (2020a) and at https://github.com/nikuehn/KBCG20. Figure 9 plots event terms (between-event residuals) against magnitude for T = 0.2 s and T = 3 s. Figure 10 plots within-event residuals for the same two periods against R RUP , indicating that attenuation is not dependent on subduction type. Overall, the residuals do not show any significant biases or trends with respect to the predictor variables that would indicate an issue with the GMM.    In both cases, R RUP = 100 km, V S30 = 400 m/s, basin depth (Z 1.0 or Z 2.5 ) is estimated from V S30 , and M B is taken from Table 5.

Example model predictions
cluttered. This figure shows that the regional differences in predicted values are larger at short periods than at long periods, which is consistent with the larger differences in the regional coefficients at shorter periods seen in Figure 7. The regional differences also become more pronounced at larger magnitudes. This is due to differences in breakpoint magnitude (see Table 5). Subduction zones with a larger value of M B have the potential to generate larger ground motions at large magnitudes. This is also seen in Figure 12 that shows the MSR of PGA for the different subduction regions and, if modeled, values of M B for different subregions within these regions. It is impractical to compare our median predicted values and standard deviations for the global and seven regional GMMs with existing subduction models in this article. Instead, the reader is referred to comparisons in the studies by Gregor et al. (2022) and Kuehn et al. (2020a).

Epistemic uncertainty
The model parameters are estimated from a finite set of data. As a result, they inherently contain epistemic uncertainty, also referred to as within-model uncertainty. This uncertainty translates into uncertainty in the median predictions. It is important to be able to estimate this uncertainty because for some event scenarios the between-model uncertainty from a logic-tree of alternative GMMs can be small compared to within-model uncertainty. Examples of within-model uncertainty models are given by Al Atik and Youngs (2014), Lanzano et al. (2019), and Kotha et al. (2020).
In a Bayesian regression, epistemic uncertainty is quantified by a posterior distribution of model parameters. In KBCG20, the posterior distribution of each regression parameter consists of 800 samples using the MCMC methodology. Thus, there are 800 sets of regression coefficients, aleatory standard deviations, and coefficient correlations, which can be used to create 800 samples of the median predicted values and standard deviations of a given scenario event from which the mean, standard deviation, and fractiles of the median prediction can be calculated. Figure 13a gives an example of the epistemic distribution of the predicted PSA spectrum for a particular scenario in JP. The results are in the form of a histogram of the median estimates of PSA at T = 0.01 s calculated from the 800 sets of coefficients from the posterior distribution. The median prediction is also shown for comparison. Figure 13b shows the mean MSR of PSA at T = 0.01 s for the scenario event calculated from the mean coefficients along with 10 randomly selected samples from the posterior distribution. Also shown are the 5% and 95% fractiles of the epistemic distribution for each magnitude. The 5% and 95% fractiles are calculated independently for each magnitude, meaning that different sets of sampled coefficients can contribute to the highest and lowest predicted values at each of the magnitudes. We highlight the predicted value from one particular sample in Figure 13b that has a steeper MSR than the mean to demonstrate the importance of within-model uncertainty in estimating ground motion. It also serves as a reminder that epistemic uncertainty permeates all aspects of the model (e.g. magnitude scaling, distance scaling, and site response).
It is especially important to account for epistemic uncertainty in partially nonergodic models so as not to underestimate the seismic hazard (Abrahamson et al., 2019). The regional adjustment terms are calculated from a relatively small number of recordings and are associated with more epistemic uncertainty than the global model in which all data are pooled. This is particularly important for regions with a small number of recordings such as CASC. This is demonstrated in Figure 14 in which we plot the epistemic standard deviation of the median predicted value of PGA (c m ) from the global model and the seven defined regions versus M, R RUP , V S30 , and T. This figure shows that epistemic uncertainty is relatively small for regions with a large number of recordings (e.g. JP and SA) and relatively large for regions with a small number of recordings (e.g. CASC and CAM). Figure 14d shows that there are larger regional differences between the values of m at short periods than at long periods. This is because the longer-period coefficients are better constrained in the regression (see Figure 7). We also show values of m for the global model, which are applicable when KBCG20 is applied to a new region. To calculate median predictions for the global model, we sample regional adjustment terms from their joint distribution and add the sampled adjustment terms to the coefficients from the posterior distribution as described in the study by Gelman and Hill (2006). This shows that the CASC values of m are closer to the values of the global model for interface events than for intraslab events, consistent with the fact that there are no CASC interface events in the database. Table 6 illustrates the effect of regionalization of epistemic uncertainty on the predicted value of PGA. Note that, although the value of f in this table is the same for all regions, it is smaller at longer periods for those regions that include a basin-response term (CASC, JP, NZ, TW). Table 6 compares the values of f, t, and total aleatory standard deviation s and the regional epistemic standard deviations c m to the standard deviations from a nonregionalized model (i.e. a model that includes only event terms). The nonregionalized model uses the same functional form and Bayesian regression methodology as KBCG20 except that none of the coefficients are regionalized. For the nonregionalized model, the coefficients are determined using all recordings in the database, which leads to a very small value of epistemic uncertainty. However, this small value of epistemic uncertainty is offset by larger values of f and t. The total predictive variability, which sums the aleatory and epistemic variances and describes the full range of all possible ground motions, is calculated as s 2 pred = t 2 + f 2 + c 2 m . This shows that the total predicted variability of the ground motion estimated from the regionalized global model when applied to a new region is very close to the total variability of the nonregionalized model, except that the epistemic uncertainty is a larger component of the total variability. This is an example of the trade-off between epistemic uncertainty and aleatory variability. The value of aleatory variability is reduced in the regionalized model, but there is a penalty for this reduction in the form of greater uncertainty for regions with less data. However, this trade-off is beneficiary for regions with a relatively large amount of data, as indicated by the reduction in the value of s pred compared to the nonregionalized model.

Engineering application
We consider KBCG20 to be generally applicable for estimating ground motions for the following range of parameters (Kuehn et al., 2020a).
The user should be aware that the epistemic standard deviation can increase (possibly considerably) if one or more of the predictor variables are near or beyond the limit of their observed values for the region of interest.
The total value of the aleatory component of variance is s 2 = f 2 + t 2 , which is the value that should be used in the hazard integral to account for aleatory variability. The most comprehensive estimate of epistemic variance is obtained from the full posterior distribution from the Bayesian regression and should be modeled using a logic tree. However, if one is only interested in the mean hazard and not its fractiles, but still wants to include epistemic uncertainty in the hazard analysis, the value of total predictive variability s 2 pred = f 2 + t 2 + c 2 m can be used in the hazard integral, thus avoiding the necessity of using a logic tree. To facilitate the use of within-model uncertainty in engineering analyses, and particularly in PSHA, we pre-calculated the epistemic within-model standard deviations (c m ) associated with the median predictions of interface and intraslab events for different magnitude and distance scenarios, holding the other predictor variables constant. These values are available at https://github.com/nikuehn/KBCG20 and can be used as a simplified version of epistemic uncertainty based on the observation that the value of c m depends only marginally on the values of Z TOR , V S30 , Z 1.0 , and Z 2.5 . They also can be used to estimate magnitude-and distance-dependent scale factors to use with the median values from KBCG20 using a scaled backbone approach similar to that proposed in the study by Atkinson et al. (2014). KBCG20 can be used to estimate ground motions for seven different subduction-zone regions. However, there are other subduction zones located throughout the world for which GMMs might be needed to perform a PSHA or deterministic seismic hazard analysis (e.g. the 79 subduction zones identified in the studies by Berryman et al., 2015 andCampbell, 2020). Some of these subduction zones might have limited or no ground-motion data. In that case the regionalized KBCG20 global model with its larger epistemic standard deviation should be used.
The epistemic uncertainty discussed above and in the Epistemic Uncertainty section of the article represents within-model uncertainty, since it only addresses uncertainty in the median estimates of the ground motion. Additional GMMs, such as those developed as part of the NGA-Sub project  or other published GMMs, should be used to account for between-model epistemic uncertainty similar to that recommended by Al Atik and Youngs (2014) for the NGA-West2 GMMs. The use of additional GMMs can address differences in basic assumptions, such as the functional form of the model and the database used to develop it, that are not accounted for in within-model uncertainty. The within-model uncertainty should be seen as a minimum representation of epistemic uncertainty.
An important aspect of KBCG20 is the use of regionalized breakpoint magnitudes (M B,if for interface events and M B,slab for intraslab events), which have a strong impact on the predicted values of ground motion at magnitudes exceeding these breakpoints. The values of breakpoint magnitude provided in this article are based on the geometry, seismogenic characteristics, and age of the subduction zone (Campbell, 2020;Ji and Archuleta, 2018) and should not be considered an inalienable parameter. KBCG20 is designed such that it is easy to replace breakpoint magnitudes in the GMM without impacting other aspects of the model, which facilitates the inclusion of epistemic uncertainty in these magnitudes for a given region or replacing them with more appropriate values for a new region (Campbell, 2020).
We did not partition the aleatory standard deviations of the within-event residuals from KBCG20 into a site term (f S2S ) and an event-and site-corrected (single-station) standard deviation (f SS ) because of the relatively narrow azimuthal range of the recordings and the potential that unmodeled path effects might be mapped into the site terms. Therefore, the aleatory variability is ergodic with respect to the site model. If a single-station standard deviation is desired, we suggest that studies available in the literature be used to estimate reduction factors to apply to f to estimate f SS . The majority of these single-station standard deviations are for crustal earthquakes, but some recent studies provide estimates of f and f SS for subduction earthquakes in Chile (Montalva et al., 2017(Montalva et al., , 2022, JP Hassani and Atkinson, 2021;Zhao et al., 2016aZhao et al., , 2016b, northern South America (Arteta et al., 2021), TW Phung et al., 2020), and globally Parker et al., 2022).

Cascadia Subduction Zone
The CASC is an example of why uncertainty in breakpoint magnitude should be considered. Campbell (2020) estimates a mean value for M B,if of 7.7 with epistemic 5% and 95% confidence limits of 7.3 and 8.2, respectively. These values are based on a relatively narrow preferred average seismogenic interface width of 68 km estimated by Berryman et al. (2015). Campbell (2020) poses the question of whether an earthquake that originates on the wider section of the CASC interface in the Puget Sound region will exhibit a larger breakpoint magnitude. To account for this uncertainty, a larger value of M B,if could be used for hypocenters located in the Puget Sound region and a smaller one for hypocenters originating along the narrower sections of the CASC interface. Alternatively, epistemic uncertainty could be modeled using a single set of breakpoint magnitudes regardless of hypocenter location. For example, Campbell (2020) suggests one possible set of values for M B,if could be 7.7 (5% confidence limit), 8.0 (mean), and 8.5 (95% confidence limit) with weights of 0.185, 0.630, and 0.185, respectively. We suggest using the mean of this distribution if only a single value of M B,if is used (see Table 5). Figure 15 shows the effect at large magnitudes of using the three-point discrete distribution of breakpoint magnitudes on the median prediction of PGA. Figure 15. Magnitude scaling for interface events in Cascadia (CASC) using the mean breakpoint magnitude and its uncertainty according to Campbell (2020).
Predictions are made for R RUP = 100 km, Z TOR = 10 km, V S30 = 400 m/s, and an estimate of Z 2.5 from V S30 .
Another potential issue with the prediction of ground motions for CASC interface events is the lack of interface earthquakes. As a result, we assumed that anelastic attenuation for CASC interface events was the same as for CASC intraslab earthquakes consistent with the CASC model of Abrahamson et al. (2018). On the contrary, Abrahamson and Gulerce (2022) assume that CASC intraslab events have a higher rate of attenuation than interface events consistent with their global GMM and Parker et al. (2022) assume that CASC interface events have a higher rate of attenuation than intraslab events consistent with their global GMM for interface events and their CASC-specific attenuation term for intraslab events. This uncertainty is not restricted to global and regional NGA-Sub models. Zhao et al. (2016aZhao et al. ( , 2016b found that intraslab events had a higher rate of attenuation than interface events, whereas Si et al. (2022) found both types of events to have the same rate of attenuation. These results indicate that our assumption of similar anelastic attenuation of interface and intraslab events is not unreasonable.
All of the NGA-Sub GMMs have stronger attenuation in the backarc region of CASC than that in the study by Abrahamson et al. (2016), but similar forearc attenuation at short periods as that by Abrahamson et al. (2016) and the M9 project simulations of Frankel et al. (2018). At long periods, the M9 simulations have less attenuation than all of the NGA-Sub models for reasons that are not known at this time. We conclude that there is a great deal of epistemic uncertainty associated with ground-motion attenuation in CASC that should be taken into account in a hazard analysis for this region.

Summary and conclusion
This article presents a partially nonergodic empirical GMM for subduction interface and intraslab events developed using the extensive global NGA-Sub database. The GMM is partially nonergodic in three aspects: (1) it includes an event term as a random effect; (2) it accounts for regional differences in overall amplitude, anelastic attenuation, linear site response, and basin response as random effects for subduction zones in AK, CAM, CASC, JP, NZ, SA, and TW; and (3) it accounts for regional differences in breakpoint magnitude between subregions of the seven regionalized subduction zones. In addition, there is a global model that can be used for sites outside of the defined regions or as an epistemic alternative to a regional model. Site terms are not evaluated in the GMM but can be estimated from published GMMs. Implementation guidelines are provided in the ''Engineering Application'' section of the article.
We recommend that our GMM be used together with other credible and appropriate GMMs for the site or region of interest to capture between-model epistemic uncertainty in a seismic hazard analysis. No single model, even if it includes within-model epistemic uncertainty as ours does, should be used to characterize seismic hazard. Our regionalized global model should be used along with its larger epistemic uncertainty when applying it to a subduction zone not characterized in this study. If some data are available for a new region, a Bayesian approach similar to that recommended by Stafford (2019) can be used to adjust the global GMM to this region. tasks of the NGA-Sub research program. Their contributions, dedication, and teamwork are greatly appreciated. We also thank three reviewers for their insightful comments that helped to improve the article.

Data and resources
The GMM developed in this study has been coded in multiple computer platforms including Visual Basic in Excel, R, MATLAB, and Python available at https://www.risksciences.ucla.edu/nhr3/gmtools. The coefficients of the model as well as an implementation of the GMM in the R software environment (R Core Team, 2023) is available at https://github.com/nikuehn/KBCG20 and the Supplemental material to the article and can be used to estimate median predictions and between-event, withinevent, and within-model standard deviations and their associated posterior distributions.