Coseismic drop of seismic velocity caused by the 2023 Turkey–Syria earthquakes

SUMMARY

The M_w 7.8 earthquake in Turkey on 6 February 2023 was extraordinary for various reasons. It originated in depth of only 10 km, ruptured along a fault plane around 300 km long and the surface was covered by an extensive network of high-quality seismic instruments. The strong motions resulted in a vast number of tragic casualties and huge material losses in Turkey and Syria. However, abundant and proximate seismic observations of this event and numerous aftershocks give an opportunity to deepen the understanding of earthquake processes. In this study, we carried out an assessment of coseismic changes of seismic velocity using passive image interferometry. We used data from one strong-motion and twenty-four broad-band sensors. We observed coseismic drops of seismic velocity, which reached up to −1.79 per cent at a location directly at the ruptured East Anatolian Fault Zone. Along the M_w 7.8 earthquake fault, we observe frequency dependence of the velocity changes. At frequencies above 0.5 Hz, the velocity drops seem to be higher at locations close to the ruptured faults than in the more distant areas.

1 INTRODUCTION

On 6 February 2023, two major earthquakes struck southern central Turkey, northwestern Syria and adjacent regions. The first earthquake with M_w 7.8 occurred at 01:17 UTC on the East Anatolian Fault (EAF, see Fig. 1). This major strike-slip fault zone separates the Arabian and Anatolian plates (DeMets et al. 1990). The fault is a 700 km long, complex tectonic structure with predominant left-lateral motion and abundant seismicity (e.g. Güvercin et al. 2022). The double couple source mechanism for the M_w 7.8 event, defined from seismic moment tensor with zero trace according to global CMT-catalogue (Dziewonski et al. 1981; Ekström et al. 2012), fits well to the southwest–northeast orientation of the EAF (Fig. 1). The second earthquake with M_w 7.5 occurred 9 hr later on another documented active fault in west–east direction (Styron & Pagani 2020). The aftershocks within 1 week after the main shocks (EMSC/CSEM, emsc-csem.org) give a rough estimate of the ruptured fault areas (Fig. 1).

In order to provide a rapid assessment of the material changes related to these events, we used a well-proven method named passive image interferometry (PII) that was introduced and first tested by Sens-Schönfelder & Wegler (2006). The authors took the advantage of possibility to approximate the Green’s function (GF) by correlating recordings of a diffuse wavefield propagating through the medium (Lobkis & Weaver 2001). In the real-world subsurface, the diffusive nature of the wavefield is ensured by omnipresent scatterers (Sato et al. 2012) and noise sources. Once the GF is repeatedly estimated over time, an assessment of relative changes of seismic velocities (dv/v) against a reference GF (e.g. an average over all the functions) can be done by stretching individual time samples against the reference (dv/v = −dt/t, see Snieder et al. 2002). This methodology was previously applied in various configurations to study coseismic velocity changes, for example in California (Brenguier et al. 2008), in Japan (Wegler & Sens-Schönfelder 2007; Wegler et al. 2009; Hobiger et al. 2012, 2014, 2016; Viens et al. 2018) and in Chile (Richter et al. 2014). Estimation of coseismic dv/v in the works of Hobiger et al. (2014, 2016) involved cross correlations of different components within individual stations. We customized their approach and used it in this study.

2 DATA AND PROCESSING

The data consisted of continuous, three-component-recordings from one strong-motion (SM) station KHMN and twenty-four broad-band (BB) stations (Fig. 1 and Table A1). Each of the stations belongs either to the KOERI network (Kandilli Observatory and Earthquake Research Institute, Boğaziçi University 1971) or Turkish National Seismic Network (Disaster and Emergency Management Authority 1990). They are abbreviated as KO and TU, respectively. Double proof was made in order to assert that the stations did not topple during the main shocks and still record reliably: (1) Rayleigh waves from the M_w 5.6 earthquake in Romania (14.2.2023) were inspected. They evidenced that the radial component preceded the vertical one by 90°; and (2) Teleseismic P-wave onsets after the M_w 6.1 event in the Philippines (15.2.2023) at the stations we used showed coherence on Z-components with the adjacent stations.

The time span of the selected data was from 1.12.2022 to 14.2.2023. There are almost no remarkable data gaps in that period, with an exception at stations MKAM, NARI and SARI (1–3 weeks long). Processing involved the YAM package (Eulenfeld 2023). Data preparation for PII included frequency filtration (0.125–0.25, 0.25–0.5, 0.5–1, 1–2 and 2–4 Hz), amplitude clipping (amplitudes exceeding daily RMS amplitude more than four times were zeroed), 1-bit normalization and spectral whitening. In addition, the SM data from KHMN were corrected for instrument response since its response to velocity is not flat. Then, individually at each station, hourly cross correlations of all three component pairs were computed (intervals with lower than 50 per cent availability were discarded). For each station, correlation functions for each pair were stacked for the time periods: (i) Prior to 6 February, 1 AM (more than 2 months) and (ii) After 6 February (8 d). This means that 23 hr on 6.2.2023 were omitted because of data gaps and deteriorated quality of the correlation functions at the day of the M_w 7.8 and 7.5 events. Pre- and post-earthquake correlation stacks from recordings at YAYL, AKCD and KHMN are shown in Fig. 2. Afterwards, the two stacks (i) and (ii) were stretched against each other in frequency-dependent coda lapse time windows (from 48–128 s at 0.125–0.25 Hz to 3–8 s at 2–4 Hz, as enumerated in the header of Table 1). The start times and lengths of the coda windows were linearly bound to the corresponding frequency bands via periods at the lowest frequency in the corresponding frequency range. Each time window started six periods after the zero lag time of the correlation function. And the time windows were always 10 periods long. This is similar to the works of Hobiger et al. (2012, 2014, 2016), in which start time of 7.5 periods after the zero lag time and length of 10 periods were successfully applied to inter- and single-station cross correlations and also to autocorrelations. After the stretching, the resulting three 1-D matrices with correlation coefficients for east–north, east–vertical, and north–vertical component pairs at each station were pointwise-averaged. Finally, the position of maximum correlation coefficient in the mean 1-D matrix (mean CC) yielded the coseismic dv/v estimate (Table 1). This was done for each of the five frequency bands. Results with mean CC lower than 0.55 were not taken into account. Note, that the results presented later in Figs 3 and 4 involve 24 hr stacks and the daily dv/v estimates came from stretching against a mean trace of the whole time span. Thus, the examples in the two figures only serve for a visualization of the data and correlation quality, while the actual coseismic dv/v estimates resulted from the two, much more solid stacks, acquired from the time periods before and after the earthquakes on 6.2.2023 shown in Table 1.

3 RESULTS

First, examples of the results are presented in Figs 3 and 4. Fig. 3 shows similarity matrices from data at SM station KHMN (2-4 Hz) with significant coseismic dv/v change of almost −2 per cent. Fig. 4 gives four examples of daily-stacked dv/v time-series, both for individual component pairs, as well as for their means. Secondly, the final and more robust estimates (as described in Section 2) are presented in Table 1 and Figs 5 and 6. The main points resulting from our analysis are summarized below.

1. Table 1 contains the actual values of coseismic velocity changes (Section 2) with their standard deviations (σ), which were calculated from three dv/v values resulting from individual component pair estimates (east–north, east–vertical and north–vertical). Fig. 5 shows maps with the velocity changes listed in Table 1. Individual subplots represent different frequency bands, with the exception of the 1–2 Hz band, which yielded no valid estimate. One of the dv/v values indicated a net velocity increase (i.e. dv/v is positive even after subtraction of the standard deviation). Namely, station CUGUR in the 0.125–0.25 Hz band with mean CC = 0.57 and dv/v = 0.19 ± 0.05 per cent. However, this is a single occurrence and the station is one of the most distant from any of the fault zones associated with the earthquake doublet. Moreover, dv/v at CUGUR at different frequency and two other estimates at nearby station DARE indicate either zero or negative dv/v of up to 0.3 per cent (including σ). Thus, we did not interpret this observation any further. On the other hand, the maximum coseismic dv/v of −1.79 ± 0.04 per cent is evident at station KHMN (2–4 Hz). This station is marked with a square in Figs 1 and 5d. This maximum drop is accompanied by the highest mean CC among the stations involved in our processing and very low σ. This SM station is located directly at the fault zone of the M_w 7.8 main shock.

2. We had a closer look at the velocity changes with mean CCs ≥ 0.55 acquired in areas at the fault zone of the M_w 7.8 main shock. These stations are shown as inverted triangles (broad-band seismometer) and as a square (strong-motion seismometer) in Figs 1 and 5. The changes are plotted against their corresponding frequency bands in Fig. 6. They indicate that the velocity drops are frequency-dependent (higher drops at higher frequencies). This seems to be true for both the overall trend in that figure, as well as for individual stations KAMA, MGND and YAYL, at which reasonable estimates in multiple frequency bands were retrieved.

3. It is apparent from Table 1 that the BB stations yielded higher mean CCs at frequencies lower than 1 Hz, while there is no observation with mean CC ≥ 0.55 in the 1–2 Hz band and only two in the 2–4 Hz band. Contrarily, the SM station KHMN yielded no sufficient mean CC value at frequencies lower than 2 Hz and the only mean CC ≥ 0.55 was observed in the 2–4 Hz band (0.87). However, it is uneasy to decide whether this difference between the BB sensors and the only SM sensor is due to the different transfer functions, electronic noise, site constructions, bedrock properties, noise sources or anything else.

4. Recordings from the SM station KHMN enabled to estimate peak ground velocity and acceleration (PGV and PGA) during the very proximate M_w 7.8 main shock and relate them to the observed velocity drop. It was not possible to include observations from the M_w 7.5 aftershock because of a data gap during that period. Our PGV and PGA estimates are based on envelopes of three-component recordings. PGV reached 1.29 ms⁻¹ and PGA was 6.53 ms⁻². Local, slope-based V_S30 estimate is 270 ms⁻¹ (Allen & Wald 2009, accessible at https://earthquake.usgs.gov/data/vs30/). V_S30 is average shear wave velocity in the shallowest 30 m of the subsurface. These values enable to calculate maximum dynamic strain |$\gamma _{d\ {\rm max}} = PGV / V_{S30}$| = 4.78 × 10⁻³ and relate the PGA to the coseismic dv/v of −1.79 per cent observed at 2–4 Hz. Both of them are generally used for characterization of site response to propagating seismic waves (e.g. Hobiger et al. 2014, 2016; Richter et al. 2014; Viens et al. 2018). Interestingly, Brenguier et al. (2014) observed volcanic fluids-related susceptibility of dv/v to PGV in Japan. Since we have a single data point, detailed comparison with previous studies is not possible in this case. Just to give an example, our data indicate considerably higher sensitivity of dv/v to PGA and to the maximum dynamic strain γ_d max than in Hobiger et al. (2014, 2016). However, the results they presented were acquired at lower frequencies (up to 2 Hz) and the coseismic velocity changes relied on borehole seismometers (they calculated PGA from recordings of vertically collocated, strong-motion sensors at the surface).

4 DISCUSSION

1. The coseismic velocity drop of 1.79 per cent at station KHMN (2-4 Hz) is one of the highest-known among the works that estimated the dv/v from scattered coda waves of ambient noise correlation functions (Wegler & Sens-Schönfelder 2007; Brenguier et al. 2008; Wegler et al. 2009; Hobiger et al. 2012, 2014, 2016; Richter et al. 2014; Viens et al. 2018). Viens et al. (2018) used single-station cross correlations of SM recordings to estimate coseismic velocity changes in the greater Tokyo area after the 2011 Tohoku-Oki earthquake in Japan. They used correlation functions filtered between 1–20 Hz and stretched them in a coda time window of 0–1 s. Such high frequencies and very early time window monitored shallow subsurface (upper ca. 100 m, as stated by the authors). With this setting, the largest-observed dv/v drop was around 11 per cent. Takagi et al. (2012) evidenced a drop of shear wave velocities of about 5 per cent due to the 2008 Iwate-Miyagi Nairiku earthquake, but this estimate resulted from direct parts or earthquake-related cross correlation functions between surface and borehole receivers, located up to 260 m apart (vertically). This also yielded shallower and more local estimates. Also, none of the cited works reported on credible coseismic velocity increase, which is consistent with our findings.

2. Our selection of frequency bands, single-station cross correlation configuration and component-pair-averaging was the same as in Hobiger et al. (2014, 2016). Regarding the temporal stacking, their approach was frequency-dependent (30–3 d with increasing frequency) and the dv/v estimate followed by curve-fitting of the final time-series. Our assessment of this coseismic change relied on direct stretching of pre- and post-earthquake temporal correlation stacks (ca. 2 months versus 8 d) and averaging of the resulting vectors of CCs. Thus, regarding the immediate coseismic changes, we stacked more prior to the earthquakes. During the post-earthquake period, we stacked fewer at frequencies below 1 Hz, while comparably and more in the higher bands. Although our correlation coefficients are generally lower, they still yielded reasonable estimates that are in the same orders of magnitude as in the other studies. Hobiger et al. (2012, 2014, 2016); Richter et al. (2014) reported on relation of coseismic d v/v with frequency (including depth extent). Namely, higher drops at higher frequencies. We observed this trend at locations in the immediate vicinity of the fault zone of the M_w 7.8 main shock, where we obtained sufficient amount of plausible dv/v estimates for such assessment (Fig. 6). This observation could be related to depth sensitivity of the dv/v, with higher frequencies imaging shallower depths. Namely, we used ambient noise wavefield that is mostly composed of surface waves whose penetration depth and depth sensitivity are frequency-dependent. However, Hobiger et al. (2012) found out that there is a trade-off between thickness of the superficial layer and the velocity drop in this layer, which does not allow for decisive conclusions without observations of dv/v in a broad frequency range and precise knowledge of the local velocity model. Also, Hobiger et al. (2016) observed the largest dv/v drops in locations close to the fault zones. This agrees well with our results at frequencies of 0.5–1 and 2–4 Hz (Figs 5c and d). However, we only obtained 5 and 3 values at that frequencies, respectively. We did not consistently observe this distance-related trend at frequencies below 0.5 Hz (Figs 5a and b).

3. We also observed a substantial difference between the GFs estimated prior to and after the earthquakes on 6.2.2023. Figs 2, 3 and 4 give a number of examples (mean CC colourscale in Fig. 4 is located on the joint dv/v curve). We suggest that the main reason is substantial bedrock damage along the fault zones. Such process also induces spatial perturbation of scatteres that affect the estimated GFs (e.g. Snieder et al. 2002). Moreover, cessation of traffic and industrial activity in heavily affected areas might explain the generally lower mean CCs at frequencies above 1 Hz (Table 1) that are usually dominated by the cultural noise wavefield. Temporal noise instability is one of the factors that unfavourably affect the GF estimate. Also, frequent aftershocks lead to a different background wavefield in comparison to the noise before the two main shocks, although this effect was reduced by the amplitude clipping (zeroing) and 1-bit normalization.

5 CONCLUSION

This report focused on imminent material changes caused by the 2023 Turkey–Syria earthquakes. They were examined with single-station cross correlations of continuous ambient noise recordings in five frequency bands ranging from 0.125 to 4 Hz. Below 1 Hz, we consistently observed frequency-dependent coseismic velocity drops in areas close to the main fault zone. The maximum change reached −1.79 per cent dv/v at 2–4 Hz in the immediate vicinity of the M_w 7.8 fault zone. This points to significant material damage in the subsurface. The results of this analysis brought promising outlook for future studies of coseismic rock damage and subsequent relaxation related to these major and well-recorded earthquakes.

ACKNOWLEDGEMENTS

We would like to thank two anonymous reviewers for their valuable remarks and suggestions during the review process. This work was cofunded by the project nr. 451271527 of the German Research Foundation (DFG).

DATA AVAILABILITY

Moment tensors are from CMT-catalogue (Dziewonski et al. 1981; Ekström et al. 2012), tectonic lines from Nuvel-1 model (DeMets et al. 1990) and Global Active Faults Database (Styron & Pagani 2020) and earthquake relocations from EMSC/CSEM (emsc-csem.org). Waveform data were provided by KOERI (Kandilli Observatory and Earthquake Research Institute, Boğaziçi University 1971) and Turkish National Seismic Network (Disaster and Emergency Management Authority 1990). Seismograms recorded at stations ANDN, DARE, GAZ, KHMN and SARI were downloaded via services provided by the International Federation of Digital Seismograph Networks, while all the other recordings were downloaded from tdvms.afad.gov.tr/continuous_data.

References

APPENDIX: STATION COORDINATES