# Rupture Dynamics and Near-Fault Ground Motion of the Mw7.8 Kahramanmaraş, Turkey earthquake of February 6, 2023

Hideo Aochi, Victor Manuel Cruz-Atienza

► To cite this version:

Hideo Aochi, Victor Manuel Cruz-Atienza. Rupture Dynamics and Near-Fault Ground Motion of the Mw7.8 Kahramanmaraş, Turkey earthquake of February 6, 2023. Seismica, 2025, 4 (1), 10.26443/seismica.v4i1.1432. hal-05454679

HAL Id: hal-05454679

https://brgm.hal.science/hal-05454679v1

Submitted on 12 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Rupture Dynamics and Near-Fault Ground Motion of the Mw7.8 Kahramanmaraş, Turkey earthquake of February 6, 2023

H. Aochi *1,2, V. M. Cruz-Atienza *3

¹Bureau de Recherches Géologiques et Minières, Orléans, France, ²Laboratoire de Géologie de l'Ecole Normale Supérieure, CNRS UMR 8538, PSL Research University, Paris, France, ³Departamento de Sismología, Instituto de Geofísica, Universidad Nacional Autónoma de México, Ciudad de México, Mexico

**Author contributions:** Conceptualization: H. Aochi, V. M. Cruz-Atienza. Data Curation: H. Aochi, V. M. Cruz-Atienza. Formal Analysis: H. Aochi, V. M. Cruz-Atienza. Funding Acquisition: H. Aochi, V. M. Cruz-Atienza. Investigation: H. Aochi, V. M. Cruz-Atienza. Methodology: H. Aochi, V. M. Cruz-Atienza. Project Administration: H. Aochi. Software: H. Aochi, V. M. Cruz-Atienza. Validation: H. Aochi, V. M. Cruz-Atienza. Visualization: H. Aochi, V. M. Cruz-Atienza. Writing – original draft: H. Aochi, V. M. Cruz-Atienza. Writing – review & editing: H. Aochi, V. M. Cruz-Atienza.

**Abstract** We studied the dynamic rupture propagation of the February 6th, 2023 (Mw 7.8, 01:17 UTC) Pazarcık (Kahramanmaraş), Turkey, earthquake by incorporating the non-planar fault structure, the regional stress field, and a data-driven friction parameterization into numerical simulations. To explain the rupture extent of 200 km and the average speed, a regional non-uniform load is necessary and was determined from the orientation and intensity of the principal stresses. Careful analysis of near-fault strong motions suggests that the critical slip-weakening distance ($D_c$) varies smoothly along the fault strike (between 0.6 - 1.2 m) with mean value of 0.86 ± 0.34 m. Such friction and prestress heterogeneities help to explain local kinematic features of the rupture process imaged by Delouis et al. (2023) (e.g., two supershear rupture transients) where the fault geometry played a major role. As expected, we found clear correlation between rupture speed and radiation efficiency ($η_r$) along the fault, both metrics with peak values near the maximum PGAs recorded. This is the first earthquake where local heterogeneity of rupture dynamics and near-fault ground motion can be studied together and the methodologies introduced will serve to generate comprehensive earthquake scenarios to assess the seismic hazard in other regions.

**要旨** 本研究では、2023年2月6日(Mw 7.8、世界標準時01:17)にトルコで発生した Pazarcık (Kahramanmaraş) 地震の動的破壊伝播について、非平面断層構造、地域応力場、データから推定されたすべり弱化摩擦法則を数値シミュレーションに統合することによって検討した。断層近傍の強震動データから慎重に解析した結果、すべり弱化の臨界距離($D_c$)は約0.6-1.2m(平均0.86m)で、走向に沿って空間的に変動することが示された。その結果、断層の形状が破壊の進行に主要な役割を果たしていること、また、破壊の広がりと破壊速度を説明するためには、地域的に不均一な応力が必要であることが示された。発表された運動学的インバージョンモデルと整合させるためには、長さ250kmにわたって最大7つの応力ゾーンが必要である。走向に沿った$D_c$変動は、その下の地震発生帯からの放射エネルギーに従って、ピーク速度分布を改善する。これは、動的破壊モデルと観測データを用いて、地震破壊と断層近傍の地震動の空間的不均質性を研究した最初の地震である。展開されたパラメタリゼーションアプローチは、包括的な地震シナリオを再現するのに有効であり、今後の地震ハザードの応用に役立つと考えられる。

Production Editor: Yen Joe Tan
Handling Editor: Alice Gabriel
Copy & Layout Editor: Kirsty Bayliss

Received:
August 21, 2024
Accepted:
April 22, 2025
Published:
June 18, 2025

<footnote>*Corresponding author: aochi.hideo@gmail.com</footnote>
**Resumen** Estudiamos la dinámica de la ruptura del terremoto del 6 de febrero de 2023 (Mw 7.8, 01:17 UTC) en Pazarcık (Kahramanmaraş), Turquía, incorporando a las simulaciones computacionales la geometría no plana de la falla, el campo regional de esfuerzos y una parametrización de la fricción desprendida de los datos. Para explicar la extensión de la ruptura de 200 km y su velocidad media, es necesaria una carga regional no uniforme, que se determinó a partir de la orientación y la intensidad de los esfuerzos principales. El análisis de los movimientos fuertes cercanos a la falla sugiere que el deslizamiento crítico de fricción ($D_c$) varía suavemente a lo largo del rumbo de la falla (entre 0,6 - 1,2 m) con un valor medio de 0,86 ± 0,34 m. Tales heterogeneidades de fricción y esfuerzos permitieron explicar características cinemáticas locales de la ruptura observadas por Delouis et al. (2023) (e.g., dos transitorios supershear) donde la geometría de la falla jugó un papel importante. Como predice la teoría, encontramos una clara correlación entre la velocidad de ruptura y la eficiencia radiativa ($\eta_r$) a lo largo de la falla, ambas métricas con valores pico cerca de los PGA máximos registrados. Este es el primer terremoto en el que se pueden estudiar conjuntamente la heterogeneidad de la dinámica de la ruptura y el movimiento del terreno próximo a la falla, de modo que las metodologías introducidas servirán para generar escenarios sísmicos que permitan evaluar con mayor confiabilidad la peligrosidad sísmica en otras regiones.

**Non-technical summary** Two strong earthquakes hit the south-east of Turkey and Syria on the February 6th, 2023. More than ten stations were installed along the Eastern Anatolian fault and recorded the strong ground motion near the fault. We study Kahramanmaraş earthquake (Mw 7.8, 01:17 UTC) by integrating the non-planar fault structure, the regional stress field, and a data-driven spatially variable slip-weakening friction into numerical simulations. A careful analysis from the near-fault strong motion data suggests a critical slip-weakening distance ($D_c$) of about 0.6-1.2 m (0.86 m on average) showing a spatial variation along the strike. The fault geometry plays a principal role in rupture progress, and non-uniform regional stress is required to explain the rupture extension and rupture speed. In order to be consistent with a published kinematic inversion model, up to seven-zone stress model is needed over 250 km long. $D_c$-variation along the strike improves the peak velocity distribution, according to the radiation energy from the seismogenic zone beneath. This is the first earthquake allowing to study the spatial heterogeneity of seismic rupture and near-fault ground motion using simulations and observations. The deployed parameterization approach is valid for reproducing comprehensive earthquake scenarios, so would be useful for further seismic hazard applications.

# 1 Introduction

On the 6th Februrary 2023, two strong earthquakes hit Eastern Turkey, an Mw 7.8 at 01:17:32 Universal Time (UTC) in Pazarcık (Kahramanmaraş) and then an Mw 7.7 at 10:24:47 UTC in Elbistan (Kahramanmaraş). Seismological information has been shared since then by the Turkish organizations AFAD (Disaster and Emergency Management Authority) and KOERI (Kandilli Observatory and Earthquake Research Institute, Boğaziçi University) in particular. The earthquakes occurred in a seismic gap previously identified for its low strain rates, i.e. for a long recurrence time of historical earthquakes (e.g. Güvercin et al., 2022; Karabulut et al., 2023). As numerous seismological/geodetic/geological studies have already shown (e.g. Melgar et al., 2023; Jia et al., 2023; Barbot et al., 2023; Delouis et al., 2023; Chen et al., 2024), these large earthquakes are related to multiple fault segments with major surface ruptures along the East Anatolian fault zone. In particular, the first event (hereafter the Kahramanmaraş earthquake) started on the Narlı normal fault before reaching the Kahramanmaraş Triple Junction where rupture propagated bilaterally with a left-lateral strike-slip mechanism for about 300 km along the main section of the East Anatolian fault (EAF).

To better understand the main rupture of the Kahramanmaraş shock in a regional context, let us examine some aspects of the 1999 Izmit earthquake, which occurred on the North Anatolian fault (i.e., 600 km northwest; Figure 1) and has been extensively studied through seismological and geodetic data, satellite interferometry and field observations. Although few near-

fault stations recorded the event, several dynamic rupture simulations were carried out to discuss the rupture transfer from one segment to another (e.g. Harris et al., 2002; Aochi and Madariaga, 2003) along the almost continuous fault trace that contained, however, some irregularities such as bends and jogs. Aochi and Madariaga (2003) tested different fault geometries and demonstrated that the dynamic rupture process of that earthquake was strongly controlled by small variations in the fault geometry. The fault structure inferred from the analysis of satellite interferograms allowed for improved earthquake models in terms of the rupture front acceleration and the resulting final slip distribution. Among the four near-field seismic stations operational during the Izmit event, the two closest within a few kilometers from the fault (SAR, YPT) recorded relatively simple velocity waveforms associated with the passage of the rupture front next to the stations. Theoretically, at such distances from the fault, the velocity waveform is close to the slip-rate at the nearby rupture front so that it was possible to quantify the fault friction in this case (Cruz-Atienza and Olsen, 2010). Dynamic rupture simulations were able to reproduce such waveforms by assuming a mechanically reasonable stress reduction (slip-weakening) process within an appropriate scale. Near-fault observations remain limited to a small number of earthquakes and observational sites, similar to station Pump Station 10 (PS10) during the 2002 Denali earthquake (Dunham and Archuleta, 2004; Eberhart-Phillips et al., 2003), where friction could also be quantified (Cruz-Atienza and Olsen, 2010). From this perspective, the Mw 7.8 Kahramanmaraş earthquake, which has
a similar faulting mechanism and was recorded by at least 11 near-fault accelerometers (i.e. within 3 km of the source), represents a globally unique opportunity to apply the similar analysis technique from the 1999 Izmit earthquake to study the dynamics of the rupture process and its impact on strong motions at the local scale, taking into account the non-planar fault geometry.

In the past, fault geometry and earthquake rupture were first examined from a geological point of view. Geometrical irregularity and fault segmentation have been shown relevant to the initiation, development and termination of the rupture process (e.g. King and Nabhêlek, 1985; Nakata et al., 1998). Dynamic rupture simulations on segmented planar faults were possible in the 1990's (Harris and Day, 1993; Kame and Yamashita, 1997; Kase and Kuge, 1998) until complex fault geometries became accessible with different methods in the 2000s (Aochi et al., 2000; Oglesby et al., 2000; Aochi and Fukuyama, 2002; Kame et al., 2003; Ando et al., 2004; Harris et al., 2009; Cruz-Atienza and Virieux, 2004; Cruz-Atienza et al., 2007). Nowadays, dynamic rupture simulations are systematically developed for many earthquakes to understand their generation process in geodynamic frames (e.g. Kaneko et al., 2010), as is also the case for seismic radiation to better estimate the seismic hazard (e.g., Guatteri et al., 2003; Olsen et al., 2009; Gallovič and Valentová, 2023). Since large earthquakes tend to occur repeatedly on known and increasingly well-characterized faults, the fault geometry is a preset condition where the governing friction law and the initial stress field represent the major challenge to achieve a better understanding of the phenomenon. For this reason, it is essential to have physically consistent methodologies to establish the prestress conditions and to extract as much information as possible about friction from the recorded seismograms, which is what is proposed in the present work.

Several studies on the dynamic rupture of the 2023 Kahramanmaraş earthquake have been conducted in two and three dimensions to explain the multiple segmentation of the rupture and emphasize the importance of the system heterogeneity (Jia et al., 2023; Ding et al., 2023; Gabriel et al., 2023; Wang et al., 2023; Abdelmeguid et al., 2023). These works focused on the mechanisms that allowed the rupture transfer from the initial splay fault to the EAF and then propagate bilaterally along the nonplanar fault that characterized the event. They also sought to explain why the Mw 7.7 Elbistan earthquake occurred nine hours later and only ~20 km to the north. In this paper, we focus on the 200 km long southwestern fault segment (Figure 1) of the EAF that ruptured in the Mw 7.8 Kahramanmaraş earthquake, because this is the first event where local heterogeneity of rupture dynamics and near-fault ground motion can be studied together from both the simulations and the near-fault seismograms, which are invaluable observations affected predominantly by the rupture process near the seismic stations. Our primary objective here is the dynamic explanation of the rupture process, described kinematically in an extraordinary way previously, and of the numerous and unprecedented near-fault strong motion records.

# 2 Earthquake Dynamic Model

## 2.1 Fault geometry

It has long been recognized that fault geometry is certainly one of the most important factors in earthquake dynamics (e.g., Aochi and Fukuyama, 2002; Aochi and Madariaga, 2003; Tago et al., 2012; Cruz-Atienza and Virieux, 2004; Cruz-Atienza et al., 2007; Adda-Bedia and Madariaga, 2008). For this reason, we built a detailed fault model based on the Line-of-Sight displacement discontinuity clearly defined in satellite interferograms (e.g. Reitman et al., 2023), where significant along-strike geometric variations are found (Figure 1). Evidence of surface rupture extends across the entire region, with offsets of up to 7.5 m in some places (e.g. Provost et al., 2024). As for the model at depth, we assumed a simple vertical fault up to 17 km depth, which is consistent with the left-lateral strike-slip focal mechanism (e.g. AFAD, Global CMT among others).

We are primarily interested in the relationship between rupture propagation and near-fault ground motions along the fault segment of the EAF shown in Figure 1b, namely the southwestern part of the Mw 7.8 rupture. Therefore, although the earthquake initiated on a secondary splay fault before reaching the EAF where rupture propagated bilaterally (e.g. Melgar et al., 2023; Barbot et al., 2023; Delouis et al., 2023), the non-planar fault model we adopted represents the main continuous segment of EAF over 250 km long without branches. Thus, in our numerical simulations, rupture nucleation is assumed around the triple junction where the splay fault meets the EAF (Figure 1b). This assumption does not undermine the generality of the model and allows us to focus the discussion on the rupture process in the target area only. The local reference frame we use is rotated 30° clockwise, so that Cartesian coordinates X (N30° E) and Y (N60° W), assumed in the analysis, roughly correspond to the fault-parallel and fault-normal directions, respectively, particularly where most of the stations of interest are located.

## 2.2 Friction Law

We assume that fault slip is governed by a linear slip-weakening law (e.g. Ida, 1972). The fault strength ($\sigma$) is thus a function of fault slip ($\Delta u$) so that

$$ \sigma(\Delta u) = \tau_r + (\tau_p - \tau_r) \left(1 - \frac{\Delta u}{D_c}\right) H\left(1 - \frac{\Delta u}{D_c}\right) \quad (1) $$

for ($\Delta u \ge 0$), where $\tau_p$ and $\tau_r$ are the peak strength and residual stresses, $D_c$ is the critical slip-weakening distance, and $H(\cdot)$ is the Heaviside step function. The breakdown strength drop is defined as $\Delta\tau_b = \tau_p - \tau_r$ and according to the Coulomb failure criterion,

$$ \tau_p = c + \mu_s \sigma_n \quad \text{and} \quad \tau_r = \mu_d \sigma_n \quad (2) $$

where $\sigma_n$ is the normal fault stress, $\mu_s$ and $\mu_d$ are static and dynamic friction coefficients, and c is the fault cohesion. The model parameters are summarized in Table 1, which are the same as those previously used by
<figure><img src="image_2.png" /></figure>

**Figure 1** Study area of the 2023 Turkish earthquake sequence. (a) Map of faults, stations and seismicity during the first 72 hours after the 01:17 Kahramanmaraş earthquake. The epicenters of the two principal events are illustrated by a star. In the upper left-hand corner, the map of Turkey is shown. (b) The detailed map of the fault model adopted for numerical simulations is shown in local coordinate (X, Y) rotated to N30°E. The map area corresponds to a rectangle in panel (a). Two open stars indicate the nucleation points selected for dynamic rupture simulations. (c) The fault model for the 1999 Izmit earthquake, for comparison, after Aochi and Madariaga (2003). The areas of panels (b) and (c) are also illustrated in a regional map on the top left of panel (a).

Aochi and Ulrich (2015). The constitutive parameters in Equation 2 are constant, but $τ_p$ and $τ_r$ are expected to vary according to $σ_n$ along both dip and strike. In addition, $D_c$ can also vary in space. We shall explain this along with the pre-stress condition in the next section.

## 2.3 Pre-Stress Condition

Although estimating the stress field prior to an earthquake is always difficult, Aochi and Madariaga (2003) and Aochi and Ulrich (2015) proposed a simulation framework where the initial and boundary conditions on the fault are consistent with generic and site-specific knowledge. In this framework, it is assumed that the optimal orientation of the fault is tangential to the great circle described by the relative motion of tectonic plates. In the region of the East Anatolian fault, the motion between the Anatolian and Arabian plates is less than half that of the North Anatolian fault region (Reilinger et al., 2006), where major earthquakes occurred over the past century, such as the Mw 7.6 Izmit earthquake in 1999 (Figure 1c). Although the horizontal velocity field in the EAF region is difficult to quantify due to its low strain rates, Aktuğ and Kiliçoğlu (2005) and Mahmoud et al. (2013) independently estimated the Euler Pole parameters associated with the relative plate motion.

Figure 2 summarizes the strike of our fault model

as well as the great circle tangential directions derived from the two Euler Pole models mentioned above. The strike of the fault varies from N70°E in the north to N20°E in the south. Moreover, since the Euler pole determined by Mahmoud et al. (2013) (49.098°N, 6.043°E) is much further away than the pole determined by Aktuğ and Kiliçoğlu (2005) (33.814°N, 38.417°E), the optimal orientation of the fault in the first case remains nearly the same at latitudes encompassing the fault (red lines), while in the second case, the optimal orientation varies considerably (blue lines) so that both models are inconsistent and thus mutually exclusive. For this reason, as we shall describe in Section 4, we decided to undertake a parametric stress analysis to find reasonable initial conditions for our earthquake model based on the following considerations.

From the strategy proposed by Aochi and Ulrich (2015), we assumed that the shear and normal stresses ($τ$, $σ_n$) on the fault plane (Equation 2) are given by the principal stresses according to the Mohr circle, as schematically illustrated in Figure 3. Considering that the Mw 7.8 earthquake occurred along a strike-slip fault, we let the axes of the maximum and minimum principal stresses ($σ_1$, $σ_3$) be in the horizontal plane and the intermediate stress axis ($σ_2$) in the vertical direction. In nature, these stresses are determined by factors at different scales such as long-term regional deformations and residual strain from local seismicity. However, since
<table><caption>**Table 1** Model parameters used in this study.</caption><thead><tr><td>Parameter</td><td>Quantity (Unit)</td></tr></thead><tbody><tr><td>Static friction coefficient μ<sub>s</sub></td><td>0.3</td></tr><tr><td>Dynamic friction coefficient μ<sub>d</sub></td><td>0.24</td></tr><tr><td>Cohesive force c</td><td>5 MPa</td></tr><tr><td>P- and S-wave velocities V<sub>P</sub> and V<sub>S</sub></td><td>6000 m/s, 3464 m/s</td></tr><tr><td>Shear Modulus G</td><td>32.4 GPa</td></tr><tr><td>Element size Δs in BIEM</td><td>500 m</td></tr><tr><td>Time step Δt in BIEM</td><td>0.0417 s</td></tr><tr><td>Grid size Δs in FDM</td><td>200 m</td></tr><tr><td>Time step Δt in FDM</td><td>0.01 s</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 2** The fault model (black) and the optimal fault plane inferred from the two different Euler pole models. We adopt the pole location at (49.098°N, 6.043°E) from Mahmoud et al. (2013) and (33.814°N, 38.417°E) from Aktuğ and Kiliçoğlu (2005). Two stars indicate the nucleation points supposed in the simulations. The triangles show the seismic station locations. On the right panel, the change in azimuth is compared along latitude.</figcaption></figure>

this study focuses on the coseismic earthquake process and the resulting ground motions, we made the simple assumption that normal tractions increase linearly with depth (Figure 3b) and that shear tractions along the fault are bound by the static and dynamic friction coefficients through the Coulomb failure criterion (straight lines in Figure 3a). For rupture to propagate spontaneously, the potential stress drop, $\Delta\tau = \tau - \tau_r$, should be positive and large enough (Das and Aki, 1977). Given the principal stresses, the optimal orientation of the fault plane is defined as the closest to the Coulomb failure. The angle for this optimal orientation, $\Phi$, is usually measured from the direction of the maximum principal stress in the mechanical framework. In this study, $\Phi$ corresponds to its azimuth in the geographical coordinate system. Thus, for such an optimally oriented fault plane, we define the parameter $T$ (Aochi and Ulrich, 2015) with respect to the Coulomb friction lines such that

$$ 
\begin{aligned}
T &\equiv \frac{\Delta\tau}{\Delta\tau_b \text{ on optimal fault } \Phi} \\
&= \frac{\tau - \mu_d\sigma_n}{c + (\mu_s - \mu_d)\sigma_n \text{ on optimal fault } \Phi}
\end{aligned} 
\qquad (3) 
$$

In this definition, $T$ is directly governed by the external principal stresses ($\sigma_1$, $\sigma_2$, $\sigma_3$) and could be negative. However, we limit our interest to $0 \le T \le 1$ because we need the rupture to start propagating spontaneously. Therefore, given a value of $T$, the initial traction vector on each point of the non-planar fault can be computed from Equation 3. We also consider that the absolute stress increases with depth due to lithostatic confining pressure as shown in Figure 3b. This condition is applied up to a depth of 12 km, below which we assume a plastic and dissipative condition where the fault strength does not increase any more ($\sigma_p = \sigma_p(z = 12\,km)$) and $D_c$ becomes much longer (Figure 3c). Above 12 km depth, based on the observations discussed in the Section 3, we initially assume $D_c = 80$ cm up to 4 km depth, where $D_c$ begins growing to 2 m at the surface ($z = 0$ km) to account for a dissipative fault zone in the shallow crust that stabilizes rupture propagation, as suggested in several previous studies of rupture dynamics (e.g. Olsen et al., 2009; Aochi and Ulrich, 2015). A larger shallow $D_c$ is also necessary to mitigate the limitations of a homogeneous medium and to emulate the effect of an upper low velocity layer that slows the rupture velocity (e.g. Mikumo et al., 1987). Along-strike variations in $D_c$ suggested by the near-fault ground motions will be discussed later.

## 2.4 Dynamic Rupture and Wave Propagation Numerical Methods

To simulate earthquake dynamic rupture, we adopt a 3D Boundary Integral Equation Method (BIEM) (Aochi et al., 2000) including the mirror source approximation for the free surface (Aochi and Fukuyama, 2002). Although the method is limited to a homogeneous half-space, the portability of this method allows the parametric stress analysis presented later in Section 3.2. In general, a suite of BIEM is useful to consider the scaling issue and fracture problem, as meshing/remeshing is flexible (e.g. Ando et al., 2004; Ide and Aochi, 2005) and to provide different rupture scenarios in different
<figure><img src="image_4.png" /><figcaption>**Figure 3** (a) Mohr-Coulomb diagram for $T = 0.75$. Mohr circles are illustrated for three different depths. The dots on the circles indicate the initial stress applied to each element of the fault model illustrated in Figure 1. It is implicitly assumed that $\sigma_2 = (\sigma_1 + \sigma_3)/2$ corresponds to lithostatic pressure minus hydrostatic pressure as a function of depth. (b) Distribution of the maximum and minimum principal stresses, $\sigma_1$ and $\sigma_3$, and the deviatoric stress $\Delta\sigma = (\sigma_1 - \sigma_3)/2$ along depth for $T = 0.75$. (c) Distribution of critical slip weakening distance $D_c$ along depth. The same parametrization as in Aochi and Ulrich (2015).</figcaption></figure>

fault systems emphasizing the role of fault geometry (e.g. Aochi and Ulrich, 2015; Ando and Kaneko, 2018), while other volumetric methods such as finite element, spectral element and discontinuous Galerkin methods allow the simulations in a more complex medium (e.g. Tago et al., 2012; Jia et al., 2023; Gabriel et al., 2023). Our interest is focused on the ground motions in the vicinity of the fault, at very short distances, where pulse-like waveforms, mostly determined by the close rupture, dominate and where therefore the assumption of a homogeneous medium is a reasonable first approximation. Our standard fault discretization consists of square subelements with a size ($\Delta$s) of 500 m, leading to 546 (along-strike) x 34 (along-depth) = 18 564 subelements. The time step is $\Delta t = \frac{\Delta s}{2V_p} = 0.0417$ s for a total simulation time of about 75 s (1820 steps). Rupture is initiated by a sudden circular crack with radius of 3 km where $\tau_p = \tau_r$ at time $t = 0$, so that a stress drop instantaneously occurs. Once an earthquake scenario is simulated, we use the slip-rate time histories on the fault to compute the ground motion in a second step by means of a 3D Finite Difference Method (FDM) (Aochi and Madariaga, 2003) that solves the elastodynamic equations in a layered half-space (Supplementary material Figure S1). As this procedure is sequential, we can test different crustal structures for the same

rupture scenario treated as a kinematic source model (Supplementary material Figure S2). In the following discussions, we adopt model (b) of Figure S1, hereafter called as our reference model, which has a slight shallow velocity variation and is close enough to the homogeneous model. This model improves the peak values at the near-source stations compared to the homogeneous model, and inhibits the generation of large surface waves behind the rupture front passage at frequencies 0.2-0.4 Hz. Based on the space and time grid sizes reported in Table 1, the maximum resolvable frequency in the FDM simulations is $f_{max} = \frac{V_{s,min}}{(5\Delta s)} = 3.2$ Hz (Levan der, 1988).

# 3 Data Analysis and Simulation Results

## 3.1 Fault Friction Constraint from Strong Motion Data

Eleven accelerometers recorded the earthquake within 3 km from the fault trace. This gives us an unprecedented opportunity to understand some aspects of the rupture front dynamics. Since the peak slip-rate at each fault point is mechanically correlated with the stress breakdown time, $T_c$ (Mikumo et al., 2003; Fukuyama et al., 2003), the peak off-fault velocities can be used to estimate the latter parameter and hence the slip-weakening distance ($D_c$, Equation 1) from displacement records, as proposed by Fukuyama and Mikumo (2007) for the 2000 Tottori and 2002 Denali earthquakes. However, the stress breakdown frequencies (lower bounded by $1/T_c$) that convey information about the dynamic process in the cohesion zone decrease exponentially with distance from the fault in sub-shear rupture earthquakes, making it difficult to estimate $D_c$ reliably (Cruz-Atienza et al., 2009). Only ground motion at fault distances less than about the width of the cohesion zone, $L_c$, is meaningful, what happened in the 2004 Parkfield earthquake (Cruz-Atienza et al., 2009) because the rupture did not reach a steady supershear rupture regime, where conical Mach waves carry such information at much longer distances, as observed for the 1999 Izmit (Figure 1c) and 2002 Denali earthquakes (Cruz-Atienza and Olsen, 2010). In the case of the Mw 7.8 Kahramanmaraş earthquake, apart from a couple of possible supershear episodes, rupture along the fault segment shown in Figure 1b maintained a sub-shear rupture propagation regime (Delouis et al., 2023).

Figure 4a shows the acceleration record at station 4616 projected into the fault-parallel (i.e., X axis) direction (N30°E). This site is located some 20 km west of the epicenter (Figure 1) and only ~2.9 km from the main fault trace. Since the actual rupture initiated on a secondary splay fault before reaching the main EAF (Melgar et al., 2023; Delouis et al., 2023), the major energy burst associated with the rupture front (with Peak Ground Acceleration (PGA) of 580 gal) arrived some 25 s after the first wave arrival, when the rupture front passed right next to the station. This feature of the seismogram repeats in all sites analyzed here (Supplementary material Figure S3), which are located to the
<figure><img src="image_5.png" /><figcaption>**Figure 4** Processing of strong motion data at station 4616 for the estimation of fault cohesive zone parameters. (a) Raw fault-parallel (N30°E) acceleration record, where Te is the estimated arrival time of the main shock wave and T0 is the initial time for further analysis. (b) Velocity window starting at T0 after one integration using an automated baseline correction algorithm and 1 s tapering. Note that Te is clearly defined in the velocity waveform. (c) Velocity and displacement (double integration by the same method) seismograms starting at Te, low-pass filtered at 0.4 Hz and unfiltered. Proxies for the stress breakdown time, Tc, and the slip-weakening distance, $D_c''$, are given at the time of peak velocity (see text).</figcaption></figure>

southwest of station 4616 (Figure 1b). To estimate the stress breakdown time, Tc, we identified the arrival time of the rupture-front shock wave in each seismogram. To this end, we first integrated the acceleration record through an automated baseline correction method (Melgar et al., 2013) to obtain velocity and displacement seismograms. Figure 4b displays the resulting velocigram cut at T0, the initial time 3 s before the main wave arrival time, denoted as Te. After 1s-Tukey tapering, we lowpass filtered the traces at 0.4 Hz. Figure 4c compares the filtered and unfiltered displacement and velocity seismograms starting at Te, where Tc corresponds to the time of the peak velocity and $D_c''$ to the displacement at that moment (Mikumo et al., 2003; Fukuyama et al., 2003). The double prime notation for $D_c''$, introduced by Cruz-Atienza et al. (2009), simply serves to differentiate the value measured on the fault, $D_c'$, from the value measured off the fault, which is subject to wave propagation and free surface effects. The values of Tc and $D_c''$ determined for the other stations with the same procedure are shown in Figure S3 and summarized in Figure 5. The blue curve in Figure 5b (left axis) gathers the $D_c''$ values measured at each site along with an error bar corresponding to an uncertainty of 40% (also valid for Tc), which is a rough estimate obtained from numerical experiments (Cruz-

Atienza et al., 2009). An average $D_c''$ value of 86 +/- 34 cm is reported in the figure legend along with the PGA (red curve, right axis) per site measured as the geometric mean of the peak values on both horizontal components.

To assess whether measured values of Tc and $D_c''$ are representative of the stress drop duration and the associated slip at the rupture front, respectively, we first estimated the width of the cohesion zone (i.e., of the rupture front), Lc, considering both, an average rupture velocity Vr of 3.5 km/s with an uncertainty of 40 %, and the 40 % uncertainty on Tc mentioned in the previous paragraph. Lc values (given by Vr times Tc) incorporating both uncertainties vary between 4 and 12 km along most of the fault (mean value of 9.2±8.3 km between -80 and 30 km), as illustrated by the blue curve in Figure 5c (left axis) with the corresponding error bars. As the rupture nears its end (stations 3131 and 3132), the width of the cohesion zone increases significantly, reaching values above 20 km. Thus, to find out whether the stations are close enough to the fault for $D_c''$ to be representative of Dc, the slip weakening distance (Equation 1), we plotted the ratio between Lc and the distance of each station to the fault trace, D, as a red curve in the same Figure 5c (right axis). Values greater than one (i.e., above the red dotted line) indicate that the sites are located
<figure><img src="image_6.png" /><figcaption>**Figure 5** Estimates of dynamic source parameters from acceleration records within 3 km from the fault trace. (a) Fault surface projection and strong motion stations. (b) Proxy of the slip weakening distance, $D_c''$ (left axis, blue curve), and peak ground acceleration (geometric mean of horizontal components) (right axis, red curve). Note the anti-correlation between the two observables. (c) Width of the rupture front cohesive zone, Lc, assuming an average rupture velocity of 3.5 km/s (left axis, blue curve). Error bars contain 40% uncertainties on rupture velocity and stress breakdown times (see text). The slip-weakening distance, $D_c$, can only be reliably estimated for distances to the fault (D) shorter than Lc. The red curve (right axis) depicts the ratio Lc / D, so sites with values greater than 1 (red dotted line) are likely at a good resolution distance for $D_c$ estimates. Note that all stations are above the resolution threshold.</figcaption></figure>

at distances from the fault less than Lc, the width of the cohesion zone, and therefore that $D_c''$ is likely representative of $D_c$ on the fault (Cruz-Atienza et al., 2009). Since all the stations are above this threshold, then the estimates of $D_c''$ reported in Figure 5b should be a reasonable proxy of the actual values of $D_c$ in fault segments close to the stations. However, Lc was not determined independently of Tc. The breakdown time, Tc, was estimated from seismograms (Figure 4c), so if not well resolved (due to wave propagation effects), then the above Lc estimates are not well resolved either and thus the above exercise is not a rigorous test of $D_c$ resolution. To mitigate such an uncertainty, Fig-

ure S4 shows the distribution of Lc along the fault determined directly from the simulation results of our preferred earthquake model discussed later in Section 4.4. Although highly variable in space (mainly due to rupture speed variations), the mean value in the upper 5 km is Lc = 6.0±4.8 km, which is close to those reported in Figure 5c (blue line) and more than twice the fault distances of all stations. From these arguments, we believe that our estimates of $D_c$ should be reasonable enough.

Possible implications of the along-strike variation of $D_c$ suggested by our results on the earthquake dynamics, along with some energy budget considerations, will be discussed in Section 4.4 from numerical simulations
in light of the observed strong motion.

## 3.2 Uniform stress field analysis

To find reasonable values for the fault prestress condition leading to sustained spontaneous rupture, we first performed a parametric analysis for the optimal fault direction ($\Phi$) and the magnitude of the Mohr circle ($T$) defined in Equation 3). The first question is whether a uniform stress field can explain the rupture extension over 250 km long. Let us focus on the southwestern fault segment. If we consider the tangential direction derived from these Euler poles as the optimal rupture direction, given the discrepancy between that direction with the strike along the fault (Figure 2), then there would be significant inconsistency/uncertainty (larger than 30° at many places) in the construction of the prestress condition. For this reason, we choose to explore systematically different values for such an optimal direction.

<figure><img src="image_7.png" /><figcaption>**Figure 6** Parameter study under a horizontally uniform stress field. Nucleation is set at (a) $X = -50$ km and (b) $X = 0$ km along the fault. The result shows the final magnitude given by simulation. White areas indicate the cases that rupture could propagate far enough beyond $X = -190$ km or $X = -40$ km in each case, respectively. (c) Rupture extension and the surface rupture in the simulation for the selected cases. The star represents the nucleation position for each simulation. Cases 2 and 6 are successful.</figcaption></figure>

We initially set the nucleation point at $X = -50$ km (in the rotated coordinate system; Figure 1b), which is west of the triple junction where the initial splay fault meets the EAF and far enough to the north to mitigate any effect of nucleation on the subsequent rupture propagation in the zone of interest, where seismic stations concentrate. Figure S5 shows an example of the initial conditions for $\Phi = N30°E$ and $T = 0.80$. Although the external principal stress is horizontally uniform, the shear and normal stresses vary along the fault as a function of fault strike because of the non-planar fault geometry. We explored values of $T \in [0.6, 1.0]$ and $\Phi \in$

$[N10°E, N70°E]$ in the parametric analysis depending on the location of the nucleation point. Figure 6a shows the simulation results in terms of the final magnitude. Since we are only looking for the model parameters that allow rupture to extend across the entire fault (i.e., beyond $X = -190$ km), the favorable model space is very limited to the white area. Outside this area, rupture either stops somewhere in the middle of the fault or fails to initiate successfully. In no case did rupture propagate to the right-side, beyond the prominent fault bend, so the prestress condition in that northern segment should be different from that in the southern segment.

To explore the northern segment ($X > -40$ km), where the initial splay fault reaches the main fault, we moved the nucleation point to $X = 0$ km and performed a similar analysis. In this case, we look for ruptures reaching $X = -40$ km, where the fault bends. Again, the favorable model space is minimal as depicted by the white area (Figure 6b). In this northern segment, the stress magnitude $T$ could be slightly lower, indicating that the fault geometry is closer to the optimal fault direction in this part. Figure S6 shows the comparison between the three conditions A, B and C, in which the stress field is extremely high ($T = 0.98$). Nevertheless, none of the conditions (with different hypocenter positions and optimal fault directions) succeeded in producing the rupture length expected between $X = 0$ km and -190 km. It should be noted that cases B and C share the same stress condition but have a different nucleation position. Since the final magnitudes of these cases are different, then the nucleation point at $X = 0$ km is not favorable for the given stress condition. In condition B, the rupture behavior is unusual, as the rupture jumps to around $X = 80$ km, which is an unrealistic scenario for this earthquake.

This parametric study allows us to conclude that a homogeneous stress field orientation across the entire fault cannot explain the rupture extension of the Kahramanmaraş earthquake, indicating that the optimal directions for the fault strike should be around $N60°E$ in the northern segment, and around $N30°E$ in the southern segment. This variation is similar to the great-circle tangential direction deduced from the model of Aktuğ and Kiliçoğlu (2005), but there remains a significant difference in the middle by a difference of more than $45°$ (Figure 2). Given the fault length of over 200 km, it is not surprising that the principal stress direction changes along the fault path, as has already been demonstrated in dynamic rupture simulations for the 1992 Mw 7.3 Landers earthquake (Aochi and Fukuyama, 2002) and the 2008 Mw 7.9 Wenchuan earthquake (Tang et al., 2021). The optimal stress magnitude ($T$) ranges between 0.70 and 0.80, which is within the limits found in previous studies (e.g. Aochi and Ulrich, 2015) and consistent with the values expected to produce near-fault ground motions in accordance with Ground Motion Prediction Equations (Aochi et al., 2017).
## 3.3 Sustained rupture propagation under non-uniform stress field

Since no combination of the model parameters explored in the previous section allowed for a complete earthquake rupture from nucleation at $X = 0$ km to -200 km, here we shall build a consistent model that allows for continuous, large, sustained rupture. To this end, we combined the two preferred models found in the parametric analysis above, namely, $\Phi = \text{N30}^{\circ}\text{E}$ and $T = 0.80$ for $X < -40$ km (Figure 6a) and $\Phi = \text{N60}^{\circ}\text{E}$ and $T = 0.75$ for $X \ge -40$ km (Figure 6b), and placed the nucleation point at $X = 0$ km. Figure 7 shows the initial stresses and fault dynamic parameters along the fault for this two-zone model. Compared to the fault parameterization under a uniform stress field (Figure S5), in this new model the potential stress drop is overall larger, particularly on the northern fault segment. As a result, the two-zone model produced a sustained and complete fault rupture as shown in Figure 8a. Furthermore, the correlation found between irregularities in fault geometry and lateral variations in the peak slip rate and final slip reveal the major role that fault geometry plays even in a simple stress tectonic setting. However, rupture speed was faster than the shear-wave velocity (supershear) over more than 190 km (i.e. between -220 and -30 km with $V_r$ close to 5 km/s), and this cannot explain the observed seismic waves as demonstrated in Figure 9a, where the model-predicted waveforms are far ahead and larger than those observed. The synthetic seismograms show a fast-propagating shock wave, a signature of supershear earthquakes, which is absent in the observations. This is consistent with previous works, which have shown that most of the rupture process of this earthquake took place in a subshear regime (e.g. Melgar et al., 2023; Delouis et al., 2023).

To slow down the rupture process, we further adapt our source model by subdividing the stress field into four zones every 80 km length. Two of them are set with lower stress levels and the southernmost zone redirecting the stress is assumed to arrest the rupture. Figure 8b shows the simulation results from this four-zone model. Although there are still some episodes of supershear rupture, the overall process maintains a sustained subshear regime with rupture velocities around 3.3 km/s that produced a remarkable ground motion prediction when compared to observed seismograms (Figure 9b), particularly for the shock wave arrival times. The amplitude itself does not change significantly between the two models. Previously, Jia et al. (2023) and Gabriel et al. (2023) were also able to reproduce the general characteristics of near-field ground motions from dynamic rupture simulations, including five common stations as our simulations (4624, 4616, 2712, 2718 and 3139). The near fault ground motions are mainly characterized by a dominant pulse corresponding to the passage of the rupture front in all simulations. It is interesting to note that their and our models are well timed at station 2712, while the simulations are generally too early at station 2718. In our simulation it is difficult to reproduce well the fault parallel component at the positions between $X = -100$ km and $X=-150$ km, although it is still good in

<figure><img src="image_8.png" /><figcaption>**Figure 7** Initial condition on the fault plane for 2-zone model. Horizontal axis presents the distance along the fault. From top to bottom, initial shear stress, initial normal stress, fault strength, stress excess required for rupturing, possible stress drop and $D_c$. The contours show the rupture times (see Figure 8b) every 1 s.</figcaption></figure>

the fault perpendicular component. Thus, our simulation could capture the macroscopic rupture propagation well, but some local conditions could be further improved, as suggested by Kidoh et al. (2024), who propose a short pulse generation on this segment near station 3139. This exercise shows that the four-zone model is globally consistent with previous knowledge of the earthquake and observational expectations. For this reason, we will consider the four-zone source model as a reference model for further discussion.

## 3.4 Generalized dynamic source model

The detailed source inversion of the Pazarcık-Kahramanmaraş earthquake introduced by Delouis et al. (2023) reveals that the rupture propagation experienced two localized supershear transients in its southwestern segment along the EAF. This kinematic model benefits from an unprecedented data set next to the fault that captured interesting properties of the rupture process on a local scale, which gives us the opportunity to look for more detailed features of the underlying dynamics. To this end, our four-zone reference model requires further complexity due, admittedly, to residual prestress heterogeneities not accounted for by our principal stress setup.

Figure 8b shows that the reference four-zone model already exhibits two main supershear rupture transients around $X = (-200, -160)$ km and $(-50, -20)$ km. However, they are spatially shifted (about 10-40 km for both cases) when compared to the supershear transients found by Delouis et al. (compare Figure 8b with their Figure 5). This indicates that our reference model is not heterogeneous enough primarily in terms of the
<figure><img src="image_9.png" /><figcaption>**Figure 8** Simulation results of dynamic rupture propagation in cases in which the nucleation at X = 0 km allows rupture to propagate until the left end. The assumption on the stress field is given by the split zones at top. The fault geometry and final slip distribution on the ground surface are illustrated in the middle. Snapshots show the spatio-temporal evolution of the slip rate on the non-planar fault, projected along the X-axis. Stars indicate the nucleation point. (a) Two-zone cases assembling the two better parameter sets from the previous parameter studies. (b) Four-zone cases, adjusted to be comparable to the observed rupture velocity.</figcaption></figure>

stress initial load and friction. After gradually and carefully increasing the number of stress zones along the

fault (see Figure S7 for five- and six-zone models), we
found that the seven-zone model shown in Figure 10a
<figure><img src="image_10.png" /><figcaption>**Figure 9** Comparison of velocity waveforms between the simulated ground motions (orange) and observations (black) at selected stations, whose locations are shown in Figure 1b. Synthetic ground motions are aligned at t = 17 s in the figure. X- and Y-components correspond briefly to fault-parallel and fault-normal components. The maximum and minimum ground velocities are indicated by red and blue dots, open marks for the simulation and solid ones for the observations. Cases (a) and (b) correspond respectively to each case in Figure 8. Low-pass filter is applied up to 1 Hz.</figcaption></figure>

best reproduces the expected overall rupture features including the two supershear transients predicted by the kinematic model (white squares); one around X = -120 km, and the other between X = -70 and -50 km, just southwest of the large, northern fault bending. The stress values *T* for the seven-zone model are indicated in Figure 10c. Comparing to Figure 8b, a relatively low stress zone around X = -20 km was indeed necessary to keep sub-shear rupture propagation in this area, which corresponds to the shadow part of the splay fault where the Mw 7.8 earthquake started. We also found neces-

sary a high stress zone surrounding the fault wrinkle at X = -130 km, which is consistent with the expected strain concentration around that fault geometric irregularity.

Source models tested so far consider an along-strike constant D_c (Figure 7 and bottom of Figure 10a) that corresponds to the average value of our D_c'' estimates determined from the strong motion records shown in Figure 5b. Although the uncertainty of these estimates is large, they feature spatially consistent along-strike variations that may be real to some extent, so let us now evaluate the effects of such D_c variations on the source
<figure><img src="image_11.png" /><figcaption>**Figure 10** Dynamic rupture simulations from the seven-zone stress model for (a) uniform horizontal $D_c$ distribution and for (b) along-strike non-uniform $D_c$ distribution estimated from $D_c''$ (Figure 5). Top, fault geometry and station locations. Second row, along-dip averaged rupture times and local rupture velocities compared with the kinematic model of Delouis et al. (2023). The third and fourth rows display the rupture velocity (normalized by Vs = 3.464 km/s) and maximum slip rate on the fault surface. Bottom, $D_c$ distributions for both models. (c) Summary of the seven-zone stress intensity along the fault for both models.</figcaption></figure>

propagation and radiation. To preserve the rupture initiation process, frictional parameters are unchanged for $X > -35$ km. As for the rest of the fault, while keeping the seven-zone stress distribution and large near-surface $D_c$ values like in all previous simulations (see also Figure S8 for further discussion), for depths between 4 and 12 km we imposed the along-strike linearly interpolated $D_c''$ values shown in Figure 5b as $D_c$ on the fault (see bottom panel of Figure 10b). Simulation result for this case is shown in Figure 10b. Although small, there are some significant differences with the along-strike constant $D_c$ model (Figure 10a). In terms of locally averaged rupture times, both models explain similarly well the inverted kinematic model of Delouis et al. (2023) (white squares). However, the two supershear transients are better captured in the variable $D_c$ model, particularly between -75 and -50 km. Rupture arrest for $X < -175$ km is also better described thanks to larger

$D_c$ estimates at the three westernmost stations, which also bound to lower values the peak slip rates (PSR) in that ending segment. We also find that two of the PSR maxima are in the supershear fault regions around -125 km and -65 km, with depths between 9 to 12 km. In fact, the correlation between average rupture speed and PSR holds along most of the fault surface, as can be seen in the phase diagram shown in Figure 11a. From this figure it is also clear that only the supershear transient around -125 km reached an averaged PSR above 3 m/s, while the other two maxima above this threshold took place in fault segments under subshear ruptures speeds.

Following Díaz-Mojica et al. (2014) and Mirwald et al. (2019), from our dynamic source model (i.e., from the evolution of the shear traction at each fault point) we estimated the radiation efficiency across the fault. Defined as $η_r = \frac{E_r}{(E_r+G)}$, where $E_r$ is the radiated energy
<figure><img src="image_12.png" /><figcaption>**Figure 11** Rupture characteristics for the seven-zone stress model with along-strike non-uniform $D_c$ distribution (Figure 10b). (a) Along-dip averaged peak slip rate as a function of rupture speed along the fault strike. (b) Along-dip averaged radiation efficiency as a function of rupture speed along the fault strike. (c) Final slip distribution and (d) absolute rupture speed.</figcaption></figure>

and G the fracture energy or breakdown work (Husseini, 1977; Venkataraman and Kanamori, 2004; Cocco et al., 2006), this source parameter quantifies how much of the energy available to propagate the rupture is radiated compared to the stress breakdown work retained in the source. Theoretical models for the three fracture modes predict that $\eta_r$ grows with rupture speed so that it is low ($\eta_r < 0.4$) for deep and tsunami earthquakes and high ($0.4 < \eta_r < 0.8$) for shallow intraplate ruptures (Venkataraman and Kanamori, 2004; Mirwald et al., 2019). Figure 11b shows the distribution of $\eta_r$ along the fault as a function of locally-averaged rupture speed normalized by the shear wave velocity. As expected

from theory and similarly to the PSR (Figure 11a), radiation efficiency is overall linearly related to the rupture velocity and spans over a wide range going from 0.1 to 0.9 for $0.3 < Vr/Vs < 1.1$ (excluding rupture arrest), with highest values above 0.7 (excluding rupture initiation) within both supershear rupture transients (i.e., around -100 and -70 km). In contrast, the fault segment exhibiting the largest PSR around -175 km (Figures 10b and 11a), while rupturing on subshear regime ($Vr/Vs \approx 0.75$), it was relatively inefficient with $\eta_r \approx 0.45$, which is explained by the low prestress level and the higher $D_c$ value in the final segment of the fault (see the lower panels of Figure 10b).
Figures 11c and 11d show the final slip and the absolute rupture speed distributions for the along-strike variable $D_c$ model already presented in Figure 10b. The slip distribution presents a segmentation controlled mainly by the fault geometry, as previously noted for our four-zone reference model (Figure 8b), with two slip deficit areas close to those determined by Delouis et al. (2023) around $X$ >-50 km and $X$ <-100 km (see blue shades in Figure 10). Figure 12 shows the synthetic seismograms at near-fault stations from the 7-zone model with and without the horizontal variation in $D_c$, corresponding to Figure 10a and 10b, respectively. Our model, though, also has a slip deficit around -130 km related to the wrinkle-like fault irregularity, where the rupture struggles to propagate. As mentioned earlier, the rupture velocity is highly variable, especially for $150 < X < -100$ km, where the wrinkle-like geometric barrier is found and where the observed PGVs (and PGAs, Figure 5b) are maximum, as shown with the blue solid curve in Figure 13 (left axis) at stations 3139, 3145 and 3144, and where the radiation efficiency overcomes 0.8 (Figure 11b). Thus, these strong-motion maxima appear to be related to the supershear transient around -120 km at station 3144 and to the fault geometry irregularity at stations 3139 and 3145. Although smaller at some sites, the model-predicted PGVs for frequencies smaller than 3 Hz (blue dotted curve) follow the same general pattern observed along the fault, with two maxima at stations 3145 and 3144, where the rupture undergoes remarkable speed changes (see orange curve). As for the PSR (Figure 11a), comparison of the average rupture velocity between 5 and 12 km depth (orange curve, right axis) with the PGVs reveals a noteworthy correlation, where the largest seismic bursts are very close to fault segments with fast rupture (e.g., stations 3144 and 2712) or where the fault undergoes a sharp geometric change (i.e., a sort of kink where large amplitude diffracted waves are expected; e.g. station 3145). In contrast, our model is unable to explain the largest observed PGV at station 3139, where rupture slows down right after clearing the wrinkle-like barrier.

The aspects of the rupture process described above can be better appreciated in Figure 14 (and Supplementary Movie S1), where the slip rate evolution is shown along with the three-dimensional fault geometry. After initiating bilaterally where the splay fault meets the EAF, the rupture propagates southwestward in a subshear regime to cross the first major fault bending, where the slip rate decreases significantly (see Movie S1). About 30 km ahead ($t = 26$ s), the rupture undergoes the first supershear transient where $D_c$ decreases (Figure 10b) just before encountering the second major fault bending, where it slows down again. Around 44 s, the earthquake reaches the second supershear transient where radiation efficiency is maximum (Figure 11b) (and where $D_c$ is minimum, see the bottom panel of Figure 10b) and where the model predicts the maximum PGVs (and PGAs, see Figure 5b) in agreement with the data (see Figure 13). The wrinkle-like fault barrier brutally slows the rupture to resume velocity in the subshear regime along the final, relatively flat segment of the fault. Thus, our dynamic source model shows

how fault geometry played a preponderant role during the Kahramanmaraş earthquake and how variations in rupture speed are responsible for the observed along-strike variations of the observed strong motions (Figures 12, 13, and S9 in Supplementary material). Along-strike variations of $D_c$ estimated directly from the rupture front shock wave improved the model predictions in terms of both the expected locations of the supershear rupture transients and the spatial distribution of the observed PGVs. Our modelling framework represents the first attempt to integrate the non-planar fault geometry, a heterogeneous regional stress field and a three-dimensional variation of $D_c$ over the source, all three factors derived from geological, tectonic and seismological observations. The resulting rupture process is thus a juxtaposition of these factors, which are in fact interrelated.

## 4 Discussion and Conclusions

We have simulated the dynamic rupture propagation and the near-fault ground motions of the February 6th 2023 01:17 UTC Pazarcık-Kahramanmaraş, Turkey, earthquake along a non-planar fault structure determined from satellite interferograms assuming a regional principal stress field and depth dependent slip-weakening friction estimated directly from near-fault strong motion records. To better understand the relationship between the dynamic rupture parameters and the near-fault ground motion along the fault strike, we focused on the ~200 km length, best-instrumented south-western segment of the earthquake. To this purpose, we adopted a modeling framework previously introduced for scenario earthquakes along the North Anatolian fault (Aochi and Ulrich, 2015). Namely, the initial stress on the fault is loaded by the external principal stresses while the depth-dependent rupture criterion and slip-weakening friction govern the rupture process. By assuming an orientation of the optimal fault plane with respect to the principal stresses of N30°E south of latitude ~37.4° and N60°E north of it, this simple framework was able to explain the rupture extent over 200 km with an average rupture velocity of 2.8 km/s, so that the arrival times and low-frequency ($f < 3$ Hz) amplitude of shock waves recorded at 11 stations along the fault strike were also explained. To this end, the intensity of the stress field should be non-uniform (i.e., should vary along the fault strike) in at least four (and optimally seven) distinct segments with lengths ranging between 20 and 80 km, conditions that produced significant variations of the rupture process. In no case was a uniform stress field along the entire source able to reproduce either the actual extent of the earthquake or the observed ground motion.

Careful analysis of strong motions next to the fault allowed us to constrain friction along the causal fault plane. Passing near the stations, the shock wave associated with the rupture front revealed a spatially consistent along-strike variation of the critical slip-weakening distance, $D_c$, ranging between 0.6 and 1.2 m with an average value of 0.86 +/- 0.34 m. The cohesion zone width also featured variations in space going from ~3 to
<figure><img src="image_13.png" /><figcaption>**Figure 12** Comparison of velocity waveforms from the simulations of 7-zone stress model without and with horizontal variation of $D_c$. The simulated ground motions (orange) and observations (black) are plotted along fault strike (See also caption of Figure 9). Cases (a) and (b) correspond respectively to each case in Figure 10. Low-pass filter is applied up to 1 Hz.</figcaption></figure>

~12 km over a ~100 km fault segment (Figure 5 bottom panel), with minimum values where recorded PGAs exceeded 640 gal (Figure 5 middle panel). Following a similar strategy, Ding et al. (2023) and Yao and Yang (2023) estimated $D_c$ from seismic records that led to higher estimates than ours. Unlike our approach, where the rupture-front shock wave was isolated prior to the double integration of acceleration via a baseline correction method, Ding et al. (2023) and Yao and Yang (2023) determined $D_c''$ from long displacement time series that suffer from the well-known baseline drift inherent in inertial accelerometers. In their procedure, these authors also ignored the effect of the free surface in the es-

timation of $D_c''$ which amplifies the motion in such a way that the factor 2 introduced by Fukuyama and Mikumo (2007) is unnecessary, as demonstrated by Cruz-Atienza et al. (2009). For these reasons, $D_c$ estimates in the aforementioned studies are likely to be significantly affected by factors unrelated to the rupture process. We acknowledge, however, that since they explored a sufficiently wide range of $D_c$ values, their main conclusions should not be significantly affected by this problem. The very same issue arising both from the double integration baseline drift and the misleading factor 2 for estimating $D_c''$ is also present in the work by He et al. (2024). Beyond the factor 2 used erroneously by
<figure><img src="image_14.png" /><figcaption>**Figure 13** Along-strike distribution of rupture speed (right axis) and off-fault observed (solid) and synthetic (dotted) fault parallel PGVs (left axis). Rupture speed is averaged between 5-12 km depth. Seismograms were low-pass filtered at 3 Hz. The station locations are shown on top.</figcaption></figure>

these and other authors in the literature, we must emphasize that the double integration of accelerations is a very delicate matter that, despite using a baseline correction, often leads to large displacement errors that grow rapidly as the record elapses (e.g. see Melgar et al., 2013). In our case, estimates of $D_c''$ derive mostly from the first 3 s (or less) of the main shock wave and within distances from the fault smaller than the local dimension of the rupture cohesion zone, which is the most reliable $D_c$ resolution criterion (Cruz-Atienza et al., 2009). Therefore, unlike the previous works mentioned and despite other sources of error intrinsic to such $D_c$ determination strategy (see Cruz-Atienza et al. (2009)), we believe that our $D_c''$ estimates (Figure 5 middle panel), which are smaller than those reported for this earthquake in previous works, should be related (and thus reliable) to some extent to the actual stress breakdown process along the fault.

The kinematic source model determined by Delouis et al. (2023) allowed us to study some details about the earthquake dynamics. Relatively small perturbations of the prestress level along the fault together with the along-strike variation of $D_c$ inferred from the strong motions, allowed us to explain satisfactorily the rupture times including the two supershear rupture transients found in the inverted model. The analysis of the energy partitioning at the rupture front revealed that the maximum PGAs observed come from a fault segment where the rupture propagated in the supershear regime. In that segment the radiation efficiency reached its maximum value above 0.8 and it is where $D_c$ is minimum around 0.6 m. The high PGAs recorded at stations 3139 and 3145, just southwest of that supershear transient, could be due to diffracted wave radiation where the fault geometry features a wrinkle-like irregularity (i.e., a sort of fault kink). Overall, the rupture velocity of our source

model is highly variable locally, with values normalized by the shear wave speed between 0.3 and 1.1, and fluctuations of $\eta_r$ between 0.1 and 0.9 that clearly correlate with rupture speed. It is because of these large local variations that the model can explain the most prominent, overall features of the earthquake, such as the kinematically inverted rupture times and the main seismic energy bursts.

In this paper, we have emphasized the importance of $D_c$ variations both vertically and horizontally on the rupture process together with the stress field, which is reasonable from the point of view of energy balance (e.g. Mikumo et al., 2003). The rupture process and seismic radiation are also influenced by the medium heterogeneiety such as the low-velocity shallow layer (e.g. Inoue and Miyatake, 1997; Abdelmeguid et al., 2025) or velocity reduction in the fault zone (e.g. Ma and Andrews, 2010; Schliwa et al., 2025). The impact of the shallow layer is very clear in the wave propagation even from the kinematic point of view, as also shown in Figure S2 (Supplementary material). The rupture process strongly depends on the condition we impose. In our configuration, seismic radiation from the shallow part is moderated by the mechanical fault consideration of a large $D_c$, which slows down the cohesive zone process. Such fault (zone) rheology in the shallow part should certainly be an important factor to be considered in seismic hazard assessment.

In summary, we could build a reasonable dynamic source model for the Mw 7.8 Pazarcık-Kahramanmaraş earthquake constrained from the near-fault seismic observations. The model is consistent with the unprecedented source inversion by Delouis et al. (2023) in terms of rupture propagation and captures the main features of the recorded strong motions at eleven stations within 3 km from the fault trace. We found that fault geome-
<figure><img src="image_15.png" /><figcaption>**Figure 14** Evolution of slip velocity along the three-dimensional fault geometry predicted by our preferred model described in Figure 10b and dissected from Figures 11 and 12. See text.</figcaption></figure>

try played a major role in rupture propagation and seismic wave radiation, and that none of the uniform prestress field assumptions can reproduce the 200 km rupture extension. At least four, slightly different prestress-intensity zones along the fault strike and two principal stress orientations are necessary, which implies that the stress field in the crust is heterogeneous at the earthquake scale. Furthermore, along-strike variations of $D_c$ estimated directly from the rupture front shock wave improved the model predictions in terms of both the

expected locations of the supershear rupture transients and the spatial distribution of the observed PGVs. Radiation efficiency and rupture speed are highly variable along the fault at local scale and correlate with each other, as expected from rupture mechanics theory, so that the largest PGA values are found where radiation efficiency is maximum (above 0.8) along one supershear transient. Observational insights into lateral friction and prestress heterogeneities may thus have important implications for further modeling scenario earthquakes
in the globe.

## Acknowledgements

We thank the editor Alice Gabriel and two reviewers, whose comments were very helpful to improve our discussion and manuscript. We thank AFAD and KOERI for the useful information on the earthquakes and data. We benefit the funding from the French ANR (Agence National de la Recherche) project E-CITY (Near-fault observation and simulation of earthquake ground motion in an urban environment) under grant ANR-21-CE22-0020 and UNAM PAPIIT project number IN111524. The simulations benefit the computational resources from the French national computing center GENCI/T-GCC and GENCI/Idris under grant A0130406700 and A0150406700.

## Data and code availability

Acceleration data are provided by Disaster and Emergency Management Authority of Turkey (AFAD), publicly available at https://tdvms.afad.gov.tr/. The numerical codes (BIEM-FDM) are available on https://doi.org/10.5281/zenodo.1472238 and https://doi.org/10.5281/zenodo.10225171. The geometry model used in this study is established with help of Bryan Raimbault and Romain Jolivet (unpublished work) and the numerical model used in this study is available in supplementary material.

## Competing interests

No competing interest is declared.

## References

Abdelmeguid, M., Zhao, C., Yalcinkaya, E., Gazetas, G., Elbanna, A., and Rosakis, A. Dynamics of episodic supershear in the 2023 M7.8 Kahramanmaraş/Pazarcık earthquake, revealed by near-field records and computational modeling. *Communications Earth & Environment*, 4(1), Dec. 2023. doi: 10.1038/s43247-023-01131-7.

Abdelmeguid, M., Elbanna, A., and Rosakis, A. Ground motion characteristics of subshear and supershear ruptures in the presence of sediment layers. *Geophysical Journal International*, 240 (2):967–987, 2025. https://doi.org/10.1093/gji/ggae422.

Adda-Bedia, M. and Madariaga, R. Seismic Radiation from a Kink on an Antiplane Fault. *Bulletin of the Seismological Society of America*, 98(5):2291–2302, Oct. 2008. doi: 10.1785/0120080003.

Aktuğ, B. and Kiliçoğlu, A. Establishment of regional reference frames for detecting active deformation areas in Anatolia, IAG-IABO-IASPO Joint Assembly, 2005. accessible on ResearchGate.

Ando, R. and Kaneko, Y. Dynamic rupture simulation reproduces spontaneous multifault rupture and arrest during the 2016 Mw 7.9 Kaikoura earthquake. *Geophysical Research Letters*, 45(23), Dec. 2018. doi: 10.1029/2018gl080550.

Ando, R., Tada, T., and Yamashita, T. Dynamic evolution of a fault system through interactions between fault segments. *Journal of Geophysical Research: Solid Earth*, 109(B5), May 2004. doi: 10.1029/2003jb002665.

Aochi, H. and Fukuyama, E. Three-dimensional nonplanar simulation of the 1992 Landers earthquake. *Journal of Geophysical Research: Solid Earth*, 107(B2), Feb. 2002. doi: 10.1029/2000jb000061.

Aochi, H. and Madariaga, R. The 1999 İzmit, Turkey, earthquake: Nonplanar fault structure, dynamic rupture process and strong ground motion. *Bulletin of the Seismological Society of America*, 93(3):1249–1266, June 2003. doi: 10.1785/0120020167.

Aochi, H. and Ulrich, T. A probabilistic earthquake scenario near Istanbul determined from dynamic simulations. *Bulletin of the Seismological Society of America*, 105(3):1468–1475, May 2015. doi: 10.1785/0120140283.

Aochi, H., Fukuyama, E., and Matsu'ura, M. Spontaneous Rupture Propagation on a Non-planar Fault in 3D Elastic Medium. *Pure and Applied Geophysics*, 157(11):2003–2027, Dec. 2000. doi: 10.1007/pl00001072.

Aochi, H., Douglas, J., and Ulrich, T. Stress accumulation in the Marmara Sea estimated through ground-motion simulations from dynamic rupture scenarios. *Journal of Geophysical Research: Solid Earth*, 122(3):2219–2235, Mar. 2017. doi: 10.1002/2016jb013790.

Barbot, S., Luo, H., Wang, T., Hamiel, Y., Piatibratova, O., Javed, M. T., Braitenberg, C., and Gurbuz, G. Slip distribution of the February 6, 2023 Mw 7.8 and Mw 7.6, Kahramanmaraş, Turkey earthquake sequence in the East Anatolian Fault Zone. *Seismica*, 2(3), Apr. 2023. doi: 10.26443/seismica.v2i3.502.

Chen, K., Wei, G., Milliner, C., Dal Zilio, L., Liang, C., and Avouac, J.-P. Super-shear ruptures steered by pre-stress heterogeneities during the 2023 Kahramanmaraş earthquake doublet. *Nature Communications*, 15(1), Aug. 2024. doi: 10.1038/s41467-024-51446-y.

Cocco, M., Spudich, P., and Tinti, E. On the mechanical work absorbed on faults during earthquake ruptures. American Geophysical Union, 2006. doi: 10.1029/gm170.

Cruz-Atienza, V. and Olsen. Supershear Mach-waves expose the fault breakdown slip. *Tectonophysics*, 493(3–4):285–296, Oct. 2010. doi: 10.1016/j.tecto.2010.05.012.

Cruz-Atienza, V. M. and Virieux, J. Dynamic rupture simulation of non-planar faults with a finite-difference approach. *Geophysical Journal International*, 158(3):939–954, July 2004. doi: 10.1111/j.1365-246x.2004.02291.x.

Cruz-Atienza, V. M., Virieux, J., and Aochi, H. 3D Finite-Difference Dynamic-Rupture Modelling Along Non-Planar Faults. *GEOPHYSICS*, 72(5):SM123 – SM137, Sept. 2007. doi: 10.1190/1.2766756.

Cruz-Atienza, V. M., Olsen, K. B., and Dalguer, L. A. Estimation of the Breakdown Slip from Strong-Motion Seismograms: Insights from Numerical Experiments. *Bulletin of the Seismological Society of America*, 99(6):3454 – 3469, Nov. 2009. doi: 10.1785/0120080330.

Das, S. and Aki, K. A numerical study of two-dimensional spontaneous rupture propagation. *Geophysical Journal International*, 50(3):643 – 668, Sept. 1977. doi: 10.1111/j.1365-246x.1977.tb01339.x.

Delouis, B., van den Ende, M., and Ampuero, J.-P. Kinematic Rupture Model of the 6 February 2023 Mw 7.8 Türkiye Earthquake from a Large Set of Near-Source Strong-Motion Records Combined with GNSS Offsets Reveals Intermittent Supershear Rupture. *Bulletin of the Seismological Society of America*, 114(2):726–740, Nov. 2023. doi: 10.1785/0120230077.

Ding, X., Xu, S., Xie, Y., Van den Ende, M., Premus, J., and Ampuero, J.-P. The sharp turn: Backward rupture branching during the 2023 Mw 7.8 Kahramanmaraş (Türkiye) earthquake. *Seismica*, 2 (3), Nov. 2023. doi: 10.26443/seismica.v2i3.1083.
Dunham, E. and Archuleta, R. Evidence for a Supershear Transient during the 2002 Denali Fault Earthquake. *Bulletin of the Seismological Society of America*, 94(6B):S256–S268, Dec. 2004. doi: 10.1785/0120040616.

Díaz-Mojica, J., Cruz-Atienza, V. M., Madariaga, R., Singh, S. K., Tago, J., and Iglesias, A. Dynamic Source Inversion of the M6.5 Intermediate-Depth Zumpango Earthquake in central Mexico: a Parallel Genetic Algorithm. *Journal of Geophysical Research: Solid Earth*, 119(10):7768–7785, Oct. 2014. doi: 10.1002/2013jb010854.

Eberhart-Phillips, D., Haeussler, P. J., Freymueller, J. T., Frankel, A. D., Rubin, C. M., Craw, P., Ratchkovski, N. A., Anderson, G., Carver, G. A., Crone, A. J., Dawson, T. E., Fletcher, H., Hansen, R., Harp, E. L., Harris, R. A., Hill, D. P., Hreinsdóttir, S., Jibson, R. W., Jones, L. M., Kayen, R., Keefer, D. K., Larsen, C. F., Moran, S. C., Personius, S. F., Plafker, G., Sherrod, B., Sieh, K., Sitar, N., and Wallace, W. K. The 2002 Denali Fault Earthquake, Alaska: A Large Magnitude, Slip-Partitioned Event. *Science*, 300(5622): 1113–1118, May 2003. doi: 10.1126/science.1082703.

Fukuyama, E. and Mikumo, T. Slip-weakening distance estimated at near-fault stations. *Geophysical Research Letters*, 34(9), May 2007. doi: 10.1029/2006gl029203.

Fukuyama, E., Mikumo, T., and Olsen, K. Estimation of the Critical Slip-Weakening Distance: Theoretical Background. *Bulletin of the Seismological Society of America*, 93(4):1835–1840, Aug. 2003. doi: 10.1785/0120020184.

Gabriel, A.-A., Ulrich, T., Marchandon, M., Biemiller, J., and Rekoske, J. 3D Dynamic Rupture Modeling of the 6 February 2023, Kahramanmaraş, Turkey Mw 7.8 and 7.7 Earthquake Doublet Using Early Observations. *The Seismic Record*, 3(4):342–356, Oct. 2023. doi: 10.1785/0320230028.

Gallovič, F. and Valentová, ☒ Broadband Strong Ground Motion Modeling Using Planar Dynamic Rupture With Fractal Parameters. *Journal of Geophysical Research: Solid Earth*, 128(6), May 2023. doi: 10.1029/2023jb026506.

Guatteri, M., Mai, P., Beroza, G., and Boatwright, J. Strong Ground-Motion Prediction from Stochastic-Dynamic Source Models. *Bulletin of the Seismological Society of America*, 93(1):301–313, Feb. 2003. doi: 10.1785/0120020006.

Güvercin, S. E., Karabulut, H., Konca, A. Ó., Doğan, U., and Ergintav, S. Active seismotectonics of the East Anatolian Fault. *Geophysical Journal International*, 230(1):50–69, Feb. 2022. doi: 10.1093/gji/ggac045.

Harris, R., Dolan, J., Hartleb, R., and Day, S. The 1999 İzmit, Turkey, earthquake: A 3D dynamic stress transfer model of intra-earthquake triggering. *Bulletin of the Seismological Society of America*, 92(1):245–255, Feb. 2002. doi: 10.1785/0120000825.

Harris, R. A. and Day, S. M. Dynamics of fault interaction: Parallel strike-slip faults. *Journal of Geophysical Research: Solid Earth*, 98(B3):4461–4472, Mar. 1993. doi: 10.1029/92jb02272.

Harris, R. A., Barall, M., Archuleta, R., Dunham, E., Aagaard, B., Ampuero, J. P., Bhat, H., Cruz-Atienza, V., Dalguer, L., Dawson, P., Day, S., Duan, B., Ely, G., Kaneko, Y., Kase, Y., Lapusta, N., Liu, Y., Ma, S., Oglesby, D., Olsen, K., Pitarka, A., Song, S., and Templeton, E. The SCEC/USGS Dynamic Earthquake Rupture Code Verification Exercise. *Seismological Research Letters*, 80(1):119–126, Jan. 2009. doi: 10.1785/gssrl.80.1:119.

He, Z., Zhang, Z., Wang, Z., and Wang, W. Slip-weakening distance and energy partitioning estimated from near-fault recordings during the 2023 Mw 7.8 Türkiye-Syria earthquake. *Tectonophysics*, 885:230424, Aug. 2024. doi: 10.1016/j.tecto.2024.230424.

Husseini, M. I. Energy balance for motion along a fault. *Geophysical Journal International*, 49(3):699–714, June 1977. doi: 10.1111/j.1365-246x.1977.tb01313.x.

Ida, Y. Cohesive force across the tip of a longitudinal-shear crack and Griffith's specific surface energy. *Journal of Geophysical Research*, 77(20):3796–3805, July 1972. doi: 10.1029/jb077i020p03796.

Ide, S. and Aochi, H. Earthquakes as multiscale dynamic ruptures with heterogeneous fracture surface energy. *Journal of Geophysical Research: Solid Earth*, 110(B11), Nov. 2005. doi: 10.1029/2004jb003591.

Inoue, T. and Miyatake, T. 3-D simulation of near-field strong ground motion: Basin edge effect derived from rupture directivity. *Geophysical Research Letters*, 24(8):905–908, Apr. 1997. doi: 10.1029/97gl00619.

Jia, Z., Jin, Z., Marchandon, M., Ulrich, T., Gabriel, A.-A., Fan, W., Shearer, P., Zou, X., Rekoske, J., Bulut, F., Garagon, A., and Fialko, Y. The complex dynamics of the 2023 Kahramanmaraş, Turkey, Mw 7.8-7.7 earthquake doublet. *Science*, 381(6661):985–990, Sept. 2023. doi: 10.1126/science.adi0685.

Kame, N. and Yamashita, T. Dynamic nucleation process of shallow earthquake faulting in a fault zone. *Geophysical Journal International*, 128(1):204–216, Jan. 1997. doi: 10.1111/j.1365-246x.1997.tb04081.x.

Kame, N., Rice, J. R., and Dmowska, R. Effects of prestress state and rupture velocity on dynamic fault branching. *Journal of Geophysical Research: Solid Earth*, 108(B5), May 2003. doi: 10.1029/2002jb002189.

Kaneko, Y., Avouac, J.-P., and Lapusta, N. Towards inferring earthquake patterns from geodetic observations of interseismic coupling. *Nature Geoscience*, 3(5):363–369, Apr. 2010. doi: 10.1038/ngeo843.

Karabulut, H., Güvercin, S. E., Hollingsworth, J., and Konca, A. Ó. Long silence on the East Anatolian Fault Zone (Southern Turkey) ends with devastating double earthquakes (6 February 2023) over a seismic gap: implications for the seismic potential in the Eastern Mediterranean region. *Journal of the Geological Society*, 180(3), May 2023. doi: 10.1144/jgs2023-021.

Kase, Y. and Kuge, K. Numerical simulation of spontaneous rupture processes on two non-coplanar faults: effects of geometry on fault interaction. *Geophysical Journal International*, 135(3):911–922, Dec. 1998. doi: 10.1046/j.1365-246x.1998.00672.x.

Kidoh, T., Nagano, M., Miyakoshi, K., Arai, K., and Sato, T. A characterized fault model with shallow subsurface slips for reproducing large-amplitude pulse-like velocity waveforms and permanent displacements observed very close to the seismic fault for the 2023 Kahramanmaraş earthquake (Mw7.8) in Türkiye. *Journal of Structural and Construction Engineering (Transactions of AIJ)*, 89(822):840–849, Aug. 2024. doi: 10.3130/aijs.89.840.

King, G. and Nábělek, J. Role of Fault Bends in the Initiation and Termination of Earthquake Rupture. *Science*, 228(4702):984–987, May 1985. doi: 10.1126/science.228.4702.984.

Levander, A. R. Fourth-order finite-difference P-SV seismograms. *GEOPHYSICS*, 53(11):1425–1436, Nov. 1988. doi: 10.1190/1.1442422.

Ma, S. and Andrews, D. J. Inelastic off-fault response and three-dimensional dynamics of earthquake rupture on a strike-slip fault. *Journal of Geophysical Research: Solid Earth*, 115(B4), Apr. 2010. doi: 10.1029/2009jb006382.

Mahmoud, Y., Masson, F., Meghraoui, M., Cakir, Z., Alchalbi, A., Yavasoglu, H., Yönlü, O., Daoud, M., Ergintav, S., and Inan, S. Kinematic study at the junction of the East Anatolian fault and the Dead Sea fault from GPS measurements. *Journal of Geodynamics*, 67:30–39, July 2013. doi: 10.1016/j.jog.2012.05.006.

Melgar, D., Bock, Y., Sanchez, D., and Crowell, B. W. On robust and
reliable automated baseline corrections for strong motion seismology. *Journal of Geophysical Research: Solid Earth*, 118(3): 1177–1187, Mar. 2013. doi: 10.1002/jgrb.50135.

Melgar, D., Taymaz, T., Ganas, A., Crowell, B., Öcalan, T., Kahraman, M., Tsironi, V., Yolsal-Çevikbil, S., Valkaniotis, S., Irmak, T. S., Eken, T., Erman, C., Özkan, B., Dogan, A. H., and Altuntaş, C. Sub- and super-shear ruptures during the 2023 Mw 7.8 and Mw 7.6 earthquake doublet in SE Türkiye. *Seismica*, 2(3), Mar. 2023. doi: 10.26443/seismica.v2i3.387.

Mikumo, T., Hirahara, K., and Miyatake, T. Dynamical fault rupture processes in heterogeneous media. *Tectonophysics*, 144(1–3): 19–36, Dec. 1987. doi: 10.1016/0040-1951(87)90006-0.

Mikumo, T., Olsen, K., Fukuyama, E., and Yagi, Y. Stress-Breakdown Time and Slip-Weakening Distance Inferred from Slip-Velocity Functions on Earthquake Faults. *Bulletin of the Seismological Society of America*, 93(1):264–282, Feb. 2003. doi: 10.1785/0120020082.

Mirwald, A., Cruz-Atienza, V. M., Díaz-Mojica, J., Iglesias, A., Singh, S. K., Villafuerte, C., and Tago, J. The 19 September 2017 (Mw7.1) Intermediate-Depth Mexican Earthquake: A Slow and Energetically Inefficient Deadly Shock. *Geophysical Research Letters*, 46(4):2054–2064, Feb. 2019. doi: 10.1029/2018gl080904.

Nakata, T., Shimazaki, K., Suzuki, Y., and Tsukuda, E. Fault Branching and Directivity of Rupture Propagation. *Journal of Geophysics (Chigaku Zasshi)*, 107(4):512–528, 1998. doi: 10.5026/j-geography.107.512.

Oglesby, D. D., Archuleta, R. J., and Nielsen, S. B. Dynamics of dip-slip faulting: Explorations in two dimensions. *Journal of Geophysical Research: Solid Earth*, 105(B6):13643–13653, June 2000. doi: 10.1029/2000jb900055.

Olsen, K. B., Day, S. M., Dalguer, L. A., Mayhew, J., Cui, Y., Zhu, J., Cruz-Atienza, V. M., Roten, D., Maechling, P., Jordan, T. H., Okaya, D., and Chourasia, A. ShakeOut-D: Ground motion estimates using an ensemble of large earthquakes on the southern San Andreas fault with spontaneous rupture propagation. *Geophysical Research Letters*, 36(4), Feb. 2009. doi: 10.1029/2008gl036832.

Provost, F., Karabacak, V., Malet, J.-P., Van der Woerd, J., Meghraoui, M., Masson, F., Ferry, M., Michéa, D., and Pointal, E. High-resolution co-seismic fault offsets of the 2023 Türkiye earthquake ruptures using satellite imagery. *Scientific Reports*, 14(1), Mar. 2024. doi: 10.1038/s41598-024-55009-5.

Reilinger, R., McClusky, S., Vernant, P., Lawrence, S., Ergintav, S., Cakmak, R., Ozener, H., Kadirov, F., Guliev, I., Stepanyan, R., Nadariya, M., Hahubia, G., Mahmoud, S., Sakr, K., ArRajehi, A., Paradissis, D., Al-Aydrus, A., Prilepin, M., Guseva, T., Evren, E., Dmitrotsa, A., Filikov, S. V., Gomez, F., Al-Ghazzi, R., and Karam, G. GPS constraints on continental deformation in the Africa-Arabia-Eurasia continental collision zone and implications for the dynamics of plate interactions. *Journal of Geophysical Research: Solid Earth*, 111(B5), May 2006. doi: 10.1029/2005jb004051.

Reitman, N. G., Briggs, R., Barnhart, W., Jobe, J. A., Duross, C. B., Hatem, A. E., Gold, R., Akçiz, S., Koehler, R. DMejstrik, J. D., and Collett, C. Fault Rupture Mapping of the 6 February 2023 Kahramanmaraş, Türkiye, Earthquake Sequence from Satellite Data (ver. 1.1, February 2024), 2023. doi: 10.5066/P98517U2.

Schliwa, N., Gabriel, A.-A., and Ben-Zion, Y. Shallow fault zone structure affects rupture dynamics and ground motions of the 2019 Ridgecrest sequence to regional distances. Jan. 2025. doi: 10.31223/x5n412.

Tago, J., Cruz-Atienza, V. M., Virieux, J., Etienne, V., and Sánchez-Sesma, F. J. A 3D hp-adaptive discontinuous Galerkin method for modeling earthquake dynamics. *Journal of Geophysical Research: Solid Earth*, 117(B9), Sept. 2012. doi: 10.1029/2012jb009313.

Tang, R., Zhu, S., and Gan, L. Dynamic Rupture Simulations of the 2008 7.9 Wenchuan Earthquake: Implication for Heterogeneous Initial Stress and Complex Multifault Geometry. *Journal of Geophysical Research: Solid Earth*, 126(12), Nov. 2021. doi: 10.1029/2021jb022457.

Venkataraman, A. and Kanamori, H. Observational constraints on the fracture energy of subduction zone earthquakes. *Journal of Geophysical Research: Solid Earth*, 109(B5), May 2004. doi: 10.1029/2003jb002549.

Wang, Z., Zhang, W., Taymaz, T., He, Z., Xu, T., and Zhang, Z. Dynamic Rupture Process of the 2023 Mw 7.8 Kahramanmaraş Earthquake (SE Türkiye): Variable Rupture Speed and Implications for Seismic Hazard. *Geophysical Research Letters*, 50(15), Aug. 2023. doi: 10.1029/2023gl04787.

Yao, S. and Yang, H. Rupture phase in near-fault records of the 2023 Turkey Mw 7.8 earthquake. Apr. 2023. doi: 10.31223/x51662.

The article *Rupture Dynamics and Near-Fault Ground Motion of the Mw7.8 Kahramanmaraş, Turkey earthquake of February 6, 2023* © 2025 by H. Aochi is licensed under **CC BY 4.0**.