Super-shear ruptures steered by pre-stress heterogeneities during the 2023 Kahramanmaraş earthquake doublet

The East Anatolian Fault (EAF) is a major intra-continental transform system with a history of destructive earthquakes. On 6 February 2023, an M7.8 earthquake on a splay of the EAF was followed ~9 hours later by an M7.5 event on the Sürgü–Çardak (S-C) fault ~90 km away. The sequence caused >50,700 fatalities across southeastern Turkey and northern Syria. Despite the EAF’s complex geometry of bends, step-overs, and sub-parallel strands, the doublet produced large, energetic ruptures. Prior studies disagree on details of rupture evolution and speeds (including possible supershear phases), particularly the transition from the initial splay rupture to propagation on the main EAF and the M7.5 rupture kinematics. The central questions are how such energetic ruptures developed on a geometrically complex system and what role ambient stress heterogeneity played in steering rupture propagation.

Previous analyses indicate the M7.8 rupture began on the Nurdağı–Pazarcık Fault (NPF) splay and then engaged the main EAF; some studies argue for an immediate bilateral rupture, while back-projection suggests a 40–50 s northeastward propagation before bilateral growth. Reported rupture speeds range from subshear (≈2.0–3.2 km/s) to possible supershear. For the M7.5 event, published estimates indicate supershear velocities from ~4.5 km/s up to ~6 km/s, though many studies lacked near-field GNSS constraints at key stations. Dynamic rupture simulations have reproduced aspects of the sequence by invoking along-strike stress rotations. These discrepancies motivate a re-evaluation using combined geodetic, strong-motion, and remote-sensing constraints and explicit estimation of pre-stress fields.

- Geodetic and remote-sensing data: Computed 2D horizontal deformation from Sentinel-1 C-band SAR and Sentinel-2 optical image pixel offsets via precise co-registration and sub-pixel correlation; inverted six independent look directions (four radar, two optical) to obtain 3D surface deformation. Validated against GNSS displacements. - Surface rupture and slip-vector extraction: Determined surface rupture traces and local strike from Rupture mapping and swath-profile analysis of 3D displacement discontinuities; used only horizontal components due to vertical uncertainty. Projected into fault-parallel and fault-normal directions; inverted fault-parallel profiles using a sum of linear and error functions to estimate surface slip and uncertainties. - Fault segmentation: Approximated the ruptures with 9 planar segments for M7.8 (A1–A9) and 6 for M7.5 (B1–B6) based on surface ruptures and deformation fields. - Bayesian inversion of geometry: Estimated subsurface dip angles, widths, and uniform strike-/dip-slip per segment using ALOS-2 ScanSAR LOS displacements (ascending/descending) and GNSS static offsets. Applied parallel sequential MCMC (1000 chains) to sample posterior distributions; dip angles for A8 (fixed 73°) and A9 (fixed vertical) constrained by aftershocks; other segment dips free within 30°–90°. Data down-sampled via quadtree; accounted for covariance. Median of posterior samples taken as preferred model; uncertainties from posterior standard deviations. - Finite-source inversion: Jointly inverted high-rate (1 Hz) GNSS, strong-motion velocity waveforms, and GNSS static offsets using a multi-time-window approach. Segments discretized into ~5×5 km subfaults (468 for M7.8; 222 for M7.5), six rows along dip; variable number along strike. Allowed rake variations (±25° around GCMT for M7.8; two slip vectors −80° and 10° for M7.5 to capture normal component). Source time function: five overlapping symmetric triangles (5 s each, 2.5 s overlap). Rupture velocities determined via grid search informed by back-projection: M7.8 split into three velocity groups (NPF A8; NE A5,A6,A7,A9; SW A1–A4). M7.5 split into supershear-capable group (B2–B4) and other (B1,B5,B6). - Waveform processing: GNSS waveforms (17 stations) and static offsets (26 stations for M7.8; 14 for M7.5) processed via PPP; filtered (GNSS 0.005–0.4 Hz; strong motion decimated to 1 Hz and filtered 0.02–0.2 Hz); 90 s windows from origin time. Green’s functions computed with layered velocity model via frequency–wavenumber integration; same filters applied. Equal weighting across data types after tests; first-order Laplacian regularization for stability. Jackknife tests by randomly removing 20% of data to assess robustness. - Pre-stress estimation: Inverted 55 pre-event focal mechanisms (2007–2020) for deviatoric stress tensors in zones along the EAF (southern A1–A3, central A4–A5,A8, northern A6–A7,A9) and separately for S-C. Compared with fault-instability metric and slip vectors from kinematic model. Also inferred stress orientation for S-C from regional strain-rate field (InSAR+GNSS, 2014–2019) assuming proportionality of stress and strain rate.

- Geometry and dips: Bayesian inversion indicates M7.8 segments are near-vertical; M7.5 segments dip more gently, especially B1, B5, B6 (~50–55°). Posterior medians align with aftershock distributions; model fits ALOS-2 LOS and GNSS data well. - M7.8 kinematics: Best-fitting models require a ~10 s delay between initial NE propagation and onset of bilateral rupture. Optimal rupture velocities are fast (3.0–4.0 km/s); improvement beyond 4.0 km/s marginal, and transient supershear cannot be resolved by grid search. Total seismic moment ≈7.76×10^20 Nm (Mw 7.86); duration ~75 s with peak moment rate around 25 s; maximum slip ~11 m; maximum slip rate ~3.0 m/s at ~40 km from epicenter. Most slip shallower than 15 km; A2 notably shallower than adjacent segments. Nucleation on NPF contributes a small fraction of total moment. Rupture propagated past a ~40° bend near the EAF–DSF intersection at high speed (A4 ~3.8 km/s, likely > Rayleigh speed). - M7.5 kinematics: Inclusion of near-field GNSS (e.g., EKZ1) enables resolving bilateral supershear over the initial ~40 km, contradicting prior unilateral-west supershear interpretations. Optimal rupture speeds reach ~5.0–6.0 km/s on B2–B4 and ~2.8 km/s on B1,B5,B6. Total seismic moment ≈5.57×10^20 Nm (Mw 7.76); duration ~35 s; peak moment rate about twice that of M7.8. Maximum slip ~11 m; slip on B3–B4 is deeper; normal-slip component becomes more prominent toward the NW termination. Estimated stress drop is large (~11–21 MPa) for the compact, high-slip rupture. - Pre-stress alignment and rotation: Inverted pre-stress tensors show a ~20° anticlockwise rotation of σ1 from SW to NE along the M7.8 rupture. Fault surfaces are nearly optimally oriented to local stresses (fault instability I > 0.7; mean ~0.98 decreasing to ~0.85 NE), explaining sustained fast rupture despite large along-strike changes in fault strike. On segment A4, σ1 is ~73° to the EAF, coinciding with fast rupture and conditions favoring dynamically increased shear/Coulomb stresses along large branch angles on the extensional side. - Limited stress transfer to S-C: Static Coulomb changes from the M7.8 event on the S-C faults were small (up to about −1.5 MPa; <1 MPa at and east of the hypocenter), and dynamic stresses moderate (up to ~7 MPa only east of the hypocenter), insufficient to explain the large, fast M7.5 rupture. - Local stress heterogeneity for M7.5: Focal-mechanism inversion along S-C suggests misalignment (σ1 ~N30°E; ~60° to S-C), but is poorly constrained (only six mechanisms). Strain-rate-derived σ1 indicates a ~30° rotation relative to NE EAF, yielding σ1 ~N47.5°E and ~34° misfit to S-C along B2–B4, placing most ruptured fault patches in high-instability regimes (I > 0.75). This local stress state can explain prolonged, supershear rupture without requiring exceptionally low friction. Observed normal slip on a N–S, west-dipping fault near the NW termination is consistent with off-fault tensional stresses induced by S-C slip. - Overall: Heterogeneous ambient pre-stress steered rupture paths and speeds in both events, enabling sustained fast and supershear propagation on geometrically complex and locally misoriented segments.

The findings address how large, energetic ruptures developed on a complex fault system. For the M7.8 event, an along-strike rotation of σ1 renders otherwise variably striking segments nearly optimally oriented, enabling fast bilateral rupture and allowing propagation through a ~40° bend near the EAF–DSF intersection. Locally high rupture speed and highly oblique σ1 (e.g., A4) amplify dynamic shear and Coulomb stresses at large branch angles on the extensional side, sustaining rupture along the bending fault. For the M7.5 event, the apparently misoriented S-C fault system can host supershear propagation if the local pre-stress substantially deviated from the regional state. Strain-rate-derived stress orientations suggest many S-C segments were actually well oriented for failure, consistent with the observed bilateral supershear and large stress drop. The limited static and moderate dynamic stress transfer from the M7.8 event cannot alone account for the M7.5 rupture, implying that pre-existing local stress heterogeneity was critical. Additional mechanisms—free-surface-induced normal stress changes on dipping strike-slip faults, bimaterial interface effects, and depth-variable friction—may further facilitate supershear transitions and high slip near the surface. Overall, integrating kinematic inversions with spatially variable pre-stress estimates reconciles previously conflicting observations and highlights the importance of local stress heterogeneity in controlling rupture dynamics on complex fault networks.

This study combines surface-rupture mapping, SAR/optical geodesy, high-rate GNSS, and strong-motion data to resolve the geometry and kinematics of the 2023 Kahramanmaraş doublet. The M7.8 rupture propagated rapidly (3.0–4.0 km/s) following a brief delay to bilateral growth, traversing a large bend owing to along-strike σ1 rotation and dynamic stress amplification. The M7.5 rupture exhibited bilateral supershear (5.0–6.0 km/s) over ~80 km, with a large stress drop and deep slip concentrations, consistent with a locally rotated pre-stress field inferred from strain-rate data. These results demonstrate that spatial heterogeneity of ambient pre-stress can steer rupture trajectories and speeds, allowing fast or supershear propagation on complex and locally misaligned faults. Future work should (i) test whether the inferred stress rotations are required to reproduce the observed kinematics via dynamic rupture simulations, (ii) refine data weighting and the role of near-field stations in inversions, and (iii) integrate strain-rate fields to estimate stress orientations where seismicity is sparse.

- Near-field station weighting and its influence on rupture-speed resolution warrant further exploration. - Bayesian geometry inversion did not use the 3D deformation maps due to lower accuracy, relying on ALOS-2 LOS and GNSS data; A8 and A9 dips were fixed from aftershocks due to weak surface constraints. - Sparse focal mechanisms along the Sürgü–Çardak faults (only six) limit robustness of stress inversions based solely on seismicity there; strain-rate-based stress inferences assume proportionality between stress and strain rate. - Finite-source grid-search inversions cannot resolve transient supershear phases; resolution of slip is better above ~15 km depth per jackknife tests. - Complex geology and material property variations introduce unavoidable model uncertainties; further dynamic simulations are needed to confirm necessity of the proposed pre-stress distributions.