Multi-scale rupture growth with alternating directions in a 1 complex fault network during the 2023 south-eastern 2 Türkiye and Syria earthquake doublet

16 A devastating doublet of earthquakes with moment magnitude M W 7.9 and M W 7.6 earth- 17 quakes contiguously occurred in SE Türkiye near the NW border of Syria. Here we per- 18 form a potency-density tensor inversion to simultaneously estimate rupture evolution 19 and fault geometry for the doublet. We find the initial M W 7.9 earthquake involved dis- 20 crete episodes of supershear rupture and back-rupture propagation

Our modeling strategy shares a similarity with that of seismic back-projection, which 164 requires very few assumptions about the fault geometry and rupture information (Ishii 165 et al., 2005;Y. Xu et al., 2009;Meng et al., 2012;Nissen et al., 2016;Satriano et al., 2012; 166 Yao et al., 2011;Taymaz et al., 2021). Our approach additionally provides kinematic in-167 formation by directly solving for the potency-rate density distribution, which should 168 enable in-depth evaluation of rupture dynamics that, for example, can be associated with 169 variable fault geometry. To perform a stable inversion with such a high degree-of-freedom 170 model without overfitting, the uncertainty of the Green's function is incorporated into 171 the data covariance matrix (Yagi & Fukahata, 2011) and the strength of smoothing is ad-172 justed using the Akaike's Bayesian Information Criterion (e.g., Akaike, 1980;Yabuki & 173 Matsu 'ura, 1992;Sato et al., 2022). We note the effect of structural heterogeneity can  gional tectonic regime (e.g., Fukahata & Wright, 2008). This exercise highlights the im-181 portance of permitting a complex rupture scenario with enough model freedom and over-182 constraining the model would fail to explain the seismic signals that are responsible for 183 the change of focal mechanism during rupture (e.g., Shimizu et al., 2020). Aftershock 184 focal mechanisms and moment tensors show a variability with a deviation from pure 185 strike-slip faulting (Fig. 1), helping to demonstrate that a flexible potency-density ten-186 sor approach is required. 187 We applied a standardized data processing workflow for our potency-density ten-188 sor approach that has been applied to earthquakes in different tectonic regimes (Shimizu 189 et al., 2020;Tadapansawut et al., 2021;Hu et al., 2021;Fan et al., 2022;Fang et al., 2022;190 Hicks et al., 2020; Yamashita, Yagi, Okuwaki, Shimizu, et al., 2022;Yagi et al., 2023). We 191 used the vertical component of teleseismic P -waveforms from a total of 39 and 37 sta-192 tions for the M W 7.9 and M W 7.6 earthquakes, respectively (Figs. S1 and S2). The data 193 were selected to ensure sufficient azimuthal coverage so that we can resolve potential 194 variations of radiation pattern during the rupture evolution and hence spatiotemporal 195 changes of fault geometry. We selected data so that we manually picked the first mo-196 tion of P -wave (e.g., Okuwaki et al., 2016). The data were then restituted to velocity at 197 1.0-s sampling interval by removing the instrumental responses. Green's functions were 198 calculated based on the method of Kikuchi and Kanamori (1991), adopting CRUST1.0 199 model (Laske et al., 2013) for the one-dimensional layered velocity structure around the 200 source region (Table S1). We further tested the robustness of our modeling against an 201 alternative structural model adopted from the ak135 model (Kennett et al., 1995) (Ta-202 ble S2), showing that the resultant pattern of potency-density tensors is less sensitive 203 to the choice of the near-source structure model (Fig. S7). The initial rupture point is 204 taken from the relocated epicenter for the M W 7.9 earthquake  and 205 on the model fault near the relocated epicenter for the M W 7.6 earthquake. We set the 206 hypocentral depth at 15 km for both earthquakes (Fig. S3). The uniformly-distributed 207 model source elements are regularly spaced 10×5 km and 5×5 km in the along-strike 208 and dip directions for the M W 7.9 and M W 7.6 earthquakes, respectively, along a ver-tically dipping non-planar model fault that aligns with the active faults (Emre et al., 2018) 210 and the relocated aftershocks (Fig. S3). Together with the curved main EAF strand, we 211 adopted a splay fault into our model fault centered on the initial rupture point, which 212 is oriented at 35 • NE, having an acute angle relative to the main EAF in NE direction 213 (Fig. S3). ). The total seismic moment is 9.6×10 20 N m (M W 7.9), which is similar to the that 220 estimated from coda waves (X. Jiang et al., 2023). The overall faulting mechanism in-221 dicated by the flexible potency density tensors is consistent with our prescribed non-222 planar model fault geometry (Fig. 2). The potency-density tensors show a largely pla-223 nar fault with depth. The space-time evolution of the rupture shows four distinct episodes 224 which we describe in the following paragraphs.

225
Rupture Episode 1. The first-motion faulting mechanism using local-regional wave-  Rupture Episode 2. After a relative quiescence for 5 s after the end of the first episode, 235 the second rupture episode starts at OT+15 s, lying 60 km NE of the epicenter. This episode 236 releases the greatest amount of seismic moment (35%; M W 7.6) of the entire rupture.

237
The rupture propagates in an asymmetric bilateral manner with a strong SW-oriented 238 direction, rupturing a total length of 120 km over 20 s duration. Most notably, the SW 239 flank of the rupture front apparently back-propagates through the hypocentral region 240 beyond 20 km SW of the epicenter (Fig. 2). The migration speed of the associated SW-241 directing back-propagating rupture signal exceeds the local S-wave velocity (Table S1;  We note that the source elements with minor potency rate may be affected by the sur-rounding major potency rate due to smoothing effects, so we do not interpret the resul-252 tant strike angle from those minor potency-rate tensors. at 120 km NE from the epicenter (Fig. 2). The strike orientation is similar to that of Episode 260 2 and remains consistent with the main EAF. We refrain from measuring rupture speeds 261 for this episode as they seem sensitive to the assumption of maximum slip duration (   The initial rupture of the M W 7.9 event has a different fault orientation than that 288 of the following main bilateral rupture that releases most (97%) of the seismic moment.

289
For example, during the peak slip of the first rupture episode (7-8 s), the strike is 36 • , 290 whilst the later bilateral rupture episode has a strike of 55 • (Fig. 3). Intense aftershock 291 activity is observed NE of the epicenter , in a lineation oriented SW to NE, seemingly connecting to the main EAF strand (Fig. 3). The alignment of these af-293 tershocks on the splay fault is consistent with the strike estimated from our inversion.

294
To the east of the epicenter, the Narlıdağ fault zone has been mapped to extend to the 295 N and NE (Perinçek & Çemen, 1990;Duman & Emre, 2013). From rapid analyses of the 296 satellite images and field measurements, surface rupture is also observed near the epi-297 center, which is elongated NE and is consistent with our estimated strike orientation (Reitman 298 et al., 2023), which is called as Nurdağı-Pazarcık fault by . Thus, the 299 first rupture episode occurred on a sub-parallel splay fault to the main EAF. Although 300 our potency-density tensor inversion finds mostly pure strike-slip faulting during the 301 first rupture episode, the first-motion mechanism from near-field waveforms suggest 302 that the rupture initiated with a weak phase of oblique-normal faulting (Fig. 3c), which 303 is likely too small to be resolved in teleseismic waveforms. From our estimated strike 304 orientations, the angle between the splay fault and the main EAF model domain is ∼18 • , 305 which is close to the peak of the splay fault angle distributions (±17 • ) that was previ-306 ously observed for active faults in California (Ando et al., 2009;Scholz et al., 2010). In 307 between the first and second rupture episodes, we only see minor moment release, which 308 may suggest a non-continuous rupture at the junction between the splay fault and main 309 EAF. However, due to the insufficient spatial resolution of the teleseismic data we used,  tions and the hypothesized rupture-front speed (Figs. S10 and S11). Such a boomerang-319 like back rupture propagation is an end-member rupture behavior that has become more 320 frequently reported with higher-resolution datasets and more detailed rupture imag-321 ing (Meng et al., 2018;Hicks et al., 2020;Vallée et 322 al., 2023). However, because the earthquakes in all of these cases studied were either 323 deep or in remote areas, there were no surface rupture observations that could have ex-324 plained the apparent back-rupture propagation. Therefore, the apparent boomerang rup-325 ture of the 2023 SE Türkiye earthquake is intriguing because we show that the rupture 326 propagated along different sub-parallel fault strands which could offer an mechanism 327 for these previously reported examples of back-propagating ruptures.

328
Although it is still difficult to find a deterministic explanation of why the initial 329 rupture occurred on the more minor bifurcated fault rather than the main EAF, the se- to be ruptured once assisted by the initial rupture on the bifurcated fault. Although our 337 sole use of teleseismic data may not rigorously discriminate the absolute location of the 338 slip on the closely located parallel faults, we favor that the apparent back-propagating 339 part of the rupture occurred on the main EAF because of the higher potency rate on the 340 main EAF model fault rather than on the splay model fault (Fig. 2c,d). This assumption 341 is supported by independent modeling using geodetic datasets that finds larger slip along

366
The strike orientation during the second rupture episode (OT+15-20 s) is slightly 367 rotated clockwise, which is also mapped in the main EAF strand west of the junction 368 ( Figs. 1 and 3). If this change in fault orientation acts as a restraining bend given the back-369 ground stress field, the rupture propagation may cause a concentration of stress at the 370 bend. This might have caused the rupture deceleration, which can be seen as the slip 371 stagnation during OT+15-20 s. Soon after this pause, dynamic stresses allowed the rup-372 ture to continue and propagate to the SW and even briefly accelerate its speed, which 373 can be consistent with the predicted behavior of a supershear rupture transition across 374 restraining bends (e.g., Bruhat et al., 2016). We emphasize here that our source model 375 does show that the Mw 7.9 earthquake is not supershear throughout the entire event, 376 but it involves discrete supershear along certain fault segments during each rupture episode.

377
Such discrete supershear pulses have been independently estimated using near-field records 378 (e.g., Delouis et al., 2023) and numerical simulations (e.g., Abdelmeguid et al., 2023). 379 We further note that the NE and SW boundaries of the second rupture episode co-380 incide with mapped fault steps near Gölbaşı and south of Nurdağı (see locations S1 and S2 in Fig. S5). Such steps may contribute to the apparent gaps of 10 s between the sec-382 ond and subsequent rupture episodes (Fig. S5). We do not have enough evidence to ex-383 plain how such gaps are physically connected, but our finding will stimulate further re-384 search to investigate how the rupture evolved across fault steps, for example, the long 385 nucleation processes or possibly inter-subevent slow deformation.

397
The strike extracted from the best-double-couple solution of our estimated potency-398 density tensors is not apparently aligned with the bulk linear trend of the active faults 399 (Fig. 2). However, because we observe non-double-couple fractions for the SW end rup-400 ture (e.g., 24% during 60-61 s; Fig. 3  We find the M W 7.6 earthquake shows a much more focused rupture process, com-414 pared with the preceding M W 7.9 event. Yet, our solution finds that the strike of the rup-415 tured fault geometry curves gradually, with a counterclockwise rotation toward the west.

416
The rotation trend can favorably be oriented to the optimal plane of the background hor-    . The blue beachballs are the GCMT solutions (Dziewonski et al., 1981;Ekström et al., 2012) Melgar et al. (2020). The stars, dots, and lines are the same as shown in Fig. 1.