Revisiting the 2020 Mw 6.8 Elaziğ, Türkiye Earthquake with Physics-Based 3D Numerical Simulations Constrained by Geodetic and Seismic Observations

Abstract

Dynamic rupture simulations of earthquakes offer crucial insights into the physical mechanisms of driving fault slip and seismic hazards. By incorporating non-planar fault models that accurately represent subsurface structures, this study provides a realistic depiction of the rupture processes of the 2020 Mw 6.8 Elazığ, Türkiye earthquake, influenced by geometric complexities. Initially, we determined its coseismic slip on the non-planar fault using near-field strong motion and InSAR observations. Subsequently, we established the heterogeneous initial stress on the fault plane based on the coseismic slip and integrated it into the dynamic rupture modeling to assess physics-based ground motion and seismic hazards. The numerical simulations utilized the curved grid finite-difference method (CGFDM), which effectively models rupture dynamics with heterogeneities in fault geometry, initial stress, and other factors. Our synthetic surface deformation and seismograms align well with the observational data obtained from InSAR and seismic instruments. We observed localized occurrences of supershear rupture during fault propagation. Furthermore, the intensity distribution we simulated closely aligns with the actual observations. These findings highlight the critical role of source heterogeneity in seismic hazard assessment, advancing our understanding of fault dynamics and enhancing predictive capabilities.

1. Introduction

On 24 January 2020, at 17:55 (UTC), a magnitude Mw 6.8 earthquake struck Elaziğ in eastern Türkiye, resulting in a total of 42 fatalities, around 1600 injuries, and severe damage to the towns of Sivrice, Doğanyol, Pütürge, and Malatya [1]. The event ruptured the Pütürge segment of the East Anatolian fault (EAF), a 600 km-long left-lateral strike-slip fault that intersects the right-lateral North Anatolian fault (NAF) [2,3] and connects the Dead Sea fault (DSF) to the Karlıova Triple Junction. Historically, the EAF experienced several strong earthquakes (Mw > 6.5) in the 1800s, but only small to moderate events occurred over the past century until the 2020 Elaziğ earthquake [4,5,6,7]. Following this event, two devastating earthquakes with magnitude 7.8 and 7.6 (Figure 1) occurred on the EAF and its side branch on 6 February 2023. The stress interactions between past earthquakes on the EAF may affect the occurrence and rupture propagation of the 2023 earthquake doublet [8], notably the most adjacent 2020 Mw 6.8 Elaziğ earthquake. A precise estimate of source characteristics of the 2020 event plays a critical role in understanding rupture behavior and strain accumulation along the EAF.

Several studies [7,10,11,12,13,14,15,16,17] have investigated the rupture kinematics of 2020 Mw 6.8 Elaziğ earthquake through joint or individual inversions of geodetic (i.e., InSAR and GNSS) and seismic (i.e., strong motion and teleseismic) observations. These kinematic models generally reveal a similar heterogeneous slip distribution with a shallow slip deficit and predominant southwestward rupture propagation. However, they differ significantly in their representations of the causative fault geometry, including single planar fault, e.g., [7,11,12,15,16,17], two planar faults, e.g., [10,14], and a curved fault [13]. In fact, the EAF is known for its complex geometry, comprising multiple main segments characterized by bends, pull-apart basins, and compressional structures, e.g., [2]. Field investigations [18] and aftershock patterns, e.g., [14], suggested that the rupture segment of the 2020 Elaziğ earthquake has complex geometry structure with two obvious bends along the strike direction, highlighting structural complexities that likely control the rupture pattern, e.g., [10,13]. The influence of intricate fault geometries on rupture dynamics has been observed in several earthquakes, such as the 2008 Mw 7.9 Wenchuan earthquake [19], the 2019 Ridgecrest earthquake sequence [20], the 2021 Mw 7.4 Maduo earthquake [21], and the 2023 Mw 7.8 and 7.6 Türkiye doublet earthquakes [22,23].

Although kinematic models can resolve the spatial and temporal evolution of slip on the predefined faults, they do not provide insight into the underlying physical mechanisms driving fault slip nor explain the variation in slip magnitudes. Physics-based 3D numerical simulations for earthquake rupture dynamics and ground motion, which incorporate complex non-planar fault systems, rough surface topography, and the heterogeneous media structures, are increasingly essential for advancing our understanding of earthquake physics and improving seismic hazard assessment, e.g., [24,25,26]. Gallovič et al. [27] performed dynamic source inversion of the 2020 Elaziğ earthquake using strong motion recordings, revealing intricate rupture dynamics such as relatively weak nucleation and small mean rupture velocity (~2 km/s). However, their model assumed a vertical planar fault, neglecting the realistic geometric complexity of the rupture segment, particularly variations. While their inversion successfully estimated stress and frictional parameters critical for physics-based hazard assessment, it required significant computational resources, with simulations taking a couple of weeks. This highlights the need for methods that address fault complexity while maintaining computational efficiency.

In this study, we first resolved the kinematic rupture model of the 2020 Mw 6.8 Elaziğ, Türkiye earthquake from a joint inversion of InSAR and strong motion data, using a linear multi-time window approach. Based on the coseismic stress drop constraints, we then conducted dynamic rupture scenario simulations that aim to the slip characteristics of the preferred kinematic model and to assess its mechanical viability. Finally, the results of the dynamic rupture simulations were imported to simulate strong ground motion using physics-based numerical methods that couple 3D complexities in topography and medium and dynamic sources. Our work highlights that physics-based 3D numerical simulations play a vital role in the seismic hazard assessment.

2. Materials and Methods

This section provides a comprehensive overview of the datasets and methodologies employed in this study. It begins with a detailed description of the data sources and preprocessing steps. Subsequently, the kinematic inversion method and the dynamic rupture simulation process are outlined.

2.1. InSAR and Strong Motion Data

Interferometric synthetic aperture radar (InSAR) is a powerful technique for detecting ground surface changes between two satellite overpasses. For this study, we collected two coseismic interferograms from the European Space Agency (ESA) Sentinel-1 satellites: ascending track A116 (20200121-20200127) and descending track D123 (20200122-20200128) (Figure 1). These SAR acquisitions were processed using the JPL/Caltech InSAR Scientific Computing Environment (ISCE) software version 2 [28]. The detailed processing procedure is described in [17]. It is important to note that surface displacements derived from the InSAR data include 3 to 7 days of postseismic deformation, which may affect the modeling of the coseismic phase. To reduce the computational burden in the inversion, the unwrapped interferograms were downsampled using the quadtree algorithm [29,30], which retains more data points in areas with rapid surface displacement changes and fewer points in areas with smoother variations. As a result, we obtained 860 and 958 data points for the A116 and D123 tracks, respectively.

Additionally, we collected 14 strong motion recordings (Figure 1) from the Disaster and Emergency Management Presidency of Turkey (AFAD). The original acceleration recordings were corrected and integrated to obtain velocity waveforms, which were bandpass filtered between 0.02 and 0.25 Hz and resampled at a 5 Hz sampling rate. The processed waveforms were then trimmed to a 90 s time window for the kinematic inversion.

2.2. Kinematic Inversion

We determined the kinematic rupture process through the joint inversion of InSAR displacements and strong motion waveforms. InSAR observations (Figure 1) reveal significant displacements northwest of the fault, suggesting a northwest-dipping fault consistent with the aftershock distribution. The 2020 Elaziğ earthquake ruptured a portion of the Pütürge segment on the East Anatolian Fault Zone [10,13]. The main surface trace of the Pütürge segment is characterized by a relatively continuous sinusoidal trend interrupted by minor bends and step-overs that do not exceed one kilometer in width [2].

Based on a detailed analysis of InSAR and pixel-offset tracking displacements, along with aftershock locations, Ragon et al. [13] built a causative fault with two structural bends. Following their results, we constructed a realistic curved fault geometry with two bends along the strike direction (Figure 2). The fault has an average strike of N239°E and dips northwest at 75°, consistent with previous studies, e.g., [12,17]. The fault plane was discretized into 35 along-strike and 14 downdip subfaults, each measuring 1.5 by 1.5 km.

The mainshock hypocenter (39.1126°E, 38.3821°N, 8.9 km) and origin time (17:55:16 UTC) were adopted from the solution from the Kandilli Observatory and Earthquake Research Institute (KOERI). Green’s functions were calculated using the frequency-wavenumber integration approach [31] based on a 1-D layered velocity model (Table 1) derived from the CRUST1.0 global crustal model. The source time function at each subfault was parameterized using five symmetric triangular functions with 50% overlap and a rise time of 2 s. The rake angle on each subfault was constrained to vary within a 90° range centered on 0° to match the left-lateral strike-slip faulting mechanism.

The rupture speed, determined through trial and error, was finalized at 2.2 km/s, consistent with Melgar et al. [11]. The linear multi-time window approach was employed with temporal and spatial smoothing constraints and solved using a non-negative least-square method implemented in the open source MudPy v1.0 software [32]. The optimal smoothing factor was determined by the L-curve criterion. After extensive trial-and-error tests, a relative weighting ratio of 5:1 was assigned to balance the fitting of strong motion and InSAR datasets.

2.3. Dynamic Rupture Simulations

We employed the curved-grid finite-difference method (CGFDM) to simulate the dynamic rupture of the 2020 Elaziğ earthquake. This method, introduced by Zhang et al. [33], utilizes a curvilinear coordinate system to discretize irregular spatial domains, solving the first-order velocity-stress equations and conducting dynamic simulations by applying fault boundary conditions based on friction laws. The CGFDM provides exceptional flexibility for simulating complex faults, making it particularly well suited for modeling the geometrically intricate faults commonly observed in nature. Moreover, it seamlessly incorporates irregular free surfaces into the simulations. Zhang et al. [34] further enhanced this method with GPU-based parallel acceleration, enabling efficient large-scale dynamic rupture simulations.

The initial stress distribution is a key parameter influencing the rupture process. In this study, the initial stress on the fault plane is divided into two components: dynamic stress (${\tau }_n {\mu }_d$), determined by the regional background stress field and fault geometry, and static stress drop ($∆\tau$). To derive the dynamic stress distribution, the traction on the fault plane is calculated based on the regional background stress and the fault’s normal vectors. The static stress drop is estimated using an analytical approach [35] and integrated with the regional background stress.

3. Results

3.1. Preferred Kinematic Model

Through the kinematic inversion of the finite fault source, we derived the fault plane slip model, the corresponding moment release rate, and the rupture evolution process, as shown in Figure 2. Starting from the hypocenter, the rupture predominantly propagated southwestward, lasting over 20 s, and the total seismic moment is 2.62 × 10¹⁹ N·m, corresponding to moment magnitude of Mw 6.88. The kinematic results reveal two main rupture zones. The first zone, located approximately 8.0 km below the surface and closer to the source, exhibits a maximum slip of about 2.35 m, while the second one, located farther and deeper (up to 15.0 km), shows a maximum slip of approximately 1.93 m. Notably, the rupture did not reach the surface, consistent with field investigations by Tatar et al. [18].

The simulation results for the InSAR data and strong motion waveforms were compared with the observed data and are shown in Figure 3 and Figure 4, respectively. The simulated InSAR displacements can generally match the observed data, yielding a root mean square error (RMSE) of 0.021 m and 0.021 m for the ascending A116 and descending D123 tracks, respectively. In terms of strong ground motion, the simulated waveforms demonstrate a good fit, successfully capturing the primary components of the signal.

3.2. Dynamic Rupture Model

To perform the dynamic rupture simulation of the Elaziğ earthquake, we required a suitable fault geometry, regional stress field, friction parameters, and velocity model. For the fault geometry, we adopted a non-planar fault model similar to that used in the kinematic inversion. The inverted slip distribution was used to constrain the initial stress.

The maximum principal compressive stress in the study area is oriented N12°E [36], the stress shape ratio, R = 0.457, is defined by the equation R = (S_v − S_H)/(S_h − S_H) [3], where S_H, S_h, and S_v denote the maximum horizontal principal stress, the minimum horizontal principal stress, and the vertical principal stress, respectively. Through a trial-and-error method, it was determined that S_H, S_h, and S_v are 60 MPa, 25 MPa, and 46 MPa, respectively. To account for stress reduction at the free surface, we implemented a gradual decrease to zero at depths less than 5 km (Figure 5a). This regional stress serves as the background heterogeneous stress induced by the non-planar fault geometry.

The friction law plays a critical role in dynamic rupture simulations. In this study, we adopted the linear slip-weakening law [37,38]. After extensive trial-and-error tests, we set the static friction coefficient to 0.40, the dynamic friction coefficient as 0.17, and the slip-weakening distance as 0.31 m. Figure 6 illustrates the distribution of normal and shear stresses on the fault plane, derived using the stress field decomposition method. The friction coefficient, stress distribution, and other specific parameters involved in the dynamic rupture are provided in Table 2.

In our dynamic simulation, the grid size is set to 50 m, the time step to 0.0027 s, and the total simulation time is 20 s. Figure 7 presents snapshots of the slip rate, illustrating the rupture propagation process. The rupture initially propagates up-dip before unilaterally advancing towards the southwest after nucleation. It slows down upon encountering the first bend and accelerates as it reaches the strongest asperity, releasing about 2.5 m of slip. Following the rupture of the largest asperity, it continues down-dip until it stops. At time t = 9.99 s and t = 12.15 s, local supershear ruptures are observed when the rupture reaches asperities. Overall, however, the rupture speed varies significantly, with a relatively small mean value of 1.85 km/s. Our results are consistent with the dynamic inversion conducted by Gallovič et al. [27].

Figure 8a presents the slip distribution obtained by the dynamic simulation, with a maximum slip of 2.5 m. The rupture does not reach the free surface due to a gradually increasing Dc and C₀ value near the free surface, which is consistent with the absence of surface rupture observed. Notably, it is evident that dynamic results differ from the kinematic inversion results. The maximum slip in the dynamic simulation is located slightly more to the right. The total released seismic moment is 1.95 × 10¹⁹ N·m, corresponding to moment magnitude of Mw 6.79, which is consistent with our kinematic inversion result and previous studies, e.g., [17,27]. The moment rate release process is more concentrated and slightly delayed (Figure 8b), which can be attributed to artificial nucleation during the dynamic rupture simulation with a peak value being slightly larger than the kinematic results.

Given the critical influence of frictional parameters on dynamic rupture simulations, we also conducted a parameter sensitivity analysis to examine how variations in the dynamic friction coefficient (${\mu }_d$) and critical slip-weakening distance ($D_c$) influence simulation outcomes. In this analysis, we maintained the static friction coefficient (${\mu }_s$) constant while systematically varying ${\mu }_d$, as the friction coefficient drop (${\mu }_s-{\mu }_d$) rather than the individual values of ${\mu }_s$ or ${\mu }_d$ serves as the primary driver of rupture propagation. The results are presented in Figures S1 and S2. Our preferred dynamic model ($D_c=0.31$ and ${\mu }_d=0.17$) shows strong agreement with the kinematic result in terms of final slip, moment rate function and overall data fittings. Models with $D_c<0.31$ or ${\mu }_d<0.17$ exhibit maximum slip exceeding 2.5 m and tend to rupture to the surface, which is not consistent with field observations. Conversely, models with $D_c\ge 0.33$ exhibit difficulty in rupture initiation and a delayed moment release rate, while those with ${\mu }_d\ge 0.18$ shows partial rupture, with rupture termination occurring prematurely. Consequently, these models were excluded through a trial-and-error method.

The simulation results of InSAR data and strong motion waveforms, presented in Figure 9 and Figure 10, provide an assessment of the selected dynamic model. The simulated displacement values can generally match the observed data, resulting in a small root mean square error (RMSE) of 0.021 m and 0.024 m for the ascending A116 and descending D123 InSAR tracks, respectively. For the strong-motion waveforms, the station locations are indicated in Figure 1. The slip rates derived from the dynamic rupture were convolved with 90 s time-window Green’s functions, precomputed using the generalized reflection and transmission coefficient method [39,40] in a 1D layered medium. As a result, synthetic seismograms were generated at the selected stations.

The observed data were preprocessed using an automatic baseline correction method [41,42] and filtered within both the observed and synthetic seismic records in the frequency range of 0.02 to 0.25 Hz. As shown in Figure 10, the synthetic waveforms exhibit high consistency with the observed waveforms for most stations. Some differences may be attributed to limitations of the current velocity model, in which discrepancies begin to appear between the two sets of waveforms. This difference may be attributed to the limitations of the current velocity model, which lacks structural details for higher frequencies exceeding 0.1 Hz (corresponding to shear wave wavelengths of approximately 20–30 km).

Additionally, we identified significant mismatches in waveform fitting for the three components of stations 2301 and 2308, particularly in the East and Up direction. This suggests that the velocity structure in this region may deviate considerably from the actual geological structure or that the source parameters used in the simulation may not fully reflect earthquake source characteristics. Further refinement of the velocity model and source parameters is necessary to address these discrepancies and enhance the simulations accuracy.

3.3. Simulated Ground Motion

After selecting an appropriate dynamic model, we utilized the CGFDM3D-EQR platform [43] to conduct numerical simulations of strong ground motion for the Elaziğ earthquake. This algorithm is based on CPU-GPU heterogeneous parallel architecture and effectively addresses three-dimensional complexities. It integrates preprocessing modules for terrain, material properties, and finite faults. Its high efficiency and robustness make it particularly suitable for earthquake simulations, intensity calculation, emergency response, and disaster assessment.

The simulations were performed on a local server equipped with 8 NVIDIA Tesla A100 GPUs. The computational domain covered an area of 200 km × 200 km × 80 km, with a spatial resolution of 200 m. The total simulation time was 90.0 s, with a time interval of 0.014 s. The 1D velocity model (CRUST1.0), and the Shuttle Radar Topography Mission (SRTM) 90 m data were used to ensure accuracy in modeling strong ground motion. Compared to empirical ground-motion predictions (e.g., ShakeMaps), snapshots of the seismic wave propagation show the significant rupture directivity effect (Figure 11).

We simulated seismic wave propagation over a total time window of 90 s to calculate the strong ground motion based on the kinematic and dynamic sources, respectively. The horizontal peak ground velocity (PGVh) distribution is primarily influenced by the fault strike. Figure 12 presents a comparison of simulated PGVh maps obtained using the dynamic source model (left panel) and the kinematic source model (right panel). The white contours represent PGVh levels, which are used to derive seismic intensity. Both models exhibit a similar overall pattern, with higher PGVh values near the fault and a gradual decrease in amplitude with increasing distance from the fault. Using the calculated PGVh distribution, we applied the Modified Mercalli Intensity (MMI) [44] scale to assess the earthquake intensity across the study area. The results indicate a diffusion of seismic intensity along the fault strike, extending predominantly from south to north.

Despite the similarities, the two source models exhibit notable differences in detail. Near the epicenter, the dynamic source model predicts a maximum intensity of VIII, concentrated along the fault zone. The overall intensity distribution shows that the region predominantly experiences intensity levels VI, VII, VIII, with levels VII and above mainly observed in the cities of Puturge and Doganyol. In contrast, the kinematic source model estimates a maximum intensity of only level VII, with a much narrower intensity distribution. This discrepancy arises because the kinematic source model is derived from filtered waveform inversion, which suppresses high-frequency information. Consequently, the absence of high-frequency subevents reduces the peak slip rate, leading to lower PGVh values and underestimated intensity levels. From this perspective, the dynamic source model provides intensity estimates that better align with observations, offering a more accurate and reliable framework for seismic intensity assessment.

4. Discussion

We developed a fault model to address the complex subsurface structure of the 2020 Turkey earthquake. In previous studies, fault models have often been simplified, with most adopting planar geometries [7,10,11,12,13,14,15,16,17]. However, natural faults commonly exhibit complex geometries, including bends, branches, and step-overs. In this study, we employed a non-planar fault model in both the kinematic inversion and dynamic rupture simulations, which allowed for a more realistic representation of the fault’s actual geometry and provided better insights into how fault geometry influences rupture propagation. The selected fault model closely aligns with the intricate structure of the East Anatolian Fault (EAF) and reveals detailed propagation characteristics, including two significant bends along the strike, consistent with the observations by [14].

Our kinematic model exhibits a predominantly unilateral rupture propagation towards southwestward direction, which aligns with most of published kinematic models, e.g., [10,11,14,17]. In contrast, Xu et al. [16] concluded that the rupture was a bilateral process with a slow rupture velocity of ~1.6 km/s. The difference may stem from their assumption of a more southwestern hypocenter. Our model indicates two asperities with a maximum slip of ~2.35 m, which is in good agreement with the results of Cheloni and Akinci [12] (~2.3 m) and Pousse-Beltran et al. [10] (~2.4 m). However, the maximum slip in our model is smaller than that reported by Melgar et al. [11] (~3.5 m) and Ragon et al. [13] (~3.5 m) but larger than that of Taymaz et al. [7] (~0.7 m) who relied exclusively on teleseismic observations. These discrepancies can be attributed not only to differences in data constraints but also to variations in the smoothing factors used in the inversion. Our dynamic model further supports that the rupture propagates from northeast to southwest at a relatively slow mean rupture velocity of ~1.85 km/s, suggesting that the event may have occurred on an immature fault. Gallovič et al. [27] conducted dynamic source inversion of the 2020 Elaziğ earthquake on a vertical planar fault, revealing rupture dynamics similar to those observed in our study. Additionally, Chen et al. [17] employed dynamic modeling to investigate the complex rupture behaviors by incorporating heterogeneous initial stresses on the simplified planar fault. Compared with their models [17,27], our study fully accounts for realistic geometric complexity along the strike, providing a more intuitive understanding of the rupture process.

Our kinematic inversion and dynamic rupture simulations complement each other, offering distinct advantages. The kinematic inversion effectively constrains the stress distribution necessary for dynamic rupture simulations, while the fully dynamic rupture model produces more physically consistent results. Due to its higher resolution, the dynamic model captures finer details of the rupture process, such as localized supershear rupture, which are challenging to resolve with kinematic inversion alone.

Moreover, the intensity distributions derived from the dynamic source model not only align more closely with observed post-earthquake intensity data but also highlight the potential for improving seismic hazard assessment. By accurately capturing the physical processes underlying rupture propagation and ground motion, dynamic models offer a pathway toward more reliable predictions of ground shaking intensity. Such advancements could enhance the design of earthquake-resistant infrastructure, improve disaster response planning, and contribute to a deeper understanding of fault mechanics, ultimately fostering a more resilient society in seismically active regions. Future studies should continue to refine these models and explore the influence of varying geological conditions on rupture dynamics, further enhancing our predictive capabilities in the face of seismic hazards.

5. Conclusions

This work presents a comprehensive investigation into the mechanics of the Mw 6.8 Elazığ earthquake in Turkey, integrating kinematic source inversion and dynamic rupture simulations. Our findings reveal that kinematic models effectively constrain plausible dynamic rupture scenarios, while dynamic simulations evaluate the mechanical feasibility of these kinematic slip characteristics. Initially, we performed a kinematic source inversion utilizing InSAR and strong ground motion data to delineate the finite fault characteristics of the earthquake. Following this, we constructed a dynamic model that incorporated the complexities of 3D topography, medium properties, and source dynamics, enabling a comprehensive simulation of the entire process of vibration generation, propagation, and arrest. Our approach, which emphasizes dynamic source inversions for determining the spatial distribution of initial stress and friction parameters, significantly enhances the reliability of seismic hazard assessments. Ultimately, this work enriches our understanding of earthquake mechanics and highlights the necessity of a multifaceted modeling approach to improve predictions of ground shaking and inform effective disaster preparedness strategies in seismically active regions.

Author Contributions

Conceptualization, Z.H. and Y.Z.; methodology, Z.H., W.W. and Z.W.; software, Z.H., W.W. and Z.W.; validation, Z.H., Y.Z. and Z.Z.; formal analysis, Z.H., Y.Z. and Z.Z.; investigation, Z.H. and Y.Z.; resources, Z.H.; data curation, Z.H.; writing—original draft preparation, Y.Z. and Z.H.; writing—review and editing, Y.Z., Z.H., W.W., Z.W., T.C.S. and Z.Z.; visualization, Z.H. and Y.Z.; supervision, Z.Z.; project administration, Z.Z.; funding acquisition, Z.Z. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Original Sentinel-1 InSAR interferograms were acquired and processed by European Space Agency Copernicus program (https://scihub.copernicus.eu, last accessed on 1 October 2024) and retrieved from Alaska Satellite Facility Distributed Active Archive Center and processed and archived online by Chen et al. 2020 (https://doi.org/10.5281/zenodo.4031927, last accessed on 1 December 2024). Strong motion data were downloaded from the Disaster and Emergency Management Authority Presidential of Earthquake Department (AFAD, https://tadas.afad.gov.tr/event-detail/8071, last accessed on 1 December 2024). Active faults were downloaded from the European Database of Seismogenic Faults (EDSF, https://doi.org/10.6092/INGV.IT-SHARE-EDSF, last accessed on 1 December 2024). The relocated seismicity and aftershock location catalogue were provided by Güvercin et al. 2022 (https://doi.org/10.5281/zenodo.5220633, last accessed on 1 December 2024). The SRTM 90m DEM data were downloaded from CGIAR Consortium for Spatial Information (CGIAR-CSI, https://srtm.csi.cgiar.org/srtmdata, last accessed on 1 October 2024). The CRUST 1.0 model was downloaded from https://igppweb.ucsd.edu/~gabi/crust1.html, last accessed on 1 October 2024. The figures were plotted using the Generic Mapping Tools version 6 [45], Python 3.10 and MATLAB R2021b software.

Acknowledgments

We are grateful to the editors and reviewers for their helpful comments and suggestions. We sincerely thank Xu Sun for proofreading the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Tectonic setting and datasets of the study area. (a) Gray stars and beachballs show epicenters and focal mechanisms of the historical earthquakes [7]. Red triangles indicate the locations of strong motion stations. The surface trace of the modeling fault is shown as the thick red line and the mapped faults are plotted as black lines [2]. The maximum and minimum horizontal principal stresses are depicted by red and blue arrows, respectively. The relocated aftershocks [7] of the first 2 months are shown as dark blue-color-filled dots. The inset box at top left shows the main fault systems around Türkiye and plate motion at a broader view (NAFZ, North Anatolian Fault Zone and EAFZ, East Anatolian Fault Zone). The brown lines denote the surface ruptured traces of the 2023 Türkiye earthquake doublet. The black arrows indicate direction of tectonic motion relative to Eurasia [9]. Coseismic InSAR line-of-sight displacements of the (b) ascending track A116 and (c) descending track D123. Positive values indicate moving away from the satellite.

Figure 2. Slip distribution and moment rate function of the preferred kinematic model. The yellow star indicates the mainshock hypocenter.

Figure 3. Observed and modeled down-sampled InSAR line-of-sight displacements from the preferred kinematic model and their residuals for the descending track (D123) and ascending track (A116). Positive values indicate moving toward the satellite.

Figure 4. Comparison between the observed (black) and synthetic (red) strong motion velocity waveforms from the kinematic inversion. All records are bandpass filtered in the frequency band of 0.02~0.25 Hz. The left numbers indicate the station names and the right numbers below waveforms denote the respective maximum amplitudes.

Figure 5. Schematic plots to illustrate the distributions of (a) tectonic principal stresses (SH [greater than] Sv [greater than] Sh), (b) critical slip-weakening distance (Dc), and (c) frictional cohesion (C0) varying with depth.

Figure 6. Heterogeneous initial stress distributions on the nonplanar fault. Shear stress in along strike direction (a) and in downdip direction (b) and normal stress (c). (d) Distribution of the ratio (μ0) between the shear stress and the normal stress.

Figure 7. Snapshots of the strike-slip rate on the non-planar fault obtained from the dynamic rupture simulation. The red star denotes the nucleation location.

Figure 8. Final slip distribution (a) and moment rate release (b) derived from the dynamic rupture simulation.

Figure 9. Observed and modeled down-sampled InSAR line-of-sight displacements from the dynamic rupture simulation and their residuals for the descending track (D123) and ascending track (A116). Positive values indicate moving toward the satellite.

Figure 10. Comparison between the observed (black) and synthetic (red) strong motion velocity waveforms from the dynamic rupture simulation. All records are bandpass filtered in the frequency band of 0.02~0.25 Hz. The left numbers indicate the station names and the right numbers below waveforms denote the respective maximum amplitudes.

Figure 11. Snapshots of the seismic wave propagation ((a): east–west velocity component, (b): north–south velocity component) obtained from physics-based 3D numerical simulations using the dynamic source model.

Figure 12. Comparison of the simulated horizontal peak ground velocity (PGVh) maps using the dynamic source model (left panel) and the kinematic source model (right panel). The white lines are the PGVh contours which are used for intensity conversion.

Table 1

One-dimensional velocity model used in this study.

Depth (km)	Vs (km/s)	Vp (km/s)	Density (g/cm3)
0.0	1.07	2.50	2.11
1.11	3.52	6.00	2.72
18.75	3.68	6.30	2.79
28.29	3.82	6.60	2.85
39.62	4.48	8.07	3.33

Table 2

Parameters used in the dynamic rupture simulations.

Parameters	Value
Static friction coefficient, ${\mu }_s$	0.40
Dynamic friction coefficient, ${\mu }_d$	0.17
Critical slip-weakening distance, $D_c$ (m)	0.31
Frictional cohesion, $C_0$ (MPa)	0.10
Element size, DH (m)	50
Nucleation patch radius, r (km)	2.0
Time step, DT (s)	0.0027
Maximum horizontal stress, SH (MPa)	60.0
Minimum horizontal stress, Sh (MPa)	25.0
Vertical stress, Sv (MPa)	46.0
SH Azimuth	N12°E

Supplementary Materials

The following supporting information can be downloaded at: www.mdpi.com/article/10.3390/rs17040720/s1. Figure S1: Sensitivity analysis of variations in the critical slip-weakening distance ($D_c$). Figure S2: Sensitivity analysis of variations in the dynamic friction coefficient (${\mu }_d$).

References

1. Çetin, K.Ö.; Ilgaç, M.; Can, G.; Çakır, E.; Söylemez, B. Sivrice-Elazig-Turkey Earthquake Reconnaissance Report. DesignSafe-CI; 24 January 2020; [DOI: https://dx.doi.org/10.17603/ds2-9jz1-e287]

2. Duman, T.Y.; Emre, Ö. The East Anatolian Fault: Geometry, segmentation and jog characteristics. Geol. Soc. Lond. Spec. Publ.; 2013; 372, pp. 495-529. [DOI: https://dx.doi.org/10.1144/SP372.14]

3. Yilmaz, H.; Over, S.; Ozden, S. Kinematics of the East Anatolian Fault Zone between Turkoglu (Kahramanmaras) and Celikhan (Adiyaman), eastern Turkey. Earth Planets Space; 2006; 58, pp. 1463-1473. [DOI: https://dx.doi.org/10.1186/BF03352645]

4. Ambraseys, N.N.; Jackson, J.A. Faulting associated with historical and recent earthquakes in the Eastern Mediterranean region. Geophys. J. Int.; 1998; 133, pp. 390-406. [DOI: https://dx.doi.org/10.1046/j.1365-246X.1998.00508.x]

5. Tan, O.; Taymaz, T. Active tectonics of the Caucasus: Earthquake source mechanisms and rupture histories obtained from inversion of teleseismic body waveforms. Postcollisional Tectonics and Magmatism in the Mediterranean Region and Asia; Dilek, Y.; Pavlides, S. Geological Society of America: Boulder, CO, USA, 2006.

6. Taymaz, T.; Eyidoğan, H.; Jackson, J. Source parameters of large earthquakes in the East Anatolian Fault Zone (Turkey). Geophys. J. Int.; 1991; 106, pp. 537-550. [DOI: https://dx.doi.org/10.1111/j.1365-246X.1991.tb06328.x]

7. Taymaz, T.; Ganas, A.; Yolsal-Çevikbilen, S.; Vera, F.; Eken, T.; Erman, C.; Keleş, D.; Kapetanidis, V.; Valkaniotis, S.; Karasante, I. et al. Source Mechanism and Rupture Process of the 24 January 2020 Mw 6.7 Doğanyol–Sivrice Earthquake obtained from Seismological Waveform Analysis and Space Geodetic Observations on the East Anatolian Fault Zone (Turkey). Tectonophysics; 2021; 804, 228745. [DOI: https://dx.doi.org/10.1016/j.tecto.2021.228745]

8. Chen, J.; Liu, C.; Dal Zilio, L.; Cao, J.; Wang, H.; Yang, G.; Göğüş, O.H.; Zhang, H.; Shi, Y. Decoding Stress Patterns of the 2023 Türkiye-Syria Earthquake Doublet. J. Geophys. Res. Solid Earth; 2024; 129, e2024JB029213. [DOI: https://dx.doi.org/10.1029/2024JB029213]

9. Nocquet, J.-M. Present-day kinematics of the Mediterranean: A comprehensive overview of GPS results. Tectonophysics; 2012; 579, pp. 220-242. [DOI: https://dx.doi.org/10.1016/j.tecto.2012.03.037]

10. Pousse-Beltran, L.; Nissen, E.; Bergman, E.A.; Cambaz, M.D.; Gaudreau, É.; Karasözen, E.; Tan, F. The 2020 6.8 Elazığ (Turkey) Earthquake Reveals Rupture Behavior of the East Anatolian Fault. Geophys. Res. Lett.; 2020; 47, e2020GL088136. [DOI: https://dx.doi.org/10.1029/2020GL088136]

11. Melgar, D.; Ganas, A.; Taymaz, T.; Valkaniotis, S.; Crowell, B.W.; Kapetanidis, V.; Tsironi, V.; Yolsal-Çevikbilen, S.; Öcalan, T. Rupture kinematics of 2020 January 24 Mw 6.7 Doğanyol-Sivrice, Turkey earthquake on the East Anatolian Fault Zone imaged by space geodesy. Geophys. J. Int.; 2020; 223, pp. 862-874. [DOI: https://dx.doi.org/10.1093/gji/ggaa345]

12. Cheloni, D.; Akinci, A. Source modelling and strong ground motion simulations for the 24 January 2020, Mw 6.8 Elazığ earthquake, Turkey. Geophys. J. Int.; 2020; 223, pp. 1054-1068. [DOI: https://dx.doi.org/10.1093/gji/ggaa350]

13. Ragon, T.; Simons, M.; Bletery, Q.; Cavalié, O.; Fielding, E. A Stochastic View of the 2020 Elazığ Mw 6.8 Earthquake (Turkey). Geophys. Res. Lett.; 2021; 48, e2020GL090704. [DOI: https://dx.doi.org/10.1029/2020GL090704]

14. Konca, A.Ö.; Karabulut, H.; Güvercin, S.E.; Eskiköy, F.; Özarpacı, S.; Özdemir, A.; Floyd, M.; Ergintav, S.; Doğan, U. From Interseismic Deformation with Near-Repeating Earthquakes to Co-Seismic Rupture: A Unified View of the 2020 Mw6.8 Sivrice (Elazığ) Eastern Turkey Earthquake. J. Geophys. Res. Solid Earth; 2021; 126, e2021JB021830. [DOI: https://dx.doi.org/10.1029/2021JB021830]

15. Lin, X.; Hao, J.; Wang, D.; Chu, R.; Zeng, X.; Xie, J.; Zhang, B.; Bai, Q. Coseismic Slip Distribution of the 24 January 2020 Mw 6.7 Doganyol Earthquake and in Relation to the Foreshock and Aftershock Activities. Seismol. Res. Lett.; 2020; 92, pp. 127-139. [DOI: https://dx.doi.org/10.1785/0220200152]

16. Xu, J.; Liu, C.; Xiong, X. Source Process of the 24 January 2020 Mw 6.7 East Anatolian Fault Zone, Turkey, Earthquake. Seismol. Res. Lett.; 2020; 91, pp. 3120-3128. [DOI: https://dx.doi.org/10.1785/0220200124]

17. Chen, K.; Zhang, Z.; Liang, C.; Xue, C.; Liu, P. Kinematics and Dynamics of the 24 January 2020 M 6.7 Elazig, Turkey Earthquake. Earth Space Sci.; 2020; 7, e2020EA001452. [DOI: https://dx.doi.org/10.1029/2020EA001452]

18. Tatar, O.; Sözbilir, H.; Koçbulut, F.; Bozkurt, E.; Aksoy, E.; Eski, S.; Özmen, B.; Alan, H.; Metin, Y. Surface deformations of 24 January 2020 Sivrice (Elazığ)–Doğanyol (Malatya) earthquake (Mw = 6.8) along the Pütürge segment of the East Anatolian Fault Zone and its comparison with Turkey’s 100-year-surface ruptures. Mediterr. Geosci. Rev.; 2020; 2, pp. 385-410. [DOI: https://dx.doi.org/10.1007/s42990-020-00037-2]

19. Zhang, Z.; Zhang, W.; Chen, X. Dynamic Rupture Simulations of the 2008 Mw 7.9 Wenchuan Earthquake by the Curved Grid Finite-Difference Method. J. Geophys. Res. Solid Earth; 2019; 124, pp. 10565-10582. [DOI: https://dx.doi.org/10.1029/2019JB018630]

20. Taufiqurrahman, T.; Gabriel, A.A.; Li, D.; Ulrich, T.; Li, B.; Carena, S.; Verdecchia, A.; Gallovič, F. Dynamics, interactions and delays of the 2019 Ridgecrest rupture sequence. Nature; 2023; 618, pp. 308-315. [DOI: https://dx.doi.org/10.1038/s41586-023-05985-x]

21. Wen, Y.; Cai, J.; He, K.; Xu, C. Dynamic Rupture of the 2021 MW 7.4 Maduo Earthquake: An Intra-Block Event Controlled by Fault Geometry. J. Geophys. Res. Solid Earth; 2024; 129, e2023JB027247. [DOI: https://dx.doi.org/10.1029/2023JB027247]

22. Jia, Z.; Jin, Z.; Marchandon, M.; Ulrich, T.; Gabriel, A.-A.; Fan, W.; Shearer, P.; Zou, X.; Rekoske, J.; Bulut, F. et al. The complex dynamics of the 2023 Kahramanmaraş, Turkey, Mw 7.8-7.7 earthquake doublet. Science; 2023; 381, pp. 985-990. [DOI: https://dx.doi.org/10.1126/science.adi0685] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37535759]

23. Wang, Z.; Zhang, W.; Taymaz, T.; He, Z.; Xu, T.; Zhang, Z. Dynamic Rupture Process of the 2023 Mw 7.8 Kahramanmaraş Earthquake (SE Türkiye): Variable Rupture Speed and Implications for Seismic Hazard. Geophys. Res. Lett.; 2023; 50, e2023GL104787. [DOI: https://dx.doi.org/10.1029/2023GL104787]

24. Infantino, M.; Mazzieri, I.; Özcebe, A.G.; Paolucci, R.; Stupazzini, M. 3D physics-based numerical simulations of ground motion in Istanbul from earthquakes along the Marmara segment of the north Anatolian fault. Bull. Seismol. Soc. Am.; 2020; 110, pp. 2559-2576. [DOI: https://dx.doi.org/10.1785/0120190235]

25. Taufiqurrahman, T.; Gabriel, A.-A.; Ulrich, T.; Valentová, L.; Gallovič, F. Broadband Dynamic Rupture Modeling with Fractal Fault Roughness, Frictional Heterogeneity, Viscoelasticity and Topography: The 2016 Mw 6.2 Amatrice, Italy Earthquake. Geophys. Res. Lett.; 2022; 49, e2022GL098872. [DOI: https://dx.doi.org/10.1029/2022GL098872]

26. Wang, W.; Li, Y.; Zhang, Z.; Xin, D.; He, Z.; Zhang, W.; Chen, X. Rapid estimation of disaster losses for the M6.8 Luding earthquake on September 5, 2022. Sci. China Earth Sci.; 2023; 66, pp. 1334-1344. [DOI: https://dx.doi.org/10.1007/s11430-022-1078-6]

27. Gallovič, F.; Zahradník, J.; Plicka, V.; Sokos, E.; Evangelidis, C.; Fountoulakis, I.; Turhan, F. Complex rupture dynamics on an immature fault during the 2020 Mw 6.8 Elazığ earthquake, Turkey. Commun. Earth Environ.; 2020; 1, 40. [DOI: https://dx.doi.org/10.1038/s43247-020-00038-x]

28. Rosen, P.A.; Gurrola, E.; Sacco, G.F.; Zebker, H. The InSAR scientific computing environment. Proceedings of the EUSAR 2012, 9th European Conference on Synthetic Aperture Radar; Nuremberg, Germany, 23–26 April 2012; pp. 730-733.

29. Jónsson, S.; Zebker, H.; Segall, P.; Amelung, F. Fault Slip Distribution of the 1999 Mw 7.1 Hector Mine, California, Earthquake, Estimated from Satellite Radar and GPS Measurements. Bull. Seismol. Soc. Am.; 2002; 92, pp. 1377-1389. [DOI: https://dx.doi.org/10.1785/0120000922]

30. Lohman, R.B.; Simons, M. Some thoughts on the use of InSAR data to constrain models of surface deformation: Noise structure and data downsampling. Geochem. Geophys. Geosyst.; 2005; 6, Q01007. [DOI: https://dx.doi.org/10.1029/2004GC000841]

31. Zhu, L.; Rivera, L.A. A note on the dynamic and static displacements from a point source in multilayered media. Geophys. J. Int.; 2002; 148, pp. 619-627. [DOI: https://dx.doi.org/10.1046/j.1365-246X.2002.01610.x]

32. Melgar, D.; Bock, Y. Kinematic earthquake source inversion and tsunami runup prediction with regional geophysical data. J. Geophys. Res. Solid Earth; 2015; 120, pp. 3324-3349. [DOI: https://dx.doi.org/10.1002/2014JB011832]

33. Zhang, Z.; Zhang, W.; Chen, X. Three-dimensional curved grid finite-difference modelling for non-planar rupture dynamics. Geophys. J. Int.; 2014; 199, pp. 860-879. [DOI: https://dx.doi.org/10.1093/gji/ggu308]

34. Zhang, W.; Zhang, Z.; Li, M.; Chen, X. GPU implementation of curved-grid finite-difference modelling for non-planar rupture dynamics. Geophys. J. Int.; 2020; 222, pp. 2121-2135. [DOI: https://dx.doi.org/10.1093/gji/ggaa290]

35. Okada, Y. Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am.; 1992; 82, pp. 1018-1040. [DOI: https://dx.doi.org/10.1785/BSSA0820021018]

36. Güvercin, S.E.; Karabulut, H.; Konca, A.Ö.; Doğan, U.; Ergintav, S. Active seismotectonics of the East Anatolian Fault. Geophys. J. Int.; 2022; 230, pp. 50-69. [DOI: https://dx.doi.org/10.1093/gji/ggac045]

37. Ida, Y. Cohesive force across the tip of a longitudinal-shear crack and Griffith’s specific surface energy. J. Geophys. Res.; 1972; 77, pp. 3796-3805. [DOI: https://dx.doi.org/10.1029/JB077i020p03796]

38. Andrews, D.J. Rupture propagation with finite stress in antiplane strain. J. Geophys. Res.; 1976; 81, pp. 3575-3582. [DOI: https://dx.doi.org/10.1029/JB081i020p03575]

39. Chen, X. A systematic and efficient method of computing normal modes for multilayered half-space. Geophys. J. Int.; 1993; 115, pp. 391-409. [DOI: https://dx.doi.org/10.1111/j.1365-246X.1993.tb01194.x]

40. Chen, X. Seismogram synthesis in multi-layered half-space Part I. Theoretical formulations. Earthq. Res. China; 1999; 13, pp. 149-174.

41. Wang, R.; Schurr, B.; Milkereit, C.; Shao, Z.; Jin, M. An Improved Automatic Scheme for Empirical Baseline Correction of Digital Strong-Motion Records. Bull. Seismol. Soc. Am.; 2011; 101, pp. 2029-2044. [DOI: https://dx.doi.org/10.1785/0120110039]

42. Melgar, D.; Bock, Y.; Sanchez, D.; Crowell, B.W. On robust and reliable automated baseline corrections for strong motion seismology. J. Geophys. Res. Solid Earth; 2013; 118, pp. 1177-1187. [DOI: https://dx.doi.org/10.1002/jgrb.50135]

43. Wang, W.; Zhang, Z.; Zhang, W.; Yu, H.; Liu, Q.; Zhang, W.; Chen, X. CGFDM3D-EQR: A Platform for Rapid Response to Earthquake Disasters in 3D Complex Media. Seismol. Res. Lett.; 2022; 93, pp. 2320-2334. [DOI: https://dx.doi.org/10.1785/0220210172]

44. Worden, C.; Gerstenberger, M.; Rhoades, D.; Wald, D. Probabilistic relationships between ground-motion parameters and modified Mercalli intensity in California. Bull. Seismol. Soc. Am.; 2012; 102, pp. 204-221. [DOI: https://dx.doi.org/10.1785/0120110156]

45. Wessel, P.; Luis, J.F.; Uieda, L.; Scharroo, R.; Wobbe, F.; Smith, W.H.F.; Tian, D. The Generic Mapping Tools version 6. Geochem. Geophys. Geosyst.; 2019; 20, pp. 5556-5564. [DOI: https://dx.doi.org/10.1029/2019GC008515]