Engineering attributes of ground motions from February 2023 Türkiye earthquake sequence

Abstract

The 2023 Türkiye earthquake sequence includes the 6 February M7.81 mainshock followed approximately 9 h later by an M7.74 event, and many smaller aftershocks including M6.81 and M6.37 events on 6 and 20 February, respectively. These events occurred in a region near the plate boundary of the East Anatolian Fault, in the proximity of which numerous ground motion recording stations had been installed north of the Türkiye–Syria border. As a result, the events were well recorded both near the fault and at rupture distances up to 582 km. We describe the available recordings and component-specific data processing performed with the aim of optimizing usable bandwidth. The resulting database includes 337, 365, 284, and 229 usable three-component recordings from the M7.81, M7.74, M6.81, and M6.37 events, respectively. We also present source, path, and site metadata that were compiled according to uniform protocols. Comparisons to a global ground motion model (GMM) for active tectonic regions and a local, Türkiye-specific model demonstrate the existence of complex path effects that result in relatively poor fits between the GMMs and observed data at large distances (generally R _JB  > 200 km). Under-predictions at some stations may be influenced by directivity and/or basin effects that affect the ground motions but that are not accounted for directly in the GMMs. We present analysis of spatial variability of several intensity measures. The residual maps produced from these analyses demonstrate that the global GMM overpredicts on the Anatolian block and underpredicts on the Arabian block, which is suggestive of distinct crustal attenuation features.

Introduction and project motivation

The 2023 Türkiye earthquake sequence occurred in a region known to have active faults and that was extensively instrumented north of the Türkiye–Syria border. Multiple large events were well recorded both near the fault and at distances up to 582 km. We describe available recordings as of the present date (June 2024); describe manual, component-specific data processing that was performed to optimize usable bandwidth; present metadata compiled using uniform protocols; compare data to a global ground motion model (GMM) for active tectonic regions and a local, Türkiye-specific model; and analyze spatial variability of several ground motion intensity measures (peak ground acceleration—PGA, peak ground velocity—PGV, 5%-damped pseudo-spectral acceleration—PSA—for 0.1, 0.3, 0.6, 1.0, 3.0, 6.0, and 10.0 s oscillator periods, Arias intensity—I _A , and cumulative absolute velocity—CAV), which are useful for ground motion estimation at sites without recordings.

Several previous ground motion compilations have been presented since the 6 February 2023 mainshock (Baltzopoulos et al., 2023; Garini and Gazeta, 2023; Gülerce et al., 2023; Sandıkkaya et al. 2024), and the work presented in this article updates the work presented in Chapter 3 (Buckreis et al., 2023b) of the Earthquake Engineering Research Institute (EERI)- Geotechnical Extreme Events Reconnaissance (GEER) Association report on this event (Cetin et al., 2023). This work is different from other previous efforts for several reasons, including the use of later-released data from the Disaster and Emergency Management Authority (AFAD) network that corrects errors from earlier releases, the application of Next-Generation Attenuation (NGA) protocols for data processing and metadata compilation, and the presentation of data for four events that were of primary interest to reconnaissance teams (6 February M7.81, 6 February M6.81, 6 February M7.74, and 20 February M6.37 events). The data presented here were compiled for use by reconnaissance teams (Cetin et al., 2023), researchers engaged in NGA projects, and other interested researchers. The data are publicly available at https://doi.org/10.17603/ds2-t115-bk16 (Buckreis et al., 2023a).

Ground motion data and metadata

Ground motion data for the February 2023 Türkiye earthquake sequence was recorded by seismic networks operating mainly in Türkiye and Cyprus. The networks operating in Syria and Lebanon report having recordings from these events (Salah Sadek, personal communication), but the data have not been released to date. The AFAD in Türkiye operates the Turkish National Strong Motion Network (TNSMN) and the Turkish National Seismic Network (TNSN). The TNSMN is a strong motion network, consisting mainly of accelerometers, while the TNSN sensors are mainly broadband seismometers. TNSMN and TNSN station locations are typically distinct, with only a few collocated stations. These networks produced the majority of recordings from the earthquake sequence, including all of the near-fault records shown in Figure 1. The principal additional network is the Kandilli Observatory and Earthquake Research Institute, which operate stations mainly in northern parts of Türkiye.

Figure 1.

Map of Türkiye depicting locations and focal mechanisms of M7.81, M6.81, M7.74, and M6.37 earthquakes, and seismic stations that recorded usable ground motion data. Surface projections of finite-fault representations for the M7.81, M7.74, and M6.37 earthquakes are represented by red, orange, and purple polygons, respectively. Sites with V _S30 assigned from measured V _S profiles or from proxy-based models are shown by green and blue triangles, respectively.

Ground motion data

Raw (unprocessed) waveforms for each event were obtained from the Earthquake Data Center System of Türkiye (TDVMS; https://tdvms.afad.gov.tr/) and the Incorporated Research Institutions for Seismology (IRIS; Trabant et al., 2012). Data were initially screened to remove duplicate records. The main content of the current data set was downloaded in late March 2023, with a few records added in February 2024; these downloads replaced records accessed in February 2023 that had been partially processed and had baseline correction problems. An initial visual inspection was performed to identify and remove records with non-usable data, which include noise-dominated, unclear event signal, early termination, late trigger, instrument malfunction, and spike records.

Many of the recordings from the M7.81 event exhibited early termination, likely a result of power failure during the strong shaking. Ground shaking exceeded the upper amplitude limit of most seismometers (instrument code “H”) during the M7.81 and M7.74 events, so recordings for these events are mostly from accelerometers (instrument code “N”). Unfortunately, there were no usable recordings near Gölbaşı, Şekeroba, and Çelikhan, and few near İskenderun and Kahramanmaraş for the M7.81 and M7.74 events. Station TK 3112 in İskenderun is located in an area with observed liquefaction; however, it did not produce usable recordings for the M7.81 and M7.74 events. Few usable records were obtained in the source regions for the M6.81 and M6.37 aftershocks (near Türkoğlu and Antakya, respectively). At the time of writing, a total of 337, 365, 284, and 229 usable three-component records at stations as far as 571, 582, 497, and 492 km were identified for the M7.81, M7.74, M6.81, and M6.37 events, respectively.

Near-fault rupture waveforms for the M7.81 event sometimes included multiple wave packets. We did not attempt to separate these waveforms and chose to use the entire ground motion time-series during processing and intensity measure computations.

Signal processing

Each of the three-component records was processed individually according to standard protocols developed for NGA projects (e.g. Goulet et al., 2021). The procedure consists of applying a baseline correction and high-pass (and sometimes low-pass) Butterworth filters in the frequency domain. Corner frequencies are selected by visual inspection of the Fourier amplitude spectra (FAS), signal-to-noise ratio (SNR), and reasonableness (lack of wobble) of displacement time-series. The lowest high-pass corner frequency that satisfies a minimum SNR threshold and produces a displacement time-series without wobble is selected.

This procedure removes any static displacement that might otherwise be present in near-field records, which is the case for records within approximately 40–70 km of the M7.81 and M7.74 events. Processing that includes preservation of the fling is not presented in this article. When using standard processing that does not allow for permanent displacement with records for which such displacements are present, unique challenges arise in the selection of high-pass corner frequencies, $f_{c,\mathrm{HP}}$ . The difficulty is well illustrated by the TK 4615 (north component) record in Figure 2, where the different panels show acceleration, velocity, and displacement. This station is located about 0.7 km from the M7.81 rupture.

Figure 2.

Example time-series plots illustrating sensitivity to the high-pass corner frequency selected for the north component of station TK 4615 from the M7.81 mainshock: (a) acceleration, (b) velocity, and (c) displacement.

In each plot, time-series are shown for the unfiltered record, and the record subjected to a relatively modest high-pass filter ( $f_{c,\mathrm{HP}}$  = 0.0159 Hz) that preserves much of the low-frequency energy and a more aggressive high-pass filter ( $f_{c,\mathrm{HP}}$  = 0.0531 Hz). These filter corners were selected considering different objectives: (1) the relatively modest filter ( $f_{c,\mathrm{HP}}$  = 0.0159 Hz) aimed to preserve the velocity pulse between 8 and 15 s, which is associated with downward (negative displacement) fling step and (2) the more aggressive filter ( $f_{c,\mathrm{HP}}$  = 0.0531 Hz) aimed to remove long-period displacement wobble that occurs when the pulse features are preserved.

As illustrated in Figure 2, there are different objectives that guide the selection of $f_{c,\mathrm{HP}}$ in the present context. The first approach has the advantage of preserving an important component of the ground motion (velocity pulse) but the disadvantage of artificial displacement wobble, while the second approach removes these features. It is not possible to preserve the velocity pulse and obtain reasonable displacements when using processing procedures that include a high-pass filter, as employed here. These considerations affect 48 and 12 ground motions from the M7.81 and M7.74 events, respectively, based on visual assessments during record processing (fling effects were not observed in recordings of the M6.37 and M6.81 aftershocks). When processing near-field fault records in the present work, we chose to prioritize preservation of the velocity pulse over the reasonableness of the displacements, which, for example, in Figure 2 produces the relatively modest value of $f_{c,\mathrm{HP}}$  = 0.0159 Hz.

The longest usable period of the ground motion is defined as T ≤ 1/(1.25 f_c,HP), where f_c,HP is the greater of the two horizontal high-pass corner frequencies selected during signal processing for combined horizontal components (e.g. RotD50 as defined in Boore, 2010). Figure 3 presents the number of usable combined horizontal components as a function of oscillator period for the four considered events. The M6.37 and M6.81 events have a more rapid rate of decay of usable records with period for T > 5 s than the M7.81 and M7.74 events, which results from the reduced long-period energy from these smaller magnitude events. The reduced long-period energy increased the likelihood that records will be noise-dominated at these long periods.

Figure 3.

(a) Distribution of records with usable data with respect to magnitude and rupture distance; (b) number of usable records versus period for each event. Longest usable period defined as T ≤ 1/(1.25 f_c,HP), where f_c,HP is the greater of the two horizontal high-pass corner frequencies.

The rate of decay for the M6.81 event also reflects challenges in selecting appropriate values for f_c,HP, because most records contained no pre-event noise due to the occurrence of an M5.7 aftershock 82 s prior to the M6.81. The relative spatial and temporal proximity of this smaller aftershock meant that many sites were still experiencing shaking when the waveforms from the M6.81 event arrived. During processing for this event, we assume the noise to be stationary and utilize noise windows prior to the M7.81 event (~10 min before the M6.81 event) to compute SNR, as described in Buckreis et al. (2023b). The displacement wobble check for recordings of this event was especially challenging due to the residual shaking from the M5.7 event; therefore, f_c,HP selections for this event erred on the side of conservatism.

Source and path metadata

Source metadata are assigned using NGA protocols (e.g. Goulet et al., 2021). The Global Centroid Moment Tensor (CMT) project (Ekström et al., 2012) origin times, seismic moments, and nodal plane solutions are preferred because they are derived using global recordings that average out radiation pattern variability and provide between-region consistency. However, CMT hypocenter locations correspond to the center of the earthquake moment distribution in time and space, which may not align well with the location of the initial slip. Therefore, preferred hypocenter locations are adopted from the USGS National Earthquake Information Center (NEIC). Table 1 summarizes the compiled moment tensor metadata for each event. The strike, dip, and rake angle provided in Table 1 are based on the solution plane that best aligns with the known orientation strike of the East Anatolian Fault.

Table 1.

Preferred source metadata for the Türkiye earthquake sequence

Origin Time (UTC)	Longitude (o)	Latitude (o)	Depth (km)	M	M0 (dyne-cm)	Strike (o)	Dip (o)	Rake (o)
2023-02-06 01:18:10	37.014	37.226	10	7.81	5.80E27	51	70	−4
2023-02-06 01:28:15	37.189	36.893	9.8	6.81	1.85E26	212	73	−13
2023-02-06 10:24:59	37.196	38.011	13.5	7.74	4.53E27	264	46	−9
2023-02-20 17:04:29	36.025	36.162	11.5	6.37	4.08E25	226	42	−16

Multiple finite-fault models (FFMs) for the M7.81, M7.74, and M6.37 events are available. We follow NGA criteria from Contreras et al. (2022) to evaluate alternate FFMs for each event, which led to the selection of the FFMs described by Xu et al. (2023) for the M7.81 event, Goldberg et al. (2023) for the M7.74 event, and the National Institute of Geophysics and Volcanology of Italy (INGV; 2023) for the M6.37 event. The FFMs for the two large events were selected because they were developed with consideration of different types of data (e.g. global teleseismic waveforms, local seismic data, and permanent crustal displacements), and produce rupture locations compatible with observed surface rupture (Kozaci et al., 2023). The Xu et al. (2023) model was applied without trimming, whereas the Goldberg et al. (2023) and INGV (2023) models were trimmed following standard NGA procedures. Simulations are used to develop a rupture surface for the M6.81 event, as outlined in Chiou and Youngs (2008) and Contreras et al. (2022). Table 2 summarizes the rupture surface representations for each event.

Table 2.

Summary of rupture surface representations for each event

M	Source	Segment	ULC Longitude (o)	ULC Latitude (o)	ULC Depth (km)	Length (km)	Width (km)	Strike (o)	Dip (o)
7.81	Xu et al. (2023)	1	36.9253	37.1198	0.0	63	21	33.5	90
		2	37.3195	37.5917	0.0	54	21	53	90
		3	36.8728	37.3949	0.0	45	21	61	90
		4	36.6052	37.0460	0.3	45	21	32	90
		5	36.1902	36.1448	0.6	108	21	20	90
		6	37.8136	37.8854	0.3	90	21	70	90
7.74	Goldberg et al. (2023) T	1	37.5969	37.9740	0.0	80.2	29.9	276	80
		2	36.6600	38.0390	0.0	25.1	23.0	250	80
		3	37.5753	37.9531	0.0	55.1	25.0	60	80
6.81	Simulated	1	37.0234	37.3113	0.7	35.8	15.0	212	73
6.37	INGV (2023) T	1	36.1689	36.2178	0.7	18	14	233	49

ULC = upper left corner.

T

Trimmed.

Source-to-site distances were evaluated using the finite rupture surfaces from Table 2 using the P4CF routine (Chiou, 2021). The P4CF routine was adapted for use with the geometrically complicated ruptures for the M7.81 and M7.74 events by computing distances for each segment separately and then selecting the minimum such distance. Recorded distance metrics include the closest distance from the site to any point on the rupture surface (R _RUP ), closest distance from the site to any point on the surface projection of the rupture surface (R _JB ), distance measured perpendicular to the fault strike from the surface projection of the top edge of the rupture surface (R _x ), distance measured parallel to the fault strike from the midpoint of the surface projection of the rupture surface (R _y ), the average distance to the fault plane (Chiou et al. 2000; Ohno et al. 1993) (R _avg ), the epicentral distance (R _epi ), and the hypocentral distance (R _hyp ).

Site metadata

Values of time-averaged shear wave velocity in the upper 30 m of the site (V _S30 ) were compiled for each distinct site. The AFAD website reports V _S30 values computed from V _S profiles measured using geophysical techniques (MASW and ReMi) for 229 of the stations, as shown in Figure 1. V _S30 values at 24 of those stations (TK 0125; TK 0201; TK 0603; TK 1201; TK 2401; TK 2518; TK 3113; TK 3116; TK 3133; TK 3143; TK 3144; TK 3301; TK 3801; TK 4614; TK 4628; TK 5001; TK 5201; TK 5505; TK 5508; TK 5801; TK 5804; TK 6004; TK 6302 and TK 6901) were updated relative to those on the AFAD website. These updates were motivated by some errors (specifically, inconsistencies with the 1D layered earth models at some stations). The corrections of these errors adjusted the V _S profiles and the V _S30 values.

Proxy-based V _S30 models were used to assign V _S30 values at locations without a measured V _S profile. For stations in Türkiye, two proxy relationships were used: (1) the Zhou (2023) proxy model based on topographic slope and (2) the Sahin et al. (2024) geology- and topography-based proxy models. In the case of the Zhou (2023) model, estimates are adjusted using kriging interpolation of V _S30 values from measured V _S profiles. The proxy relationship applied for each site is recorded in the station-site metadata file provided by Buckreis et al. (2023b). Outside of Türkiye, V _S30 values are assigned using the topographic slope-based proxy model proposed by Wald and Allen (2007).

Analysis of ground motion data

Comparisons of data to global and local GMMs

We compare recorded data to the following GMMs for shallow crustal earthquakes in tectonically active regions: (1) a global model—Boore et al. (2014; hereafter BSSA14) and (2) a Türkiye-specific model by Kale et al. (2015; hereafter KAAH15). The BSSA14 model can be applied with or without regional adjustments of the rate of anelastic attenuation, which is relatively low in this GMM for Türkiye and China (slower attenuation) and relatively high for Italy and Japan (faster attenuation). Without a regional adjustment, the model is considered applicable in other active crustal regions for which validation was performed, which are California, New Zealand, and Taiwan.

The GMMs provide predictions of ground shaking intensity given M, rupture mechanism/style-of-faulting (strike-slip for all events considered), R _JB , and V _S30 . BSSA14 has a model component that adjusts ground motions for the effects of V _S iso-surface depths (i.e. z_1.0); however, given that this site parameter is unavailable for the vast majority of sites, we choose not to apply any basin adjustments during this analysis (i.e. the differential depth δz_1.0 = 0).

Figure 4 shows the distance-dependence of RotD50 PGA, PGV, and PSA (1.0 s) for the M7.81, M7.74, M6.81, and M6.37 events. Recorded data are subdivided into three site groups: (1) rock (V _S30  ≥ 760 m/s), (2) stiff soil (360 ≤ V _S30  < 760 m/s), and (3) soft soil (V _S30  < 360 m/s). The median predictions shown for each of the GMMs are for a V _S30  = 360 m/s site condition. All GMMs fit the observed PGA data reasonably well out to distances of about 200 km. For PGV and PSA (1.0 s), the global and Italian regionalized versions of BSSA14 are nearly identical and fit the data better than the Türkiye-specific model (KAAH15) up to distances of approximately 100 km. At larger distances, the regional components of BSSA14 significantly affect predictions, with Italian and Türkiye versions attenuating faster and slower, respectively. At these larger distances, the global and Italian models fit the data better than the Türkiye model, which overpredicts. The KAAH15 GMM significantly underpredicts PGV in the near-field. Although the model includes an anelastic attenuation term in its functional form, the amount of anelastic attenuation is small, causing it to over-predict PGV and PSA (1.0 s) at far distances (100–200 km).

Figure 4.

Observations versus distance colored by site condition (rock, stiff soil, or soft soil) for (a) M7.81, (b) M7.74, (c) M6.81, and (d) M6.37 events. Predictions shown for global BSSA14, regional BSSA14, and KAAH15 GMMs for a V _S30 of 360 m/s.

The poor predictions of the BSSA14 model with the Türkiye regional adjustment may initially be surprising. However, the Turkish data considered in the development of that model were mainly from the northern part of Türkiye, near the North Anatolian Fault. This region corresponds to the Eurasian and Anatolian plates, whereas the present data are in the south of Türkiye near the interface of the Anatolian and Arabian plates (Figure 1). The results suggest that different crustal properties are encountered in the two parts of Türkiye, which is analyzed further in the “Spatial analysis” section of the article. For subsequent analyses, we use the global version of the BSS14 model, which predicts average rates of attenuation.

To better evaluate GMM performance relative to the data, we compute total residuals (R _ij ), considering the appropriate path and site conditions as follows:

$$R_{\mathrm{ij}}=\ln \left(Y_{\mathrm{ij}}\right)-\ln ({\mu }_{\mathrm{ij}})$$

(1)

where Y _ij is the value of the observed ground motion intensity measure (e.g. PGA) for event i from station j, and ${\mu }_{\mathrm{ij}}$ is the median GMM prediction. Figure 5 presents plots of R _ij versus distance for all GMMs. Binned means are also shown for R _ij , along with their 95% confidence intervals. Only data from free-field sites are shown; data from 18 instrumented dam sites and 2 borehole stations are screened out due to the potential for impactful soil–structure interaction effects.

Figure 5.

Total residuals (R _ij ) versus distance for (a) M7.81, (b) M7.74, (c) M6.81, and (d) M6.37 events. Residuals computed for PGA, PGV, and PSA (1.0 s) for the global BSSA14 and KAAH15 GMMs.

The residuals shown in Figure 5 reflect the combined impact of potential differences between the data and GMMs related to source, path, and site effects. Residual amplitudes at short distances are generally indicative of source effects, whereas the general trend of R _ij against R _JB reflects path effects. Site effects cannot be easily interpreted from plots of R _ij versus R _JB . The trends observed in Figure 5 support the initial observations drawn from Figure 4, namely the existence of complex path effects that result in relatively poor fits between the GMMs and data at large distances, generally R _JB  > 200 km (e.g. downward- or upward-trending residuals with distance suggesting higher or lower rates of attenuation, respectively). In general, BSSA14 performs the best at distances less than approximately 200 km, beyond which uncaptured path effects are apparent. KAAH15 performs comparably, except for a tendency to under-predict PGV over a wide distance range, and the other intensity measures at short-to-moderate distances.

Event terms ( ${\eta }_{E,i}$ ), which represent the event-specific bias, can be approximated as the average R _ij for each event. This approach approximates the value that would be obtained from mixed-effects analyses since the number of recordings per event is large. To avoid mapping potential bias from the uncaptured path effects into the estimation of ${\eta }_{E,i}$ , only observations within 200 km are used. Figure 6 presents plots of ${\eta }_{E,i}$ and their standard errors for PGA, PGV, and PSA at periods ranging from 0.01 s to 4 and 10 s for KAAH15 and BSSA14, respectively. These results show that the BSSA14 model underpredicts long-period motions (>0.5 s) and PGV for the three largest events (positive ${\eta }_{E,i}$ ). On the contrary, it generally captures shorter period motions and PGA with little bias. The KAAH15 model is relatively effective at long periods, but underpredicts at short-to-moderate periods (T < 0.2 s). For the M6.37 event, BSSA14 overpredicts all intensity measures whereas KAAH15 underpredicts at short-to-moderate periods (T < 0.5 s) and over-predicts to a similar degree as BSSA14 at longer periods.

Figure 6.

Event terms ( ${\eta }_{E,i}$ ) computed from stations within 200 km versus oscillator period for (a) M7.81, (b) M7.74, (c) M6.81, and (d) M6.37 events computed using global BSSA14 and KAAH15 GMMs.

Spatial analysis

Correlation models are developed to characterize the spatial pattern of within-event residuals ( $\delta W_{\mathrm{ij}}$ ) for each event from the BSSA14 (global) model. To this end, δW _ij at the station locations are computed as:

$$\delta W_{\mathrm{ij}}=R_{\mathrm{ij}}-{\eta }_{E,i}$$

(2)

The $\delta W_{\mathrm{ij}}$ values are then normalized by their standard deviation ${\varphi }_j$ as:

$$\delta {\stackrel{\prime }{W}}_{\mathrm{ij}}=\frac{\delta W_{\mathrm{ij}}}{{\varphi }_j}$$

(3)

The $\delta \stackrel{\prime }{W}ij$ values are used to develop the correlation models. Our approach for developing correlation models accounts for between-station closeness in Euclidean distance ( ${\rho }_E$ ) and azimuthal distance ( ${\rho }_A$ ), as proposed by Bodenmann et al. (2023):

$${\rho }_{\mathrm{EA}}(d_E,{\gamma }_E)={\rho }_E(d_E,{\gamma }_E)·{\rho }_A(d_A)$$

(4)

where:

$${\rho }_E(d_E,{\gamma }_E)=\exp [-1·{\left(\frac{d_E}{{\ell }_E}\right)}^{{\gamma }_E}]$$

(5)

$${\rho }_A\left(d_A\right)=(1+\frac{d_A}{{\ell }_A})·{(1-\frac{d_A}{180})}^{\frac{180}{{\ell }_A}}$$

(6)

and $d_E$ is the Euclidean distance (km) between two seismic stations, ${\gamma }_E$ is a model coefficient, ${\ell }_E$ is the Euclidean correlation length (km), $d_A$ is the azimuthal angle (degrees) between two seismic stations, and ${\ell }_A$ is the azimuthal correlation angle (degrees).

Bayesian inference is used to estimate the model coefficients ${\ell }_E$ , ${\gamma }_E$ , and ${\ell }_A$ , which allows for establishing prior distributions for the model parameters and updating them using the available ground motion data. Bayesian inference provides multiple models that in aggregate describe the model uncertainty. A set of 1000 model coefficients is sampled from the posterior joint distributions, resulting in 1000 correlation models. The maximum a posteriori (MAP) model (i.e. the model most supported by the data) is also estimated and used to develop δW _ij maps. Example correlation models for PGA for the M7.81 mainshock are presented in Figure 7.

Figure 7.

Bayesian correlation models for peak ground acceleration (PGA) within-event residuals from the M7.81 mainshock. (a) Correlation models dependent on the Euclidean distance ( ${\rho }_E$ ) and (b) correlation models dependent on the azimuthal angle ( ${\rho }_A$ ).

The correlation models are applied with ordinary Kriging to estimate the mean δW _ij at unsampled locations ( $\delta {\tilde{W}}_{\mathrm{ij}}$ ) and the associated variance (or standard deviation) of the interpolation error. The interpolation errors depend on the distances from a given uninstrumented location to the closest seismic stations. Locations close to a seismic station have low variance in the estimated ground motions, while locations further away from seismic stations have higher variance, which at the limit approach the within-event variance from the GMM. The variance of the interpolation error at seismic stations is zero (for more details about the ground motion estimation procedure, please refer to Pretell et al. (2024a)). The spatial correlation models and raster files with within-event residual across Türkiye for ground motion intensity measures of the four events are accessible at a data repository (Pretell et al., 2024b).

Figures 8 and 9 present the derived $\delta {\tilde{W}}_{\mathrm{ij}}$ maps that demonstrate spatial patterns of PGA, PGV, PSA (0.3 s), and PSA (1.0 s) residuals for the M7.81 and M7.74 events, respectively. The colorization in Figures 8 and 9 indicates areas of weaker- and stronger-than-average ground motions (blue and red, respectively). Interestingly, there are discontinuities in the blue and red shading at the East Anatolian Fault plate boundary, which are most acute for high-frequency intensity measures. We generally find negative residuals (blue) on the Anatolian plate (northwest of the fault) and positive residuals (red) on the Arabian plate (southeast of the fault). These results suggest that the global BSSA14 model overpredicts on the Anatolian plate and underpredicts on the Arabian plate. The results shown in Figure 8 are fairly consistent across the other events, as shown in Figure 9 for the M7.74 event and supporting figures provided in the Electronic Supplement. These findings suggest the presence of tectonic block-dependent anelastic attenuation in southern Türkiye, with faster attenuation (lower seismic quality factor, Q; Aki and Richards, 1980) on the Anatolian plate and slower attenuation (higher Q) on the Arabian plate, consistent with previous studies (e.g. Aydin, 2020; Aydin and Şahin, 2011).

Figure 8.

Mean within-event residual for (a) PGA, (b) PGV, (c) PSA (0.3 s), and (d) PSA (1.0 s) from the M7.81 mainshock developed using the maximum a posteriori correlation models. Blue and red shading indicate areas of weaker- and stronger-than-average ground motions, respectively. Non-shaded areas are beyond the usable distance range of the BSSA14 GMM.

Figure 9.

Mean within-event residual for (a) PGA, (b) PGV, (c) PSA (0.3 s), and (d) PSA (1.0 s) from the M7.74 event developed using the maximum a posteriori correlation models. Blue and red shading indicate areas of weaker- and stronger-than-average ground motions, respectively. Non-shaded areas are beyond the usable distance range of the BSSA14 GMM.

Ground motion estimation procedure

Mean ground motion estimates ( $\hat{y}$ ) can be obtained at sites without instrumentation using GMM median predictions, event terms, and within-event residuals maps as follows:

$$\ln ({\widehat{y}}_{ij})=\ln ({\mu }_{ij})+{\eta }_{E,i}+\delta {\tilde{W}}_{ij}$$

(7)

where the computation of ${\mu }_{\mathrm{ij}}$ requires estimation of V _S30 and R _JB for the site (in the case of BSSA14). Maps of $\delta {\tilde{W}}_{\mathrm{ij}}$ using the BSSA14 model, with ${\eta }_{E,i}$ from Figure 6, are presented in Figures 8 and 9 and in the Electronic Supplement. Uncertainties in the ground motions are also provided in the Electronic Supplement. For any given location, these uncertainties result from the aleatory variability in the Kriged within-event residual and epistemic uncertainty in $\delta {\tilde{W}}_{\mathrm{ij}}$ that is related to the multiple possible semi-variogram models. This procedure has been used to estimate RotD50 PGA, PGV, and 5%-damped PSA at 0.1, 0.3, 0.6, 1.0, 3.0, 6.0, and 10.0 s, I _A , and CAV for sites of interest to reconnaissance and other research teams (Cetin et al., 2023; Moss et al., 2024; Moug et al., 2024; Numanoglu et al., 2024personal communication Hariri Mohammadamin, NIST, September 2023).

Summary and conclusions

This article presents the ground motion data and spatial correlation models for four events of the February 2023 Türkiye earthquake sequence. The processed ground motion data and metadata are accessible at Buckreis et al. (2023b), and correlation models are available at Pretell et al. (2024b). This data will be influential for future GMM development efforts and used as validation data for ground motion simulations, among other applications. Comparisons of ground motion data to global and local GMMs suggest the existence of event-specific effects and complex path effects in the region, uncaptured by the existing models. Spatial analysis of the within-event residuals suggests that the complex attenuation features are likely related to differences in crustal properties between the Anatolian and Arabian plates.

Supplemental Material

sj-pdf-1-eqs-10.1177_87552930241259024 – Supplemental material for Engineering attributes of ground motions from February 2023 Türkiye earthquake sequence

Supplemental material, sj-pdf-1-eqs-10.1177_87552930241259024 for Engineering attributes of ground motions from February 2023 Türkiye earthquake sequence by Tristan E Buckreis, Renmin Pretell, M Abdullah Sandikkaya, Özkan Kale, Aysegul Askan, Scott J Brandenberg and Jonathan P Stewart in Earthquake Spectra