GS24b and GS24bc Ground Motion Models for Active Crustal Regions Based on Non-Traditional Modeling Approach

1. Introduction

Ground motion models (GMMs) also called ground motion prediction equations (GMPEs), or originally attenuation relations, use datasets of recorded ground motion parameters at multiple seismic stations during different earthquakes and in various seismic source regions in a similar tectonic environment to generate equations. These equations are later used to estimate site-specific ground motions that may shake a site if an earthquake of a certain magnitude occurs at a nearby location. These models describe the distribution of ground motion in terms of a median and a logarithmic standard deviation and are crucial in assessing seismic hazard, thereby providing estimates of the loading that a structure may undergo during a future earthquake.

Typically, Ground Motion Models are developed using an empirical regression of observed amplitudes against an available set of predictor variables. Joyner and Boore [1,2] proposed performing analyses of data using the two-stage regression based on the algorithm [3,4] for applying the random effects model to regression analyses. Abrahamson and Youngs [5] presented an alternative mixed-effects algorithm that is more stable according to the authors although computationally less efficient. In these approaches, coefficients were obtained separately for each period resulting in response spectra that demonstrates jaggedness and consequently requiring smoothing (e.g., [6,7,8,9]).

In contrast to the above-described traditional approaches, we are using the non-traditional approach first introduced by Graizer and Kalkan [10,11,12] and later expanded in [13,14]. In this method at the first stage the closed form expression of the two functions (peak ground acceleration$PGA$ and spectral shape $SA$) composition is developed separately for the active crustal and stable continental regions (ACR and SCR) approximating ground motion attenuation. This type of modeling was originally based on the expanded NGA-West1 dataset and is developed using the Nelder-Mead method of nonlinear minimization [15]. However, this closed form approximation is not flexible enough to describe all variations of response spectral attenuation and to produce lower standard error. This is why we incorporated the second stage when previously developed closed form generic model is adjusted based on residuals for each period.

In [16] the new set of strong motion data from the two strongest 2023 earthquakes in Turkey were used to test the ergodic GK17 GMM [14]. The GK17 model developed using the NGA-West2 database [17] for the active crustal regions was applied to the dataset of recordings from the two moment magnitude $\left(\mathbf{M}\right)$ 7.8 and 7.5 earthquakes in Turkey [18,19]. The GK17 model demonstrates acceptable performance while mostly underpredicting spectral accelerations at near fault up to ~100 km and far-field more than ~400 km rupture distances for short periods $T<1$ s. In [16] GK17 model was modified by applying additional $R_{rup}$ distance and $V_{S30}$ residuals corrections creating an updated GK model tuned for Turkey and called the GK_T model. Turkey specific GK_T partially non-ergodic model shows better agreement with recorded data than the ergodic GK17 model, especially at short periods and most importantly short rupture distances.

We recently developed two new GMMs [20]:

• The global backbone GS24b model that uses the closed form approximation of the spectral acceleration as a multiplication of the $PGA$ and spectral shape (normalized spectral acceleration spectrum) functions. This model can be later used for adjusting to the specific ACR region (e.g., Southern and Northern California) creating partially non-ergodic models.
• The ergodic GS24 model representing the backbone GS24b GMM adjusted for the depth to the shear-wave velocity horizon of 2.5 km/s ($Z_{2.5}$), style of faulting ($SoF$) and also for the moment magnitude $\mathbf{M}$, time-averaged shear-wave velocity in the upper 30 meters $V_{S30}$ and closest distance to the fault rupture $R_{rup}$ residuals.

In this paper we are presenting the backbone model GS24b and the GS24bc model created by adjusting backbone model for residuals. We also tested previously developed GK17 GMM [14] against the NGA-West2 dataset expanded with recordings from the three 2023 strong Turkish earthquakes.

2. Datasets

The new GS24b model is based on Pacific Earthquake Engineering Research (PEER) Center Next Generation Attenuation Phase 2 (NGA-West2) data expanded with recordings from the three 2023 Turkish earthquakes with moment magnitudes of $\mathbf{M}$ 6.3, 7.5 and 7.8 [18,19] ground motion database compiled from shallow crustal earthquakes in active crustal regions to develop a GMM for the “average” RotD50 (50th percentile or median response [21] horizontal components of peak ground acceleration ($PGA$), peak ground velocity ($PGV$) and 5% damped elastic pseudo-absolute acceleration response spectral ordinates ($PSA$) at 21 oscillator periods (T) ranging from 0.01 to 10 sec [17]. The specific spectral periods are 0.01, 0.02, 0.03, 0.05, 0.075, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.75, 1, 1.5, 2, 3, 4, 5, 7.5 and 10 sec. This set of periods is considered to be representative in current US engineering practice. The number of predictors used in the current model is limited to a few measurable parameters: moment magnitude ($\mathbf{M}$), closest distance to the fault rupture plane ($R_{rup}),$ time-averaged shear-wave velocity in the upper 30 m of the geological profile ($V_{S30}$), and apparent (intrinsic and scattering) attenuation factor $\left(Q_{SA}\right(f,\mathbf{M})$) of 5% damped spectral acceleration.

The currently used dataset consists of the original subset of the NGA-West2 dataset with 13,241 recordings with the addition of 685 Turkish recordings with the total number of 401 earthquakes in the combined dataset. Figure 1 demonstrates the distribution of chosen recordings with respect to moment magnitude $\mathbf{M}$ and closest distance to the rupture $R_{rup}$ (upper panel), and with respect to the $PGA$ and $V_{S30}$ (middle panel). NGA-West2 data are shown with open circles and additional Turkish data with red circles in Figure 1. We limited the dataset by potentially damaging earthquakes with 4.0≤ $\mathbf{M}$ ≤ 7.9 and rupture distances $R_{rup }$≤ 400 km. The GS24 models are based on recordings from California, Alaska (crustal events), Taiwan, Turkey, Italy, Greece, New Zealand and Northwest China ACRs characterized by similar tectonic environments. We did not include data from Japan because most of the sites in Japan are characterized by a subsurface geology significantly different from the site conditions in other ACRs. Similar to our previous GMMs GS24b should be considered a global model since it includes recordings from multiple ACRs.

Aftershock records are treated differently in recent studies than data from main shocks because of concern that the median ground motions from aftershocks are systematically lower than those from the mainshocks. Existing publications provide conflicting results, in some cases finding different magnitude scaling for aftershocks relative to main shocks (e.g., [7]), while in other cases finding similar motions for similar magnitudes [22]. As a result, some GMM developers exclude aftershock data from regressions [7], exclude aftershocks closest to the mainshock [8,9], include a separate coefficient for the attenuation term for the aftershocks data used [6], or simply include all well recorded aftershocks [14,23]. Following our previous approach, we are including all well recorded aftershocks in the datasets.

The record processing procedures applied to all records in the NGA-West2 flatfile include the selection of record-specific corner frequencies to optimize the usable frequency range. The most important filter applied to the data is the low-cut filter, which removes low frequency noise. For each record, the maximum usable period applied in our analysis was taken as the inverse of the lowest usable frequency given in the NGA-West2 flatfile. The lowest usable frequency is usually equals 1.25 of the high-pass (equivalent to low-cut) corner frequencies used in the processing of the two horizontal components. Figure 1 (lower panel) demonstrates the number of data points used in our models for different periods with decreasing amount of data for longer periods.

The current versions of the GS24b model includes testing that uses the following three datasets:

• The 1st dataset of 4 ≤ $\mathbf{M}$ ≤ 7.9 and distances $R_{rup }$ ≤ 400 km called M4_R400 includes all 13,926 data points (Figure 1, lower panel, and Figure 2).
• The 2^nd dataset (subset of the 1^st dataset) includes 5,063 data points for $\mathbf{M}$ ≥ 5.0 and $R_{rup }$ ≤ 150 km (called M5_150), and the final dataset:
• The 3rd dataset (also subset of the 1^st dataset) includes 6,045 data points covering the range of 4≤ $\mathbf{M}$ ≤7.9 and $R_{rup }$≤ 250 km (called M4_R250).

The 1^st dataset consists of the same NGA-West2 dataset [17] as used in the development of GK17 [14] model (total of 13241 data points) enhanced by the data from the three recent 2023 Turkish earthquakes with $\mathbf{M}$ 6.3, 7.5 and 7.8 [18,19]. The 2^nd dataset is a subset of the 1^st dataset limited by magnitude and rupture distance. The 3^rd dataset is also a subset of the 1^st group. This dataset was created based on the recordings of earthquakes that can potentially produce structural damage. It consists of the seven half-magnitude bins with varying distance thresholds depending upon magnitude correlated to potentially damaging rupture distance:

• 4.0 ≤ $\mathbf{M}$ < 4.5 and $R_{rup }$≤ 25 km
• 4.5 ≤ $\mathbf{M}$ < 5.0 and $R_{rup }$ ≤ 50 km
• 5.0 ≤ $\mathbf{M}$ < 5.5 and $R_{rup }$ ≤ 75 km
• 5.5 ≤ $\mathbf{M}$ < 6.0 and $R_{rup }$ ≤ 100 km
• 6.0 ≤ $\mathbf{M}$ < 7.0 and $R_{rup }$ ≤ 150 km
• 7.0 ≤ $\mathbf{M}$ < 7.5 and $R_{rup }$ ≤ 200 km
• 7.5 ≤ $\mathbf{M}$ ≤ 7.9 and $R_{rup }$ ≤ 250 km

We consider the 3^rd dataset to be the most important from the engineering application point of view and will call it “engineering” dataset.

As shown in Figure 2, the first dataset is heavily dominated by the relatively low magnitude (4 ≤ $\mathbf{M}$ ≤ 5) earthquakes less likely to produce damage, while the effect of lower magnitude earthquakes is absent in the second dataset. The 3^rd dataset is the most magnitude-distance balanced and consequently most important from the engineering and seismic hazard assessment point of view, practically limiting data to the events that can potentially produce any damaging effects. This dataset is also much more uniform in terms of magnitude distribution and mostly limited to more uniform set of accelerographs recordings (Figure 2).

3. GS24b Backbone Ground Motion Model

The backbone model development included multiple steps, like those described in [11,12]. As a first step exploratory analyses using the first largest dataset were performed and demonstrated the demand to update some of the previously developed relations coefficients. However, all previously described formulations remain valid.

3.1. PGA Scaling

We started with testing and adjusting the backbone model presented in Graizer (2018) and originally developed by [10,11,12] to the described datasets. As a reminder, spectral acceleration is a closed form combination of the two functions $PGA$ and normalized spectral shape${ SA}_{norm}$ functions:

$$SA=PGA x{ SA}_{norm}$$

(1)

$$\log \left(PGA\right)=F_1\left(M\right)+F_2\left(R\right)$$

(2)

Where $F_1\left(M\right)$ is PGA magnitude scaling and $F_2\left(R\right)$ is distance scaling. $PGA$ magnitude scaling function has a linear scaling in logarithmic space for small magnitudes and the same style of saturation approximation function as in our previous model for larger magnitudes [24]:

$$F_1\left(\mathbf{M},k_{scale}\right)=\left\{\begin{array}{c}log\left(\left[c_{21}\mathit{exp}(c_{22}\mathbf{M}\right)]k_{scale}\right) 4\le M<5.0\\ \\ log\left(\left[c_1\mathit{arctan}(\mathbf{M}+c_2\right)+c_3]k_{scale}\right) M\ge 5.0\end{array}\right.$$

(3)

where $c_n$are coefficients and $k_{scale}$ is a scaling factor (Figure 3). In the current GS24b model we modified the turning point in magnitude scaling from M = 5.5 in GK17 to 5.0 effectively increasing the scaling for $\mathbf{M}$ < 5.5.

In the GS24b backbone model, like in [13,14] we are assuming attenuation of ground motion in the near-source region to be associated with the shear-waves geometric spreading of R_rup^{- 1} while at distances larger than about 50 km maximum ground motion to be associated with surface waves attenuating of R_rup^-0.5. This average transition distance was found by fitting NGA-West2 data [14].

$$F_2(R_{rup})=\left\{\begin{array}{c}log{\displaystyle \frac{1}{\sqrt{{\left[1-\frac{R_{rup}}{R_2}\right]}^2+4D_2^2(R_{rup}/R_2)}}}{,R}_{rup}<50km\Rightarrow ~{\displaystyle \frac{1}{R_{rup}}}\\ log{\displaystyle \frac{w}{\sqrt{[1-\sqrt{\frac{R_{rup}}{R_2}}{]}^2+4D_2^2\sqrt{\frac{R_{rup}}{R_2}}}}},R_{rup}\ge 50km\Rightarrow ~{\displaystyle \frac{1}{\sqrt{R_{rup}}}}\end{array}\right.$$

(4)

$$w=log\sqrt{{\displaystyle \frac{(1-\sqrt{\frac{50}{R_2}}{)}^2+4D_2^2\sqrt{\frac{50}{R_2}}}{{\left[1-\frac{50}{R_2}\right]}^2+4D_2^2(50/R_2)}}}$$

Parameter $w$ is a scaling factor applied to Equation 4 to avoid step in PGA attenuation slope, and R₂ is the corner distance in the near-source defining the plateau without significant PGA attenuation. Parameter $R_2$ is directly proportional to the moment magnitude of an earthquake; the larger is $\mathbf{M}$, the wider is the plateau defined by R₂. The scaling law of the corner distance was originally developed for the active tectonic environment [10], and also used in previous models [14,25]:

$$R_2=c_4\mathbf{M}+c_5$$

(5)

Equations (4 and 5) imply that for larger magnitudes, the turning point on the attenuation curve occurs at larger distances varying from 1.4 km for $\mathbf{M}$ = 4.0 to 10.4 km for $\mathbf{M}$ = 8.0 [10]. We assigned parameter $D_2$= 0.5, which is equivalent to minor “bump”: Increase in amplitude of ground motion at certain distances from the fault assuming over saturation of PGA and a smooth transition from a plateau to the R^-1 or R^-0.5 geometrical spreading.

3.2. Spectral Shape Model

Response spectral shape is modeled using closed form approximation function developed in [11]. This function is a combination of a single-degree-of-freedom oscillator and a modified log-normal probability density function. The updated spectral shape model ($S{A}_{norm,0}$) is a continuous function and is formulated as shown in Equation 6:

$$\mathit{S}{\mathit{A}}_{\mathit{n}\mathit{o}\mathit{r}\mathit{m},\mathbf{0}}\left(\mathit{T},\mathit{M},{\mathit{R}}_{\mathit{r}\mathit{u}\mathit{p}},{\mathit{V}}_{\mathit{S}\mathbf{30}}\right)=$$

$$\mathit{I}(\mathit{M},{\mathit{R}}_{\mathit{r}\mathit{u}\mathit{p}})\mathit{exp}\left(-\mathbf{0.5}{\left({\displaystyle \frac{\mathit{ln}(\mathit{T})+\mathit{\mu }(\mathbf{M},{\mathit{R}}_{\mathit{r}\mathit{u}\mathit{p}},{\mathit{V}}_{\mathit{S}\mathbf{30}}}{\mathit{S}(\mathit{M},{\mathit{R}}_{\mathit{r}\mathit{u}\mathit{p}})}}\right)}^{\mathbf{2}}\right)+{\left[{\left(\mathbf{1}-{\left({\displaystyle \frac{\mathit{T}}{{\mathit{T}}_{\mathit{s}\mathit{p},\mathbf{0}}}}\right)}^{\mathit{\zeta }}\right)}^{\mathbf{2}}+\mathbf{4}{{\mathit{D}}_{\mathit{s}\mathit{p}}}^{\mathbf{2}}{\left({\displaystyle \frac{\mathit{T}}{{\mathit{T}}_{\mathit{s}\mathit{p},\mathbf{0}}}}\right)}^{\mathit{\zeta }}\right]}^{-{\displaystyle \frac{\mathbf{1}}{\mathbf{2}}}}$$

$$\mathit{\mu }(\mathit{M},{\mathit{R}}_{\mathit{r}\mathit{u}\mathit{p}},{\mathit{V}}_{\mathit{S}\mathbf{30}})={\mathit{m}}_{\mathbf{1}}{\mathit{R}}_{\mathit{r}\mathit{u}\mathit{p}}+{\mathit{m}}_{\mathbf{2}}\mathit{M}+{\mathit{m}}_{\mathbf{3}}$$

$$\mathit{I}(\mathit{M},{\mathit{R}}_{\mathit{r}\mathit{u}\mathit{p}})=({\mathit{a}}_{\mathbf{1}}\mathit{M}+{\mathit{a}}_{\mathbf{2}})\mathit{exp}({\mathit{a}}_{\mathbf{3}}{\mathit{R}}_{\mathit{r}\mathit{u}\mathit{p}})$$

$$\mathit{S}(\mathit{M},{\mathit{R}}_{\mathit{r}\mathit{u}\mathit{p}})={\mathit{s}}_{\mathbf{1}}{\mathit{R}}_{\mathit{r}\mathit{u}\mathit{p}}-({\mathit{s}}_{\mathbf{2}}\mathit{M}+{\mathit{s}}_{\mathbf{3}})$$

$${\mathit{T}}_{\mathit{s}\mathit{p},\mathbf{0}}(\mathit{M},{\mathit{R}}_{\mathit{r}\mathit{u}\mathit{p}},{\mathit{V}}_{\mathit{S}\mathbf{30}})=\mathit{max}\left\{\mathbf{0.3},\mathit{a}\mathit{b}\mathit{s}({\mathit{t}}_{\mathbf{1}}{\mathit{R}}_{\mathit{r}\mathit{u}\mathit{p}}+{\mathit{t}}_{\mathbf{2}}\mathit{M}+{\mathit{t}}_{\mathbf{3}})\right\}$$

(6)

m1	m2	m3	a1	a2	a3	Dsp
-0.0012	-0.38	*	0.017	1.27	0.0001	0.75

$$\mathbf{*}{\mathit{m}}_{\mathbf{3}}=\mathbf{0.0006}{\mathit{V}}_{\mathit{S}\mathbf{30}}+\mathbf{4.15}$$

t1	t2	t3	s1	s2	s3
0.001	0.59	**	0.00	0.077	0.3251	***

$$\mathit{*}\mathit{*}{\mathit{t}}_{\mathbf{3}}=-\mathbf{0.0005}{\mathit{V}}_{\mathit{S}\mathbf{30}}-\mathbf{2.9}$$

$${\mathbf{*}\mathbf{*}\mathbf{*}(\mathrm{M},\mathit{R}}_{\mathit{r}\mathit{u}\mathit{p}}) =-\mathbf{0.25} \mathrm{M}+\mathbf{3} - \mathbf{0.0004} {\mathbf{R}}_{\mathbf{r}\mathbf{u}\mathbf{p}}$$

The $S{A}_{norm,0}(T,\mathbf{M},R_{rup},V_{S30}$) is a continuous function of spectral period (or frequency), moment magnitude, closest distance to the fault rupture surface and shear-wave velocity in the upper 30 m of profile $V_{S30}$. Spectral shape function is anchored at $PGA$ = 1g (Figure 4a). Spectral maximum shifts toward longer periods with an increase of magnitude (Figure 4b) and distance (Figure 4c) and also shifts toward shorter periods with an increase of shear-wave velocity (Figure 4d). Spectral slope at longer periods decreases with an increase of magnitude and rupture distance (Figure 4b).

4. GS24b Backbone Model

We developed a backbone model which is a combination of the spectral shape (Equation 6), $PGA$ (Equation 2), site amplification and apparent anelastic attenuation functions:

$$S{A}_{backbone}=G_1\left(V_{S30},f\right)\times G_2\left(Q_{SA}\left(f,M\right)\right)\times S{A}_{norm,0}\times PGA$$

(7)

where $f$ is frequency. In Equation 7, $S{A}_{backbone}$ is response spectral acceleration, $G_1(V_{s30},f)$ is for shallow site amplification and $G_2\left(Q_{SA}\left(f,M\right)\right)$ is for apparent anelastic (intrinsic and scattering) attenuation correction [26].

4.1. Site Response Term

The $G_1(V_{S30},f)$ term (referred to these items as filters in our previous publications) in Equation 7 is for shallow site amplification due to geological conditions in the upper section of geological profile as characterized by the parameter${ V}_{S30}$. Site correction was developed in [25] based on multiple runs of different representative shear-wave velocity profiles through SHAKE-type [27] 1-D equivalent-linear (EQL) programs using time histories and random vibration theory (RVT) approaches [28] and on EQL RVT type code developed by the staff of the U.S. Nuclear Regulatory Commission. Soil profiles were chosen from the set of California profiles collected by the U.S. Geological Survey, California Geological Survey, UC Santa Barbara and other organizations (e.g., [29,30]) and the NRC library of appropriate profiles. For the $G_1(V_{S30},f)$ filter we used same functional form as that used in [14] for ACR:

$$G_1(V_{S30},f)=1+{\displaystyle \frac{k_{Vs30}}{\sqrt{{\left[1-\frac{f_{Vs30}}{f}\right]}^2+1.96\left(\frac{f_{Vs30}}{f}\right)}}}$$

(8)

$$k_{Vs30}=-0.5\mathit{ln}(V_{S30}/1100)180\le V_{S30}\le 1100m/s$$

$$f_{Vs30}=\left\{\begin{array}{c}{\displaystyle \frac{275}{120}}-2.0 180\le V_{S30}\le 275 m/s\\ {\displaystyle \frac{V_{S30}}{120}}-2.0 275\le V_{S30}\le 1100 m/s\end{array}\right.$$

Like in [14] we characterized our reference rock (also sometimes called western US hard rock) as $V_{S30}$ = 1100 m/s. This value of reference rock is identical to that of [8] and close to $V_{S30}$ = 1130 m/s in [9] and to $V_{S30}$= 1180 m/s in [6]. Site amplification functions are calculated for different $V_{S30}$ relative to the reference rock of $V_{S30}$ =1100 m/s. The limits of model saturation are 180 ≤ $V_{S30}$ ≤ 1100 m/s covering most of the velocity profiles in the database. Figure 5 shows $V_{S30}$ site amplification from Equation 8 representing 1-D site amplification derived from a collection of S-wave velocity profiles basically demonstrating a semi-empirical approach.

Our approach to implementing nonlinearity is discussed in detail in Graizer (2017 and 2018). Before implementing nonlinearity in our model, we reviewed current approaches for modeling it in site response (e.g. [8,31,32] and performed independent analyses. Our assessment showed significant variability in results and approaches. We can speculate that there is no simple nonlinear correlation between site amplification and $V_{S30}$. For example, [33,34] demonstrated that the addition of kappa is needed to get reasonable soft rock to hard rock spectral amplification ratios. Based on the analysis of strong motion recordings from the 2014 $\mathbf{M}$ 6.0 South Napa earthquake at Carquinez Bridge geotechnical array [35] demonstrated that the apparent S-wave velocity decreased when the $PGA$ was greater than 0.07g, confirming that soil nonlinear behavior was observed due to strong shaking. Another argument about the deficiency in current approaches to nonlinearity is that they are actually capping on the possible level of strong ground motion, contradicting a number of strongest records: e.g., more than 2 g recorded at Parkfield Fault Zone 16 station during the 2004 $\mathbf{M}$ 6.0 Parkfield earthquake (Shakal et al., 2006), 1.9 g Tarzana record from the $\mathbf{M}$ 6.7 Northridge earthquake or a 2.2 g Heathcote Valley School record from the $\mathbf{M}$ 6.3 Christchurch, New Zealand earthquake [36]. As a result of the tests and analysis performed, we decided to implement nonlinearity in a simple way by putting cap on the site amplification. As can be seen from Equation 8 and Figure 5 the coefficient of amplification reaches its maximum at$V_{S30}$ = 180 m/s and is constant for lower velocities. Similarly, coefficient of amplification reaches its minimum at $V_{S30}$ = 1100 m/s and is constant and equal to 1 for higher S-wave velocities (Figure 5).

4.2. Apparent Attenuation of Spectral Accelerations

Geometric attenuation shown in Equation 4 is representing power function type attenuation (R^-α) while apparent anelastic attenuation is responsible for the ground motion exponential decay (e^(-kR)). Apparent attenuation represents combined intrinsic absorption and scattering dissipation responsible for the beyond geometrical spreading attenuation of spectral acceleration. We avoid calling it “anelastic” because it includes elastic scattering.

The $G_2(Q_{SA},f,\mathbf{M})$ filter adjusts the distance attenuation rate by including the apparent attenuation of the response spectra given as:

$$G_2(Q_{SA},f,\mathbf{M})=\mathit{exp}(-{\displaystyle \frac{\pi f{R}_{rup}}{\beta Q_{SA}\left(f,\mathbf{M}\right)}})$$

(9)

in which $Q_{SA}(f,\mathbf{M})$ is a frequency dependent apparent attenuation quality factor of SA amplitudes, and $\beta$ is apparent wave velocity. It is important distinguishing between the well-known “seismological” $Q\left(f\right)$ measured using Fourier spectra of S-, Lg- or coda-waves (e.g., [37,38,39,40]), and the quality factor of response spectral acceleration amplitudes $Q_{SA}(f,\mathbf{M})$. Estimates of apparent anelastic attenuation and the corresponding response spectra $Q_{SA}(f,\mathbf{M})$ quality factors were performed using inversions similar to those applied to the Fourier spectral amplitudes by for example in [41] but replacing Fourier with response spectral amplitudes. For the ACR inversions were performed at 21 frequencies for 7 different averaging distance intervals from 50 to 400 km (Figure 6). An average over frequency range of 0.1 to 100 Hz approximation of apparent quality factor based on our inversion results [14,26]:

$$Q_{SA}(f,\mathbf{M})=Q_{SA0}{f}^{0.96}$$

(10)

where $Q_{SA0}$ = 120 for the ACRs corresponding to the average $\mathbf{M}$ = 5.25 in the currently used NGA-West2 dataset of 4 ≤ $\mathbf{M}$ ≤ 7.9 and distances $R_{rup}$ ≤ 400 km. In [14] we used approximation of $Q_{SA}\left(f,\mathbf{M}\right)$ as a combination of three frequency intervals. To avoid jaggedness when combining these three intervals’ approximations we are approximating $Q_{SA}\left(f,\mathbf{M}\right)$ with a smooth Equation in the interval 50 ≤$R_{rup}$ ≤ 400 km:

$$Q_{SA}(f,\mathbf{M})={cs}_2+{cs}_4{Q}_1\mathbf{M}$$

(11)

$$Q_1=Q_{SA0}{f}^{\alpha }$$

$$\alpha ={cs}_1{f}^{{cs}_2}$$

Comparisons of the seismological and response spectra apparent attenuations demonstrate the differences in slope and amplitudes of these functions. In the meantime, the values of $Q_{SA0}\approx Q_0$ (seismological quality factor at a frequency of 1 Hz) apparent quality factors of response spectra and seismological $Q_S$ are similar at a frequency near 1 Hz. As was shown in [14], another important consequence of Equation 10 is that apparent attenuation is almost frequency independent at rupture distances of more than 50 km.

$$Q_{SA}\left(f,\mathbf{M}\right)\approx Q_{SA0}{f}^1 {\Rightarrow G}_2\left(Q_{SA},f\right)\approx \mathit{exp}\left(-{\displaystyle \frac{{\pi R}_{rup}}{ Q_{SA0}}}\right); R_{rup}>50 km$$

(12)

Figure 6 also demonstrates magnitude dependence of the resulting effective apparent attenuation factor. We found that magnitude dependent $Q_{SA}$ better fits existing data. A reason for this may be that larger earthquake faults generally penetrate deeper in the crust than the smaller ones, and consequently a significant portion of waves generated by the source are propagating through the deeper layers with higher $Q_{SA}$.

If $Q_{SA}\left(f,\mathbf{M}\right)\approx Q_{SA0}{f}^1$the anelastic attenuation rate $\gamma \left(f\right)$ becomes frequency independent as shown in [26], and the exponential factor in Equation 12 may be included in the distance-dependent power law R^-α. Frequency independent attenuation may result in an alternative to the exponential power-law attenuation and correspondently instead of R^-0.5 apparent geometrical attenuation rate can become of ~R^-09 at distance more than 50 km exceeding geometrical spreading rate of surface waves. This approach applied to Fourier spectra was previously suggested in [42].

These results are empirical and applicable to shallow crustal earthquakes in active crustal regions associated with the geometrical spreading of surface waves. Most importantly, apparent attenuation of response spectral amplitudes is different from that of the “seismological” $Q_S$ -factor and should be estimated based on actual spectral acceleration data and not transferred from seismological measurements. More details about apparent attenuation of spectral accelerations in different tectonic environments are presented in [26].

4.3. Examples of Backbone Model

Figure 7 shows examples of response spectral accelerations at rupture distances of 10, 50 and 200 km for moment magnitude earthquakes in the range of 4.5 ≤ $\mathbf{M}$ ≤ 8 calculated using the backbone model. As expected, all curves are smooth with maxima shifting toward longer periods and long period slope decreasing for larger magnitudes (Figure 7). Figure 8 shows examples of different periods attenuation with distance to the fault. As prescribed by Equation 4 attenuation rate changes from that of shear to surface waves at $R_{rup}$= 50 km.

The backbone model GS24b was tested against all three above-described subsets. As expected, standard error is higher for the first largest subset (M4_R400) (Figure 9). If compared to the GK17 model, GS24b mostly demonstrates higher sigma except for short period range of 0.01-0.1 s for the M4_R250 dataset considered the most important. In the meantime, we did not expect the backbone model to perform as well as the final GK17 model adjusted for residuals as well as for the style of faulting and additional deep sediment (basin) effect. Another important conclusion is that GK17 model demonstrates acceptable performance against all the three datasets confirming results achieved in Graizer (2024) using Turkish data.

The backbone GS24b model includes magnitude $F_1\left(\mathbf{M},k_{scale}\right)$, distance geometrical spreading$F_2\left(R_{rup}\right)$, site response term $G_1(V_{S30},f)$ and apparent (anelastic) attenuation $G_2(Q_{SA},f,\mathbf{M})$ scaling. We referred to them as filters in our previous publications. We intentionally did not incorporate style of fault $G_3(SoF,T)$ and an additional deep sediment $G_4(Z_{2.5},T)$ effect corrections into GS24b backbone model to allow future adjustments to the specific areas.

As expected GS24b is completely smooth (Figure 7 and Figure 8) representing closed form approximation. However, this model still results in residual biases and should be considered as a first-step approximation. Table 1 presents natural logarithmic standard deviation sigma for the GS24b model.

Peak Ground Velocity

The current GS24b backbone model includes PGV calculation not provided in the previous GK-17 model. To compute PGV we are converting 5% damped SA into spectral velocity and calculating maximum of 5% damped SV.

$$max\left(S{V}_{backbone}\right)=max(980 {\displaystyle \frac{S{A}_{\mathit{backbone}}}{2\mathit{\pi f}}})=max(980 {\displaystyle \frac{S{A}_{backbone} T}{2 \pi }})$$

(13)

assuming that spectral acceleration is presented in units of g, and PGV is in cm/s. In the backbone model we are calculating maximum amplitude of spectral velocity SV which may not be equal to the recorded PGV. However, the M_, Rrup and V_S30 residuals don’t show significant bias with an acceptable natural logarithmic standard deviation sigma ${\sigma }_{PGV}$ = 0.684.

4.4. GS24bc Model

The purpose of the next step is to demonstrate what can be done in a relatively simplified way by adjusting for residuals the generic active crustal region backbone GS24b model to the target region having a specific ground motion database. Similar approach was previously introduced in Graizer (2024) by creating a partially non-ergodic model for Turkey.

As compared to the GS24b backbone model, the GS24bc model (GS24b backbone model corrected for moment magnitude M, shear-wave velocity $V_{S30}$ and fault rupture distance Rrup residuals) is an example of how backbone model can be used to create a non-ergodic model specific for certain region or database. The residuals from the GS24b demonstrate bias, while Figure 10, Figure 11 and Figure 12 demonstrating residuals of the GS24bc model show no bias after the backbone model was corrected for residuals.

As another step we are correcting maximum spectral velocity SV for residuals to get PGV. Figure 13 demonstrates PGV residuals versus moment magnitude M, shear-wave velocity $V_{S30}$ and fault rupture distance Rrup . As shown in Figure 13 there is no bias after $max\left(S{V}_{backbone}\right)$) was corrected for residuals. Compared to [43] showing the average ratio between PGV and $max\left(SV\right)$~2.2 our average is ~ 1.7.

Figure 14 demonstrates examples of response spectral accelerations at rupture distances of 10, 50 and 200 km for moment magnitudes range of 4.5 ≤ $\mathbf{M}$ ≤ 8 calculated using the GS24bc model. As expected, GS24bc curves are not as smooth as for the backbone GS24b model shown in Figure 7. Figure 15 shows examples of different periods attenuation with distance to the fault differing in smoothness from those shown in Figure 8. Similar to the backbone model (Figure 8) attenuation rate changes from that of shear to surface waves at rupture distance $R_{rup}$ around 50 km.

As discussed before, we are considering the 3^rd dataset that includes 6045 data points, covering the range of moment magnitudes 4 ≤ $\mathbf{M}$ ≤ 7.9 and rupture distances $R_{rup}$ ≤ 250 km (M4_R250) to be the most important from the engineering applications point of view since it was created based on the recordings of earthquakes that can potentially produce damage. Figure 16 and Table 1 demonstrate the natural logarithmic standard deviation sigma ($\sigma$) for the GS24b, GS24bc and GK17 model predictions applied to the 3^rd dataset M4_R250. As expected, the newest GS24bc model demonstrates the lowest total sigma. However, previously developed GK17 sigma also demonstrates acceptable results.

We also split the total standard deviation ($\sigma$) into within-event ($\varphi$) and between-event ($\tau$) standard deviations.

$$\sigma =\sqrt{{\varphi }^2+{\tau }^2}$$

(14)

The within-event standard error ($\varphi$) is calculated by averaging within-event standard errors of earthquakes represented well enough in the database by more than 10 data points. Between-event standard deviation (τ) is calculated as

$$\tau =\sqrt{{\sigma }^2-{\varphi }^2}$$

(15)

Table 1 shows total, within-event and between-event standard deviations of the models using the 3^rd engineering dataset.

Figure 17, Figure 18 and Figure 19 demonstrate comparisons of the GS24b, GS24bc and GK17 ground motion models with empirical spectral accelerations $SA$ data from the total M4_R400 dataset for the three magnitude ranges:

• 4.77 ≤ $\mathbf{M}$ ≤ 5.21 with an average $\mathbf{M}$ = 4.99
• 5.70 ≤ $\mathbf{M}$ ≤ 6.24 with an average $\mathbf{M}$ = 6.06
• 7.28 ≤ $\mathbf{M}$ ≤ 7.9 with an average $\mathbf{M}$ = 7.64.

The GK17 model’s sigma recommended for use in hazard calculations published in Graizer (2018, Table 1) is lower, having been developed using subset of data with $\mathbf{M}$ ≥ 5.5 and $R_{rup}$ ≤ 140 km following the practice of NGA-West2 GMM developers. For example, Campbell and Bozorgnia (2014) recommended using sigma for the subset of data with $\mathbf{M}$ ≥ 5.5 and $R_{rup}$ ≤ 80 km. We believe that it is more appropriate for engineering applications in seismic hazard calculations to use the 3^rd subset of data shown in Figure 1 and Figure 2 for estimating sigma. Figure 16 also demonstrates $PGV$ sigma shown on the plot by a circle beyond the GMM range.

5. Results

The new GS24b backbone model is developed for the average RotD50 horizontal component of ground motion for shallow crustal earthquakes in active crustal regions. The model is derived based on a subset of the same NGA-West2 dataset [17] used in the development of the GK17 [14] model enhanced by the data from the three strong 2023 Turkish earthquakes with moment magnitudes $\mathbf{M}$ 6.3, 7.5 and 7.8 with a total of 13926 data points [20]. We did not use data from earthquakes with moment magnitudes $\mathbf{M}$ < 4.0 and limited the rupture distance range to 400 km. Records from earthquakes in Japan were also not used, because most of sites in Japan are characterized by subsurface geology significantly different from site conditions in other ACRs.

We consider the subset of 4 ≤ $\mathbf{M}$ ≤ 7.9 and $R_{rup}$ ≤ 250 km (M4_R250) of 6045 data points to be the most important from the engineering application point of view since it is limited to the data from earthquakes with magnitude-distance combinations that can potentially produce structural damage.

Similar to the GK17 [14], the new model has a bilinear attenuation slope of R_rup^-1 representing geometrical spreading of body waves for the closest 50 km from the fault, and R_rup^-0.5 at larger distances representing geometrical spreading of surface waves. The 50 km fault distance shift in the geometrical spreading is supported by the NGA-West2 data.

The number of input predictors in the GS24b model are limited to a few measurable parameters: moment magnitude $\mathbf{M}$, closest distance to fault rupture plane $R_{rup}$, time-averaged shear-wave velocity in the upper 30 m of the profile $V_{S30}$ and apparent anelastic attenuation SA quality factor $Q_{SA}\left(f,\mathbf{M}\right)$. This GS24b backbone GMM is applicable for earthquakes with 4.0 ≤ $\mathbf{M}$ ≤ 8.5, at rupture distances from 0 < $R_{rup}$ ≤ 400 km, at sites having 150 ≤ $V_{S30}$ ≤1500 m/s, and for spectral periods of 0.01 ≤ $T$ ≤ 10 sec.

A distinctive feature of the GS24b global backbone model is that it uses the closed form approximation of the spectral acceleration as a multiplication of the $PGA$ and spectral shape functions. This model can be later used for adjusting to the specific active crustal region, for example, Southern or Northern California, Italy or Turkey creating region specific partially non-ergodic models, requiring only residual and apparent anelastic attenuation adjustment.

As an example, we created GS24bc model by adjusting the generic active crustal region backbone GS24b model for moment magnitude M, shear-wave velocity $V_{S30}$ and fault rupture distance Rrup residuals using the NGA-West2 enhanced database. The GS24bc is an example of an easily created partially non-ergodic model demonstrating no residual biases. Compared to GK17 currently developed models also calculates PGV.

The previously developed GK17 [14] model tested against the engineering dataset M4_R250 demonstrated acceptable performance with standard error lower than the GS24b backbone model for periods T ≥ 0.1 s.

The GS24b and GS24bc ground motion models for spectral acceleration and peak ground velocity are developed using MATLAB software.

6. Data and Resources

The newly developed ground motion models discussed in this paper are based on the PEER Center Next Generation Attenuation Phase 2 (NGA-West2) database of processed ground motions from shallow crustal earthquakes in active tectonic regimes (http://ngawest2.berkeley.edu/). The NGA-West2 flatfile and the associated website were last accessed in July 2016. Turkish strong motion accelerometer data were retrieved from the Department of Earthquake, Disaster and Emergency Management Authority (AFAD; https://deprem.afad.gov.tr/stations) of Türkiye through the Turkish National Strong Motion Network (DOI: 10.7914/SN/TK), and Kandilli Observatory and Earthquake Engineering Research Institute of Boğaziçi University (KOERI, Istanbul; DOI: 10.7914/SN/KO) and processed by Buckreis et al. (2023a, 2023b) [18,19]. Other data used in this study came from the published sources listed in the text or in the references. All data processing was performed offline using a commercial software package (MATLAB R2024b, The MathWorks Inc., Natick, Massachusetts, 2024).