Kinematics and Dynamics of the 24 January 2020 Mw 6.7 Elazig, Turkey Earthquake

We determine rupture kinematics of the 2020 Mw 6.7 Elazig, Turkey earthquake from joint inversion of interferometric synthetic aperture radar (InSAR) measurements, regional 1 Hz Global Navigation Satellite System (GNSS), strong motion, and teleseismic waveforms, and we also use dynamic modeling to assess the faulting properties to explain the observed kinematics. Our work shows that this event predominantly ruptured unilaterally toward the SW along the East Anatolian Fault Zone at a speed as slow as 2.0 km/s for ~20 s, and three main asperities are formed with a depth ranging from 20 km to the surface, but the surface rupture seems negligible. Besides, the dynamic model reveals an initial heterogeneous stress distribution with variations up to 30 MPa, which has been probably built up during the interseismic period. While this event does not seem to promote the failure of Pazarcık seismic gap, it remains elusive to evaluate the disturbed seismic potential between Elazig and Bingol regions.


Introduction
Historical ruptures and plate kinematics (e.g., fault locking depth and interseismic loading rate) are essential to assess seismic hazards (e.g., Avouac, 2015;Dolan et al., 2007;Hubert-Ferrari et al., 2020;Satake & Atwater, 2007). Information on past earthquakes are commonly inferred from the available paleoseismic field studies. Regional deformation rates, on the other hand, are nowadays routinely estimated based on geodetic measurements, for example, interferometric synthetic aperture radar (InSAR) and Global Navigation Satellite System (GNSS) time series, which have been well developed and archived over the past decades (Blewitt et al., 2018).
As one of the major intracontinental transform faults in the Eastern Mediterranean region, the left-lateral East Anatolian Fault Zone (EAFZ) forms an~580 km plate boundary between the Arabian and Anatolian plates and has been responsible for a series of damaging historical earthquakes (Duman & Emre, 2013;Taymaz et al., 1991) (see Figure 1). Interestingly, after a seismic burst with a series of M w > 6.8 events between 1871 and 1905, the EAFZ has been relatively quiescent until 2020. Under such circumstances, many studies have attempted to assess seismic hazards along the EAFZ from various aspects. For example, Nalbant et al. (2002) calculated the stress evolution caused by the modeled Ms ≥6.6 earthquakes since 1822 and tectonic loading to determine the likely location and magnitude of future earthquakes along the fault zone. Very recently, Hubert-Ferrari et al. (2020) investigated the seismic cycle along EAFZ using paleoseismic record. However, paleoseismic records are only limited to a few places and precise locations, magnitudes, and fault geometries of past events may not be well recognized from paleoseismological investigations (Grant Ludwig, 2015).
A number of previous studies have also tried to estimate the fault locking depth and slip rate using InSAR and GNSS velocity fields (e.g., Bletery et al., 2020;Cavalié & Jónsson, 2014;Walters et al., 2014), which contribute to evaluate the strain building and identify the slip deficits. With regard to the geodetic approach, a long outstanding problem is that there exists a trade-off in the estimation of slip rates and locking depth (Savage & Burford, 1973). As a result, we can see a significant discrepancy for seismic hazard assessment along EAFZ among different studies. For instance, despite both use InSAR line of sight (LOS) velocity field, the very shallow locking depth 4.5 km estimated by Cavalié and Jónsson (2014)) differs significantly from 15 ± 5 km by Walters et al. (2014), and the identified seismic gaps also locate from the east Karlıova triple junction (Bayrak et al., 2015) to the west Pazarcık and Amanos segments (Duman & Emre, 2013).
The 24 January 2020 M w 6.7 Elazig earthquake, which struck the Puturge segment in Elazig at 17:55 UTC and caused 41 deaths (Cetin et al., 2020), marked the end of seismic calm period since 1971 M6.8 Bingol earthquake. The M w 6.7 Elazig event could provide us valuable insights in reinterpreting seismic potential along EAFZ and surrounding regions. In this study, we investigate the rupture kinematics of the 24 January 2020 M w 6.7 from joint inversion of InSAR interferograms, high-rate (1 Hz) GNSS, strong motion, and broadband teleseismic P waveforms. We also perform dynamic rupture simulations based on the stress drop inferred from the kinematic model to ensure a physically plausible rupture process. Furthermore, we discuss and analyze how our results shed light on locking depth estimation, historical rupture identification, and seismic potential elevation along EAFZ.

Kinematic Source Modeling
We use 961 downsampled InSAR measurements from an ascending and a descending pair, 1 Hz displacement waveforms at 6 GNSS stations, velocity waveforms at 10 strong motion sensors, and P wave records at 22 broadband seismic stations to invert the rupture kinematics. GNSS, strong station distributions, and InSAR interferograms are shown in Figure 1 (note that the GNSS station ERGN is colocated with strong motion 2104); more information about the data processing and preparation is supplemented in Text S1in the supporting information (Chen & Zebker, 2002;Dziewonski & Anderson, 1981;Farr et al., 2007;Jónsson et al., 2002;Rosen et al., 2012). With respect to the fault geometry, we define a fault plane of 60 × 22.5 km 2 that is sufficiently large enough to cover the deformation zone as revealed by the InSAR data, test different hypocenter locations, origin time from several institutions (see Table S1), run interations of strike, dip angles with initial values which are based on global centroid moment tensor (GCMT) solutions (Ekström et al., 2012), and vary the rupture speed from 1.6 to 3.4 km/s. We adopt the frequency-wavenumber integration method (Zhu & Rivera, 2002) to compute Green's functions for near-field observations using 1-D layered velocity model (see Table S2) CRUST 2.0 (https://igppweb.ucsd. edu/~gabi/crust2.html), and a thorough description about the inversion procedure is provided in Text S2 (Chen et al., 2018;Dziewonski & Anderson, 1981;Hartzell & Heaton, 1983;Konca et al., 2010).
The preferred strike and dip angles after iterations are 240°and 75°, respectively. The fault plane follows approximately the fault trace of EAFZ (see Figure 2). We try to refine the fault geometry in order to fit better the curved fault trace but find a sharp drop of the InSAR fits, and we believe that our geometry should be an adequate, though idealized, representation of fault structure at the scale of the seismogenic zone. Furthermore, the rupture speed that fits all the data sets best is as low as 2.0 km/s (see Figure S1), consistent with what is estimated based on back projection using a teleseismic array (Pousse-Beltran et al., 2020).

Earth and Space Science
The low rupture speed is usually considered as an indicator of fault immaturity, which reconfirms the EAFZ as a young strike-slip fault (e.g., Bulut et al., 2012). Our favorable slip model, the corresponding moment release rate, and rupture evolution with snapshots at 2 s interval are shown in Figure 2; slip uncertainty estimation through jackknife test for each patch is appended in Figure S2; and data fits for InSAR, high-rate GNSS, strong motion, and teleseismic waveforms are supplemented in Figures S3-S6. The peak slip is slightly above 2 m at a depth of 4.5 km, and the total seismic moment (1.62 × 10 19 Nm, equivalent to a magnitude of M w 6.7) is almost released in 20 s. Along the strike, the earthquake propagates~35 km toward the southwest and~10 km toward the northeast, mainly unilaterally as clearly demonstrated. Along with the dip, this event ruptured a depth ranging from~20 km to the surface, and the derived~10 km centroid depth is close to the GCMT solution (12 km). Only small shallow slip (generally less than 40 cm) is revealed by this model, which are actually consistent with the absence of primary surface rupturing observed in preliminary field investigations (Cetin et al., 2020).
Four main asperities might be distinguished (denoted as Asperities A, B, C, and D in Figure 2). Asperity A locates NE to the hypocenter, while the other three are on the SW side. Rupture evolution snapshot shows that Asperity A is the first to rupture. However, the rupture front reached its NE limit after just~7 s and since then the propagation is dominantly to the SW. The moment release rate peaked at~10 s after the origin time with a maximum of slip rate of~1.5 m/s, the sharp increase implies a pulse-like rupture as also seen in the high-rate GNSS waveforms, and Asperities B and C ruptured during this period. After that, we find a significant drop in the moment release rate over the following~3 s, but a rupture zone that extended further along the strike at a depth between Asperities B and C accelerated the moment release rate again.
The data fits are fairly well for InSAR and P waveforms with variance reductions (VR) as 70.9% and 52.6%, respectively. Note that the size of this event is on the small side for typical teleseismic finite source inversion; the good fits of P waveform imply that the~20 s duration of the moment release rate model should be reasonably well resolved. By contrast, VR of high-rate GNSS, strong motion is just 40.2% and 33.4%. Particularly, due to a directivity pulse, the peak displacement waveforms at high-rate GNSS stations MALY and ADY1 are severely underestimated, which is also reported in Melgar et al. (2020). Further improvements of the strong motion fits would probably require a more precise 3-D local velocity structure.

Dynamic Source Modeling
To investigate the possible causes of the source complexity of the 2020 Turkey earthquake, we run the dynamic rupture simulation of this earthquake. We use a curved grid finite difference method (Zhang et al., 2014) assuming a slip-weakening law (Ida, 1972). The fault geometry is the same as the planar one used in the kinematic inversion, and the other dynamic rupture parameters used in our dynamic rupture simulations are listed in Table S3. A constant critical slip-weakening distance Dc of 0.3 m is assumed at depth larger than 2.5 km. The Dc value decreases linearly from 0.3 m at a depth of 2.5 km to 1.5 m at the free surface so as to match the observed slip deficit near the free surface.
The initial stress is an important factor governing the rupture behavior. To reproduce the complex rupture behaviors of the Elazig earthquake, we implement heterogeneous initial stresses on the simplified planar fault. Considering that the strike slipping dominates the Elazig earthquake, we only evaluate the initial stress along the strike direction, leaving zero values for the initial stress along the dip direction. Moreover, to make our dynamic rupture modeling more reliable, we approximate the initial shear stress on the fault plane according to the slip pattern from the inversion work by (1) calculating the static stress drop Δτ using the slip distribution of the kinematic inversion results (Figure 3a), (2) setting a constant stress drop τ at the place where having a positive value of s 0 , and (3) approximating the initial stress τ 0 that was used in the dynamic rupture modeling by τ 0 = τ + μ d · σ n , in which the dynamic friction coefficient μ d and the normal stress σ n are listed in Table S3. After trial and error tests, the τ value is set to 7 MPa to generate an M w 6.8 earthquake. Moreover, to trigger the spontaneous dynamic rupture, a slightly higher stress (0.1%) than fault strength is set at the hypocenter, which is represented by the circle-shaped nucleation patch with a radius of 1 km. Finally, the initial stress on the fault plane is shown in Figure 3a. Figure 3b illustrates the slip distributions on the fault plane from the dynamic rupture simulation, and Figure S6 shows associated synthetic waveforms at GNSS and strong motion stations. Generally, the slip pattern of the dynamic model (Figure 3b) is close to that derived from the kinematic inversion results. The large 10.1029/2020EA001452

Earth and Space Science
CHEN ET AL. slip in Asperities A, B, and D can be observed in the dynamic model ( Figure 3b). However, due to the high stress within the nucleation patch as it is a mandatory condition to trigger the rupture in the dynamic modeling, the slip near the hypocenter is large comparing to that from the kinematic inversion. Moreover, the visible patch of the coseismic slip of the inversion work at the bottom (Asperity C in Figure 2a) of the fault plane is not shown in our dynamic model (Figure 3b). Heterogeneous slip pattern indicates the complexities in the rupture process, which agrees with the kinematic inversion results (Figure 2). Waveform fits are acceptable and comparable to that of the kinematic slip model.

Discussions
We compare our kinematic model (Figure 2c) with those produced by USGS based on teleseismic waveforms (https://earthquake.usgs.gov/earth quakes/eventpage/us60007ewc/finite-fault), Cheloni and Akinci (2020) and Pousse-Beltran et al. (2020) based on InSAR observations, and Melgar et al. (2020) based on joint InSAR and high-rate GNSS records. Furthermore, we supplement rupture models from each individual data set (see Figure S7) to check their constrains on slip features. Overall, all the published slip models agree with others quite well on the shallow slip deficit and dominating SW rupture propagation. Despite that teleseismic-only inversion from USGS shows a very similar bimodal moment rate function as our model, the location of USGS model apparently contradicts the InSAR observations, which is caused by the improper epicenter selection. Pousse-Beltran et al. (2020) use two colinear fault planes with a step over to characterize the rupture; we find that while the two fault planes do increase, the InSAR fits a bit but at the expense of decreasing the waveform fits. In fact, Cheloni and Akinci's (2020) fault geometry inversion using InSAR observations also favors one single fault plane. Besides, both the models of Cheloni and Akinci (2020) and Pousse-Beltran et al. (2020) look a bit oversmoothing without evident asperities as in our model due to the lack of incorporating waveforms for inversion. Our model differs from Melgar et al. (2020) mainly in the moment release rate and Asperity D. There does not exist sharp down and up in moment release rate after~10 s in Melgar et al. (2020), and interestingly, their maximum moment release rate is lower than our model in spite of a larger peak slip. Asperity D is constrained by both strong motion and P waveforms (see Figure S7) that are not used by Melgar et al. (2020), and we find that removing Asperity D will underestimate the strong motion and P waveforms (see Figure S8). However, ignoring Asperity C almost has no influence on the data fits; we speculate that Asperity C could be an artifact induced by high-rate GNSS waveforms rather than a true feature, as also indicated by the dynamic model.
As mentioned in section 1, previous InSAR studies (Cavalié & Jónsson, 2014;Walters et al., 2014) show an apparent inconsistency about the fault-locking depth estimation. The scattered slip along dip in our study implies that this inconsistency, instead of being a result from different inversion approaches, may represent significant spatial variations in locking depth along EAFZ, which has been identified along NAFZ as caused by heterogeneous friction properties along the fault (Kaneko et al., 2013). Furthermore, the shallow slip deficit indicates possibly near-surface creeping, as also confirmed by a recent survey by Dogan et al. (2020). However, all of the InSAR tracks used by Cavalié and Jónsson (2014) and Walters et al. (2014) do not show any sharp displacement discontinuities, indicating that the fault creep should not have reached the surface during the InSAR observation period and is thus possibly a transitory phenomenon. Taking all the above into account, we believe that besides imparted by past events (Nalbant et al., 2002), heterogeneous stress distribution is being built up constantly during the interseismic period on locked patches along EAFZ as revealed by a recent study by Bletery et al. (2020), which can then fail separately or collectively during earthquakes (Kaneko et al., 2010), and this is probably the main reason for the irregularity of earthquake intervals along EAFZ.

Earth and Space Science
Areas of particular seismic risk (see Figure 4) before the 2020 M6.8 event have been identified by several previous studies (e.g., Bayrak et al., 2015;Duman & Emre, 2013;Nalbant et al., 2002). While both Duman and Emre (2013) and Nalbant et al. (2002) suggested Pazarcık segment to be one of the most important seismic gaps along the EAFZ and expected it to yield a large event in the future; they have different views as to the exact rupture segments of the 1874 M7.1 and 1875 M6.7 earthquakes. Nalbant et al. (2002) inferred the two events to have occurred on the Puturge segment (see the approximated regions in Figure 4). However, following the paleoseismological exploratory trenching study of Cetin et al. (2003), Duman and Emre (2013) proposed that these two earthquakes should have occurred on the Palu-Lake Hazar segment (see Figure 4), just to the northeast of the Puturge segment. Given the~9 mm/yr slip rate as adopted by Nalbant et al. (2002), the~1.3 m accumulated slip deficit (ignoring possible creeping effects) fails to match our kinematic model. In this regard, the seismic potential between Elazig and Bingo (EB) regions estimated by Nalbant et al. (2002) may need to be revised.
Last but not least, historical catalog (Ambraseys & Jackson, 1998) shows seismic episodes with bursts every few hundred years separated by quiescence periods along the strike-slip faults in the Eastern Mediterranean region, and the 2020 M w 6.7 event raises an open question of whether this event and the earlier 2010 M6.1 Kovancılar earthquake will mark the beginning of an earthquake cluster, similar to the nineteenth century sequence. Essentially, the seismic cluster results from the fact that a large earthquake produces stress perturbations on the surrounding faults, which are likely to rupture soon after. The 2010 M6.1 Kovancılar earthquake was considered to have increased <0.1 bar Coulomb stress levels on the fault of the 2020 event (Akkar et al., 2011), quite unlikely to be a causative factor (King et al., 1994). We then calculate the Coulomb stress failure Δσ f associated with the 2020 event: where Δτ and Δσ n are the changes in shear and normal stress on the receiving fault plane and μ′ is the apparent friction coefficient, here we set μ′ to be 0.4 as adopted in the study of (Nalbant et al., 2002). The obtained Δσ f (see it in inset map in Figure 4) is negligible on the Pazarcık segment, and thus, the 2020 M w 6.7 earth-  1875 M6.8, 1893 M7.1, 1905 M6.8 events are also outlined, red and blue dashed rectangles denote rupture areas adopted from Duman and Emre (2013) and Nalbant et al. (2002), respectively. Inset map shows Coulomb Stress Failure induced by the 2020 M w 6.7 earthquake, note color bar is saturated. quake will probably not advance ruptures on this seismic gap. The seismic moments accumulated on the EB region, however, if were only partly released during the 1874 M7.1 event, could potentially be accelerated to a damaging event in the context of increasing stress levels from 1971 M6.8, 2010 M6.1, and 2020 M6.8 earthquakes. Seismic risk estimation is sensitive to the previous history of large earthquakes in the region, and we need more detailed investigations to constrain the exact rupture geometries of previous earthquakes on these segments.

Conclusions
We have retrieved the rupture kinematics of the 2020 M w 6.7 Elazig, Turkey earthquake using both geodetic and seismic measurements, and our preferred slip model is also proven to be physically defensible as demonstrated by dynamic rupture simulation. This earthquake propagates mainly unilaterally toward SW along EAFZ at a slow speed 2.0 km/s; the~20 s duration might result in three asperities: two of them at between 2 and 10 km depths and a deeper slip region extending down to 20 km depth, but there is an apparent absence of surface rupturing. The slip distribution indicates a significant variation of locking depth and heterogeneous stress building during the interseimsic period, which accounts for the irregularity of historical earthquake intervals. Our work also challenges the view that the 1874 M7.1 earthquake occurred on the Puturge segment. We conclude that the 2020 M w 6.7 earthquake will not advance ruptures of the Pazarcık seismic gap.

Data Availability Statement
Original Sentinel-1 data were acquired and processed by European Space Agency Copernicus program (https://scihub.copernicus.eu/) and retrieved from Alaska Satellite Facility Distributed Active Archive Center; these data are stored online (https://zenodo.org/record/4031927). High-rate GNSS displacement waveforms were provided by Melgar et al. (2020); strong motion records were provided by Disaster And Emergency Management Presidency of Turkey (AFAD) (https://tadas.afad.gov.tr/event-detail/8071). Teleseimic waveforms were obtained through the Data Management Center of the Incorporated Research Institutions for Seismology (https://ds.iris.edu/wilber3/find_stations/11175173). All of the links are last accessed 2020 September. Most of the figures in this paper were prepared using Generic Mapping Tools (Wessel et al., 2013).