Sub-and super-shear ruptures during the 2023 M w 7.8 and M w 7.6 earthquake doublet in SE Türkiye

Abstract

Abstract An earthquake doublet (Mw 7.8 and Mw 7.6) occurred on the East Anatolian Fault Zone on February 6 th , 2023. The events produced significant ground motions, in excess of 150%g, and caused major impacts to life and infrastructure throughout SE Türkiye and NW Syria. Here we show the results of earthquake relocations of the first 11 days of aftershocks and rupture models for both events inferred from the joint kinematic inversion of HR-GNSS and strong motion data considering a multi-fault and 3D rupture geometry. We find that the first event nucleated on a previously unmapped fault before transitioning to the East Anatolian Fault (EAF) rupturing for~350 km and that the second event ruptured the Sürgü fault for~160 km. Maximum rupture speeds were estimated to be 3.2 km/s for the Mw 7.8 event. For the Mw 7.6 earthquake, we find super-shear rupture at 4.8 km/s westward but sub-shear eastward rupture at 2.8 km/s. Maximum slip for both events were as large as~8 m and~6 m, respectively. Özet (Turkish) 6 Şubat 2023 tarihinde Doğu Anadolu Fay Zonu'nda Mw 7.8 ve Mw 7.6 büyüklüklerinde bir deprem çifti meydana geldi. Depremlerin ürettiği kuvvetli yer hareketleri,1.5g değerini aşarak, güneydoğu Türkiye ve kuzeybatı Suriye'de yaşam ve altyapı üzerinde önemli etkilere ve yıkımlara neden oldu. Bu çalışmada, deprem sonrasında ilk 11 günlük artçı depremlerin yeniden lokasyonlarını ve deprem çifti için çoklu fay ve bütünleşik üç boyutlu (3B) fay geometrisi kullanarak, HR-GNSS ve kuvvetli yer hareket verilerinin (SGM) birlikte analiz edildiği kinematik ters çözüm sonuçlarını göstermekteyiz. Anaşok, Doğu Anadolu Fayı (DAF) ile doğrudan ilişkilendirilmemiş ve daha önce haritalanmamış fay parçası üzerinde gelişerek yaklaşık 350 km uzunluğunda ve ikinci deprem Sürgü Fayı üzerinde yaklaşık 160 km boyunca kırılmıştır. Mw 7.8 depremi için maksimum yırtılma hızı 3.2 km/s olarak hesaplanmıştır. Mw 7.6 depremi için batıya doğru 4.8 km/s yüksek ve doğuya doğru 2.8 km/s düşük yırtılma hızlarını bulduk. Deprem çifti için maksimum yer değiştirme miktarı sırasıyla 8m ve 6m olarak hesaplanmıştır. Non-technical summary Two very large earthquakes occurred in south-eastern Türkiye on February 6 th 2023. In this paper we calculated kinematic models of how much the faults moved during both events and found very large displacements of as much as 6-8 m. We further calculated how fast the faults broke and found a "normal" behavior for the magnitude 7.8 earthquake. Meanwhile, the magnitude 7.6 broke extremely quickly in one direction (west) but at normal speed in the other direction (east). This fact is scientifically interesting and important to explain why ground shaking was so strong in the region.

Overview of the events
On February 6 th , 2023 at 01:17:35 UTC the M w 7.8 Nurdağı-Pazarcık earthquake nucleated~15 km southeast of the mapped trace of the East Anatolian Fault Zone (EAFZ, Figure 1A). Relocations (Figures 1B,1C) place the hypocenter at (37.0234°E, 37.2444°N, depth=12 km) and analyses of teleseismic data show a left-lateral source mechanism on a vertical or near vertical fault. A vigorous aftershock sequence followed and a little over 9 hours after the first event, at 10:24:49 UTC, the M w 7.6 Ekinözü earthquake occurred with a hypocenter at (37.2756°E, 38.0900°N, depth=15 km). It locates close to the mapped trace of the Sürgü fault (SF), and, as the event is of large magnitude and on a separate structure, we consider it as part of a "doublet" rather than a traditional mainshock/aftershock sequence (see Taymaz et al., 2022).
Ground motions recorded by a dense network of strong motion stations and inferred from the ShakeMap product showed intensities as high as MMI 8 or 9 for both events (U.S. Geological Survey, 2023a,b). At the time of writing this article, reports in the news media indicate at least 55,000 fatalities and over 5 million displaced people in Türkiye and Syria. The two earthquakes represent the largest in the EAFZ system and produced the largest ground motions in instrumental times including widespread liquefaction phenomena (Taftsoglou et al., 2023). They have been catastrophic for the entire region.
The EAFZ is one of the most seismically active areas in Türkiye and the Middle East. Its tectonics are complex and are still being studied to fully understand the geologic history of the region. The EAFZ is part of a major fault zone that runs through eastern Türkiye as it accommodates the tectonic movement between the Arabian and Anatolian microplates (Ambraseys, 1989). This shear deformation zone is represented by a 580-km long plate boundary and is associated with frequent shallow seismicity in the top~20-25 km of the crust (Taymaz et al., 1991;Tan and Taymaz, 2006;Taymaz et al., 2021;Melgar et al., 2020b). Relative plate motion is accommodated primarily by leftlateral strike-slip faulting at slip rates of 10±1 mm/yr (Reilinger et al., 2006) and has caused a series of destructive earthquakes in eastern Türkiye and northwest Syria as documented by historical records (Ambraseys and Jackson, 1998;Taymaz et al., 1991;Tan and Taymaz, 2006). Recent geological and geomorphic data indicate that the EAFZ has displaced the Euphrates River by 12 km since the mid-Quaternary (Trifonov et al., 2018) thus attaining a mean geological slip rate of 12-15 mm/yr. Yet, despite the dramatic effects of this fault's activity, a lack of high-resolution geodetic displacement data (e.g., achievable with continuous high-rate GNSS observations, HR-GNSS) has limited the capacity of constraining fault segmentation patterns, slip rate variations, earthquake recurrence intervals, and rupture dynamics.  Within this context, the earthquake doublet is of keen scientific interest for the region and for the study of large strike-slip systems generally. Here we will present the results of aftershock relocations and of kinematic slip inversions on a multi-fault 3D geometry using HR-GNSS and strong motion data. We will show that, for the M w 7.8 event, the kinematics are complex -it nucleates on a previously unmapped structure and propagates to the EAF which then triggers and slips bilaterally with a maximum rupture speed of 3.2 km/s. Likewise the M w 7.6 event ruptures bilaterally on the curved, roughly E-W striking Sürgü fault (SF) at super-shear speeds westward, likely as high as~4.8 km/s, but sub-shear eastward at 2.8 km/s. The slip is then partitioned between a splay parallel to the EAF and the continuation of the SF to the intersection with the EAF.

Available Data and Methods
We used regional geodetic and seismological data to produce an aftershock catalog and slip model as follows.

Double-Difference Hypocenter Relocations
We relocate a total of 5077 earthquakes, including, the mainshocks of the doublet, and 9 large aftershocks with magnitudes between M w 5.5 to 6.6. The phase data for this were acquired from the Disaster and Emergency Management Presidency of Türkiye (AFAD). It includes Pand Sarrivals from available stations selected by an automatized earthquake detection based on LTA/STA algorithm and initial locations estimated by the Hypoinverse algorithm (Klein, 2014). Most of the hypocentral depth estimates for these auto-located earthquakes range from 6.9 km to 7.1 km, i.e., more than 60% percent of aftershocks in this limited catalog.
To improve on this, we applied a relative earthquake location algorithm, hypoDD (Waldhauser and Ellsworth, 2000) using absolute Pand S-wave travel-time phase readings published in the AFAD bulletin. The algorithm makes use of earthquake pairs with very small hypocentral differences compared to event to inter-station distances. This allows direct association of the spatial offsets between the pairs to time delays between two events observed at a single station (Waldhauser and Ellsworth, 2000). hypoDD minimizes the difference between observed and calculated travel time residuals using relative hypocenter locations and origin times for all observed event-station pairs in an iterative manner. This approach overcomes potential bias originating from insufficient knowledge of structural complexities (e.g., velocity heterogeneities) along the sourcereceiver path, and, in this way, provides high-resolution hypocenter locations.
Travel-time differences are estimated for event pairs with less than 10 km of interevent distances and with a minimum of 8 connections between stations to define up to 10 neighbors at all 177 stations located within 200 km distance from the center of cluster. Initially 4756 out of 5077 aftershocks within the first 11 days were located following the Nurdağı-Pazarcık and Ekinözü earthquake doublet. Relative locations and origin times (OT) were obtained by a single set of 15 iterations in which large residuals were underestimated to suppress potential bias in the solution. We employed a 1-D initial velocity model that was updated through the relocation process of the 24 January 2020 M w 6.7 Doğanyol-Sivrice earthquake and its aftershock sequence (Taymaz et al., 2021;Melgar et al., 2020b). Our final database (see Data and Code Availability) consists of 2909 relocated events that had the highest resolution solutions ( Figure 1A).

Source inversion
Given the complexity of the rupture process of both events, defining the 3D geometry for inversion (Figures 1, 2) is critical for the success of the modeling. We combined information from the aftershocks, mapped traces of all known structures (EAFZ, and SF) and mapped surface ruptures from remote sensing (Reitman et al., 2023) to decide on the geometry as follows. We inferred there is a structure, which we hence call the Nurdağı-Pazarcık Fault (NPF) offset from the main strand of the EAF. We used the general trend of the aftershocks and a small, mapped surface rupture from remote sensing data (Reitman et al., 2023) to define its strike. As we will discuss later, this fault is necessary to fit the data. Further, the large (~15 km) offset between the hypocenter and the trace of the EAF provides additional support for its existence. For the junction of the SF with the EAF we used the mapped trace which connects the two faults. We also extended the SF into a small splay parallel to the EAF which is clearly visible in mapped surface offsets (Reitman et al., 2023). We used a vertical dip for the EAF southwest of the junction with the SF and used a vertical dip for the NPF as well. For the EAFZ northwest of the junction with the SF, and for the SF itself, we used a northward trending dip of 80°. This is supported by observations by Taymaz et al. (1991), Melgar et al. (2020b), and Taymaz et al. (2021) that reports a northward dipping geometry during the M w 6.7 Doğanyol-Sivrice earthquake in the segment of the EAF immediately northeast of where rupture for the M w 7.8 arrests. Additionally, the aftershocks are offset from the mapped surface traces and suggest a gentle northward deviation from vertical. We extended these geometries to a seismogenic depth of 20 km; this is supported by general observations of seismicity in the region from (Türkelli et al., 2003) and from the aftershocks ( Figure 1B, C). The 3D surface is meshed into triangles of mean vertex length of~5 km, resulting in 482 subfault elements and 256 subfault elements for the M w 7.8 and M w 7.6, respectively.
Next, we processed the geodetic and geophysical data as follows. HR-GNSS solutions were calculated at 1 Hz sampling rate using the precise point positioning method (PPP) as implemented in GipsyX (Bertiger et al., 2020). We used Jet Propulsion Laboratory rapid clocks and orbits (Noll, 2010) and rotated the solutions from geodetic coordinates to topocentric north, east, and up (vertical) coordinates. The displacement waveforms were low-pass filtered to 0.4 Hz prior to inversion. Likewise, the strong motion data were processed by remov- ing the instrument gain, removing the pre-event mean, and integrating to velocity for the M w 7.6 and to displacement for the M w 7.8. They were then bandpass filtered between 0.05 and 0.4 Hz, a total of 60 waveforms extracted from 12 three-component GNSS and 8 three-component strong motion sites contributed to the source inversion of the M w 7.8 event. For the M w 7.6 we used 10 three-component HR-GNSS and 5 threecomponent strong motion stations for a total of 45 waveforms. Locations of the stations are in Figure 1 and the station codes in Figure S1.
For the kinematic inversion, we employed the opensource MudPy code (Melgar and Bock, 2015), which implements the linearized multi-time window method. Elastodynamic Green's functions for both data sets were computed using the frequency-wavenumber integration approach of Zhu and Rivera (2002) with the sum of point sources used to represent each subfaults finite extent (see Koch et al., 2019). We assumed the 1-D layered model of Taymaz et al. (2021), which is appropriate for the region. The Green's functions (GFs) were filtered in the same passbands as the data before inversion. Rupture is allowed to nucleate at the hypocenter and a maximum rupture speed, v r max , is imposed. Note that in a multi-time window inversion this rupture speed is the upper bound allowed, slower rupture speeds are possible with subsequent time windows. We tested for both ruptures a range of values from 2.4 to 3.8 km/s for the M w 7.8 and 2.0 to 6.0 km/s for the M w 7.6 ( Figure 2C).
Each subfault is allowed to slip on one of five triangular slip rate functions. Each subfault has a fixed rise time: we used 5 s and 3 s for the M w 7.8 and M w 7.6, respectively. These values are obtained from the measurements of average rise times by Melgar and Hayes (2017) for large events worldwide. Each window has an overlap of 50% with the previous one, such that in total, at any given subfault, slip is possible for as long as 15 s for the M w 7.8 and 10 s for the M w 7.6. A nonnegative least squares solver is used, and we restrict the rake vector for all subfaults to a 90°window between -45°and 45°. The inversion is stabilized using Tikhonov regularization; no smoothness constraint (e.g., such as a Laplacian) is imposed. The regularization parameter is chosen using the L-curve criterion. Each of the two types of data are weighted according to their individual L2 norms as explained in Melgar et al. (2020b) and the vertical component of the HR-GNSS is down weighted by a factor of 3 to account for its higher noise levels (e.g., Melgar et al., 2020a).

Mainshock hypocenters and aftershock relocations
A careful inspection of time sequence of aftershock activity reveals three large aftershocks ranging from M w 5.6 to 6.6 that occurred within 18 minutes of the first main-shock with locations respectively southwest of it and a M w 5.6 to northeast 46 minutes later. The second main-shock occurred 9 hours after the first on the Sürgü fault and it had a M w 5.9 aftershock after 10 hours on the western end of the same fault system ( Figure 1B).
Within the entire aftershock sequence, the distribution of our event relocations indicates a spatially elongated set of events throughout the southwestern segment of the EAFZ. This includes epicenter of Nurdağı-Pazarcık earthquake and along the E-W oriented Sürgü fault following Ekinözü earthquake ( Figure 1C). Our relocations for the two mainshocks show 12.3 km of hypocenter depth falling within the upper crust for the M w 7.8 Nurdağı-Pazarcık earthquake whereas the Ekinözü earthquake is deeper at 15.2 km corresponding to the mid-crustal depth range. The depth distribution of the relocated aftershocks suggests the entire crust between 3 km to 25 km underwent deformation, mainly along major fault zones.

Kinematics of the M w 7.8 Nurdağı-Pazarcık earthquake
The event hypocenter is offset~15 km towards south from from the trace of the EAF. Additionally, there is a distinct cloud of aftershocks offset from the EAF. It is not feasible to associate it to the EAF given the good confidence in the hypocenter's location, and indeed, inversions that do assume this have very poor fits to the data. We infer thus that a secondary structure, the NPF, hosts the rupture initiation.
Consider Figure 3 where we show HR-GNSS site ONIY and strong motion site TK.2712 which are located 68 and 79 km away from the hypocenter (Figure 1). Here from these time-series, we find that there are clearly two stages of ground motion. This must be considered during the kinematic inversion. Thus, we tested two scenarios for how rupture transfers from the NPF to the EAF. First, we allowed rupture on the EAF that starts at a time equivalent to the moment when S-waves from the NPF reach it. In this case, the fits to the GNSS and strong motion were poor ( Figure S2), particularly regarding the early stages of the waveforms at near-field HR-GNSS sites ANTE and ONIY and strong motion sites 2712 and 2718 (locations in Figure S1) which were hard to model (e.g., Figures S2, 3, 4). In a second scenario, we delayed the onset of slip on the EAF until the time the rupture front originating at the NPF reaches the intersection of the two faults. Here we finally see the fits to the data improve significantly (Figure 4). Snapshots of rupture propagation ( Figure 5) and an animation (Supplementary S1, see Data and code availability) show that once the rupture reaches the EAF, at~10s after origin time, it spreads bilaterally across the fault. Slip rates reach as high as 1.5 m/s in the model, the total length of rupture on the EAF is~350 km and peak slip is 9 m -this yields a final moment of M 0 = 6.51x10 20 N-m (M w 7.8). The apparent complexity of the source time function is identified by many peaks reflecting the interaction of these two faults ( Figure 2B). Finally, we find that fits to the data are highest for v r max = 3.2 km/s which corresponds to about~90% of shear wave speed at the depths where most of the slip takes place.

Kinematics of the M w 7.6 Ekinözü earthquake
For the M w 7.6 rupture nucleates on the SF, spreads bilaterally ( Figure 5) and tapers at both ends of the fault (Figure 2A). The event has high peak slip,~7 m and the total length of rupture is~160 km. Fits to the waveforms are also good ( Figure 6) and have similar RMS ( Figure 2C) although there are later arrivals at strong motion sites TK.0205 and TK.4404 that cannot be modeled smoothly. These could reflect path or site-specific conditions that lead to amplifications that cannot be explained within our simple 1D approach. Nonetheless, the fits are good and the final model (shown in Figure 2) has a seismic moment of M 0 =3.64x10 20 N-m (M w 7.6). The most interesting aspect of this event is that the joint modeling of HR-GNSS and strong motion required two different v r max values in order to fit the data ( Figure 2C). We tried several values of a single v r max and quickly noticed that stations east or west of the rupture prefer different values ( Figure 2C). There is a broad plateau of low RMS between 4.6 and 5.4 km/s for sites west of the rupture. As an example, the ground motions recorded at stations TUF1 and FEEK (locations in Figure 1A) cannot be explained by sub-shear speeds. This preferred v r max is much larger than the~3.7 km/s shear wave speed at the depths where most of the slip takes place. This is compelling evidence that the event had super-shear rupture to the west. Interestingly the sites to the east display high misfits when v r max is high and prefer much lower values closer to~2.8 km/s. For our best-fitting model (Figure 2A,6) we imposed a combination of super-shear to the west and sub-shear to the east. This can be seen clearly in Animation S1 in the supplementary data as well.

Discussion and outstanding questions
The results shown here are a brief "first-look" analysis into two complex events and point to several important open questions, which will warrant further investigation. The M w 7.8 earthquake ruptured the southern three segments of the EAFZ which last broke in 1513, 1872, and 1893 (see Taymaz et al., 2021) and arrested at the source zone of the recent 2020 M w 6.7 Doğanyol-Sivrice earthquake (e.g., Melgar et al., 2020b). Meanwhile, the M w 7.6 likely broke the entire Sürgü fault which had not hosted a significant earthquake since 1544 (Taymaz et al., 2021). Understanding the timing, stress interactions between these events, and further implications for other neighboring structures will be important.
Regarding the ruptures, the strong evidence provided by the near-field HR-GNSS and strong motion data supports the conclusion that the second event involved a super-shear rupture, based on the relatively high estimate of v r max . The rapid finite-fault model published by the U.S. Geological Survey (2023b) similarly shows zones of super-shear rupture. Our preferred westward v r max of 4.8 km/s is very high but has been seen before in other super-shear strike slip events such as during the 1999 M w 7.4 Izmit, 2004 M w 7.8 Denali, 2013 M w 7.5 Craig, Alaska and other earthquakes (Bouchon et al., 2001;Frankel, 2004;Yue et al., 2013). Addition- ally, we notice again that this is the maximum allowed speed: slower speeds are possible with the multi-time window approach and indeed, in Figure 5, we observe that to the east the initial stage of rupture has very modest slip rates and the slower rupture speeds correspond to larger slip rates. The area where stations TUF1 and FEEK are located, towards the west, is where the slip pulse exhibits significant slip at v r max . Understanding the contributions of these kinematics to ground motion will be of great importance. Finally, a remaining open question is why there is no obvious super-shear rupture in the first event. Rosakis et al. (2023), from analysis of a strong motion record, suggest there must be super-shear rupture during the initial stage of the event on the Nurdağı-Pazarcık fault. Our data do not require this. However, we note that because our inversion relies on long period waveforms, it is possible for this early super-shear process to exist but to not be as obvious in the records. This is possible if the rupture transitions later to sub-shear velocities on the EAFZ. Indeed, on the EAFZ, rupture seems to prefer propagation right at Rayleigh wave speeds and, while increasing v r max to slightly above shear-wave speed still produces low RMS ( Figure 2C, S3), the result is not nearly as obvious or dramatic as for the M w 7.6 earthquake.
At a more granular level, a few structural questions remain as well. What is the exact nature of the NPF and how frequently does it participate in significant events? Additionally, the intersection of the SF and the EAF is structurally complex, mapped traces and the aftershocks hint at a secondary structure, sub-parallel to the EAF and immediately north of the SF. We find that rupture also branched out onto this structure. Here, remote sensing observations of crustal deformation will hold important clues. Using these data will not be without its challenges, as most observations, from InSAR for example, will have captured both events and many, if not all, the large aftershock. Separating the contributions to crustal deformation from individual events has been done for other similarly complex earthquake sequences (e.g., Taymaz et al., 2007;Fielding et al., 2013;Ganas et al., 2018Ganas et al., , 2021Goldberg et al., 2020;Taymaz et al., 2022) but it will require significant effort.

Conclusions
Here we have shown kinematic rupture models from joint inversion of HR-GNSS and strong motion data-sets and relocated aftershocks for the two events in the 2023 SE Türkiye earthquake doublet. We have used a com-plex multi-fault 3D geometry for inversion. We find that rupture speed is very close to the sub-to super-shear transition for the M w 7.8 event and that it is super-shear for the westward rupture of the M w 7.6 earthquake but sub-shear to the east. Peak slip exceeds 8 m for both events and slip rates as high as~1.5 m/s are pervasive throughout. Rupture lengths where~350 km for the M w 7.8 event and~160 km for the M w 7.6 earthquake. Sciences (TÜBA) in the framework for Young Scientist Award Program (TÜBA-GEBİP), and the Alexander von Humboldt Foundation Research Fellowship Award for financial support and for further providing computing facilities and other relevant computational resources through Humboldt-Stiftung Follow-Up Program. AG acknowledges funding from H2020 project European Plate Observing System Sustainability Phase. VT was funded by a Hellenic Foundation for Research and Innovation (ELIDEK) doctoral scholarship. SV received funding by the project "Risk and Resilience Assessment Center -Prefecture of East Macedonia and Thrace -Greece." (MIS 5047293) which is implemented under the Action "Reinforcement of the Research and Innovation Infrastructure", funded by the Operational Programme "Competitiveness, Entrepreneurship and Innovation" (NSRF 2014-2020) and co-financed by Greece and the European Union (European Regional Development Fund). We have further benefited from fruitful discussions with Tamer Y. Duman and Selim Özalp for interpretation of active neotectonics features observed in the region, with Ercan Yüksel and Oğuz C. Çelik for structural damages in the catastrophic area near the apocalypse, and with D. Goldberg for the kinematic modeling. TT thanks Dr. Beyza Taymaz for her phenomenal support during hectic days dealing with global media requests and organizing national and international scientific collaborations throughout four weeks of sleepless nights. We thank Kiran Kumar Thingbaijam, Ezgi Karasözen, an anonymous reviewer, and editor Théa Ragon for constructive comments.

Data and code availability
The hypoDD earthquake relocation code is freely available from the developers at https: //www.ldeo.columbia.edu/~felixw/hypoDD.html. The MudPy slip inversion code is available on GitHub at https://github.com/dmelgarm/MudPy and a snapshot of the code can be found on Zenodo at Melgar et al. (2021). Data and models produced in this paper can be found at https://zenodo.org/record/7699971#.ZATVUhPMLnwthis includes the catalog of relocated events, waveforms and station locations used for inversion, 3D geometry definitions of the fault, and the final best fitting rupture models. GNSS solutions were processed with JPL's GipsyX software package, licensed to BWC at University of Washington.