The Destructive Earthquake Doublet of 6 February 2023 in South‐Central Türkiye and Northwestern Syria: Initial Observations and Analyses

The two devastating earthquakes of 6 February 2023, and their associated aftershock sequences, in south‐central Türkiye and northwestern Syria are sobering reminders that earthquakes can neither be predicted nor prevented, but can only be prepared for. The earthquakes were the deadliest ones in Türkiye for centuries. The two strong earthquakes occurred in rapid succession but on different faults. The epicenter of the first shock of magnitude $Mw$ 7.8 was located ∼15 km east of the east Anatolian fault (EAF) at 37.288° N, 37.043° E, 8.6 km depth, (origin time 01:17 a.m. UTC). Only 09:07 hr later, the second event (magnitude $Mw$ 7.6) initiated 90 km north of the EAF on the east–west‐trending Sürgü fault (38.089° N, 37.239° E, 7.0 km depth, origin time 10:24 a.m. UTC).

The EAF defines the active plate boundary between the Arabian and the Anatolian plates (Fig. 1a). Over its ∼500 km length, the left‐lateral EAF has an estimated slip rate of ∼10 mm/yr (Aktug et al., 2016). Together with the right‐lateral north Anatolian fault (NAF), the EAF bounds the westward extrusion of the Anatolian plate from the Arabia–Eurasia collision zone (Pousse‐Beltran et al., 2020, and references therein). The section of the EAF that broke in the first mainshock extends into the Hatay triple junction between the EAF, the Cyprus arc, and the Dead Sea fault branching to the south (Fig. 1b).

During the last ∼100 yr, both the NAF and EAF varied in terms of releasing tectonic stress in large earthquakes. The NAF produced a sequence of large earthquakes in the twentieth century that initiated with the 1912 $Mw$ 7.2 Ganos earthquake at the western end of the Marmara Sea and then continued with the devastating $Mw$ 7.8 Erzincan earthquake in 1939 on the eastern NAF that killed over 30,000 people. This tragic quake was followed by 10 moderate‐to‐large events (1942–1967; Barka, 1996) and the 1999 $Mw$ 7.6 Izmit and $Mw$ 7.2 Düzce earthquakes east of the Marmara Sea, leaving major fault segments near Istanbul unbroken since 1766. In contrast, only three moderately sized earthquakes occurred on the EAF in the last ∼50 yr (1971 $Mw$ 6.7; 2020 $Mw$ 6.1; and 2020 $Mw$ 6.8), located on the northwestern section of the EAF. However, large historical earthquakes are documented along the southern EAF (Ambraseys, 2009; Meghraoui, 2015), such as in 1114 B.C.E. $Mw∼7.8$⁠, 1872 B.C.E. $Mw∼7.2$⁠, and 1822 B.C.E. $Mw∼7.5$⁠. The second rupture of the 2023 earthquake sequence occurred on the Sügür fault—a side branch of the EAF strand that is thought to have last ruptured in 1544 B.C.E. (Fig. 1a). The seismic activity of the NAF and EAF is reflected in corresponding seismic hazard maps (e.g., Akkar et al., 2018; Demircioğlu et al., 2018; Pagani et al., 2018; Şeşetyan et al., 2018; Fig. 1b).

In the past ∼100 yr, only a few continental strike‐slip earthquakes with magnitudes exceeding $Mw$ 7.5 occurred, limiting the available data and hence detailed studies of such large earthquakes. Examples are the 1990 $Mw$ 7.8 Luzon, Philippines, the 2001 $Mw$ 7.8 Kunlun, China, the 2002 $Mw$ 7.9 Denali, Alaska, the 2013 $Mw$ 7.7 Balouchistan, Pakistan, and the 2016 $Mw$ 7.8 Kaikoura, New Zealand, events. Each of these earthquakes revealed an intricate rupture process related to geometrical fault complexities (Klinger, 2022, and references therein): rupture lengths exceeding 100 km and reaching up to 300 km (2001 Kunlun earthquake); average horizontal surface slip of 3–4 m, reaching locally up to 7–9 m, and strongly varying along‐strike in relation to the fault‐trace geometry.

The unique character of the 2023 sequence is that two large‐magnitude earthquakes occurred only 9 hr apart on nearby faults. Pairing of large continental earthquakes over such a short time had not been observed before. Previous pairing always involved longer separation in time, like 14 days between the $Mw$ ∼8 Tsetserleg and $Mw$ ∼8 Bulnay events (Choi et al., 2018) and 4 months for the $Mw$ 7.4 Izmit–$Mw$ 7.2 Düzce sequence (Konca et al., 2010). In addition, the spatial dimensions of the two main 2023 quakes estimated from real‐time aftershock locations (Fig. 1b) reach lengths of ∼350 km on the EAF for the initial $Mw$ 7.8 earthquake and ∼170 km for the $Mw$ 7.6 earthquake on the Sürgü fault. For the regional seismogenic width of 20 km (Ozer et al., 2019), these source dimensions are consistent with the events’ magnitudes based on source‐scaling relations (Thingbaijam et al., 2017).

We used pixel‐offset tracking of Sentinel‐1 radar images to map coseismic surface displacements around the two faults and the extent of surface fault rupturing (e.g., Fialko et al., 2001; Wang and Jónsson, 2015). Based on ascending and descending orbit images, as well as along‐track (azimuth) and across‐track (range), we derived pixel offsets (Fig. S1, available in the supplemental material to this article), yielding four different offset images from which we inverted for 3D surface displacements (Liu et al., 2022; Fig. S2). The resulting horizontal surface displacements and their spatial pattern exhibit left‐lateral motion across the two main faults (Fig. 2a). Vertical displacements are small in comparison (Fig. S2c), confirming the almost pure strike‐slip mechanism of both events. The length of the main surface rupture along the EAF in the first earthquake is ∼320 km, whereas the surface rupture of the second mainshock is markedly shorter (∼150 km). Hence, for both cases the mapped surface rupture is 20–30 km shorter than indicated by aftershock locations (Fig. 1d).

Both earthquakes produced large surface‐fault offsets (exceeding 4 m) over extended sections along the faults (Fig. 2). From their horizontal surface displacement fields, we measured fault offsets at 5 km intervals along the two main faults. The results show that surface fault slip along the EAF has 2–3 slip maxima, with the largest slip found northeast of the epicenter (6–7 m), ∼30 km east of the city of Kahramanmaras. Another slip maximum (∼4 m) is found farther southwest, near Islahiye, with fault slip abruptly decreasing near Antakia at the southwestern end of the rupture. The maximum surface offset of the second fault is even larger than for the first rupture, exceeding 7 m near the epicenter and 6 m over a fault length of ∼60 km (Fig. 2b). This large surface offset may compensate the relatively short rupture length for the measured magnitude of $Mw$ 7.6.

Using far‐field teleseismic data, we estimated the space–time rupture evolution via backprojection (e.g., Ishii et al., 2005; Koper et al., 2011; Li and Ghosh, 2017) and finite‐fault modeling (e.g., Mendoza and Hartzell, 2013, and references therein). For this purpose, we compiled two datasets. For the backprojection, we used seismic stations in Alaska and selected only stations with average cross correlation (CC) above 0.6 for the first 25 s around the P‐phase arrival filtered in the range 0.1–2 Hz. This results in 205 and 201 stations for the $Mw$ 7.8 and 7.6 events, respectively. The targeted backprojection region extends from 35° to 39.5°, both in latitude and longitude, with 0.01° grid spacing in both the directions. We calculated theoretical travel times based on the Preliminary Reference Earth model from the source grid to each seismic station, with source depths fixed at catalog depths. In addition, we applied time shifts and relative polarity estimation from the peak cross‐correlation (CC) coefficients of the first arrival P phases for a shorter 5 s time window, relative to a reference station with maximum average CC coefficients, as empirical time and polarity correction. To then image the rupture evolution, we deployed a 6 s sliding time window and 0.1 s time step to the continuous waveform data.

Backprojection results show bilateral rupture of the $Mw$ 7.8 earthquake, with an average rupture speed ∼2.5 km/s to the east and ∼2 km/s to the west, estimated using epicentral distance (which underestimates average rupture velocity along the fault itself if fault geometry changes; Fig. S3). However, with a priori knowledge of the fault traces and assuming nearly vertical dip angles, we were able to backproject the radiated seismic energy directly to the fault (Fig. 3). This better illustrates the complex rupture process of the $Mw$ 7.8 event (inset in Fig. 3f). The backprojection suggests bilateral rupture on a small fault east of the EAF where the hypocenter is located (Disaster and Emergency Management Authority [AFAD] catalog). Rupture to the southwest stopped after a few seconds but continued to the northeast until it reached the intersection with the EAF at ∼10 s. The backprojection then locates strong radiation to the east, but near the intersection of the nucleation branch and the EAF. This correlates with the strongest recorded ground shaking at the station TK4614. The rupture continued to the northeast with an average rupture speed of ∼3.1 km/s along the EAF until ∼55 s rupture duration, with strong radiation from a short east–west branch of the fault, the junction between the EAF and the eastern extension of the Sürgü fault before it stopped ∼25 km farther east of the junction (inset in Fig. 3f). The rupture to the southwest along the main EAF appears delayed with limited seismic radiation west of the epicentral region. However, continued rupture to the east may have altered the stress state on this segment, thereby promoting further rupture to the southwest until ∼80 s, when rupture suddenly terminated near Antakia with strong observed seismic radiation.

Backprojection results of the $Mw$ 7.6 event reveal a frequency‐dependent rupture process due to rupture directivity (Fig. 3e–h and Fig. S3; Li et al., 2022). The 0.1–0.5 Hz results capture rupture to the east of the Sürgü fault, which then changed its direction toward the northeast. The 0.5–1 Hz results, on the other hand, mainly track rupture to the west on the Sürgü fault. The two strong radiation sources are located where the fault geometry changes.

In addition to backprojection imaging, we estimated two sets of finite‐fault source models for the two earthquakes, one from satellite radar data‐derived coseismic horizontal surface displacements (Fig. 2a) and the other from teleseismic observations. For the geodetic source model, we used the fault traces mapped from satellite radar offsets, extended the fault lengths a few kilometers beyond the mapped surface ruptures, extended the fault widths to 25 km, and discretized them into 5 km × 5 km fault patches. The first fault is vertical, whereas the second mainshock fault dips 78° to the north. We then estimated spatially variable slip on the faults (e.g., Jónsson et al., 2002) from the coseismic horizontal surface displacements using an appropriate degree of spatial fault‐slip smoothing.

For inverting the teleseismic data to derive kinematic finite‐fault source models, we downloaded seismic waveforms for stations situated at teleseismic distances of 30°–90°, ensuring good azimuthal coverage. By visual inspection, we selected 18 stations for the first mainshock (Fig. S4a) and 17 stations for the second event (Fig. S5a) with high signal quality, using initially only the P wavetrain. Waveforms were then band‐pass filtered (5–20 s, Butterworth filter) to remove high‐frequency noise. To infer kinematic finite‐fault parameters, we only used vertical components and applied covariance matrices to account for errors in both measurements and theory (e.g., Vasyura‐Batke et al., 2020, and references therein). Because the kinematic finite‐fault source‐parameter estimation suffers from nonunique solutions, we explored the model space using Bayesian inference with sequential Monte Carlo sampling implemented in the Bayesian earthquake analysis tool (BEAT) code (Vasyura‐Batke et al., 2020; Figs. S4b and S5b).

The kinematic rupture model for the first main event comprises four major segments with uniform dip angle of 89°. Each segment is subdivided into 5 km × 5 km subfaults. Segment 1‐a, on which the hypocenter is located, has dimensions of 55 km × 25 km, segment 1‐b and 1‐c along the EAF expand over 90 km × 25 km, and 80 km × 25 km, respectively, whereas the southernmost segment 1‐d covers an area of 140 km × 25 km (Fig. 4). In total, the four fault segments form a 365 km long rupture plane that extends to 25 km depth, parameterized by 365 subfaults.

For the 2nd mainshock, we discretized the fault rupture model into three segments, each subdivided into 5 km × 5 km subfaults. Segment 2‐a, on which rupture nucleated, spans 70 km length; segment 2‐b (to the east) and segment 2‐c (to the west) are 70 km and 50 km long, respectively (Fig. 4a). Each segment extends for 25 km along the fault‐dip direction (78° dip angle). The assumed rupture plane is thus 190 km long.

On each subfault, we solved for slip, rupture onset time, and rise time. Furthermore, we allowed for variable rupture velocity, searching in the range of 2.5–4.0 km/s. The local source time function was set to a half cosine; the slip direction (rake angle) was allowed to vary in the range −90° to 90° with respect to a reference rake angle of 0° (parallel to the strike direction).

Figure 4 summarizes the finite‐fault rupture models inverted from coseismic surface displacements and teleseismic P‐wave data. We did not perform any joint inversion, so each model is constrained by a single dataset. However, both the models reveal consistent slip distributions. The first $Mw$ 7.8 earthquake is characterized by three main areas of fault slip, showing up to 7 m of slip near the surface on segment 1‐b, up to 6 m on segment 1‐c, and 4–7 m on segment 1‐d (in which we find the largest difference between the seismic and the geodetic models). Segment 1‐a, on which the rupture started, had less fault slip but still up to ∼3 m (Fig. S4c). The southernmost segment 1‐d shows an area of high slip before rupture abruptly stops, in agreement with the backprojection results, creating a strong stopping phase. Together with rupture directivity along this 140 km fault segment, strong seismic radiation was generated toward the south into the Hatay‐Antakia region that combined with local site effects, created severe local shaking and extensive damage. The total finite‐fault seismic moment for this rupture models is $1.03×1021N·m$ (⁠$Mw$ 7.97) and $6.13×1020N·m$ (⁠$Mw$ 7.84) from the teleseismic and geodetic data, respectively consistent with other estimates (i.e., Jiang et al., 2023).

The backprojection imaging and teleseismic source inversion thus reveal a “T‐Bone” geometry with rupture propagating backward relative to the initial direction, seen only in few previous cases, that is, Kaikoura (Klinger et al., 2018), Romanche (Hicks et al., 2020), and to a lesser extent Landers event (e.g., Fliss et al., 2005; Wollherr et al., 2019).

Finite rupture models of the second mainshock show very large near‐surface fault slip with the maximum slip exceeding 8 m on segment 2‐a (on which rupture nucleated, see Fig. S5c) along the Sürgü fault, with slip values reaching 6 m over an extended stretch along that segment. Slip values then decreased toward northeast and southwest along segments 2‐b and 2‐c, respectively. These inferred slip values are in good agreement with surface displacement derived from geodetic data (Fig. 2). The inferred seismic moment of this rupture model is $5.03×1020N·m$ (⁠$Mw$ 7.77) and $3.32×1020N·m$ (⁠$Mw$ 7.65) for the seismic and geodetic model, respectively.

We collected strong‐motion recordings from the AFAD (see Data and Resources) for first $Mw$ 7.8 earthquake at 254 stations based on the following selection criteria: (1) instrument response removed, band‐pass filter applied (low cutoff frequency: 0.025–0.1 Hz and high cutoff frequency: 25–40 Hz), and baseline corrected; (2) no abnormal recordings (e.g., no pre‐event signals and no obvious peaks); (3) three‐component recordings available; and (4) $VS30$ values available. Moreover, we obtained regional‐distance seismic waveforms from several sites located along the southward extension of the EAF, the Dead Sea fault, and the Gulf of Aqaba, including stations in Saudi Arabia. We removed the instrument response from waveforms, and then filtered them between 0.01 and 50 Hz.

Figure 5 presents an overview of locations of strong‐motion sites (Fig. 5a), selected examples of near‐source recordings that illustrate pulse‐like motions due to directivity effects and long‐duration shaking (Fig. 5b), as well as regional‐scale waveforms with well‐developed surface waves (Fig. 5c). The peak ground acceleration (PGA) ShakeMap in Figure 5a documents PGAs exceeding 0.5g in mainly three areas: near Adiyaman, around the wider epicentral region, and over a large area in the Hatay‐Antakia region. Locally, PGA values reached 1g, with one site (TK4614) even showing ∼2g horizontal ground acceleration. In addition, we collected strong‐motion recordings from 150 stations for the second mainshock, applying the same criteria as for the first event. First‐order analysis of PGA values of the second event reveals overall lower ground motions (the maximum recorded PGA 0.56g at site TK4612, the closest to the epicenter at ∼67 km distance; Fig. S12a). However, due to the lack of stations closer to the fault, even higher shaking levels may have occurred but were not recorded.

We further compared observations with two empirical ground‐motion models (GMM) used in the 2018 Turkish probabilistic seismic hazard assessment (PSHA; Akkar et al., 2018) (Figs. S6–S11). Whereas the first GMM is specific for Türkiye (Akkar and Çağnan, 2010), the second one applies to Europe and the Middle East (Akkaret al., 2014). Our preliminary analyses suggested that observed ground motions exceed median GMM predictions not only for these two GMMs, both for “raw” observations, and if site‐specific corrections for $VS30$‐based site amplification are applied (Figs. S7, S9, and S11). These observations are consistent with Gülerce et al. (2016) who modified the Next Generation Attenuation‐West1 (NGA‐W1) GMMs using the Turkish strong‐motion database (so‐called “TR‐adjusted models”). These TR‐adjusted GMMs better replicate recorded ground motions, because they adopted the well‐constrained large‐magnitude scaling of the global dataset in the NGA‐W1 models. ShakeMaps for two spectral periods (T = 0.2 s and T = 1.0 s) of spectral acceleration SA(T) reveal the concentrated strong shaking in several areas (Fig. S13a–d). The regions of particularly high shaking levels correspond approximately to fault areas with high slip, whereby extended strong shaking in the Hatay‐Antakia region can be explained by a combination of strong seismic radiation and local site effects. At several sites (e.g., Antakia, Iskenderun, and Arsuz), spectral accelerations exceeded the current building code of Türkiye (TBEC‐18) at periods T >1 s relevant for tall structures (Figs. S13e,f).

We conducted a first‐order analysis of the rupture process of the magnitude $Mw$ 7.8 and 7.6 earthquake doublet of 6 February 2023 in south‐central Türkiye using both satellite and seismic data. Both earthquakes are large, predominantly bilateral strike‐slip ruptures. The $Mw$ 7.8 earthquake initiated on a side branch to the EAF and transitioned onto the main EAF with bilateral rupture into the northeast and southwest directions. Although the event stopped abruptly in the northeast (after ∼55 s), rupture continued to the south where it then terminated after ∼80 s. Directivity effects due to rupture propagation along extended straight fault segments as well as stopping phases due to sudden rupture cessation at fault extremities led to locally strong seismic radiation for the $Mw$ 7.8 earthquake. The $Mw$ 7.6 earthquake initiated on the Sürgü fault, which is 90 km north of the $Mw$ 7.8 epicenter and ruptured bilaterally for about 150 km. Given its magnitude, the $Mw$ 7.6 rupture is shorter and more compact. Considering the length of both ruptures and their strike‐slip mechanisms, supershear rupture propagation may locally be expected. We found evidence for such behavior in the backprojection imaging for $Mw$ 7.6 event, but refined analyses based on strong‐motion records is needed to confirm this initial observation.

PGAs during the $Mw$ 7.8 earthquake locally reached 2g and exceeded 0.5g over a wide area in the Hatay‐Antakia region. Directivity effects and strong stopping phases are partially responsible for the observed strong‐motion characteristics. Site effects further amplified ground motions locally. An initial analysis reveals that shaking levels exceeded median predictions from GMMs used in the most recent regional PSHA (Figs. S6–S11). Locally, observed spectral accelerations exceeded the design spectra of the current building code (Figs. S12a and S13). Ground motions of the second mainshock then hit already weakened or partially collapsed buildings and infrastructure, further increasing damage and destruction. In combination, these effects may provide partial explanations for the widespread damage and large destruction of these two earthquakes.

Although the occurrence of two such large earthquakes as a “doublet” is uncommon, the second event can be physically explained by stress changes in its epicentral area imposed by the first mainshock that brought the fault closer to failure (Stein et al., 2023). Given size and location, we consider the $Mw$ 7.6 earthquake therefore a second mainshock and not an unusually large aftershock. Large strike‐slip earthquakes like the $Mw$ 7.8 and 7.6 ruptures of 6 February 2023 are rare but not uncommon, because they have been observed in the past. Such multisegment ruptures forming “compounded” events on geometrically complicated fault structures are a challenge in standard PSHA.

Real‐time aftershock locations are provided by the Disaster and Emergency Management Authority (AFAD) available at https://deprem.afad.gov.tr/event-catalog. Focal mechanism solutions of the two mainshocks are available at https://geofon.gfz-potsdam.de/old/eqinfo/list.php?mode=mt. Focal mechanisms of significant seismicity during the period 1 January 2018–6 February 2023 can be downloaded at https://www.globalcmt.org/. Satellite data are available via the Sentinel data hub at https://scihub.copernicus.eu. Teleseimic waveforms for backprojection were downloaded from the Incorporated Research Institutions for Seismology (IRIS) at https://ds.iris.edu/wilber3/find_stations/11654089 for the $Mw$ 7.8 event and https://ds.iris.edu/wilber3/find_stations/11654205 for the $Mw$ 7.6 event. Teleseismic waveforms for finite‐fault inversion were obtained from the IRIS at https://www.iris.edu, Geoforschungsnetz (GEOFON) at https://geofon.gfz-potsdam.de, and Observatories and Research Facilities for European Seismology (ORFEUS) at https://www.orfeus-eu.org data centers, respectively. Broadband waveforms of Figure 5 are from King Abdullah University of Science and Technology (KAUST) seismic network (COLA, ASCO, available upon request to the authors), GEOFON (EIL, MSBI, and UJAP) available at https://geofon.gfz-potsdam.de, and Lebanese Centre National de la Recherche Scientifique‐École Normale Supérieure (CNRS; BHL, available upon request to the authors). Strong‐motion data are available at https://tadas.afad.gov.tr/event-detail/15499. The 2018 seismic hazard map of Türkiye is available at https://www.resmigazete.gov.tr/eskiler/2018/03/20180318M1-2-1.pdf. High‐resolution aftershock locations are available at “A. Lomax (2023). Precise, NLL‐SSST‐coherence hypocenter catalog for the 2023 Mw 7.8 and Mw 7.6 SE Turkey earthquake sequence. (v2.0) [Data set]. Zenodo, doi: 10.5281/zenodo.7727678”. The supplemental material provides additional material to explain results and support discussions in the main paper. Specifically, the supplemental material includes figures on the processing of Sentinel‐1 radar images, further information on the teleseismic back‐projection and finite‐fault inversions, and detailed comparisons of the observed ground motions with empirical ground‐motion models and the design spectra of the Turkish building code. All websites were last accessed in February 2023.

The authors acknowledge that there are no conflicts of interest recorded.

Geçtiğimiz bu zor günlerde tüm kalbimizle Türkiye ve Suriye vatandaşlarının acılarını paylaşıyoruz. Bu yıkıcı depremlerden etkilenen, ailelerini, arkadaşlarını, komşularını, evlerini ve eşyalarını kaybeden herkese başsağlığı diliyoruz.

قلوبنا متضامنة مع الشعب التركي والسوري في هذه الأوقات العصيبة.

نقدم تعازينا الخالصة لكل من فقد صديق أو قريب أو جار أو منزل أو ممتلكات بسبب الزلازل.

(We share with all our hearts the suffering of the citizens of Türkiye and Syria during these difficult days. Our condolences to all those affected by these devastating earthquakes, who lost their families, friends, neighbors, homes, and belongings.)

We are indebted to the Saudi Geological Survey (SGS) and the National Center for Geophysics from Centre National de la Recherche Scientifique‐École Normale Supérieure (CNRS) of Lebanon for sharing regional distance seismic waveforms. Comments and constructive criticism by two anonymous reviewers and the Associate Editor helped to improve our work. The research presented in this article is supported by King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia, by Research Grant Numbers BAS/1/1339‐01‐01 and BAS/1/1353‐01‐01.