Linearly combined ground motion model using quadratic programming for low- to mid-size seismicity region: South Korea

This study suggests a linearly combined ground motion model (GMM) optimized for the within-rock conditions of South Korea. Ground motions recorded by accelerometers positioned in within-rock layers are used as target intensity measures (IMs), and five local GMMs and four global GMMs are used as candidates for combination. Optimization, which seeks to find a weight vector minimizing the uncertainty of the combined model, is performed using the quadratic programming (QP) technique, which provides very fast and solid results for the linear combination problem. This study illustrates how to use the QP technique for the linear combination problems. Also, we suggest optimized weight vectors for GMM combinations for two conditions: 1) the IM prediction of a scenario event without observations and 2) the IM prediction of a past event with observations. Among the local and global GMMs considered, JB03, Eea15, JH21, and BSSA14 are selected as the best four GMMs for the first condition, and Jea02, JB03, JH21, and CB14 are selected as the best four GMMs for the second condition. The combined model reduces the standard deviation of residuals in natural logarithms by 10% and 8% for the first and second conditions, respectively, compared to the best individual GMM at each period. Among the GMMs considered, the prediction by Eea15 is only applicable for magnitudes less than 5. Hence, for large magnitude (Mw > 5) prediction, CB14 is recommended instead of Eea15 for the best four models for the first condition.

1 Introduction

After two damaging earthquakes, 2016 M_L 5.8 earthquake in Gyeongju and 2017 M_L 5.4 earthquake in Pohang, there is a strong need to decide which ground motion model (GMM) or what combination of GMMs is best for the engineering problems (e.g., earthquake early warning and the post-earthquake intensity distribution map). For the development of region-specific GMMs in South Korea, a stochastic method simulating ground motions for the moment magnitude (M_W) range of 4–7 and regressing a GMM using simulated intensity measures (IMs) was popular (Noh and Lee, 1995; Jo and Baag, 2001; Jo and Baag, 2003; Park et al., 2001; Junn et al., 2002). There are only limited ground motion records with less than mid-size events (M_W ≤ 5.4) in South Korea, so the development of the simulation-based GMM was a reasonable approach to consider the large M_W range. Cumulating ground motion records over time with an increased number of seismic stations, some GMMs have been developed empirically (Emolo et al., 2015), but still, the largest M_W is low to be considered for the seismic design practice. A hybrid approach, combining empirical data set and ground motion from stochastic simulation, was another way for the GMM development (Jee and Han, 2021).

In the probabilistic seismic hazard assessment, multiple GMMs are combined using a logic-tree approach that assigns appropriate weights to each GMM and linearly combines them (Kramer, 1996). The linear combination of multiple prediction models with appropriate weights reduces the uncertainty of model estimation (Kwak et al., 2018). These weights could be defined based on the rank of GMMs using several ranking procedures (Scherbaum et al., 2009; Kale and Akkar, 2013; Akkar and Kale, 2014; Nizamani and Park, 2021). The ranking provides the order of preference, but it does not suggest the most efficient weights for minimizing the standard deviation of predictions of the combined model. For example, Scherbaum et al. (2009) suggest weights based on the amount of the log-likelihood, but this does not necessarily provide the minimum standard deviation of predictions.

This study suggests the best combination of GMMs with optimized weights minimizing the standard deviation of predictions using the quadratic programming (QP) optimization technique for the IM of ground motion in South Korea. The optimized weights for multiple models were previously found using the Monte Carlo (MC) simulation (Kwak et al., 2018), but the MC simulation is not efficient in terms of execution time and repeatability of the results. QP solves an optimization problem involving a quadratic function (Goldfarb and Idnani, 1983). The QP technique is adopted in this study to find the optimized weights by adjusting the inequality constraints to satisfy the weight conditions of linear combination. The QP is much faster than the MC simulation and provides a solid result.

For the linear combination, five local (JB01, Jo and Baag, 2001; Jea02, Junn et al., 2002; Yea08, Yun et al., 2008; Eea15, Emolo et al., 2015; and JH21, Jee and Han, 2021) and four global (ASK14, Abrahamson et al., 2014; BSSA14, Boore et al., 2014; CB14, Campbell and Bozorgnia, 2014; and CY14, Chiou and Youngs, 2014) GMMs are considered as candidates. The local GMMs were used for the seismic hazard map development (Hong, 2017) and to predict spectral acceleration up to 5 s. To consider IMs for large M_W events for the combined GMM, which are out of range for ground motion records in Korea, we also considered four global GMMs from the NGA West2 project (Bozorgnia et al., 2014). We tested the performance of the combined model by comparing it to the seismic records observed during 2000–2018 by the National Seismic Network operated by the Korea Meteorological Administration (KMA).

This study suggests two optimized weight sets: one for a scenario event (where it is anticipated that no records are available) and the other for past events (where it is anticipated that records are available). It should be noted that both weight sets are estimated using ground motion records observed in South Korea. For the scenario event, weights are found based on the total residuals of GMMs. If records are available, the between-event residuals can be calculated, which makes it possible to compute the within-event residuals. Therefore, for the past event, the optimized weight set is suggested, minimizing standard deviation of within-event residuals.

The subsequent sections describe the procedure for solving the weight optimization problem for linear combination, introduce the ground motion data set from South Korea and GMMs, and evaluate the performance of the combined model. Suggestion and conclusion are followed.

2 Quadratic programming defining optimized weights

This section illustrates the linear combination and how to adjust inequality constraints in the QP technique to determine optimized weights in the linear combination problem. This method becomes the base of GMM combinations and provides optimized weight vectors.

The weighted combination of n number of models is expressed as:

${\widehat{Y}}_{comb}={\displaystyle \sum }_{i=1}^n{w}_i{\widehat{Y}}_i,(1)$

where ${\widehat{Y}}_{comb}$ is the combined prediction, w_i is the weight, and ${\widehat{Y}}_i$ is the prediction by the ith model. w_i has the summation condition and lower and upper bound conditions as:

${\displaystyle \sum }_{i=1}^n{w}_i=1,(2)$

$0\le w_i\le 1.(3)$

If we express the residual of the ith model as ε_i, then the residual of the ${\widehat{Y}}_{comb}$ (i.e., ε_comb) can be expressed as:

${\epsilon }_{comb}=Y-{\widehat{Y}}_{comb}=Y-{\displaystyle \sum }_{i=1}^N{w}_i{\widehat{Y}}_i=Y-{\displaystyle \sum }_{i=1}^N{w}_i\left(Y-{\epsilon }_i\right)={\displaystyle \sum }_{i=1}^N{w}_i{\epsilon }_i.(4)$

This indicates that ε_comb is the weighted mean of ε_i. Then, the variance of ε_comb (${\sigma }_{comb}^2$) can be expressed as:

${\sigma }_{comb}^2=Var\left({\epsilon }_{comb}\right)=Var\left({\displaystyle \sum }_{i=1}^N{w}_i{\epsilon }_i\right)=w_1^2{\sigma }_1^2+w_2^2{\sigma }_2^2+\cdots +2w_1{w}_2{\sigma }_{12}+\cdots ,(5)$

where ${\sigma }_i^2$ is the variance of ${\epsilon }_i$, and ${\sigma }_{ij}$ is the covariance of ${\epsilon }_i$ and ${\epsilon }_j$. Using the matrix notation, Eq. 5 becomes:

${\sigma }_{comb}^2=\left[w_1\cdots w_n\right]\left[\begin{array}{ccc}{\sigma }_1^2& \cdots & {\sigma }_{1n}\\ ⋮& \ddots & ⋮\\ {\sigma }_{n1}& \cdots & {\sigma }_n^2\end{array}\right]\left[\begin{array}{c}{w}_1\\ ⋮\\ w_n\end{array}\right]={\mathbf{w}}^{\mathrm{T}}{\mathbf{\sigma }}^2\mathbf{w},(6)$

where w is the w_i vector and σ² is the covariance matrix. w minimizing ${\sigma }_{comb}^2$ in Eq. 6 is called the optimized weight vector.

Equation 5 is a quadratic function. Hence, finding the object variable (i.e., w) by minimizing the value of Eq. 5 (i.e., ${\sigma }_{comb}^2$) can be solved using the QP optimization technique. QP solves an optimization problem involving the quadratic function, where the general form is expressed as follows (Goldfarb and Idnani, 1983):

$\begin{array}{c}\underset{x}{\min }:\frac{1}{2}{\mathbf{x}}^{\mathrm{T}}\mathbf{G}\mathbf{x}+{\mathbf{a}}^{\mathrm{T}}\mathbf{x}\\ s.t.:{\mathbf{C}}^{\mathrm{T}}\mathbf{x}\ge \mathbf{b}\end{array},(7)$

where x and a are n vectors, G is an n×n symmetric positive definite matrix, C is an n×m matrix, and b is an m vector. Equation 7 contains the quadratic part ($\frac{1}{2}{\mathbf{x}}^{\mathrm{T}}\mathbf{G}\mathbf{x}$), the linear part (${\mathbf{a}}^{\mathrm{T}}\mathbf{x}$), and the inequality constraints (${\mathbf{C}}^{\mathrm{T}}\mathbf{x}\ge \mathbf{b}$). Goldfarb and Idnani (1983) suggested a numerical solution for finding the optimized x, minimizing the function value. We can use this solution to find optimized w_i, replacing variables in Eq. 7 to be compatible with Eq. 6 and constraints to satisfy the weight conditions (Eq. 2 and Eq. 3). The adjusting process is shown as follows:

1) Because Eq. 6 does not have a linear part, ${\mathbf{a}}^{\mathbf{T}}\mathbf{x}$ is ignored by setting a as a zero vector;

2) x and 0.5 G were replaced as w and σ², respectively;

3) The C^T matrix is set as follows (using three models, n = 3 as an example):

${\mathbf{C}}^{\mathrm{T}}=\left[\begin{array}{c}\begin{array}{ccc}1& 1& 1\end{array}\\ \begin{array}{ccc}1& 0& 0\end{array}\\ \begin{array}{ccc}0& 1& 0\end{array}\\ \begin{array}{ccc}0& 0& 1\end{array}\end{array}\right].(8)$

4) The b vector is set as follows (using three models, n = 3 as an example):

$\mathbf{b}=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right].(9)$

It should be noted that m becomes n+1 in Eqs 8, 9. For general n, the first row must be set as 1 and an n×n identify matrix from the second for C^T, and the first row as 1 and 0 for the last for b. By adjusting C^T and b as mentioned previously, Eq. 7 becomes (n = 3 as an example):

$\begin{array}{c}\underset{w}{\min }:{\mathbf{w}}^{\mathbf{T}}{\mathbf{\sigma }}^{\mathbf{2}}\mathbf{w}\\ s.t.:\left[\begin{array}{c}\begin{array}{ccc}1& 1& 1\end{array}\\ \begin{array}{ccc}1& 0& 0\end{array}\\ \begin{array}{ccc}0& 1& 0\end{array}\\ \begin{array}{ccc}0& 0& 1\end{array}\end{array}\right]\left[\begin{array}{c}{w}_1\\ w_2\\ w_3\end{array}\right]\ge \left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right].\end{array}(10)$

The inequality constraint of Eq. 10 indicates:

${\displaystyle \sum }_{i=1}^n{w}_i\ge 1,(11)$

$w_i\ge 0,i=\mathrm{1,2,3}.(12)$

Since ${\sigma }_i^2$ is a positive real number, to minimize w^Tσ²w in Eq. 10, the sum of w_i in Eq. 11 should be minimized. This means that the optimization solution of Eq. 10 results in the sum of w_i as unity, which satisfies the weight condition of Eq. 2. Also, w_i is greater than zero due to Eq. 12 constraint and the sum is unity, so w_i cannot be greater than 1, which satisfies the other weight condition in Eq. 3.

The QP algorithm is available in several programming languages (e.g., MATLAB, Python, or R). We used the function solve.QP in the quadprog library in R to find the optimized weights (R Core Team, 2022). QP provides the unique w resulting in the smaller ${\sigma }_{comb}^2$ than the solution found from the Monte Carlo simulation approach (Kwak et al., 2018), and the computing speed is much faster. We applied this QP technique to find the best combined GMM for within-rock conditions of South Korea in subsequent sections.

3 Ground motion records and models

3.1 Ground motion records

At the current time (September 2022), KMA maintains 261 seismic observatory stations, and there are 169 outdated stations that were replaced with new sensors or moved to other locations (Korea Meteorological Administration, 2022). Figure 1 shows the locations of all 430 KMA stations. For the stations, there are three stratigraphic levels of sensor locations: 1) surface (0 m depth; called first level hereafter); 2) depth at the top of the bedrock layer, which is defined as the soft rock layer or stiffer (20–70 m depth; called second level hereafter); and 3) 100 m depth (called third level hereafter). The KMA had installed sensors at the first level, but they have been switching the sensor locations to the second or third levels for the purpose of broadband detection of microseisms with minimized ambient noise. For the second-level location, the installation depth is intentionally kept at 20 m even though the bedrock depth is shallower than 20 m. Currently, there are 203 stations that only have sensors at downhole (second and/or third levels), 50 stations that have surface sensors only (first level), and eight stations that have both surface and downhole sensors.

FIGURE 1

FIGURE 1. Location of KMA seismic stations. Symbol shapes indicate the installation depth of sensors per station, and hollowness distinguishes the operation status.

Among three levels of sensor locations, our interest is in the second level. For practical purposes, predicting IMs at the surface (i.e., first level) is necessary; the seismic hazard analysis of superstructures or geotechnical failure events that occurred at the Earth’s surface (i.e., landslide, lateral spreading, liquefaction, and/or ground settlement) are related to the intensity of the ground motion at the surface. However, the major portion of the KMA sensors is at the second or third levels without sensors at the first level. Moreover, the KMA keeps switching the first level sensors to the second or third levels. This leads to a lack of first-level recordings, which limits the development of a GMM for the surface. To predict first-level IMs, alternatively, an amplification model from the second level to the first level can be coupled with the second-level IM prediction. The third level would have less variability of ground characteristics than the second level, so the GMM targeting the third level would result in less ground motion prediction uncertainty, but the number of records and stations for the third level are less than the second level, and availability of geophysical investigations down to 100 m depth is limited at this time. It should be noted that the data available for the second level are 3,001 records for 207 stations, and the data for the third level are 1,979 records for 65 stations. Hence, in this study, we targeted the second level and used records from the second-level sensors for the target data of the combined GMM.

For the local magnitude (M_L) range of 3.5–5.8, the epicentral distance (R_epi) range of <210 km, and the time interval between 2000 and 2018, 674 ground motion records have been collected at the second-level sensors. It should be noted that the official magnitude type of earthquake announcement in South Korea is M_L, which is used by the KMA. Among these records, we selected 614 ground motions by filtering out poorly recorded ground motions containing very high ambient noise or instrumental errors. We set the lower limit of the M_L as 3.5 because it is rare for an event with M_L < 3.5 to cause damage to structures, and M_L 3.5 is the KMA’s low threshold for the earthquake early warning service. The upper limit of the M_L is 5.8 because it is the greatest magnitude recorded in South Korea. The M_L versus R_epi distribution for the selected records is shown in Figure 2. The major portion of the data is located at the long R_epi (>100 km) and small M_L (<4).

FIGURE 2

FIGURE 2. Distribution of the local magnitude (M_L) and epicentral distance (R_epi) of ground motions recorded at the second-level sensors during 2000–2018 by the KMA network.

3.2 Ground motion models

Five region-specific GMMs for South Korea are summarized in Table 1. Three GMMs were developed using stochastic simulation (i.e., random vibration theory; Boore, 1983), one GMM was developed using the empirical data set, and one GMM was developed using the hybrid approach. The stochastic simulation method extracts the stress drop ($\Delta \sigma$) and Q factor from records, generates ground motions based on point-source event generation for variable magnitude and source-to-site distances, and finds regression models fitting generated IMs. The empirical method finds the best-fit regression model to the empirical IMs collected from records. The hybrid method combines simulated IMs and recorded IMs and finds the best-fit regression model. The $\Delta \sigma$ and Q factors used for the ground motion simulation per study, the dataset used, and the applicable ranges of magnitudes and source distances are shown in Table 1. All GMMs considered in this study have no site effect considerations for general sites, and the reference site condition is the hard-rock outcrop. It should be noted that simulation-based and hybrid GMMs used the moment magnitude (M_W) and hypocenter distance (R_hypo) or rupture distance (R_rup) as the model input parameters. The GMM based on the empirical data (Emolo et al., 2015) used M_L and R_epi as the model input parameters.

TABLE 1

TABLE 1. Comparison of parameters of five Korean-specific GMMs.

The ground motion records collected in this study have limited information: M_W or M_L, R_epi, and hypocentral depth (Z_hyp). The M_L is available for all ground motions, but M_W is only available for large events. To use records with M_L as an input parameter for GMM, we transformed M_L to M_W using the M_L–M_W relationship suggested for South Korea events (combination of Eq. 8 and Eq. 10 in Sheen et al. (2018)):

$M_W=\left(0.9234M_L+0.8034\right)/1.086.(13)$

It should be noted that the M_L in Eq. 13 is the local magnitude suggested by the KMA, and the KMA uses the vertical component for the M_L estimation.

Figure 3 compares the selected IMs with respect to varying M_W (3–7) at R_hypo = 10, 50, and 200 km. For Eea15, the maximum applicable M_L (=4.9) is less than other GMMs, so predictions at large M_W are excessively high. Except Eea15, all GMMs have similar magnitude scaling with the non-linear functional shape.

FIGURE 3

FIGURE 3. Comparison of five Korean GMMs (Jea02, JB03, Yea08, Eea15, and JH21) and four global GMMs (ASK14, BSSA14, CB14, and CY14) for PGA, SA_0.2, SA_1.0, and SA_2.0 with varying M_W at R_hypo=10, 50, and 200 km.

Four global GMMs from the NGA West2 project (ASK14, BSSA14, CB14, and CY14) are also considered as candidates for the model combination. NGA West2 GMMs require more detailed input parameters for source, path, and site effects (e.g., fault type, hanging wall effect, depth to the top of the fault, and region) than local GMMs. Table 2 shows the list of input parameters in the NGA West2 project GMMs. However, the data set collected for South Korea only has the parameters M_W, R_epi, and Z_hyp. Hence, we need to assume some parameters for the NGA West2 GMMs with default values. We set up input parameters as “unknown” if applicable, otherwise assigning general values. The fault dip angle, rake angle, and depth to the seismogenic crust (Z_bot) need default values, for which 90, 0, and 20 km are used, respectively. We assumed no hanging wall effect, site condition with V_S30 = 760 m/s, and an unknown basin depth term. The trends of IMs for NGA West2 GMMs are also shown in Figure 3. It should be noted that the global GMMs have distinct moment saturation points comparing to region-specific GMMs.

TABLE 2

TABLE 2. Variation of input parameters adopted for the model development in the NGA West2 project and parameters selected to use for Korean event prediction.

4 Residual analysis

4.1 Site condition adjustment

The target site conditions of five local GMMs are rock outcrop, and those of global GMMs are surfaces with varying V_S30. However, the target site condition in this study is the within-rock layer (i.e., second level), and within-rock layer records are collected. Hence, we expect that direct comparison of GMM predictions with the collected records is not appropriate—predictions (surface) would be larger than observations (within-rock). Figure 4 shows boxplots of total residuals (res_tt) per GMM. Generally, the medians (centerline of the box) of res_tt are smaller than zero for local GMMs for a wide range of periods, but those for global GMMs are period-dependent. For long periods (SA_1.0), medians are less than zero, but for low periods (PGA and SA_0.2), medians are greater than zero. This means that for local GMMs, most ground motion observations are lower than predictions as expected, but the global GMMs predict too low IMs for certain periods. This might be attributed to the fact that the global GMMs were developed using ground motion records from different regions.

FIGURE 4

FIGURE 4. Total residual boxplots of four IMs (PGA, SA_0.2, SA_1.0, and SA_2.0) binned by M_L for nine GMMs.

To remove this discrepancy arising from the difference of the target site condition and regional difference, we omitted the global bias of res_tt by performing mixed-effect linear regression. The mixed-effect linear regression provides random effects and fixed effects, where the fixed effects are equivalent to the global bias of res_tt, the random effects are between-event residuals, and the remaining are the within-event residuals (Abrahamson and Youngs, 1992; Boore et al., 2014). Using the mixed-effect linear regression, res_tt can be disintegrated as:

$re{s}_{tt}=c_k+\delta B_e+\delta W_{es},(14)$

where c_k is the global bias for model k, $\delta B_e$ is the between-event residual for event e, and $\delta W_{es}$ is the within-event residual for event e and site s. $\delta W_{es}$ includes path variability, site-to-site variability, and within-site variability (Al Atik et al., 2010). Table 3 provides c_k resulting from mixed-effect linear regression for selected IMs for five local and four global GMMs. Omitting c_k does not completely remove the effects of the actual target site condition discrepancy, but at least it makes an unbiased mean estimation for the target site condition.

TABLE 3

TABLE 3. Global bias of the mixed-effect linear regression for five local and four global GMMs.

4.2 Comparison of GMM performance

Figure 5 shows the standard deviation of res_tt (denoted σ_tt), $\delta B_e$ (denoted σ_b), and $\delta W_{es}$ (denoted σ_w) for PGA and spectral accelerations with varying periods (SA_period). For the σ_b, global GMMs perform better than local GMMs especially at short periods (<1 s). At long periods (>1 s), σ_b is similar among all GMMs except Eea15 which has a higher σ_b. Eea15 predicts biased estimation for records with M_L > 4.9, especially for long-period spectral ordinates, which is the reason for the higher σ_b. CY14 has the lowest σ_b in average. For the σ_w, local GMMs have lower σ_w for periods <1 s than global GMMs. Among local GMMs, JB03 has the lowest σ_w in average, and among global GMMs, BSSA14 has the lowest σ_w. For σ_tt, JH21 has the lowest σ_tt over a wide range of periods and JB03 is the second lowest. Local GMMs generally perform better than global GMMs, while Eea15 has high σ_tt at periods >0.5 s. Global GMMs have relatively high σ_tt at periods of 0.05–0.3 s, and among them, CY14 has the lowest σ_tt in average. Because global GMMs have less uncertainty about between-event residuals and local GMMs have less uncertainty about within-event residuals, global GMMs have better constraints for the source term, while local GMMs perform better for the path and site term in South Korea’s records.

FIGURE 5

FIGURE 5. Standard deviation of (A) total residuals, (B) between-event residuals, and (C) within-event residuals for nine individual GMMs and combined GMMs suggested in this study.

5 Performance of combined GMM

5.1 Selection of best combination

For the sake of the IM prediction for future earthquakes, we cannot specify $\delta B_e$. Hence, σ_tt combining the uncertainty of both $\delta B_e$ and $\delta W_{es}$ is the appropriate metric for the evaluation of the model’s performance. However, if we apply a GMM for an event with records available, we can calculate $\delta B_e$, that is, equivalent to the mean misfit of res_tt (Abrahamson and Youngs, 1992). In this case, $\delta B_e$ in Eq. 14 is defined, so lowering σ_w would only contribute to the reduction of the total uncertainty. GMMs are often used to evaluate the intensity distribution at the post-event stage, so it is also of interest to find the GMM combination to minimize σ_w. This study suggests two types of linear combinations: one for σ_tt minimization and the other for σ_w minimization.

In practice, for simplicity, using two to four GMMs is preferable for the combination rather than use of several to tens of GMMs. Hence, in this section, we examine the efficiency of the model combination using different numbers of GMMs. Among the nine GMMs (five local and four global), we picked one to five GMMs as combination candidates and found the optimized weights using the QP process. The QP process is summarized as follows:

1) The covariance matrix (σ² in Eq. 6 and Eq. 10) of the res_tt and δW_es between GMMs selected was calculated;

2) σ² is inserted into Eq. 10;

3) The optimized w is calculated.

It should be noted that the QP should have σ² as the positive definite matrix, so we slightly altered σ² to be the nearest positive definite matrix using the nearPD function in R (Higham, 2002) if σ² is not positively definite. Figure 6 shows the average of σ_tt and σ_w of IMs including PGA and 15 kinds of SAs ranging from 0.01 to 5 s of the best combined models for different number of model combination. The standard deviation is the highest for a single GMM case, and it decreases steeply for the case of a combination of two models. The standard deviations are decreasing with the greater number of GMM combinations, but the reduction rate slows down. Table 4 shows the GMMs selected for the best combination for different numbers of models and indicates the percentage difference in average standard deviation compared to the combination of all nine GMMs. The single best GMM has 10.2% larger σ_tt than one from a nine-GMM combination and 8.45% larger σ_w than one from nine GMM combination, but the difference becomes small if the number of GMMs for combination increases. The differences of average standard deviation between a four-GMM combination and a nine-GMM combination become negligible (less than 0.2%), so we selected a four-GMM combination as the suggested GMM model in this study.

FIGURE 6

FIGURE 6. Percentage difference of standard deviation between one from the combination of nine GMMs and the other from a combination of different numbers of GMMs. (A) σ_tt and (B) σ_b.

TABLE 4

TABLE 4. Best GMMs for different numbers of candidates.

Table 5 describes the optimized weights suggested for res_tt and $\delta W_{es}$. Three local GMMs (JB03, Eea15, and JH21) and one global GMM (BSSA14) are recommended for res_tt, and other three local GMMs (Jea02, JB03, and JH21) and one global GMM (CB14) are recommended for $\delta W_{es}$.

TABLE 5

TABLE 5. Suggested weights of four GMMs for res_tt and $\mathit{\delta }{\mathit{W}}_{\mathit{e}\mathit{s}}$ per IM.

5.2 Performance of combined GMM compared to records

To evaluate the performance of the combined GMM, we compared the median predictions of the GMMs with the recorded distribution. Figure 7 shows observed IMs (PGA, SA_0.2, SA_1.0, and SA_3.0), median IMs of six GMMs (Jea02, JB03, Eea15, JH21, BSSA14, and CB14) adjusted by c_k (Table 3), and median IMs of two combined GMMs minimizing σ_tt and σ_w using optimized weights (Table 5). For the comparison, we made two groups: 1) observed IMs with M_L 3.5–4.5 and median prediction of GMMs using M_L = 3.6, that is, the median M_L of the observed IMs and 2) observed IMs with M_L 4.5–5.8 and median prediction of GMMs using M_L = 5.1, that is, the median M_L of the observed IMs. Binned means of observed IMs are also shown in Figure 7. The combined GMM minimizing σ_tt (denoted Comb_tt) fits well with the observations for all IMs considered. Because $\delta B_e$ are not corrected in the observed IMs in Figure 7, the combined GMM minimizing σ_w (denoted Comb_w) does not fit well with the observations as good as Comb_tt; but still, Comb_w reasonably follows the data trend. This proves that the combined model with appropriate weights provides more accurate predictions across a wide range of IMs than prediction from a single model.

FIGURE 7

FIGURE 7. Comparison of four modified Korean GMMs and two global GMMs (Jea02, JB03, Eea15, JH21, BSSA14, and CB14) with combined GMMs minimizing σ_tt and σ_w (Comb_tt and Comb_w, respectively) for PGA, SA_0.2, SA_1.0, and SA_2.0 with M_L = 3.6 and 5.1 varying R_epi. Records and binned means for M_L 3.5–4.5 and M_L 4.5–5.8 are overlapped.

5.3 GMMs for large magnitude

The suggested GMMs (Comb_tt and Comb_w) mentioned earlier were tested with ground motion records with M_L ≤ 5.8. In South Korea, however, the magnitude range of 6.0–7.0 is used for seismic design practice (Earthquake Engineering Society of Korea, 2018), for which local records are not available. To check the effectiveness of the suggested GMM for this extrapolated range of magnitude, this section compares the trend of Comb_tt and Comb_w with the NGA West2 GMMs, which were developed using large magnitude ground motion records—the minimum of the upper limit M_W is 7.5.

Figure 8 shows selected IMs including PGA, SA_0.1, SA_0.2, SA_0.5, SA_1.0, and SA_2.0 with an increase of M_W from 5 to 7.5. For this calculation, R_epi and R_jb are set to 20 km, Z_hyp is set to 10 km, and other parameters for NGA West2 GMMs are set as shown in Table 2. The Comb_tt and Comb_w are extrapolated up to M_W 7.5 and compared with the NGAW2_eq (a combined GMM with equal weights for the four NGA West2 models with mean bias correction using coefficients from Table 3). At low periods (PGA and SA ≤ 0.5 s), the IMs of Comb_tt increase linearly, so they become excessively high at large magnitudes. This is attributed to the fact that Eea15 has higher weights at low periods (Table 5), and the magnitude application range of Eea15 is less than 5. Hence, including Eea15 in the Comb_tt is not appropriate if Comb_tt is used for M_W > 5. For Comb_w, on the other hand, the trend of magnitude saturation is comparable with the NGAW2_eq. The best four GMMs selected for Comb_w do not include Eea15 (Table 4), and the M_W range is up to 7 (Table 1). It means that Comb_w could be used for large M_W events (>5) even though it is not validated with observed IMs.

FIGURE 8

FIGURE 8. Comparison of magnitude scaling of PGA, SA_0.1, SA_0.2, SA_0.5, SA_1.0, and SA_2.0 for combined GMMs (Comb_tt, Comb_tt-woE, and Comb_w) with equally weighted NGA W2 models (NGAW2.eq).

For the suggestion of Comb_tt for large M_W, we omitted Eea15 and found the best four GMMs and their weights as indicated in Table 6. The trend of magnitude saturation of this model (Comb_tt-woE) is shown in Figure 8, which is comparable to NGAW2_eq for all spectral periods. However, the difference in average standard deviation of the Comb_tt-woE is 1.62% compared to all nine models (Figure 4), which is higher than the difference of Comb_tt (0.15%). Hence, it is recommended to use Comb_tt for small M_W, which is the optimized model, and use Comb_tt-woE for large M_W, which avoids excessive extrapolation at large M_W.

TABLE 6

TABLE 6. Suggested GMMs and weights minimizing σ_tt per IM for large M_W (>5).

6 Conclusion

This study suggests the best combination of GMMs for within-rock layer ground motion in South Korea. Local and global GMMs were used for model combination, and ground motions recorded in the within-rock layer were used for validation. The optimized weights were found using the QP technique, which provides fast and solid results in the model combination problem and effectively selects GMMs for the best combination, which minimizes the uncertainty of the combined model. Two cases of optimized weights were suggested: the first case for the prediction without ground motion observations (e.g., a future or scenario event) and the second case for the prediction with observations (e.g., a past event with ground motion records). Among the five local and four global GMMs considered, a combination of JB03, Eea15, JH21, and BSSA14 is suggested for the first case with small M_W (M_W < 5), and a combination of JB03, JH21, BSSA14, and CB14 is recommended for large M_W (M_W ≥ 5). The combination of Jea02, JB03, JH21, and CB14 is evaluated as the best for the second case. The performance was evaluated by comparing the standard deviation of the model residuals. The combined model reduces the average standard deviation of residuals by up to 10% compared to the best single GMM. The combination of even two models reduces the standard deviation significantly. The combined GMMs were also compared with equally weighted global GMMs for larger magnitudes where ground motion records for South Korea are not available. This comparison revealed that the suggested combined GMMs are applicable up to M_W = 7.5. The linear combination using the optimized weights suggested in this study minimizes the standard deviation of the residuals of the combined model.

The limitation of this study is that the suggested method is based on the empirical data, so records are required for the estimation of optimized weights and for the removal of the mean misfit. Also, the model uncertainty would be increased for the ground motions, for which the magnitude and distance were not covered in this study. In addition, the suggested model is for the within-rock condition, so coupling with site amplification is required to estimate ground motion amplitudes at the surface. The suggested model provides PGA and spectral accelerations for periods from 0.01 to 5 s given earthquake magnitude and hypocentral distance, but a recent study reveals that the near-fault effect such as a pulse-like seismic wave can result in more serious damage for slopes for the same magnitude and distance (Zhu et al., 2020; Zhu et al., 2022). The future studies will find optimized weights without records, estimate model uncertainty for ground motions, which are out of the range of the dataset, couple with a site amplification model and validate with surface ground motion records, and find the GMMs that consider near-fault effects, as promised.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

DK, J-KA, and T-WK contributed to the conception and design of the study. DJ organized and analyzed the database and prepared figures. DJ and DK wrote the first draft of the manuscript. JK-A and T-WK reviewed and analyzed the results and recommended corrections. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This work was supported by the Hanyang University (grant number HY-2019-N) and the Korea Meteorological Administration (grant number KMA2022-02121).

Acknowledgments

Authors greatly appreciate Hanyang University and KMA's supports. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect those of the Hanyang University and KMA.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Abrahamson, N. A., Silva, W. J., and Kamai, R. (2014). Summary of the ASK14 ground motion relation for active crustal regions. Earthq. Spectra 30, 1025–1055. doi:10.1193/070913eqs198m

CrossRef Full Text | Google Scholar

Abrahamson, N. A., and Youngs, R. R. (1992). A stable algorithm for regression analyses using the random effects model. Bull. Seismol. Soc. Am. 82, 505–510. doi:10.1785/bssa0820010505

CrossRef Full Text | Google Scholar

Akkar, S., and Kale, Ö. (2014). Reply to “comment on ‘A new procedure for selecting and ranking ground-motion prediction equations (GMPEs): The Euclidean distance-based ranking (EDR) method’ by Özkan Kale and Sinan Akkar” by Sum Mak, Robert Alan Clements, and Danijel Schorlemmer. Bull. Seismol. Soc. Am. 104, 3141–3144. doi:10.1785/0120140169

CrossRef Full Text | Google Scholar

Al Atik, L., Abrahamson, N., Bommer, J. J., Scherbaum, F., Cotton, F., and Kuehn, N. (2010). The variability of ground-motion prediction models and its components. Seismol. Res. Lett. 81, 794–801. doi:10.1785/gssrl.81.5.794

CrossRef Full Text | Google Scholar

Boore, D. M., Stewart, J. P., Seyhan, E., and Atkinson, G. M. (2014). NGA-West2 equations for predicting PGA, PGV, and 5% damped PSA for shallow crustal earthquakes. Earthq. Spectra 30, 1057–1085. doi:10.1193/070113eqs184m

CrossRef Full Text | Google Scholar

Boore, D. M. (1983). Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra. Bull. Seismol. Soc. Am. 73, 1865–1894. doi:10.1785/BSSA07306A1865

CrossRef Full Text | Google Scholar

Bozorgnia, Y., Abrahamson, N. A., Al Atik, L., Ancheta, T. D., Atkinson, G. M., Baker, J. W., et al. (2014). NGA-West2 research project. Earthq. Spectra 30, 973–987. doi:10.1193/072113eqs209m

CrossRef Full Text | Google Scholar

Campbell, K. W., and Bozorgnia, Y. (2014). NGA-West2 ground motion model for the average horizontal components of PGA, PGV, and 5% damped linear acceleration response spectra. Earthq. Spectra 30, 1087–1115. doi:10.1193/062913eqs175m

CrossRef Full Text | Google Scholar

Chiou, B. S.-J., and Youngs, R. R. (2014). Update of the Chiou and Youngs NGA model for the average horizontal component of peak ground motion and response spectra. Earthq. Spectra 30, 1117–1153. doi:10.1193/072813eqs219m

CrossRef Full Text | Google Scholar

Earthquake Engineering Society of Korea (EESK) (2018). Commentary of general seismic design code. KDS 17 (10 00), 2018.

Google Scholar

Emolo, A., Sharma, N., Festa, G., Zollo, A., Convertito, V., Park, J. H., et al. (2015). Ground-motion prediction equations for South Korea Peninsula. Bull. Seismol. Soc. Am. 105, 2625–2640. doi:10.1785/0120140296

CrossRef Full Text | Google Scholar

Goldfarb, D., and Idnani, A. (1983). A numerically stable dual method for solving strictly convex quadratic programs. Math. Program. 27, 1–33. doi:10.1007/bf02591962

CrossRef Full Text | Google Scholar

Higham, N. (2002). Computing the nearest correlation matrix - a problem from finance. IMA J. Numer. Anal. 22, 329–343. doi:10.1093/imanum/22.3.329

CrossRef Full Text | Google Scholar

Hong, T. K. (2017). Establishment of standardized procedure and development of national seismic hazard maps Report R&D 2017-MPSS31-107. Sejong, Korea: Ministry of Interior and Safety, 635. (in Korean).

Google Scholar

Jee, H. W., and Han, S. W. (2021). Regional ground motion prediction equation developed for the Korean Peninsula using recorded and simulated ground motions. J. Earthq. Eng. 26, 5384–5406. doi:10.1080/13632469.2021.1871682

CrossRef Full Text | Google Scholar

Jo, N., and Baag, C. (2003). Estimation of spectrum decay parameter k and stochastic prediction of strong ground motions in southeastern Korea. J. Earthq. Eng. Soc. Korea 7, 59–70. (in Korean). doi:10.5000/EESK.2003.7.6.059

CrossRef Full Text | Google Scholar

Jo, N., and Baag, C. (2001). Stochastic prediction of strong ground motions in southern Korea. J. Earthq. Eng. Soc. Korea 5, 17–26. (in Korean).

Google Scholar

Junn, J.-G., Jo, N.-D., and Baag, C.-E. (2002). Stochastic prediction of ground motions in southern Korea. Geosci. J. 6, 203–214. doi:10.1007/bf02912691

CrossRef Full Text | Google Scholar

Kale, Ö., and Akkar, S. (2013). A new procedure for selecting and ranking ground-motion prediction equations (GMPEs): the Euclidean distance-based ranking (EDR) method. Bull. Seismol. Soc. Am. 103, 1069–1084. doi:10.1785/0120120134

CrossRef Full Text | Google Scholar

Korea Meteorological Administration (KMA) (2022). National earthquake comprehensive information system (NECIS). Available at: https://necis.kma.go.kr/ (Accessed September 22, 2022).

Google Scholar

Kramer, S. L. (1996). Geotechnical earthquake engineering. Upper Saddle River, NJ: Prentice-Hall.

Google Scholar

Kwak, D. Y., Seyhan, E., and Kishida, T. (2018). “A method of linear combination of multiple models for epistemic uncertainty minimization,” in 11th National Conference on Earthquake Engineering, Los Angeles, CA, 25–29 June, 2018.

Google Scholar

Nizamani, Z. A., and Park, D. (2021). Testing ground-motion prediction equations against moderate magnitude earthquake data recorded in Korea. Bull. Seismol. Soc. Am. 111, 321–338. doi:10.1785/0120200095

CrossRef Full Text | Google Scholar

Noh, M. H., and Lee, K. H. (1995). Estimation of peak ground motions in the Southeastern part of the Korean peninsula (II): Development of predictive equations. J. Geol. Soc. Korea 31, 175–187. (in Korean).

Google Scholar

Park, D. H., Lee, J. M., Baag, C.-E., and Kim, J. K. (2001). Stochastic prediction of strong ground motions and attenuation equations in the southeastern Korean Peninsula. J. Geol. Soc. Korea 37, 21–30. (in Korean).

Google Scholar

R Core Team (2022). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing.

Google Scholar

Scherbaum, F., Delavaud, E., and Riggelsen, C. (2009). Model selection in seismic hazard analysis: An information-theoretic perspective. Bull. Seismol. Soc. Am. 99, 3234–3247. doi:10.1785/0120080347

CrossRef Full Text | Google Scholar

Sheen, D.-H., Kang, T.-S., and Rhie, J. (2018). A local magnitude scale for South Korea. Bull. Seismol. Soc. Am. 108, 2748–2755. doi:10.1785/0120180112

CrossRef Full Text | Google Scholar

Yun, K.-H., Park, D.-H., Chang, C.-J., and Sim, T.-M. (2008). “Estimation of aleatory uncertainty of the ground-motion attenuation relation based on the observed data,” in Proceedings of EESK conference, Seoul, Korea, 21 March, 2008, 116–123.

Google Scholar

Zhu, C., Cheng, H., Bao, Y., Chen, Z., and Huang, Y. (2022). Shaking table tests on the seismic response of slopes to near-fault ground motion. Geomech. Eng. 29 (2), 133–143.

Google Scholar

Zhu, C., Huang, Y., and Sun, J. (2020). Solid-like and liquid-like granular flows on inclined surfaces under vibration – implications for earthquake-induced landslides. Comput. Geotech. 123, 103598. doi:10.1016/j.compgeo.2020.103598

CrossRef Full Text | Google Scholar