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

The 2023 M7.8 and M7.5 earthquake doublet near Kahramanmaraş, Turkey, provides insight regarding how large earthquakes rupture complex faults. Here we determine the faults geometry using surface ruptures and Synthetic Aperture Radar measurements, and the rupture kinematics from the joint inversion of high-rate Global Navigation Satellite System (GNSS), strong-motion waveforms, and GNSS static displacement. The M7.8 event initiated on a splay fault and subsequently propagated along the main East Anatolian Fault with an average rupture velocity between 3.0 and 4.0 km/s. In contrast, the M7.5 event demonstrated a bilateral supershear rupture of about 5.0–6.0 km/s over an 80 km length. Despite varying strike and dip angles, the sub-faults involved in the mainshock are nearly optimally oriented relative to the local stress tensor. The second event ruptured a fault misaligned with respect to the regional stress, also hinting at the effect of local stress heterogeneity in addition to a possible free surface effect.

To determine the 2D horizontal deformation field, we used a subpixel correlation technique on optical images acquired before and after the M7.8 and M7.5 earthquakes.These orthorectified images were obtained from the European Space Agency's Sentinel-2 satellite sensor that acquires images at 10 m pixel resolution.Specifically, the deformation maps were estimated from the images using the COSI-Corr software which uses a phase correlation scheme to estimate the sub-pixel shift of a subset of pixels under the assumption that a translation in the spatial domain is equivalent to a shift in phase 1,2 .Here we used a correlation window with a dimension of 64×64 pixels and a step size of 16 pixels.This produces a 2D deformation map of the horizontal surface motion with a uniform pixel resolution of 160 m.
We used right-looking C-band Sentinel-1 SAR images acquired before and after the M7.8 and M7.5 events to estimate radar amplitude pixel offsets.In total we used data from three tracks, two ascending (Track 14 and 116) and one descending (Track 21) that were acquired in TOPS mode.Together these tracks cover the entire rupture constraining surface motion in the range and azimuthal (i.e., in-flight) directions.To estimate the pixels offsets we used the InSAR Scientific Computing Environment (ISCE-2) developed at the Jet Propulsion Lab to first co-register the reference and secondary single-look complex images (SLCs) 3 .For cross-correlation we used a window of 64×64 pixels, a search distance of 20 pixels and a skip size of 32 pixels, using the denseOffsets.pyroutine in ISCE.Following correlation, the pixel offsets are geocoded from radar to geographic co-ordinates.The variability for each of the radar offset results is listed in Supplementary Table 1, which is estimated from the variability in surface motion of a far-field stable region.

Supplementray Note 2. 3D surface deformation from inversion of radar and optical pixel offsets
To align the various datasets into the same co-ordinate system, the offsets are taken relative to a far-field stable point that is common to all tracks (37.0209°E, 37.0172° N) instead of using GNSS as there is insufficient station coverage located across each track.This gives a surface deformation result that is in a reference frame relative to the stable region.
The 3D surface deformation () is then estimated at each pixel from inverting the satellite pixel offsets () with at least three independent look directions using a weighted least squares (WLS) approach,  = &′) !"  ′  .The diagonal weighting matrix () contains values that are the inverse of the variances of each dataset (with 1 uncertainty up to 0.2 m and 0.5 m for Sentinel-1 range and azimuths respectively, while we find uncertainty of 1 = 0.25 m and 0.31 cm for Sentinel-2 east-west and north-south components respectively, see Supplementary Table 1 for details).Measurements of the surface fault slip show very strong statistical agreement with measurements from other studies using similar pixel tracking geodetic imaging data (Supplementary Figure 3).

Supplementray Note 3. Prestress tensor inversion
The 3D deviatoric stress tensor is estimated by inverting the focal mechanisms from Güvercin et al. 4 under the Wallace-Bott assumption that slip is parallel to the shear stress 5 .The stress inversion gives the orientation and shape of the 3D deviatoric stress tensor but not its magnitude.The principal deviatoric stresses  " ,  # ,  $ are ordered from most to least compressive.To resolve spatial variations of stress along the rupture we distinguish three zones along the M7.8 rupture (Supplementary Figure 21), making sure that each contains sufficient diversity of fault orientation to resolve the stress tensor according to the criteria given by Hardebeck and Hauksson 6 .Where for all zones the angle of the slip vectors have an RMS > 30° from the average (see Supplementary Figure 22 which shows the distributions of each zone).To determine the nodal plane for our inversion for the stress tensors, we follow the approach of Vavryčuk 7 by using an iterative approach to pick nodal planes that are optimally aligned to the stress field that have the highest fault instability value.
To invert the unit slip vectors for stress we use a L2 least squares inversion 8 .To minimize overfitting of the data and to constrain stress to be spatially smooth along the rupture we use a damping constraint to the inversion that penalizes large changes of the stress orientation between neighboring zones (i.e., the gradients of the model vector between cells) 9 .The uncertainties (see Supplementary Figure 22) of the stress model are then estimated from bootstrapping via random replacement of the original unit slip vectors (Supplementary Figure 23), and the angular misfits are shown in Supplementary Figure 24.

A m a n o s s e g m e n t Supplementary Figure 8 .
Quadtrees for down-sampling ALOS-2 LOS displacement maps.Quadtrees applied for down-sampling ALOS-2 LOS displacement maps for the ascending track (a) and the descending track (b), respectively.The longitude and latitude have been shifted with respect to the reference point (35°287 E, 36°N).

Supplementary Figure 9 . 3 -
D view of fault geometry and uniform slip from different perspectives.3-D view of fault geometry and uniform slip on each segment, estimated by the medians of their posterior samples from the Bayesian inversion.N denotes the north direction.
fault geometry, slip distribution, and aftershocks from different perspectives.3-D view of the fault geometry color-coded with slip and the distribution of aftershocks (black dots, Ding et al., 2023).All segments are labeled.The black arrow indicates the north direction.Red stars represent hypocenters determined by the U.S. Geological Survey (USGS), GEOFOrschungsNetz (GEOFON), the Incorporated Research Institutions for Seismology (IRIS) and the Turkey Disaster and Emergency Management Authority (AFAD).

Supplementary Figure 11 .Supplementary Figure 12 .Supplementary Figure 13 .
The fit to the data for the Bayesian inversion of the fault geometry.(a) The horizontal GNSS observations (blue) versus the corresponding model predictions (red).The ALOS-2 lineof-sight observations (b and e), the corresponding model predictions (c and f), and the unmodeled residuals (d and g) from the ascending and descending tracks are color-coded in the same scale.Fault segments are marked with black solid lines.AZ, azimuthal direction.LOS, line of sight.Positive pixel values indicate motion towards the satellite.Data fit against rupture speeds of the M7.8 event.Variance reduction (%) to high-rate GNSS waveforms, strong motion and GNSS offsets as a function of the maximum allowed rupture speed for the M7.8 event.Data fit to the waveform data of the M7.8 event.GNSS (a) and strong motion (d) waveform observations (black) and fits (red) for the M7.8 event.For each plot, the numbers on the right show the peak amplitude values of the waveform.

Table 3
Information of Sentinel-1 radar and Sentinel-2 optical images The 1D velocity model used in this study

Table 5
Mean fault strike from the slip model and the SHmax