Complex multi-fault rupture and triggering during the 2023 earthquake doublet in southeastern Türkiye

The Anatolian microplate is being extruded westward between the converging Arabian, African, and Eurasian plates, producing major strike-slip systems: the right-lateral North Anatolian Fault Zone (NAFZ) and the left-lateral East Anatolian Fault Zone (EAFZ). The EAFZ comprises geometrically complex, segmented strike-slip faults that bifurcate into a northern strand and a main strand, with variable slip rates (~10 mm/yr in the northeast to ~4 mm/yr in the southwest). Although less active than the NAFZ during the 20th century, geodetic coupling indicated accumulated strain capable of large events, which was realized on February 6, 2023, in a destructive doublet (Mw 7.8 and Mw 7.7). This study aims to resolve the detailed kinematic rupture processes on the complex fault network using joint inversion of dense seismic and geodetic data and to assess static stress transfer and potential triggering between the two events and onto neighboring faults. Understanding these processes is crucial for seismic hazard assessment in southeastern Türkiye and adjacent regions.

Prior work documents the geometry, segmentation, and jog characteristics of the EAFZ, including conjugate fractures, pull-apart basins, bends, and stepovers that control rupture size and occurrence. Historical large events occurred in 1513, 1795, and multiple earthquakes between 1822 and 1905. Recent geodetic studies quantified interseismic coupling and slip rates along the EAFZ. The 2020 Doğanyol-Sivrice Mw 6.7 event ruptured the central EAFZ main strand and provided context for strain release patterns. Early studies of the 2023 sequence mapped complex surface ruptures from satellite data and produced initial finite-fault models, identifying mixed sub- and supershear rupture phases. The literature also notes that supershear ruptures often exhibit paucity of aftershocks in high-speed segments and tend to radiate lower moment-scaled energy compared to the global mean for strike-slip events.

- Data sets: Joint inversion of multiple data types for both events: regional strong-motion accelerograms (corrected to obtain ground velocities, displacements, and static offsets), static GNSS coseismic offsets, high-rate GNSS time series, and teleseismic P and SH broadband waveforms with good azimuthal coverage at 30°–90° distances. - Strong-motion processing: Introduced an updated baseline correction scheme improving on a bi-linear approach by using an iteratively smoothed natural curve correction. Steps included: pre-seismic baseline removal, integration to velocity, estimation of post-seismic linear trend, construction of a starting correction curve, iterative smoothing under constraints, final correction curve determination, and baseline correction to obtain velocity and displacement. Applied to 52 stations (Mw 7.8) and 26 stations (Mw 7.7), yielding stable coseismic displacements at 21 and 4 near-fault stations, respectively. Regional waveforms bandpass filtered 0.02–0.5 Hz and sampled at 0.2 s; 300 s time windows; onsets hand-picked. - GNSS data: High-rate PRIDE PPP-AR solutions for selected stations (0.2 s sampling; 300 s windows); static offsets from JPL rapid-orbit 5-min series at 29 (Mw 7.8) and 7 (Mw 7.7) GNSS sites; vertical static included only at EKZ1 for Mw 7.7 due to precision limits. - Fault geometry: Multi-segment fault models constrained by post-event satellite rupture mapping and relocated aftershocks. Mw 7.8: six segments (aF1–aF6) spanning the main EAFZ and an initiating splay fault. Mw 7.7: five segments (bF1–bF5) along the northern strand with curved western extent. Segment parameters in Supplementary Tables. - Inversion: Nonlinear finite-fault inversion in the wavelet domain simultaneously fitting seismic waveforms (L1+L2 norm across wavelets) and static displacements (sum-squared residuals). Green’s functions computed with a regional 1D velocity model. Initial models randomized with total moment fixed to GCMT. Weighting: static vs waveform errors equal; strong-motion-derived statics weighted at half of GNSS statics. Multiple simulated annealing runs (10 seeds per event) to assess model stability and uncertainty. - Validation and fits: Compared synthetic to observed strong-motion, high-rate GNSS, static GNSS, and teleseismic waveforms; assessed residuals and noted misfits possibly due to 1D Earth model and site effects. - Coulomb stress: Computed ΔCFS = Δτ + μΔσn with μ = 0.4 using PSGRN/PSCMP. Tested multiple receiver fault geometries for the Mw 7.7 initial segment (USGS, AFAD, GCMT, and this study) at 10 km depth to assess mainshock-induced loading and potential triggering. Calculated regional ΔCFS for three target areas (A1: NE EAFZ; A2: DSF; A3: diffuse faults west of source) combining stress from the 2023 doublet and the 2020 Mw 6.7 event. - Teleseismic radiated energy: Retrieved from IRIS EQEnergy and normalized by seismic moment to evaluate E_r/M_0 relative to global strike-slip events.

- Near-fault static displacements: Strong-motion derived coseismic offsets robustly constrained the slip, with largest horizontal permanent displacement ~2.8 m at station 4614 and vertical ~0.6 m at station 4615 for the Mw 7.8 event. Results generally agree with Sentinel-1 pixel tracking, with some near-fault differences due to methodological uncertainties. - Mw 7.8 Pazarcık earthquake (main EAFZ): - Nucleation on a short, previously unmapped splay fault, then cascading onto the steeply dipping main strand, propagating bilaterally (~160 km NE, ~180 km SW; total ~340 km). - Predominantly left-lateral strike-slip with minor normal/thrust components; strong spatial heterogeneity. - Seismic moment M0 ≈ 7.1 × 10^20 N·m (Mw 7.82); rupture duration ~90 s; peak slip ~8.1 m near the splay-main strand intersection. - Rupture velocities: NE average ~3.5 km/s (locally up to ~4.5 km/s between 30–40 s, supershear) and SW average ~2.5–3.2 km/s with a supershear phase up to ~3.8 km/s between 55–70 s. Supershear segments coincide with aftershock paucity. - Mw 7.7 Ekinözü earthquake (northern strand of EAFZ): - Bilateral rupture along E–W striking, curved segments with complex SW geometry; predominantly strike-slip with significant shallow slip. - Largest slip concentrated on bF1–bF2; complementary aftershock patterns (dense where slip is low; sparse in large-slip/supershear zones). - Maximum slip ~11 m (bF1); rupture duration ~65 s; M0 ≈ 5.0 × 10^20 N·m (Mw 7.7). - Early (<8 s) bilateral supershear ~4.0 km/s followed by subshear speeds: NE ~3.0 km/s (~0.88 Vs), SW ~2.7 km/s; SW deceleration linked to fault discontinuity. - Triggering and stress transfer: - Mainshock-induced ΔCFS at the Mw 7.7 source increased by ~0.014–0.189 MPa across four tested receiver geometries at 10 km depth, exceeding typical ~0.01 MPa triggering threshold, consistent with delayed triggering of the back-branch rupture. - Combined ΔCFS from the 2023 doublet and 2020 Mw 6.7 event shows increased loading (up to ~0.1 MPa) on: (A1) northeastern EAFZ; (A2) Dead Sea Fault south of the Amanos segment; and positive changes in (A3) diffuse faults west of the northern strand. - Radiated energy: - Mw 7.8: E_r ≈ 1.36 × 10^16 J; Mw 7.7: E_r ≈ 7.29 × 10^15 J; moment-scaled E_r/M_0 ≈ 1.9 × 10^−5 and 1.5 × 10^−5, respectively, both lower than the global mean for large strike-slip events, consistent with tendencies for supershear cases. - Damage implications: Strong directivity along the SW Amanos fault with average ~3.2 km/s during Mw 7.8 likely amplified shaking toward western Syria, contributing to severe damage despite slip being within Türkiye.

The study resolves how complex multi-segment ruptures with variable rupture velocities, including supershear phases, developed during the 2023 doublet. The detailed kinematic models, constrained by near-fault static displacements and diverse seismic/geodetic data, demonstrate that the Mw 7.8 mainshock’s rapid, patchy slip and supershear stages produced strong directivity and stress redistribution favorable to triggering. The static ΔCFS analysis indicates that the mainshock likely advanced failure on the receiver geometry of the Mw 7.7 event and increased loading on adjacent fault systems (NE EAFZ and DSF), offering a plausible mechanistic link for the 9-hour delayed doublet. Variations in rupture speed and slip heterogeneity correlate with geometric complexities (bends, stepovers, branching) and potential lateral variations in stress and frictional properties. The findings align with other studies identifying supershear phases in this sequence but further resolve the distinct temporal staging of supershear for each event. The relatively low E_r/M_0 values compared to global strike-slip averages are consistent with smoother high-speed propagation on straight segments, whereas rougher slip would radiate more short-period energy. Overall, the results elucidate how geometric complexity, stress heterogeneity, and prior loading set the stage for cascading, partially supershear ruptures with significant hazard implications.

By integrating dense strong-motion, GNSS (static and high-rate), and teleseismic data with satellite-constrained fault geometries and relocated aftershocks, the study establishes robust kinematic rupture models for the 2023 Türkiye earthquake doublet. Key contributions include: (1) development and application of an enhanced strong-motion baseline correction enabling reliable near-fault static offsets; (2) resolution of complex multi-fault cascading ruptures with pronounced bilateral propagation and temporally distinct supershear phases; (3) quantification of rupture speeds, peak slips, durations, and spatial slip heterogeneity; and (4) demonstration of positive ΔCFS at the Mw 7.7 source and increased loading on neighboring faults, supporting a triggering relationship and indicating future seismic potential on the NE EAFZ, DSF, and diffuse western faults. Future work should refine 3D structural models and fault geometries (e.g., dips, bends, offshore segments), expand very-near-fault instrumentation to reduce model uncertainty, integrate site effects for improved waveform fits at higher frequencies, and pursue dynamic rupture modeling to test physical mechanisms governing supershear onset and cascading across complex networks.

- Earth model simplification: Green’s functions computed with a regional 1D velocity model likely contribute to waveform misfits, particularly regionally where 3D structure and site effects are important. - Model geometry uncertainties: Segment dips and detailed geometries vary among studies; uncertainties in receiver fault parameters affect ΔCFS estimates. - Data limitations: High-frequency content of strong-motion records is not fully explained; limited very near-fault observations in some areas increase slip-model standard deviations (notably on aF3, aF6, bF1, bF2). - Inversion non-uniqueness: Simulated annealing shows sensitivity to random seeds when multiple near-optimal solutions exist; though large-slip features are stable, local variations persist. - Static offsets: Strong-motion-derived statics have inherent uncertainties due to baseline correction, reflected in reduced weighting relative to GNSS statics.