Complex rupture dynamics on an immature fault during the 2020 Mw 6.8 Elazığ earthquake, Turkey

Physical laws governing friction on shallow faults in the Earth and spatial heterogeneity of parameters are critical to our understanding of earthquake physics and the assessment of earthquake hazards. Here we use a laboratory-derived fault-friction law and high-quality strong-motion seismic recordings of the 2020 Elazığ earthquake, Turkey, to reveal the complex rupture dynamics. We discover an initial Mw 5.8 rupture stage and explain how cascading behavior of the event, involving at least three episodes, each of M > 6, caused it to evolve into a large earthquake, contrarily to other M5+ events on this part of the East Anatolian Fault. Although the dynamic stress transfer during the rupture did not overcome the strength of the uppermost ~5 kilometers, surface ruptures during future earthquakes cannot be ruled out. We foresee that future, routine dynamic inversions will improve understanding of earthquake rupture parameters, an essential component of modern, physics-based earthquake hazard assessment. The 2020 Mw 6.8 Elazığ earthquake displayed cascading rupture behavior along an immature fault with most of its energy spent fracturing the wall rock rather than displacing the surface, according to a detailed dynamic source inversion approach.

W hile kinematic modeling of earthquakes for their space-time slip history on faults has become almost routine task nowadays 1,2 , applications of the fault friction-laws in dynamic source inversions are rare and challenging [3][4][5] . Dynamic inversions, involving causative stresses and strengths have the potential to surmount standard kinematic approaches that are considered to provide strongly non-unique results 6,7 . Instead, dynamic inversion can obtain a betterconstrained picture of the rupture evolution, because rupture dynamics strongly couple the complete rupture history in terms of the energy release, and the mechanical fault features along the propagating rupture front via the employed friction law. Although the feasibility of dynamic rupture inversions was demonstrated in the 1990's 5,8,9 , there are just a handful of realdata discoveries in coseismic dynamics, with none yet from the East Anatolian Fault Zone (EAFZ).
The EAFZ is a major intra-continental left-lateral strike-slip structure between the Arabia plate and the Anatolian Block of the Eurasia plate, slipping at~1 cm/y [10][11][12][13][14] . Despite being known as the birthplace of significant historical events 15 , the first earthquake illuminating rupture details of the EAFZ occurred just recently 16,17 , on January 24, 2020 (Mw 6.8, Global Centroid Moment Tensor, GCMT), see Fig. 1. The Elazığ earthquake killed 41 people and heavily damaged over 700 buildings 18 . Coseismic surface faulting has not been observed.
The EAFZ is considered an immature fault zone, 2-5 Ma old, with 22-33 km offset 19,20 . Despite its historical seismicity and geomorphological evidence, some studies 12,19 claim that the seismic quiescence of the EAFZ during the last~100 years suggests the EAFZ is locked, accumulating elastic strain energy. Contrarily, recent studies employing InSAR, GPS and creepmeter data report creep within a zone from surface to the seismogenic depth, along the segments to the north-east of Lake Hazar [21][22][23] .
Yet, many open questions about the creep segmentation (location, rate, origin, etc.) need to be addressed in more detail 24 .
Here, based on seismic data of the Turkish networks, and using innovative methodology, we decipher stress and frictional parameters of the Elazığ earthquake fault rupture in a detail allowing new tectonic and mechanical interpretations. In particular, we explain why the recent earthquake did not produce surface faulting, and we show that moderate earthquakes of the interseismic period are undeveloped rudiments of potentially large events.

Results
Bayesian dynamic model of the fault rupture. For our recently developed Bayesian dynamic source inversion 25 , we assume a planar fault of 75 × 20 km based on the preceding Multiple-point source modeling (see Methods and Fig. 1), adjusted to conform to the aftershocks. The fault plane reaches the surface. The spontaneous rupture is governed by linear slip-weakening friction 26 ( Fig. 2a) with spatially variable left-lateral prestress and frictional parameters. We model strong-motion records at stations shown in Fig. 2c, d. The frequency range considered is 0.05-0.30 Hz in the case of the closest stations (2301, 2302, 2308, 4401, 4404), and 0.05-0.15 Hz elsewhere. Fig. 2c shows the waveform fit for our estimated maximum a posteriori (MAP) source model with variance reduction of 70%. Some of the stations were deliberately omitted in the inversion and serve as an independent check of the predictive power of the dynamic model ( Supplementary Fig. S1). The main kinematic and dynamic parameters are displayed in Fig. 3. The seismic moment is 1.7 × 10 19 Nm (Mw 6.8), the average stress drop is 12 MPa, and the mean slip is 1.5 m. Despite the considerable uncertainty of the dynamic parameters (see Supplementary Figs. S2-S5), the key characteristics of the rupture remain stable. In particular, the rupture starts from a patch at the located hypocenter H. It initially propagates up-dip and towards the north-east (Fig. 2b), creating~5 s long minor episode in the moment rate release (Fig. 2d). The rupture reaches its final extent towards the northeast close to H', which we located as a secondary hypocenter based on strong P' phases identified in the seismograms (see Methods and Fig. 1). The rupture then turns to south-west and produces the first major moment release for another~5 s. After minor weakening, the rupture continues unilaterally towards the south-west between 10 and 15 s after the origin time and releases about 4 meters of slip at the strongest asperity with the local stress drop of 40 MPa. After almost ceasing, the rupture reactivates in the third major episode breaching bilaterally a smaller asperity of similar stress drop within the last~5 s. The rupture velocity varies significantly, having a relatively small mean value of~2 km/s. We point out that none of the accepted rupture scenarios reaches the surface, in agreement with the fact that there was no co-seismic surface rupture observed.
The three main episodes revealed here entirely agree with the three-point-source models ( Fig. 1 and Methods) and the SCARDEC time function 27 . The slip distribution and mean rupture speed agree with smooth images inferred from geodetic inversions 16,17 . As we show in Supplementary Fig. S6, although the episodic character is not very visible in the rather smooth slip distribution, it is clearly expressed in the slip rate functions, exhibiting multiple peaks due to the partial reactivation after each major stress drop episode. The ceasing of the rupture is dynamically due to the relatively large values of D c and strength excess (difference between strength and prestress) along the fault, resulting in locally large fracture energy expended in the rupture process.
The main rupture parameters 28 are the slip-weighted mean stress drop 12.0 ± 0.7 MPa, mean D c 1.04 ± 0.07 m, total fracture energy 3.3 ± 0.4 × 10 15 J, and radiated energy 0.5 ± 0.1 × 10 15 J, suggesting that relative to other shallow earthquakes the radiation efficiency of~0.1-0.2 is low 28 . The mean fracture energy, also known as breakdown work density, is 10 ± 1 MJm −2 . The latter value is in agreement with inferences for past earthquakes 29 .
The birth and growth of the rupture. The earthquake records reveal two P onsets ( Fig. 1), weak and strong, and the source of the latter (H') was located on the fault about 10 km to the northeast from the hypocenter (H). The inferred model explains this preparatory part as weak nucleation lasting 4-5 s, including a fast-slow-fast transition of the rupture from H to H' (Fig. 3). We performed a simple experiment where we reduced the initial stress in the top 12 km by 10% (corresponding to an increase of the strength excess) to simulate earlier state of the fault. We obtained an Mw 5.8 earthquake, which failed to develop into a larger event ( Supplementary Fig. S7). During 1964-2020, seven events of such size occurred at the investigated fault segment, despite being considered a seismic gap resisting significant rupture in-between the strong earthquakes of the last 150 years 30 , see Supplementary Fig. S8. These seven events ranged from Mw 5.1 to 5.7 and included an Mw 5.3 of the 4th April 2019, situated almost at the epicenter of the 2020 event. This suggests that these previous events failed to rupture large fault segments, analogously to the above simple experiment with the reduced initial stress, perhaps because they occurred too early when the segment was not yet ready to go 31,32 . Contrarily, by January 2020, the whole segment was ready to fail and the initial Mw 5.8 nucleation event spontaneously evolved to the Mw 6.8 rupture. Therefore, the final magnitude of the 2020 earthquake was not determined at the time of its nucleation. Although it may not necessarily be a typical behavior, these findings might contribute to the explanation of the wide scatter in empirical correlations between the observed Pwave onsets and the earthquake magnitude [33][34][35][36] , or perhaps even indicate a lack of any causal relation. This represents an important contribution to early warning systems.
The area of the dominant slip in our model is characterized by a paucity of aftershocks ( Fig. 3 and Supplementary S2). Most of the aftershocks appeared below the major slip, mainly in the nucleation area, and at the north-east and south-west edges of the slipped area. The spatial distribution of aftershocks is in good agreement with the Coulomb stress change 37 calculated from the inverted slip model ( Supplementary Fig. S9). The rupture seems to have stopped by reaching other distinct segments of the EAFZ, namely restraining bend near the Yarpuzlu city to the south-west and releasing bend near Lake Hazar to the north-east.
The largest slip (4 m) was released at the patch with the largest initial stress (and thus the largest stress drop) of 40 MPa, see Fig. 3. The patch is also associated with large D c and . The GCMT focal mechanism is shown with gray beachball. Three major subevents of multiplepoint source inversion of broadband and strong motion records, MPS-BB-SM, are plotted with circles sized with their seismic moment and colored with their rupture time. The multiplicity of the subevents (and their nodal lines) reflects the uncertainty revealed by jackknifing. The three best-fitting subevents, shown with black-white beachballs, are labeled by their moment magnitudes. Green and blue stars refer to epicenters H and H' calculated from the first (weak) and the secondary (strong) P and P' onsets, respectively, delimiting the source nucleation. Relocated aftershocks are color-coded relative to their occurrence time. Positions of the Ms~7 historical earthquakes are shown schematically. b Velocity records filtered between 0.05-2.5 Hz for strong-motion stations 4401 and 2302 (plotted in panel a), rotated into the fault-parallel and fault-normal system. Real data (black) are compared with synthetics (red) calculated using the dynamic model of Fig. 2 (see also Supplementary Fig. S12). The waveforms illustrate the ground motion complexity, including directivity effect. The zoomed initial parts of the waveforms demonstrate the P and P' onsets due to H and H', respectively (for other stations see Supplementary Fig. S12).
strength, resulting in locally 5-times larger fracture energy than average. The rupture speed reaches supershear values of~4 km/ s when approaching the asperity, and then the rupture slows down to propagate over the asperity at subshear speed (2-3 km/ s). This asperity is mainly responsible for strong directivity pulses at stations in the south-west direction (e.g., 4401 in Fig. 1b). After breaking this strong asperity, the rupture experienced a subsequent short inhibition due to lower prestress and large fracture energy. This can also be understood as a proxy of some (unmodeled) geometrical fault complexity [38][39][40][41] , perhaps related to the bending of the Euphrates River (Fig. 2d). The rupture then proceeded into another area of a 20-MPa stress drop. This time the entered area had small fracture energy, which resumed the rupture. In most of the accepted models, the crack propagated over this patch bilaterally and possibly at supershear speed locally, but always with strong radiation expressed by the later pulses in the seismograms (Fig. 2c).    Supplementary Fig. S1. The bottom-right inset is the moment rate function of the rupture model.

Shallow slip deficit and implications for East Anatolian Fault.
No surface rupture was created by the Elazığ earthquake. Fig. 3 suggests that the stress transfer from the rupture did not overcome the strength at the subsurface~5-km fault strip, and thus the earthquake failed to reach the surface. One of the explanations of the subsurface slip deficit is the creep in Paleozoic-Mesozoic metamorphic rocks, Mesozoic serpentiniterich ophiolites, and volcanoclastics with moderate to very low strength, characterizing the EAFZ 42,43 . This would resemble vertical segmentation of the fault hosting the 2011 Mw 7.1 Van, Turkey, earthquake, at 10-20 km depths, suggested to explain the lack of a shallower rupture 44 . This slip deficit was later at least partially ameliorated by 1.5 years of postseismic deformation 45 . In our case, a persistent creep seems to be contraindicated by the weak and rapidly decaying afterslip observed from the InSAR measurements 16 . A counteracting of the slip deficit by inelastic off-fault deformation during the event 46,47 seems also unlikely because field observations 48 and dynamic rupture simulations with fault zone plasticity 49 suggest that this effect is limited to much shallower depths (< 0.1-1 km). Although a significant slip deficit can be a general property of immature faults 50 , we propose that the energy release rate of a future earthquake may, at shallow depths, surpass the abundant fracture energy or switch to brittle rheology 51 . This way, future earthquake will eventually reach the surface, conciliating the shallow slip deficit. Examples of such surface-breaking ruptures on faults considered immature are the recent earthquake at Ridgecrest, California, Mw 7.1, 2019 52 , Norcia, Central Italy, Mw 6.5, 2016 53 , and Fandoqa, southeast Iran, Mw 6.6, 1998 54 .

Discussion
Our Bayesian dynamic source inversion revealed information imprinted in the near-fault strong motion seismograms, suggesting the following picture of the Elazığ earthquake and its broader role in the fault evolution. The event occurred on the Pütürge segment of the EAFZ, situated between strong earthquakes that happened 100-150 years ago, and, since that time, the segment remained with no Mw 6+ occurrence. The fault is loaded by aseismic slip at the base of the brittle (velocity-weakening) part of the fault at~15-20 km depth, which produces earthquakes sporadically reaching magnitude Mw above 5. Most of these Mw 5+ events occurred too early and thus failed to rupture the whole Pütürge segment. In 2020, relatively weak nucleation at 15 km depth corresponding in terms of its duration and size to an Mw5.8, succeeded to advance further. After initial up-dip and north-east propagation, the rupture stopped at the Lake Hazar releasing bend and triggered off-fault aftershocks, suggesting a distributed fault network. The rupture then continued mostly unilaterally towards the south-west. It ruptured a very strong asperity of 40 MPa at a moderate speed, that was then inhibited by locally large fracture energy. The rupture resumed and broke another 20-MPa asperity at high speed due to relatively low fracture energy, prolonging the rupture duration and increasing the earthquake seismic moment. This episodic, partially bilateral rupture propagation resulted in a relatively small mean rupture velocity (~2 km/s). Up to 90% of the total available energy was expended in the fracture process, while just 10% was radiated by seismic waves, which might be a signature of the fault zone immaturity. The earthquake failed to propagate through a 5-km subsurface fault strip and did not rupture the Earth's surface. The slip deficit and weak afterslip suggest that the shallow portions of the fault remain primarily locked and could release its accumulated strain in a future earthquake.

Methods
Mainshock relocation. We manually picked the first weak (yet identifiable) Ponset from 10 strong-motion (SM) and 8 broadband (BB) records at epicentral distances <110 km (Supplementary Table S1). During picking we estimated the inaccuracy for each station independently (P-picks: 0.03-0.6 s, S-picks: 0.2-1.7 s). The S-picks were uncertain due to emergent onset of S waves from complex preceding wavegroups. We located hypocenter H by the NonLinLoc probabilistic method 55 considering only P-picks with their uncertainty as a data error. To assess the error due to structural effects, we repeated location in several velocity models. The effect was significant mainly upon the source depth. As detailed in the next paragraph (and Supplementary Fig. S10), we prefer location in model VM1 56 (Supplementary Fig. S11a). While the epicenter is well constrained at 38.360°N, 39.091°E, with +/−2 km uncertainty in the NNW-SSE direction ( Fig. 1 and Supplementary S11), the depth remained poorly constrained between 10 and 20 km. The S-minus-P travel time difference of 3.9 s at the nearest SM station (code 2308) points to the hypocenter depth of 12-14 km and origin time 17:55:10.8 (UTC).
Note that we identified large time residuals for certain SM stations (not shown in Supplementary Table S1), obviously due to GPS time error, and hence we excluded them from the relocation. Later, in the waveform modeling, two important SM stations with the large residuals were used after forward-shifting their records by the residuals (6.6 and 4.6 s at 2301 and 2302, respectively).
Recordings at many stations exhibit a second (strong) P' onset with small amplitudes between P and P' (Fig. 1 and Supplementary S12). We picked P' (Supplementary Table S1) and located respective hypocenter H' at 38.439°N 39.188°E, 0-10 km depth, and~4.5 s after the origin time (see Fig. 1 and Supplementary S11).
Relocation of the sequence. We located the whole Elazığ sequence using manual picks adopted from the DDA Catalog of AFAD, the HYPOINVERSE code 57 , and velocity model VM1. Using other models of the region, VM2 58 , VM3 59 , and VM4 60 , spurious concentrations of hypocenters occurred at certain depths (Supplementary Fig. S10). Such artifacts are quite common: For events situated closely above a sharp velocity discontinuity, observed first arrivals are interpreted by the location code as head waves; travel times are then weakly sensitive to actual depth differences, and hence hypocenters collapse to a narrow depth range. Therefore, for all modeling in this paper, we use solely VM1. Relocation with the HypoDD code 61 efficiently focused the hypocenters, clearly revealing three aftershock clusters ( Fig. 1 and Supplementary S13). The central cluster (no. 1 in Supplementary  Fig. S13), near the mainshock hypocenter, is quite tight. The north-east (no. 2) and south-west (no. 3) clusters reveal a fault dip of 70-80°to the north-west, with a lower dip angle of the third cluster.
Multiple-point source (MPS) model. The GCMT solution of the mainshock with strike/dip/rake (°) = 246/67/-9, representing a point-source model at low frequencies and large distances, corresponds to a left-lateral strike-slip motion, consistent with the local strike of EAFZ. We improved the space-time characterization of the source by employing 12 non-clipped near-regional BB stations and 7 local SM stations ( Supplementary Fig. S14), and jointly inverted them for an MPS-BB-SM model in the complementary frequency ranges of 0.01-0.05 and 0.05-0.10 Hz, respectively. Using iterative deconvolution of the established ISOLA software [62][63][64][65][66][67] we calculated full and deviatoric moment tensors (MTs). Still, the double-couple (DC) part of those models was low, DC < 50%, due to the multiplicity of subevents. Therefore, we adopted a DC-constrained MPS modeling (DC = 100%). The multiple point sources were grid-searched at several fixed depths, 5-15 km, with just a weak preference of 10 km, hereafter assumed as a source depth in the MPS modeling. Preliminary tests have revealed that the subevent positions are better resolved along the EAFZ strike than across. Therefore, a linear grid of trial sources (azimuth 246°) is a justified approximation. Our preferred model consists of three subevents, providing the total magnitude Mw = 6.7 and fitting waveforms with global variance reduction (VR) of 0.77 ( Fig. 1 and Supplementary S14). One or two (largest) subevents would produce only Mw 6.5, VR 0.34, and Mw 6.6, VR 0.59, respectively. The improvement by considering two or three subevents is significant, see Supplementary Fig. S15. Additional subevents are smaller and have unstable properties, and do not improve the fit considerably.
Although the DC-constrained MTs of the three subevents were free in the inversion, their focal mechanisms are remarkably similar. Scalar moments and the subevent timing are only weakly varying with various perturbations of the model setup, for example, with the depth and azimuth of the grid, frequency range, used stations, etc. Fig. 1 shows the 3-point MPS-BB-SM model accompanied by a jackknife test (repeatedly removing one station). The MPS models are independent of the assumed hypocenter position. Nevertheless, the source process starts between H and H' with a subevent of Mw 6.4 at~8 s after origin time. Then the process continues with a major subevent of Mw 6.5 at~12 s, located towards the south-west, near the GCMT centroid. The process ends at~17 s with a third subevent of Mw 6.3, further displaced to the south-west. The focal mechanism of the latest subevent appears to have a gentler fault dip (in agreement with aftershocks) and a thrust-faulting component. We note that the first-motion polarities of phase P support the left-lateral strike-slip faulting mechanism of initial nucleation, which is too weak to be captured by the earliest subevent of our model.

Back-projection.
We also back-project the normalized high frequency (HF, 2-8 Hz) S-waveforms from the strong-motion recordings following the previously developed methods 68,69 . The result (in Supplementary Fig. S16) shows that~5 s after the origin time, S-wave HF energy appears close to H'. Similarly to the paper 16 , the HF sources then migrate towards the south-west along the fault strike, concentrating at similar locations as the three MPS subevents (Fig. 1).
Dynamic source inversion. We apply Bayesian dynamic source inversion 25 . The method has been validated via synthetic tests, inversion of the 2016 Amatrice, Central Italy earthquake 70 , and via modeling of rupture scenarios to reconcile the stress drop ambiguity 71 . We assume a planar fault of 75 × 20 km based on the previous Multiple-point source modeling with strike 246°, dip 80°and rake 0°, adjusted by a 2 km northward shift of hypocenter (i.e., within its inaccuracy) to conform to the aftershocks. The assumed fault plane reaches the surface. We calculate waveforms following the representation theorem, i.e., by convolving Green's functions calculated in velocity model VM1 with slip rates calculated by very efficient 3D finite-difference (FD) dynamic rupture simulator FD3D_TSN 72 .
The spontaneous rupture is governed by the linear slip-weakening friction 26 ( Fig. 2a) with spatially variable left-lateral prestress τ 0 and frictional parameters. The latter include characteristic slip-weakening distance Dc and strength τ s = µ s σ n , where µ s is the static coefficient of friction, and σ n is the normal stress assumed to increase linearly with depth (Supplementary Table S2). Dynamic friction coefficient µ d is considered constant. Since the Elazığ earthquake did not rupture the surface, the total value of stress is unimportant, and the rupture is controlled by (i) initial stress τ i , i.e., the difference between prestress τ 0 and dynamic stress τ d = µ d σ n , and (ii) friction drop Δµ = µ s − µ d from the static to dynamic friction. The parameters are bilinearly interpolated from 31 × 9 regularly distributed control points to the FD grid. At the control points, the parameter values are sought using Monte Carlo Markov Chain (MCMC) approach, assuming Gaussian data errors with a constant standard deviation of the observed displacements (0.05 m) and prior uniform parameter distribution in broad ranges (Supplementary Table S2). Rupture nucleation (a process in the area with negative strength excess over the initial stress) is permitted only within a 5-km radius from the located hypocenter.
The MCMC exploration has been initiated from a relatively smooth kinematic model obtained by the recently developed parametric inversion method 73 . The initial model had a negative VR, so it is not shown here. Significant improvement of the MCMC performance was achieved by adding a constraint on the seismic moment, formally introduced as a prior Gaussian with the mean 1.8 × 10 19 Nm (Mw 6.8) and standard deviation 0.3 × 10 19 Nm. We also added a constraint on the maximum overstress in the nucleation patch (2 MPa) to limit the effect of the nucleation on the subsequent rupture propagation.
The calculations ran on an in-house farm of 10 GPUs, where a single rupture simulation took less than 1 min. With the total number of the simulations exceeding one million, the computing time was~3 weeks. Supplementary Figs. S2-S5 demonstrates the estimated uncertainty as captured by 2800 accepted models within 2% of the estimated maximum of the posterior probability density function.