Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.
Communications Earth & Environment volume 5, Article number: 431 (2024 ) Cite this article fault recording
Rupture velocity affects the shaking damage potential of an earthquake, but it is hard to estimate. Here we track the rupture front of the first 25 s of the 2023 Kahramanmaraş earthquake using the P-wave polarization vector vs time as measured by the dense three-component accelerometer network surrounding the slipping fault. We find a coherent increase of the rupture velocity along portions of the fault zone where high coseismic slip has been previously mapped. This is detected in a domain of the Pazarcık segment of the East Anatolian fault at 40 km NE of the mainshock epicenter where the rupture becomes super-shear. This domain spatially matches the part of the Pazarcık segment where most micro-earthquakes nucleate in the inter-seismic phase and edges at north a nearly aseismic domain. We postulate the existence of a velocity-weakening subsegment on the Pazarcık fault sandwiched between velocity-strengthening domains and/or mechanic barriers to the fault slipping.
On February 6th, 2023 at 01:17:35 UTC the Mw 7.8 Kahramanmaraş earthquake struck a region of Turkey close to the border with Syria causing about 60,000 fatalities. The seismotectonic setting in the region is ruled by a complex fault network that accommodates the stress built up from the relative motion at the triple junction of Anatolian, Arabian, and African plates1,2,3,4. The Kahramanmaraş earthquake nucleated on a so far unknown splay fault of the East Anatolian fault (EAF) and ruptured a segment of about 350 km through an asynchronous bilateral rupture that caused a maximum fault slip of about 8 m concentrated on the eastern branch of the fracture5,6. The EAF is a transform fault that separates the Arabian and Anatolian plates and accommodates a relative motion between the plates of up to about 10 cm per year7. About 9 h later a second large magnitude shock (Mw 7.6) occurred northeast of the first mainshock, rupturing a 200 km long, EW-oriented, secondary branch of the EAF system8.
The overall space and time evolution of the doublet earthquake ruptures has been studied by several authors using different approaches, such as inversion and modeling of the strong motion waveforms, GNSS and SAR static displacement, and aftershock relocation5,6,8,9,10,11,12,13,14. The majority of models agree that the Mw 7.8 mainshock fracture initiated along a minor branch of the EAF (on the Nurdaği-Pazarcık Fault (NPF)—within the so-called Narli Fault Zone (NFZ)) and propagated during the first 25 s mainly north-eastwards along the EAF, while at later times it propagated both in NE and SW directions. However, all models are consistent in showing that most of the slip occurred in the NE part of the fractured fault. Despite this general agreement of the overall fracture propagation, the details of the slip dynamics and kinematics depend on the fault model, inversion method, and data.
A critical aspect of the earthquake rupture process is the spatiotemporal evolution of the fracture front and its spatially variable speed15,16. The seismology community has put much attention on this aspect because the potential damage caused by an earthquake may be strictly related to changes in the rupture velocity due to the possible generation of Mach wave fronts17,18.
The rupture front migration along the fault surface is modeled as a crack propagation19,20. In this theory, an equation relates the rupture velocity with the dynamic stress drop and the slip velocity21,22,23. In24,25 it was shown that abrupt changes of rupture velocity and/or slip velocity during the rupture propagation can affect the amplitude of high-frequency radiated body wave field thus potentially increasing the damaging effect of earthquakes as related to source-controlled strong ground motion. While for large earthquakes the rupture speed is generally estimated with a relatively high uncertainty from the kinematic inversion modeling of seismic and GPS records, for the vast majority of studies on microearthquake source parameter determination the value of the rupture velocity is fixed a priori based on given rupture mode developments (circular, elliptical, rectangular, ...) and theoretical relations between source parameters (rupture radius/length) and measured quantities on seismograms (body wave duration or corner frequency). In turn, accurate estimates of the mean rupture speed during the earthquake faulting allow for a reliable determination of the static stress drop and related fault slip, which are generally estimated with large uncertainty due to the poor determination of the rupture length and area. In that sense, a reliable and independent estimate of the rupture velocity is very important for a good characterization of the seismic source.
Until the sixties of the last century, theoretical models have excluded that during an earthquake the fracture front can propagate with a speed larger than the Rayleigh wave velocity in materials in case of mode I (tensile) crack propagation. The same limiting velocity has been assumed for mode II (inplane) and III (antiplane) cracks that are best appropriate to describe the shear faulting mechanism of crustal earthquakes26,27. Subsequent theoretical works28,29 and experimental evidence30 confirmed that, although in most materials rupture propagates at sub-shear velocities, it may happen that it overcomes this limit and becomes super-shear, although this condition makes the crack propagation unstable31.
In more recent years, observations of “super-shear” seismic ruptures both from earthquakes and laboratory experiments have increased32,33,34,35,36,37, and the models of fracture propagation have accounted for a combination of different fracture modes and frictional conditions19,29,38,39,40, that theoretically allow rupture speeds larger than the shear velocity and upper limited by the P-wave velocity (subsonic waves).
Here we focus on the rupture velocity profile during the first phase of the earthquake process, i.e., during the first 25 s (T) of the Kahramanmaraş earthquake fracture, following the approach proposed in41. We first calculate the polarization vectors from the recordings of the three-component accelerometers surrounding the slipping area (Fig. 1), and for each station, the polarization vector is computed in a time window sliding over the signals, i.e., running over time. At each time, the position of the radiating point is determined by averaging the solutions over the stations. In this way, we track the rupture front from the origin time (OT) to T. From the position of the rupture front over time, we infer the instant rupture velocity (see the section “Methods” for a detailed description of the numerical workflow). The time history of the rupture propagation is finally critically discussed in relation to the spatiotemporal properties of the co-seismic slip, of the density of the aftershocks, and of the known seismotectonic segmentation of the EAF.
a Map of the strong-motion stations (black triangles), the epicenter of the Mw 7.8, February 6, 2023, Kahramanmaraş earthquake (star), and projection of the slipped fault on the surface (black line)6. b Graphical sketch displaying main fault segments of the EAF activated during the mainshock.
As pointed out in the original paper41, this approach provides stable and reliable solutions for unilateral rupture propagation, for which waves emitted by a source moving along a specific near-horizontal source rupture direction dominate the radiated wavefield. For this reason, here we will track the rupture front and consequently determine the rupture speed from the OT until a characteristic time T. The value of T was chosen critically revising the current knowledge of the kinematics of the Kahramanmaraş earthquake and based on the evidence arising from the kinematic modeling of the high-frequency, large amplitude arrivals vs epicentral distance at strong motion sites located along the EAF (see “Methods”). Most of the studies so far published show evidence for asynchronous bi-lateral rupture propagation, e.g., the nucleation along the two opposite branches of the bi-lateral fault rupture is not simultaneous, with a later activation of the southern segment relative to the northern one. The time delay at which the southern rupture nucleates is variable between 30 s and 50 s in all asynchronous bilateral models5,6,8,10,11,12,14,42, except for9,13,43, which determined a <15-s delay, mostly associated with the rupture of the Nurdaği-Pazarcık segment. In the latter cases, it cannot be excluded that low-frequency, broad-band teleseismic and geodetic GNSS modeling might not have the time and spatial resolution to determine small-scale source details such as the about 10 s delay of the asynchronous starting of the rupture along northern and southern branches of the EAF. Given these literature observations and modeling and our present analysis of high-frequency large amplitude arrivals, in our work, we assumed the hypothesis for an asynchronous bilateral rupture propagation along the EAF and therefore a dominant unilateral rupture at its initial stage (first 25 s).
Figure 2 shows the location of the radiating points over time (filled-squares). The square color scales with the time from OT. In the same figure, the small colored dots display the aftershocks during the first 11 days after the mainshock, the black dots show the points sampling the fault, while the red numbers on the right mark the along-strike distance from the epicenter of some radiating points. The size of the position error bars of the radiating points accounts for the uncertainty produced by the 10-s time window used to infer the polarization vectors.
Location of the radiating point of the fault at every second (squares) with color-codes radiation time. Horizontal error bars show the uncertainty associated to the position of the radiating points and are mostly an effect of the 10-s time window adopted to compute the polarization vector, as described in the “Methods” section. The error bars displaying this uncertainty are rotated from along the fault for better readability. Black dots mark the modeled fault planes (as reported in ref. 6) projected over the surface. Blue shaded dots show the epicenters of the aftershocks of the first 10 days (size is scaled by magnitude and light blue points mean shallower events). Black crosses mark the edges of the Gölbaşi basin, while violet crosses display the along-strike limits of the segment of the Pazarcık fault capable of generating seismic sequences and earthquakes with a magnitude larger than 2.54. The numbers in red show the distances in km measured along the strike direction (sketched as a red arrow) from the epicenter of some radiating points. The star marks the epicenter position. The red circle close to the epicenter shows the position of the minimum of the error function (measured as the difference between theoretical azimuths and measured polarization azimuth averaged over the stations) fixing potential radiating points on a planar grid with 0.1°-step both on latitude and longitude. This point is plotted as a black diamond in the Supplementary Fig. 13.
The solutions confirm a unilateral north-eastward migration of the radiating points over time. Moreover, we remark that the radiating point located at time zero (nucleation of the fracture) well matches (within the uncertainty) the position of the epicenter, as one would expect.
In Fig. 3 the instant rupture velocities are plotted as a function of time from the rupture onset together with the source time function (STF—red line) and the slip rate (as determined in ref. 6—dashed black line) averaged over the longitude range 37°–37.5°, which broadly encloses the fault segment activated during the first 25 s of the rupture. While STF averages the effects of the slip along the entire rupture process over the whole slipping fault, the dashed black line should better represent the effects of the slip localized in the segment of the slipping fault activated until T. The uncertainty on the rupture velocity has been experimentally estimated by applying our approach to synthetic signals produced by a rectangular fault rupturing unilaterally under controlled conditions (see “Methods” for details). We remind that the time of each solution marks the center of the 10 s separating the positions of the two radiating points used for the calculation of the instant velocity.
Each point (square) represents the rupture velocity estimated from the location of the radiating points 10 s apart and is located on the x scale at the center of the 10 s. So, the first point displays the rupture propagation velocity during the first 10 s of the earthquake. Y-uncertainties have been fixed after applying our approach to synthetic signals and computing the oscillations of the solutions around the expected mean value, as described in the “Methods” and in Supplementary Note 1. The red line shows the STF, while the black dotted line displays the mean slip rate averaged over the longitude range of 37°–37.5°, which broadly encloses the fault segment activated during the first 25 s of the rupture. Both STF and slip rates are shown after6. A good agreement between the time history of the rupture velocity and of the cumulative slip and its rate is clearly visible.
We note an overall positive correlation between rupture velocity and STF/local slip rate. Moreover, two relative maxima of the velocity can be identified: one at the very beginning of the fracturing (during the first 10 s), and another at about 18 s (±5 s) from rupture onset, where the rupture velocity ranges in 3.4–3.7 km/s. Since the S-wave velocity in the depth range where the fault slipped is around 3.6 km/s44, this implies that the second maximum could indicate a possible super-shear fracturing. Consistently, other authors found super shear velocities at similar times by the inversion of the near-source acceleration waveforms and GNSS static displacements, combined with evidence of surface breakages5,45.
In Fig. 4a the velocities have been plotted as a function of the distance along the fault from the epicenter together with the cumulative slip (over the whole slipping area—dashed red line) and the mean slip rate averaged over the first 25 s of the rupture (dashed black line). The ratio between the instant rupture velocity and the shear wave velocity has been here plotted for comparison. We find an overall good spatial correlation between the rupture velocity and both the cumulative slip and the local slip rate: the segments of the fault where the rupture front propagates faster spatially match the segments where the fault mostly slips and where it slips faster, and the highest rupture speeds are located in fault segments that largely slipped during the earthquake. Moreover, we note that a correlation between the slip rate and rupture front velocity during the earthquake faulting process has been predicted by numerical simulations where the fault friction was governed by a slip weakening regime46. Similarly, numerical simulations show that the fault slip peak correlates with the rupture velocity47, in line with our observations.
a Rupture velocity as a function of the distance from the epicenter (squares) normalized to the S-wave velocity in the areas at shallow depths (3.6 km/s44). Y error bars display the uncertainty of the rupture velocity (shown in Fig. 3) normalized to the S-wave velocity; x error bars show the standard deviation of the along-strike distances between the two radiating points used for the calculation of the velocity. The red line and black dotted line display, respectively, the cumulative slip and the slip rate averaged over the first 25 s of the rupture. The magenta triangles mark the subsegment of the Pazarcık fault (PAZ-SUB) where sequences and M2.5+ earthquakes nucleate, as indicated in4. The green triangles mark the conjunction of the Narli fault with the EAF. The two black triangles show the position of the Gölbaşi basin. b Location of the aftershocks during the first 11 days after the mainshock; c as (b), but for the earthquakes occurring in March–May 2023. The red line shows the density of the seismicity along the strike. A minimum in the seismicity density appears at distances from the epicenter broadly ranging between 45 km and 70 km.
The results displayed in Fig. 4a suggest the evidence for a segment of the EAF located at along-strike distances of about 35–50 km from the epicenter with peculiar coupling and/or frictional properties that allow high slip and high rupture velocity. Such a segment has been already indicated as originating a super-shear fracture propagation5. This domain borders at the north with a segment with a very low aftershock density zone, located at a distance range from the epicenter of about 45–75 km (Fig. 4b, c)5,6. In fact, the aftershock density of the first 110 days after the main shock in the along-strike distance range 45–75 km is 0.07 earthquakes per day and kilometer, while this value is 0.62 for the distance ranges 0–45 km from the epicenter (Fig. 4c). Moreover, this segment shows low earthquake production also in inter-earthquake periods, i.e., in the order of few micro-earthquakes (magnitude lower than 2.5) per year4.
Our analysis of the rupture tracking of the early 25 s of the Kahramanmaraş earthquake in comparison with coseismic slip and aftershock distributions leads to the following main conclusions:
The rupture speed changes during the first 25 s along the NSF and the EAF faults well match the coseismic slip variations, in that an along-strike high-velocity segment is spatially co-located with the high-slip fault segment;
The rupture propagates with a speed close to or above the shear velocity in two segments during the first 25 s: just after the nucleation on the NPF and after about 18 s from the nucleation on the EAF, where we inferred super-shear migration;
Low-velocity rupture segments spatially correlate with low-slip segments along the NSF;
Low-slip segments along the EAF in the along-strike distance range 45–75 km from the epicenter spatially correspond to the fault segments with very low aftershock and inter-event micro-earthquake production.
As shown in Figs. 2–4, the mainshock rupture reaches the connection between NSF and the Pazarcık segment of the EAF around 10 s after the nucleation, when it has traveled about 30 km from the epicenter. This implies that the segment where the fracture propagates with super-shear velocity is located on the Pazarcık fault4,48. In general, we observe a correlation between fault slip and rupture velocity (and a super-shear fracture velocity) broadly in the time interval 15–25 s from the nucleation.
Geodetic measurements show that the Pazarcık segment of EAF slips at a rate of 4–9 mm per year4,48. Moreover, it is separated north from the Erkenek segment by the Gölbaşi Releasing Bend (GRB) that corresponds to the Gölbaşi basin, which constitutes a discontinuity and a step-over of the EAF4,48,49. Only a subsegment of the Pazarcık fault (PAZ-SUB) was able to generate sequences and M2.5+ earthquakes in the period 2007–2019, which include most of the earthquakes occurring along the Pazarcık fault4. The edges of both fault segments corresponding to GRB and PAZ-SUB are shown in Figs. 2 and 4, by black and magenta markers, respectively.
We find a maximum in the rupture velocity that spatially matches the position of PAZ-SUB. The co-seismic slip reaches a maximum in the same subsegment close to its western edge. Moving eastwards from the western edge of PAZ-SUB, that is from the fault segment showing high coseismic slip, a sharp decrease of the inter- and aftershock production appears4,6,48 (Fig. 4a). Moving farther eastwards, GRB marks the structural separation between the Pazarcık and Erkenek fault and spatially corresponds to an increase of the aftershock production (Fig. 4).
All these elements are consistent with the hypothesis of a multiple segmentation of the Pazarcık fault with heterogeneous capability to generate seismic radiation and edged at its eastern limit by GRB, which can act as a barrier for slip propagation. This hypothesis has also confirmation from historical seismicity since the coseismic slip areas of the M7.0 1795 and M7.1 1893 earthquakes ruptured separately the Pazarcık and Erkenek segments of the EAF, respectively4,50. This “barrier” could not stop the fracture propagation towards the NE in the case of the 2023 Kahramanmaraş earthquake, which ruptured additional already identified segments of the EAF, as occurred also for other large earthquakes capable of skipping low-strength “barrier” zones51.
Summarizing, we observed along the PAZ-SUB: 1. a correlation between co-seismic slip and rupture velocity; 2. that sequences and M2.5+ earthquakes occurred, as from the instrumental catalog 2007–2019; 3. that the mainshock fracture propagates with super-shear velocity. All these observations suggest that this subsegment could be controlled by a velocity-weakening frictional behavior. Differently, on the segment of Pazarcık fault north of PAZ-SUB low coseismic slip and earthquake production are observed, suggesting that this segment can behave under a different frictional regime. In this scheme, PAZ-SUB would broadly represent a transition between different fault-slipping properties. In that sense, interestingly the hypocentres of the early aftershocks deepen when crossing the PAZ-SUB6 (Fig. 4b), compatible with a sharp along-strike change of the fault conditions.
West of PAZ-SUB along NSF we observe a moderate coseismic slip, large aftershock production, and a dominance of small rupture velocities, although an almost super-shear velocity has been measured also during the very first seconds of the rupture, as observed also by other authors5. So, a possible segmentation also of this fault cannot be excluded, but it is beyond the resolution of our analysis.
A possible interpretative scheme of the results is sketched in Fig. 5. We hypothesize a segment of the Pazarcık fault (blue-shaded, partially including PAZ-SUB) governed by a velocity-weakening frictional regime that allows high coseismic slip that in turn favors an increase of fracture speed. East of this segment (yellow-shaded segment), a velocity-strengthening frictional regime can explain the low earthquake production and coseismic slip, and it is potentially compatible with aseismic slip. PAZ-SUB could represent a transition between two frictional regimes, with an asperity generating the micro-earthquakes embedded in an otherwise aseismically slipping fault, similar to the mechanism generating repeating earthquakes52,53.
In this paper, we can infer the rupture velocity of the segments plotted in red and blue; the blue segment shows high rupture velocity and high coseismic slips, and edges at the north with a segment with low aftershock density and low co-seismic slip. We postulate the existence of a velocity-weakening segment surrounded by domains with different frictional regimes. At the northern edge of the Pazarcık fault, GRB can constitute a mechanic barrier potentially able to inhibit fracture propagation, as observed for historical earthquakes.
Along the NSF (red-shaded segment) the observed high earthquake production combined with low coseismic slip and rupture velocities could be the effect of the activation of splay faults that do prevent continuous fracture propagation, but do not limit the aftershock production. Finally, at the eastern edge of the Pazarcık fault, the basin of the GRB (gray-shaded segment) can act as a mechanic barrier for rupture propagation during moderate earthquakes and marks the separation between the Pazarcık and the Erkenek fault. GRB marks the transition between low to high coseismic slip and earthquake production when moving from the eastern edge of the Pazarcık fault to the Erkenek fault (Fig. 4).
We have tracked the rupture front migration during the first 25 s of the Kahramanmaraş earthquake of February 2023 by calculating the azimuth direction from the polarization vector using 240 three-component accelerometers surrounding the slipping fault. From the position of the rupture front every second, we have inferred the speed of fault fracturing northeastwards.
We have found that the rupture propagated along the first 30 km from the epicenter with moderate velocity in a region with high aftershock production along the NFZ, while at an along-strike distance from the epicenter of about 40 km on the Pazarcık fault segment of the EAF we have detected a super-shear rupture velocity culminating about 18 s after the nucleation. Moreover, we have observed a positive correlation between rupture velocity and slip and a sharp decrease of the inter- and aftershock earthquake production when moving eastwards from the segment with high coseismic slip and rupture velocity.
We argue that the rupture propagates at the beginning through a highly fractured splay fault zone that limits the propagation velocity. When the front reaches the EAF, it propagates instead with high speed along a segment of the Pazarcık fault supposedly ruled by a velocity-weakening frictional regime and bounded at its eastern edge by a segment dominated by a velocity-strengthening slipping regime potentially slipping aseismically. Eventually, the GRB, which separates the Pazarcık and the Erkenek faults, can act as a further factor inhibiting the fault slip; however, in the case of the Mw 7.8 Kahramanmaraş earthquake, this barrier was surpassed by the rupture propagation.
The observation of elements compatible with the existence of a subsegment able to accommodate super-shear rupture propagation velocities, as well as of a sub-segmentation within a segment of EAF can have very important implications for the seismic risk assessment and the definition of shaking maps in the region. As well, insights on possible aseismic slip over small segments of the Pazarcık fault provide a potential research target for future observations.
For this work, the acceleration time histories recorded at the strong motion instruments of the AFAD network (https://tdvms.afad.gov.tr/) have been utilized. We selected 240 three-component accelerograms sampled at a frequency of 100 Hz, from stations located at an epicentral distance ranging between 70 km and 570 km.
In the following, the different stages of the data workflow are individually described. Data are the three-component accelerometric recordings at 240 receivers (see Fig. 1) sampled at 100 Hz. The first arrival has been manually picked on all vertical traces. Signals at all stations have been synchronized on the P arrivals, so that the common reference time is the rupture onset.
Figure 6a shows the accelerograms at stations along or near the EAF line6 and is shown in Fig. 6b. The y-scale represents the time from the OT. The x-scale displays the offset calculated from the junction point along EAF at the crossing point of the NPF with the EAF fault. Negative distances are associated with stations located SW of the epicenter (hereinafter SW stations), while stations NE of the epicenter (NE stations) correspond to positive distances. We emphasize that stations on the NE side are off-fault and sparse, while on the SW side, the stations are well deployed along the fault trace and denser distributed (see Fig. 6b).
a Section of the accelerograms (EW component) recorded at the along-fault stations displayed in (b). x- and y-scale show, respectively the time from OT and the distance from NPF-EAF junction point. The green line marks the NE propagation of the rupture along the EAF at times around 10 s from OT, and the blue dotted line shows the propagation (with a speed close to the S-wave velocity) of the radiation mostly from the junction between NPF and EAF at similar times, the red line shows the SW rupture propagation at about 30 s from OT; b map of the along-strike stations. The star and the square show respectively the earthquake epicenter and the NPF-EAF junction point. c Sketch of the rupture dynamics: during the first about 10 s, the rupture propagates along the NPF. Between about 10 s and 30 s from the OT, the rupture propagates NE and emits high-frequency radiation mostly at times 10–15 s from OT when it reaches the junction between NPF and EAF, where moreover large co-seismic slip and high rupture speed have been estimated. At later times, the rupture starts propagating southwestwards from the junction.
Seismograms at SW stations show two main high-frequency, large amplitude wave packets at 10–15 s (wave packet 1—WP1) and 30–50 s (wave packet 2—WP2) after the OT; the blue dotted line and the red dashed line broadly fit the arrival times at SW stations, respectively, of the first and the second packet. The slope of the lines fitting WP1 and WP2 are, respectively, 3.9 km/s and 3.1 km/s. Noteworthily, the maximum amplitude of WP1 strongly attenuates moving from near-epicenter stations to farther stations; differently, WP2 shows nearly uniform high-amplitude arrivals along the SW stations. Moreover, we remark that the best-fitting lines of WP1 and WP2 intercept the zero of the x-axis, respectively at about 10 s and 30 s from OT.
The different amplitude vs offset patterns and apparent velocity allow us to discriminate the nature of the WP1 and WP2 phases. WP1 is mostly formed by S-wave trains originated by a source located near the NPF-EAF junction, which justifies both the relatively higher wave speed and strong distance attenuation. Instead, WP2 behaves as high-frequency S-waves radiated from the closest patches of the fault as the rupture propagates toward SW. This explains well both the lower apparent velocity (controlled by the average sub-shear rupture speed) and the higher amplitudes vs offset.
Along the northern branch of EAF, a high-frequency, large amplitude arrival (WP3) is visible at the few available NE stations and it is fit by a green dashed line in Fig. 6a with a slope of 2.8 km/s, similar to the mean SW rupture speed - red dashed line. Although showing smaller amplitudes than the southern wave packet, WP3 still shows relatively large amplitudes that can be associated with a northeastward propagating rupture. The lack of waveform similarity between SW and NE stations is imputable to the combined effect of the location of the NE sensors (which are sparse and off-fault) and of the rapid attenuation of high frequencies. For the same reason, it is hard to unequivocally identify later arrivals at NE stations as the effect of a unique propagating radiation.
According to our travel time modeling (Fig. 6), the timing of WP1/WP3 and WP2 can be straightforwardly interpreted as the high-frequency radiation originated by an asynchronous bilateral rupture once the rupture front (nucleated on the NPF) has reached the EAF. While the OT of WP1/WP3 is consistent with the time at which the rupture reached EAF (about 10 s after OT), the time delay between WP2 and WP1/WP3 arrivals constrains the activation of the southern rupture about 15–20 s later, e.g., 25–30 s after the OT.
Figure 6c summarizes the proposed scheme. The source of radiation S1 corresponds to the nucleation of the Kahramanmaraş earthquake. After about 10 s from the nucleation, the main radiating source (S2) is located at the junction between NPF and EAF and triggers the NE rupture propagation and the S-wave radiation visible on the two sides of the fault as WP3 and WP1, respectively. About 30 s after the nucleation, the SW rupture propagation (S3) starts from the NPF-EAF junction generating the WP2 arrivals. So, the activation of the SW rupture propagation would occur about 20 s after the activation of the NE rupture and about 30 s after the OT. This scheme is in line with the interpretation of most of the literature studies that identify the SW rupture nucleation at about 30 s after the OT, while the NE rupture would start immediately once the rupture reaches the EAF (see “Introduction”).
The starting point of our approach is the calculation of the polarization vector individually for each station in the selected time window sliding over time. We use the Kanasewich algorithm based on the diagonalization of the covariance matrix54. The polarization vector is eventually identified as the eigenvector of the covariance matrix corresponding to the largest eigenvalue. The direction of this vector is identified by two angles: the dip and the azimuth angles. The former measures the angle between the polarization vector and the vertical direction, and the latter measures clockwise the angle between the polarization vector and the north direction. In this work, we only consider the azimuth angle, which implies that we will focus on the projection of the particle motion on the horizontal plane. A polarization vector is defined with an uncertainty of 180° because the algorithm can only determine the direction of oscillation but not the sense.
A rectilinearity (RL) is associated with the polarization vector; it measures the degree of linearity of the ground particle motion. A RL equal to 1 means perfectly linear oscillations, RL = 0 means that the three eigenvalues of the covariance matrix are identical and no preferential direction of oscillation can be identified. To analyze only meaningful polarization vectors, we discarded solutions with RL lower than 0.7.
Before inferring the polarization vector along the traces, the frequency filter and the time window for the computation of the Kanasewich algorithm have been fixed. Different bands in the range between 0.01 Hz and 40 Hz and different time windows and sliding overlapping fractions at each time step have been tested: for each of them, the stability of the polarization solutions for the signals at each station has been evaluated together with the capability to track the epicenter at t = 0 of the rupture. At the same time, the optimal time window and its overlapping at each time step have been selected as a trade-off between the time resolution and the stability of the polarization solutions. In this way, we selected an optimal frequency band of 0.1–10 Hz and a time window of 10 s sliding by 1 s at each time step (see Supplementary Figs. 9–12 for examples of polarization).
To mitigate potential instability effects of the solutions, the polarization azimuth measured at one receiver in a time frame is calculated as the median of three consecutive solutions corresponding thus to the polarization vector of three 10-s time windows starting within 2 s.
We consider the points of the fault activated during the earthquake as potential radiating points and thus potential source positions. In that sense, our approach needs some knowledge of the geometry of the slipping fault. For each station and each fault point, the theoretical azimuth is computed. The position of the radiating point for one station and the one-time frame is eventually fixed as the point of the fault that minimizes the difference between theoretical azimuth and observed polarization azimuth, accounting also for the 180° uncertainty of the polarization azimuths. The final source position in a one-time frame is defined as the median in longitude and latitude of the position of the radiating points over all stations in the same time frame. In this step, only stations for which the S-wave is expected to arrive after the processed time frame are considered. The expected S arrival relative to the P-wave is estimated using a P-wave velocity of 6 km/s and an S-wave velocity of 3.6 km/s44.
Once the position of the radiating fault point at one time is determined, the procedure is repeated for all time frames. In this way, a time history of the position of the radiation is calculated. For the first 15 s of the rupture, we constrain the source to be radiated from the fault segment NPF, where the earthquake nucleates; this must avoid artificial scattering of the solutions between the NPF and EAF due to very similar theoretical azimuths between the points of the two branches at most stations.
Error in the location of the source for each station has been fixed assigning an error to the observed azimuths of 2.5° as in ref. 41 and projecting this cone along an arc with a radius equal to the source-station distance. Remarkably, this azimuth uncertainty is similar to the minimum of the azimuth error function at the epicenter for a planar grid of sources not limited to the fault points (see Supplementary Fig. 13) and can thus be considered reasonable also in our case. In this way, the farther is the station from the epicenter, the larger the uncertainty associated with the azimuths. As stations at short distances from the source are discarded soon after the earthquake nucleation because they may have already recorded the S-wave, the error on source location increases with time (Supplementary Fig. 14). This is also a good argument to limit the analysis to the first phase of the rupture. The typical position error derived from the azimuth uncertainty is in the order of 15–18 km (Supplementary Fig. 14); this value is about half of the distance traveled by the rupture front in 10 s, i.e., the time window used for computing the polarization vector (assuming a typical migration velocity in the order of 3–3.5 km/s), which can be thus considered the main source of uncertainty. Based on that, we eventually fixed the error on the position of the radiating point at each time step as +/− the position error derived by the azimuth uncertainty averaged over the stations (plotted in Supplementary Fig. 14).
The center of the 10-s window is used as the reference time for each time frame. Starting from the source position at every second, the rupture speed is calculated by dividing by 10 s the along-strike distance between two solutions 10 s apart. The along-strike position of a radiating point is calculated by integrating the distance along the fault from the epicenter. The position associated with a rupture velocity estimate is the average between the along-strike positions of the two solutions (10 s apart) used for the rupture velocity calculation. The error on this position is fixed as the standard deviation of the along-strike distance between the two solutions—the x-uncertainty indeed increases for solutions with larger velocities.
We have tested the capability of this approach to track a unilateral rupture propagation generating synthetic signals produced by a 150-km long rectangular fault fracturing with a uniform velocity55. As described in Supplementary Note 1, we can correctly track the rupture front and retrieve on average the expected rupture velocity. We used this simulation to assign an uncertainty to the velocity estimate: using a conservative approach, we fixed the velocity error \({\sigma }_{{v}_{r}}\) as +/− three times the standard deviation of the rupture velocity oscillations around the mean value in the synthetic case. Based on this, it was eventually fixed \({\sigma }_{{v}_{r}}\) = 0.36 km/s.
As a reference fault line, we used the fault model presented in6. This fault geometry model has been built using the combined information from the relocated aftershocks with the double-difference technique, mapped traces of all known structures (EAF Zone and Sürgü Fault), and mapped surface ruptures from remote sensing56. Although other rupture models have been proposed in the successive literature5,8 there are no substantial changes in the mainshock rupture geometry as proposed by the original paper6 that is therefore adopted in our study.
For the location of the radiating points, we only considered the horizontal position of the points discarding their depths—our approach is based on the polarization azimuths and can thus only locate the radiation on the map, while it cannot constrain the depth of the radiation.
The points sampling the faults are displayed in Figs. 1 and 2. The mean spacing along longitude and latitude between points of the fault is about 0.006° and 0.004°, respectively, corresponding to distances in the order of 500 m.
Data from strong motion instruments have been downloaded from the website of the AFAD Institute (https://tdvms.afad.gov.tr/). A list of the used instruments can be found at https://doi.org/10.5281/zenodo.12731462. Fault position, STF, co-seismic slip model, and catalog of the aftershocks of the first 11 days after the mainshock have been extracted from the supplementary material of the paper6. A catalog of the earthquakes in the period March–May 2023 has been downloaded from the KOERI/RETMC (Kandilli Observatory and Earthquake Research Institute/Earthquake-Tsunami Monitoring Center) website (http://udim.koeri.boun.edu.tr/zeqmap/hgmmapen.asp).
The numerical codes for the analysis have been developed using Matlab (R2022b, academic license https://matlab.mathworks.com/, last access 03/2024). The same software was used to produce Figs. 1–5 of the main text and all figures of the Supplementary Material; Fig. 6 of the main text was produced using the package matplotlib (v. 3.1.3, https://matplotlib.org/) and some utilities of Obspy57 (https://docs.obspy.org/) for Python 3.8. All codes are available at https://doi.org/10.5281/zenodo.12731462.
Şengör, A. M. C., Görür, N. & Şaroğlu, F. Strike-slip faulting and related basin formation in zones of tectonic escape: Turkey as a case study1. In Strike-Slip Deformation, Basin Formation, and Sedimentation (SEPM Society for Sedimentary Geology, 1985); https://doi.org/10.2110/pec.85.37.0211.
Dewey, J. F., Hempton, M. R., Kidd, W. S. F., Saroglu, F. & Şengör, A. M. C. Shortening of continental lithosphere: the neotectonics of eastern anatolia-a young collision zone. Geol. Soc. Lond. Spec. Publ. 19, 1–36 (1986).
Lyberis, N., Yurur, T., Chorowicz, J., Kasapoglu, E. & Gundogdu, N. The east anatolian fault: an oblique collisional belt. Tectonophysics 204, 1–15 (1992).
Güvercin, SE, Karabulut, H., Konca, A. Ö., Doğan, U. & Ergintav, S. Active seismotectonics of the east anatolian fault. Geophys. J.Int. 230, 50–69 (2022).
Wang, Z. et al. 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. https://doi.org/10.22541/essoar.168614604.40356164/v2 (2023).
Melgar, D. et al. Sub- and super-shear ruptures during the 2023 Mw 7.8 and Mw 7.6 earthquake doublet in se türkiye. Seismica 2, 1–10 (2023).
Reilinger, R. et al. Gps constraints on continental deformation in the africa-arabia-eurasia continental collision zone and implications for the dynamics of plate interactions. J. Geophys. Res. Solid Earth 111, B05411 (2006).
Jia, Z. et al. The complex dynamics of the 2023 kahramanmaraş, turkey, Mw 7.8–7.7 earthquake doublet. Science 381, 985–990 (2023).
Liu, C. et al. Complex multi-fault rupture and triggering during the 2023 earthquake doublet in southeastern türkiye. Nat. Commun. 14, 5564 (2023).
Petersen, G. M. et al. The 2023 southeast türkiye seismic sequence: rupture of a complex fault network. Seismic Rec. 3, 134–143 (2023).
Okuwaki, R., Yagi, Y., Taymaz, T. & Hicks, S. P. Multi-scale rupture growth with alternating directions in a complex fault network during the 2023 south-eastern türkiye and syria earthquake doublet. Geophys. Res. Lett. 50, e2023GL103480 (2023).
Abdelmeguid, M. et al. Dynamics of episodic supershear in the 2023 m7. 8 kahramanmaraş/pazarcik earthquake, revealed by near-field records and computational modeling. Commun. Earth Environ. 4, 456 (2023).
Ren, C. et al. Supershear triggering and cascading fault ruptures of the 2023 kahramanmaraş, Türkiye, earthquake doublet. Science 383, 305–311 (2024).
Ding, X. et al. The sharp turn: backward rupture branching during the 2023 Mw 7.8 kahramanmaraş (türkiye) earthquake. Seismica https://doi.org/10.26443/seismica.v2i3.1083 (2023).
Das, S. The need to study speed. Science 317, 905–906 (2007).
Das, S. Supershear earthquake ruptures–theory, methods, laboratory experiments and fault superhighways: an update. Perspectives on European Earthquake Engineering and Seismology: Volume 2. Geotech. Geol. Earthq. Eng. 39, 1–20 (2015).
Bernard, P. & Baumont, D. Shear mach wave characterization for kinematic fault rupture models with constant supershear rupture velocity. Geophys. J. Int. 162, 431–447 (2005).
Mello, M., Bhat, H. S., Rosakis, A. J. & Kanamori, H. Identifying the unique ground motion signatures of supershear earthquakes: theory and experiments. Tectonophysics 493, 297–326 (2010).
Andrews, D. J. Dynamic growth of mixed-mode shear cracks. Bull. Seismological Soc. Am. 84, 1184–1198 (1994).
Ben-Zion, Y. Dynamic ruptures in recent models of earthquake faults. J. Mech. Phys. Solids 49, 2209–2244 (2001).
Madariaga, R. High-frequency radiation from crack (stress drop) models of earthquake faulting. Geophys. J. Int. 51, 625–651 (1977).
Madariaga, R. High frequency radiation from dynamic earthquake fault models. Ann. Geophys. 1, 17–23 (1983).
Kanamori, H. Mechanics of earthquakes. Annu. Rev. Earth Planet. Sci. 22, 207–237 (1994).
Spudich, P. & Frazer, L. N. Use of ray theory to calculate high-frequency radiation from earthquake sources having spatially variable rupture velocity and stress drop. Bull. Seismol. Soc. Am. 74, 2061–2082 (1984).
Bernard, P. & Madariaga, R. A new asymptotic method for the modeling of near-field accelerograms. Bull. Seismol. Soc. Am. 74, 539–557 (1984).
Kostrov, B. Some general problems of mechanics of brittle fracture. Arch. Mech. Stos. 6, 749–776 (1970).
Freund, L. Energy flux into the tip of an extending crack in an elastic solid. J. Elast. 2, 341–349 (1972).
Andrews, D. Rupture velocity of plane strain shear cracks. J. Geophys. Res. 81, 5679–5687 (1976).
Das, S. & Aki, K. A numerical study of two-dimensional spontaneous rupture propagation. Geophys. J. Int. 50, 643–668 (1977).
Archuleta, R. J. A faulting model for the 1979 imperial valley earthquake. J. Geophys. Res. Solid Earth 89, 4559–4585 (1984).
Freund, L. The mechanics of dynamic shear crack propagation. J. Geophys. Res. Solid Earth 84, 2199–2209 (1979).
Rosakis, A., Samudrala, O. & Coker, D. Cracks faster than the shear wave speed. Science 284, 1337–1340 (1999).
Bouchon, M. & Vallée, M. Observation of long supershear rupture during the magnitude 8.1 kunlunshan earthquake. Science 301, 824–826 (2003).
Xia, K., Rosakis, A. J., Kanamori, H. & Rice, J. R. Laboratory earthquakes along inhomogeneous faults: directionality and supershear. Science 308, 681–684 (2005).
Bouchon, M. et al. Faulting characteristics of supershear earthquakes. Tectonophysics 493, 244–253 (2010).
Chounet, A., Vallée, M., Causse, M. & Courboulex, F. Global catalog of earthquake rupture velocities shows anticorrelation between stress drop and rupture velocity. Tectonophysics 733, 148–158 (2018).
Bao, H. et al. Global frequency of oceanic and continental supershear earthquakes. Nat. Geosci. 15, 942–949 (2022).
Xu, S. et al. Strain rate effect on fault slip and rupture evolution: Insight from meter-scale rock friction experiments. Tectonophysics 733, 209–231 (2018).
Kocharyan, G. G., Budkov, A. M. & Kishkina, S. B. Effect of slip zone structure on earthquake rupture velocity. Phys. Mesomech. 25, 549–556 (2022).
Pampillón, P., Santillán, D., Mosquera, JC & Cueto-Felgueroso, L. The role of pore fluids in supershear earthquake ruptures. Sci. Rep. 13, 398 (2023).
Bayer, B., Kind, R., Hoffmann, M., Yuan, X. & Meier, T. Tracking unilateral earthquake rupture by p-wave polarization analysis. Geophys. J. Int. 188, 1141–1153 (2012).
Čejka, F., Zahradník, J., Turhan, F., Sokos, E. & Gallovič, F. Long-period directivity pulses of strong ground motion during the 2023 Mw 7.8 kahramanmaraş earthquake. Commun. Earth Environ. 4, 413 (2023).
Delouis, B., van den Ende, M. & Ampuero, J.-P. Kinematic rupture model of the 6 february 2023 Mw 7.8 türkiye earthquake from a large set of near-source strong-motion records combined with gnss offsets reveals intermittent supershear rupture. Bull. Seismol. Soc. Am. 113, 1–15 (2023).
Taymaz, T. 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 804, 228745 (2021).
Rosakis, A., Abdelmeguid, M. & Elbanna, A. Evidence of early supershear transition in the Mw 7.8 Kahramanmaraş earthquake from near-field records. Preprint at https://doi.org/10.48550/arXiv.2302.07214 (2023).
Bizzarri, A. Rupture speed and slip velocity: What can we learn from simulated earthquakes? Earth Planet. Sci. Lett. 317-318, 196–203 (2012).
Schmedes, J., Archuleta, R. J. & Lavallée, D. Correlation of earthquake source parameters inferred from dynamic rupture simulations. J. Geophys. Res. Solid Earth 115, 1–12 (2010).
Khalifa, A., Çakir, Z., Owen, L. A. & şinasi, K. Morphotectonic analysis of the east anatolian fault, turkey. Turkish J. Earth Sci. 27, 110–126 (2018).
Duman, T. Y. & Emre, Ö. The east anatolian fault: geometry, segmentation and jog characteristics. Geol. Soc. Lond. Spec. Publ. 372, 495–529 (2013).
Nalbant, S. S., McCloskey, J., Steacy, S. & Barka, A. A. Stress accumulation and increased seismic risk in eastern Turkey. Earth Planet. Sci. Lett. 195, 291–298 (2002).
Robinson, D. P., Das, S. & Watts, A. B. Earthquake rupture stalled by a subducting fracture zone. Science 312, 1203–1205 (2006).
Uchida, N. & Bürgmann, R. Repeating earthquakes. Annu. Rev. Earth Planet. Sci. 47, 305–332 (2019).
Palo, M., Picozzi, M., De Landro, G. & Zollo, A. Microseismicity clustering and mechanic properties reveal fault segmentation in southern italy. Tectonophysics 856, 229849 (2023).
Kanasewich, E. R. Time Sequence Analysis in Geophysics (University of Alberta, 1981).
Zollo, A. & Bernard, P. How does an asperity break? new elements from the waveform inversion of accelerograms for the 2319 ut, october 15, 1979, imperial valley aftershock. J. Geophys. Res. Solid Earth 96, 21549–21573 (1991).
Reitman, N.G. Fault Rupture Mapping of the 6 February 2023 Kahramanmaraş, Türkiye, Earthquake Sequence From Satellite Data (U.S. Geological Survey Data Release, 2023); https://doi.org/10.5066/P985I7U2.
Beyreuther, M. et al. Obspy: a python toolbox for seismology. Seismological Res. Lett. 81, 530–533 (2010).
We acknowledge the KOERI/RETMC institution for making available the Turkey earthquake catalog. This work has received financial support under the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.3—Project RETURN—CUP E63C22002000002.
Department of Physics “Ettore Pancini”, University of Naples Federico II, Naples, Italy
Mauro Palo & Aldo Zollo
You can also search for this author in PubMed Google Scholar
You can also search for this author in PubMed Google Scholar
M.P. co-led the study; prepared the datasets; implemented and tested the whole numerical workflow; contributed to the design of validation tests; co-participated in the design of the figures and tables; co-participated in the interpretation of the results; analyzed the synthetic signals in support of the rupture velocity estimate validation; co-elaborated the literature review and data analysis in support of the asynchronous bi-lateral rupture mechanism; and wrote the original draft. M.P. is the corresponding author. A.Z. conceived and co-led the study; contributed to the design of validation tests; co-participated in the design of the figures and tables; participated in the interpretation of the results; conceived the simulation experiment and generated the synthetic signals in support of the rupture velocity estimate validation; conceived and co-elaborated the literature review and data analysis in support of the asynchronous bi-lateral rupture mechanism; and acquired the financial support for the project leading to this publication. Both authors contributed to finalizing the manuscript.
The authors declare no competing interests.
Communications Earth & Environment thanks František Gallovič and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor: Joe Aslin. A peer review file is available.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Palo, M., Zollo, A. Small-scale segmented fault rupture along the East Anatolian fault during the 2023 Kahramanmaraş earthquake. Commun Earth Environ 5, 431 (2024). https://doi.org/10.1038/s43247-024-01597-z
DOI: https://doi.org/10.1038/s43247-024-01597-z
Anyone you share the following link with will be able to read this content:
Sorry, a shareable link is not currently available for this article.
Provided by the Springer Nature SharedIt content-sharing initiative
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.
Communications Earth & Environment (Commun Earth Environ) ISSN 2662-4435 (online)
cable fault distance locator Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.