On February 6, 2023, a doublet earthquake sequence comprising two major events of moment magnitude M_w 7.8 and M_w 7.6 struck southeastern Türkiye and northern Syria (Melgar et al., 2023; Okuwaki et al., 2023). The sequence occurred within a complex fault system at the triple junction of the Anatolian, Arabian, and African Plates (Figure 1). The first event (M_w 7.8) nucleated at 37.2444°N, 37.0234°E (Melgar et al., 2023) on the Narlı fault, a left-lateral strike-slip branch of the East Anatolian Fault (EAF). Rupture initiated on the Narlı fault, propagated along the branch, and then transferred to the EAF, where it continued bilaterally for ~90 s, producing an ~350-km-long surface rupture (Barbot et al., 2023; Mai PM et al., 2023; Melgar et al., 2023; Okuwaki et al., 2023). Nine hours later, the second event (M_w 7.6) nucleated on the Sürgü–Misis Fault, ~90 km north-northeast of the M_w 7.8 epicenter. This event also ruptured bilaterally, generating an ~150-km-long surface rupture (Jia Z et al., 2023; Xi X et al., 2025). Aftershocks of the doublet are shown in Figure 1. The combined ruptures caused extensive structural damage across 11 major cities in Türkiye and northwestern Syria, resulting in approximately 50,000 fatalities, 100,000 injuries, and direct economic losses exceeding USD 100 billion. The sequence is one of the largest recorded continental strike-slip earthquakes in recent history.

The East Anatolian Fault Zone (EAFZ) is one of the most seismically active regions in Türkiye and the Middle East. Situated at the triple junction of the Anatolian, Arabian, and African Plates, the EAFZ exhibits elevated seismicity driven primarily by interactions between the Arabian and Anatolian Plates (Ambraseys, 1989). Northward indentation of the Arabian Plate induces predominant left-lateral strike-slip displacement along the EAFZ at nearly 10 mm/yr (Reilinger et al., 2006). The EAFZ consists of segmented strike-slip faults with associated bends and branch structures; it has generated multiple historical devastating earthquakes (Taymaz et al., 1991; Ambraseys and Jackson, 1998; Tan O and Taymaz, 2006). Notably, a pronounced seismic gap persists along the southern EAFZ segment, where the last documented rupture before the 2023 M_w 7.8 event was the M_s 7.4 Earthquake in 1513. Long-term strain accumulation within this gap may be the principal driver of the 2023 Kahramanmaraş Earthquake.

Since the 2023 Türkiye Earthquake doublet, numerous studies have been published (Gabriel et al., 2023; Jia Z et al., 2023; Liu CL et al., 2023; Melgar et al., 2023; Rosakis et al., 2023; Ren CM et al., 2024; Li B et al., 2025). Most results indicate a complex rupture process for the M_w 7.8 event, with ongoing debates over the nucleation site on the main fault and rupture velocity on the branch fault. Using InSAR (interferometric synthetic aperture radar) data, seismic records, and aftershock locations to constrain the fault model and joint inversion to resolve the rupture process, Liu CL et al. (2023) proposed that rupture on the main fault nucleated at the triple junction between the main fault and the branch fault and propagated bilaterally; Melgar et al. (2023) also supports this interpretation. Integrating multi-scale seismic and geodetic records with kinematic inversion and dynamic rupture simulations, Jia Z et al. (2023) inferred that rupture on the main fault first propagated toward the northeast, and after a delay extended southwestward. Similar findings were reported by Wang ZJ et al. (2023) and Li B et al. (2025). In contrast, Ren CM et al. (2024), through joint inversion of near- and far-field seismic waveforms, GNSS (Global Navigation Satellite System) data, and the waveform polarization direction, concluded that the main fault nucleated ~9 km west of the triple junction and subsequently ruptured bilaterally. Considering the complex rupture process of the M_w 7.8 Kahramanmaraş event, resolving the nucleation site on the main fault is critical for understanding the triggering mechanism involving the Narlı branch fault.

The interpretation of rupture velocity along the Narlı branch fault also remains inconsistent. Rosakis et al. (2023) analyzed the ratio and waveforms of fault-parallel and fault-normal components from two near-fault stations (4615 and NAR, as shown in Figure 1) and inferred that supershear rupture occurred ~19 km from the epicenter on the Narlı branch fault. In contrast, Melgar et al. (2023), through kinematic multiple-time-window inversion, characterized the overall Narlı branch fault as subshear. Li B et al. (2025) also indicated through dynamic simulations that supershear velocity is not necessary to extend a rupture to the main fault. For the main fault, Jia Z et al. (2023) indicated a rupture velocity of 3.2 km/s along the northeastern segment, consistent with subshear rupture, whereas Ren CM et al. (2024) and Li B et al. (2025) concluded that the northeastern segment exhibited supershear characteristics.

In this study, we integrate relocated aftershock data, surface rupture traces, GNSS observations, and near-field strong ground motion records to construct a three-dimensional (3D) dynamic rupture model of the 2023 M_w 7.8 Kahramanmaraş Earthquake. Finite fault solutions are used to constrain the stress drop on the fault plane. Through 3D dynamic rupture simulations, we compare synthetic results with recorded seismic station data and GPS data to constrain the nucleation location where rupture propagated from the branch fault to the EAF during the 2023 M_w 7.8 Kahramanmaraş Earthquake.

Dynamic rupture simulations provide valuable insight into the triggering mechanism of Narlı branch faults on the main fault. However, key simulation parameters, such as fault geometry, background stress distribution, and friction parameters, are still difficult to constrain directly through in situ measurements.

From the surface fault trace data provided by Reitman et al. (2023) and the aftershock relocation results from Melgar et al. (2023), we construct a 3D nonplanar fault model for the 2023 M_w 7.8 Kahramanmaraş Earthquake event. We retain only the Narlı fault where the mainshock occurred and exclude minor branch faults (Melgar et al., 2023). The M_w 7.8 Kahramanmaraş Earthquake ruptured more than 350 km along three primary faults—the Amanos, Pazarcık, and Erkenek Faults—and the Narlı branch faults. Our fault dataset comprises approximately 150 scattered points distributed along a roughly 380-km-long fault system. These points are projected onto specific depth levels (0, 10, 15, and 20 km) according to their depth ranges, and smooth curves are fitted at each depth to generate four constraining lines. We then apply the spline-constraint method in Gmsh to construct the fault model. The overall fault dip angle is 85°, decreasing to 80° in the northeastern Erkenek segment, consistent with the results of Melgar et al. (2023). The fault width is 20 km.

In our study, we adopt the slip-weakening law as follows (Andrews, 1976):

Here, τ represents the frictional stress at a cumulative slip δ; $ {\tau }_{p}={\mu }_{s}{\sigma }_{n} $ denotes the yield strength, where $ {\mu }_{s} $ is the static friction coefficient and $ {\sigma }_{n} $ is the normal stress; $ {\tau }_{d}={\mu }_{d}{\sigma }_{n} $ is the residual stress, where $ {\mu }_{d} $ the dynamic friction coefficient; and D_c is the critical slip-weakening distance. Khalifa et al. (2018) examined rock strength near the EAF and noticed substantial spatial variations along its segments. Accordingly, we assign different static and dynamic friction coefficients along the EAF. We set $ {\mu }_{s}=0.4 $ and $ {\mu }_{d}=0.21 $ for the Narlı fault segment, and $ {\mu }_{s}=0.412 $ and $ {\mu }_{d}=0.12 $ for the EAF segment. The characteristic slip distance D_c is 0.2 m within the depth range of 0–20 km.

Initial stresses, a critical parameter controlling rupture dynamics, remain challenging to constrain for historical earthquakes. Li B et al. (2025) integrated seismotectonic observations, assumptions of fault-zone fluid overpressure, and the Mohr–Coulomb failure criterion to constrain segment-specific orientations of maximum principal stress in the EAF region. Following Li B et al. (2025), we assign maximum principal stress directions ranging from N170°E and N205°E, divided into four segments according to the fault location (as shown in Figure 2b). Notably, the maximum principal stress orientations on branch faults deviate substantially from those in adjacent EAF segments, which we attribute to localized tectonic rotation.

We define the initial triaxial principal stresses. They vary with depth within the upper 3 km but remain constant below a depth of 3 km to prevent an excessive stress drop at the fault base. The vertical principal stress is given by

where $ \rho $ denotes the rock density and $ \gamma $ is the pore pressure factor, adopted as 0.71 in this study. The stress regime ratio R, defined as $ R=({S}_{v}-{S}_{H})/({S}_{h}-{S}_{H}) $, is also important to constrain the initial stress, where $ {S}_{h} $, $ {S}_{H} $, $ {S}_{v} $ are the minimum, maximum, and vertical principal stress, respectively. Using the focal mechanism inversion results, Yilmaz et al. (2006) proposed an average R value of 0.715 in the north part of the EAF. In our model, R is set to 0.7 for the northeastern segment of the EAF and the Narlı fault, consistent with the dominant transpressional tectonic regime in this region. In contrast, for the southwestern segment of the EAF, which is characterized by a transtensional tectonic environment (Lyberis et al., 1992), R is assigned as 0.3. Another controlling factor, $ k={S}_{H}/{S}_{v} $, is set to 1.5 for the northeastern segment of the EAF and the Narlı fault and to 1.1 for the southwestern segment of the EAF.

We project the triaxial principal stresses onto the fault plane to obtain the initial normal stress $ {\sigma }_{n} $ distribution (Figure 2a). Heterogeneous initial shear stresses τ₀ based on the kinematic slip model are further introduced. Ren CM et al. (2024) provided a reference distribution of fault slip. We calculate the stress drops on the fault by using the slip model and the Okada model (Okada, 1985, 1992). The residual stress $ {\tau }_{d} $ is derived from the normal stress and the dynamic friction coefficient. Superimposing the residual stress and the calculated stress drop yields the initial stress distribution on the fault (Figure 3a). The stress intensity parameter S (defined as $ ({\tau }_{p}-{\tau }_{0})/({\tau }_{0}-{\tau }_{d}) $) is distributed as shown in Figure 2c. The distribution of the S-value indicates that large asperities predominantly occur along geometrically continuous and planar fault segments (Zhang YJ et al., 2023). By comparing simulated ground motion with recorded station data and GPS records through trial and error, we develop three distinct stress models. The stresses illustrated in Figure 3a are designated as Model 1. On the basis of Model 1, we introduce an additional stress perturbation in the localized regions to the southwest of the triple junction of the EAF to define Model 2 (Figure 3b). As shown in Figure 3c, Model 3 shares the same initial stresses as Model 2 except for the Narlı branch fault, where the initial shear stress is higher. The hypocenter location is assigned based on high-precision hypocenter relocation results from Melgar et al. (2023) and is situated ~35 km from the triple junction along a branch fault at a depth of 12 km. The radius of the nucleation zone is 3.5 km. The initial shear stress of the nucleation zone is set 0.1 MPa higher than the yield strength to trigger rupture. Our model incorporates the 3D velocity structure in the EAFZ derived by Wang Z et al. (2024) using double-difference tomography integrated with multi-parameter seismic tomography. The density is calculated following the empirical equation (Brocher, 2005).

We utilize the open-source software drdg3d, developed by Zhang WQ et al. (2023), to simulate the dynamic rupture processes of the 2023 M_w 7.8 Kahramanmaraş Earthquake. The drdg3d software utilizes a discontinuous Galerkin (DG) method with tetrahedral meshes (Hesthaven and Warburton, 2008) and has been widely adopted for dynamic rupture simulations of earthquakes (Ramos et al., 2021; Biemiller et al., 2022; Wang ZJ et al., 2023). The drdg3d method incorporates a hybrid numerical flux combining upwind and central fluxes, thereby reducing mesh dependency and enhancing computational efficiency. The method has demonstrated excellent performance across multiple benchmark models in the Statewide California Earthquake Center (SCEC)/ U.S. Geological Survey (USGS) Spontaneous Rupture Code Verification Project (https://strike.scec.org/cvws/; Harris et al., 2009). To raise computational efficiency, our model utilizes a grid size of 1 km with third-order spatial accuracy. The time step is 0.07 s, and the total simulation time is 85 s.

We construct a 3D fault geometry for the 2023 M_w 7.8 Kahramanmaraş Earthquake based on surface fault trace data (Reitman et al., 2023) and relocated aftershocks (Melgar et al., 2023). We then develop heterogeneous initial stresses along the fault, constrained by the kinematic slip inversion (Ren CM et al., 2024), near-field seismogram, and GPS observations. We present three different dynamic rupture models to investigate the nucleation behavior near the triple junction of the main EAF.

Model 1 reproduces the inverted slip distribution, observed near-field seismic waveforms, and GPS observations. Following nucleation on the Narlı branch fault, where supershear rupture occurs, rupture nucleates near the triple junction and propagates northeastward first. Nearly 4 s later, the rupture also propagates southwest toward the Amanos Fault. To facilitate nucleation at 9 km southwest of the triple junction of the main EAF, as hypothesized by Ren CM et al. (2024), an artificially elevated initial shear stress is required in that region, as shown in Model 2. However, this nucleation difference itself does not greatly influence the overall characteristics of the final slip distribution, near-field ground motion, and surface displacement, demonstrating that the complex 3D geometry and rupture process mask early nucleation details. Thus, both nucleation locations, the triple junction and 9 km southwest of the triple junction, are possible for the 2023 M_w 7.8 Kahramanmaraş event. However, the triple junction is more likely to nucleate because of the positive dynamic Coulomb stress perturbation there, whereas the condition for the southwest nucleation is stricter—a higher initial shear stress perturbation is necessary. Moreover, we elevate the initial shear stress along the Narlı branch fault in Model 3, which increases the rupture velocity before the rupture transfers to the main EAF, and the P-wave approaching supershear rupture velocity results in a narrower Mach cone, making the southwest nucleation on the main EAF more difficult, although the comparison of source time function of the three models and the USGS result indicates such a P-wave approaching supershear speed is less likely for the 2023 M_w 7.8 event. Our models reconcile kinematic interpretation of the 2023 M_w 7.8 Kahramanmaraş Earthquake and highlight the critical role of branch-fault rupture velocity and initial stress perturbation in controlling rupture transfer across complex fault junctions.