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

Domain: earth_science
Published: 2023-09-09
Editorial outcome: major_revision

Abstract
Two major earthquakes (M 7.8 and M 7.7) ruptured left-lateral strike-slip faults of the East Anatolian Fault Zone (EAFZ) on February 6, 2023, causing >59,000 fatalities and ~$119B in damage in southeastern Türkiye and northwestern Syria. Here we derived kinematic rupture models for the two events by inverting extensive seismic and geodetic observations using complex 5-6 segment fault models constrained by satellite observations and relocated aftershocks. The larger event nucleated on a splay fault, and then propagated bilaterally ~350 km along the main EAFZ strand. The rupture speed varied from 2.5-4.5 km/s, and peak slip was ~8.1 m. 9-h later, the second event ruptured ~160 km along the curved northern EAFZ strand, with early bilateral supershear rupture velocity (>4 km/s) followed by a slower rupture speed (~3 km/s). Coulomb Failure stress increase imparted by the first event indicates plausible triggering of the doublet aftershock, along with loading of neighboring faults.

Introduction
The crust of Türkiye is fragmented by escape tectonics, with the Anatolian microplate displacing westward as the Arabian and African plates move northward toward the Eurasian plate (Fig. 1a). This produces active continental faulting along the right lateral strike-slip North Anatolian Fault Zone (NAFZ) and the obliquely intersecting left-lateral strike-slip East Anatolian Fault Zone (EAFZ) 1–9.

The EAFZ is a suite of primarily strike-slip faults formed by the transpressional collision between the Anatolian microplate and the Arabian Plate. The EAFZ bifurcates into a northern strand and a main strand near Celikhan (Fig. 1b). The main strand extends ~700 km with a strike averaging N60°E from the northeastern Karlıova triple junction to near the southwestern gulf of İskenderun 9–13 (Fig. 1a). Some tectonic interpretations identify a candidate Maras triple junction located near Türkoglu, joining the African, Anatolian and Arabian plates 14 with the EAFZ extending relatively straight along the Karatas-Osmaniye Fault Zone, while others extend the EAFZ southwestward along the Amanos fault with strike N35°E 13 to a candidate Amik triple junction 15 near the northern end of the north-south striking Dead Sea Fault (DSF) (Fig. 1b). The DSF bounds the African and Arabian Plates, extending through Syria, Lebanon, Israel, and Jordan. Geodetic and geological studies indicate that the main strand of the EAFZ is divided into several distinct geometric segments by conjugate fractures, parallel faults, pull-apart basins, bends, and stepovers that may govern the size and occurrence of large earthquakes 15–17. The northern strand of the EAFZ involves the Sürgü, Cardak, Savrun-Toprakkale, and Yumurtalik-Düzici-Iskenderun fault segments with varying strikes from E-W to N30°E 15.

The EAFZ was less seismically active than the NAFZ during the twentieth century 18,19. However, inter-seismic geodetic coupling analysis clearly indicated that strain accumulation along the EAFZ 19–21 was capable of producing significant earthquakes, like those that struck the fault in 1513, 1795, and other events between 1822 and 1905 6,19 (Fig. 1a). Estimates of the slip rate along the main EAFZ segment range from ~10 mm/y in the northeast to ~4 mm/y in the southwest where the fault system connects to the DSF 7,22–24. On January 24, 2020, an M 6.7 rupture struck the Doğanyol-Sivrice region of the central EAFZ main strand, northeast of the bifurcation; this was the largest event on the fault in the last 50 years prior to 2023 25–31. The cumulative slip on the main strand of the EAFZ is modest, with estimates of 22–27 km 15,32.

On February 6, 2023, an M 7.8 rupture (denoted the Pazarcık earthquake) initiated on a short, previously unmapped splay fault extending southward from the main of the EAFZ (Fig. 1b), with hypocentral parameters reported by the USGS National Earthquake Information Center (USGS-NEIC) being (37.226°N, 37.014°E, 10 km deep, at 01:17:34.332 UTC). The USGS-NEIC W-phase moment tensor had 81% double couple solution with a best double couple with near-vertical left lateral strike-slip with strike 228°, dip 89°, rake −1° with a seismic moment of 5.389×1020 N-m (M 7.75), while the Global Centroid Moment Tensor (GCMT) solution has a best double couple with strike 51°, dip 70°, rake −4°, M = 5.8×1020 N-m (M 7.8). At 10:24:49.640 UTC, a second large event (denoted the Ekinözü earthquake) with hypocentral parameters (38.011°N, 37.196°E, 7.4 km deep), struck along the northern strand of the EAFZ with 34% double couple W-phase solution with best double couple strike 277°, dip 78°, rake 4° and seismic moment of 2.637×1020 N-m (M 7.55). The GCMT solution for the Ekinözü earthquake has a best double couple with strike 264°, dip 46°, rake −9°, M 0 = 4.53×1020 N-m (M 7.7). This pair of major earthquakes, designated as a doublet because of their similar size (M 7.8 and M 7.7) and close space-time proximity, produced devastating ground motions across southeastern Türkiye and northwestern Syria, responsible for >59,000 fatalities and ~$119B in damage. The events ruptured complex fault networks, involving multiple fault segments resolved by satellite images 33–39. The ground motions for both events were extensively recorded by regional strong-motion accelerometers, GNSS stations, and global broadband seismic stations (Fig. 2, Supplementary Figs. 1, 2), and the recorded signals are herein inverted for kinematic rupture models of the two major events to shed light on the ground motion generation that resulted in regional catastrophe.

Results
Near-fault coseismic displacements from the strong-motion data
The availability of dense near-fault strong-motion observations presents an excellent opportunity to study the detailed rupture processes of the 2023 Türkiye earthquake doublet. Near-field static deformation provides robust constraints on the slip distributions due to its insensitivity to the rupture process and precise Earth structure. To better constrain the coseismic slip distributions of the doublet events, we developed a new baseline correction method to determine static displacements from the near-fault strong-motion data, enhancing a prior approach 40 (see Methods). The robustness and effectiveness of this new method are demonstrated by comparing estimates from strong-motion data with static coseismic GNSS displacements for several historical great earthquakes (see Code Availability). Ultimately, we successfully derived coseismic displacements at 19 near-fault strong-motion stations for the M 7.8 event and four stations for the M 7.7 event (Supplementary Table 1), thus compensating for a paucity of near-fault GNSS-based coseismic displacements (Fig. 2 and Supplementary Fig. 2). The recovered coseismic displacements show prominent sinistral strike-slip characteristics, with the largest horizontal and vertical permanent values for the M 7.8 event being 2.8 m at station 4614 and 0.6 m at station 4615, respectively (Fig. 2). These derived displacements provide valuable constraints on the coseismic slip distributions. The strong-motion derived coseismic displacements are generally consistent with horizontal displacements derived from pixel-tracking offsets of Sentinel-1 satellite radar images 37. However, there are some differences at near-fault stations, as illustrated in Supplementary Fig. 3. Given the inherent uncertainties associated with both approaches, such as the resolution of the pixels and the orientation error of the strong motion stations, these uncertainties inevitably contribute to differences in both magnitude and direction of the derived horizontal displacements.

Kinematic slip model of the M 7.8 Pazarcık earthquake
The February 6, 2023 earthquake doublet in southeastern Türkiye involves the most complex rupture evolution recorded in Türkiye throughout the last century, with backward branching 41. The geometrical complexity of the fault ruptures was well captured by post-earthquake satellite data 33, resolving uncertainty in which fault segments ruptured and constraining the absolute location of the significant faulting with a precision of less than 1 km. Using the satellite data together with the relocated aftershock distribution 42 (Supplementary Fig. 4), we constructed a six-segment fault model (aF1–aF6) for the M 7.8 event, with the model parameters listed in Supplementary Table 2. Based on this constrained geometry, using a well-established nonlinear finite fault inversion method (see Methods), we determined a robustly constrained space-time slip model of the M 7.8 earthquake by joint inversion of seismological and geodetic measurements, including strong-motion data, teleseismic waveforms, static GNSS, high-rate GNSS, and the coseismic displacements derived from strong-motion observations (see data processing in the Methods section).

Our finite fault model of the M 7.8 event indicates that the rupture began on a small fault extending southwestward from the main branch of the EAFZ, then spread onto the main branch, dipping steeply toward the northwest, with slip extending over 160 km to the northeast and terminating southwest of the rupture region of the 2020 Doğanyol-Sivrice M 6.9 event 29 (Fig. 3). Simultaneously, rupture propagated about 180 km toward the southwestern end of the main southeast-dipping EAFZ strand, manifesting a strong bilateral rupture process. The slip distribution exhibits significant spatial heterogeneity, with the largest slip concentrated on aF1 and aF2, presenting a complementary distribution with aftershocks (Fig. 3a). The slip on the other fault segments is relatively low, and these have dense aftershock distributions. Some available finite-fault slip models 36–38,42 show relatively smooth slip variations across much of their fault models. Despite the differences, all models are characterized by a peak slip near the epicenter while showing a minor slip along the northeastern fault segment. The maximum slip is about 11 m, located on bF1, and the total rupture duration is ~65 s. The computed seismic moment M 0 = 5.0×1020 N·m (M 7.7), which is comparable with the GCMT solution (4.97×1020 N·m). The early (<8 s) bilateral rupture propagation had a high rupture velocity of ~4 km/s; subsequent rupture was slightly faster towards the northeast (3.0 km/s, 88% of the local shear wave velocity) than to the southwest (2.7 km/s).

Kinematic slip model of the M 7.7 Ekinözü earthquake
Following the same procedure as for the M 7.8 earthquake, we constructed a 5-fault (bF1-bF5) segment model (Supplementary Table 2) to investigate the coseismic slip model of the M 7.7 event on the northern strand again by joint inversion of strong-motion data, teleseismic waveforms, static GNSS, high-rate GNSS and coseismic displacements derived from strong-motion observations (see Data processing in the Methods section). The preferred model shows that the rupture of the M 7.7 event propagated bilaterally along the east-west strike direction and is primarily characterized by strike-slip offsets with significant shallow motion (Fig. 5a), and the rupture in the multiple southwest segments was more complicated than the northeastern rupture, which propagated parallel to the main strand rather than continuing on a fault extending to the main strand. The slip distribution has substantial spatial heterogeneity, with the largest slip concentrated on bF1 and bF2, presenting a complementary distribution with aftershocks (Fig. 5a). The slip on the other fault segments is relatively low, and these have dense aftershock distributions. Some available finite-fault slip models 36–38,42 show relatively smooth slip variations across much of their fault models. Despite the differences, all models are characterized by a peak slip near the epicenter while showing a minor slip along the northeastern fault segment. The maximum slip is about 11 m, located on bF1, and the total rupture duration is ~65 s. The computed seismic moment M 0 = 5.0×1020 N·m (M 7.7), which is comparable with the GCMT solution (4.97×1020 N·m). The early (<8 s) bilateral rupture propagation had a high rupture velocity of ~4 km/s; subsequent rupture was slightly faster towards the northeast (3.0 km/s, 88% of the local shear wave velocity) than to the southwest (2.7 km/s).

The fits of static displacements and seismic waveforms of all datasets are shown in Fig. 5b and Supplementary Figs. 10, 11, and 12. The 5-fault segment model can satisfactorily explain most observations. However, despite the complexity of the fault model, there is inevitable model oversimplification and neglect of detailed 3D site effects (there is limited available information about the shallow crustal structure near the stations), and the high-frequency content of strong-motion seismic recordings are not all well explained as a result.

Simulated annealing inversions frequently exhibit slight dependence on the chosen random seeds, mainly when multiple optimal solutions exist within the model space, exhibiting indistinguishable objective function values 44. Moreover, the varying random seeds result in distinct initial fault models and Markov chains. To address this uncertainty and explore its impact, we conducted ten inversions for each event in the earthquake doublet using different random seeds in each case. The tests indicate that large-slip distributions of the ten models for the M 7.8 event exhibit relatively stable behavior, with consistency among the models (Supplementary Fig. 13a). In general, the standard deviation (STD) across most fault segments is negligible, with the exception of segments aF3 and aF6 (Supplementary Fig. 13b). Similarly, the STD for the M 7.7 event is typically small compared with the average slip (Supplementary Fig. 13c), but exceptions are found in the western bF1 and bF2 fault segments (Supplementary Fig. 13d). It is suspected that the higher STD in some parts of the fault model is caused by the absence of corresponding very near-fault observations, suggesting the need for further investigation in these areas.

Coseismic Coulomb stress changes and earthquake-triggering effects
This is a rare strike-slip major earthquake doublet with a separation interval of only 9 h; the first M 7.8 earthquake ruptured the main branch of the EAFZ, and the second M 7.7 earthquake ruptured the northern branch of the EAFZ. To investigate the triggering mechanism of the M 7.7 event, we analyzed the Coulomb stress changes induced by the M 7.8 earthquake (see Methods). Due to the significant variability in the estimated dip angle for the larger event, with faulting geometries dipping to the northwest or to the southeast at angles from 42° to 86° being reported by different seismological institutes (Fig. 6), as well as the sensitivity of the results to the receiver fault parameters, we performed analyses using four different receiver fault models with varying parameters (Fig. 6). This provides a more comprehensive exploration of the triggering process by comparing the loading patterns on the initial geometry and location of the M 7.7 event. The calculated results for 10 km depth indicate that allowing for the uncertainty in the precise geometry and slip distribution of the M 7.8 earthquake, the Coulomb stress at the source of the M 7.7 event increased by ~0.014–0.189 MPa for four different receiver target fault geometries, in all cases exceeding the minimal earthquake triggering threshold of ~0.01 MPa 45 (Fig. 6) for favorably oriented static stress change. This is compatible with direct, albeit delayed, triggering of the back-branch rupture. It is also important to remember that triggering is complex, and larger dynamic stresses during the passage of elastic waves from the first event did not immediately trigger failure. Accumulation of pre-stress near the failure limit and favorable orientation of the fault relative to the Coulomb stress perturbation is essential for triggering failure, with large doublet events being relatively rare as a result.

To assess the future effect on seismicity in southeastern Türkiye, we calculated the coseismic stress changes resulting from the combination of the 2023 Türkiye earthquake doublet and the 2020 Doğanyol-Sivrice M 6.7 event 29. We targeted three potential regions that Coulomb stress perturbations may impact: the main northeastern strand of the EAFZ (A1), the DSF to the south (A2), and the region of the Anatolian microplate around the EAFZ (A3). The calculated Coulomb failure stress in A1 increased by up to 0.1 MPa for EAFZ receiver geometry (as shown in Fig. 7), which is concerning given that the most recent large events in the area northeast of Lake Hazar occurred in 1874 and 1866 46 (Fig. 1a). In A2, the receiver geometry of the left-lateral strike-slip Dead Sea fault, located just south of the main shock rupture along the Amanos Fault, is calculated to have a loading increase (up to 0.1 MPa), suggesting an advance toward the next rupture. Numerous parallel strike-slip faults exist in region A3 as a result of distributed tectonic activity in the transpressional regime. For a receiver geometry given by the average orientation of these faults, Coulomb stress changes were computed, revealing a positive stress change zone towards the west of the northern strand. Given that significant delays ranging from years to decades between mainshocks and major aftershocks are frequently observed worldwide 45, identifying possible future rupture zones based on the stress perturbations from the recent faulting can assist in directing mitigation efforts toward these regions.

Methods
Dataprocessing
Teleseismic data. We selected 40 P wave and 26 SH wave broadband waveforms for the M 7.8 earthquake and 40 P wave and 36 SH wave for the M 7.7 earthquake from the IRIS data management center based on high signal-to-noise ratio and well-distributed azimuthal coverage at teleseismic (30°–90°) distances (Supplementary Fig. 1). We then removed instrument responses to obtain ground displacements with durations of 100 s for the M 7.8 event and 60 s for the M 7.7 event, in the passband 1 s–300 s. Finally, we precisely aligned all the P and SH wave initial motions manually.

Geodetic observations. We chose the displacement time series at six GNSS stations for the M 7.8 event and five GNSS for the M 7.7 event, respectively, from Türkiye Ulusal Sabit GNSS Ağı (TUSAGA-Aktif) (see Data Availability), which were computed by PRIDE PPP-AR 54. All data were re-sampled at 0.2 s intervals, and a time window of 300 s was used for the joint inversion. The first-motion arrivals of all ground displacement waveforms are hand-picked.

We also selected coseismic displacements at 29 GNSS sites (Fig. 2) for the M 7.8 earthquake and seven GNSS sites (Supplementary Fig. 2) for the M 7.7 earthquake from 5-min sample rate time series derived with rapid orbits by the Jet Propulsion Laboratory (see Data Availability). Due to the relatively low precision of GNSS positions for vertical components, only the near-fault vertical component recorded at station EKZ1 was utilized in the joint inversion for the M 7.7 earthquake.

Strong-motion data. Strong-motion data used in this study were recorded and provided by the Disaster and Emergency Management Presidency of Türkiye (AFAD-TK) (see Data Availability). Usually, strong-motion records include different sources of noise affecting the information to be retrieved. The most well-known problem is caused by shifts in the reference baseline, which prevent accurate ground velocity and displacement recovery through integration. The 2023 Türkiye earthquake doublet was well captured by dense strong-motion stations. To obtain high-quality velocity waveforms and coseismic permanent displacements, here, we present an updated scheme for the bi-linear baseline correction approach of Wang et al. 40. To minimize the uncertainty of this approach, we replace the broken-line correction with a natural curve correction obtained through iterative smoothing of the uncorrected velocity seismogram to better recover the ground velocities, ground displacements, and permanent coseismic offsets. The new scheme consists of the following three steps.

Step 1. Integrate raw accelerograms to the uncorrected velocity seismogram
First, for a raw accelerogram a raw(t) within the time window t 0 to t end with the known P wave arrival t 0 within the time window. In practice, we suggest a pre-seismic window t pre - t 0 between 5 s and 30 s, and a generous signal window t end - t pre.

Second, estimate the pre-seismic baseline offset
Δa pre = 1/(t pre - t 0) ∫ t 0 t pre a raw(τ)dτ, (1)
and remove it from the raw accelerogram to get an accelerogram including only seismically induced baseline errors, a 0(t) = a raw(t) - Δa pre. (2)

Third, to estimate the time when the co-seismic baseline shift is stabilized to a constant post-seismic offset, we introduce function
E(t) = ∫ 0 t pst |a 0(τ)|dτ, (3)
and time t γ satisfying E(t γ) = γE(t pst), and assume that t pst = t γ = 85% can be regarded as the time when the co-seismic baseline shift has been stabilized.

Finally, integrate a 0(t) to velocity seismogram,
v 0(t) = ∫ t 0 t a 0(τ)dτ: (4)

Step 2. Estimate post-seismic baseline shift and the starting velocity correction.
First, calculate the post-seismic linear trend of v 0(t) via least-squares regression
f(t) = v pst + (v end - v pst)/(t end - t pst) * (t - t pst), (5)
where v pst and v end are the start and end value of f(t) at t = t pst and t end, respectively.

Second, define another function
g(t) = { v 0(t), t 0 ≤ t < t pre, w(t), t pre ≤ t < t pst, f(t), t pst ≤ t ≤ t end } (6)
as the starting correction curve, where the function w(t) is a weighted average of v 0(t) and f(t), i.e., the sum of right-tapered v 0(t) and left-tapered f(t).

Step 3. Get final velocity correction via iterative smoothing
First, smooth g(t) iteratively using a small moving window, but fixing g(t pre) = 0 and g(t end) = v end,
g(t) := 1/2Δt * ∫ t-Δt t+Δt g(τ)dτ, t pre < t < t pst, (7)
where Δt is the time sample and a := b means updating a by b. The smoothing process is terminated when g(t) has no extrema after t pre and at maximum only one extremum before t pst. So g(t) becomes a smooth and, in most cases, monotonic curve.

Second, set the final velocity correction curve v err(t) = 0 for t 0 ≤ t ≤ t pre and v err(t) = g(t) for t pst ≤ t ≤ t end. For the remaining co-seismic period t pre < t < t pst, v err(t) need to be estimated specially. In the case that the co-seismic and post-seismic baseline shift have the same sign, i.e., v pst * (v end - v pst) ≥ 0, and the former is smaller than the latter, i.e., |v pst| < |v end - v pst|, we shift t pre rightward to max[t pre, (t pre + 2 * t pst - t z)/3], where t z is the time of zero-crossing of the post-seismic trend f(t).

Third, construct a monotonic v err(t) in the co-seismic period t pre < t < t pst, which best fits v 0(t) in this period in the least-squares sense.

Finally, make the baseline correction on the velocity seismogram
v(t) = v 0(t) - v err(t) (8)
and integrate it into the corrected displacement seismogram
u(t) = ∫ 0 t v(τ)dτ: (9)

Using the correction procedure outlined above, we successfully corrected strong-motion waveforms at 52 stations for the M 7.8 event and 26 stations for the M 7.7 event, and obtained stable coseismic displacements at 21 near-fault stations for the M 7.8 event and four near-fault stations for the M 7.7 event, respectively (Supplementary Table 1). Before the joint inversion, all regional seismic waveforms were filtered with a bandpass filter of 0.02–0.5 Hz and sampled at 0.2 s intervals. The first-motion arrivals of all ground velocity waveforms are hand-picked, and a time window of 300 s was used for the joint inversion. Verification of the ground displacement estimation from the strong-motion recordings using the foregoing procedure is provided by applications to large datasets of strong-motion and GNSS observations for the 2014 M 8.2 Pisagua earthquake, 2014 M 7.6 Iquique earthquake, and 2011 M 9.0 Tohoku-Oki earthquake (see Code Availability). Favorable recovery of both horizontal and vertical displacements is achieved.

Finite-fault inversion
A combined analysis of seismic and geodetic data is very effective in understanding the rupture process of large earthquakes. So, we utilized both data types to invert the rupture process of the M 7.8 and M 7.7 events using a multi-segment fault model with geometries determined by surface rupture 33 and relocated aftershocks 42 (Supplementary Fig. 4). A nonlinear finite fault inversion method is employed 55,56, which can simultaneously invert geodetic and seismic observations in the wavelet domain. The sum of L1 and L2 norms of the seismograms in different wavelets quantifies the misfit between the recorded and synthetic waveforms. Sum-squared residuals have been adopted as the evaluation criteria to measure the difference between observed and synthetic static displacements. All inversions commence with a randomly generated initial model with a total moment equal to the GCMT solution. The weight assigned to the static error is set to be equal to the waveform error, but for the statics, the weight on the coseismic displacements derived from the strong-motion data is taken as half of GNSS statics accounting for the inherent uncertainties associated with baseline correction. All inversion parameters of this earthquake doublet are presented in Supplementary Tables 2, 3. All Green’s functions for both statics and waveforms are computed using a regional 1D velocity model 57.

Coulomb failure stress
The Coulomb failure stress (ΔCFS) change can be defined as 58: ΔCFS = Δτ + μ, Δσ N, where Δτ and Δσ N are changes in the shear stress and normal stress on a receiver fault, respectively, caused by the earthquake. In this study, the friction coefficient (μ,) was set to 0.4 as a common choice. The values of Δτ and Δσ N are defined with respect to the slip and normal directions of the receiver fault, respectively. Hence, a positive value of ΔCFS indicates that the earthquake-induced stress changes push the receiver fault closer to rupture, while a negative value of ΔCFS suggests the opposite.

Using the code PSGRN/PSCMP 59, we calculated the coseismic Coulomb stress change at the location of the M 7.7 event caused by the M 7.8 earthquake, as well as the evolution of ΔCFS on the surrounding main faults caused by the combined stress contributions from the 2023 Türkiye earthquake doublet and the 2020 Doğanyol-Sivrice M 6.9 event 29.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper presents detailed slip models for the destructive earthquake doublet that struck\nsouthwestern Turkiye in early 2023. A few models have been published recently, but we are very\nmuch in the early stages of assessing the source of these events. A key feature of the models in\nthis paper is that they have used an extensive set of near-field data from GPS displacements and\nstrong motion records. These data are critical for localizing the slip in the model in space and time.\nOverall, the models look good to me. The figures are mostly good, but the arrows in the panels\nlike Figure 5a are too small to see. Overall, I think this is a valuable contribution that will advance\nour knowledge of the earthquake source for these events, and I think publication after minor\nrevisions is reasonable.\nI marked an annotated manuscript. There are some awkward phrasings noted, and some English\ncorrections suggested. Awkward phrasing markings should be interpreted to mean that the\nauthor’s intended meaning could not be uniquely interpreted from their words, so they need to\nrephrase the text to be more clear. However, overall the writing and organization are clear and the\npaper is easily understood.\nLines 46-51. This discussion is not so clear, and the potential triple junction points are not clearly\nindicated on the figures. Perhaps show the alternative geometries with different colors or line\nstyles in Figures 1 or 2?\nLines 354-357. The use of “postpone” and “prepone” is technically correct, but “prepone” is a\nrarely used word. It might be simpler to use “advance” and “delay”. Shifting t_pre to an earlier\ntime should not cause any issues, but since t_pre is supposed to be the P-wave arrival time,\nshifting t_pre to a later time means that the P waves would be included in the pre-earthquake\nbaseline window, which seems problematic. Please clarify. Also, there is a supplemental figure that\nshows this method applied to the Tohoku earthquake, which does not appear to be referenced in\nthe text. The method should end with a comparison of the estimated static displacements from\nstrong motion records to observed GPS displacements, which I think they need to use the Tohoku\nearthquake for. So say that in the main text, and refer to the figure.\nLines 390-398. Please add the EOS article reference provided for the UNR products, in addition to\nthe specific URLs here.", "reviewer_2": "+++++++++++++++++++++++\nREVIEW FOR 427543\nComplex multi-fault rupture and triggering 1 during the February 6, 2023, earthquake doublet in\nsoutheastern Türkiye\nBy Chengli Liu and co-authors\n++++++++++++++++++++\nEarthquake source inversions are key to understand the physics of the earthquake rupture\nprocess, which in turn helps to unravel the underlying causes (tectonics; acting stresses; fault-\nzone complexities) of this particular event, but also allow to make inferences on what to possibly\nexpect in future quakes. The earthquake sequence in Turkey of February 06, 2023, was a\nparticular violent one in terms of shaking & damage, but also in terms of magnitudes of the two\nevents (M 7.8 and M 7.6 only 9 hrs later), and hence studying these earthquakes in detail is of\ngreat importance.\nIn this study, Liu et al use a combined data set of geodetic and seismic data for finite-fault\ninversions of this earthquake doublet. They also calculate the Coulomb Failure Stress (CFS)\nresulting from these ruptures (assuming alternative fault geometries) to understand (i) how the\nfirst earthquake may have facilitated the triggering of the 2nd event, and (ii) how the stresses\nhave changed in the surroundings of these ruptures to potentially bring nearby faults closer to\nfailure.\nThe paper is well written, the data analysis/interpretation solid, and the inverted finite-fault\nrupture model are in agreement with recent findings in several other studies. In this context, a few\nreferences are missing (Barbot et al, 2023, in Seismic; Mai et al, 2023, in The Seismic Record;\nGoldberg et al, 2023, in The Seismic Record; Petersen et al, 2023, The Seismic Record) which the\nauthors should look at and compare their results with.\nI have moderate-to-major comments, mostly editorial and for clarification, which I list below (sort\nof sequential).\nThe finite-fault inversion approach is not well documented/explained, uncertainties are not stated,\nthe weighting of the different data sets is not explained, the fault discretization etc. This makes it\nimpossible to asses the robustness of the inferred models.\nAbstract\n+ L 23 — perhaps better “>59,000 fatalities;\n+ the abstract should tell which data are used and what has been done, before giving key results.\nThis should be stated at the current end of Line 24, before rupture nucleation and propagation is\npresented. In some sense, Line 31-34 need to be moved forward, and then need to be expanded a\nbit\nIntroduction\n+ L 82: See above comment to Line 23\n+ L 92: earth —> Earth\n+ L 101: the authors here refer to the 2nd event to be of Mw 7.7. However, this has not been\nmeasured/inferred yet. In Line 80, they mention the USGS NEIC W-phase estimate of Mww 7.55,\nso this value of Mw 7.7 comes at a surprise and is confusing. Many other studies report M 7.6 for\nthe 2nd event. Please clarify.\nKinematic Slip Models\n+ L 131: “patchy heterogeneity” — somewhat awkward expression … perhaps “spatial\nheterogeneity”?\n+ L 133: unclear why variations in ’minor normal and thrust faulting’ would be related to ‘lateral\ninhomogeneities of crustal rocks’ (assuming you mean the velocity-density structure here).\nPerhaps rephrase …\n+ L 154: remove ‘oversimplified’\nTriggering and Seismic Hazard Evaluation\n+ This section header is misleading, as no seismic hazard calculation/evaluation is done. At most,\nthere could be a “clock advance/delay”. Please rephrase\n+ L 190: add some references to the statement of ‘significant variability in dip angles’,\nand\nperhaps quantify this variability (5 degrees; 25 degrees?)\n+ Line 195-201: Here it should be stated at which depth the CFS-values are calculated / reported.\nIs a stress increase of ~0.014 “significant”? And why is “~0.01 MPa” the minimal triggering\nthreshold? This should be explained , and backup with references\n+ Line 202: seismic “risk” is not assessed in this study. Not even seismic hazard. See comment\nabove.\n+ Line 212: please quantify “incremental loading increase”\nDiscussion\n+ Line 253: I find the argument on ‘significant non-uniformity of stress accumulation’ somewhat\nconfusing, in particular when it is then linked to “high stress buildup” two lines below. Non-\nuniformity due to small-scale fault segmentation or fault roughness? Any evidence of this due to,\nsay, variations in moment-tensor solutions in the local background seismicity there?\n+ Lines 260-280: this is interesting argument, which could be further quantified by computing the\nseismic-radiation efficiency. It is actually quite sobering and counterintuitive if these very\ndamaging earthquakes had low radiation efficiency …\n+ Line 285 - 298: This seems rather speculative, in my opinion, suggesting that super-shear\nrupture may even be common (or the norm?) for large earthquakes on the NAF and EAF. I don’t\nthink the case of prevalence of supershear event can be made, or even should be made, based on\nonly four documented observations.\nFinite Fault Inversion\n+ there is no mention of how the different datasets are weighted;\n+ there is no mention of uncertainties at all;\n+ there is no mention of how misfits are computed (variance reduction?)\n==> this section needs to be expanded. The readers need information on robustness and\nuncertainties of the solutions\nCoulomb Stress\n+ Line 418: “we calculated the triggering” … this is actually not correct. Only the static stress\nchanges were computed. It is then inferred/assumed that somewhat higher stress near the\nhypocenter may have initiated / facilitated the 2nd event. But there is dynamic rupture simulation\ncarried out that actually shows that …\nReferences\n+ several recent publications that also conducted finite-fault inversions are not listed, but should\nbe added (see above)\nFigures\n+ Figure 3 and 4, panels: from which fault-normal distance were aftershocks projected onto the\nfault plane? Up to 5 km on either side of the fault plane? More or less?\n+ Figures 8 and 9 are confusing: they show time-slices (snapshots) of the rupture process, but in\neach time slices, the rupture-time contours of the entire rupture is shown. Shouldn’t the time-slice\n0-10 sec only contain contours for the first 10 second?\nAlso, I find the color-range choice suboptimal, as tiny variations in shades of blue cannot be\ndistinguished. Regions of zero slip should be shown in white, or light gray, and thus be clearly\ndistinct from regions with 0.5-1.0 m. This is not case here.\nSupplementary Material\nFigure 3 is confusing here and should be removed. Maps related to the 2011 Tohoku earthquake\nare not needed here", "reviewer_3": "Authors presented a thorough analysis on the rupture processes and possible impacts of the\ndevastating Turkey earthquake doublet. The results suggest that both events ruptured multiple\nfault segments, featured supershear rupture episodes, and the first earthquake possibly aided the\ntriggering of the second event. I like the idea of enriching the static displacement measurements\nusing the near-field strong motion accelerograms. However, I have some concerns with the\nvalidity of the data product and results, as well as some unclear delivery of the method approach.\nAcknowledging that the implication of this work is important, here I suggest a major revision, see\nthe following points for details.\nMajor comments:\n1. The fault geometry is one major assumption and is critical to the finite fault inversions\nconducted in this study. How the fault geometries, particularly the dip angles, are constructed, is\nunclear, which raises concerns about the validity of the findings. Additionally, the absence of\nsatellite data and sources further compounds the issue.\n2. One of the major claims is the supershear stages found during both the northeast and\nsouthwest rupture episodes of the first event. However, how to validate (quantify the uncertainty)\nof such estimates? To date, most finite-fault inversions and back projection-based studies suggest\nsubshear ruptures along the EAF during the first event. These includes:\nMeng et al. (https://doi.org/10.21203/rs.3.rs-2747911/v1)\nMai et al. (https://doi.org/10.1785/0320230007)\nGoldberg et al. (https://doi.org/10.1785/0320230009)\nMelgar et al. (https://doi.org/10.26443/seismica.v2i3.387)\nAmong them, Meng et al. conducted a Mach Cone analysis arguing that the rupture is\npredominantly subshear.\n3. Authors’ slip models have significant slip gaps along the northeast segment of the East\nAnatolian Fault Zone, and the west of the Sürgü Fault. However, this appears to contradict the\nsatellite surface displacement data which are continuous for both events (see Mai et al. Fig. 2;\nGoldberge et al. Fig. 4). How certain are these slip gaps? The quickest way is to do a forward\nprediction of the satellite displacement field and compare with the data.\n4. The conversion from strong motion accelerograms to static displacements impressively enriched\nthe near-fault static observations. However, I am concerned with the robustness. For example, the\nstatic slip of station 4615 during the first earthquake has a 1.5m amplitude and is perpendicular to\nthe dominant strike slip direction of the East Anatolian Fault Zone, whereas the satellite radar\nimages do not show any significant southeast motions (Mat et al, Goldberg et al.). Authors need to\nmake sure that such conversions are valid.\n5. 0.1 MPa is still small for static triggering of the Mw 7.7 second earthquake. Is there a source of\nminimal earthquake triggering threshold of ~0.01 MPa for M7 earthquakes?\nMinor comments:\n1. The colorbars for slip amount are severely saturated in Figures 3, 5, 8, 9. I can’t distinguish\nanything between 6 and 12 m.\n2. Still for these figures (especially Fig. 3 and 5), the slip models (subpanel a) are inversely\noriented compared with the map (subpanel b). I had a hard time to realize that the left-hand side\nindicates the faults to the east."}, "author_response": "Response to reviews of \"Complex multi-fault rupture and triggering during the\nFebruary 6, 2023, earthquake doublet in southeastern Türkiye\" by Liu et al.,\nsubmitted to Nature Communications. The review comments are reproduced below\nin black type, with our responses and indications of how we have revised the\nmanuscript to address the reviews indicated in blue.\n\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nThis paper presents detailed slip models for the destructive earthquake doublet that\nstruck southwestern Turkiye in early 2023. A few models have been published\nrecently, but we are very much in the early stages of assessing the source of these\nevents. A key feature of the models in this paper is that they have used an extensive\nset of near-field data from GPS displacements and strong motion records. These\ndata are critical for localizing the slip in the model in space and time.\nOverall, the models look good to me. The figures are mostly good, but the arrows in\nthe panels like Figure 5a are too small to see. Overall, I think this is a valuable\ncontribution that will advance our knowledge of the earthquake source for these\nevents, and I think publication after minor revisions is reasonable.\nWe modified Figure 5a to make the arrows clearer.\nI marked an annotated manuscript. There are some awkward phrasings noted, and\nsome English corrections suggested. Awkward phrasing markings should be\ninterpreted to mean that the author's intended meaning could not be uniquely\ninterpreted from their words, so they need to rephrase the text to be more clear.\nHowever, overall the writing and organization are clear and the paper is easily\nunderstood.\nWe appreciate the reviewer's comments and suggestions, and we reworded the\nawkward phrasings in the revised manuscript. All comments in the annotated\nmanuscript were addressed.\nLines 46-51. This discussion is not so clear, and the potential triple junction points\nare not clearly indicated on the figures. Perhaps show the alternative geometries\nwith different colors or line styles in Figures 1 or 2?\nWe clarified the discussion and revised Figure 1b by adding two diamond symbols\nto locate the potential triple junctions.\nLines 354-357. The use of \"postpone\" and \"prepone\" is technically correct, but\n\"prepone\" is a rarely used word. It might be simpler to use \"advance\" and \"delay\".\nShifting t_pre to an earlier time should not cause any issues, but since t_pre is\nsupposed to be the P-wave arrival time, shifting t_pre to a later time means that the\nP waves would be included in the pre-earthquake baseline window, which seems\nproblematic. Please clarify. Also, there is a supplemental figure that shows this\nmethod applied to the Tohoku earthquake, which does not appear to be referenced\nin the text. The method should end with a comparison of the estimated static\ndisplacements from strong motion records to observed GPS displacements, which I\nthink they need to use the Tohoku earthquake for. So say that in the main text, and\nrefer to the Figure.\nThanks for the suggestion. A significant baseline shift does not necessarily start with\nthe P-wave arrival. To possibly avoid over-corrections within the coseismic period,\nit is sometimes necessary to delay t_pre, particularly when the average of transient\nbaseline shift appears with the same sign as but smaller in magnitude than the\npermanent one. Of course, the way to delay t_pre that we suggested is empirically-\nbased. We revised the descriptions of the method for strong-motion baseline\ncorrection to make it more straightforward. We removed the supplemental Fig. 3\nshowing an application to Tohoku event data to avoid confusion, but added a new\ndatabase in the section of Code Availability (https://zenodo.org/record/8058010)\nto demonstrate the effectiveness of the new correction method.\nLines 390-398. Please add the EOS article reference provided for the UNR products,\nin addition to the specific URLs here.\nWe added it as suggested.\nWe thank the reviewer for their comments, and the manuscript is improved by our\nrevisions to address them. This is noted in the revised acknowledgments.\n\nReviewer #2 (Remarks to the Author):\n+++++++++++++++++++++++\nREVIEW FOR 427543\nComplex multi-fault rupture and triggering 1 during the February 6, 2023,\nearthquake doublet in southeastern Türkiye\nBy Chengli Liu and co-authors\n++++++++++++++++++++\nEarthquake source inversions are key to understand the physics of the earthquake\nrupture process, which in turn helps to unravel the underlying causes (tectonics;\nacting stresses; fault-zone complexities) of this particular event, but also allow to\nmake inferences on what to possibly expect in future quakes. The earthquake\nsequence in Turkey of February 06, 2023, was a particular violent one in terms of\nshaking & damage, but also in terms of magnitudes of the two events (M 7.8 and M\n7.6 only 9 hrs later), and hence studying these earthquakes in detail is of great\nimportance.\nIn this study, Liu et al use a combined data set of geodetic and seismic data for finite-\nfault inversions of this earthquake doublet. They also calculate the Coulomb Failure\nStress (CFS) resulting from these ruptures (assuming alternative fault geometries)\nto understand (i) how the first earthquake may have facilitated the triggering of the\n2nd event, and (ii) how the stresses have changed in the surroundings of these\nruptures to potentially bring nearby faults closer to failure.\nThe paper is well written, the data analysis/interpretation solid, and the inverted\nfinite-fault rupture model are in agreement with recent findings in several other\nstudies. In this context, a few references are missing (Barbot et al, 2023, in Seismic;\nMai et al, 2023, in The Seismic Record; Goldberg et al, 2023, in The Seismic Record;\nPetersen et al, 2023, The Seismic Record) which the authors should look at and\ncompare their results with.\nWe added the recently published references and made the comparisons (see below)\nas suggested in the revised manuscript.\n“The slip distribution exhibits significant spatial heterogeneity, characterized by\npredominant strike-slip motion with minor occurrences of normal or thrust faulting\n(Fig. 3a). This pattern aligns closely with published models36-38, highlighting the\npresence of lateral variations in tectonic stress, frictional properties within the crust,\nand intricate fault zone structures along the rupture.”\n“. Some available finite-fault slip models36-38,42 show relatively smooth slip variations\nacross much of their fault models. Despite the differences, all models are characterized\nby a peak slip near the epicenter while showing a minor slip along the northeastern\nfault segment.”\nI have moderate-to-major comments, mostly editorial and for clarification, which I\nlist below (sort of sequential).\nThe finite-fault inversion approach is not well documented/explained, uncertainties\nare not stated, the weighting of the different data sets is not explained, the fault\ndiscretization etc. This makes it impossible to asses the robustness of the inferred\nmodels.\nWe added more text (see below) to clarify the finite-fault inversion method,\nuncertainties, the weighting of the different data sets, and fault discretization in the\nrevision.\nResults:\n“Simulated annealing inversions frequently exhibit slight dependence on the chosen\nrandom seeds, mainly when multiple optimal solutions exist within the model space,\nexhibiting indistinguishable objective function values44. Moreover, the varying random\nseeds result in distinct initial fault models and Markov chains. To address this\nuncertainty and explore its impact, we conducted ten inversions for each event in the\nearthquake doublet using different random seeds in each case. The tests indicate that\nlarge-slip distributions of the ten models for the M 7.8 event exhibit relatively stable\nW\nbehavior, with consistency among the models (Supplementary Fig. 13a). In general, the\nstandard deviation (STD) across most fault segments is negligible, with the exception\nof segments aF3 and aF6 (Supplementary Fig. 13b). Similarly, the STD for the M 7.7\nW\nevent is typically small compared with the average slip (Supplementary Fig. 13c), but\nexceptions are found in the western bF1 and bF2 fault segments (Supplementary Fig.\n13d). It is suspected that the higher STD in some parts of the fault model is caused by\nthe absence of corresponding very near-fault observations, suggesting the need for\nfurther investigation in these areas.”\nMethod:\n“The sum of L1 and L2 norms of the seismograms in different wavelets quantifies the\nmisfit between the recorded and synthetic waveforms. Sum-squared residuals have\nbeen adopted as the evaluation criteria to measure the difference between observed\nand synthetic static displacements. All inversions commence with a randomly\ngenerated initial model with a total moment equal to the GCMT solution. The weight\nassigned to the static error is set to be equal to the waveform error, but for the statics,\nthe weight on the coseismic displacements derived from the strong-motion data is\ntaken as half of GNSS statics accounting for the inherent uncertainties associated with\nbaseline correction.”\nAbstract\n+ L 23 — perhaps better \">59,000 fatalities;\nWe changed the text as suggested.\n+ the abstract should tell which data are used and what has been done, before giving\nkey results. This should be stated at the current end of Line 24, before rupture\nnucleation and propagation is presented. In some sense, Line 31-34 need to be\nmoved forward, and then need to be expanded a bit\nWe revised the abstract as suggested.\nIntroduction\n+ L 82: See above comment to Line 23\n+ L 92: earth —> Earth\nWe changed the text, as suggested.\n+ L 101: the authors here refer to the 2nd event to be of Mw 7.7. However, this has\nnot been measured/inferred yet. In Line 80, they mention the USGS NEIC W-phase\nestimate of Mww 7.55, so this value of Mw 7.7 comes at a surprise and is confusing.\nMany other studies report M 7.6 for the 2nd event. Please clarify.\nThe Mw 7.7 for the second event is based on the revised GCMT solution. We added a\ndescription in the introduction.\nKinematic Slip Models\n+ L 131: \"patchy heterogeneity\" — somewhat awkward expression … perhaps\n\"spatial heterogeneity\"?\nWe corrected it, as suggested.\n+ L 133: unclear why variations in' minor normal and thrust faulting' would be\nrelated to 'lateral inhomogeneities of crustal rocks' (assuming you mean the\nvelocity-density structure here). Perhaps rephrase …\nWe rephrased the statement (see below) that predominantly strike-slip motion\naccompanied by minor normal or thrust faulting indicates lateral inhomogeneities\nof tectonic stress and friction in the crust and fine fault zone structure along the\nrupture.\n“The slip distribution exhibits significant spatial heterogeneity, characterized by\npredominant strike-slip motion with minor occurrences of normal or thrust faulting\n(Fig. 3a). This pattern aligns closely with published models36-38, highlighting the\npresence of lateral variations in tectonic stress, frictional properties within the crust,\nand intricate fault zone structures along the rupture.”\n+ L 154: remove 'oversimplified'\nWe removed the word, as suggested.\nTriggering and Seismic Hazard Evaluation\n+ This section header is misleading, as no seismic hazard calculation/evaluation is\ndone. At most, there could be a \"clock advance/delay\". Please rephrase\nWe changed the title of this section to “Coseismic Coulomb stress changes and\nearthquake triggering effects”\n+ L 190: add some references to the statement of 'significant variability in dip\nangles', and perhaps quantify this variability (5 degrees; 25 degrees?)\nFigure 6 in the main text indicates the dip angle estimates from different\nseismological measures (used as target fault geometries), ranging from 42° to 86°,\nwith either NW or SE plunge. We expand the text to discuss this in the revision.\n“Due to the significant variability in the estimated dip angle for the larger event, with\nfaulting geometries dipping to the northwest or to the southeast at angles from 42° to\n86° being reported by different seismological institutes (Fig. 6),”\n+ Line 195-201: Here it should be stated at which depth the CFS-values are\ncalculated / reported. Is a stress increase of ~0.014 \"significant\"? And why is \"~0.01\nMPa\" the minimal triggering threshold? This should be explained , and backup with\nreferences\nThe 10 km depth for the calculation is added to the text and indicated in Fig. 6.\nPrevious studies have confirmed that a change in Coulomb stress value between\n0.01 and 0.1 MPa can be enough to trigger a subsequent earthquake, with 0.01 MPa\nbeing a threshold value for earthquake-triggering. We add a reference to Ross\nStein’s work in the revision.\n+ Line 202: seismic \"risk\" is not assessed in this study. Not even seismic hazard. See\ncomment above.\nYes, we agree that we do not formally address risk or hazard per se, but we do\nquantify stress perturbations that are likely to influence future activity. This is\nclarified in the revision. It has been demonstrated by many studies that significant\nCoulomb stress changes can promote or inhibit subsequent earthquake occurrence.\nThis is significant for earthquake interactions and seismic hazard evaluations, and\ncan be utilized for assessing aftershock migration in the future. In our study, we\nidentified notable increases in Coulomb stress along the Dead Sea Fault and along\nthe northeastern segment of the EAF, which is important for future seismic\nrisk/hazard assessments to consider.\n+ Line 212: please quantify \"incremental loading increase\"\nWe added a specific value in the revision, as suggested.\n“In A2, the receiver geometry of the left-lateral strike-slip Dead Sea fault, located just\nsouth of the mainshock rupture along the Amanos Fault, is calculated to have a\nloading increase (up to 0.1 MPa), suggesting an advance toward the next rupture.”\nDiscussion\n+ Line 253: I find the argument on 'significant non-uniformity of stress\naccumulation' somewhat confusing, in particular when it is then linked to \"high\nstress buildup\" two lines below. Non-uniformity due to small-scale fault\nsegmentation or fault roughness? Any evidence of this due to, say, variations in\nmoment-tensor solutions in the local background seismicity there?\nNo major earthquakes have occurred on the EAF since the 19th century, and the\nbackground seismicity has produced relatively few recent focal mechanisms directly\nalong the relevant section of the EAF. Historical studies referenced in the\nmanuscript (Duman & Emre, 2013, https://doi.org/10.1144/SP372.14; Güvercin et\nal., 2022, https://doi.org/10.1093/gji/ggac045) indicate that the EAF exhibits\ncomplex behavior, which is likely influenced by various factors, including oblique\nplate motion, variations in seismic coupling and fault maturity, and geometric\ncomplexities. Similar to determinations in other studies, our slip models clearly\ndemonstrate a significantly non-uniform distribution in slip and rupture velocity,\nhighlighting the uneven release of stress along the southwest section of the EAF. As\nwe discuss, the EAF exhibits variations in stress orientations and the strain rate\nfield, which is consistent with the diverse solutions for focal mechanisms. Our slip\nmodels show a non-uniformity of stress release along the fault, so we modify our\nlanguage to \"non-uniformity of stress release\" instead of \"significant non-uniformity\nof stress accumulation\" in the revision.\n+ Lines 260-280: this is interesting argument, which could be further quantified by\ncomputing the seismic-radiation efficiency. It is actually quite sobering and\ncounterintuitive if these very damaging earthquakes had low radiation efficiency\nYes, the modest moment-scaled radiated energy is an important attribute of these\nevents. Quantifying radiation efficiency also requires very reliable estimates of\nstress drop along with apparent stress, and for such a complex rupture, our\nconfidence in stress drop values is low, as there is much data- and model-\ndependence of the estimates. At this time, we prefer not to pursue radiation\nefficiency estimation.\n+ Line 285 - 298: This seems rather speculative, in my opinion, suggesting that\nsuper-shear rupture may even be common (or the norm?) for large earthquakes on\nthe NAF and EAF. I don't think the case of prevalence of supershear event can be\nmade, or even should be made, based on only four documented observations.\nWe modify the wording to make it less speculative, with only about a dozen\nconvincing cases of supershear faulting, having 4 of them in Anatolia is quite a\nconcentration, but we agree it is not a lot of cases.\nFinite Fault Inversion\n+ there is no mention of how the different datasets are weighted;\n+ there is no mention of uncertainties at all;\n+ there is no mention of how misfits are computed (variance reduction?)\n==> this section needs to be expanded. The readers need information on robustness\nand uncertainties of the solutions\nAs suggested, we expanded this section in the revision (see below), to clarify the\nweights, uncertainties (revised Supplementary Fig 12), and misfits calculations in\nour finite fault inversions.\nResults:\n“Simulated annealing inversions frequently exhibit slight dependence on the chosen\nrandom seeds, mainly when multiple optimal solutions exist within the model space,\nexhibiting indistinguishable objective function values44. Moreover, the varying\nrandom seeds result in distinct initial fault models and Markov chains. To address this\nuncertainty and explore its impact, we conducted ten inversions for each event in the\nearthquake doublet using different random seeds in each case. The tests indicate that\nlarge-slip distributions of the ten models for the MW 7.8 event exhibit relatively stable\nbehavior, with consistency among the models (Supplementary Fig. 13a). In general, the\nstandard deviation (STD) across most fault segments is negligible, with the exception\nof segments aF3 and aF6 (Supplementary Fig. 13b). Similarly, the STD for the MW 7.7\nevent is typically small compared with the average slip (Supplementary Fig. 13c), but\nexceptions are found in the western bF1 and bF2 fault segments (Supplementary Fig.\n13d). It is suspected that the higher STD in some parts of the fault model is caused by\nthe absence of corresponding very near-fault observations, suggesting the need for\nfurther investigation in these areas.”\nMethod:\n“The sum of L1 and L2 norms of the seismograms in different wavelets quantifies the\nmisfit between the recorded and synthetic waveforms. Sum-squared residuals have\nbeen adopted as the evaluation criteria to measure the difference between observed\nand synthetic static displacements. All inversions commence with a randomly\ngenerated initial model with a total moment equal to the GCMT solution. The weight\nassigned to the static error is set to be equal to the waveform error, but for the statics,\nthe weight on the coseismic displacements derived from the strong-motion data is\ntaken as half of GNSS statics accounting for the inherent uncertainties associated with\nbaseline correction.”\nCoulomb Stress\n+ Line 418: \"we calculated the triggering\"… this is actually not correct. Only the\nstatic stress changes were computed. It is then inferred/assumed that somewhat\nhigher stress near the hypocenter may have initiated / facilitated the 2nd event. But\nthere is dynamic rupture simulation carried out that actually shows that …\nYes, we only calculated the static stress changes caused by the Mw 7.8 event, and\nfound that the Mw 7.8 event stress change promoted the occurrence of the Mw 7.7\nevent. We rephrased the text.\nReferences\n+ several recent publications that also conducted finite-fault inversions are not\nlisted, but should be added (see above)\nWe added several recent publications that have appeared since our paper was\nsubmitted in the revised manuscript.\nFigures\n+ Figure 3 and 4, panels: from which fault-normal distance were aftershocks\nprojected onto the fault plane? Up to 5 km on either side of the fault plane? More or\nless?\nAftershocks less than 20 km on either side of the fault plane are projected onto the\nfault plane, but only the closest aftershocks are projected onto the fault plane for the\nintersecting faults. This is noted in the revised Fig. 3 caption.\n+ Figures 8 and 9 are confusing: they show time-slices (snapshots) of the rupture\nprocess, but in each time slices, the rupture-time contours of the entire rupture is\nshown. Shouldn't the time-slice 0-10 sec only contain contours for the first 10\nsecond?\nAlso, I find the color-range choice suboptimal, as tiny variations in shades of blue\ncannot be distinguished. Regions of zero slip should be shown in white, or light gray,\nand thus be clearly distinct from regions with 0.5-1.0 m. This is not case here.\nWe modified the figures as suggested in the revision.\nSupplementary Material\nFigure 3 is confusing here and should be removed. Maps related to the 2011 Tohoku\nearthquake are not needed here\nWe removed the original Supplemental Fig. 3 and added a link to a related database\nin the section of Code Availability.\nWe thank the reviewer for their very constructive comments, and the manuscript is\nimproved by our revisions to address them. This is noted in the revised\nacknowledgments.\n\nReviewer #3 (Remarks to the Author):\nAuthors presented a thorough analysis on the rupture processes and possible\nimpacts of the devastating Turkey earthquake doublet. The results suggest that both\nevents ruptured multiple fault segments, featured supershear rupture episodes, and\nthe first earthquake possibly aided the triggering of the second event. I like the idea\nof enriching the static displacement measurements using the near-field strong\nmotion accelerograms. However, I have some concerns with the validity of the\ndata product and results, as well as some unclear delivery of the method approach.\nAcknowledging that the implication of this work is important, here I suggest a major\nrevision, see the following points for details.\nMajor comments:\n1. The fault geometry is one major assumption and is critical to the finite fault\ninversions conducted in this study. How the fault geometries, particularly the dip\nangles, are constructed, is unclear, which raises concerns about the validity of the\nfindings. Additionally, the absence of satellite data and sources further compounds\nthe issue.\nWe elaborate on the choices made regarding fault dip, which does vary amongst\npublished models. Preliminary inversions with various strike/dip combinations for\nsparse fault models indicate typical limited direct waveform-based constraints on\nthe dip for steeply-dipping strike-slip geometries. Thus, it is desirable to constrain\nthe faulting geometry with a priori information to the extent possible. As discussed,\nwe use the InSAR mapping to prescribe the active subfault locations and outcrop\nstrikes along both ruptures. The distribution of seismicity with respect to the\nsurface ruptures indicates that the dip direction of the EAF changes along strike. A\ndetailed study of the aftershock activity of the Sivrice earthquake, just northeast of\nthe Mw 7.8 rupture, indicates an NW steeply-dipping fault (Güvercin et al., 2022).\nThis also holds for the northeastern part of the Mw 7.8 rupture based on cross-\nsections (D-D’, and possibly C-C’) of relocated aftershock as shown in new\nSupplementary Fig. 4. The aftershocks of the Mw 7.8 earthquake in the\nsouthwestern region (cross-sections A-A’ and B-B’) are mainly located southeast of\nthe surface fault trace, indicating a SE dip, in agreement with the aftershock focal\nmechanisms in this zone. The along-strike change in fault dip may have influenced\nthe bilateral rupture expansion. These relocated aftershock locations guide us in our\nchoice of dip direction and angle, but there is clearly some uncertainty.\n2. One of the major claims is the supershear stages found during both the northeast\nand southwest rupture episodes of the first event. However, how to validate\n(quantify the uncertainty) of such estimates? To date, most finite-fault inversions\nand back projection-based studies suggest subshear ruptures along the EAF during\nthe first event. These includes:\nMeng et al. (https://doi.org/10.21203/rs.3.rs-2747911/v1)\nMai et al. (https://doi.org/10.1785/0320230007)\nGoldberg et al. (https://doi.org/10.1785/0320230009)\nMelgar et al. (https://doi.org/10.26443/seismica.v2i3.387)\nAmong them, Meng et al. conducted a Mach Cone analysis arguing that the rupture is\npredominantly subshear.\nWe utilize a substantial amount of near-field strong-motion data, surpassing the\ndata coverage of previous studies, to constrain the slip evolution of the Mw 7.8\nevent. Extensive testing was conducted on the rupture speed, revealing a notable\ncontrast in the overall velocity of rupture propagation between the northeast (~ 3.5\nkm/s) and southwest (~2.5 km/s, but with the average rupture velocity on the aF5\nsegment being ~3.3 km/s, as depicted in Fig. 3a). This difference in rupture velocity\nwas inferred by waveform fitting analysis, as documented in the Supplement. The\nrelatively high rupture velocity observed in both directions is consistent with the\npreprint describing Meng's analysis and the teleseismic analysis of Mai et al. (2023).\nHowever, the BP analysis only provides an average rupture velocity for the\nnortheast and southwest directions, and does not resolve local variations in rupture\nvelocity. This is especially true if there are delays in rupture between some\nsegments, as appears to be the case for the Mw7.8 rupture, or if only a portion of the\nrupture is supershear (Mach wave coherence for surface waves will only be\npronounced if there are long intervals of supershear rupture, which we do not think\nis the case). The excellent strong motion data coverage right along the fault allows\nus to constrain local intervals with variable rupture velocity with reasonable\nconfidence.\n3. Authors' slip models have significant slip gaps along the northeast segment of the\nEast Anatolian Fault Zone, and the west of the Sürgü Fault. However, this appears to\ncontradict the satellite surface displacement data which are continuous for both\nevents (see Mai et al. Fig. 2; Goldberge et al. Fig. 4). How certain are these slip gaps?\nThe quickest way is to do a forward prediction of the satellite displacement field and\ncompare with the data.\nTo address the uncertainty and investigate its impact, we used different random\nseeds to perform ten inversions for the earthquake doublet. The results indicate that\nthe large-slip distributions among the ten models for the Mw 7.8 event display\nrelatively consistent behavior, demonstrating agreement among the models.\nGenerally, the standard deviation (STD) across most fault segments is negligible,\nexcept for the aF3 and aF6 segments (see revised Supplementary Fig. 12), where\nhigher variations are observed. Similarly, for the Mw 7.7 event, the STD is typically\nsmall compared to the average slip, with some exceptions found in the western\nparts of the bF1 and bF2 fault segments (revised Supplementary Fig. 12). The higher\nSTD in these cases corresponds to slip gaps. We suspect that the elevated STD in\ncertain areas of the fault segment is likely due to the lack of near-fault observations.\nThis finding suggests the need for future investigation in these specific regions. We\nnote that the slip models from teleseismic and geodetic data separately shown in\nMai et al. (2023) differ significantly, but both display segment variability\nqualitatively similar to our inversions for the Mw 7.8 event. Our model is also\nsimilar to the Goldberg et al. (2023) model, which has substantial non-uniformity\nalong strike. The along-strike variation in fault-parallel displacement inferred from\npixel-tracking indeed has a more continuous distribution, possibly influenced by\nconditioning effects, but there are still rapid spatial slip variations of ~3-4 m over\n40 km long sections, so there is substantial variability. Ongoing detailed mapping of\nsurface ruptures along the fault zone may shed light on the true spatial\nheterogeneity. Besides, we did no include the pixel-tracking offsets in the joint\ninversions due several factors, including low resolution near the fault, the\nindistinguishability of coseismic displacements from two earthquakes, and potential\nearly afterslip effects.\nOverall, our analysis highlights the stability and consistency of the large slip\ndistributions for this earthquake doublet, as inferred from the abundant seismic\nobservations and static offset measurements. We also point out the presence of\nlocalized variations and the need for additional research in certain areas where\nnear-fault observations are lacking.\n4. The conversion from strong motion accelerograms to static displacements\nimpressively enriched the near-fault static observations. However, I am concerned\nwith the robustness. For example, the static slip of station 4615 during the first\nearthquake has a 1.5 m amplitude and is perpendicular to the dominant strike slip\ndirection of the East Anatolian Fault Zone, whereas the satellite radar images do not\nshow any significant southeast motions (Mat et al, Goldberg et al.). Authors need to\nmake sure that such conversions are valid.\nWe compared the localized coseismic displacements derived from strong-motion\ndata and the horizontal displacements derived from pixel-tracking offsets of\nSentinel-1 satellite radar images (Mai et al., 2023), finding overall reasonable\nagreement. However, there are exceptions observed at near-fault stations, as\nillustrated in new Supplementary Fig. 3. Given the inherent uncertainties associated\nwith both approaches, such as the resolution of the pixels and the orientation error\nof the strong motion stations, these uncertainties inevitably contribute to\ndifferences in both magnitude and direction of the derived horizontal\ndisplacements. As for station 4614 and 4615 specifically, they exhibit high-quality\nstrong-motion recordings without significant baseline shift other than correctable\ndrift (see the following Figure), and has a strong fault-normal component for station\n4615 (relative to the initial splay fault orientation) which is very similar to that of\nthe pixel-tracking estimate (see revised Supplementary Fig. 3). A larger discrepancy\nbetween the estimates is seen at station 4614, but that site is likely influenced by\ncontributions from slip on both the splay and the main EAF, so it is not unreasonable\nto have oblique motion. Therefore, using our new processing method, we have\nconfidence in the stability and reliability of the static displacements derived from\nthe strong-motion data.\nFigure. Displacements integrated from strong-motion data at station 4614 and 4615\n(black) with baseline corrections (red). The green curves indicate the estimated\nbaseline effect.\n“The strong-motion derived coseismic displacements are generally consistent with\nhorizontal displacements derived from pixel-tracking offsets of Sentinel-1 satellite\nradar images37. However, there are some differences at near-fault stations, as\nillustrated in Supplementary Fig. 3. Given the inherent uncertainties associated with\nboth approaches, such as the resolution of the pixels and the orientation error of the\nstrong motion stations, these uncertainties inevitably contribute to differences in both\nmagnitude and direction of the derived horizontal displacements.”\n5. 0.1 MPa is still small for static triggering of the Mw 7.7 second earthquake. Is\nthere a source of minimal earthquake triggering threshold of ~0.01 MPa for M7\nearthquakes?\nPrevious studies have confirmed that the change in Coulomb stress value between\n0.01 and 0.1 MPa is thought to be enough to trigger an earthquake in the future,\nwith 0.01 MPa being the threshold value for earthquake triggering. We added a\nreference to Ross Stein’s work in the revision. Triggerability is, of course, influenced\nby proximity to failure stress, so that very small peak dynamic strains can trigger\nfailures (e.g., van der Elst and Brodsky, JGR, 2010).\nMinor comments:\n1. The colorbars for slip amount are severely saturated in Figures 3, 5, 8, 9. I can’t\ndistinguish anything between 6 and 12 m.\nWe changed the colorbar to clarify the slip magnitude.\n2. Still for these figures (especially Fig. 3 and 5), the slip models (subpanel a) are\ninversely oriented compared with the map (subpanel b). I had a hard time to realize\nthat the left-hand side indicates the faults to the east.\nWe revised Figs. 3 and 5 to clarify the fault orientation.\nWe thank the reviewer for their comments, and the manuscript is improved by our\nrevisions to address them. This is noted in the revised acknowledgments.\n\nREVIEWERS' COMMENTS\nReviewer #3 (Remarks to the Author):\nI would like to acknowledge that the authors have made a significant effort to address most of my\nconcerns. The added comparison between aftershock distribution and fault dip justifies the\nselection of fault geometry. The additional test on multiple random seeds reflects the variation of\nslip distribution in the current settings and assumptions. Comparison between strong motion-\nderived displacements and the pixel tracking offsets illustrates the overall robustness of the\nderivation. Authors also improved their figures. Therefore, I recommend publication of the\nmanuscript after evaluation of my following minor comments.\n1. The overall high rupture velocity of the first event is consistent with other studies. The part I\nfeel slightly concerned about is, the methodology from your slip inversion is essentially the same\nas Meng et al. and Goldberg et al., but the final kinematics has some (nontrivial) difference. Is this\ndifference a reflection of the model uncertainty as well? Is it due to the imposed difference of\nassumptions on fault geometries and delays? Author may add one sentence in the discussion as a\ncomment.\n2. Although 0.1 MPa can induce seismicity, triggering of large destructive earthquakes are still\nmuch more rare than the occurrence of stress perturbation of 0.1 MPa. Is it somewhat related to\nthe stochasticity of large earthquake nucleation? Author may add one sentence in the discussion\nas a comment.\n\nResponse to reviews of \"Complex multi-fault rupture and triggering during the February 6,\n2023, earthquake doublet in southeastern Türkiye\" by Liu et al., submitted to Nature\nCommunications. The review comments are reproduced below in black type, with our\nresponses and indications of how we have revised the manuscript to address the reviews\nindicated in blue.\n\nREVIEWERS' COMMENTS\nReviewer #3 (Remarks to the Author):\nI would like to acknowledge that the authors have made a significant effort to address most\nof my concerns. The added comparison between aftershock distribution and fault dip\njustifies the selection of fault geometry. The additional test on multiple random seeds\nreflects the variation of slip distribution in the current settings and assumptions.\nComparison between strong motion-derived displacements and the pixel tracking offsets\nillustrates the overall robustness of the derivation. Authors also improved their figures.\nTherefore, I recommend publication of the manuscript after evaluation of my following\nminor comments.\n1. The overall high rupture velocity of the first event is consistent with other studies. The\npart I feel slightly concerned about is, the methodology from your slip inversion is\nessentially the same as Meng et al. and Goldberg et al., but the final kinematics has some\n(nontrivial) difference. Is this difference a reflection of the model uncertainty as well? Is it\ndue to the imposed difference of assumptions on fault geometries and delays? Author may\nadd one sentence in the discussion as a comment.\nYes, there is variation among the published kinematic models, including one presented in\nthe new paper by Jia et al. that we now cite. It is challenging to specifically state exactly\nwhy there are differences, but in general, these will occur due to differences in precise\nmodel geometries (different dips are used for some fault segments between studies),\ndifferent data (for example, our study is the only to be using the expanded near-fault static\ndisplacement observations obtained from strong-motion instruments to constrain the\ncoseismic slip), along with differences in inversion algorithms and model parameterization.\nWe added a sentence to note these issues, as suggested.\n2. Although 0.1 MPa can induce seismicity, triggering of large destructive earthquakes are\nstill much more rare than the occurrence of stress perturbation of 0.1 MPa. Is it somewhat\nrelated to the stochasticity of large earthquake nucleation? Author may add one sentence\nin the discussion as a comment.\nIndeed, triggering intuitively requires pre-stress accumulation that has approached failure,\nand only then can ~0.01 MPa stress increments that are favorably oriented initiate rupture\nonset. The empirically estimated lower stress increment level of ~0.01 MPa increment may\nreflect observational limits, but is a threshold consistent with many observations. We add\na brief discussion of this issue as suggested."}]

Decision letter