An offshore non-ergodic ground motion model for subduction earthquakes in Japan Trench area

Abstract

With the improvement of the world’s largest seafloor observation network (S-net) and the increase in the quantity and quality of records, the ergodic assumptions can be further relaxed in the modeling of offshore ground motion models (GMMs). This allows accounting for systematic and repeatable source, site, and path effects to further understand the characteristics of offshore ground motion in the Japan Trench region. We developed an offshore ergodic backbone GMM for subduction earthquakes and classified the sites into four categories using horizontal–vertical response spectral ratio to investigate site amplification. The offshore ergodic GMM is applicable for subduction earthquake scenarios with moment magnitudes ranging from 4.0 to 7.4 and rupture distances ranging from 10 to 300 km. Comparing offshore ergodic GMMs with onshore GMMs for subduction earthquakes, we found that offshore GMMs were significantly different from onshore GMMs, especially in the long-period and unburied states. Then a new offshore non-ergodic GMM was developed based on the offshore ergodic GMM. The systematic and repeatable source and site effects were captured by the spatially varying coefficients represented by Gaussian processes, while the systematic and repeatable path effects were captured by cell-specific anelastic attenuation proposed by Dawood and Rodriguez-Marek (2013), calculated with the Cohen-Sutherland computer graphics algorithm. The non-ergodic GMM revealed systematic and repeatable source, site, and path effects that were not captured by the ergodic GMM. Moreover, the non-ergodic GMM showed reduced aleatory variability and epistemic uncertainty on ground motion estimation compared to ergodic GMM. The reduction of aleatory variability and epistemic uncertainty had a significant impact on probabilistic seismic hazard analysis. Quantifications of these results are contributed to conduct reasonable seismic design and seismic risk assessment for marine engineering in offshore regions vulnerable to strong subduction earthquakes.

Introduction

The current research indicates that the offshore ground motion is significantly different from the onshore ground motion in response spectrum (Chen et al., 2015; Dhakal et al., 2021; Hu et al., 2020; Zhang and Zheng, 2019), temporal characteristics (Zhang and Zheng, 2020), elastic relative and absolute input energy (Hu and Tan, 2022), vertical-to-horizontal spectrum ratio (Boore and Smith, 1999; Chen et al., 2017; Tan and Hu, 2021), damping modification factors (Hu et al., 2022b), and seismic ground motion simulation (Lan et al., 2021; Li et al., 2011), especially for subduction earthquakes. Furthermore, Tan and Hu (2023) compared the differences in horizontal peak ground acceleration (PGA), spectral acceleration, frequency content, and significant duration between offshore and onshore ground motion characteristics using two large offshore earthquakes (February 13, 2021, M_w 7.1 and March 17, 2022, M_w 7.4), wherein the offshore records were from the S-net network and the onshore records were from the K-NET network. Cui et al. (2023) employed the same method to analyze offshore and onshore ground motion characteristics for a large earthquake. Both studies revealed significant differences and factors influencing the characteristics of offshore and onshore ground motion.

Strong offshore ground motion could cause the destruction of structures (Li et al., 2022), such as the collapse and damage of cross-sea bridges and submarine tunnels (Zhang et al., 2016), the leaking of marine oil drilling platforms (Borrero, 2005), the fracturing of submarine communication cables (Bi and Hao, 2016), and the large-scale liquefaction of artificial island backfills (Tanaka, 2000). For example, several earthquakes, such as the 2011 M_w 9.0 Tohoku-Oki subduction interface earthquake (Simons et al., 2011) and the 2021 M_w 7.1 Fukushima subduction intraslab earthquake (Kiser and Kehoe, 2021; Kubota et al., 2021), have revealed substantial economic losses and secondary disasters (e.g. tsunami, mudslide, etc.) caused by ground motions in Japan trench area. Therefore, with the development of performance-based earthquake engineering, it is necessary to conduct reasonable seismic design and seismic risk assessment of marine engineering in offshore regions that are easily affected by strong subduction earthquakes, such as the east and south coasts of China, Chile, the west coast of the United States, New Zealand, and Japan.

Under the current concept of performance-based seismic engineering, probabilistic seismic hazard analysis (PSHA) is based on empirical ground motion models (GMMs). In PSHA, the slope of the hazard curve is sensitive to aleatory standard deviation (SD) of empirical GMMs (Bommer and Abrahamson, 2007). The ergodic assumption is convenient because it combines ground motion data from similar tectonic environments around the world to develop ergodic GMMs based on magnitude, distance, and site conditions. For partially non-ergodic GMMs, there are two descriptions: one can only capture the systematic site effects, known as single station GMMs; the other is that GMMs have different coefficient values for different regions, such as constant term, anelastic attenuation term, site term, and basin amplification term. Fully non-ergodic GMMs include non-ergodic terms to capture the systematic effects (e.g. source, path, and site). For example, the NGA-west1 (e.g. Abrahamson and Silva, 2008; Boore and Atkinson, 2008; Campbell and Bozorgnia, 2008; Chiou and Youngs, 2008; Idriss, 2008) GMMs are fully ergodic GMMs for shallow crustal tectonic settings. In addition, four of the NGA-west2 (e.g. Abrahamson et. al., 2014; Boore et al., 2014; Campbell and Bozorgnia, 2014; Chiou and Youngs, 2014) GMMs are partially non-ergodic GMMs for shallow crustal tectonic settings because it contains the regional terms in their models. However, there are still regional differences in ergodic GMMs over a wide range of regions, which has led to the development of GMMs for specific regions and non-ergodic GMMs. For instance, many regional ergodic GMMs have already been developed for the Japanese subduction zone (e.g. Si et al., 2020; Zhao et al., 2016a, 2016b; Hassani and Atkinson, 2021). Furthermore, the NGA-Subduction (NGA-Sub) project (Bozorgnia et al., 2022; Mazzoni et al., 2022) established a large subduction ground motion database containing a large number of strong earthquake records from subduction zones around the world and developed new GMMs for subduction earthquakes (e.g. Abrahamson and Gülerce, 2020, 2022; Gregor et al., 2022; Kuehn et al., 2020; Parker et al., 2020, 2022). Note that the NGA-Sub GMMs are also partially non-ergodic GMMs that capture the average differences in ground motion scaling among the seven broad regions by accounting for site term, the linear distance term, basin term, constant term, and break point for large-magnitude scaling. In addition, several recent studies (e.g. Abrahamson et al., 2019; Kuehn, 2021; Kuehn et al., 2019; Landwehr et al., 2016; Lavrentiadis et al., 2021; Lavrentiadis and Abrahamson, 2021; Liu et al., 2022a, 2022b; Macedo and Liu, 2022; Sung et al., 2021) have shown regional differences in the scaling of ground motion represented by the non-ergodic GMMs, revealing systematic and repeatable source, path, and site effects. These studies indicate that the non-ergodic model can capture epistemic uncertainty of spatially varying effects. If the trade-off between the aleatory variability and epistemic uncertainty for seismic risk is insufficient, ergodic GMMs may lead to biased estimates of intensity measures (IMs) hazard curves when performing PSHA (Abrahamson et al., 2019).

With the improvement of the offshore strong-motion network (S-net) and the quantity and quality of subduction earthquake records in the Japan Trench region, the ergodic assumption can be further relaxed in the modeling of offshore GMMs to account for systematic and repeatable source, site, and path effects. In this study, we introduced the S-net seafloor ground motion database and utilized a horizontal-to-vertical response spectral ratio (H/V) with a damping ratio of 5% for each station to classify the site. We then developed the offshore ergodic backbone GMM for subduction earthquakes in the Japan Trench area. Furthermore, we compared our offshore GMM with onshore GMM for subduction earthquakes to gain further insight into offshore ground motion characteristics. We further developed the offshore non-ergodic GMM based on the offshore ergodic backbone GMM. The systematic and repeatable source, site, and path effects of the offshore non-ergodic GMM were modeled using Gaussian processes and cell-specific attenuation approach, respectively. Finally, we discussed the residuals and aleatory SDs of the non-ergodic GMM and the systematic and repeatable source, site, and path effects captured by the non-ergodic GMM, and drew some conclusions. The results of this article are expected to provide a basis for revealing the source, site, and path effects and offshore ground motion characteristics for subduction earthquakes in the Japan Trench region.

S-net seafloor ground motion database

Introduction of S-net seafloor observation network

The S-net is situated in the northeast region of Japan, aligning with the arc boundary where the Pacific plate subducted beneath the Eurasian plate. Its purpose is to enhance the Japan Meteorological Agency’s earthquake early warning (EEW) and tsunami early warning systems (Aoi et al., 2020).

Initiated in 2011 following the catastrophic 9.1 Tohoku-Oki earthquake, the construction of S-net was completed in March 2017. The seafloor observation network consists of six subsystems: S1, S2, S3, S4, S5, and S6. Figure 1 illustrates their geographical distribution. It comprises 150 seafloor observation stations, spanning an area of approximately 1000 km by 300 km from off-Hokkaido to off-Kanto. Stations situated in bathymetry below 1500 m (as shown in Figure 1), numbering 41, are buried in grooves approximately 1 m beneath the seafloor to circumvent fishery activities, minimize ground noise, and enhance ground coupling. The remaining 109 stations are placed freely on the seafloor in areas with bathymetry deeper than 1500 m. These stations are spaced at approximately 30 (50–60) km intervals in the east-west (north-south) direction, perpendicular to (parallel with) the Japan Trench axis. Each subsystem contains 22–28 stations connected in series to land stations via ocean-bottom fiber optic cables, with lengths ranging from 800 to 1500 km. The total length of these cables is approximately 5500 km. The S-net observation device contains built-in sensors, including accelerometers, velocimeters, and pressure gauges. The device has a cylindrical appearance and is placed horizontally along the direction of the cable.

Figure 1.

The geographical distribution of the S-net seafloor observation network and earthquake epicenters. The color bar indicates the bathymetry at each observation station. The blue and red hypocentral mechanism symbols represent interface events and intraslab events, respectively. The hypocentral depth of the interface event is smaller than 50 km and that of the intraslab event is greater than 50 km.

When an earthquake occurs, S-net transmits continuous waveform data to the National Research Institute for Earth Science and Disaster Resilience (NIED) data center via the ocean-bottom fiber optic cables (Dhakal et al., 2021).

Data selection and processing

First, we gathered earthquake events in the Japan Trench area from the Full Range Seismograph Network of Japan (F-net), with the geographic distribution of earthquake epicenters displayed in Figure 1. Based on the earthquake catalog, we downloaded 10 min continuous waveform data recorded by accelerometers from the S-net seafloor observation network, starting 1 min before the earthquake occurred. These data include X components parallel to the cable’s long axis and Y and Z components perpendicular to it.

S-net sensors are prone to rotations during strong ground motion, complicating the recovery of accurate strong ground motion data (Takagi et al., 2019). We estimated tilt and rotation angles using the matrix operation proposed by Takagi et al. (2019) and converted the X, Y, and Z components of continuous waveform data into east, north, and up (ENU) components. Figure 2a displays rotation angles of original acceleration records greater than 0.1°, while Figure 2b shows acceleration records before and after matrix operation.

Figure 2.

The M_w7.1 off-Fukushima earthquake in 2021: (a) the original waveform of acceleration at station S2N01; (b) The comparison between the original waveform and the converted waveform at station S3N17.

We employed the short-time average over long-time average (STA/LTA) algorithm proposed by Trnkoczy (2002) to extract P-waves and determined noise windows. The filter corner was determined by comparing the Fourier amplitude spectra (FAS) of ground motion records with those of noise windows. Figure 3 illustrates the plot for selecting filter corner frequencies. For the low-pass/high-pass corner frequency (fc_LP/fc_HP), we used a minimum/maximum signal-to-noise ratio (SNR) equal to 0.5 (Bahrampouri et al., 2021; Boore and Bommer, 2005; Douglas and Boore, 2011).

Figure 3.

Illustrative plot for selecting the filter corner frequencies.

To achieve the research objectives of this article and account for the model’s influencing factors, the following criteria were used to exclude certain ground motion records in the model’s application:

1. Records with rotations greater than 0.1° (as shown in Figure 2a).

2. Records with an SNR less than three for the frequency band of 1∼5 Hz (as shown in Figure 3).

3. Records with obvious low-quality, such as those with multiple wave packets or missing pre-event memory.

4. Records with rupture distances greater than 300 km.

5. Records that do not capture the start of the P-wave.

6. Records with PGA lower than 1e−3 g.

7. Records that do not belong to subduction earthquake events, as determined by the method proposed by Zhao et al. (2015). This means that earthquake events are not located near the subduction zone (as shown in Figure 1).

8. An earthquake with fewer than three records.

The final data set used in this study comprised 5061 sets of three-component records from 194 intraslab earthquakes and 5270 sets of three-component records from 214 interface earthquakes between 2016 and 2023. Figure 4 illustrates the distribution relationship between moment magnitude (M) and rupture distances (R_rup), with the density map of moment magnitude and rupture distance representing the centralized distribution of ground motion data. Rupture distance was calculated using the method proposed by Kaklamanos et al. (2011).

Figure 4.

The scatter distribution of ground motion data between moment magnitude (M_w) and rupture distances (R_rup).

The NGA-west2 data processing flow was used for record processing (Ancheta et al., 2014), which included steps such as removing pre-event mean from acceleration, tapering, and padding. The ground motion records used in this article were uniformly processed by baseline correction and fourth-order acausal Butterworth band-pass filter. These records were processed using empirical values for high-pass corner frequencies proposed by Puglia et al. (2018). In contrast, a fixed frequency of 30 Hz was used for the low-pass corner frequency. Douglas and Boore (2011) noted that low-frequency noise (<1 Hz) significantly impacted strong motion intensity parameters, and filtering high-frequency noise was only necessary in specific situations.

Furthermore, the two horizontal components of ground motion used in this study are RotD50 (Boore, 2010), which represents the median ground motion at all rotation angles and is orientation-independent.

Site classification of offshore stations

The estimation of site effect is commonly done using the time-averaged shear-wave velocity over the top 30 m (V_s30, m/s) (e.g. Abrahamson and Gülerce, 2020, 2022; Gregor et al., 2022; Kuehn et al., 2020; Parker et al., 2020, 2022). However, an alternative parameter for site effect estimation is the site predominant period (T_n), which has been increasingly used in recent years (e.g. Kwak and Seyhan, 2020; Zhao and Xu, 2013). V_S30 reflects the stiffness of the shallow soil layer and T_n depends on the stiffness of the soil layer and the depth of sediment to the bedrock. Both measurements are associated with site amplification (Arteta et al., 2021; Zhao and Xu, 2013; Zhu et al., 2020).

Due to the lack of V_s30 data for S-net stations, we utilized the horizontal-to-vertical response spectral ratio (HVRSR) with a damping ratio of 5% for each station to classify the sites. To ensure accuracy, we excluded records with PGA greater than 0.05 g, as high-amplitude records may cause nonlinear behavior in the underlying soil (Arteta et al., 2021; Kawase et al., 2019).

The sites were then categorized based on the predominant period (T_n) of each station, resulting in four categories (SC1, SC2, SC3, and SC4). The classification criteria were modified based on the criteria proposed by Idini et al. (2017) and Zhao (2006). The T_n range of site classification was determined using the method proposed by Idini et al. (2017), and the corresponding approximate range of V_s30 was based on the work proposed by Zhao (2006). However, the method of predominant period does not consider the factor of spectral shape. If there are not many records at a station or there are many peaks in the HVRSR, it is difficult to judge according to the predominant period (T_n). To improve the classification accuracy, Zhao (2006) proposed a classification method based on spectral shape index, and Ghasemi et al. (2009) proposed a classification method based on Spearman rank correlation coefficient. The site classification index (SI) was calculated using the method proposed by Ghasemi et al. (2009) to improve the classification accuracy, as follows:

$$SI=1-6\sum \frac{d_i^2}{n(n^2-1)}$$

(1)

in which n is the total number of periods, d_i is the difference between the mean HVRSR of the target category and the mean HVRSR for a station. Figure 5 shows the HVRSR curves of some selected stations from S-net. The gray line and black line represent the HVRSR curves of each record and the mean HVRSR curve for a station, respectively.

Figure 5.

The HVRSR curves of some selected stations: (a) Individual and mean HVRSR curves at station S1N07; (b) Individual and mean HVRSR curves at station S2N09; (c) Individual and mean HVRSR curves at station S3N06; (d) Individual and mean HVRSR curves at station S3N11.

Site categories SC1, SC2, SC3, and SC4 were classified depending on T_n and SI according to Table 1. Table 1 provides the site categories, predominant period ranges, corresponding approximate ranges of V_s30, and the number of stations and records. Figure 6 shows the comparison results of the number of four site categories for buried and unburied stations. The figure shows a more complicated site condition for buried stations than unburied stations (ie., There are four site categories for buried stations and three for unburied stations). The number of unburied stations is mainly distributed in site categories SC1and SC2, whereas the number of buried stations is mainly distributed in site categories SC3 and SC4. This means that most of the unburied stations are located at stiff and hard site conditions, while most of the buried stations are located at medium and soft site conditions.

Table 1.

Classification criteria of every category and the number of stations and records

Category	Predominant period, Tn(s)	Site description	Corresponding approximate range of Vs30 (m/s)	No. of buried stations	No. of unburied stations	No. of records (Interface)	No. of records (Intraslab)
SC1	Tn ≤ 0.2	Rock/stiff soil	Vs30 ≥ 600	7	40	1640	1649
SC2	0.2 < Tn ≤ 0.4	Hard soil	300 ≤ Vs30 < 600	7	59	2261	2191
SC3	0.4 < Tn ≤ 0.8	Medium soil	200 ≤ Vs30 < 300	16	10	1042	940
SC4	Tn > 0.8	Soft soil	Vs30 < 200	11	0	327	281
Total				41	109	5270	5061

Figure 6.

The comparison results of the number of four site categories for buried and unburied stations.

Figure 7 shows the mean HVRSR and the standard deviation (SD) of all stations classified with the predominant period for the site categories SC1, SC2, SC3, and SC4. Figure 7a clearly illustrates that the amplitudes and shapes of the mean HVRSR of the four site categories are significantly different. The mean HVRSR amplitude of unburied stations is greater than that of buried stations. Meanwhile, Figure 7b shows that each of the four site categories has a large SD near the peak amplitude of the mean HVRSR.

Figure 7.

(a) Mean H/V response spectral ratio and (b) corresponding standard deviations for all data from the S-net stations under buried and unburied states for the site categories SC1, SC2, SC3, and SC4.

Figure 8 presents the geographical distribution for the site categories SC1, SC2, SC3, and SC4. Site categories SC1 and SC2 are predominantly located near the Japan Trench, with most of stations being unburied states. On the contrary, site categories SC3 and SC4 are mostly distributed near the onshore, with most of stations being buried states.

Figure 8.

The geographical distribution for the site categories SC1, SC2, SC3, and SC4.

Offshore ergodic backbone GMM for subduction earthquakes

The functional form of the ergodic backbone GMM we adopted here is modified based on the functional form of the GMM of Si et al. (2022) for subduction earthquakes. The modified functional form is as follows:

$${\log }_{10}A(T)=b(T)+g(T,R_{\mathrm{rup}},M,D)-k(T)R_{\mathrm{rup}}+G(T)+\delta B_e+\delta W_{\mathrm{es}}$$

(2)

in which A(T) represents the observed ground motion IM, such as PGA and the RotD50 values of 5%-damped spectral acceleration for different periods; M is moment magnitude; R_rup is rupture distance; D is hypocentral depth. δB_e is the inter-event residual for eth event and δW_es is the intra-event residual for the sth record from the eth event. The residuals δB and δW are uncorrelated and conform to the normal distribution with variances τ² and ϕ², respectively.

The b(T) in Equation 2 is the source term, and the functional form adopted is shown in Equation 3. The magnitude scaling contains a linear term and a nonlinear term. Coefficient a₀(T) captures the linear scaling associated with magnitude, and coefficient a₁(T) captures the nonlinear scaling associated with magnitude. Coefficients d₀(T) and d₁(T) are the calibrations of earthquake types and are used to reflect differences between interface and intraslab events in GMM. Coefficient e(T) is the constant term.

$$b(T)=a_0(T)M+a_1(T)M^2+d_0(T)S_0+d_1(T)S_1+h(T)D+e(T)$$

(3)

where

$$S_0=\{\begin{array}{c}1\begin{array}{cc}\mathrm{Interface}\hfill & \mathrm{earthquake}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{0\begin{array}{cc}\mathrm{Intraslab}\hfill & \mathrm{earthquake}\hfill \end{array}}\end{array}$$

(4)

$$S_1=\{\begin{array}{c}0\begin{array}{cc}\mathrm{Interface}\hfill & \mathrm{earthquake}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{1\begin{array}{cc}\mathrm{Intraslab}\hfill & \mathrm{earthquake}\hfill \end{array}}\end{array}$$

(5)

The path scaling includes the effect of geometric spreading (log R term) and the anelastic attenuation term (linear R term). The g(T,R_rup, M,D) in Equation 2 is the geometric spreading term, and the k(T)R_rup is the anelastic attenuation term. C(T, M) represents the near-source saturation term associated with period and magnitude changes. The anelastic attenuation coefficient k(T) is modified to capture the offshore ground motion characteristic associated with period. These two terms are correlated in most GMM regressions (Arteta et al., 2021). The geometrical spreading coefficient is used to avoid positive value of the anelastic attenuation coefficient. The data of Moho depth obtained from CRUST 1.0. ¹

$$\begin{array}{c}\multicolumn{1}{c}{g(T,R_{\mathrm{rup}},M,D)=\{\begin{array}{c}-{\log }_{10}[R_{\mathrm{rup}}+C(T,M)]\begin{array}{c}D\le \mathrm{Moho}\begin{array}{cc}\mathrm{depth}\hfill & \begin{array}{c}\mathrm{or}\hfill \end{array}\hfill \end{array}{R}_{\mathrm{rup}}<1.7D\hfill \end{array}\hfill \\ \multicolumn{1}{c}{0.6{\log }_{10}[1.7D+C(T,M)]-1.6{\log }_{10}[R_{\mathrm{rup}}+C(T,M)]\begin{array}{c}D>\mathrm{Moho}\begin{array}{cc}\mathrm{depth}\hfill & \begin{array}{c}\mathrm{and}\hfill \end{array}\hfill \end{array}{R}_{\mathrm{rup}}\ge 1.7D\hfill \end{array}}\end{array}}\end{array}$$

(6)

where

$$C(T,M)=\{\begin{array}{c}0.0055\times {10}^{0.5M}\begin{array}{c}T<0.3\mathrm{s}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{[0.00081-0.00897{\log }_{10}\left(T\right)\times {10}^{0.5M}\begin{array}{c}0.3\mathrm{s}\le T\le 0.6\mathrm{s}\hfill \end{array}}\\ \multicolumn{1}{c}{0.0028\times {10}^{0.5M}\begin{array}{c}T>0.6\mathrm{s}\hfill \end{array}}\end{array}$$

(7)

$$k(T)=\{\begin{array}{c}0.00135\begin{array}{c}T<0.3s\hfill \end{array}\hfill \\ \multicolumn{1}{c}{-0.00055-0.00153{\log }_{10}(T)\begin{array}{c}0.3s\le T\le 1s\hfill \end{array}}\\ \multicolumn{1}{c}{0.00055\begin{array}{c}T>1s\hfill \end{array}}\end{array}$$

(8)

For the estimation of site effects, we introduced a dummy variable (W) with values of 0 and 1 for buried and unburied stations, respectively, as well as each site category in section “Site classification of offshore stations.”Tan and Hu (2023) analyzed the characteristics of offshore ground motion using two largest offshore earthquakes recorded by the S-net. They examined the influence of water column and sediment thickness using residual analysis. Their results showed that the residuals of unburied stations were greater than those of buried stations. Hu et al. (2022a) also found that whether stations are buried or not has a significant influence on the ground motion IM. After residual analysis (see section “Residual analysis of ergodic backbone GMM”), the residuals of spectral acceleration can be evenly distributed around zero with the thickness of sediment layer (D14) and water column, respectively. Therefore, the difference between spectral acceleration and the sediment thickness can be well captured by dividing S-net stations into buried and unburied stations. This means that whether stations are buried or not can be considered as a predictor variable in the ergodic GMM to reflect the influence of water column and the sediment thickness. In addition, several studies showed that site classification can capture the resonant characteristics of a group of sites, allowing adequate description of site amplification on the establishment of GMM (Arteta et al., 2021; Idini et al., 2017; Zhao and Xu, 2013). Based on the above findings, we propose the following functional forms to estimate site effect:

$$G(T)=w(T)W+s_i(T)S{C}_i$$

(9)

where s_i is a regression coefficient and i refers to the ith site category.

The regression coefficients (a₀, a₁, d₀, d₁, h, e, w and s_i) were determined by mixed-effect regression (Abrahamson and Youngs, 1992; Bates et al., 2015). The regression coefficients and SDs of residuals are presented in Table 2. σ is the total SD of the total residual (δB_e + δW_es) in Table 2, which can be calculated as $\sqrt{{\tau }^2+{\varphi }^2}$ .

Table 2.

Regression coefficients and standard deviations of the offshore ergodic backbone GMM for subduction earthquakes

Period(s)	a0	a1	d0	d1	h	e	s1	s2	s3	s4	w	τ	ϕ	σ
0	−0.7103	0.1177	0.1651	0.2122	0.0023	0.1882	0.1540	0.2301	0.1595	0.1628	0.2937	0.146	0.249	0.289
0.01	−0.7103	0.1177	0.1651	0.2122	0.0023	0.1882	0.1540	0.2301	0.1595	0.1628	0.2937	0.146	0.249	0.289
0.02	−0.7103	0.1177	0.1651	0.2122	0.0023	0.1882	0.1540	0.2301	0.1595	0.1628	0.2937	0.146	0.249	0.289
0.03	−0.7135	0.1179	0.1701	0.2176	0.0024	0.1934	0.1606	0.2353	0.1631	0.1668	0.2915	0.146	0.249	0.289
0.05	−0.7153	0.1182	0.1786	0.2271	0.0025	0.2023	0.1748	0.2456	0.1647	0.1678	0.2771	0.147	0.25	0.29
0.08	−0.6056	0.1076	0.1091	0.1652	0.0030	0.1366	0.1317	0.1753	0.0793	0.0904	0.2321	0.153	0.256	0.298
0.1	−0.5097	0.0972	0.0552	0.1168	0.0033	0.0854	0.0998	0.1155	0.0137	0.0471	0.2385	0.157	0.261	0.304
0.15	−0.6557	0.1075	0.2701	0.3306	0.0035	0.2998	0.3405	0.3266	0.1990	0.2538	0.2617	0.158	0.274	0.317
0.2	−0.7708	0.1188	0.3899	0.4358	0.0029	0.4124	0.3909	0.4754	0.3155	0.3908	0.3208	0.152	0.275	0.314
0.25	−0.6475	0.1102	0.2611	0.2987	0.0024	0.2795	0.2154	0.3502	0.2337	0.2523	0.3329	0.147	0.272	0.309
0.3	−0.4804	0.0985	0.0762	0.1122	0.0019	0.0939	0.0246	0.1579	0.0793	0.0479	0.3211	0.145	0.266	0.303
0.4	−0.2136	0.0791	−0.2278	−0.1987	0.0008	−0.2136	−0.3062	−0.1588	−0.1899	−0.2202	0.3138	0.143	0.265	0.301
0.5	−0.0135	0.0649	−0.4678	−0.4345	−0.0001	−0.4514	−0.5631	−0.4029	−0.3883	−0.4501	0.2954	0.14	0.264	0.298
0.6	0.1845	0.0506	−0.7105	−0.6744	−0.0007	−0.6928	−0.8121	−0.6480	−0.6143	−0.6807	0.3066	0.141	0.264	0.299
0.75	0.4307	0.0327	−1.0163	−0.9685	−0.0018	−0.9929	−1.1076	−0.9606	−0.9005	−0.9666	0.3199	0.149	0.265	0.304
1	0.6466	0.0183	−1.3192	−1.2641	−0.0033	−1.2922	−1.4090	−1.2659	−1.1951	−1.2401	0.3383	0.147	0.262	0.3
1.5	0.9036	0.0028	−1.7173	−1.6548	−0.0041	−1.6867	−1.7975	−1.6668	−1.5895	−1.6234	0.3441	0.142	0.26	0.296
2	1.0030	−0.0021	−1.9217	−1.8579	−0.0043	−1.8904	−1.9905	−1.8725	−1.8029	−1.8331	0.3537	0.135	0.25	0.284
2.5	0.9561	0.0050	−1.9643	−1.9147	−0.0041	−1.9400	−2.0241	−1.9322	−1.8571	−1.8774	0.3495	0.129	0.236	0.269
3	0.7982	0.0216	−1.8886	−1.8421	−0.0041	−1.8658	−1.9442	−1.8617	−1.7823	−1.8058	0.3242	0.124	0.225	0.257
4	0.4750	0.0522	−1.6834	−1.6383	−0.0040	−1.6613	−1.7268	−1.6651	−1.5815	−1.5980	0.2843	0.124	0.212	0.245
5	0.2138	0.0743	−1.5013	−1.4516	−0.0040	−1.4769	−1.5450	−1.4753	−1.4005	−1.4232	0.2780	0.129	0.21	0.247
6	−0.0216	0.0944	−1.3339	−1.2726	−0.0041	−1.3039	−1.3767	−1.2939	−1.2327	−1.2632	0.2701	0.133	0.207	0.246
8	−0.3428	0.1210	−1.1103	−1.0633	−0.0032	−1.0872	−1.1680	−1.0739	−1.0183	−1.0320	0.2554	0.136	0.212	0.252
10	−0.5717	0.1376	−0.9339	−0.8926	−0.0026	−0.9137	−1.0000	−0.8996	−0.8457	−0.8409	0.2639	0.14	0.218	0.259

Residual analysis of ergodic backbone GMM

Figure 9 shows residual analysis with respect to moment magnitude (M_w), rupture distance (R_rup), the thickness of sediment up to Vs = 1.4 km/s (D14), and water column for PGA and spectral accelerations at T = 0.5 and 1.0 s. Figure 9a1 to a3 shows the plots of inter-event residuals versus moment magnitude. We note that most of the data correspond to earthquakes with magnitudes less than about 6.5. A possible bias at lager magnitudes is very close to 0; however, the low number of data points within that magnitude range implies that ergodic GMM can capture them. Therefore, these residuals can be considered as epistemic uncertainty for the ergodic GMM due to the lack of data. As may be observed from Figure 9b1 to b3 and Figure 9c1 to c3, the ergodic GMM, in general, is unbiased for all considered magnitude and distance ranges. Ergodic GMM considers whether stations are buried or not as a site predictor variable, and the distribution of the intra-event residuals with respect to the sediment thickness (e.g. Figure 9d1 to d3) and the water column (e.g. Figure 9e1 to e3) are uniformly distributed around zero. This indicates that there is no need to consider water column and sediment thickness as site predictor variables in the ergodic GMM. Figure 10 shows the intra-event residuals with respect to tectonic regime, buried state, and each site category for PGA and spectral accelerations at T = 0.5 and 1.0 s. The residuals and the mean of residuals are centered at zero. This indicates that the ergodic GMM can well reflect the differences of tectonic regime, buried state, and each site category.

Figure 9.

Residual analysis with respect to moment magnitude (M_w), rupture distance (R_rup), the thickness of sediment up to Vs = 1.4 km/s (D14), and water column for PGA and spectral accelerations at T = 0.5 and 1.0 s. (a1)–(a3) the inter-event residuals versus moment magnitude (M_w); (b1)–(b3) the intra-event residuals versus moment magnitude (M_w); (c1)–(c3) the intra-event residuals versus rupture distance (R); (d1)–(d3) the intra-event residuals versus thickness of sediment up to Vs = 1.4 km/s (D14); and (e1)–(e3) the intra-event residuals versus water column. Orange circles and green error bars represent the mean and standard deviation of the residuals, respectively.

Figure 10.

The intra-event residuals with respect to tectonic regime ((a1)–(a3)), buried state ((b1)–(b3)), and each site category ((c1)–(c3)) for PGA and spectral accelerations at T = 0.5 and 1.0 s.

Comparison of offshore ergodic GMM for subduction earthquakes

Figure 11 presents the comparison of median response spectra at different buried states for different site categories with M = 7 and R_rup = 100 km for interface and intraslab earthquakes. For buried stations, the four site categories exhibit varying site amplifications, with peaks occurring at different periods. The generic soil sites (category SC2) exhibit maximum amplification values at T = 0.25 s. Categories SC3 and SC4 spectra exhibit less amplification and a major period shift at the peak, in comparison to category SC1, which corresponds to stiff soil or rock sites. The figure shows that unburied stations show a similar site amplification to the buried stations and a large difference and amplitude for unburied stations.

Figure 11.

The comparison of median response spectra at different buried states for different site categories with M = 7 and R = 100 km for (a) interface and (b) intraslab earthquakes. D = 30 km for interface earthquakes and D = 50 km for intraslab earthquakes, Moho = 30 km.

Figure 12 shows the comparison of the median response spectra predicted by the offshore ergodic backbone GMM for different hypocentral depths with M = 7 and buried states at a distance of 150 km for interface and intraslab earthquakes, considering a reference site category SC1. The comparison of the median response spectra values for different hypocentral depths shows a large amplitude for deep earthquakes and a stronger depth dependence at shorter periods than at longer periods.

Figure 12.

Comparison of the median response spectra predicted by the offshore ergodic backbone GMM for different hypocentral depths with M = 7 and buried states at a distance of 150 km for (a) interface and (b) intraslab earthquakes, considering a reference site category SC1. D = 30 km for interface earthquakes and D = 50 km for intraslab earthquakes; Moho = 30 km.

Figure 13 shows the comparison results of PGA and Sa_RotD50 (T = 1 s) attenuation curves at different M values and buried states predicted by the offshore ergodic backbone GMM for interface and intraslab earthquakes with a reference site category SC1. This figure shows a clear amplitude-saturation for lager magnitude events at short distance and a strong M dependence at different periods. The amplitude of intraslab events at the same magnitude is greater than that of interface events and unburied stations show a larger amplitude than buried stations. The comparison of PGA and Sa_RotD50 (T = 1 s) attenuation curves for smaller earthquake events shows a steeper attenuation with distance than for larger earthquake events.

Figure 13.

The comparison results of ((a) and (b)) PGA and ((c) and (d)) Sa_RotD50 (T = 1s) attenuation curves at different M values and buried states predicted by the offshore ergodic backbone GMM for interface (left panel) and intraslab (right panel) earthquake events with a reference site category SC1. D = 30 km for interface earthquakes and D = 50 km for intraslab earthquakes; Moho = 30 km.

Figure 14 shows the comparison results of median response spectra at different M values and buried states predicted by the offshore ergodic backbone GMM at R_rup = 60 km and R_rup = 120 km for interface and intraslab earthquakes with a reference site category SC1. The figure shows the intraslab earthquakes have lager amplitude than interface ones at the same rupture distance, and the long-period attenuation become slower as M increases. The amplitude of median response spectra decreases as R_rup increases.

Figure 14.

The comparison results of median response spectra at different M values and buried states predicted by the offshore ergodic backbone GMM at ((a) and (b)) R_rup = 60 km and ((c) and (d)) R_rup = 120 km for (left panel) interface and (right panel) intraslab earthquakes with a reference site category SC1. D = 30 km for interface earthquakes and D = 50 km for intraslab earthquakes; Moho = 30 km.

In recent years, two of the largest offshore subduction-slab earthquakes (February 13, 2021, M 7.1 and March 17, 2022, M 7.4) have caused significant damage to onshore structures. Figure 15a and b compare the offshore ground motion of the 2022 M7.4 off-Fukushima earthquake recorded by S-net with the median values predicted by the offshore ergodic backbone GMM for PGA and Sa_RotD50 (T = 1 s). The PGA and Sa_RotD50 (T = 1 s) of offshore ground motion are significantly different between the buried and unburied states. The offshore ergodic backbone GMM captures the characteristics of the offshore ground motion well in terms of distance scaling, both in buried and unburied states. Figure 15c and d compare the offshore ground motion of the 2021 M 7.1 off-Fukushima earthquake with the median values predicted by the offshore ergodic backbone GMM for PGA and Sa_RotD50 (T = 1 s). Note that the offshore ergodic backbone GMM is developed based on data less than 300 km, but the ergodic GMM in Figure 15 demonstrates that it is still somewhat predictive for data greater than 300 km. In addition, the PGA attenuation curves of buried and unburied states exhibit consistent attenuation with observation data greater than 300 km. However, there are some differences between the Sa_RotD50 (T = 1 s) attenuation curves of buried and unburied states and observed ground motion data greater than 300 km. For Sa_RotD50 (T = 1 s) observation data, ground motion data of unburied stations exhibits faster attenuation beyond 300 km, while ground motion data of buried stations exhibits slower attenuation. This may be due to the increased ground coupling caused by buried stations, resulting in slow attenuation at long-period and long-distance.

Figure 15.

The comparison of offshore ground motion recorded by S-net with median values predicted by the offshore ergodic backbone GMM at PGA (left panel) and Sa_RotD50 (T = 1 s) (right panel) for ((a) and (b)) 2022 M7.4 off-Fukushima earthquake and ((c) and (d)) 2021 M 7.1 off-Fukushima earthquake with a reference site category SC1.

Comparison of the proposed GMM with previous studies for subduction earthquakes

We compared the ergodic GMM with NGA-Sub GMMs (e.g. Abrahamson and Gülerce, 2020, referred to as AG20; Parker et al., 2020, referred to as PSBAH20; Si et al., 2020, referred to as SMK20) and some subduction GMMs from Japan (e.g. Zhao et al., 2016a, referred to as Zea16; Zhao et al., 2016b, referred to as Zea16; Hu et al., 2020, referred to as HTZ20). Note that the HTZ20 is a single-station offshore GMM based on onshore and six offshore stations, and we choose the mean values of the six offshore stations for comparison. For the NGA-Sub GMMs, the magnitude scaling breakpoint values for the Japanese Pacific tectonic plates are plotted. The median values predicted by the offshore ergodic backbone GMM are presented at buried and unburied states for interface and intraslab earthquakes with a reference site category SC2. The comparison results of PGA and spectral acceleration at period of 1 s (Sa(T = 1 s)) attenuation curves are presented for interface and intraslab events with magnitude 7 for the defined site condition with V_s30 values of 500 m/s, which corresponds to the V_s30 range for the site category SC2. For the response spectra, the comparison results are provided at two distances of 75 and 200 km for interface and intraslab events with magnitude 7. Note that the Zea16 and HTZ20 are developed using the geomean horizontal component, whereas the NGA-Sub GMMs and our ergodic GMM are developed using the RotD50 horizontal component. For the distance scaling, the Zea16 and HTZ20 use the hypocentral distance, while our ergodic GMM and the NGA-Sub GMMs use the rupture distance. The comparison is a general comparison between our offshore ergodic GMM and other GMMs for subduction earthquakes, without using spectral period-dependent scaling (Boore and Kishida, 2017) to adjust the geometric horizontal component to the RotD50 component (Gregor et al., 2022), and without adjusting on the distance scaling.

Figure 16 shows the comparison results of PGA and Sa (T = 1 s) attenuation curves for interface and intraslab events with magnitude 7. Figure 16a and b shows the comparison results of PGA attenuation curves for interface and intraslab events. The offshore ergodic GMM and HTZ20 show a smaller difference in amplitude because the HTZ20 is developed based on data from six offshore stations and is able to represent the offshore ground motion of PGA. For the interface events, the amplitude of the PGA attenuation curve of the offshore GMM is larger than that of the NGA-Sub GMMs. In contrast, for intraslab events, the differences between the PGA attenuation curves of offshore GMM for the buried states and the PGA attenuation curves of onshore GMMs are relatively small. However, there are still significant differences between the PGA attenuation curves of offshore GMM for the unburied states and the PGA attenuation curves of onshore GMMs. Figure 16c and d shows the comparison results of Sa (T = 1 s) attenuation curves for interface and intraslab events. The Sa (T = 1 s) attenuation curves of the offshore GMM for buried and unburied states show a lager amplitude and a steeper attenuation than those of the NGA-sub GMMs.

Figure 16.

The comparison results of ((a) and (b)) PGA and ((c) and (d)) Sa (T = 1 s) attenuation curves for interface (left panel) and intraslab (right panel) events with magnitude 7 and a reference site category SC2. For NGA-Sub GMMs, Ztor = 10 km for interface events and Ztor = 50 km for intraslab events; hypocentral depth (Z_HYP) is 30 km for interface events and hypocentral depth is 50 km for intraslab events; V_s30 = 500 m/s.

Figure 17a and b shows the comparison results of response spectra for interface and intraslab events with magnitude 7 at a distance of 75 km. The offshore ergodic GMM with unburied states exhibits similar trends to the HTZ20 at short periods, but still has some differences at long periods. For interface events, the difference between offshore GMMs and onshore GMMs becomes larger as period increases. For intraslab events, the offshore ergodic GMM with buried states shows a similar trend to NGA-sub GMMs at short periods. However, the difference between offshore GMMs and onshore GMMs increases as period increases, which is similar to the trend of interface events. Figure 17c and d shows the comparison results of response spectra for interface and intraslab events with magnitude 7 at a distance of 200 km. The offshore ergodic GMM and HTZ20 show great similarity in amplitude and trend, especially for the unburied states. For interface events, the offshore ergodic GMM are still dramatically different from the NGA-sub GMMs, with differences similar to the response spectra at 75 km. For intraslab events, the comparison results are similar to those of the response spectra at a distance of 75 km. These comparisons demonstrate that offshore GMMs are significantly different from onshore GMMs for subduction earthquakes, especially in the long-period and unburied states. Furthermore, the larger period of the peak response spectra of offshore GMMs compared to onshore GMMs implies that offshore ground motions are characterized by rich long-period ground motion.

Figure 17.

The comparison results of response spectra for interface (left panel) and intraslab (right panel) events with magnitude 7 at a distance of ((a) and (b)) 75 km and ((c) and (d)) 200 km with a reference site category SC2. For NGA-Sub GMMs, Ztor = 10 km for interface events and Ztor = 50 km for intraslab events; hypocentral depth (Z_HYP) is 30 km for interface events and hypocentral depth is 50 km for intraslab events; V_s30 = 500 m/s.

The significant differences in attenuation curves and response spectra between offshore and onshore models highlight the impracticality of directly applying onshore GMMs to offshore areas. Therefore, it is necessary to develop reasonable prediction models for offshore ground motion based on offshore ground motion observed data, particularly for subduction earthquake events that can cause significant damage. The large differences between onshore and offshore GMMs may be attributed to various factors, such as the deployment method of the seafloor stations, the degree of coupling between the stations and the seafloor, the topography of the seafloor, the thickness of the seafloor sediment, and the depth of seawater (Cui et al., 2023; Dhakal et al., 2021; Tan and Hu, 2023). Furthermore, the significant differences between buried and unburied stations suggest that the attenuation curves and response spectra of unburied stations are more susceptible to the influence of complex seafloor environments. In comparison, offshore buried stations increase the coupling with the seafloor, resulting in less difference in attenuation curves and response spectra between buried stations and onshore stations. However, some differences still exist due to the influence of seafloor topography and other factors.

Development of offshore non-ergodic GMM for subduction earthquakes

The development of non-ergodic GMM is based on an ergodic GMM to account for systematic and repeatable source, site, and path effects. We briefly describe the development of the non-ergodic GMM, using functional forms similar to Abrahamson et al. (2019), Kuehn et al. (2019), and Macedo and Liu (2022). We use the offshore ergodic backbone GMM with adjustment terms to describe the offshore non-ergodic GMM, which can be written as:

$${\log }_{10}A(T)={\mu }_{\mathrm{ergodic}}(x)+\delta c_0+\delta L2L\left(t_e\right)+\delta S2S\left(t_s\right)+\delta P2P(t_e,t_s)+\delta B_e^0+\delta W_{\mathrm{es}}^0$$

(10)

$$\delta L2L\left(t_e\right)=\delta c_{\mathrm{eq}}\left(t_e\right)$$

(11)

$$\delta S2S\left(t_s\right)=\delta c_{\mathrm{stat}}\left(t_s\right)+\delta S2S_0$$

(12)

$$\delta P2P(t_e,t_s)=\sum _1^n{c}_i{R}_i-{\theta }_{\mathrm{atten}}R$$

(13)

where μ_ergdic(x) represents the median value of ground motion IMs estimated from the ergodic backbone GMM (e.g. Equation 1). δc₀ is the constant intercept indicating the change in weight of each record due to the non-ergodic model considering spatial correlation. δL2 L(t_e), δS2 S(t_s) and δP2 P(t_e, ts)) are systematic and repeatable source, site and path effects, respcetively, which are epistemic uncertainties that cannot be represented by ergodic GMM; and t_e and t_s are the earthquake (epicenter) and station coordinates in the geography. δc_eq(t_e) term is the spatially varying source-specific adjustment term, which is intended to capture spatially correlated and repeatable source effects. δc_stat(t_s) is the spatially varying site adjustment term, which is intended to capture the repeatable regional site effects. δS2S₀ is the spatially location-specific site adjustment, which captures the repeatable location-specific site effects. δB⁰ _e , δW⁰ _es , and δS2S₀ are the inter-event, intra-event, and site-to-site residuals of non-ergodic GMM, which conform to the normal distribution of N(0, τ₀²), N(0, ϕ₀²), and N(0, ϕs₂s₀²), respectively. θ_attenR is the anelastic attenuation function in the offshore ergodic backbone GMM; $\sum _1^n{c}_i{R}_i$ is the cell-specific anelastic attenuation which can be seen as a alternative term for the anelastic attenuation function. c_i is the anelastic attenuation coefficient for the ith cell, which should be negative to avoid nonphysical attenuation coefficients (Abrahamson et al., 2019; Lavrentiadis et al., 2021). R_i is the length of the path within the ith cell. n is the total number of cells in the grid.

The cell-specific anelastic attenuation method

The cell size of 20 km × 20 km provides a good coverage for the Japan Trench region using the cell-specific anelastic attenuation method proposed by Dawood and Rodriguez-Marek (2013). The Cohen–Sutherland algorithm is used to calculate the path length R_i through each cell, which is a line clipping algorithm used in computer graphics (Foley et al., 1994; Liu et al., 2022a; Macedo and Liu, 2022). The algorithm was first proposed by Foley et al. (1994) and then extended by Liu et al. (2022a). The algorithm effectively determines whether a specific propagation path passes through a cell, avoiding unnecessary R_i calculations when the path does not pass through certain cells. Figure 18a shows the schematic diagram of the propagation path of seismic wave from the epicenter to the station (ray path diagram (green line)). A total of 2052 cells of 20 km × 20 km size was covered in the Japan Trench region using the cell-specific anelastic attenuation method, of which 931 cells contained propagation paths as shown in Figure 18b. In addition, the number of paths in cells contains at most 1285 paths and at least 1 path.

Figure 18.

(a) The schematic diagram of the propagation path and (b) the generated grid cells covered the Japan Trench region and the number of paths passed by each cell. The epicenters and stations are shown with red dots and blue triangles, respectively.

Bayesian inference with the integrated nested Laplace approximation

Bayesian inference provides a method to solve the uncertainty of the model by using probability theory, which can model the uncertainty of the non-ergodic GMM as probability distribution (Liu et al., 2022a, 2022b; Macedo and Liu, 2022). The integrated nested Laplace approximation (INLA) method can be implemented using the R-INLA package (www.r-inla.org) in the R computing environment (Lindgren and Rue, 2015; Rue et al., 2017). The INLA method, proposed by Rue et al. (2009), is a computationally effective alternative to Markov chain Monte Carlo designed for latent Gaussian models (Lindgren and Rue, 2015). For an introduction to INLA theory and example models, you can refer to works such as Blangiardo and Cameletti (2015), Wang et al. (2018), Krainski et al. (2019), and Gómez-Rubio (2020). The specific process of INLA is as follows: (1) a latent Gaussian random field (GRF) is used for modeling in a region of interest. (2) The finite-element method is used to divide the spatial domain into a group of disjoint triangles, and the latent GRF is represented as the latent Gaussian Markov random field (GMRF). (3) The latent GMRF with Matérn covariance matrices is used as solutions to stochastic partial differential equations (SPDEs) (Lindgren and Rue, 2015).

The spatially varying source and site adjustment (δc_eq(t_e) and δc_stat(t_s)) in the non-ergodic GMM are modeled as a Gaussian process with mean zero and Matérn covariance function (also known as the kernel function) (Landwehr et al., 2016; Macedo and Liu, 2022), as follows:

$$\delta c_{\mathrm{eq}}\left(t_e\right)~\mathrm{GP}(0,k(t_e,{t_e}^{\prime }))$$

(14)

$$\delta c_{\mathrm{stat}}\left(t_s\right)~\mathrm{GP}(0,k(t_s,{t_s}^{\prime }))$$

(15)

$$k(t,t\prime )=\frac{{\psi }^2}{2^{\upsilon -1}\Gamma \left(\upsilon \right)}{(\kappa \Vert t-t\prime \Vert )}^v{K}_{\upsilon }(\kappa \Vert t-t\prime \Vert )$$

(16)

where k(·,·) is a covariance function that determines the spatial correlation of systematic and repeatable source and site effects at any two locations; ψ is a hyperparameter that represents the marginal SD of the space field, Г(υ) is the gamma function, K_υ is the modified Bessel function of the second kind, κ is a scale parameter and υ is a smoothness parameter. $\rho =\sqrt{8\nu }/\kappa$ is another hyperparameter, which represents the length scale associated with the spatially varying term. The value of υ determines the mean square differentiability of the process, which is usually fixed because it is difficult to identify in applications. For GP with two-dimensional coordinates, the INLA usually takes the value of υ as 1, which is consistent with previous papers (e.g. Bakka et al., 2018; Li et al., 2011; Lindgren and Rue, 2015; Rue et al., 2009). Figure 19 shows the triangular mesh generated by INLA for constructing the SPDE model using the finite-element method, which covers the observed epicenter and station locations in the Japan Trench region. To account for the anelastic attenuation term of the non-ergodic GMM using modeling in km, we convert geographic coordinate systems with longitude and latitude values to objects with UTM projections. The mesh contains 5443 vertices and 10,762 triangles with a maximum edge length of 10 km and an outer boundary of 50 km. In general, the complexity of the model increases as the number of vertices increases (and the edge length decreases). Grids with too large edge lengths can lead to biases in inference (for advice on grid generation, please refer to the website). ²

Figure 19.

The triangular grid generated by the integrated nested Laplace approximation (INLA) method. The epicenters and stations are shown with red dots and blue triangles, respectively.

INLA is a Bayesian method, which means that a prior distribution for the estimated parameters needs to be set first. The parameters include: fixed effects (δc₀); SDs of offshore non-ergodic GMM (τ₀, ϕ₀ and ϕs₂s); and hyperparameters for the spatially varying term (e.g. spatial range ρ and SD ψ). To represent the prior distribution of fixed effects, the following normal distribution is used:

$$\delta c_0~N(0,0.2)$$

The prior distributions of the SD τ₀, ϕ₀, and ϕs₂s₀ use a log-gamma distribution with a shape parameter of 0.9 and a rate parameter of 0.01:

$$\begin{array}{c}\ln ({\tau }_0)~\log -\mathrm{Gamma}(0.9,0.01)\hfill \\ \ln (f_0)~\log -\mathrm{Gamma}(0.9,0.01)\hfill \\ \ln (f{s}_2{s}_0)~\log -\mathrm{Gamma}(0.9,0.01)\hfill \end{array}$$

The prior distribution of hyperparameters for spatially varying terms is based on the penalized complexity (PC) prior proposed by Fuglstad et al. (2019) for the SPDE model. The PC prior aims to shrink the parameters toward a simpler base model where the spatially varying terms have zero marginal variance and infinite range. For further details, you can refer to works such as Fuglstad et al. (2019) and Franco-Villoria et al. (2019). The spatial range ρ of the spatially varying term is specified with a prior probability of 50% for a spatially varying term smaller than 30 km; the SD ψ of the spatially varying term is specified with a prior probability of 1% for a spatially varying term larger than 0.3:

$$\begin{array}{l}P({\rho }_{eq/stat}<30)=50\%\\ P({\psi }_{eq/stat}>0.3)=1\%\end{array}$$

In addition, the scale of anelastic attenuation coefficient c_i is very different, and we also use the PC prior for the prior distribution of its SD:

$$P(\psi >0.01)=1\%$$

Results of offshore non-ergodic GMM

Spatially varying coefficient term

Figures 20 and 21 show maps of spatially varying source and site coefficient term (ie., δc_eq and δc_stat values) calculated by the INLA method, and the output results contain the mean, quartile values, and SD of the marginal posterior distribution. Positive values in the map indicate that the systematic and repeatable source and site effects in the region are higher than the average level; negative values indicate that the source and site effects are lower than the average level.

Figure 20.

(a) the mean values of source effects for PGA; (b) the standard deviations of source effects for PGA; (c) the mean values of source effects for Sa_RotD50 (0.5 s); (d) the standard deviations of source effects for Sa_RotD50 (0.5 s); (e) the mean values of source effects for Sa_RotD50 (1 s); (f) the standard deviations of source effects for Sa_RotD50 (1 s); (g) the mean values of source effects for Sa_RotD50 (8 s); (h) the standard deviations of source effects for Sa_RotD50 (8 s). In plots, the red dots represent the epicenters.

Figure 21.

The mean and standard deviation (SD) of the spatially varying site effects for PGA and spectral accelerations at T = 0.5, 1.0, and 5.0 s: (a) the mean values of site effects for PGA; (b) the standard deviations of site effects for PGA; (c) the mean values of site effects for Sa_RotD50 (0.5 s); (d) the standard deviations of site effects for Sa_RotD50 (0.5 s); (e) the mean values of site effects for Sa_RotD50 (1 s); (f) the standard deviations of site effects for Sa_RotD50 (1 s); (g) the mean values of site effects for Sa_RotD50 (5 s); and (h) the standard deviations of site effects for Sa_RotD50 (5 s). In plots, the blue triangles represent the distribution of stations.

Figure 20 shows maps of the mean and SD of the spatially varying source effects (δc_eq) for PGA and spectral accelerations at T = 0.5, 1.0, and 8.0 s. Figure 20a, c, e and g reveal that the positive values of the spatially varying source effects are predominantly concentrated on the side adjacent to the coast, while the negative values are primarily dispersed on the offshore side. Moreover, both positive and negative values are clustered near fixed locations and do not vary greatly with the period. This indicates that the subduction earthquakes happened on the coastal side are endowed with a greater potential for destruction in the Japan trench area. The INLA method takes the spatial correlation length into account and controls the attenuation of the correlation with the separation distance. Therefore, the map shows non-zero mean values for some locations where no earthquakes occurred. In addition, the spatial pattern exhibited by the spatial varying source effects indicates that the spatial correlation length varies greatly in different periods. The Table 3 shows the values of hyperparameter ρ for different periods controlling the range of spatial correlation length of spatially varying terms. Figure 20 and Table 3 show a lager spatial correlation length of spatially varying source effect at shorter periods than at longer periods. The spatially varying source effects exhibit different spatial patterns due to the nature of the ground motion parameters. The short-period parameters are sensitive to high-frequency ground motion, while the long-period parameters are sensitive to low-frequency ground motion. This difference in frequency content can lead to different spatial patterns of spatially varying source effects. The shorter correlation length for long-period parameters might be related to the attenuation of high-frequency ground motion with distance. High-frequency ground motion attenuates more rapidly than low-frequency ground motion, which can result in shorter correlation lengths for long-period parameters. This is because the energy of high-frequency waves is absorbed more quickly by the Earth’s crust, resulting in more localized effects. Figures 20b, d, f, and h indicates that the locations where historical earthquakes occurred exhibit a smaller SD. Conversely, the locations where no historical earthquakes occurred display near-zero mean values and large SDs, implying that there is insufficient data to constrain the systematic and repeatable source effect.

Table 3.

The mean values of the posterior distributions of the hyperparameters and the inter-event, intra-event, site-to-site, and total aleatory standard deviations for the non-ergodic GMM

Period(s)	δc0	ρeq (km)	ρstat (km)	ψeq	ψstat	τ0	ϕ0	ϕs2s0	σ0
0	0.002	126.22	85.59	0.173	0.129	0.135	0.177	0.067	0.232
0.01	0.003	126.17	85.62	0.173	0.129	0.135	0.177	0.068	0.232
0.02	0.003	126.12	85.57	0.173	0.129	0.135	0.177	0.068	0.232
0.03	0.004	126.60	85.90	0.173	0.130	0.134	0.177	0.068	0.232
0.05	0.011	302.11	92.60	0.190	0.137	0.137	0.176	0.081	0.237
0.08	0.030	207.86	103.18	0.203	0.158	0.137	0.174	0.083	0.236
0.1	0.026	215.13	87.73	0.212	0.164	0.138	0.173	0.081	0.236
0.15	0.007	174.48	90.50	0.203	0.182	0.136	0.174	0.108	0.246
0.2	−0.005	98.50	85.91	0.172	0.147	0.134	0.180	0.129	0.258
0.25	−0.026	106.52	79.93	0.164	0.138	0.139	0.184	0.114	0.257
0.3	−0.023	197.65	68.88	0.173	0.147	0.132	0.187	0.074	0.240
0.4	−0.015	201.45	65.79	0.164	0.141	0.136	0.185	0.079	0.243
0.5	−0.013	213.45	62.73	0.159	0.131	0.139	0.185	0.083	0.246
0.6	−0.019	242.97	59.49	0.165	0.134	0.139	0.183	0.078	0.243
0.75	−0.023	259.03	51.18	0.159	0.133	0.140	0.181	0.076	0.241
1	−0.007	130.01	73.33	0.143	0.126	0.141	0.176	0.086	0.241
1.5	−0.041	88.89	83.53	0.136	0.119	0.123	0.167	0.116	0.238
2	−0.062	59.92	54.82	0.141	0.137	0.111	0.160	0.070	0.207
2.5	−0.084	80.16	74.71	0.125	0.110	0.114	0.153	0.107	0.219
3	−0.094	82.99	82.14	0.086	0.085	0.117	0.147	0.112	0.218
4	−0.078	61.47	58.12	0.079	0.079	0.118	0.146	0.092	0.209
5	−0.068	61.00	47.68	0.108	0.100	0.126	0.149	0.067	0.206
6	−0.053	56.64	53.66	0.088	0.087	0.124	0.153	0.088	0.216
8	−0.035	58.70	50.24	0.104	0.096	0.123	0.164	0.059	0.213
10	0.000	52.14	48.34	0.106	0.088	0.116	0.171	0.064	0.216

Figure 21 shows maps of the mean and SD of the spatially varying site effects (δc_stat) for PGA and spectral accelerations at T = 0.5, 1.0, and 8.0 s. Figure 21a, c, e, and g show that the positive and negative values of the spatially varying site effects change as the period changes. The negative values of spatially varying site effects are distributed mainly near the coast at short periods and transition to be distributed near the trench as the period increases. Conversely, the positive values of spatially varying site effects are distributed mainly near the trench at short periods and transition to be distributed near the coast as the period increases. Figure 21b, d, f, and h shows the SDs of the spatially varying site effects, which is similar to the spatially varying source effects. In regions without station distribution, the site effects have near-zero means and large SDs. The spatial pattern exhibited by the spatially varying site effects indicates that the spatial correlation length does not change much in different periods. The hyperparameter ρ that control the range of spatial correlation length of spatial varying terms is shown in Table 3. The reason why the spatial correlation length of the spatially varying site effects remain relatively stable may be attributed to many factors, such as station density (ie., the 150 offshore stations are uniformly distributed near the Japan Trench), consistency of site conditions and geological stability (ie., each station maintains a stable site condition at a fixed location for different periods).

These non-zero means of spatially varying source and site effect coefficients represent source and site characteristics that cannot be captured by the ergodic backbone GMM. However, these systematic and repeatable source and site effects can be constrained using the ground motion records in the non-ergodic GMM.

Figure 22 shows the estimated mean and SD of the site-specific adjustment term δS2S₀ for PGA. The mean values of δS2S₀ do not exhibit spatial correlation (i.e. the spatial correlation length of δS2S₀ is zero, and there are no systematic positive or negative regions), indicating that spatially varying site effects δc_stat correctly capture the spatial correlation component of the site effects.

Figure 22.

The estimated (a) mean and (b) standard deviation of the site-specific adjustment term δS2S₀ for PGA.

The cell-specific attenuation coefficient

Figure 23 shows the mean and epistemic SD of the cell-specific anelastic attenuation coefficient c_i of each cell for PGA and Sa_RotD50 (1 s). The cell-specific anelastic attenuation coefficients represent systematic and repeatable path effects. For cells without paths, the mean cell-specific anelastic attenuation coefficients c_i are considered to be consistent with the anelastic attenuation coefficient k(T) of the ergodic backbone GMM. Furthermore, most of the cell-specific anelastic attenuation coefficients are consistent with the anelastic attenuation coefficient k(T) of the ergodic GMM. However, some cell-specific anelastic attenuation coefficients exhibit nonphysical behavior, with a small amount of positive cell-specific attenuation coefficients in the region. Figure 23c shows nonphysical behavior for Sa_RotD50 (1 s) in the region of buried stations due to the increased ground coupling of buried stations. Cells containing a large number of paths have smaller epistemic SDs, while cells containing fewer paths have larger epistemic SDs. The large epistemic SD implies that the current data set is insufficient to provide a good estimate of the cell-specific anelastic attenuation coefficient in this cell.

Figure 23.

The mean anelastic attenuation coefficient and epistemic standard deviation of each cell estimated by the cell-specific anelastic attenuation method: (a) the mean cell-specific anelastic attenuation coefficients for PGA; (b) the epistemic standard deviations for PGA; (c) the mean cell-specific anelastic attenuation coefficients for Sa_RotD50 (1 s); and (d) the epistemic standard deviations for Sa_RotD50 (1 s).

Residuals and standard deviations

Figure 24 shows residual analysis of non-ergodic GMM with respect to moment magnitude (M), rupture distance (R_rup), thickness of sediment up to Vs = 1.4 km/s (D14), and water column for PGA and spectral accelerations at T = 0.5 and 1.0 s. Compared with Figure 9, the distribution of residuals with respect to each explanatory variable (ie., M, R_rup, D14, and water column) has less scattering and become more concentrated near zero. This is due to the fact that the non-ergodic GMM incorporates repeatable and systematic source, site, and path effects into the ergodic GMM, which makes the non-ergodic GMM have less epistemic uncertainty.

Figure 24.

Residual analysis of non-ergodic GMM with respect to moment magnitude (M), rupture distance (R_rup), thickness of sediment up to Vs = 1.4 km/s (D14), and water column for PGA and spectral accelerations at T = 0.5 and 1.0 s. (a1)–(a3) the inter-event residuals versus moment magnitude; (b1)–(b3) the intra-event residuals versus moment magnitude; (c1)–(c4) the intra-event residuals versus rupture distance; (d1)–(d3) the intra-event residuals versus thickness of sediment up to Vs = 1.4 km/s; (e1)–(e3) the intra-event residuals versus water column. Orange circles and green error bars represent the mean and standard deviation of the residuals, respectively.

Figure 25 shows the inter-event, intra-event, and total aleatory SD terms of the ergodic and non-ergodic GMM for each period. The non-ergodic GMM exhibits a significant reduction in aleatory variability, with the values of the total aleatory SDs reduced by about 22%–39%. This reduced aleatory variability is important for explaining the epistemic uncertainty in the non-ergodic terms and for performing non-ergodic PSHA. The mean values of the posterior distributions of the hyperparameters, as well as the inter-event, intra-event, site-to-site, and total aleatory SDs for the non-ergodic GMM, are listed in Table 3 for periods from 0 to 10 s. The hyperparameters of spatially varying source effects vary greatly, ranging from 52 to 302 km in different periods; while the hyperparameters of spatially varying site effects vary slightly, ranging from 47 to 103 km.

Figure 25.

The inter-event, intra-event, and total aleatory standard deviation terms of ergodic backbone and non-ergodic GMM for each period.

Conclusion and discussion

An offshore ergodic backbone GMM for subduction earthquakes in the Japan Trench region was developed in this study, and its functional form was implemented by modifying the GMM of Si et al. (2022). Site amplification was investigated by classifying S-net stations into four site categories using horizontal-to-vertical response spectral ratio (HVRSR). The amplitudes and shapes of the mean HVRSR of four site categories were significantly different, especially for unburied stations. The offshore ergodic GMM can accurately capture the observed offshore ground motion in terms of distance scaling for different subduction earthquake scenarios. It is applicable for subduction earthquake scenarios in the Japan Trench area with moment magnitudes ranging from 4.0 to 7.4 and rupture distances ranging from 10 to 300 km. In addition, we compared the offshore ergodic GMM with NGA-Sub GMMs and some GMMs from Japan in terms of attenuation curves and response spectra. The results showed that offshore GMMs are significantly different from onshore GMMs for subduction earthquakes, especially in the long-period and unburied states, which highlights the impracticality of directly applying onshore GMMs to offshore regions. The large differences between onshore and offshore GMMs may be related to the deployment method of the seafloor stations, the degree of coupling between the stations and the seafloor, the topography of the seafloor, the thickness of the seafloor sediment, and the depth of seawater. The significant differences between buried and unburied stations indicated that the unburied stations are more susceptible to the influence of complex seafloor environments. Therefore, it is necessary to develop reasonable prediction models for offshore ground motion, particularly for subduction earthquakes that can cause significant damage.

The non-ergodic GMM for subduction earthquakes was developed based on the offshore ergodic backbone GMM. In the non-ergodic GMM, the systematic and repeatable source and site effects were captured by the spatially varying coefficients represented by Gaussian processes; the systematic and repeatable path effects were captured by cell-specific anelastic attenuation proposed by Dawood and Rodriguez-Marek (2013), calculated with the Cohen–Sutherland computer graphics algorithm. The spatial pattern exhibited by the spatial varying source effects varied slight for different periods. The positive values of the spatially varying source effects were mainly concentrated on the side near the coast, while the negative values were mainly scattered on the offshore side. However, the spatial pattern exhibited by the spatial varying site effects varied continuously with the period. The negative values of site effects were mainly distributed near the coast in short periods, transitioning to be distributed near the trench in long periods. Conversely, the positive values of site effects were mainly distributed near the trench in short periods, transitioning to be distributed near the coast in long periods. The spatial varying source and site effects had near-zero means and large epistemic SDs for regions without epicenter and station distributions, implying that there are not enough data to constrain the systematic and repeatable source and site effects. The mean anelastic attenuation coefficient and epistemic SD of each cell calculated by Bayesian inference with INLA can represent systematic and repeatable path effects. Cells containing a small number of paths had larger epistemic SD. The epistemic uncertainty of spatially varying effects is consistent with intuitive data availability. This implies that the epistemic uncertainty of estimated spatially varying effects is greater for regions without sufficient data.

Compared to the residuals of ergodic GMM, the distribution of residuals of non-ergodic GMM with respect to each explanatory variable had less scattering and was more concentrated around zero. This is because the non-ergodic GMM incorporates systematic and repeatable source, path and site effects into GMM. In addition, the values of the total aleatory SDs of non-ergodic GMM were reduced by about 22%–39%. The non-ergodic GMM developed in this study can better capture the aleatory variability and epistemic uncertainty of ground motion in the Japan trench region. Such aleatory variability and epistemic uncertainty should be incorporated into PSHA calculations because a small variation in SDs can lead to an underestimation of seismic risk. These results can provide a basis for rational seismic design and seismic risk assessment of marine engineering in offshore subduction zone, as well as the prediction of ground motions for EEW and seismic demand analysis of offshore facilities in the Japan Trench region.