Introduction

The destruction that unfolded in southeast Turkey and northwest Syria after the 6 February 2023, earthquake doublet (Fig. 1) was devastating. The first M w 7.8 earthquake was more than twice as large as the most significant known regional earthquakes, which had reached up to M w 7.4 (Fig. 1a; Duman and Emre, 2013), and is the most powerful earthquake recorded in Turkey since 1939. It initiated south of the Eastern Anatolian fault (EAF) on a splay fault known as the Nurdağı-Pazarcık fault (NPF) before branching northeast and southwest and bilaterally rupturing > 300 km of the EAF (Fig. 1a; Barbot et al., 2023;Goldberg et al., 2023;Jia et al., 2023;Karabacak et al., 2023;Mai et al., 2023;Melgar et al., 2023;Okuwaki et al., 2023). The EAF did not host significant earthquakes during the last century until the 2020 M w 6.8 Elaziğ earthquake (e.g., Cakir et al., 2023), located north of the first rupture. Nine hours after the first earthquake, a second M w 7.7 earthquake ruptured the geometrically complex Çardak fault network ≈60 km north of the NPF-EAF junction. The Çardak fault is a part of the predominantly strike-slip Sürgü-Çardak-Savrun fault (SCSF) system, although the Sürgü segment, which connects the Çardak fault to the EAF, did not rupture.

The sequence ruptured segments of the northeast-striking EAF-a major intracontinental left-lateral strike-slip fault that accommodates northward convergence of the Arabian of the Arabian plate and westward motion of the Anatolian plate, Goldberg et al., 2023). Shaded areas show the inferred extent of historic surface ruptures labeled by year and magnitude (Duman and Emre, 2013). Red and blue numbers correspond to fault segments modeled in this study named following Duman and Emre (2013). The first earthquake is modeled using six segments of the EAF: 1 and 2, Amanos segment; 3, Pazarcık segment; 4, Nurdağı-Pazarcık fault (NPF); 5, unnamed Erkenek splay; 6, Erkenek segment. The second earthquake ruptures four segments of the SCSF: 7, Çardak fault; 8, Göksun bend segment; 9, Malatya fault; and 10, unnamed Göksun splay. The Sürgü fault (segment 11) is shown in Figure 5. Inset shows regional tectonic map modified from Barbot and Weiss (2021). The black rectangle outlines the area shown in the bottom panel. (b) Bottom: zoomed view of East Anatolian fault zone principal strain rate directions in purple (first component) and pink (second component) from Weiss et al. (2020). In dark and light gray, we show the seismologically inferred maximum and minimum principal horizontal stress components from Güvercin et al. (2022), as well as in dark and light blue, the maximum and minimum principal horizontal stress orientations used in this study. (c) Initial conditions for 3D dynamic rupture modeling of both large earthquakes. SH max ° is the orientation of the maximum horizontal compressive stress from a new stress inversion we perform (based on Güvercin et al., 2022, Fig. S5), D c is the critical slip-weakening distance in the linear slip-weakening friction law, R 0 is the maximum relative prestress ratio, and R < R 0 is the fault-local relative prestress ratio modulated by varying fault geometry and orientation. Although the assumed SH max is the same in distinction to the dynamic rupture models in Jia et al. (2023), no additional smaller scale initial prestress or fault strength heterogeneity is prescribed.

resulting in geologic and geodetic slip rates of up to 10 mm/yr across the complex EAF system (inset of Fig. 1a, Taymaz et al., 2007;Duman and Emre, 2013;Weiss et al., 2020). Ongoing approximately north-south Arabian-Eurasian collision squeezes and extrudes the Anatolian plate westward (Fig. 1a), driving westward migration and counterclockwise rotation of eastern Anatolia relative to Eurasia, along with transpression that is accommodated partly by intraplate distributed deformation and complex faulting, but primarily by slip along its major boundary faults: the right-lateral North Anatolian fault (NAF) to the north and the EAF to the southeast (e.g., Barbot and Weiss, 2021). The eastern SCSF branches to the west from the EAF and curves south to follow the EAF subparallel into the Gulf of Iskenderun to the southwest (Fig. 1a), where the SCSF links into the Cyprean Arc and a strand of the EAF continues south into the Dead Sea Transform fault (DSTF). Before hosting surprisingly large left-lateral coseismic slip on 6 February, the sense of offset along the SCSF remained debated (e.g., Koc and Kaymakcı, 2013;Duman and Emre, 2013), with fault-bounded geomorphologic features and slip orientations measured from fault scarps leading some to infer dextral motion (Koc and Kaymakcı, 2013). The tectonic and structural complexity of the associated fault systems reflects the complex modern and paleotectonics of the region (Fig. 1a,b), highlighting how strain partitioning across distributed networks of nonuniformly oriented fault segments can accommodate sharp lateral variations in local tectonic loading (Duman and Emre, 2013;Weiss et al., 2020;Barbot and Weiss, 2021;Güvercin et al., 2022). The EAF is considered relatively immature compared to the NAF, initiating 2-5 Ma ago and since accruing 22-33 km of offset (Saroglu, 1992). Interferometric Synthetic Aperture Radar (InSAR) and Global Navigation Satellite System (GNSS) data show slightly elevated principal strain rates around the SCSF-EAF that are above interior Anatolian background levels but far lower than those along the NAF (Fig. 1b; Weiss et al., 2020). Principal strain rate directions rotate ≈30°along the EAF from north-northwest near the Gulf of Iskenderun to north-northeast around the SCSF-EAF junction and northeast further along strike (Fig. 1b). Principal horizontal stress directions inverted from focal mechanisms of nearby earthquakes show similar ≈20°-30°clockwise rotations along the SCSF-EAF, from primarily north-northeasttrending in the southwest to northeast-trending near their intersection and further northeast (Fig. 1b). Similar to other recent significant earthquakes, such as the 2016 Kaikōura, New Zealand, earthquake and the 2019 Ridgecrest, California, sequence, the 2023 Turkey earthquakes activated more fault segments than expected from geodetic slip rates and historical earthquakes, due to the geometrically complex fault structures interacting across space and time scales. These earthquakes illustrate the difficulties in reliable estimates of expected earthquake magnitudes due to short-term records, irregular cycles, and multifault rupture dynamics (Goldfinger et al., 2013;Milner et al., 2022).

Multiple slip episodes occurring close in time and activating fault segments nearby challenge data-driven analysis of large earthquakes. Initial imaging and data-driven modeling efforts based on strong motion, teleseismic, and high-rate GNSS data reveal dynamic and structural complexity of both the events, highlighted by opposing interpretations of fault system interaction and the characteristics of each earthquake, for example, the inferred local rupture speeds, activated fault segments, and seismic moment release. Joint data-driven interpretation using various geophysical and geologic datasets can illuminate the spatiotemporal evolution of the rupture sequence (e.g., Jia et al., 2023) but typically cannot probe dynamically viable pre-and coseismic mechanical conditions. Detailed, physicsbased interpretations can pose a unifying approach capable of explaining independent data sets and unifying unexpected field observations but are often only available on timescales of months to years after large earthquakes (e.g., Taufiqurrahman et al., 2023). The complex interactions between various fault segments during multifault earthquake sequences can lead to more significant damage, complex seismic and geodetic data signatures, and modeling challenges, because they require considering the fault geometry and slip characteristics of multiple faults simultaneously.

We present 3D dynamic rupture models based on seismic and geodetic observations available within hours to days after the earthquakes. Our models are not initialized using solutions to inverse problems, such as a static slip model in the different simulations presented in Jia et al. (2023). We provide a new comparison of the fault offsets measured from Sentinel-2 displacement fields with the dynamic rupture models and compare the timing of the observed peak ground motions at unclipped strong-motion stations. Our models also explain the northward deviation of the second earthquake's western rupture and the minimal slip observed on the Sürgü segment. Our physics-based simulations disentangle the complex set of observations from the earthquake doublet and their interrelationship with general implications for the often-underestimated hazard from multifault rupturing earthquakes. Furthermore, we highlight how quickly developed 3D dynamic rupture modeling constrained by observations available before and soon after complex ruptures can provide timely insight into the postseismic stress, strength, and rheological conditions of such fault systems, complementing postevent data collection and informing hazard estimation efforts.

Methods

In data-driven kinematic models, a large number of free parameters enables close fitting of detailed observations, often at the expense of mechanical consistency and uniqueness of the solution. Dynamic rupture modeling involves simulating how earthquakes nucleate, propagate, and arrest. 3D models can directly reproduce geophysical and geologic observables, such as seismic waves, geodetic deformation patterns, and surface rupture patterns, in a physically self-consistent manner. Computational advances enable large-scale 3D dynamic rupture modeling informed by observations in various tectonic and scientific contexts. Such models require prescribed initial conditions, including fault geometry, fault strength, initial stress distribution, and material properties (Ramos et al., 2022). Our 3D dynamic rupture simulations use the same spatially variable fault geometries, ambient tectonic loading, and fault strengths in a mechanically linked model of both the earthquakes that captures the respective effects of static and dynamic stress transfers during and between both the earthquakes. The basics of the dynamic rupture model setups and their computational cost (Fig. S7, available in the supplemental material to this article) are equivalent to the refined, more heterogeneous models of Jia et al. (2023). Our simulations allow us to analyze how the earthquakes are related by analyzing the conditions that contributed to the doublet occurrence, including fault geometry and mechanical properties. We base our model's initial conditions on a few parameters for which we can constrain their spatial variability from those observations available before and soon after the earthquakes, including fault geometry from space-geodesy and seismicity (see the Fault geometries section), fault loading from regional seismotectonics (see the Initial stresses and fault loading section), and dynamic parameters from observed earthquake kinematics such as moment rate release (see the Fault friction, offfault plasticity, and dynamic fault strength section).

Fault geometries

Our fault geometry includes 10 curved, intersecting segments with variable dip (Fig. 1c). We constrain the geometry of the fault system from rupture traces mapped from coseismic horizontal surface displacements from pixel correlation of 10 m resolution Sentinel-2 satellite imagery (Table S3; Fig. S9; see Data and Resources). We show the coseismic east-west and north-south displacement fields and measured fault offsets in Figure 2. Our fault geometries capture large-scale geometrical complexities of the fault system, including fault bends, stepovers, and secondary segments. We extend the mapped surface fault traces to a depth of 20 km with varying dip angles ranging between 90°for the first earthquake's segments and 70°N for the main segments of the second earthquake, which is likely simplified but to first-order consistent with relocated aftershocks (Lomax, 2023). The small Göksun splay (segment 10) dips 90°. The revisit frequency of the Sentinel-2 optical or Sentinel-1 Synthetic Aperture Radar (SAR) satellite constellations is 5 and 6 days, respectively, implying that the surface rupture trace of a large continental earthquake can be constrained first order within this timeframe (see Fig. S9 for the optical displacement field we obtained 3 days after the doublet). We explore alternative models incorporating the Sürgü fault, which connects the SCSF to the EAF, with geometries constrained from geologic mapping (Emre et al., 2018) in the Discussion section.

Initial stresses and fault loading

We expose all fault segments to depth-dependent and laterally rotating initial stresses resembling the regional state of stress in the upper crust and regional stress inversion (Fig. 1c, Fig. S5). We combine the Mohr-Coulomb theory of frictional failure with dynamic parameters to reduce the ample parameter space of dynamic rupture modeling (e.g., Taufiqurrahman et al., 2023). We construct a 3D Cartesian stress tensor S ij within the model domain (Fig. S7) from a new stress inversion we perform (based on Güvercin et al., 2022, Fig. 1b, Fig. S5) for the orientation of the maximum horizontal compressive stress SH max and the stress shape ratio ν. The latter balances the amplitude of principal stress components and is here ν = 0.5, reflecting a strike-slip regime. We assume depth-dependent effective normal stress and overpressurized pore fluids. Above 6 km depth, we assume a pore fluid pressure ratio of γ ρ water =ρ 0:66. At a larger depth, the pore fluid pressure gradient mirrors the lithostatic stress gradient, leading to constant effective normal stress below 6 km depth (Rice, 1992). There are few constraints on SH max near the seismically quiet Çardak fault (Fig. 1b). However, using a few trials of dynamic rupture simulations, we find that assuming along-strike rotating loading with Sentinel-2 displacement fields across the (c) M w 7.8 and (d) M w 7.7 ruptures with the fault offsets measured from the dynamic rupture models. See also Figure S9. (g,h) Comparison of observed (orange) and dynamic rupture modeled (blue) horizontal components of Global Navigation Satellite System (GNSS) displacements for the (g) M w 7.8 and (h) M w 7.7 earthquakes. The vector error ellipses represent a confidence interval of 95%.

close-to-optimal SH max close to the hypocenter of the second earthquake (Fig. 1c) is required to dynamically generate the observed large slip and surface displacements.

The relative prestress ratio R of potential dynamic stress drop to full frictional breakdown strength is a crucial parameter controlling earthquake dynamics and dynamic triggering potential. We define

in which τ 0 and σ n are the initial shear and normal tractions on a fault plane, and μ d and μ s are the static and dynamic friction coefficients, respectively. R governs the ratio of strain energy and fracture energy and, thus, determines local stress drop and acceleration or deceleration of the rupture front. Prescribing the maximum R 0 ≥ R, as the relative prestress ratio of an optimally oriented fault portion, allows us to constrain deviatoric stress amplitudes by balancing the amplitude of prestress relative to the frictional strength drop. We find that the observed complexities of the first and the second earthquake require R 0 to be variable within the model domain.

The relative prestress ratio of an optimally oriented fault in the complex stress field, which is first-order aligned with regional rotations of principal stress directions, ranges between R 0 0:3…0:8 (Fig. 1c, Fig. S5). Local fault geometry further modulates fault prestress and strength, resulting in heterogeneous fault-local R ≤ R 0 , and implying that locally more optimally oriented fault portions are closer to critically prestressed. Our parameterization satisfies the dynamic constraints that the second earthquake's fault structure is not dynamically triggered during the first earthquake rupture in our linked simulations. As we show in the Discussion section, these tectonically constrained prestress and fault strength conditions also explain why the Sürgü fault did not coseismically slip during the second earthquake, without requiring choosing locally different conditions. Importantly, we must prescribe larger fracture energy and larger nucleation energy in addition to dynamically stronger main faults for the second earthquake to capture their distinct rupture dynamics, as explained in the next section.

Fault friction, off-fault plasticity, and dynamic fault strength We demonstrate that simple friction parameters can give rise to the distinct slip characteristics of both the earthquakes, and promote dynamic and static multifault earthquake cascading in complex tectonic contexts. We use the widely used linear slip-weakening friction law (Andrews, 1976), which we parameterize with static friction coefficient μ s 0:6 and dynamic friction coefficient μ d 0:2 (Table S2). The critical slip distance varies between D c 0:5 m for faults hosting the first earthquake and larger D c 1:0 m, implying larger fracture energy, for the main faults hosting the second earthquake. All faults will begin to slip when shear stresses locally exceed the high-static frictional fault strength. Fault strength then decreases linearly from static to dynamic levels over the critical slip distance. We assume a depth-dependent, nonassociated Drucker-Prager elasto-viscoplastic rheology to model coseismic off-fault plastic deformation. Off-fault plasticity is parameterized by bulk internal friction coefficient and 1D variable plastic cohesion. We use a uniform bulk friction coefficient of 0.6, matching our on-fault static friction coefficient, and define plastic cohesion C plast as proportional to the 1D depth-dependent shear modulus (Table S1).

We link both the dynamic rupture earthquake models in the same simulation to account for the dynamic and static stress changes. To this end, we initiate the second earthquake at 150 s simulation time after initiating the first earthquake to ensure that no significant wave-transmitted dynamic stresses of the first earthquake interfere with the second dynamic rupture simulation. We do not account for potential interseismic aseismic slip or healing. We use nucleation patches that grow smoothly in time and across a minimal-sized perturbation area, reproducing rupture kinematics determined in several trial dynamic rupture simulations. The nucleation radius must be chosen larger for the second earthquake than the first, despite the second fault system being closer to critical prestress levels: Nucleation patch sizes are 2 km and 3 km for the M w 7.8 and 7.7 earthquakes, respectively.

Results

Our linked dynamic rupture scenarios of the M w 7.8 and 7.7 doublet are each dynamically consistent with respective seismic and geodetic observations and the observed rupture patterns and fault-system interaction. We summarize the modeled rupture dynamics of both the earthquakes in Figure 3, and compare them with geodetic and near-field seismic observations in Figures 2 and 4. Although the slip distributions resulting from these simpler models are smoother, rupture dynamic results are generally consistent with more heterogeneous models (Fig. S6; Jia et al., 2023), demonstrating the robustness of our results and the extent to which broad crustal mechanical conditions established by regional seismotectonics may control notable features of the rupture sequence evolution.

Rupture dynamics of the M w 7.8 earthquake The M w 7.8 dynamic rupture scenario shows three distinct rupture phases. Rupture is artificially nucleated (Fig. 3) and starts as a subshear crack-like rupture on the NPF-an north-northeast-striking splay fault south of the EAF. After this first rupture episode, the geometric barrier formed by the NPF-EAF intersection nearly arrests the rupture (Fig. S2), leading to temporarily negligible seismic moment release rates. This delay is stable for different initial conditions (i.e., is also observed in alternative rupture models in Jia et al., 2023 and Fig. 5) and corresponds to a pronounced trough in teleseismically inferred moment rate release at 18 s rupture time (Fig. 3d) Time (s) (a) (b) (c) (d) unclipped strong-motion stations that recorded the (c) M w 7.8 and (d) M w 7.7 earthquakes. The size of the circles indicates the PGV value, and the colors indicate the time at which the PGV occurred.

The simulated waveforms in the dynamic rupture models resolve frequencies of at least 1 Hz close to the fault system (Fig. S7). PGV is here computed as PGV x PGV y p . For a color-coded comparison of PGV amplitudes and quantification of the differences in PGV timing and amplitudes, see Figure S8. We account for topography, viscoelastic attenuation, and off-fault plasticity but use a 1D model of subsurface structure (Table S1). and a delay in observed peak ground velocities (PGVs) after rupture onset (Fig. 4c, Fig. S8).

Crack-like rupture branching in the forward direction (toward the northeast) dominates the second rupture phase, whereas rupture branching to the southwest is dynamically unfavorable and causes a pronounced delay before the third rupture phase. The forward direction branch of the EAF is oriented at ≈145°to the NPF, which is close to optimal rupture branching angles (Kame et al., 2003). Although forward-direction branching to the northeast is expected to be mechanically favorable, direct backward branching of the same slip style as the main rupture (i.e., here, left-lateral) is unexpected and theoretically unfavorable; the main rupture is expected to induce shear stress in the opposite rake direction on the backward branch, impeding rupture there (e.g., Poliakov et al., 2002;Fliss et al., 2005). Therefore, backward branching at ≈35°is discouraged under all prestress conditions in theoretical and numerical analysis accounting for the dynamic stress field at the rupture front. However, backward branching has been observed in real earthquakes such as the Landers earthquake, the Hector Mine earthquake, and the Romanche transform fault earthquake. To explain these observations, complex mechanisms, including rupture jumping, delayed reactivation of weakened faults, free-surface interactions, and 3D prestress or locking heterogeneities have been suggested (Kame et al., 2003;Oglesby et al., 2003;Fliss et al., 2005;Wollherr et al., 2019;Hicks et al., 2020). Our findings agree with these previous static and dynamic analyses in that immediate bilateral dynamic rupture branching to both sides of the EAF is dynamically unfavorable. However, we show that 3D dynamic effects, including progressive unclamping, transient shear stressing, and static stress build-up at the fault intersection and the southwest branch itself (Fig. S2) due to the propagating northeast rupture, eventually lead to delayed and self-sustained branching toward the southwest. This mechanism is simpler than rupture jumping, free-surface effects, or 3D prestress variability.

In the third phase of the dynamic model, two bilateral pulse-like rupture fronts, with strong directivity and median rupture speeds of 3 km/s, unzip the EAF in both the directions. Figure 3b illustrates the modeled fault slip reaching up to 8.2 m on the Erkenek segment. We observe shallow dip-slip components due to dynamic rake changes enhanced near the free surface where confining stresses are low. Rupture speed remains overall subshear during the first earthquake.

Our dynamic rupture scenario reproduces various geodetic observations, including SAR and GNSS offsets, and optical correlation data of the first earthquake (Fig. 2a,b,c,e,g). In particular, the predicted fault offsets fit observations remarkably well (Fig. 2c), demonstrating that our model matches the surface slip amplitude variations due to fault system segmentation. The model also reproduces the timing and amplitude of the impulsive strong ground motion signals recorded at stations along the fault to first order (Fig. 4a), which are generated due to a combination of surface rupture, partially pulse-like rupture, and strong directivity. Our models show that intense near-field pulses (Fig. 5, e.g., station 4615) can be caused by surface-breaking subshear strike-slip rupture as well as by supershear rupture (see the Discussion section).

Rupture dynamics of the M w 7.7 earthquake Rupture of the M w 7.8 earthquake does not coseismically trigger sustained coseismic slip along the second earthquake's fault system (Fig. S1a). However, the static stress changes due to the first earthquake bring its central segments closer to failure while considerably shadowing its eastern regions (Fig. S1b). Despite the compound rupture extent, our modeled rupture dynamics are strongly asymmetric, and feature westward supershear rupture speeds and eastward subshear rupture (Fig. 3c). In comparison to the first earthquake, the modeled fault slip (up to 9.4 m, Fig. 3b) and stress drop (Fig. S3) are considerably larger, whereas peak slip rates are lower, reflecting differences in prestress and fracture energy (D c , Fig. 1c) between these scenarios. Figure S4 shows an alternative simulation of the second earthquake without accounting for the change in Coulomb failure stress (ΔCFS, Fig. S1b) of the first earthquake. In the east, where the stress shadow of (negative) ΔCFS is the largest, the second earthquake's rupture speed, slip amplitude, and peak slip rate are considerably increased (Fig. S4c). In distinction, fault slip is smaller on the central segment for which ΔCFS is positive (9.2 m). The static stress changes due to the first earthquake moves the second earthquake's eastern segments away from failure, aiding asymmetry between westward and eastward rupture dynamics, and increasing peak slip.

As in the preceding M w 7.8 earthquake, the combined dynamic rupture model synthetics closely match the Sentinel-2 surface displacement field (Fig. 2a) and the eastwest fault offsets (Fig. 2d) of the second earthquake. Both the dynamic rupture earthquake scenarios also largely match the amplitude of the horizontal displacement direction at site EKZ1 (the GNSS station closest to the second earthquake) with slight differences in orientation (Fig. 2g). The mismatch for the eastward north-south offsets of the second earthquake may imply that either the assumed dip angle of this fault segment is too shallow, the modeled dip-slip component is too large, or both.

Few near-source recordings are available for the second earthquake, challenging model verification with strong ground motion. Our synthetics capture the shape and amplitude of recorded ground-motion pulses reasonably well (Fig. 4b). The fit with the timing and amplitudes of observed PGV is good (Fig. S8), which supports our fault stress and strength assumptions differing for the second earthquake: assuming larger fracture energy and more critically stressed faults lead to larger stress drop (Fig. S3) but does not cause larger PGVs. In distinction, the high but subcritical EAF prestress, governing the first earthquake results in delayed backward branching, and rupture are unable to immediately activate the fault system of the second earthquake.

Discussion

Unexpected doublet dynamics explained by stress and strength heterogeneity The dynamic stressing around a propagating rupture is enhanced with increasing rupture velocity and may drive rupture across steep branching angles, as discussed in Kame et al. (2003), for varying sub-Rayleigh rupture speeds. Motivated by this, Figure 5a-e explores dynamic rupture scenarios featuring supershear versus subshear rupture speeds on the NPF during the M w 7.8 earthquake. Both overcome the geometric barrier (Fig. S2) and are difficult to distinguish in near-fault ground-motion synthetics.

Although the supershear Mach front promotes rupture jumping, the dominant driving factor aiding subshear rupture branching from the NPF to the EAF is the abrupt stopping of surfacebreaking strike-slip rupture at the geometric barrier (Fig. S2, Fliss et al., 2005). In both the cases, the abrupt stopping of the first rupture at the intersection of the hypocentral NPF with the EAF causes a clear delay in the moment release rate and aids rupture branching to the northeast of the EAF. The varying mechanisms involve equally delayed activation of a shallower part of the western EAF. We do not explore models featuring local eastward supershear rupture in the second phase of the first earthquake, at the EAF after the NPF-EAF junction (Wang et al., 2023), which may improve moment rate release match at 25 s rupture time.

The M w 7.7 rupture deviated northward in this model (and in Jia et al., 2023), breaking the Malatya fault (Fig. 5g), instead of connecting across the Sürgü or Doğanşehir fault to the EAF (as suggested by Melgar et al., 2023). We include this fault as an 11th segment in an alternative scenario of the second earthquake, keeping other model parameters unchanged. Surprisingly, the second rupture does not trigger the now included Sürgü fault coseismically (Fig. 5f-h), and the dynamics of the sharp deviation are unaffected, despite static stress changes favoring a straight rupture path (Fig. S1b, Jia et al., 2023), and without added heterogeneity. Segment 11 is poorly aligned to our regional stress field resulting in a lower subcritical prestress (R ≤ 0.2) than the Malatya fault. The large coseismic stresses from the propagating rupture and radiating seismic waves (Fig. S1) cannot overcome these barriers nor trigger even minor-to-moderate slip.

Larger fracture energy and more critically stressed faults of the SCSF system (Fig. 1, Fig. S5) are required to model the compound, supershear rupture, and higher stress drop of the second earthquake (Fig. 3, Fig. S3). This may reflect dynamically triggered fault valving due to upwelling fluids (Sibson, 1992) or locally more immature fault structures than the EAF, however, supershear rupture is generally assumed to be favored by more mature faults (Perrin et al., 2016). Distributed faulting and deformation across subparallel strands and splays of the SCSF-EAF suggest that these complex fault systems may have initiated and developed in response to relatively recent changes in regional tectonic loading, compared to potentially steadier loading responsible for long-lived and highly localized slip on the northeast portion of the EAF. Although the SCSF includes reactivated Miocene structures (e.g., Duman and Emre, 2013), its diffuse and complex geometry suggests that it may be effectively less mature than the well-developed EAF.

Regional strain partitioning across complex faults may predispose local and complex fault-to-fault stress interactions that prevent steady, long-term slip (Fletcher et al., 2016) but favor cascading earthquake dynamics and sequences. An evident lack of seismicity at the SCSF after the 2020 Elaziğ earthquake as well as after the 2023 M w 7.8 earthquake (Lomax, 2023) supports differences in dynamic fault strength. As a result, both the rupture scenarios feature surprisingly different and distinct earthquake dynamics that are difficult to incorporate into earthquake hazard assessment: The M w 7.8 grew larger than expected due to the dynamic breaching of geometric barriers, including branching into steep backward direction considered mechanically unfavorable. The M w 7.7 scenario features local supershear rupture on a dynamically stronger fault segment than those hosting subshear rupture during the first earthquake.

Unlocking the keystones of multifault ruptures

Because our dynamic rupture simulations necessitate relatively high but subcritical initial prestress on the EAF to correspond with the observed slip patterns in the first earthquake, we hypothesize that the EAF could be an "imperfect keystone" fault for southeast Turkey (Fletcher et al., 2016), governing the overall strength and eventual multisegment failure of a complex fault system. Smaller intersecting or neighboring faults, such as the NPF, may reach critical stress states from tectonic loading before the keystone fault, subsequently dissipating strain energy through minor slip during the keystone fault's late interseismic period and remaining "pinned" and unable to rupture until the keystone fault nears failure.

Our models illustrate the predisposition of complex fault geometries, prevalent in tectonically complex immature fault systems, for cascading multifault and multievent earthquake sequences. The Turkey-Syria doublet and other recent large earthquakes involving multifault rupture sequences highlight how understanding their underlying mechanics is crucial to improving earthquake hazard assessment and mitigation. Geometrically complex fault systems are often immature or located in regions of complex or recently reoriented tectonic loading, implying that the complexity of our presented scenarios may not be unusual. Identifying and characterizing potential keystone faults conjointly with local prestress may help improve hazard estimates for complex fault networks and constrain key inputs to dynamic rupture simulations.

Toward "rapid" dynamic rupture modeling Our results demonstrate that regionally constrained and simple dynamic rupture scenarios are capable of matching important characteristics relevant for earthquake engineering: our physics-based earthquake scenarios resolving up to 1 Hz seismic-wave propagation (Fig. S7) match the timing of observed near-fault PGVs and near-fault ground-motion velocity pulses (Fig. 4, Fig. S8). We do not solve inverse problems nor use solutions to inverse problems to initialize our models but validate 3D dynamic rupture forward simulations with regional observations retrospectively. Our modeled slip distributions (Fig. 3b) are relatively smooth, result in a larger moment of the first earthquake, and have slightly worse geodetic misfit (by 2.2%, Fig. S6) than in Jia et al. (2023). A potential explanation is that our initial conditions do not account for smaller scale heterogeneities, such as local prestress or fault strength asperities, along-strike variations in seismogenic depth, higher variability of pore fluid pressure, and 3D subsurface structure. Both dynamic rupture models, simple and refined, do not match all observed waveform complexities (Figs. 4a,b, 5), which may be due to unaccounted local acceleration and deceleration of the rupture front due to fault segmentation or roughness, frictional or stress heterogeneities, off-fault damage and scattering by structural heterogeneities.

This study models both earthquakes together, which limits the typically vast parameter space of initial conditions for dynamic rupture simulations (Harris, 2004;Ramos et al., 2022) in a self-consistent manner (Taufiqurrahman et al., 2023). For instance, this approach enforces the same regional background stress conditions at the connection of both fault systems and constraints fault strength (relative prestress ratio) to avoid the dynamic triggering of segments that did not slip while reproducing rupture dynamics where slip occurred. Directly constraining remaining uncertainties would require denser observations (e.g., Ben-Zion et al., 2022). However, the models can, based on their self-consistency and the achieved fit to observational data, serve as strong first-order constraints of the doublet's mechanics. They can also serve as a base for refined dynamic rupture simulations to which additional data inferences can be added as they become available. These can be from static, kinematic, or dynamic slip inversion models, for example, using a Bayesian approach as for the 2020 Elaziğ earthquake (Gallovič et al., 2020).

Joint observing and modeling of the long-and short-term, dynamic, and static earthquake interactions are crucial for seismic hazard assessment of active multifault systems to overcome the resulting challenges to empirical hazard assessment and rapid response efforts. For example, joint mapping and seismological analyses may help validate our imperfect keystone hypothesis by determining pre-and postsequence crustal stresses. Combined with dynamic rupture modeling, static and dynamic mechanical inferences could inform rapid assessments of the likelihood of triggering of fault segments within and near a main ruptured fault system by helping estimate the spatial distributions of maximum dynamic stress perturbations from both the earthquakes, the net static stress changes from the rupture sequence, and the relative strengths of different fault segments subject to the local tectonic stress field.

Conclusions

We show that structurally and geodetically informed fault geometries, regional seismotectonics, and early observations typically available within hours to days after large crustal earthquakes (Lacassin et al., 2020) are sufficient constraints to enable meaningful postevent dynamic rupture modeling. Our results have implications for seismic hazard assessment: large multifault earthquakes commonly exhibit seemingly unexpected rupture behaviors, which may be explained by dynamic and static stress transfers between variably oriented fault segments in 3D variable tectonic environments. Unexpected slip patterns and their associated hazards are dynamically plausible, and may be a typical, expected behavior specifically for ruptures of complex immature multifault systems. High stress-drop events may be triggered in geometrically complex fault systems. A single fault system with variable relative strength can host dynamic rupture propagating at highly variable sub-to-supershear rupture speeds across different segments and in different directions, including branching and delayed triggering in backward direction.

Data and Resources

The three-component coseismic displacements recorded by Global Positioning System (GPS) stations of the CORS network (https://geodesy.noaa.gov/CORS/) have been downloaded from the Nevada Geodetic Laboratory website (http://geodesy.unr.edu/). Sentinel-1 Synthetic Aperture Radar (SAR) offsets are processed by the National