# The overall-subshear and multi-segment rupture of the 2023 Mw7.8 Kahramanmaraş, Turkey earthquake in millennia supercycle

Liuwei Xu, Saeed Mohanna, Lingsen Meng, Chen Ji, Jean-Paul Ampuero, Zhang Yunjun, Masooma Hasnain, Risheng Chu, Cunren Liang

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

## To cite this version:

Liuwei Xu, Saeed Mohanna, Lingsen Meng, Chen Ji, Jean-Paul Ampuero, et al.. The overall-subshear and multi-segment rupture of the 2023 Mw7.8 Kahramanmaraş, Turkey earthquake in millennia supercycle. Communications Earth & Environment, 2023, 4 (1), pp.379. 10.1038/s43247-023-01030-x. hal-04885069

HAL Id: hal-04885069

https://hal.science/hal-04885069v1

Submitted on 15 Jan 2025

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

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

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

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

# The overall-subshear and multi-segment rupture of the 2023 Mw7.8 Kahramanmaraş, Turkey earthquake in millennia supercycle

Liuwei Xu ID <fnref n="1" />, Saeed Mohanna ID <fnref n="1" />, Lingsen Meng ID <fnref n="1" />✉, Chen Ji ID ²✉, Jean-Paul Ampuero ID ³, Zhang Yunjun ID ⁴, Masooma Hasnain ID <fnref n="1" />, Risheng Chu ID ⁵ & Cunren Liang⁶

On February 6, 2023, an Mw7.8 earthquake hit the East Anatolian Fault (EAF) and Narlı Fault (NF), followed by an Mw7.5 event on the Sürgü Fault. We combine multiple seismic datasets, global navigation satellite system recordings, and radar satellite images with finite fault inversion and slowness enhanced back-projection to study the rupture kinematics. Our analysis reveals that the rupture originated on the NF, propagating 120 km northeast at 3.05 km/s and 200 km southwest at 3.11 km/s after reaching the EAF junction, exhibiting overall subshear speeds. Further Mach wave analysis confirms the subshear rupture, matching the prediction using close-Rayleigh speeds. The unexpectedly-large slip on some EAF segments suggests a supercycle lasting ≥900 years. The EAF geometry is similar to the San Andreas-San Jacinto Fault system, while the latter has higher slip rates but without large earthquakes on its southern segments since 1857, carrying the potential of an M8 earthquake.

<footnote n="1">¹Department of Earth, Planetary and Space Sciences, University of California Los Angeles, Los Angeles, CA, USA. ²Department of Earth Science, University of California Santa Barbara, Santa Barbara, CA, USA. ³Observatoire de la Côte d'Azur, Université Côte d'Azur, IRD, CNRS, Valbonne, France. ⁴National Key Laboratory of Microwave Imaging Technology, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing, China. ⁵State Key Laboratory of Geodesy and Earth's Dynamics, Innovation Academy for Precision Measurement Science and Technology, Chinese Academy of Sciences, Wuhan, China. ⁶Institute of Remote Sensing and Geographical Information System, School of Earth and Space Sciences, Peking University, Beijing, China.
✉email: lsmeng@g.ucla.edu; ji@geol.ucsb.edu</footnote>
On February 6, 2023, two powerful earthquakes of moment magnitude (Mw) 7.8 and 7.5 hit south-central Turkey and northern Syria (Fig. 1). The earthquakes, which occurred only 9 h apart, shook the regions around Gaziantep, Kahramanmaraş, Malatya, and Hatay, severely affecting the local population. These earthquakes produced ruptures that are among the longest ever recorded in continental strike-slip earthquakes, with fault lengths of ~368 km and ~133 km for the Mw 7.8 and Mw 7.5 ruptures, respectively¹. The earthquakes triggered a basin-wide tsunami alert, and a small tsunami was generated in the Eastern Mediterranean Sea. The two earthquakes resulted in severe damage to buildings and infrastructure in Turkey and Syria, with a reported loss of at least 50,000 lives and many more injured².

The earthquake doublet occurred in the East Anatolian Fault Zone (EAFZ), and the rupture reached the northern end of the Dead Sea Fault Zone (DSFZ). The DSFZ is a major tectonic boundary between the African and Arabian plates in the Middle East. It has a complex geometry and has caused devastating earthquakes in the past, including the 1138 M7.1 Aleppo, 1837 M7.1 Galilee, and 1927 M6.3 Jericho earthquakes. The EAFZ is a major left-lateral strike-slip fault that runs for approximately 700 km (ref. ³) along the eastern edge of the Anatolian Plateau. It has been responsible for numerous large earthquakes in the past, including the 995 (M7.4), 1114 (M ≥ 7.8), 1789 (M7.4), 1893 (M7.1), 1905 (M6.8) events⁴⁻⁶. More recently, there has been a series of damaging earthquakes across southeastern Turkey including the 1998 M6.3 Adana-Ceyhan, 2003 M6.4 Bingol, 2011 M7.1 Van, 2010 M6.1, and 2020 M6.7 Elazig earthquakes. Fault mapping of the EAFZ shows several fault segments and splay faults that have accommodated varying degrees of seismicity throughout their histories³,⁷. Prior to the 2023 earthquakes, the 2020 Elazig earthquake was the most recent M > 6.5 event to

rupture the EAF. It occurred further east, on the Pütürge segment of the EAF, and it was not estimated to increase the Coulomb stress on the Pazarçık segment⁸,⁹. The 2023 Mw 7.8 earthquake started close to the junction of the Pazarçık and Amonos segments of the EAF. The Disaster and Emergency Management Authority of Turkey (AFAD) placed its epicenter close to a splay fault referred to as the Narlı Fault (NF). The Mw 7.5 aftershock occurred about 90 km north, on the Cardak fault segment of the Sürgü Fault.

To analyze the source processes associated with the Kahramanmaraş earthquakes in Turkey, we utilize two methodologies: seismic array Slowness Enhanced Back-Projection (SEBP)¹⁰ and joint Finite Fault Inversion (FFI)¹¹,¹². SEBP is applied to both the Mw 7.8 mainshock and the Mw 7.5 aftershock, while FFI is specifically applied to the Mw 7.8 mainshock and resolves the slip distribution on the NF and EAF. SEBP and joint FFI indicate that during the Mw 7.8 mainshock, the rupture first initiates and propagates on the NF. After reaching the junction with the EAF, it propagates bilaterally at overall subshear speeds, extending 120 km to the northeast at 3.05 km/s and 200 km to the southwest at 3.11 km/s. We further validate subshear speeds by local seismic recordings and far-field Rayleigh waves analysis. Our investigations reveal the presence of an earthquake supercycle on the EAF whose complete period is ≥900 years. We discover geometric similarities between the EAF system and the San Andreas Fault system. These findings hold crucial implications for seismic hazard evaluations in California and other global strike-slip fault systems.

## Results and discussion

**Rupture kinematics of the Mw 7.8 and 7.5 events.** We apply the SEBP and joint FFI to analyze the source processes associated with the 2023 Mw 7.8 Kahramanmaraş, Turkey earthquake and

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

**Fig. 1** Spatiotemporal distribution of high-frequency (High-Freq) radiators and multiple datasets used in joint FFI. Diamonds denote the High-Freq radiators for the Mw 7.8 earthquake, color-coded by rupture time relative to the origin time of the event and with size proportional to the normalized BP power. Radiators imaged by the China array are shown in **a** and by the Alaska array are shown in **b**. All China and Alaska stations are shown in the lower inset of **a**. Color circles denote the same as diamonds but for the Mw 7.5 earthquake. Blue to red background shows the ground displacement in east-west direction (**a**) and in north-south direction (**b**) from radar satellites (Methods). **a** The gray dots denote the seismicities occurring from Jan 1st, 2023 to Mar 14th, 2023, from the AFAD catalog. The white triangles denote the strong motion stations used in SEBP validation (Methods). The diamonds with magenta edges correspond to the 5 s radiator (see Results section). **b** S1-S6 denote the vertical fault planes adopted for FFI. The red arrows indicate the coseismic deformation measured by GNSS stations. The white and magenta triangles denote the strong motion stations and teleseismic stations used by joint FFI, respectively. The segments of the EAF by Duman and Emre³² are denoted by the brown dash lines and labeled names. **c** The tectonic setting of the 2023 Mw 7.8 Kahramanmaraş earthquake. The red rectangle indicates the study region in **a** and **b**. The fault trace data are from Emre et al.¹⁴, Mahmoud et al.⁹⁵, and Styron et al.⁹⁶.
apply SEBP to the Mw 7.5 event. The slip planes of the Mw 7.8 Kahramanmaraş earthquake are constructed on six vertical fault segments (segments S1–S6, Fig. 1b) extracted from the SAR image-based surface traces (Fig. 1a, b) and aftershock distributions. During the joint FFI study, we divide these fault segments into subfaults and invert for the rupture histories of individual subfaults simultaneously¹¹,¹² using seismic waveforms observed in both close-fault (<50 km) strong motion and teleseismic (3000–10,000 km) stations as well as the ground deformation captured by high-rate Global Navigation Satellite System (GNSS), and Synthetic Aperture Radar (SAR) speckle tracking from Sentinel-1 satellite, InSAR and multiple aperture interferometry (MAI) from ALOS-2 and LuTan-1 satellites (Supplementary Figs. 1–8). The SEBP analysis is based on the technique developed by Meng et al.¹⁰ (see Methods for more details). We collect the broadband seismograms from two large-aperture arrays at teleseismic distance, located in China (CH) and Alaska (AK; the lower inset in Fig. 1a). Slowness correction terms for the mainshock are derived using nearby aftershocks (Supplementary Table 1 and Supplementary Figs. 9 and 10) to further mitigate the spatial bias caused by path effects (Methods).

To account for the influence of random seed selection in simulated annealing inversions, we conduct 10 FFIs for the Mw 7.8 event using different random seeds. This approach helps us explore the uncertainty arising from multiple optimal solutions within the model space. All 10 inversions result in similar solutions as expected, with negligible standard deviation on the objective function values (~0.6% of the average), coseismic slip distributions and moment rate functions (Supplementary Fig. 11). Considering that all 10 models are plausible solutions, we present the model with the smallest objective function value. Its coseismic slip distribution is shown in Fig. 2 and Supplementary Figs. 12

and 13, the rupture evolution on fault planes is presented in Supplementary Movie 1, and the coseismic slip rate is shown in Supplementary Fig. 13b. The final coseismic moment is $M_0 = 7.67 \times 10^{20}$ Nm (Mw 7.85), with an average slip of 3.1 m. The preferred slip model features three large slip asperities north and south of the hypocenter, one located on the EAF near the junction with the splay fault, another at 70 km northeast of the junction, and one at the south Pazarcik and north Amanos segments. The peak slip is 9.2 m and occurs on the Pazarcik segment. The south rupture on the Amanos segment of the EAF has relatively uniform and shallow slip, with an average value of ~3.1 m. The along-fault averaged slip as a function of depth is shown in Supplementary Fig. 14. The largest average slip occurs at ~4 km depth, with a value of 4.4 m. According to the FFI, the rupture duration is ~100 s, with the southwest rupture tip reaching the southern Amanos fault, close to the Karasu Fault (Fig. 1). The total rupture length is estimated to be ~370 km, based on our SEBP results. This value is very close to that estimated by field and satellite data (~368 km, ref. ¹).

Our FFI analysis reveals that the 2023 Mw 7.8 Kahramanmaras mainshock has a heterogeneous rupture process. As shown in Supplementary Fig. 13a and Movie 1, the rupture starts on the splay NF, but 92% of the seismic moment is released on the Erkenek, Pazarcik, and Amanos segments of the EAF as three high slip patches (Fig. 2 and Supplementary Fig. 15). The rupture first propagates on the ~N33.5°E oriented NF for ~50 km and reaches the NF’s northern end at ~20 s, consistent with the SEBP results of the AK and CH arrays (Fig. 3b). At the junction of the NF and the EAF, our slip model shows spatially and temporally complex rupture fronts. The earliest slip on the EAF occurs at ~15 s, then propagates bilaterally towards the northeast and southwest. The northeast rupture lasts ~40 s and extends to the

<figure><img src="image_5.png" /><figcaption>**Fig. 2 3D view of the preferred slip model.** Segment indices and lengths are labeled at bottom. Fault names are annotated on the top. The yellow star denotes the hypocenter.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Spatial and temporal evolution of High-Freq radiators and slip rates. **a** The moment rate function (MR, black line) and the BP power for SEBP results of CH array (purple line) and AK array (blue line). **b** Along-fault distance and time of BP radiators and slip rates. The junction of the splay fault and the main fault is set as the origin of the distance axis, and the northeast is the positive direction. The origin time of the Mw 7.8 event is set as the origin of the time axis. Unless otherwise stated, subsequent references to distance and time use this definition. The diamonds and circles denote the High-Freq radiators imaged by the CH and AK arrays, respectively. The colormap denotes the slip rate averaged across depth. The red and magenta slants show the fitted rupture speeds for the SW and NE fronts delineated by High-Freq radiators, respectively. **c** Vertically exaggerated envelopes of detrended, demeaned, bandpass filtered (1-4 Hz) fault-parallel components of strong motion stations shown in Fig. 1a. The red line indicates the SW speed fit using the timing and distance of second envelope peaks at southern stations. The red dots indicate the time at which these peaks occur for a particular station's distance along the fault.</figcaption></figure>

Pütürçe segment. The southwest rupture continues on the Pazarcık and Amanos segments for ~80 s after the junction. To better compare the rupture evolution imaged by the SEBP and FFI, we project the high-frequency (High-Freq) radiators imaged by SEBP onto the FFI fault planes (Fig. 3b). The subfaults' slip rates, defined as the slip amplitude divided by the rise time, are averaged across the depth and plotted as a function of their curvilinear distance along the FFI faults. The High-Freq radiators imaged by SEBP show an overall subshear northeast rupture process consistent with the FFI slip model, with an average propagation speed of ~3.05 km/s. The first 80 km of the southwest rupture is only partially imaged by SEBP and is obscured by the strong northeast front; this interference is expected when the separation distance between radiators of a bilateral rupture is shorter than the BP resolution¹³. The interference of the bilateral fronts also slightly affects BP imaging between −80 and −150 km, during which the slip in FFI appears 1–5 s earlier than High-Freq radiators. The overall south rupture speed resolved by SEBP is ~3.11 km/s, which also indicates a subshear rupture. A secondary BP energy peak appears in the power curves between 120 and 140 s, after the trough at 100 s (Fig. 3a, Supplementary Movies 2–5). This peak corresponds to the radiators in the Eastern Mediterranean Sea (Fig. 1a, b) and could be generated by an immediate aftershock or a subevent on the offshore part of the Karasu fault, close to the Cyprus arc (see Rupture Speed and Extent section for more discussion).

The SEBP results are further validated by local strong motion recordings (Fig. 3c). The recordings show two pulses of energy that form easily distinguishable peaks at most displayed stations. The first pulse attenuates with distance from the epicenter as expected from the geometrical spreading of seismic waves. We trace S-phases radiated from selected High-Freq radiators obtained from the BP results of the CH array to test whether the arrival times of these S-phases at station locations coincide with time frames of high energy release. The magenta vertical bars indicating the theoretical S-wave arrival from the 5 s High-Freq radiator (the diamond with a magenta edge in Fig. 1a) coincide well with this first energy pulse. However, the second

pulse of energy does not attenuate with distance. In Fig. 3c, we notice an abundance of High-Freq radiators, indicated by the cyan diamonds, at southern stations at the time of the second energy pulse. They correspond to the SW-ward rupture front passing through these stations along the Amonos segment. The amplitude of the second pulse of energy is thus maintained by the contribution of source areas closest to each station. We do not see two envelopes of energy at stations in the north, which we attribute to the rupture pattern observed from SEBP: the north rupture travels at speed comparable to the local shear wave speed (Vs = 3.39 km/s, ref. ³), so that S waves from the hypocenter and the rupture front arrive near northern stations almost simultaneously. In contrast, the south stations see two envelope peaks separately because the SW rupture starts at least 15 s after the earthquake's initiation and north from the epicenter.

For the SEBP of the Mw 7.5 event (Fig. 1 and Supplementary Fig. 10, Supplementary Movies 6–9), the CH array result shows that the rupture propagates bilaterally along the East-West trending Sürügü Fault. Further east, the rupture traveled along a curved path, which corresponds well with bends seen in the fault mapping of the Sürügü Fault Zone¹⁴. The rupture fronts then terminate at ~65 km east and west of the epicenter, spanning a rupture length of ~130 km, which is also close to the length estimated from field and satellite data¹. The rupture speeds are ~2 km/s eastward and ~3.6 km/s westward. The rupture seems to stagnate in the first 8 s of the event, likely due to the aforementioned signal interference in bilateral ruptures¹³. The westward rupture speed of 3.6 km/s is very close to the local shear wave speed of 3.3–3.5 km/s (refs. ³,¹⁵), and supershear rupture might occur on parts of the west Sürügü Fault. Faster western rupture is reported in the slip model by Melgar et al. ¹⁶, in which a maximum rupture speed of 4.8 km/s is preferred.

**Rupture speed and extent of the mainshock.** The Mw 7.8 event ruptured multiple faults and involved a complex propagation process. A key question is whether the rupture speeds during the mainshock exceed local shear wave speed (Vs). Supershear
ruptures, which propagate faster than Vs, can generate stronger ground motion¹⁷ and cause more severe damage than subshear ones. Supposing a local Vs of 3.39 km/s at 4–12 km depth where most slip occurs³, the overall speeds for the northeast ($V_{NE}$ = 3.05 km/s) and southwest ($V_{SW}$ = 3.11 km/s) are both subshear. Contemporary studies show highly debated rupture speeds (Supplementary Table 2). Mai et al.¹⁸ perform beam-forming BP and report a $V_{NE}$ of 3.1 km/s and a $V_{SW}$ of 2 km/s. Petersen et al.¹⁹ utilize finite source inversion and BP techniques, revealing an average velocity range of 1.8–2.6 km/s and a maximum velocity of 3.4 km/s. Yao & Yang²⁰ determine a $V_{SW}$ range of 3.1–3.4 km/s and a $V_{NE}$ of ~3.19 km/s through near-field data analysis and numerical rupture simulations. Melgar et al.¹⁶ conduct FFI based on high-rate GNSS and local strong-motion data, identifying a preferred rupture speed of 3.2 km/s for both the northeast and southwest branches. Delouis et al.²¹ perform finite source inversion and find that the rupture speeds on the EAF are globally subshear, with several transiently supershear fault portions. Apart from Delouis et al.²¹, supershear ruptures are also declared in some studies, ranging from 5.1 km/s at the splay fault²² to 6 km/s at the southwest EAF²³. Rosakis et al.²² analyze three strong motion stations near the splay fault and identify that the rupture propagates at speed near 1.55 Vs (~5.1 km/s) on the NF after traveling 19.45 km at subshear speeds. Okuwaki et al.²³ perform a potency-density tensor inversion and find a 5–6 km/s speed for the southwest branch of the rupture. Wang et al.²⁴ simulate the dynamic rupture process of the mainshock, and report that $V_{NE}$ is supershear while $V_{SW}$ varies repeatedly between supershear and subshear. Abdelmeguid et al.²⁵ report the occurrence of supershear ruptures at the end of the southwest

segment and the northeast segment of the EAF based on their observation of larger fault-parallel amplitudes compared to fault-normal amplitudes in near-field seismic recordings. However, it should be noted that the interpretation of component amplitudes alone as strong evidence for supershear ruptures is a subject of debate²¹.

We resolve the discrepancy in reported rupture speeds between our results and other studies by showing that Mach wave signatures of supershear rupture are absent in Rayleigh waves recorded globally (Methods). Long and persistent ruptures propagating faster than Rayleigh and Love waves produce surface wave Mach cones with distinguishing signatures in seismograms. At stations located near a pair of azimuths around the rupture direction, surface waves emitted by all portions of the source arrive simultaneously and interfere constructively²⁶,²⁷, leading to waveform shapes highly similar to those of a nearby smaller event with shorter rupture duration (e.g., an M5-6 aftershock or foreshock). The waveform amplitude ratio between the mainshock and the smaller event should equal their seismic moment ratio. Figure 4a, b shows the cross-correlation coefficients (CCs) between the 15–25 s period mainshock waveforms and a nearby M5.3 aftershock (Supplementary Table 3) waveforms, and Fig. 4c, d shows their amplitude ratios over moment ratios. We find that the azimuthal distributions of CCs and amplitude ratios are distinct from those predicted by a supershear rupture. In simulations²⁸ and observations of supershear events²⁶,²⁷, the CC and amplitude ratios feature two maxima at azimuths of 25°–50° relative to the rupture direction and a local minimum at the rupture direction. But here the CC and amplitude ratios peak at the two rupture directions rather than on their sides. This

<figure><img src="image_7.png" /><figcaption>**Fig. 4** Far-field Rayleigh wave analysis. **a** Spatial distribution of stations. Triangles denote the stations used for analysis, color-coded by the cross-correlation coefficients (CCs) between 15 to 25 s Rayleigh waves of the mainshock waveforms and its M5.3 aftershock (Supplementary Table 3) waveforms. The dark green line and word indicate the direction of the northeast rupture. The blue line and word indicate the direction of the southwest direction. **b** CCs distribution as a function of the station azimuth. The vertical black dashed lines indicate the northeast and southwest rupture directions. The pink dash curve indicates the envelope delineating the distribution pattern of CCs. **c** Spatial distribution of stations. Triangles denote the stations used for analysis, color-coded by [waveform amplitude ratio]/[seismic moment ratio] between the mainshock and the aftershock. **d** [Amplitude ratio]/[moment ratio] distribution as a function of the station azimuth.</figcaption></figure>
distribution pattern is consistent with simulation results of a
subshear rupture traveling at Rayleigh wave speed (V<sub>Rayleigh</sub> =
0.92 Vs) (ref. ²⁸). For the southwest rupture, several stations
record highly-similar waveforms (CCs > 0.9; Fig. 4a and Supple-
mentary Fig. 16), indicating a rupture speed around V<sub>Rayleigh</sub>.
Using a local Vs of 3.39 km/s and V<sub>Rayleigh</sub> = 0.92 Vs, the inferred
rupture speed of 3.12 km/s remarkably matches our SEBP results
(3.11 km/s). For the northeast rupture, the maximum CCs are
slightly below 0.9, indicating its speed should be marginally lower
than V<sub>Rayleigh</sub>, consistent with our SEBP results.

The rupture process of the mainshock involves bilateral rupture propagation across multiple phases. To examine the effects of this bilateral rupture and the rupture kink on the cross-correlation coefficient of the empirical Green's function (EGF) and mainshock waveform, we conduct synthetic tests (Methods, Supplementary Figs. 17 and 18). Based on the results of our bilateral rupture synthetic tests, it is observed that when a supershear rupture occurs, there is a distinct peak couple of CCs and amplitude ratios observed at specific azimuths around the strike of the supershear branch (Supplementary Fig. 17c, d, f, g). These findings are consistent with observations from previous supershear earthquakes such as the 2017 Komandorski and the 2016 Romanche events, where similar peak couples were observed within Mach cones²⁸. The synthetic waveform's amplitude contributed by the subshear branch is found to be low, indicating its negligible influence on the observation of Mach waves and Mach cones. The results from the third synthetic case (Supplementary Fig. 17h-j) where rupture speeds are set to be identical to those resolved by our SEBP (V<sub>NE</sub> = 0.9Vs and V<sub>SW</sub> = 0.92Vs) align with the observations: when both V<sub>NE</sub> and V<sub>SW</sub> are close to the Rayleigh wave speed, one peak of CCs and amplitude ratios appears at the northeast rupture strike, and another one occurs at the southwest rupture strike (Fig. 4, Supplementary Fig. 17i, j). The consistency between the observed and synthetic CCs and amplitude ratio distributions further verify our conclusion that both V<sub>NE</sub> and V<sub>SW</sub> are close to the Rayleigh wave speeds and are overall subshear. Additionally, we conduct a synthetic test to simulate a supershear rupture with a 20° kink on the northeast of the EAF (Supplementary Fig. 18). The results indicate that the presence of a small angle kink has little impact on the distribution patterns of CC and amplitude ratio, and the Mach cone pair remains clearly identifiable. The results of our synthetic tests confirm the overall subshear speeds observed by our SEBP and joint FFI. However, we acknowledge the possibility of local and transient supershear behavior, as found by Delouis et al.²¹. For supershear ruptures with lengths comparable to the fault width, Mach waves could not be identified due to the long-period nature of surface waves²⁶.

The small event corresponding to the secondary BP power
peak and occurring around 120 s after the origin time (Fig. 3b)
could be a subevent of the mainshock on the offshore portion of
the Karasu Fault or an immediately triggered aftershock. Similar
short-period emissions are also observed in 108–117 s of BP
results by Petersen et al.¹⁹ and are identified as the last stage of
the mainshock by them. We observe scattered and discontinuous
BP radiators and a minimum BP power during 90–120 s (Figs. 1
and 3a), as well as an offshore aftershock gap near Karasu Fault
and the coast (Fig. 1a). In contrast, many aftershocks concentrate
near the southern end of the Amanos segment. Since aftershocks
tend to happen on the edge of major coseismic slip areas to
release the increased stress caused by the mainshock²⁹,³⁰, we
interpret that the mainshock rupture ended at the Amanos
segment and the offshore event is a triggered aftershock. To
further validate our hypothesis, we calculate the theoretical S
wave travel time from the last High-Freq radiator on the Amanos
segment to the first radiator in the sea (Supplementary Fig. 19)

using the local velocity model from Guvercin et al.³. The
theoretical travel time (30.84 s) matches the time interval between
these two radiators (31.33 s) very well. The time consistency
suggests that the offshore event was dynamically triggered by the
S wave emitted by the stopping phase at the south end of the
mainshock rupture.

**Earthquake supercycle.** During the Mw 7.8 event, most of the slip and seismic moment took place on the Pazarcık, Amanos, and Erkenek segments of the EAF (Supplementary Fig. 15). During the instrumental period between 2007 and 2020, these segments had far less seismic activity than the northeast part of the EAF³. Two major patches with coseismic slip amplitude greater than 8 m on the Pazarcık and Erkenek segments (Fig. 2) coincide with the seismicity voids on them (Supplementary Fig. 20). Ambraseys⁴ investigated the historical earthquakes in SE Turkey and reported temporary seismic quiescence in the 20th century, though violent events occurred within the preceding few centuries. The lack of contemporary seismicity on the main fault contributes to the stress accumulation and gives the main fault potential to host large coseismic slip during earthquakes. Historical documents and field investigations demonstrate that the most recent events on or near the Pazarcık segment occurred in 1114, 1513, and 1795, with magnitudes of 7.8+, 7.4, and 7.0, respectively (Fig. 5b; refs. 4–6). On the Amanos segment, researchers have identified an ~M7 event in 601, an M7.5 event in 1822, and an M7.2 event in 1872 (Fig. 5b) (ref. 5). On the Pazarcık and Amanos segments, the time intervals since the reported last large earthquakes are 218 and 141 years, respectively. The long-term slip rates inferred from GPS data are 7 mm/year and 3 mm/year, respectively³¹,³². The maximum earthquake magnitudes are 7.0–7.3 every 237–772 years (Pazarcık segment) and 7.2–7.4 every 414–915 years (Amanos segment)³.

During the Mw 7.8 earthquake, the equivalent coseismic Mw
on the Pazarcık and Amanos segments are both greater than 7.49
(Supplementary Fig. 15). Such large magnitudes, released within
relatively short time periods since the last large earthquakes
(100–200 years), exceed the expectation based on the character-
istic earthquake cycle model. The pattern can be explained within
the framework of the earthquake supercycle model on segmented
faults, in which a series of events breaking fault segments
separately can be followed by a larger earthquake that breaks
multiple segments at once³³. In the supercycle model, the
magnitude of one earthquake is not only determined by the
cumulative strain since the last earthquake, but also related to the
residual stresses left by earlier events³⁴. Such supercycle models
might also feature clusters of events separated by quiescence
periods³⁵. Supercycle behaviors have been identified on the San
Andreas Fault, the Sumatra subduction zone, the Cascadia
margin, and the Tohoku region³⁴,³⁶⁻³⁹, and could be common,
as suggested by Wallace³⁵.

The latest very large earthquake (M ≥ 7.8) on the EAF before
2023 took place in 1114 (ref. 5), demonstrating a possible
supercycle return period of ~900 years. To evaluate the validity of
the 900-year period, we compute the seismic moment accumula-
tion history since 1114 and check if the balance could be achieved
between build-ups and releases. We first compute the moment
build-up rates on the Amanos, Pazarcık, and Erkenek segments
where coseismic slip occurred during the Mw 7.8 mainshock. The
moment build-up rate: m' = μ·L·H·r, where μ is the rock rigidity
inferred from the 1D model³, L is the segment length, H is the
seismogenic zone depth, and r is the long-term fault slip rates. We
set H as 20 km based on the coseismic slip distribution in the
preferred model and the seismicity distribution presented by
Güvercin et al.³. The long-term slip rates are from Duman and
<figure><img src="image_8.png" /><figcaption>Fig. 5 Historical earthquakes on East Anatolian Fault and San Andreas</figcaption></figure>

**Fault.** a History of potential moment accumulation on the EAF segments involved in the Mw 7.8 mainshock. The blue curve is the moment history of the low accumulation case, and the orange curve is the moment history of the high accumulation case. See text for details. The last big drop is the moment released by the 2023 Mw 7.85 event. b Map of historic earthquakes<sup>4–6,8</sup> and optimal slip orientation. The blue circles and numbers indicate the location, year and magnitude of historic earthquakes close to the faults, and included in the moment history in a. The dark-green circles denote historic events far from the FFI fault planes and are not included in a. The red bars indicate the optimal slip orientation computed based on the regional strain rate field<sup>49</sup> by assuming a rock friction coefficient of 0.6 and the local strain direction being the same as the regional horizontal strain rate direction. The lower inset shows the geometry of San Andreas Fault and San Jacinto Fault and historical events. The black arrows near San Andreas Fault indicate the relative motion of two blocks. Dark blue bars represent the paleoseismic sites of the 1812 earthquake, adapted from Lozos<sup>58</sup>. Dark green bars represent the paleoseismic sites of the latest 1726 earthquake, adapted from refs. 52,53. The pink line denotes the ruptured fault during the 1857 M7.9 event<sup>39</sup>. Red dots and indices denote all M ≥ 6.5 events since 1812. Information of events 1–21 in the lower inset can be found in Supplementary Table 6.

Emre³², and the segment lengths are determined according to our FFI fault planes and fault maps from Duman and Emre³² (Supplementary Table 4).

We assume all moment build-ups are released seismically. For major earthquakes (M ≥ 7), historical investigations⁴,⁵ provide information on their magnitudes, locations, and rupture lengths (Supplementary Table 5). Assuming these catalogs are complete, we can compute the moment release by those major events. As depicted in Fig. 5b and Supplementary Table 5, the events in 1872 (M7.2), 1795 (M7), and 1893 (M7.1) are believed to have occurred on the Amanos, Pazarcık, and Erkenek segments, respectively (refs. ³⁻⁵). The 1513 event (M7.4) is likely to have taken place on a secondary fault rather than the EAF⁴. The occurrence of the 1822 event (M7.5) on the EAF⁴ remains a

subject of debate due to the absence of fresh fault-related topographical features³². On the other hand, it is hard for these studies to include all smaller events (M < 7) before the instrumental period. We utilize the Gutenberg–Richter (G-R) law⁴⁰ to evaluate the seismicity rates of events between Mw 1–7 (excluding 7). These background seismicity rates are influenced by the b values in the G-R relations, where a lower b value corresponds to a higher rate. Güvercin et al.³ reports two sets of b values, differing by approximately 0.2. Thus, we explore two scenarios: a low accumulation case with b values around 0.9 (ref.³), and a high accumulation case with b values around 1.1 (ref.³; see Table 1 and Supplementary Table 4). Concurrently, we include the previously mentioned 1822 event (M7.5) in the low accumulation case but exclude it from the high accumulation case (Table 1). Our calculations show that in the low accumulation scenario, the net moment accumulation between 1114 and 2022 is approximately 32% of the coseismic moment of the 2023 Mw 7.85 mainshock (Fig. 5a), which is equivalent to an Mw 7.52 event. In the high accumulation case, the percentage and equivalent Mw are 90% and 7.82, respectively. Although these estimates are subject to many uncertainties, such as historic earthquake magnitude estimations⁴¹ and temporal variations of geodetic slip rates⁴²,⁴³, it suggests that the 900-year period is close to the low bound of the supercycle duration. Additional estimates of the interseismic coupling rate and the fraction of moment released seismically⁴⁴ are required to explore the upper bound of the super-cycle duration. A contemporary study⁴⁵ also conducted similar moment calculations and estimated cycle lengths, employing lower rates of geodetic strain accumulation and a smaller dataset of historical events. They also observe that the accumulated moment since 1114 is insufficient to account for the Mw 7.8 event in 2023 and attribute this deficit to potential overestimations of historical events.

**Nucleation on splay faults and implication for San Andreas-San Jacinto Fault system.** SEBP and FFI results show that the rupture initiates on the splay Narlı Fault and continues on the main EAF. Similar phenomena have been observed in other fault systems. The 2001 Kokoxili earthquake⁴⁶, the 2002 Mw 7.9 Denali⁴⁷, and the 2016 Mw 7.8 Kaikoura⁴⁸ earthquakes also nucleated on secondary faults and continued onto the main faults. The reason why Mw 7.8 Kahramanmaraş earthquake started on the splay NF deserves investigation. One possibility is that current stress fields make the NF more favorable for earthquake nucleation, which is supported by more abundant seismicity observed on the NF than on its nearby EAF between 2007 and 2020 (Supplementary Fig. 20; ref. ³). To validate the hypothesis, we calculate the optimal orientation for the failure planes on and near EAF (Fig. 5b) based on the strain rate field⁴⁹, assuming a rock friction coefficient of 0.6 and that the strain shares the same direction as the horizontal strain rate. Near the hypocenter, the optimal orientations align better with the strike of the NF than that of the Pazarcık segment. This is consistent with the active seismicity on the NF³. The regional optimal failure orientations also align well with the strike of the Amanos and the Pütürge segments of the EAF, but the former shows quiescent seismicity, and the latter is seismically active. This contrast might be attributed to the discrepancy in long-term slip rates. From the east to the west along the EAF, the slip rates substantially decrease from 10 mm/yr to 1–4 mm/yr.

There is a high resemblance between the San Andreas Fault (SAF) and the EAF systems (Fig. 5b and Supplementary Fig. 21) as similar main fault-splay fault geometry is observed at many sites along the SAF, such as its junction with the San Jacinto Fault (SJF), the San Gabriel Fault, and the San Gregorio Fault
<table><caption>Table 1 The *b* values and moment accumulation rates on the East Anatolian Fault.</caption><thead><tr><th>Scenarios</th><th>Geodetic cumulative rate (Nm/year)</th><th>Background seismicity (Mw 1-7) rate (Nm/year)</th><th>b values in GR relations</th><th>Historic M ≥ 7 events (yyyy, M)</th></tr></thead><tbody><tr><td rowspan="3">Low accumulation</td><td rowspan="3">9.628*10<sup>17</sup></td><td rowspan="3">2.55*10<sup>17</sup></td><td>Erkenek: 0.88</td><td>1795, M7</td></tr><tr><td>Pazarcik: 0.94</td><td>1822, M7.5</td></tr><tr><td>Amanos: 0.92</td><td>1872, M7.2<br>1893, M7.1</td></tr><tr><td rowspan="3">High accumulation</td><td rowspan="3">9.628*10<sup>17</sup></td><td rowspan="3">2.41*10<sup>16</sup></td><td>Erkenek: 1.30</td><td>1795, M7</td></tr><tr><td>Pazarcik: 1.10</td><td>1872, M7.2</td></tr><tr><td>Amanos: 1.10</td><td>1893, M7.1</td></tr></tbody></table>

The geodetic cumulative rate = [long-term slip rate]*[fault length]*[fault width]*[rigidity]. The background seismicity rate is a sum of seismicity moments for Mw 1-7 events. The net moment accumulation rate = [geodetic cumulative rate] – [background seismicity rate].

(Supplementary Fig. 22). Among them, the SJF has the largest slip rate thus should be paid close attention. At Cajon Pass, the N40-45°W oriented SJF branches out from the old and mature SAF, which strikes ~N60°W at this site. Both the SJF and NF are relatively immature compared to the main faults (the SAF and EAF), and feature many stepovers³²,⁵⁰ and dense off-fault seismicity³,⁵¹. The long-term slip rates on the SAF (19–34 mm/ year) and SJF (16 mm/year; ref. ⁴³) are much higher than on the EAF (10 mm/year), and there has been a lack of major earthquakes (M > 7) on the ~500-km-long SAF segment from Parkfield to the Salton Sea since 1857 (the lower inset in Fig. 5b). In the southernmost section of the SAF, specifically near the Salton Sea and Palm Springs, the most recent surface ruptures are attributed to the 1726 earthquake (Fig. 5b; refs. ⁵²,⁵³). On the SJF further south of the 1812 event, the time elapsed since the last great earthquake (Mw > 7, Supplementary Table 6) is also more than 250 years³⁹,⁵⁴,⁵⁵, which is longer than the time interval between the 2023 Mw 7.8 event and the 1800s ones on the EAF. An ~1000 years supercycle on the SAF was proposed based on observations at Wrightwood site, where the strain accumulation is currently approaching the highest level seen in the past 1500 years³⁹. The recent Kahramanmaraş earthquakes, therefore, have important implications for California. Resemblant to the NF in Turkey, the SJF also has very active seismicity⁵¹ and the rupture on SJF could potentially trigger the rupture on the SAF. In fact, the 1812 M7.5 San Juan Capistrano earthquake on the SJF and SAF might to some extent resemble the rupture propagation process of the Mw 7.8 Kahramanmaraş earthquake: surface ruptures in the early 1800s are identified along the SAF and SJF⁵⁶,⁵⁷, and a dynamic rupture model built by Lozos⁵⁸ shows a plausible scenario in which the rupture began near Mystic Lake on the south SJF, and propagated towards the northwest. After arriving at the junction of SJF and SAF, the rupture jumped through the stepover and went further northwest along SAF. If a “Kahramanmaraş earthquake” scenario happens: a rupture starts on the SJF and branches bilaterally on the SAF and releases most moment accumulation at once, the higher slip rate and longer time interval suggest a potentially more severe earthquake in California, with a magnitude of ~8. In-depth investigations into fault interactions during the 2023 Mw 7.8 Kahramanmaraş earthquake sequence will yield valuable insights applicable to the SAF and other major fault systems. One aspect worth exploring is how the rupture on the splay NF triggers the bilateral failure on the EAF. This investigation is particularly relevant because the NF and the southwest EAF have an orientation difference of only approximately 28°. Consequently, slip on the NF is expected to decrease the static Coulomb failure stress on the southwest EAF, creating a “stress shadow” that hinders fault failure⁵⁹. Detailed examinations are needed to understand how the southwest rupture surpasses this stress shadow and its correlation with the background stress conditions. Such inspections will contribute to

a better estimation and mitigation of seismic hazard risks associated with the global strike-slip faults.

## Methods

**SAR data processing.** We apply speckle tracking to the Sentinel-1 SAR images for the surface deformation in slant range (radar line-of-sight, or cross-track) and ground azimuth (along-track) directions via the InSAR Scientific Computing Environment (ISCE) version 2 software⁶⁰ using the “topsApp.py” application⁶¹ (Supplementary Fig. 6a–f). The two level 1 Single Look Complex (SLC) data are coregistered with each other using the geometrical co-registration with refinement from the enhanced spectral diversity technique⁶². Then we applied the speckle tracking (also known as pixel tracking, amplitude cross-correlation) analysis to pairs of Sentinel-1 SLC images (detailed pairs information are listed in Supplementary Table 7) to estimate the coseismic surface deformation. We oversampled the coregistered SLC by a factor of 32, used a cross-correlation window or chip with a relatively large size of 128 pixels in range direction and 64 pixels in azimuth direction (about 585 ×900 m on the ground) to reduce the estimation noise⁶³. To further reduce the offset estimation noise, we applied a median filter with a ground size of 1.5 ×2.2 km; masked out pixels on the water body or with abnormal offset value (>5 m); for range offsets masked out pixels with offset standard deviation (STD) > 0.5 m, and for azimuth offsets masked out pixels with offset STD > 3.0 m. Then we geocoded all offsets to a grid in size of 0.00277778° (~300 m) in latitude and longitude.

We apply InSAR and MAI to the ALOS-2 ScanSAR images for the surface deformation in slant range and ground azimuth directions via the ISCE-2 software⁶⁰ using the “alos2App.py” and “alos2burstApp.py” applications⁶⁴⁻⁶⁶ (Supplementary Fig. 6g–j); and apply InSAR to the LuTan-1 stripmap images for the surface deformation in the slant range direction via a combination of Gamma⁶⁷ and ISCE-2⁶⁰ software packages (Supplementary Fig. 6k). For all interferograms from ALOS-2 and LuTan-1, we use the SRTM DEM⁶⁸ (SRTMGL1, 1 arc second with void-filled, ~30 m), apply the Goldstein filter with a strength of 0.8, unwrap the interferogram using the minimum cost flow method⁶⁹ with near-fault pixels masked out using the surface rupture traces from Sentinel-1 speckle tracking¹. For ALOS-2 InSAR, we multilook the interferogram by 5 and 28 looks in range and azimuth directions, respectively; correct for ionospheric delay using the range split-spectrum technique⁶⁶; geocode into a grid in size of 0.000833334° (~90 m; Supplementary Fig. 23a, d). For LuTan-1 InSAR, we coregister the single look complex (SLC) images using DEM and precise orbits with refinements from cross correlation, multilook the interferogram by 24 and 28 looks in range and azimuth directions, respectively; correct for phase unwrapping errors by adding an integer number of phase cycles (2π) to one side of the Mw 7.8 rupture using the ALOS-2 descending
interferogram as reference; geocode to a grid in size of 0.00055556° (~60 m; Supplementary Fig. 24). The ionospheric delay is not corrected due to the expected weak impact at 6:30 am local time. For ALOS-2 MAI, we extract burst SLC images from ALOS-2 full aperture ScanSAR products. The burst SLCs are used to form burst interferograms. We then compute the MAI interferograms using the burst interferograms. We multilook the MAI interferograms by 14 and 4 looks in range and azimuth directions, respectively. Subsequent processing is the same as regular InSAR processing. The ALOS-2 ScanSAR system is a four-look burst system, so we can get three MAI interferograms with burst separations one, two and three. The final MAI result is a weighted mean of the three MAI interferograms<sup>65</sup>, and geocoded to a grid in size of 0.000833334° (~90 m).

For all measurements in the range direction (Sentinel-1 range offsets, ALOS-2 and LuTan-1 interferograms), we correct for tropospheric delay using the ERA5 global atmospheric model<sup>70</sup> via PyAPS<sup>71</sup> (Supplementary Fig. 23b, e), for solid Earth tides following the 2010 IERS convention<sup>72</sup> via PySolid<sup>73</sup> (Supplementary Fig. 23c, f).

We combine all observations above to derive the three dimensional co-seismic deformation (Supplementary Fig. 7), including three range offsets and three azimuth offsets from Sentinel-1, two regular interferograms and two MAI interferogram from ALOS-2, one regular interferogram from LuTan-1. Built upson Fialko et al. (ref. 74), we use a weighted approach to balance the different uncertainties among speckle tracking, InSAR and MAI techniques, and propagate the uncertainties into the east, north and up components for each pixel. More specifically, we use the pixelwised speckle tracking STD for Sentinel-1, a constant STD of 0.32 and 0.30 m for ALOS-2 ascending and descending MAI (estimated from the non-deforming region in the far-field) to account for the ionospheric impact, and a constant STD of 0.05 m for ALOS-2 and LuTan-1 interferograms to account for residual tropospheric delays. The estimated east, north and up components are masked using a STD threshold of 0.4, 0.8 and 0.3 m, respectively.

**Joint Finite Fault Inversion.** The Finite Fault Inversion (FFI) method was developed after the 1979 Imperial Valley earthquake<sup>75,76</sup> to study the rupture complexities of earthquakes. With decades of development, current FFI algorithms could include global and local data to simultaneously determine the slip amplitude, rupture time, and rise time<sup>11,77,78</sup>. The method has been widely applied to large earthquakes and provides important information for hazard assessment and mitigation (e.g., 2004 Sumatra earthquake<sup>79</sup>; 2011 Tohoku earthquake<sup>80</sup>).

We use the joint FFI to image the fault's rupture process and slip distribution<sup>11,12</sup>. The seismic waveform and geodetic data are adopted to constrain the slip model. The seismic waveform data include 20 broadband teleseismic P waves (bandpass filtered to 2–333 s, Supplementary Fig. 1), 13 broadband teleseismic S waves (bandpass filtered to 2–333 s, Supplementary Fig. 1), 23 long-period teleseismic Rayleigh waves (bandpass filtered to 166–250 s, Supplementary Fig. 2), 17 long-period teleseismic Love waves (bandpass filtered to 166–250 s, Supplementary Fig. 2), 57 near-field strong motion recordings (P, SH, SV for 19 stations, bandpass filtered to 4–40 s and integrate to velocity seismograms, Supplementary Figs. 3 and 5). The geodetic data include 24 high-rate GNSS recordings (UD, EW, NS for 8 stations, displacement waveforms, Supplementary Figs. 4 and 5), static ground displacements derived from Sentinel-1 SAR data, ALOS-2 and LuTan-1 InSAR and MAI data (Fig. 1a, b, Supplementary Figs. 6–8; see SAR Data Processing for more details). The combination of multiple datasets could provide better constraints

on the slip amplitude, slip distribution, rupture speed, and rise time than any single dataset<sup>81,82</sup>. The seismic and high-rate GNSS data are transformed into wavelet domains. Compared with conventional inversion only in frequency or time domain, the wavelet transform allows us to simultaneously capture the characteristics in the high-frequency and low-frequency and preserve the information in the time-domain<sup>11,12</sup>. The misfits for static ground deformation measured by SAR, InSAR and MAI are calculated separately since they could not be included in the wavelet transform. We use the weighted sum-squared residuals to measure the difference between observed and synthetic data. The weights are inversely proportional to the standard deviation at each data point. We identify the fault surface traces from the SAR images. Six vertical fault planes (segments S1–S6 in Fig. 1b) are set to model the rupture. We parameterize the fault planes with 45 and 7 subfaults in the along strike and dip directions, respectively. We adopt the Laplacian smoothing constraints for fault slip<sup>11</sup> and perturbations of rupture initiation time at individual subfaults<sup>83</sup>. The subfault size is 9 km (along strike) by 3 km (downdip; Supplementary Fig. 12). The hypocenter is set at 37.218°N, 37.007°E, 7.5 km depth (adapted from a relocated hypocenter shared by Sezim Ezgi Güvercin on Twitter). The rupture velocity is allowed between 1.5 km/s to 4.9 km/s with a reference speed of 3.0 km/s, and the rake angle is allowed ±30° from the reference angle of 11° provided by Global Centroid Moment Tensor Catalog. The Green's functions for local strong motion waveforms and static displacements are computed using the f-k integration approach by Zhu and Rivera<sup>84</sup> and the 1-D crust model by Güvercin et al.<sup>3</sup>. Subfaults are interpolated to 1 km × 1 km when plotting, to improve the visualization (Fig. 2 and Supplementary Fig. 13).

**Slowness-Enhanced Back-Projection (SEBP) and validation.** A method that allows for clear source images to be obtained as soon as data is available and without accurate fault parameters is back-projection (BP). This source imaging technique tracks the propagation of earthquake rupture fronts based on coherent seismic recordings from dense networks of seismometers<sup>85</sup>. Results of BP show the spatiotemporal evolution of ~1 Hz radiators and are useful in obtaining rupture parameters such as length and speed<sup>86,87</sup>.

In our study, we use the Multitaper-MULTiple SIgnal Classification (MUSIC) type of BP enhanced by the use of 'reference windows'. This removes an apparent drift in the radiators towards stations in the array that would otherwise be displayed in the image results, allowing for sharper images to be obtained<sup>88</sup>. In standard BP, 1-dimensional velocity models (e.g., IASP91, ref. 89) are used to obtain travel times of P-waves from the mainshock hypocenter to station locations. In reality, these travel times deviate from those predicted using a 1D model due to ray path effects caused by Earth's 3D structure. To account for this, BP corrects the travel times using the hypocenter alignment method<sup>86,90</sup>. This method assumes that the first arrival on observed seismograms comes from the hypocenter location obtained by a local catalog (Disaster and Emergency Management Authority of Turkey; AFAD for our application) and then cross-correlates the first 8 s of the P-wave recordings. The key assumption of this standard approach is that these travel time errors are uniform over the whole rupture area, which is only valid for locations close to the hypocenter.

To accurately account for travel time errors further away from the hypocenter, the BP method we use utilizes the spatial derivatives of travel times (slowness) with respect to source locations for a given station, and is hereafter referred to as Slowness-Enhanced back-projection (SEBP<sup>10</sup>). This additional
correction term is computed using locations of aftershocks within the rupture area. For each aftershock, the travel times relative to the hypocenter are compared with the predictions without the slowness correction term. Then, any differences are mapped to the slowness correction term by dividing them by the distance between the mainshock and the aftershock. Several studies have shown that SEBP provides more precise rupture location and speed estimates¹⁰,²⁷,⁹¹,⁹².

We collect the vertical components of the P wave broadband seismograms from two large-aperture arrays at teleseismic distance, located in China (CH) and Alaska (AK, the lower inset in Fig. 1a). The regular-size AK array contains ~200 stations, and the frequency band with high resolution and enough signal-to-noise ratio (SNR) is 0.5–2 Hz. The vast CH array includes more than 700 stations available for BP, which ensure high SNR in 1–4 Hz BP and provide better resolution than regular-size arrays. For the SEBP of the Mw 7.8 event, 718 seismograms from the CH array in 1–4 Hz and 182 seismograms from the AK array in 0.5–2 Hz are included. For the SEBP of the Mw 7.5 event, we process 719 and 153 seismograms from the CH and AK arrays, respectively. We carried out SEBP with a 10-s-long sliding window with a 1-s-long step. Slowness correction terms are derived using 11 and 7 nearby aftershocks for the CH and AK applications to the Mw 7.8 event, with 6 and 5 aftershocks for the CH and AK array applications to the Mw 7.5 event, respectively. The mainshock and aftershocks’ hypocenter locations and magnitudes based on the AFAD catalog are shown in Supplementary Tables 1a–d, along with comparisons of our aftershock locations before and after slowness calibration (Supplementary Figs. 9 and 10). The root-mean-square distance error relative to the AFAD locations decreases by at least 61% after calibration for all our applications.

The dense network of accelerometer stations along the Pazarlık and Amanos segments (white triangles in Fig. 1a) provides us with a good benchmark to validate the SEBP results. We collect recordings from 33 near-field accelerometer stations. After detrending and demeaning accelerometer recordings, we bandpass-filter them using a frequency band of 1–4 Hz to match the BP band of the CH array application. We then calculate and smooth the envelopes of each recording using the Hilbert transform method. We do not include several northern stations on our plot in Fig. 3c since, at the time of this study, they did not have available data in the time window of interest. We trace S-phases radiated from the fixed High-Freq radiators obtained from CH array BP results (Fig. 1a) to test whether the arrival times of these S-phases at station locations coincide with time frames of high energy release. We use the TauP toolkit for ray tracing, where we assumed a 1-D velocity model from Guvercin et al.³. Since that velocity model does not provide a Moho depth estimate, we adopt a Moho depth of 37 km from the CRUST1.0 global velocity model at the location of the mainshock. The High-Freq radiators resolved by CH array BP are also plotted in Fig. 3c according to their along-fault distance and time.

**Mach wave searching and synthetic tests.** For ruptures propagating at speeds lower than Rayleigh and Love waves speeds, surface waves generated by all portions of the source arrive at far-field stations at different time. If the rupture speeds (Vr) increase to faster than Rayleigh and Love speeds, at stations located near a pair of azimuths around the rupture direction, surface waves emitted by all portions of the source arrive simultaneously and interfere constructively²⁶. The azimuth pair is called Mach cones, and the surface waves recorded within Mach cones are Mach waves. Due to the simultaneous arrival within Mach cones, the waves from all sub-sources interfere constructively.

Theory predicts that the shapes of bandpassed seismograms should be identical or highly similar to those from a smaller reference event (hereinafter referred to as EGF event) with analogous focal mechanisms²⁶. The EGF event usually has a magnitude of 5–6. The amplitude ratios of waveforms between the large supershear and EGF events equal their seismic moment ratios. Since the Rayleigh and Love waves have speeds very close to the shear wave speeds ($V_{\text{Rayleigh}} = 0.92$ Vs), Mach waves and Mach cones are also believed to be unique features of supershear ruptures.

The similarities between waveforms are measured by standard cross-correlation coefficients (CCs): $CC(U, u) = \frac{\sum_{i=1}^{N} U_i u_i}{\sqrt{\sum_{i=1}^{N} U_i U_i \sum_{i=1}^{N} u_i u_i}}$, where U and u are the main-shock and EGF event displacement waveforms, respectively. We use $\Phi$ to denote the angle between the rupture direction and the source-station azimuth. Theory predicts Mach cones at $\Phi_M$: $\Phi_M = \arccos(V_{\text{Rayleigh}}/Vr)$, where $V_{\text{Rayleigh}}$ is the Rayleigh wave speed and $Vr$ is the rupture speed. The Mach angle is then positively related to $Vr$. For the rupture with a speed of $V_{\text{Rayleigh}}$, the Mach cone appears in the frontal direction of the rupture with $\Phi_M=0$. When $Vr$ increases to $1.5V_{\text{Rayleigh}}$ (1.38Vs), $\Phi_M$ increases to 49°.

We collect 46 vertical components of broadband seismic recordings (BHZ). To limit the effect of surface wave dispersion, seismograms are bandpass filtered to 15–25 s. The period is shorter than the mainshock duration but longer than the EGF duration, keeping the point source effect of the EGF event and the finite source effect of the mainshock rupture. In actual observation, due to the dispersion of surface waves and the source directivity effect, highly-similar waveforms can be found in an azimuth range around $\Phi_M$. The apparent source duration is $T_a = T_0 \times D(\Phi)$, where $D(\Phi) = 1-\cos(\Phi) \times Vr/V_{\text{Rayleigh}}$. Similar waveforms appear in angle ranges satisfying $T_a \ll T_c$, where $T_c$ is the center period of waves and is the 20 s in our case. Considering the possible travel time errors caused by the 3D path effect, we manually pick the EGF Rayleigh wave arrivals. We calculate the CCs in the 200–300 s window to include the entire envelope of EGF Rayleigh waveforms. Apart from CCs, we also compute the ratio: $R = R_{\text{amp}}/R_{\text{mom}}$, where $R_{\text{amp}}$ is the amplitude ratio between the waveforms of the mainshock and EGF event, $R_{\text{mom}}$ is the seismic moment ratio between the mainshock and EGF event.

To better understand how the distribution patterns of CC and R vary with rupture speeds, Bao et al.²⁸ performed synthetic tests to explore scenarios with different rupture speeds and present CC and R's distribution as a function of $\Phi$. When $Vr \approx V_{\text{Rayleigh}}$, CC and R are maximized in the rupture direction. For supershear ruptures, CC and R reach their maximum on the two Mach cones, whose $\Phi_M$ varies from 25° to 50° depending on the rupture speeds. The 2023 Mw 7.8 Kahramanmaraş earthquake ruptured bilaterally, with the northeast branch striking ~N40°E and the southwest branch striking ~S20°W. The maximum peaks of CC and R appear near these two directions (Fig. 4b, d), indicating that both branches should propagate by speeds close to $V_{\text{Rayleigh}}$ (see the main text for more detailed discussions).

Vallée and Dunham²⁶ derive the theory and present observations of Mach waves based on the assumption of a long unilateral rupture, where the rupture length is substantially greater than the fault width. Moreover, the selected stations should be in far-field, where the station epicenter distance is much larger than the rupture length. To meet the far-field assumption, 44 of our stations are located beyond 20° (>2200 km), with 2 closer stations (6–8°, approximately 660–880 km) included as well. In a compound rupture comprising segments with different rupture speeds and directions, waveforms from subshear segments interfere destructively at stations along the Mach cone, contributing minimally to the overall recordings²⁶. It is confirmed by
the observations that Mach waves from long and sustained supershear segments can be distinctly identified in many large strike-slip events involving multiple rupture phases. These events include the following examples: the 2001 Mw 7.8 Kokoxili earthquake, where a supershear segment is enclosed between two subshear segments²⁶; the 2016 Mw 7.1 Romanche earthquake, characterized by a westward supershear back-propagation following an eastward subshear rupture⁹³; and the 2021 Mw 7.4 Maduo earthquake, exhibiting simultaneous westward subshear and eastward supershear propagations⁹⁴.

The rupture process of the 2023 Mw 7.8 Kahramanmaraş earthquake is highly complex, involving bilateral rupture propagation across multiple phases. To examine the effects of this bilateral rupture and the rupture kink on the cross-correlation coefficient of the empirical Green's function (EGF) and mainshock waveform, we conduct synthetic tests. In these tests, the synthetic mainshock rupture propagation is simulated as a series of sub-sources. The waveforms of the synthetic mainshock rupture are modeled as the superposition of EGF waveforms, considering the rupture times and locations of different sub-sources. The slip distribution on the East Anatolian Fault is used to establish a bilateral rupture model (Supplementary Fig. 17a). The northeast rupture has a strike of N50°E and a length of 140 km, while the southwest rupture has a strike of S20°W and a length of 200 km. The interval between sub-sources is 9 km, evenly distributed along the rupture path. We tested three rupture scenarios: (1) northeast rupture speed $V_{NE} = 1.2$ Vs, southwest rupture speed $V_{SW} = 0.7$ Vs; (2) $V_{NE} = 0.7$ Vs, $V_{SW} = 1.2$ Vs; (3) $V_{NE} = 0.9$ Vs, $V_{SW} = 0.92$ Vs. Case (3) represents the rupture speeds resolved by our SEBP and joint FFI (as discussed in the main text). The synthetic waveforms and EGF waveforms are filtered to 15–25 s, and we assume a Rayleigh wave speed of 0.92 Vs at the center frequency (1/20 Hz). Cross-correlation coefficients are then measured between the synthetic waveforms and EGF waveforms (Supplementary Fig. 17c, f, i), as well as the amplitude ratios of the synthetic waveforms to the seismic moment ratios (Amp_ratio/Mom_ratio, Supplementary Fig. 17d, g, j). We also conduct a synthetic test to simulate a 140-km-long supershear rupture ($V_r = 1.2$ Vs) with a 20° kink in the middle (Supplementary Fig. 18). This case simulates the kinked fault we observe on the northeast East Anatolian Fault.

**Moment accumulation history.** We include the Amanos, Pazarcik, Erkenek segments of the EAF, where the coseismic slip of the Mw 7.8 quake occurred, in the computation. The moment accumulation history consists of two parts: the moment build-up and moment release. The build-up rate:

$m' = \mu \cdot L \cdot H \cdot r$, where $\mu$ is the rock rigidity, $L$ is the segment length, $H$ is the seismogenic zone depth, and $r$ is the long-term fault slip rates. We set $H$ as 20 km based on the coseismic slip distribution in the preferred model (Fig. 2) and the seismicity distribution presented by Güvercin et al.³. The long-term slip rates are from Duman and Emre³². The segment lengths are determined according to our FFI fault planes and fault maps in Duman and Emre³² (parameters are listed in Supplementary Table 4).

We assume all moment build-ups are released seismically and divide the release into two parts. The first part is the moment unloading by major events (M ≥ 7). Historical investigations⁴⁻⁶ provide information about their magnitudes, epicenter locations, and rupture lengths. Assuming these catalogs are complete, we can compute the moment release by major events. However, it is hard for these studies to include all smaller events (M < 7) before the instrumental period. We utilize the Gutenberg–Richter (G-R) law⁴⁴ to evaluate the seismicity rates of events between Mw 1 and

7. The G-R law states that the number (N) of earthquakes with a magnitude larger than M over a certain period on a given fault obeys:

$log_{10}(N) = a - b \cdot M$, where a and b are G-R law constants determined empirically by regional seismicity data.

Güvercin et al.³ investigated the seismicity of the EAF between 2007 and 2020. They evaluated the a and b values for each segment of EAF (listed in Supplementary Table 4). With these values, we can calculate the seismicity rate (R) for magnitudes 1–7 on each segment. The annual seismic moment (*Mo')* of earthquakes with magnitude M is then:

$$Mo' = R \cdot 10^{(1.5M+9.1)}$$

Summing *Mo'* for M1–7 earthquakes, we get their annual moment release rate on the fault. The final moment accumulation history is the moment build-up minus moment release (Fig. 5a).

## Data availability

The Copernicus Sentinel-1 data are provided by the European Space Agency and obtained from the Alaska Satellite Facility (https://search.asf.alaska.edu/) via Seamless SAR Archive (SSARA), a service provided by the EarthScope Consortium. The ALOS-2 data are provided by the Japan Aerospace Exploration Agency through Sentinel Asia and the International Disaster Charter. The LuTan-1 data are provided by the China Centre for Resources Satellite Data and Application. The fault rupture traces are obtained from USGS¹. The high-rate GNSS data are provided by TUSAGA-Aktif (https://www.tusaga-aktif.gov.tr/). The strong motion data are available at the Disaster and Emergency Management Authority of Turkey (AFAD; https://tadas.afad.gov.tr/event-detail/15499). The moment tensor solutions come from the Global Centroid Moment Tensor project (CMT; http://www.globalcmt.org). The seismic data are provided by the IRIS (https://ds.iris.edu/wilber3/) and the Data Management Center of China National Seismic Network at Institute of Geophysics, China Earthquake Administration (http://www.seisdmc.ac.cn). The local earthquake catalog is available at AFAD (https://deprem.afad.gov.tr/event-catalog). The SAR data used in inversions, the station list used in Mach wave searching, SEBP results, and the FFI slip model are available at: https://doi.org/10.5281/zenodo.8402190.

## Code availability

The SAR processing is performed using the ISCE-2 software⁶⁰, available at https://github.com/isce-framework/isce2 with postprocessing recipe available at https://doi.org/10.5281/zenodo.8098260. The MATLAB code of SEBP is available at https://github.com/Ismeng/MUSICBP/tree/SEBP. The Python code of Mach wave analysis is available at https://github.com/Ismeng/MUSICBP/tree/MachWave. The finite fault inversion codes used in this study are available upon request.

Received: 29 June 2023; Accepted: 29 September 2023; Published online: 17 October 2023

## References

1. Reitman, N. G. et al., Preliminary fault rupture mapping of the 2023 M7.8 and M7.5 Türkiye Earthquakes. Turkey Earthquake Emergency Response (2023). https://doi.org/10.5066/P98517U2. Accessed 16 May 2023.

2. Gocümlü, B. Ç. “Sağlık Bakanı Koca: 10 ilde 17 bin 929'u hekim olmak üzere 143 bin 829 personelimiz hizmet veriyor” [Health Minister Husband: 143 thousand 829 personnel, 17 thousand 929 of whom are physicians, provide service in 10 provinces] (in Turkish). Anadolu Agency. Retrieved 10 February (2023).

3. Güvercin, S. E., Karabulut, H., Konca, A. Ö., Doğan, U. & Ergintav, S. Active seismotectonics of the East Anatolian Fault. *Geophys. J. Int.* **230**, 50–69 (2022).

4. Ambraseys, N. Temporary seismic quiescence: SE Turkey. *Geophys. J. Int.* **96**, 311–331 (1989).

5. Ambraseys, N. & Jackson, J. Faulting associated with historical and recent earthquakes in the Eastern Mediterranean region. *Geophys. J. Int.* **133**, 390–406 (1998).

6. Palutoglt, M. & Sasmaz, A. 29 November 1795 Kahramanmaras Earthquake. *Southern Turkey. Bull. Mineral Res. Explor.* **155**, 187–202 (2017).

7. Westaway, R. 1994. Present-day kinematics of the Middle East and eastern Mediterranean. *J. Geophys. Res.* **99**, 12071–12090 (1994).

8. Chen, K., Zhang, Z., Liang, C., Xue, C. & Liu, P. Kinematics and dynamics of the 24 January 2020 Mw6.7 elazig, Turkey earthquake. *Earth Space Sci.* **7**, 11 (2020).
9. Lin, X. et al. Coseismic slip distribution of the 24 January 2020 Mw 6.7 Doganyol Earthquake and in relation to the foreshock and aftershock activities. *Seismol. Res. Lett.* **92**, 127–139 (2020).

10. Meng, L., Zhang, A. & Yagi, Y. Improving back projection imaging with a novel physics-based aftershock calibration approach: a case study of the 2015 Gorkha earthquake. *Geophys. Res. Lett.* **43**, 628–636 (2016).

11. Ji, C., Wald, D. J. & Helmberger, D. V. Source description of the 1999 Hector Mine, California, Earthquake, Part I: wavelet domain inversion theory and resolution analysis. *Bull. Seism. Soc. Am.* **92**, 1192–1207 (2002).

12. Ji, C., Helmberger, D. V., Wald, D. J. & Ma, K.-F. Slip history and dynamic implications of the 1999 Chi-Chi, Taiwan, earthquake. *J. Geophys. Res. Solid Earth.* **108**, 2412 (2003).

13. Xu, L. et al. Understanding the rupture kinematics and slip model of the 2021 Mw 7.4 Maduo earthquake: a bilateral event on bifurcating faults. *J. Geophys. Res. Solid Earth.* **128**, e2022J025936 (2023).

14. Emre, Ö. et al. Active Fault Map of Turkey with an Explanatory Text 1:250,000 Scale. *Special Publication Series.* **30** (General Directorate of Mineral Research and Exploration, 2013).

15. Acarel, D., Cambaz, M. D., Turhan, F., Mutlu, A. K. & Polat, R. Seismotectonics of Malatya Fault, Eastern Turkey. *Open Geosci.* **11**, 1098–1111 (2019).

16. Melgar, D. et al. Sub- and super-shear ruptures during the 2023 Mw 7.8 and Mw 7.6 earthquake doublet in SE Türkiye. *Seismica* **2**, 3 (2023).

17. Dunham, E. M. & Bhat, H. S. Attenuation of radiated ground motion and stresses from three-dimensional supershear ruptures. *J. Geophys. Res. Solid Earth.* **113**, B08319 (2008).

18. Mai, P. M. et al. The destructive earthquake doublet of 6 February 2023 in South-Central Türkiye and Northwestern Syria: initial observations and analyses. *Seism. Record* **3**, 105–115 (2023).

19. Petersen, G. M. et al. The 2023 Southeast Türkiye seismic sequence: rupture of a complex fault network. *Seism. Record* **3**, 134–143 (2023).

20. Yao, S., & Yang, H. Rupture phase in near-fault records of the 2023 Turkey Mw 7.8 earthquake. *EarthArXiv* https://doi.org/10.31223/X51662 (2023).

21. Delouis, B., van den Ende, M., & Ampuero, J.-P. Kinematic rupture model of the February 6th 2023 Mw 7.8 Turkey earthquake from a large set of near-source strong motion records combined by GNSS offsets reveals intermittent supershear rupture. *Authorea* https://doi.org/10.22541/essoar.168286647.71550161/v1 (2023).

22. Rosakis, A., Abdelmeguid, M., & Elbanna, A. Evidence of Early Supershear Transition in the Mw 7.8 Kahramanmaraş Earthquake From Near-Field Records. Preprint at https://arxiv.org/abs/2302.07214 (2023).

23. Okuwaki, R., Yagi, Y., Taymaz, T. & Hicks, S. P. Multi-scale rupture growth with alternating directions in a complex fault network during the 2023 south-eastern Türkiye and Syria earthquake doublet. *Geophys. Res. Lett.* **50**, e2023GL0103480 (2023).

24. Wang, Z. et al. Dynamic rupture process of the Mw 7.8 Kahramanmaraş earthquake (SE Türkiye): variable rupture speed and implications for seismic hazard. *Authorea* https://doi.org/10.22541/essoar.168614604.40356164/v1 (2023).

25. Abdelmeguid, M. et al. Revealing the dynamics of the Feb 6th 2023 M7.8 Kahramanmaracs/Pazarıkçı Earthquake: near-field records and dynamic rupture modeling. *arXiv* https://doi.org/10.48550/arXiv.2305.01825 (2023).

26. Vallée, M. & Dunham, E. M. Observation of far-field Mach waves generated by the 2001 Kokoxili supershear earthquake. *Geophys. Res. Lett.* **39**, L05311 (2012).

27. Bao, H. et al. Early and persistent supershear rupture of the 2018 magnitude 7.5 Palu earthquake. *Nat. Geosci.* **12**, 200–205 (2019).

28. Bao, H. et al. Global frequency of oceanic and continental supershear earthquakes. *Nat. Geosci.* **15**, 942–949 (2022).

29. Hartzell, S. H. & Heaton, T. H. Rupture history of the 1984 Morgan Hill, California, earthquake from the inversion of strong motion records. *Bull. Seism. Soc. Am.* **76**, 649–674 (1986).

30. Mendoza, C. & Hartzell, S. H. Aftershock patterns and main shock faulting. *Bull. Seism. Soc. Am.* **78**, 1438–1449 (1988).

31. Reilinger, R. et al. GPS constraints on continental deformation in the Africa-Arabia-Eurasia continental collision zone and implications for the dynamics of plate interactions. *J. geophys. Res.* **111**, B05411 (2006).

32. Duman, T. Y. & Emre, Ö. The East Anatolian Fault: geometry, segmentation and jog characteristics. *Geol. Soc. Spec. Publ.* **372**, 495–529 (2013).

33. Philibosian, B. & Meltzner, A. J. Segmentation and supercycles: a catalog of earthquake rupture patterns from the Sumatran Sunda Megahrust and other well-studied faults worldwide. *Quat. Sci. Rev.* **241**, 106390 (2020).

34. Salditch, L. et al. Earthquake supercycles and long-term fault memory. *Tectonophysics* **774**, 228298 (2020).

35. Wallace, R. E. Grouping and migration of surface faulting and variation in slip rates on faults in the Great Basin province. *Bull. Seismol. Soc. Am.* **77**, 868–877 (1987).

36. Goldfinger, C., Ikeda, Y., Yeats, R. S. & Ren, J. Superquakes and supercycles. *Seismol. Res. Lett.* **84**, 24–32 (2013).

37. Satake, K. Geological and historical evidence of irregular recurrent earthquakes in Japan. *Phil. Trans. R. Soc. A* **373**, 20140375 (2015).

38. Sieh, K. et al. Earthquake supercycles inferred from sea-level changes recorded in the corals of West Sumatra. *Science.* **322**, 1674–1678 (2008).

39. Weldon, R., Scharer, K., Fumal, T. & Biasi, G. Wrightwood and the earthquake cycle: what a long recurrence record tells us about how faults work. *GSA Today.* **14**, 4–10 (2004).

40. Gutenberg, B. & Richter, C. F. Frequency of earthquakes in California. *Bull. Seism. Soc. Am.* **34**, 185–188 (1944).

41. Kokum, M. & Özcelik, F. An example study on re-evaluation of historical earthquakes: 1789 Palu (Elazığ) earthquake, Eastern Anatolia, Turkey. *Bull. Min. Res. Explor.* **161**, 157–170 (2020).

42. Savage, J. & Prescott, W. Asthenosphere readjustment and the earthquake cycle. *J. Geophys. Res. Solid Earth.* **83**, 3369–3376 (1978).

43. Tong, X., Smith-Konter, B. & Sandwell, D. T. Is there a discrepancy between geological and geodetic slip rates along the San Andreas Fault System? *J. Geophys. Res. Solid Earth.* **119**, 2518–2538 (2014).

44. Avouac, J. P. From geodetic imaging of seismic and aseismic fault slip to dynamic modeling of the seismic cycle. *Ann. Rev. Earth Planet. Sci.* **43**, 233–271 (2015).

45. Karabulut, H., Güvercin, S. E., Hollingsworth, J. & Konca, A. Ö. Long silence on the East Anatolian Fault Zone (Southern Turkey) ends with devastating double earthquakes (6 February 2023) over a seismic gap: implications for the seismic potential in the Eastern Mediterranean region. *J. Geol. Soc.* **180**, jgs2023–021 (2023).

46. Lasserre, C. G. et al. Coseismic deformation of the 2001 Mw = 7.8 Kokoxili earthquake in Tibet, measured by synthetic aperture radar interferometry. *J. Geophys. Res.* **110**, B12408 (2005).

47. Eberhart-Phillips, D. et al. The 2002 Denali Fault earthquake, Alaska: a large magnitude, slip-partitioned event. *Science.* **300**, 1113–1118 (2003).

48. Cesca, S. et al. Complex rupture process of the Mw 7.8, 2016, Kaikoura earthquake, New Zealand, and its aftershock sequence. *Earth Planet. Sci. Lett.* **478**, 110–120 (2017).

49. Weiss, J. R. et al. High-resolution surface velocities and strain for Anatolia from Sentinel-1 InSAR and GNSS data. *Geophys. Res. Lett.* **47**, e2020GL087376 (2020).

50. Wesnousky, S. Seismological and structural evolution of strike-slip faults. *Nature.* **335**, 340–343 (1988).

51. Ross, Z. E., Hauksson, E. & Ben-Zion, Y. Abundant off-fault seismicity and orthogonal structures in the San Jacinto fault zone. *Sci. Adv.* **3**, e1601946 (2017).

52. Scharer, K. M. & Yule, D. A maximum rupture model for the Southern San Andreas and San Jacinto Faults, California, derived from paleoseismic earthquake ages: observations and limitations. *Geophys. Res. Lett.* **47**, e2020GL088532 (2020).

53. Rockwell, T. K., Meltzner, A. J. & Haaker, E. C. Dates of the two most recent surface ruptures on the Southernmost San Andreas fault recalculated by precise dating of Lake Cahuilla dry periods. *Bull. Seism. Soc. Am.* **108**, 2634–2649 (2018).

54. Fialko, Y. Interseismic strain accumulation and the earthquake potential on the southern San Andreas fault system. *Nature.* **441**, 968–971 (2006).

55. Olsen, K. B. et al. Strong shaking in Los Angeles expected from southern San Andreas earthquake. *Geophys. Res. Lett.* **33**, L07305 (2006).

56. Seitz, G., Weldon, R. II & Biasi, G. P. The Pitman Canyon paleoseismic record: a re-evaluation of southern San Andreas Fault segmentation. *J. Geodyn.* **24**, 129–138 (1997).

57. Kendrick, K. J. & Fumal, T. E. Paleoseismicity of the northern San Jacinto fault, Colton and San Bernardino, southern California; preliminary results. *Geol. Soc. Am. Abstr. Prog.* **37**, 559 (2005).

58. Lozos, J. C. A case for historic joint rupture of the San Andreas and San Jacinto faults. *Sci. Adv.* **2**, e1500621 (2016).

59. Yamashita, T. & Umeda, Y. Earthquake rupture complexity due to dynamic nucleation and interaction of subsidiary faults. *Pure Appl. Geophys.* **143**, 89–116 (1994).

60. Rosen, P. A., Gurrola, E., Sacco, G. F., & Zebker, H. The InSAR scientific computing environment. In *EUSAR 2012; 9th European Conference on Synthetic Aperture Radar*. Nuremberg, Germany, 23-26 April 2012, 730–733 (VDE, 2012).

61. Fattahi, H., Agram, P. & Simons, M. A network-based enhanced spectral diversity approach for TOPS time-series analysis. *IEEE Trans. Geosci. Remote Sens.* **55**, 777–786 (2017).

62. Yagüe-Martinez, N. et al. Interferometric processing of sentinel-1 TOPS data. *IEEE Trans. Geosci. Remote Sens.* **54**, 2220–2234 (2016).

63. De Zan, F. Accuracy of incoherent speckle tracking for circular Gaussian signals. *IEEE Geosci. Remote Sens. Lett.* **11**, 264–267 (2014).

64. Liang, C. & Fielding, E. J. Interferometry with ALOS-2 full-aperture ScanSAR data. *IEEE Trans. Geosci. Remote Sens.* **55**, 2739–2750 (2017a).

65. Liang, C. & Fielding, E. J. Measuring azimuth deformation with L-band ALOS-2 ScanSAR interferometry. *IEEE Trans. Geosci. Remote Sens.* **55**, 2725–2738 (2017b).
66. Liang, C., Liu, Z., Fielding, E. J. & Bürgmann, R. InSAR time series analysis of L-band wide-swath SAR data acquired by ALOS-2. *IEEE Trans. Geosci. Remote Sens.* **56**, 4492–4506 (2018).

67. Werner, C., Wegmüller, U., Strozzi, T., & Wiesmann, A. Gamma SAR and interferometric processing software. *Proceedings of the ERS-Envisat symposium*. (Gothenburg, Sweden, 2000).

68. Farr, T. G. et al. The shuttle radar topography mission. *Rev. Geophys.* **45**, 2 (2007).

69. Chen, C. W. & Zebker, H. A. Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization. *J. Opt. Soc. Am. A* **18**, 338–351 (2001).

70. Hersbach, H. et al. The ERA5 global reanalysis. *Q. J. R. Meteorol. Soc.* **146**, 1999–2049 (2020).

71. Jolivet, R., Grandin, R., Lasserre, C., Doin, M. P. & Peltzer, G. Systematic InSAR tropospheric phase delay corrections from global meteorological reanalysis data. *Geophys. Res. Lett.* **38**, L17311 (2011).

72. Petit, G., & Luzum, B. IERS Conventions (IERS Technical Note No. 36). Retrieved from Frankfurt, Germany: http://www.iers.org/TN36/.

73. Yunjun, Z. et al. Range geolocation accuracy of C-/L-band SAR and its implications for operational stack coregistration. *IEEE Trans. Geosci. Remote Sens.* **60**, 5227219 (2022).

74. Fialko, Y., Simons, M. & Agnew, D. The complete (3-D) surface displacement field in the epicentral area of the 1999 MW7.1 Hector Mine Earthquake, California, from space geodetic observations. *Geophys. Res. Lett.* **28**, 3063–3066 (2001).

75. Olson, A. H., & Apsel, R. Finite fault and inversion theory with applications to the 1979 Imperial Valley earthquake, *Bull. Seism. Soc. Am.* **72**, 1969–2001 (1982).

76. Hartzell, S. H. & Heaton, T. H. Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California earthquake. *Bull. Seism. Soc Am.* **73**, 1553–1583 (1983).

77. Hartzell, S. H., Liu, P. & Mendoza, C. The 1994 Northridge, California, earthquake: investigation of rupture velocity, rise time, and high-frequency radiation. *J. Geophys. Res. Solid Earth.* **101**, 20091–20108 (1996).

78. Yue, H. et al. Depth varying rupture properties during the 2015 Mw 7.8 Gorkha (Nepal) earthquake. *Tectonophysics.* **714**, 44–54 (2017).

79. Ammon, C. J. et al. Rupture process of the 2004 Sumatra-Andaman Earthquake. *Science* **308**, 1133–1139 (2005).

80. Minson, S. E. et al. Bayesian inversion for finite fault earthquake source models - II: the 2011 great Tohoku-oki, Japan earthquake. *Geophys. J. Int.* **198**, 922–940 (2014).

81. Kanamori, H. & Given, J. W. Use of long-period surface waves for rapid determination of earthquake-source parameters. *Phys. Earth Planet. Inter.* **27**, 8–31 (1981).

82. Wald, D. J. & Graves, R. W. Resolution analysis of finite fault source inversion using 1D and 3D Green's functions. II Combining seismic and geodetic data. *J. Geophys. Res. Solid Earth.* **106**, 8767–8788 (2001).

83. Shao, G. et al. Focal mechanism and slip history of the 2011 Mw 9.1 off the Pacific coast of Tohoku Earthquake, constrained with teleseismic body and surface waves. *Earth Planet. Sp.* **63**, 559–564 (2011).

84. Zhu, L. & Rivera, L. A. A note on the dynamic and static displacements from a point source in multilayered media. *Geophys. J. Int.* **148**, 619–627 (2002).

85. Kiser, E. & Ishii, M. Back-projection imaging of earthquakes. *Annu. Rev. Earth Planet. Sci.* **45**, 1 (2017).

86. Ishii, M. et al. Extent, duration and speed of the 2004 Sumatra-Andaman earthquake imaged by the Hi-Net array. *Nature.* **435**, 933–936 (2005).

87. Kiser, E. & Ishii, M. The 2010 Mw 8.8 Chile earthquake: triggering on multiple segments and frequency-dependent rupture behavior. *Geophys. Res. Lett.* **38**, L07301 (2011).

88. Meng, L., Inbal, A. & Ampuero, J. P. A window into the complexity of the dynamic rupture of the 2011 Mw 9 Tohoku-Oki earthquake. *Geophys. Res. Lett.* **38**, L00G07 (2011).

89. Kennet, B. L. N. I. A. S. P. E. I. 1991 seismological tables. *Terra Nova.* **3**, 122 (1991).

90. Ishii, M., Shearer, P. M., Houston, H. & Vidale, J. E. Teleseismic P Wave imaging of the 26 December 2004 Sumatra-Andaman and 28 March 2005 Sumatra earthquake ruptures using the Hi-net array. *J. Geophys. Res.* **112**, B11307 (2007).

91. Meng, L. et al. Double pincer movement: encircling rupture splitting during the 2015 Mw 8.3 Ilapel earthquake. *Earth Planet. Sci. Lett.* **495**, 164–173 (2018).

92. Meng, L., Huang, H., Xie, Y., Bao, H. & Dominguez, L. A. Nucleation and kinematic rupture of the 2017 Mw 8.2 Tehuantepec earthquake. *Geophys. Res. Lett.* **46**, 3745–3754 (2019).

93. Hicks, S. P. et al. Back-propagating supershear rupture in the 2016 Mw 7.1 Romanche transform fault earthquake. *Nat. Geosci.* **13**, 647–653 (2020).

94. Zhang, X., Feng, W., Du, H., Samsonov, S. & Yi, L. Supershear rupture during the 2021 MW 7.4 Maduo, China, earthquake. *Geophys. Res. Lett.* **49**, e2022GL097984 (2022).

95. Mahmoud, Y. et al. Kinematic study at the junction of the East Anatolian fault and the Dead Sea fault from GPS measurements. *J Geodyn.* **67**, 30–39 (2013).

96. Styron, R. & Pagani, M. The GEM global active faults database. *Earthq. Spectra.* **36**, 160–180 (2020).

Acknowledgements

We thank Gilles Peltzer for valuable comments and discussions, Qingyue Yang, Junfeng Li and Yachao Wang for help with the LuTan-1 data processing. We thank the editor and anonymous reviewers for helping improve the quality of the manuscript. We thank the Disaster and Emergency Management Authority of Turkey (AFAD) for providing strong motion data and local earthquake catalog. We thank TUSAGA-Aktif for providing high-rate GNSS data. L.X., S.M., and L.M. were supported by NSF grant no. EAR-1848486 and Leon and Joanne V.C. Knopoff Fund. C.J. was supported by grants from the Southern California Earthquake Center (SCEC)-funded by NSF cooperative agreement EAR-0109624 and USGS cooperative agreement 02HQAG0008. J.-P.A. was supported by the French government through the UCAJEDI Investments in the Future project (ANR-15-IDEX-01) managed by the French National Research Agency (ANR). R.C. was supported by NSFC grants U20A2095 and 42325401. C.L. was supported by NSFC grant 42274026. Figures were produced using Generic Mapping Tools (GMT), MATLAB, and Python. The Python software package ObsPy was used for data requesting, waveform filtering, and CCs calculation.

Author contributions

L.X.: performed the SEBP, the joint FFI, and the Rayleigh Mach wave analysis; designed the figures and tables; participated in the interpretation of the results; wrote the original draft. S.M.: performed the SEBP and the SEBP validation; designed the figures and tables; participated in the interpretation of the results; wrote the original draft. L.M.: conceived and led the study; performed the SEBP and the SEBP validation; designed the figures and tables; participated in the interpretation of the results; wrote the original draft. C.J.: conceived and led the study; performed the joint FFI; designed the figures and tables; participated in the interpretation of the results; wrote the original draft. J.-P.A.: performed the SEBP validation; participated in the interpretation of the results; wrote the original draft. Z.Y.: processed the SAR images; wrote the original draft. M.H.: performed the SEBP. R.C.: processed the China array seismic data. C.L.: processed the SAR images. All authors contributed to finalize the manuscript.

Competing interests

The authors declare no competing interests.

Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-023-01030-x.

**Correspondence** and requests for materials should be addressed to Lingsen Meng or Chen Ji.

**Peer review information** *Communications Earth & Environment* thanks the anonymous reviewers for their contribution to the peer review of this work. Primary Handling Editor: Joe Aslin. A peer review file is available

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023