# Source Scaling and Ground-Motion Variability along the East Anatolian Fault

Dino Bindi, Riccardo Zaccarelli, Fabrice Cotton, Graeme Weatherill, Sreeram

Reddy Kotha

► To cite this version:

Dino Bindi, Riccardo Zaccarelli, Fabrice Cotton, Graeme Weatherill, Sreeram Reddy Kotha. Source Scaling and Ground-Motion Variability along the East Anatolian Fault. The Seismic Record, 2023, 3 (4), pp.311-321. 10.1785/0320230034. hal-04501989

HAL Id: hal-04501989

https://hal.science/hal-04501989

Submitted on 3 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Source Scaling and Ground-Motion Variability along the East Anatolian Fault

Dino Bindi<fnref n="1" />, Riccardo Zaccarelli<fnref n="1" />, Fabrice Cotton<fnref n="1" />, Graeme Weatherill<fnref n="1" />, and Sreeram Reddy Kotha²

## Abstract

We investigate the source scaling and ground-motion variability of 1585 earthquakes with $M_w$ > 3 occurring along the East Anatolian fault since 2010. We compile a dataset of 17,691 Fourier amplitude spectra of S waves recorded by 186 stations. A spectral decomposition is applied to isolate the source contribution from propagation and site effects. Source spectra are fit with Brune's model to estimate seismic moment and corner frequency and to compute the stress drop $\Delta\sigma$. The 10th, 50th, and 90th percentiles of the $\Delta\sigma$ distribution are 0.18, 0.51, and 1.69 MPa, respectively, and the average $\Delta\sigma$ increases with earthquake magnitude. For the two mainshocks of the 2023 sequence, the estimated $\Delta\sigma$ is about 13 MPa, significantly larger than the $\Delta\sigma$ of the smaller events. At intermediate and high frequencies, the interevent residuals are correlated with $\Delta\sigma$. When recorded peak ground accelerations and velocities for $M_w$ < 6 are compared with the predictions from ground-motion models proposed in the literature, the negative value of the average interevent residuals is consistent with low values of $\Delta\sigma$. Contrariwise, the average residuals for the peak parameter of the $M_w$ 7.8 and 7.5 mainshocks of the 2023 sequence are almost zero, but with distance dependencies.

**Cite this article as** Bindi, D., R. Zaccarelli, F. Cotton, G. Weatherill, and S. Reddy Kotha (2023). Source Scaling and Ground-Motion Variability along the East Anatolian Fault, *The Seismic Record*. **3**(4), 311–321, doi: 10.1785/0320230034.

## Introduction

On 6 February 2023, a seismic sequence struck Türkiye and Syria along the East Anatolian fault (EAF) zone in the Kahramanmaraş province. The sequence was initiated by an $M_w$ 7.8 earthquake along the Nurdağı-Pazarcık fault, located south of the EAF, and followed approximately 9 hr later by an $M_w$ 7.5 event involving the Sürgü and Çardak faults (Goldberg et al., 2023; Melgar et al., 2023). According to the Disaster and Emergency Management Authority (AFAD), the combined death toll in Türkiye and Syria exceeded 52,000, and the extension of the affected area contributed to generate a widespread damage that exceeded expectations. The two mainshocks were recorded by a large number of accelerometers belonging to the AFAD network (see Data and Resources), allowing the characteristics of the ground shaking to be studied even at close source distances. The aim of this study is to relate the differences between observed and predicted ground motions to the distribution of source parameters, and in particular to the scaling of the stress drop with the seismic

moment. The strategy we follow is the following: first, we estimate the source parameters of events with magnitude greater than 3 occurring since 2010 using a spectral decomposition approach; then, we calibrate a regional ground-motion model (GMM) that allows us to capture the spectral scaling with distance and magnitude in eastern Türkiye, partitioning the residuals into event- and station-specific components; finally, to relate the results of the ground-motion residual analyses to the source parameter distributions, we quantify the correlation

<footnote n="1">1. German Research Centre for Geoscience GFZ, Potsdam, Germany, https://orcid.org/0000-0002-8619-2220 (DB); https://orcid.org/0000-0001-9133-1750 (RZ); https://orcid.org/0000-0002-9242-3996 (FC); https://orcid.org/0000-0001-9347-2282 (GW); 2. Université Grenoble Alpes, Institut des Sciences de la Terre, Grenoble, France, https://orcid.org/0000-0002-4874-3730 (SRK)</footnote>

<footnote>*Corresponding author: bindi@gfz-potsdam.de</footnote>

<footnote>© 2023. The Authors. This is an open access article distributed under the terms of the CC-BY license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.</footnote>
<figure><img src="image_3.png" /></figure>

between interevent residuals and stress drop at intermediate and high frequencies.

## Data

Waveforms are downloaded from AFAD and European Integrated Data Archive repositories (see Data and Resources) using the *stream2segment* software (Zaccarelli et al., 2019). Download is guided by the AFAD seismic catalog, setting the minimum magnitude to 3 for events occurring in the study area (Fig. 1) before February 2023 and to 4 for events occurring later. The maximum selected depth is 30 km, and only recordings at hypocentral distances smaller than 250 km are considered. The dataset is mostly composed of strong motion recordings, but broadband recordings are also collected for events with a magnitude below 4. The processing workflow includes pass-band filtering with magnitude-dependent high-pass ($f_{hp}$) and low-pass ($f_{lp}$) corner frequencies (i.e., $f_{hp}$ = 0.3 Hz for magnitude below 6 and $f_{hp}$ = 0.08 Hz for larger events; $f_{lp}$ = 40 Hz for all events), correction for the instrumental response, computation of Fourier amplitude spectra (FAS) for windows selected using an energy criteria, and smoothing with a triangular window. The energy criterion is applied by selecting the following distance-dependent percentiles for the cumulative of the squared velocity: 0.025–0.95 for hypocentral distances $R < 25$ km; 0.05–0.85

<caption>**Figure 1.** Analyzed dataset for the East Anatolian fault zone. Locations (circles) are taken from the Disaster and Emergency Management Authority catalog; moment magnitudes $M_w$ and stress drop $\Delta\sigma$ are derived in this study. Black triangles indicate the location of the considered stations. The focal mechanisms are downloaded from the Global Centroid Moment Tensor Project. See Data and Resources, for details. The inset map shows the location of the study area (black box) within the Euro-Mediterranean region</caption>

for $25 \le R < 50$ km; 0.05–0.8 for $50 \le R < 100$ km; 0.05–0.75 for $R \ge 100$ km; the minimum duration is fixed to 5 s and the maximum to 30 s for $R < 20$ km, to 40 s for $20 \le R < 80$ km, and to 50 s for $R \ge 80$ km. Only spectra with a signal-to-noise ratio larger than 5 for frequencies within the range 1–6 Hz and larger than 2 outside this range are considered for the analysis. We consider events and stations with at least three recordings. The final dataset, shown in Figure 1, is composed of 1585 earthquakes, 186 stations, and 17,691 FAS, computed as the square root of the sum of the two horizontal spectra squared.

## Ground-Motion Residuals

We use the ground shaking intensity model (eGSIM) service (see Data and Resources) to compute the observed minus predicted residuals for peak ground acceleration (PGA) and
velocity (PGV), considering the previous derived GMM for Europe. For events with a magnitude smaller than 7, predictions are computed using the (Bindi et al., 2014) GMM. We select a GMM based on a simple functional form to avoid the need to use a priori values for unknown variables. We also select a GMM based on the hypocentral distance. We set the shear-wave velocity averaged in the uppermost 30 m to 420 m/s, that is, the median value for stations that recorded the main-shock at distances within 300 km. The mean ± one standard deviation of the residuals are (−0.94 ± 0.83) and (−0.47 ± 0.70) for the natural logarithm of PGA and PGV, respectively, indicating that PGA and PGV values were overestimated, on average, by predictions. To investigate the magnitude dependence of the overestimation, Figure 2a,b shows the interevent standardized residuals against magnitude. The negative mean values obtained for both PGA and PGV are mostly controlled by magnitudes smaller than 5, whereas residuals for larger magnitudes are approaching zero. Therefore, we also evaluate the residuals for the M<sub>w</sub> 7.8 and 7.5 events. For the two main-shocks, we analyze the peak parameters and use the site information provided for these events by the Engineering Strong Motion service (see Data and Resources) that disseminates processed data revised by experts. To account for the extension of the fault, residuals are computed using the European model developed by Kotha et al. (2016), here applied considering the regional attributes for Türkiye and the Joyner–Boore distances R<sub>JB</sub> (i.e., the minimum horizontal distance to the surface projection of the fault, and set to zero for those stations located within the surface projection). Figure 2c,d shows the PGA and PGV standardized residuals versus distance for the M<sub>w</sub> 7.8 event; Figure 2e,f for the M<sub>w</sub> 7.5 event, considering R<sub>JB</sub> < 300 km. For M<sub>w</sub> 7.8, the average residual for both PGA and PGV are close to zero but show distance dependencies. For R<sub>JB</sub> < 100 km, the residuals are mostly positive (i.e., observations larger than prediction), whereas the residuals decrease and become negative at larger distances. A similar pattern is also observed for M<sub>w</sub> 7.5, although in this case the average standardized residuals for PGV are positive (0.37). The rapid attenuation observed for the EAF region is consistent with the adaptation of the European backbone model (Kotha et al., 2020), as proposed by Weatherill et al. (2020) for probabilistic hazard assessment in Europe. It is worth noting that more recent European models (e.g., Kotha et al., 2020) may better capture the regional variability of ground shaking and near-source ground motion. Because the main objective of the residual analysis performed in this

study is to confirm whether the ground motions generated by small and large events behave differently compared to the predictions, a comprehensive residual analysis considering more complex models requiring regional adaptations is beyond the current scope of the study.

In summary, the ground-motion residual analyses show that although the ground motions generated by the main-shocks, measured in terms of PGA and PGV, are on average consistent with expectations from GMMs, the peak parameters predicted for the small-to-moderate events overestimate the observations. To relate the results of the residual analysis to source effects, in the next section we use the source spectra isolated from propagation and site effects to evaluate the scaling of stress drop with seismic moment.

## Spectral Decomposition

We apply spectral decomposition (Bindi et al., 2020) to split source, propagation, and site amplification effects. The decomposition is performed for each frequency separately using a nonparametric approach described by the following:

$$ \log\text{FAS}_{es} = \log S_e + \log P(R_{es}) + \log Z_s, \quad (1) $$

in which FAS<sub>es</sub> is the spectral amplitude of event e recorded by station s, P(R<sub>es</sub>) describes the propagation effects between event e and station s; R<sub>es</sub> is the hypocentral distance; S<sub>e</sub> is the source spectral amplitude of earthquake e; Z<sub>s</sub> represents the site amplification effects at station s; and log indicates the logarithm in base 10. Considering all available e and s combinations, equation (1) generates an overdetermined linear system that we solved in a least-squares sense. To remove two unresolved degrees of freedom, we constrain the propagation term to one at 10 km, and we assume a reference site amplification model.

We first perform the spectral decomposition focusing on the propagation term. Therefore, we select events with magnitudes in the AFAD catalog above 4, and we constrain the average of site amplification over the whole set of stations to one. The distance range is discretized, with nodes placed every 4 km for R<sub>es</sub> < 102 km; every 6 km for 102 < R<sub>es</sub> < 220 km; and every 10 km for R<sub>es</sub> > 220 km. The spectral attenuation is linearized between two consecutive nodes and smoothed by requiring its second derivative with distance to be small. The obtained spectral attenuation P(R<sub>es</sub>,f) for the considered frequencies f is shown in Figure 3b and provided as a table in the supplemental material. After correcting the spectral
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 2.** (a) Interevent residuals for peak ground acceleration (PGA) for *M*<sub>w</sub> < 7 events and using the (Bindi et al., 2014) ground-motion model (GMM). (b) The same as in panel (a), but for peak ground velocity (PGV). Residuals for (c) PGA and (d) PGV for the *M*<sub>w</sub> 7.8 event using the parametric information provided by the Engineering Strong</figcaption></figure>

Motion database (see Data and Resources) and considering the Kotha *et al.* (2016) model. (e,f) The same as in panels (c) and (d), but for the *M*<sub>w</sub> 7.5 earthquake. In all the panels, dotted lines indicate the mean residual, and the trend lines in black are the result of local regressions.
<figure><img src="image_10.png" /><figcaption>**Figure 3.** Results of the spectral decomposition. (a) Source spectra (events with $M_w > 6$ are shown with different colors). (b) Spectral attenuation curves (the average attenuation computed over different frequency ranges is shown using different colors as in the legend); models $1/R$ and $1/R^2$ are shown for comparison. (c) Site amplification terms: site amplifications at selected reference stations are shown in blue, with their average constrained to the amplification model shown in red; the green spectrum is the amplification obtained for station 4628 of the TK network (see the Spectral Decomposition section, for more details).</figcaption></figure>

amplitudes for $P(R_{es}, f)$, the linear system generated by equation (1), modified by considering $FAS_{es}/P(R_{es})$ on the left-hand side, is solved to separate the source spectra $S_e$ from the site amplifications $Z_s$. This second decomposition is performed considering all events but limiting the distances to 110 km to increase the signal-to-noise ratio for events with magnitudes below 4 and to mitigate the impact of propagation complexities affecting longer paths.

The assumed reference amplification model (Fig. 3c, red line) corresponds to a generic crustal amplification model multiplied by a negative exponential term controlled by $k_0 = 0.035$ s (Bindi et al., 2021), in which the $k_0$ value is empirically chosen to obtain, on average, flat source acceleration spectra as predicted by the considered $\omega^2$ source model (Brune, 1970). The reference amplification model is used to constrain the average amplification of a set of reference stations. The reference stations (Fig. 3c, blue lines), which are selected as stations with small amplification compared to the network average, are nearly frequency independent. To choose these, we inspect the results of a preliminary inversion performed by setting the average amplification over all stations to one, regardless the frequency, and select stations with an overall amplification median < 0.5 and 95th percentile < 1.5.

The obtained site amplifications $Z_s$ are shown in Figure 3c (black); the green curve exemplifies the result for station 4628 (Kahramanmaraş) belonging to the AFAD network TK (see Data and Resources). The AFAD site characterization performed for station 4628 indicated that the shear-wave velocity averaged over the uppermost 30 m is $V_{S30} = 186$ m/s. The resonant peak of amplification at 1.8 Hz, as evidenced by the microtremor horizontal to vertical spectral ratios (H/V) for response spectra (see Data and Resources), is in agreement with the peak also shown by the spectral decomposition result. Contrariwise, the decomposition result shows a broad range of amplifications above the fundamental resonance not captured by the H/V ratios for microtremors, suggesting possible 2D–3D amplification effects also affecting the vertical component.

The nonparametric source spectra (Fig. 3a) are fit to the Brune model (Brune, 1970):

$$S(f) = \frac{\Omega^0}{[1 + (f/f_c)^2]}, \qquad (2)$$

in which $f_c$ is the corner frequency, and the low-frequency spectral level $\Omega_0$ is proportional to the seismic moment $M_0$. To make the moment magnitudes compatible with those in
<figure><img src="image_11.png" /><figcaption>**Figure 4.** Scaling of source parameters. (a) Seismic moment versus corner frequency color with colors proportional to the root-mean-square error (RMSE) of the source spectral fit. (b) Stress drop values versus seismic moment color coded with hypocentral depth; the trend line is determined by performing a local regression (loess). (c) Stress drop values against hypocentral depth, with colors representing the moment magnitude. In all the panels, the red squares indicate the results obtained for the $M_w$ 7.8 and 7.5 mainshocks using equation (9).</figcaption></figure>

the AFAD catalog, we first fit the source spectra to get the values of the low-frequency plateau level $\Omega_0$ and $f_c$ for events with a magnitude below 6; then, the $\Omega_0$ values are transformed to $M_0$ and the average offset with respect to the values obtained from the catalog; finally, the source fit is repeated over the spectra corrected for the offset and constraining $M_0$ to the values obtained from the catalog when available (i.e., for those events with $M_w$ provided by the AFAD catalog only $f_c$ is determined). The $M_0$ and $f_c$ values (Fig. 4a; supplemental information) are used to compute the stress drop $\Delta\sigma$ for a circular rupture with uniform $\Delta\sigma$ (Eshelby, 1957):

$$ \Delta\sigma = \frac{7 M_0}{16 r^3}, \qquad (3) $$

$$ r = \frac{0.37\beta}{f_c}, \qquad (4) $$

in which we considered the shear-wave velocities $\beta$ at the hypocentral depths extracted from the 1D velocity model proposed by Güvercin et al. (2022). The scaling of $\Delta\sigma$ with $M_0$ and depth is shown in Figure 4b and 4c, respectively.

Regarding the mainshocks, the extension and complexity of the rupture processes cannot be described by a simple circular rupture model. Nevertheless, to provide a first-order estimate of $\Delta\sigma$ also for the two main events consistent with the values provided for smaller events, we consider an elliptical rupture model with semiaxes $a$ and $b$, that is,

$$ \Delta\sigma = \frac{3}{4} \eta \frac{M_0}{\pi a b^2}, \qquad (5) $$

in which $\eta$ is a function of Poisson's modulus (here assumed equal to 0.25) of the aspect ratio $a/b$ and of the complete elliptical integrals of the first and second kinds (Eshelby, 1957). In equation (5), we compute the semiaxis as $a = L/2$ and $b = W/2$, in which $L$ and $W$ are the length and width of the rupture, respectively. Using $L = 300$ km and $W = 20$ km for the event with magnitude $M_w$ 7.8, and $L = 110$ km and $W = 20$ km for the $M_w$ 7.5 event (Goldberg et al., 2023), we get $\Delta\sigma = 13$ MPa for both the events (shown as red squares in Fig. 4). Although these estimates are based on an oversimplified model, and alternative estimates of $L$ and $W$ would produce different values, the obtained values are consistent with the trend shown by the scaling in Figure 4b, which suggests $\Delta\sigma$ values for the two
main events of the order of 10 MPa. Stress drop values of the same order, that is, between 8 and 10 MPa, were also estimated for the two large strike slip events that occurred in 1999 along the North Anatolian fault—the $M_w$ 7.6 Izmit and $M_w$ 7.2 Düzce earthquakes (Akinci et al., 2006).

## GMM

We perform a mixed effects regression (Bates et al., 2015) to quantify source-related deviations from the median predictions of a GMM capturing the regional scaling with magnitude and distance in eastern Türkiye. For each frequency, we consider the following model:

$$ \text{LogFAS} = F(M) + G(R) + Q(R) + \delta B_e + \delta S_s + \epsilon, \quad (6) $$

in which the magnitude scaling $F(M)$ is given by

$$ F(M) = e_1 + b_1(M - 4.5) + b_2(M - 4.5)^2, \quad (7) $$

and the distance-dependent attenuation scaling is defined by the geometrical spreading term:

$$ G(n_1, n_2, n_3; \text{Log}(R), R_a, R_b) = \begin{cases} n_1 \text{Log}(\frac{R}{10}) & \text{if } R \le R_a \\ n_1 \text{Log}(\frac{R_a}{10}) + n_2 \text{Log}(\frac{R}{R_a}) & \text{if } R_a < R \le R_b, \\ n_1 \text{Log}(\frac{R_b}{10}) + n_2 \text{Log}(\frac{R_b}{R_a}) + n_3 \text{Log}(\frac{R}{R_b}) & \text{otherwise} \end{cases} \quad (8) $$

and by the anelastic attenuation term:

$$ Q(k; R, R_a, R_b) = \begin{cases} 0 & \text{if } R < R_a \\ k_1(R - R_a)/100 & \text{if } R_a \le R < R_b, \\ k_1(R_b - R_a)/100 + k_2(R - R_b)/100 & \text{otherwise} \end{cases} \quad (9) $$

in which the fixed effects $e_1, b_1, b_2, c_1, c_2, c_3, k_1, k_2$ define the median predictions for different $M_w$ and R scenarios, and the hinge distances are set to $R_a = 8$ km and $R_b = 60$ km. The random effects $\delta S_s$ and $\delta B_e$ (Atik et al., 2010) are introduced to partition the residuals according to station (interstation) and event (interevent) grouping factors, respectively, whereas the leftover distribution $\epsilon$ absorbs contributions to the residuals that are not event or station specific (e.g., path or radiation pattern effects). The interstation and interevent

residuals approximate normal distributions with zero mean and standard deviation $\phi_s$ and $\tau$, respectively. Because the random effects are independent distributions, the standard deviation $\sigma$ of the complete residual distribution is given by $\sigma = \sqrt{\tau^2 + \phi_s^2 + \phi_0^2}$, in which $\phi_0$ is the standard deviation of $\epsilon$. The values of the standard deviations for the analyzed frequencies are shown in Figure 5a. In agreement with the previous studies (e.g., Bindi et al., 2020; Kotha et al., 2022), the largest contribution to ground-motion variability $\sigma$ is given by the site component $\phi_s$, in particular at high frequencies.

## Results and Discussions

### Source parameters

The source scaling (Fig. 4b) shows that the average $\Delta\sigma$ increases with $M_0$. For events with magnitudes below 4, the 10th, 50th, and 90th percentiles of the $\Delta\sigma$ distribution are 0.15, 0.37, and 1.09 MPa, respectively, and they increase to 0.32, 0.79, and 2.08 MPa for magnitudes between 4 and 5. The exponent $\epsilon = 0.13$ found for the $M_0 \propto f_c^{-(3+\epsilon)}$ relationship (Kanamori and Rivera, 2004) confirms the weak departure from self-similarity. The increase of $\Delta\sigma$ with $M_0$ with the tendency of saturating for large magnitudes has already been observed for sequences in other regions, such as the 2019 Ridgecrest in southern California (Bindi et al., 2020; Trugman, 2020) and in central Italy (Morasca et al., 2022). The average of the $\Delta\sigma$ distribution is 0.51 MPa; the standard deviation of 0.38 for $\log \Delta\sigma$ is within the range of 0.35–0.45 found in the recent studies (e.g., Baltay et al., 2013; Trugman, 2020). Although comparing absolute values obtained by different studies is hindered by the impact of different methodologies and model assumptions (e.g., Shearer et al., 2019; Bindi et al., 2023), $\Delta\sigma$ for magnitudes below 5 are, on average, lower than the values obtained for the 2019 Ridgecrest sequence in California. Using the same methodology applied in this study, Bindi et al. (2020) obtained 1.2, 2.9, and 7.2 MPa for the 10th, 50th, and 90th $\Delta\sigma$ percentiles of events with $M < 4$, and 2.6, 5.0, and 9.3 MPa for $4 < M < 5$, that is, values about one order of magnitude larger than those found in the present study for the EAF zone.

Considering the competing effects between source and depth-dependencies of attenuation, the possibility of detecting trends showing an increase of $\Delta\sigma$ with hypocentral depth is debated in the literature (Abercrombie et al., 2021; Bindi et al., 2021; Pennington et al., 2021). The distribution of $\Delta\sigma$ shown in Figure 4c does not show any trend with depth, and events with similar magnitude but occurring at different depths show similar stress drop values. The cluster of events
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

with depth at 7 km is probably an artifact of the location procedure, and more accurate locations and velocity models are needed to better characterize the $\Delta\sigma$ distribution with depth.

## Ground-motion residuals

For a given seismic moment, $\Delta\sigma$ determines the value of the corner frequency and of the source spectral level at high frequencies. Using the random vibration theory (Hanks and McGuire, 1981), it has been shown that for an $\omega^2$ source model, the natural logarithm of PGA is expected to increase with the logarithm of $\Delta\sigma$ with slope 0.8 and with the logarithm

**Figure 5** (a) Standard deviations of the residual distributions computed for a GMM derived in the Fourier domain (see the GMM section, for details); red circles indicate the standard deviation $\tau$ of the model derived by (Bindi et al., 2020) for the 2019 Ridgecrest sequence; blue circles indicate the standard deviation $\tau_e$ of the model derived by (Kotha et al., 2022) for crustal earthquakes in Europe. (b) Frequency dependence of Pearson's correlation coefficients computed between interevent residuals $\delta B_e$ and stress drop $\Delta\sigma$. (c) Interevent residuals for 0.52 Hz against stress drop. (d) Interevent residuals for 5.2 Hz against stress drop. The trend lines in panels (c) and (d) are the results of a loess regression.

of the rupture velocity with slope 2.4 (see equation 3 in Causse and Song, 2015, and references therein). Because the median model (equation 6) is describing the scaling with distance and
magnitude, it is not capturing the impact on high-frequency spectral values of $\Delta\sigma$ variations for events with similar $M_0$. Therefore, it is expected that the interevent residuals $\delta B_e$ absorb the $\Delta\sigma$ variability, although a formal connection between $\delta B_e$ and source parameters is not straightforward (Picozzi et al., 2021, their equation 8).

Because ground-shaking variability is mostly influenced by local site effects, to highlight the source contribution, we partitioned the residuals into interevent and interstation components (equation 6). Figure 5a confirms that the standard deviation $\phi_s$ of the interstation residuals is the largest contributor to the overall variability at high frequencies, as expected by the large variability of site amplifications shown in Figure 3c. Figure 5a also compares the standard deviation $\tau$ of the interevent residual with the result obtained by Bindi et al. (2020) for the 2019 Ridgecrest sequence and by Kotha et al. (2022) for Europe, where $\tau_e$ is also accounting for the between-locality variability $\tau_{L2L}$. The three datasets show a similar spectral dependency of $\tau$, with the highest values obtained for the models derived at the continental spatial scale and the lowest ones for the very localized Ridgecrest dataset.

To quantify the correlation between $\Delta\sigma$ and the interevent residuals $\delta B_e$, we compute the Pearson correlation coefficient considering $\delta B_e$ for different frequencies. The correlation increases with frequency, assuming values greater than 0.5 for frequencies >2 Hz (Fig. 5b). The strongest correlation is reached between 6 and 10 Hz, in which the correlation coefficient is about 0.75. For frequencies higher than 10 Hz, the correlation coefficient decreases, suggesting that parameters other than stress drop also contribute to the variability of the source spectra at high frequencies, resulting in event-dependent deviations from an $\omega^2$ spectral shape. The frequency dependency of the correlation between $\delta B_e$ and $\Delta\sigma$ distributions is exemplified in Figure 5c,d, in which the results for 0.52 (weak correlation) and 5.2 Hz (strong correlation) are visualized.

## Conclusions

We determined the source parameters of 1585 earthquakes with $M_w > 3$ occurring in the East Anatolian fault region since 2010. We estimated an average stress drop $\Delta\sigma$ equal to 0.51 MPa, with the 10th and 90th percentiles equal to 0.18 and 1.69 MPa, respectively. The $\Delta\sigma$ distribution shows a positive trend with seismic moment and a significant correlation with interevent residuals evaluated at intermediate and high frequencies (Pearson's correlation coefficient > 0.6 above 2 Hz). When observed PGAs and PGVs are compared with predictions from

the previously derived GMM for Europe and the Middle East, the negative average interevent residuals obtained for magnitude < 6 are consistent with $\Delta\sigma$ values lower than those estimated for the European events contributing to the GMM (Bindi and Kotha, 2020, their fig. 11, median $\Delta\sigma = 2.8$ MPa for magnitudes in the range 3.5–5.5). The $\Delta\sigma$ values for $M_w < 6$ are also lower than those obtained for the 2019 Ridgecrest sequence in California by Bindi et al. (2020) using the same approach adopted in this study (median $\Delta\sigma = 5$ MPa for $4 \le M_w \le 5$). The average residuals for PGA and PGV for the $M_w$ 7.8 and 7.5 mainshocks are almost zero but show distance dependencies. Although the recent models developed to better capture the regional variability of ground shaking in Europe (e.g., Kotha et al., 2020) could improve the overall precision of the predictions compared to the results obtained this study, the differences in residual distributions for small and large events agree with the main conclusions drawn for the stress drop. The $\Delta\sigma = 12$ MPa estimated for the two mainshocks is in the same order as the values between 8 and 10 MPa estimated for the 1999 $M_w$ 7.6 Izmit and $M_w$ 7.2 Düzce earthquakes estimated by Akinci et al. (2006), both strike-slip events occurring over the North Anatolian fault. For the 2020 $M_w$ 6.8 Sivrice-Elazığ earthquake that occurred about 200 km northeast of the 2023 $M_w$ 7.8 earthquake (Cheloni and Akinci, 2020) estimated $\Delta\sigma = 9$ MPa; for the 2011 $M_w$ 7.1 Van earthquake (eastern Türkiye) that featured reverse faulting (Akinci et al., 2014) estimated $\Delta\sigma = 20$ MPa. The consistency of the stress drop values and the fact that the interevent residuals for these two events are close to zero confirm that the seismic motion generated by the two mainshocks of the 2023 sequence was not, on average, unusually high for this size of earthquake.

## Data and Resources

Softwares used in this study are R language (R Core Team, 2020), Generic Mapping Tools (Wessel et al., 2013), stream2segment (https://github.com/rizac/stream2segment), and ground shaking intensity model (eGSIM) (https://github.com/rizac/eGSIM). Data repositories used in this study are the Engineering Strong Motion (ESM) database (https://esm-db.eu/#/home), the European Integrated Data Archive (EIDA; https://www.orfeus-eu.org/data/eida/), the Disaster and Emergency Management Authority (AFAD) catalog (https://deprem.afad.gov.tr/event-catalog), and AFAD waveforms (https://deprem.afad.gov.tr/fdsnws/dataselect/1/). The data about the AFAD station report are available at https://tadas.afad.gov.tr/station-detail/385. Seismic networks used in
the study are Turkish National Strong Motion Network TK (DOI: 10.7914/SN/TK), Kandilli Observatory and Earthquake Research Institute (KOERI) KO (DOI: 10.7914/SN/KO), and Turkish National Seismic Network TU (DOI: 10.7914/SN/TU). Focal mechanisms shown in Figure 1 are available at the Global Centroid Moment Tensor (Global CMT) Project (https://www.globalcmt.org/). The supplemental material to this article include the following files: table reporting the spectral attenuation model; table including the nonparametric source spectra; table including the site amplifications; table with the source parameters computed in this study; list of the selected reference stations and the reference amplification model; and table with the coefficients of the calibrated ground-motion model (GMM). Files are distributed as Bindi (2023). Supplements to source scaling and ground-motion variability along the East Anatolian fault, Zenodo, DOI: 10.5281/zenodo.10044523. All websites were last accessed in September 2023.

**Declaration of Competing Interests**

The authors declare no competing interests.

**Acknowledgments**

The authors thank Annemarie Baltay and Aybige Akinci for their valuable comments.

**References**

Abercrombie, R. E., D. T. Trugman, P. M. Shearer, X. Chen, J. Zhang, C. N. Pennington, J. L. Hardebeck, T. H. W. Goebel, and C. J. Ruhl (2021). Does earthquake stress drop increase with depth in the crust? *J. Geophys. Res.* **126**, no. 10, e2021JB022314, doi: 10.1029/2021JB022314.

Akinci, A., L. Malagnini, and R. B. Hermann (2014). High-frequency attenuation in the Lake Van region, eastern Turkey, *Bull. Seismol. Soc. Am.* **104**, no. 3, doi: 10.1785/0120130102.

Akinci, A., L. Malagnini, R. B. Hermann, R. Gok, and M. Sorensen (2006). Ground motion scaling in the Marmara region, Turkey, *Geophys. J. Int.* **166**, 635–651.

Atik, L. A., N. Abrahamson, J. J. Bommer, F. Scherbaum, F. Cotton, and N. Kuehn (2010). The variability of ground-motion prediction models and its components, *Seismol. Res. Lett.* **81**, no. 5, 794–801.

Baltay, A. S., T. C. Hanks, and G. C. Beroza (2013). Stable Stress-drop measurements and their variability: Implications for ground-motion prediction, *Bull. Seismol. Soc. Am.* **103**, no. 1, 211–222.

Bates, D., M. Mächler, B. Bolker, and S. Walker (2015). Fitting linear mixed-effects models using lme4, *J. Stat. Softw.* **67**, no. 1, 1–48.

Bindi, D., and S. R. Kotha (2020). Spectral decomposition of the Engineering Strong Motion (ESM) flat file: Regional attenuation, source scaling and Arias stress drop, *Bull. Earthq. Eng.* **18**, 2581–2606.

Bindi, D., M. Massa, L. Luzi, G. Ameri, F. Pacor, R. Puglia, and P. Augliera (2014). Pan-European ground-motion prediction equations for the average horizontal component of PGA, PGV, and 5 %-damped PSA at spectral periods up to 3.0 s using the RESORCE dataset, *Bull. Earthq. Eng.* **12**, 1573–1456.

Bindi, D., H. N. T. Razafindrakoto, M. Picozzi, and A. Oth (2021). Stress drop derived from spectral analysis considering the hypocentral depth in the attenuation model: Application to the Ridgecrest region, California, *Bull. Seismol. Soc. Am.* **111**, no. 6, 3175–3188.

Bindi, D., D. Spallarossa, M. Picozzi, A. Oth, P. Morasca, and K. Mayeda (2023). The community stress-drop validation study—Part II: Uncertainties of the source parameters and stress drop analysis, *Seismol. Res. Lett.* **94**, no. 4, 1992–2002.

Bindi, D., R. Zaccarelli, and S. R. Kotha (2020). Local and moment magnitude analysis in the Ridgecrest region, California: Impact on interevent ground-motion variability, *Bull. Seismol. Soc. Am.* **111**, no. 1, 339–355.

Brune, J. N. (1970). Tectonic stress and the spectra of seismic shear waves from earthquakes, *J. Geophys. Res.* **75**, no. 26, 4997–5009.

Causse, M., and S. G. Song (2015). Are stress drop and rupture velocity of earthquakes independent? Insight from observed ground motion variability, *Geophys. Res. Lett.* **42**, 7383–7389.

Cheloni, D., and A. Akinci (2020). Source modelling and strong ground motion simulations for the 24 January 2020, Mw 6.8 Elazığ earthquake, Turkey, *Geophys. J. Int.* **223**, no. 2, 1054–1068, doi: 10.1093/gji/ggaa350.

Eshelby, J. D. (1957). The determination of the elastic field of an ellipsoidal inclusion, and related problems, *Proc. R. Soc. Lond.* **241**, no. 1226, 376–396.

Goldberg, D. E., T. Taymaz, N. G. Reitman, A. E. Hatem, S. Yolsal-çevikbilen, W. D. Barnhart, T. S. Irmak, D. J. Wald, T. Öcalan, W. L. Yeck, et al. (2023). Rapid characterization of the February 2023 Kahramanmaraş, Türkiye, earthquake sequence, *Seismic Record 3*, no. 2, 156–167.

Güvercin, S. E., H. Karabulut, A. O. Konca, U. Doğan, and S. Ergintav (2022). Active seismotectonics of the East Anatolian fault, *Geophys. J. Int.* **230**, no. 1, 50–69.

Hanks, T. C., and R. K. McGuire (1981). The character of high-frequency strong ground motion, *Bull. Seismol. Soc. Am.* **71**, no. 6, 2071–2095.

Kanamori, H., and L. Rivera (2004). Static and dynamic scaling relations for earthquakes and their implications for rupture speed and stress drop, *Bull. Seismol. Soc. Am.* **94**, no. 1, 314–319.

Kotha, S. R., D. Bindi, and F. Cotton (2016). Partially non-ergodic region specific GMPE for Europe and Middle-East, *Bull. Earthq. Eng.* **14**, 1573–1456.

Kotha, S. R., D. Bindi, and F. Cotton (2022). A regionally adaptable ground-motion model for Fourier amplitude spectra of shallow crustal earthquakes in Europe, *Bull. Earthq. Eng.* **20**, 1573–1456.

Kotha, S. R., G. Weatherill, D. Bindi, and F. Cotton (2020). A regionally-adaptable ground-motion model for shallow crustal earthquakes in Europe, *Bull. Earthq. Eng.* **18**, 4091–4125.
Melgar, D., T. Taymaz, A. Ganas, B. Crowell, T. Öcalan, M. Kahraman, V. Tsironi, S. Yolsal-çevikbilen, S. Valkaniotis, T. S. Irmak, et al. (2023). Sub- and super-shear ruptures during the 2023 Mw 7.8 and Mw 7.6 earthquake doublet in SE Türkiye, *Seismica 2*, no. 3, doi: 10.26443/seismica.v2i3.387.

Morasca, P., D. Bindi, K. Mayeda, J. Roman-Nieves, J. Barno, W. R. Walter, and D. Spallarossa (2022). Source scaling comparison and validation in central Italy: Data intensive direct S waves versus the sparse data coda envelope methodology, *Geophys. J. Int.* **231**, 1573–1590.

Pennington, C. N., X. Chen, R. E. Abercrombie, and Q. Wu (2021). Cross validation of stress drop estimates and interpretations for the 2011 Prague, OK, earthquake sequence using multiple methods, *J. Geophys. Res.* **126**, no. 3, e2020JB020888, doi: 10.1029/2020JB020888.

Picozzi, M., F. Cotton, D. Bindi, A. Emolo, G. M. Adinolfi, D. Spallarossa, and A. Zollo (2021). Spatiotemporal evolution of ground-motion intensity at the Irpinia near-fault observatory, southern Italy, *Bull. Seismol. Soc. Am.* **112**, no. 1, 243–261.

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

Shearer, P. M., R. E. Abercrombie, D. T. Trugman, and W. Wang (2019). Comparing EGF methods for estimating corner frequency and stress drop from P wave spectra, *J. Geophys. Res.* **124**, no. 4, 3966–3986.

Trugman, D. T. (2020). Stress-drop and source scaling of the 2019 Ridgecrest, California, earthquake sequence, *Bull. Seismol. Soc. Am.* **110**, no. 4, 1859–1871.

Weatherill, G., S. Kotha, and F. Cotton (2020). A regionally-adaptable “Scaled Backbone” ground motion logic tree for shallow seismicity in Europe: Application to the 2020 European Seismic hazard model, *Bull. Earthq. Eng.* **18**, 5087–5117.

Wessel, P., W. H. F. Smith, R. Scharroo, J. Luis, and F. Wobbe (2013). Generic mapping tools: Improved version released, *EOS Trans. AGU* **94**, no. 45, 409–410.

Zaccarelli, R., D. Bindi, A. Strollo, J. Quinteros, and F. Cotton (2019). Stream2segment: An open-source tool for downloading, processing, and visualizing massive event-based seismic waveform datasets, *Seismol. Res. Lett.* **90**, no. 5, 2028–2038.

Manuscript received 2 September 2023

Published online 15 November 2023