Rupture segmentation on the East Anatolian fault (Turkey) controlled by along-strike variations in long-term slip rates in a structurally complex fault system

The East Anatolian fault (EAF) system (Turkey) forms a major plate boundary between the Arabian and Anatolian plates (Fig. 1), with long-term slip rates varying from 1 mm/yr to 10 mm/yr among seven segments on the main strand (Fig. 1) that are defined by structural complexities and inferred historical earthquake rupture extents (Duman and Emre, 2013). The EAF is a relatively immature fault exhibiting major structural complexity, including multiple step-overs, releasing and restraining bends, branch faults, and parallel structures like the Bitlis-Zagros fold-and-thrust belt (Duman and Emre, 2013). These structures accommodate transpression along the Erkenek, Pütürge, and Palu segments that parallel the Bitlis-Zagros suture and transtension along the Amanos segment that bounds the Karasu trough. As a result of slip trade-off within a complex fault network, the long-term slip rates along the main strand of the EAF may vary along strike. Inferred fault slip rates are usually associated with large uncertainties (Duman and Emre, 2013; Table S1 in the Supplemental Material¹), especially along the southern segments. For example, the slip rate on the Pazarcık segment is inferred to be 6.5–7.0 mm/yr (Reilinger et al., 2006) to ~10 mm/yr (Aktug et al., 2016) from geodetic studies depending on modeling assumptions (Text S6; Figs. S7 and S8), 9 mm/yr from paleoseismological data (Karabacak et al., 2011), and 4.0–4.6 mm/yr from geomorphology (Westaway, 2003).

Despite the fault structural immaturity, many large earthquakes (magnitude >7) have occurred on the EAF during the historical and instrumental periods. These M 6.7–7.5 earthquakes, with the exception of the 1114 CE M > 7.8 earthquake (Ambraseys and Jackson, 1998), occupy different parts of the EAF main strand, showing a clear pattern of rupture segmentation (Fig. 1). These ruptures are thought to terminate at structural complexities along the EAF (Duman and Emre, 2013), although the precise rupture extents of these historical earthquakes, or even the faults they ruptured, are still debated (Duman and Emre, 2013; Güvercin et al., 2022).

However, on 6 February 2023, a devastating Mw 7.8 earthquake hit the region near Kahramanmaraş, unzipping multiple segments over 300 km. Nine hours later, another Mw 7.6 earthquake ruptured most of the Çardak-Sürgü fault. The official death toll reached 60,000, making it one of the deadliest earthquakes in history. The total estimated economic loss exceeds 100 billion U.S. dollars (Dal Zilio and Ampuero, 2023). The cause of such a large earthquake crossing multiple structural barriers is still unknown. Rupture segmentation is likely linked to fault geometry (Wesnousky, 2006) and variations in frictional properties (Kaneko et al., 2010), normal stress (Molina-Ormazabal et al., 2023), or long-term slip rate (Perez-Silva et al., 2022). Understanding rupture segmentation and seismic supercycles on continental transforms is vital to seismic hazard evaluation.

We investigated rupture segmentation by modeling seismic cycles on the EAF and the Çardak fault with three-dimensional (3-D) fault geometry. With increasing variation of long-term slip rate among segments of the EAF, earthquake cycles transition from systematic multi-segment ruptures to supercycles with single- and multi-segment ruptures (Fig. 2A; Fig. S2; Table S3). Some of these cycles mimic the historical rupture patterns of the EAF, with mostly single-segment ruptures and occasional multi-segment earthquakes, and some individual events reproduce the crustal deformation of the February 2023 Kahramanmaraş doublet. Our physics-based seismic cycle simulations highlight the importance of heterogeneous long-term slip rates and shed light on the long-term variability of earthquake size that shapes seismic hazards.

We conducted quasi-dynamic seismic cycle simulations using the boundary integral method, employing the cross-validated code Unicycle (Barbot, 2018, 2019; Jiang et al., 2022). Due to the technical difficulty of simulating long earthquake cycles on disconnected fault segments (Ozawa et al., 2023), we focused on exploring the role of heterogeneous long-term slip rate (Fig. 2A). We use a smooth fault geometry of the EAF and the Çardak fault without abrupt bends or stepovers constrained by surface fault traces and seismicity. Both faults are modeled down to 20 km depth. The EAF is modeled from Antakya to the south of Lake Hazar, covering the southernmost four segments: the Amanos, Pazarcık, Erkenek, and Pütürge segments. The Amanos and Pazarcık segments are divided into two subsegments with different fault strike. The faults are discretized into 1-km-square vertical subfaults with smoothly varying strike, leading to a total of 10,580 subfaults (Fig. 2B). Each subfault is assigned a set of parameters, including normal stress, friction parameters, and long-term slip rates via back slip (Table S2). The long-term loading accumulates only shear stress on the faults. We assume the normal stress loading is accommodated by off-fault structures. The models resolve the cohesive zone width and the nucleation size (Text S2) by a factor of at least 4–5. We first explore the influence of variability in long-term slip rate within the range of reported geodetic and geologic slip rates on the isolated EAF. The long-term slip rates come from a compilation of various geological and geodetic studies in the literature (Duman and Emre, 2013; Styron and Pagani, 2020; Table S1). Several simulations readily align with the historical earthquake records (Fig. S3; Tables S3 and S4). Starting from models that produce both single- and multi-segment ruptures, we then add the Çardak fault to the model. Due to the low long-term slip rate and the rare occurrence of earthquakes on the Çardak fault, it does not disrupt the rupture segmentation pattern on the EAF (Fig. S10). We vary the slip rate and normal stress to reproduce both the historical earthquake record and the seismo-geodetic data of the 2023 doublets (Barbot et al., 2023; Jia et al., 2023; Mai et al., 2023; Ren et al., 2024; Xu et al., 2023). For brevity herein, we mainly discuss one of the models that best reproduces the pattern of paleoearthquakes on the EAF and the geodetic observations of the 2023 Kahramanmaraş and Elbistan earthquakes. Some alternative models fit the constraints equally well (Fig. S3; Tables S3 and S4).

The preferred model produces supercycles of single- and multi-segment ruptures (Fig. 2) with magnitudes of partial ruptures approximating those inferred from the historical earthquake record on the EAF (Figs. 2B and 2D). The recurrence time of the simulated earthquakes at the center of each fault segment varies from ~200 yr to ~1300 yr (Fig. 2D; Fig. S5), consistent with constraints from historical records (Fig. 1B) (Duman and Emre, 2013), the paleoseismic record at Lake Hazar (Hubert-Ferrari et al., 2020), and those inferred from seismic catalogs and geodetic strain rate (Güvercin et al., 2022). A Mw 6.9 earthquake comparable to the 2020 Mw 6.8 Elazig earthquake that ruptured part of the Pütürge segment is also reproduced in the model (green areas in Figs. 2B and 2D). We do not attempt to reproduce every detail of historical earthquakes (e.g., exact timing of occurrence, hypocenter locations) due to the nonlinear nature of the modeled system (e.g., Gauriau et al., 2023) and aseismic creep near segment boundaries as a result of changing slip rate or fault orientation, but focus on the overall patterns. The simulated ruptures show consistent endpoints that coincide with transitions between different slip rates and changes in fault orientation. Alternative models with intermediate variations in along-strike slip rates also produce both single- and multi-segment ruptures with roughly consistent endpoints (Fig. 2A; Fig. S3).

Large earthquakes that rupture multiple segments are occasionally observed with recurrence times of thousands of years (Fig. 2D), consistent with estimations from the moment budget of earthquake cycles (Xu et al., 2023). Although some ruptures cross only one barrier, others overcome multiple barriers and result in several-hundred-kilometer-long ruptures, as occurred during the Mw 7.8 2023 Kahramanmaraş and the 1114 CE large earthquakes.

We select events from the simulations that individually resemble the 2023 doublet (Fig. 2D; Fig. S3), but they are separated by several hundreds of years in time (Table S5) instead of 9 hours as observed in nature. The slip distribution of the simulated earthquakes is in overall agreement with that inferred from kinematic inversions of the 2023 doublet using geodetic data (Barbot et al., 2023), although some differences persist (Fig. 3A). For both the observed 2023 and simulated EAF multi-segment rupture, large slip is observed on the Amanos, Pazarcık, and Erkenek segments, with peak slip reaching ~8 m, concentrating on the Pazarcık segment. In both observation and simulation, a larger peak slip of ~10 m is produced on the Çardak fault. The larger slip and smaller rupture area of the Mw 7.6 Elbistan earthquake is consistent with a higher stress drop on the Çardak fault. Some notable differences between the simulated and observed 2023 earthquakes include a slight mismatch in the extent of the mainshock rupture near the southern and northern tips, a larger slip amplitude on the northeast portion of the Çardak-Sürgü fault, and a more homogeneous simulated slip distribution for both earthquakes. Our slip distribution is also consistent with 3-D dynamic rupture models (e.g., Wang et al., 2023), although the quasi-dynamic model prevents discussions on supershear ruptures observed on Narlı fault, Pazarcık segment, and Çardak fault (Abdelmeguid et al, 2023; Ren et al., 2024).

We compare the modeled surface displacement field against the geodetic observation from high-rate GNSS stations and optical and synthetic aperture radar data. The model successfully reproduces the first-order patterns of the surface displacement field due to the 2023 doublet (Figs. 3B–3E). The variance reduction (Text S3) for GNSS observations for the Mw 7.8 Kahramanmaraş earthquake and Mw 7.6 Elbistan earthquake reaches 74% and 97%, respectively. The model also explains the synthetic aperture radar data with a variance reduction ranging from 70% to 80% in the range direction and 32%–55% in the azimuthal direction for the Sentinel-1 pixel-tracking data, 80%–83% for the Advanced Land Observing Satellite-2 (ALOS-2) interferograms, and 26% for Sentinel-2, which is likely affected by snow coverage (Barbot et al., 2023) (Table S5). Notably, simulated event pairs with rupture extents like the 2023 doublet within alternative models also fit the geodetic observations well, with variance reduction comparable to the model shown here in the main text (Table S4).

Heterogeneous long-term slip rate distribution along a fault is common in nature (e.g., Weldon et al., 2004). On strike-slip faults, a non-uniform slip rate distribution can originate from structural complexities such as a change in fault orientation (e.g., Li et al., 2023) or the presence of parallel or branching faults (e.g., Gauriau and Dolan, 2021). This potentially results in neighboring structural segments with distinctively different slip rates. In contrast, subduction zones often exhibit smoothly varying slip rates along strike (DeMets et al., 1990), due to varying proximity to the Euler pole that characterizes relative plate motion.

Our simulation highlights the significant role of long-term slip rates on rupture dynamics. Variations in long-term slip rates along a fault desynchronize the loading of neighboring segments, placing them in varying states of readiness to rupture at any given time. Spatially varying levels of initial stress favor complex rupture dynamics, with earthquakes of varying sizes and irregular recurrence patterns. Better constraints on spatially and temporally varying long-term slip rate on natural faults may therefore improve physics-based estimates of seismic hazard.

Fault geometry is well known to correlate with rupture segmentation for strike-slip faults. Earthquake rupture endpoints often coincide with fault stepovers or bends (King and Nábělek, 1985). When a fault stepover surpasses 5 km in width, it acts as a strong barrier that inhibits rupture propagation (Biasi and Wesnousky, 2016). Numerical studies also find fault bends are favorable sites for rupture nucleation and termination (Duan and Oglesby, 2005; Ozawa et al., 2023; Sathiakumar and Barbot, 2021).

In our simulations, a smooth non-planar fault geometry alone is insufficient to account for the observed rupture segmentation (Fig. 2C). Instead, we observe that the along-strike variation in long-term slip rates is required to reproduce the observed rupture segmentation (Fig. 2D; Fig. S10). However, the influence of geometrical complexities is somewhat subdued in our model due to the use of smooth fault bends and the absence of stepovers (Delogkos et al., 2023; Howarth et al., 2021). Fault step-overs may occasionally prevent rupture propagation, setting the conditions for sequences of single- and multi-segment ruptures. However, the role of geometrical complexities as potential barriers to ruptures may be limited if the step-overs only persist within shallow layers. The effects of structural complexity and along-strike variations in long-term slip rate may be strongly linked, as variations in long-term slip rates are associated with fault geometry and the complexity of the fault network.

The along-strike variation of long-term slip rates, which arises from kinematic and structural complexity within a fault network, plays a crucial role in shaping the dynamics of earthquakes, giving rise to cycles of single- and multiple-segment ruptures. The interplay of fault geometry with heterogeneous loading within an inter-connected fault network shapes a complex recurrence pattern with earthquakes of diverse sizes (Fig. 4). Physics-based simulations of the seismic cycles may reproduce some first-order characteristics including the rupture extent, magnitude, recurrence times, and surface displacements of earthquakes. Future models with additional physics such as wave-mediated stress, migration of fluids, and enhanced dynamic weakening may reconcile the yet unexplained delayed triggering of the Çardak event. Integration of realistic constitutive behavior and structural setting is key to reproducing the range of earthquake sizes and to better estimating the resulting seismic hazards. More near-fault geodetic and geologic constraints on the heterogeneous loading and high-quality paleoseismic trench sites that provide robust recurrence records may further advance understanding of earthquake cycles.

We thank Sezim Ezgi Güvercin, Judith Gauriau, and Baoning Wu for helpful suggestions on the manuscript. We thank three anonymous reviewers for their constructive comments. We thank editor Tracy Rushmer and the Geology staff for their support throughout the publication process.