Directivity analysis of the 2017 December Kerman earthquakes in Eastern Iran

Abstract

Using an empirical Green’s function (EGF) approach and data from local to regional distances we analyzed rupture propagation directivity in the three mainshocks (ML 6.0–6.1) and in six of the largest aftershocks (ML 5.0 – 5.5) of the 2017 Kerman, Iran, seismic sequence. The EGF procedure was based on data from smaller events (ML 4.0 – 4.8). Deconvolution was applied separately to P and S phases. Using the P-wave data, we calculated relative source-time functions and examined azimuthal variations in rupture duration. In the S-wave analysis, we investigated along strike rupture directivity of the mainshocks and the largest aftershocks by evaluating azimuthal variation of the amplitude spectra. Two of the mainshocks and four of the aftershocks clearly showed rupture propagation from the south-east toward the north-west. The third mainshock and one of the aftershocks suggested almost bilateral rupture propagation, and one aftershock showed rupture directivity to the southeast. It seems that the rupture propagation direction in the area is generally to the north-west and the events which have different propagation directions are located within the NW and SE ends of the faulting area. We suggest that the general rupture propagation direction in the area is steered by regional tectonic stress field regarding the faulting orientations which have been affected by stress redistribution around a restraining bend.

1 Introduction

On 1 December 2017, a magnitude (ML) 6.2 event occurred in Eastern Iran, located about 50 km north of the Kerman province (Fig. 1a). Eleven days later, on the 12th of December, two events with ML 6.2 (at 08:43 UTC) and 6.0 (at 21:41 UTC) ruptured almost the same area. Based on the Iranian Seismological Center (IRSC; http://irsc.ut.ac.ir) locations, maximum epicentral separation between the three events was about 10 km. The close spatial and temporal clustering of these evenly sized earthquakes characterizes them as a triplet (Kagan and Jackson 2005).

The triplet is located in the boundary region between the low-lying Lut-block to the east and mountainous topography to the west (Fig. 1a). The area is surrounded by major right-lateral strike-slip faults to the southeast (the Gowk fault), northwest (Kuh-Banan), and north (the Lakar-Kuh, and Nayband faults). The three events seem to be located along the southern continuation of the Lakar-Kuh fault. According to the ISC catalog from 1964 to November 2017 (Bondár and Storchak 2011), instrumental seismicity is mostly limited to the Gowk and Kuh-Banan faults with only a few recordings near the Lakar-Kuh and Nayband faults (Fig. 1a). The most recent destructive event in the area was the 2005 Mw 6.4 at Zarand which occurred on a reverse almost E-W trending fault plane (gray beach ball in Fig. 1a). The three mainshocks of the 2017 triplet sequence also exhibit a predominantly reverse faulting mechanism (Table 1) despite their proximity to the major strike-slip faults.

Fig. 1

a Seismo-tectonic map of the study area. Major active faults are in red. Faults are from Walker et al., (2010). Blue circles are events with magnitude ≥ 5, between 1964 to Nov 2017, from the ISC catalog. Focal mechanisms are from global-cmt. The gray beach ball specifies the 2005 Zarand event. Stars are the 2017 triplet events located by IRSC. Squares mark the cities. Inset: study area shown in a larger scale. b Station distribution used for the study.Blue stations are used for averaging in the SE direction from the fault and red stations in the NW direction in the S-wave spectral analysis. Inverted triangles mark the stations used for MS1 analysis, triangles used for MS2 and circles used for MS3. The black stations in each mainshock analysis were considered as perpendicular to fault strike

Triplet earthquakes are uncommon features and are of interest since they may provide valuable insights regarding the stress transfer between the adjoining faults. Rupture directivity is one of the key issues for understanding rupture propagation and can be important in assessing hazard. Along-strike rupture directivity has been often observed and estimated by several methods. Warren and Shearer (2006) determined rupture directivity from analysis of long-period P-wave spectra in a large, globally distributed database and showed that the estimated rupture direction could be used to distinguish the fault plane among the two nodal planes of the CMT solution. Cesca et al. (2011) detected directivity by inversion of amplitude spectra from regional P-wave seismograms. They obtained apparent rupture duration at each station assuming a spatial point source approximation. Kane et al. (2013) examined rupture directivity along the San-Andreas Fault by looking at azimuthal variations of the P-wave amplitude spectra. They applied an iterative correction process using a large cluster of events to eliminate the effects of wave propagation from the displacement spectra. Other studies have used empirical Green’s function (EGF) deconvolution to eliminate path effects and assess source parameters including rupture directivity of earthquakes. Viegas et al. (2010) used velocity seismograms of both P and S-waves at regional and local distances to calculate spectral ratios at individual stations. Calderoni et al. (2013) analyzed acceleration spectra of S-waves as a function of azimuth to assess directivity of the 2009 L’Aquila seismic sequence. Later, Calderoni et al. (2015) performed a systematic detection of along-strike rupture directivity of 70 earthquakes (3.0 ≤ Mw ≤ 6.1) in the 2009 L’Aquila seismic sequence. They analyzed azimuthal variation of source spectra by looking into EGF spectral deconvolution of S-waves in broadband seismograms. All these studies detected an almost unilateral rupture propagation behavior for the seismic sequences they have considered.

Table 1 list of target events and EGF events used in this study

Savidge et al. (2019) studied the 2017 triplet using regional and teleseismic data to determine epicenter and focal mechanism of the mainshocks and 15 of the largest aftershocks. They also used InSAR to deduce the fault geometries and estimated the slip distributions of the mainshocks. They concluded that the three mainshocks occurred on shallow reverse faults which have not been recognized before. Here, we investigate rupture directivity of the three Kerman mainshocks and six of the largest aftershocks to estimate the prevailing rupture propagation direction of the seismic sequence. Applying an empirical Green’s function method, we analyzed P and S-waves separately and with different methods, and compared the results of these two, to elucidate and evaluate the consistency of the employed methods and the reliability of the results.

2 Data

Data from the three mainshocks (MS1, MS2, MS3), six of the largest aftershocks (EV 4-9) and nine smaller aftershocks (EGF1-9) have been used for this study (Table 1). The waveform data are velocity seismograms sampled at 50 Hz, collected from broadband stations in a distance range of 80 to 400 km (Fig. 1b). The initial P and S arrival times were estimated using the TauP toolkit (Crotwell et al. 1999) with the IASPEI91 1D velocity model (Kennett 1991). For the P-wave deconvolution analysis, vertical seismograms were band-pass filtered 0.5 to 1.2 Hz using a zero-phase Butterworth filter. A window of 10 s, bracketing the maximum P-wave amplitude, was then used for analyses. To estimate the signal to noise ratio (SNR) of the P-wave, a noise window of length 10 s was defined to end just 5 s before the first P arrival. The SNR value was then measured as the relative amplitude spectra of the noise and signal windows around 1 Hz, which shows SNR values of at least 2 and 5 for the aftershocks and mainshocks, respectively. For the S-wave spectral analysis, a 30-s window bracketing the maximum amplitudes of the observed S-wave was chosen on both transversal and radial components. Displacement amplitude spectra were computed through FFT on unfiltered displacement seismograms. The earthquake locations used are from the IRSC catalog which provides location uncertainties of 0.87 km in latitude and 1.1 to 2.2 km in longitude direction. A maximum epicentral separation of 3 km between the larger event and EGFs was chosen. In order to make sure we have selected appropriate EGF events, we examined the epicentral separations according to the calibrated locations of Savidge et al. (2019). The maximum separation obtained was 4.4 km for mainshock-EGF pairs and 8.7 km for largest aftershocks (EV 4-9)-EGF pairs. Calibrated epicentral uncertainties are all below 3 km at the 90% confidence level. The maximum separation between the events is at most one-tenth of the source to receiver distance and is about three wavelengths at the highest frequency used (2.8 km at 1.2 Hz). Thus, for this geometry and these frequencies, the propagation paths of the target and EGF events should be sufficiently similar to avoid major bias due to path effects. Fault plane solutions are taken from global-cmt catalog (Dziewonski et al. 1981; Ekström et al. 2012) as listed in Table 1.

3 Method

3.1 Theory

While earthquakes are very complex phenomena, an apparently clear unilateral rupture behavior is often seen for both large events (e.g., McGuire et al. 2002; Mai et al. 2005) and small earthquakes (Calderoni et al. 2015), well-observed by station networks. Unilateral rupture on a fault plane produces characteristic azimuthal variations in source signature (Ben-Menhaten 1961), including shorter duration and higher amplitude of the source time function (STF) along the rupture propagation direction, and longer duration lower amplitude signals in the opposite direction. Following, e.g., Haskell (1964), we illustrate the concept using the simplest possible mathematical description of the phenomenon, and assume a line (long, thin) fault in a homogeneous Earth where rupture initiates at one end and propagates to the other, and all slip on each part of the fault occurs instantaneously as the rupture passes. The apparent (observed) rupture duration is as follows (Lay and Wallace 1995):

$$ T(\theta)=\frac{L}{V_{r}}-\frac{Lcos\theta}{c}=\frac{L}{V_{r}}(1-\frac{v_{r}}{c}cos\theta) $$

(1)

Here, T(𝜃) is the apparent duration (s) of the rupture along a line fault measured at 𝜃 which is the azimuth between the direction of rupture and recording stations. L is the rupture length (km), v_r is the rupture velocity, and c is the seismic wave velocity (km/s). Similar simple expressions exist if the fault has significant width as well as length. In this case, the dip of the fault enters the equation, and an assumption of the form of rupture must be made (e.g., a line rupture, or a rupture pattern which is initially circular). Additionally, slip duration on each segment of a fault is not instantaneous. If slip duration is long compared with total rupture time, then essentially the whole fault radiates simultaneously and rupture direction and duration has negligible effect. If slip duration is smaller relative to rupture time (e.g., Heaton 1990) then details of the rupture process may significantly effect the radiation pattern. Further complications may exist including, e.g., if rupture velocity and amount of slip or slip duration vary significantly along the fault, if there are repeated episodes of slip on some segments as the earthquake develops, if adjacent faults with different orientations and slip directions rupture simultaneously, or if wave propagation in the inhomogeneous Earth leads to major directionally dependent perturbations of the waveforms. Such complications suggest that analysis based on simple rupture, slip and wave propagation models may not be adequate to deduce fully reliable descriptions of the source. However, experience from several studies (mentioned above) has shown that the general character of the slip process means that by concentrating on frequencies which are relatively low compared with the rupture- and slip-duration, analysis based on the simple concepts related to Eq. (1) may be adequate.

For the analysis, we must consider the propagation of the signal from the source to each observing station. Ray paths are in general curved rather than straight lines. The wave trains observed are the sum of energy which may have followed significantly different routes through the Earth and may be significantly affected by scattering, phase conversion, etc. Accounting for these effects can be approached by using a velocity model based on previous data to model wave propagation. Alternatively, a more direct empirical Green’s function approach may be used.

The type of analysis we conduct in this study fulfills various criteria suggesting that all of the complications above are manageable. The signal frequencies which we use are relatively low compared with the size of the faults involved; we use relative methods to reduce the effects of wave propagation; and focusing on signals arriving close to the P and S first arrivals reduces the complications related to multi-pathing and phase conversion. For the P-waves, this is because rays traveling by a very different route than the direct wave from source to receiver have traveled significantly farther and will arrive outside the time window used for analyses, as will converted waves which have propagated more than a limited distance as S-waves. Similar logic is valid for the S-waves, though there is the additional complication of possible interference from late P-wave arrivals. In general, this problem is limited because S-wave amplitudes are higher than P.

Different methods exist to suppress the effects of wave propagation in this type of analysis. One approach is to deconvolve the target event with a smaller event, i.e., using the smaller event as an empirical Green’s function, EGF (e.g., Hartzell 1978; Mueller 1985). Assumptions include that the smaller event is sufficiently close to the larger event so the propagation paths are very similar, that the same EGF is adequate for signals originating from all segments of the larger fault, and that the source of the small event is sufficiently small that this event can be regarded as a “point source” at the frequencies used. In practice, the latter implies that the smaller events should be at least one order of magnitude smaller than the larger event. When using such an EGF approach, especially one focusing on signal amplitudes as opposed to time signature (observed duration), there may also be an implicit assumption that the focal mechanisms of the larger and smaller events, and thereby the P and S-wave radiation patterns, are similar. If these requirements are appropriately met for a given data set, deconvolution can provide a reliable relative source time function (RSTF) which includes interesting features about source properties such as rupture extent and propagation direction. Deconvolution can be performed either in the time or frequency domains. Extracting detailed source properties needs careful data analysis and application of appropriate deconvolution methods (e.g., Hough and Dreger 1995; Vallee 2007), but according to, e.g., Vallee (2004), the width of the retrieved RSTF is fairly stable and insensitive to the deconvolution process.

Deconvolution in the frequency domain is straightforward, being based on spectral division together with some stability constraints. In the frequency domain, the directivity manifests itself as an azimuthally dependent corner frequency, which appears as higher spectral amplitudes and corner frequencies for signals aligned with the direction of rupture (Calderoni et al. 2015).

3.2 Application of the method to the dataset

The observed signal amplitudes at a station are affected by the size of the source, the radiation pattern, and the wave propagation path. We focus first on relative observed initial durations as a function of geographical position with respect to the source and later compare this to a separate analysis based on observed spectral amplitudes. In an EGF approach, we deconvolve the waveform of the larger event using that of the smaller event. With the underlying assumptions about source behavior (discussed above), at the frequencies we use (0.5–1.2 Hz), the radiation pattern for the smaller event will be close to that of a point source. That of the larger event will be slightly different in character because of interference between time-delayed signals from different parts of the fault, which will distort the radiation pattern from that of a point source, making it asymmetric. If the slip orientation of the large and small events are similar, the normalization achieved via deconvolution should reveal any significant such asymmetry which is present. As the large and slightly smaller events used both reflect the stress regime in the source region, we might expect the slip directions to be similar. However, the orientation of the radiation patterns may sometimes differ significantly. Since for the studied seismic sequence no focal mechanism solution is available for the aftershocks in the global-cmt catalogue, nor in the published studies (except EGF9 in Savidge et al. 2019), so we choose three EGF events for each larger target event and considered waveform similarities between them (Fig. 2a). To further reduce the possible risk of bias due to differences in radiation patterns of the large and smaller events, we choose to initially focus on observed (deconvolved) relative P-wave duration between the stations. This should be less sensitive to possible radiation pattern mismatch than an approach based on signal amplitudes.

Fig. 2

Waveform deconvolution procedure at station TVBK with azimuth 215^∘ for MS1. a Waveform similarities between the mainshock and EGF events on the vertical component. The boxes mark the P and S-wave windows used for analysis. b Normalized waveforms showing similarity of the mainshock with three EGF events, in 10-s P-wave window, filtered 0.5 – 1.2 Hz. c 3s of the EGF events was chosen to be deconvolved from the mainshock. d Waveforms obtained after the deconvolution. The dotted waveform is the absolute value of the three waveforms summed up together. The green dotted trapezoid marks the retrieved RSTF

For the P-wave analysis, after considerable experimentation, a filter band 0.5 to 1.2 Hz was applied prior to extraction of the data segments. This choice of filtering appeared to provide the greatest similarities between the waveforms of the mainshocks and the associated EGFs. A window length of 10 s was chosen for all the stations. This 10 s is the whole P-window length before the direct S phase arrival for the closest station (80-km epicentral distance). Then 3 s of the EGF waveform at each station was chosen in a way to include the most similar part of the EGF events. A conventional Hanning taper was applied and the signal was zero-padded to keep the 10-s length (Fig. 2c). After transformation to the frequency domain, deconvolution was performed by spectral amplitude division with water-level (1%) regularization and then transformed back to the time domain. The derived RSTF waveforms (Fig. 2d) provide information on the rupture directivity of the mainshock in the form of azimuthal variations in the width (duration) of the observed pulses.

For the subsequent S-wave analysis, we applied a different method, to use both a different subset of data (S rather than P) and an alternative approach which would allow an additional assessment of the consistency of results. A window length of 30 s was chosen for both the radial and transverse components of the seismograms (Fig. 2a). This length was chosen so as to include the significant part of the S-wave energy at the more distance stations ($\sim 400 km$). The data segments for the large and small events were the same length (i.e., no zero padding). As for the P-wave analysis, the EGF deconvolution was performed using spectral division, with 1% water-level regularization for numerical stabilization. Unlike in the P-wave analysis, we did not transform back to the time domain but instead considered spectra averaged separately over stations at one end of the fault, compared with stations at the other end (almost the same approach as used by Calderoni et al. 2015). We chose stations within ± 45^∘ aligned with the fault strike. As shown in Eq. (1) and assuming a constant rupture velocity of 50% of the P-wave velocity, this window would include apparent relative rupture durations of 0.5(L/vr) to 0.65(L/vr) in the direction of rupture and 1.35(L/vr) to 1.5(L/vr) in the opposite direction. Slip duration is not included in Eq. (1). Longer slip durations will lead to smaller differences in relative apparent rupture durations, but for unilateral rupture, the directionality pattern should remain. We present results from the averaged spectral amplitude comparison at frequencies 0.5-1.2 Hz, the same frequency band used for P-wave analysis. This frequency band is sufficiently high to avoid point source approximation for the larger events, and is adequately low that the EGF events will effectively be point sources (note the spectra and the box that marks the mentioned frequency band on Figs. 5 and 6).

3.3 Synthetic analysis

In order to assess, for our particular cases, the reliability of the method for estimating the width of the RSTF, a number of synthetic tests were conducted using realistic EGFs and synthetic main events. For the main shock, we assumed a simple line source with 10-km length, Vp 6.0 km/s, and Vr = 0.5Vp. The rupture was assumed to start from one end of the line source with strike 310^∘ and propagate through the other end. A constant slip history of 1 s was also assumed and included in the model. The source time function (STF) at each azimuth was calculated based on Eq. (1). The synthetic mainshock was produced by convolution of the source time function with real waveforms of an EGF event (EGF4 in Table 1) at every station. Another EGF event (EGF5 in Table 1) was then selected to be deconvolved with the synthetic mainshock. The same deconvolution procedure as used for the empirical P-wave data was then used. To estimate the width of the obtained RSTF (Fig. 3a), we calculated the envelope of the waveforms and examined different proportions of amplitude reduction from the peak in the envelope as measures of the width of the RSTF. We found that a 35% reduction from the peak value of the envelope was the most reliable measure of the width of the RSTF (Fig. 3a).

Fig. 3

a Retrieved width of the RSTF in the synthetic test. The waveforms are normalized. Dotted waves are the results of the deconvolution, solid waves are the absolute values of the waveforms and gray curves show the envelopes. Dashed lines mark the beginning of the RSTF and empty triangles mark the ends, both determined based on 35% drop from the peak value of the envelope. Arrows and solid triangles mark the real beginning and end of the RSTF, as defined in the test, where it is not equal to the calculated one. b Retrieved width of the RSTF for the first mainshock, marked by the dashed line in the beginning and triangles at the ends

Comparison between the predefined width of the STF and the calculated width shows that the exact values of the rupture duration are not precisely recovered at some stations, but the relative ordering of durations regarding the azimuth is stable. The tests revealed that the width of the input STF might be overestimated at stations in the rupture direction and underestimated at the stations located in the opposite direction. That might be due to the limited frequency band and the simple deconvolution method that has been used here. In spite of these deficiencies, a valid pattern of the azimuthal dependency is recovered, so we use the same duration estimation procedure on the real data.

4 Results

4.1 P-wave deconvolution of the mainshocks

For the first mainshock, MS1, nine suitable stations were available in a distance range of 80 to 300 km (Fig. 1b, inverted triangles). Three EGF events (EGF1, EGF2, EGF3 in Table 1) were chosen based on their estimated epicentral distances from the mainshock (less than 3 km). Results of waveform deconvolution with different EGFs were similar, as we would expect for the data segments used, because the waveforms of the three EGFs were similar. After visual comparison of results from deconvolving with the individual EGFs, the three deconvolved waveforms (one from each EGF) were summed at each station, in order to help suppress noise which appears in the early and late tails of the retrieved RSTFs. For presentation, the resulting waveforms are then normalized to ease the comparison between the different stations. The envelope procedure described above was used to estimate duration at each station. The retrieved width (Fig. 3b) clearly shows shorter length for stations oriented toward NW, than those toward SE.

For the second mainshock, MS2, we had 12 stations (Fig. 1b, triangles), which improved the azimuthal coverage. Events EGF4, EGF5, EGF6 in Table 1, were used for the waveform deconvolution (Fig. 4a). A shorter RSTF width at stations around azimuth 300^∘ is distinguishable.

Fig. 4

Retrieved width of the RSTF for a the second and b the third mainshock. Lines and marks as in Fig. 3b

For the third mainshock, MS3, we had the same station distribution as the MS2, thus a good azimuthal coverage compared with the first event. Unlike the MS1 and MS2, this time the waveform deconvolution results do not show a clear preferred rupture direction. The width of the RSTF is almost equal in the NW and SE directions, suggesting a more complex rupture pattern, consistent with bilateral rupture propagation along the fault strike.

4.2 Spectral analysis of mainshocks with the S-waves

Analyzing the relative displacement spectra of the S-waves provides a simple approach to indicate along-strike directivity, based on spectral amplitude rather than time duration. For non-strike-slip events, up-dip or down-dip rupture propagation may dominate changes in rupture duration with direction. However, we must also consider the aspect ratio of the fault segment which has slipped. If the along-strike length is greater than the along-dip length, then rupture duration distortion will be more apparent along strike, irrespective of slip direction. In our case, InSAR investigations of the triplet (Savidge et al. 2019) indicated a highly elongate slip distribution, much narrower downdip ($\sim$5 km) than along strike ($\sim$15 km), for the first two mainshocks. For the third mainshock, their obtained slip model is less stable compared with the first two, but suggested a $\sim$7 km along strike with $\sim$3 km downdip fault plane. These values imply an aspect ratio of 2 and higher for these reverse mainshocks and implying that the along-strike rupture directivity investigations are appropriate.

As mentioned above, EGF methods assume a similar slip orientation for the target (larger) and EGF events. In some directions relative to the slip, the S-wave radiation pattern implies low amplitudes, while radiated energy in a directions parallel and perpendicular to the fault slip is higher. Close to these directions, the rate of change with azimuth of energy in the radiation pattern is relatively small. Thus, if the slip direction of the small event does not differ that much with the larger event, then the normalized (deconvolved) energy observed in directions close to parallel or perpendicular to slip should reveal the asymmetry in the radiation pattern of the larger event. For both pure strike-slip and pure dip-slip events, one of these directions is along strike. In our analysis, we choose to normalize separately with data from several smaller events, allowing assessment of the stability of results, and increasing the likelihood that the final averaged result will not be significantly biased by an inappropriate slip direction for the small event. Furthermore, for this S-wave data, in addition to averaging the spectra of stations aligned with the fault strike (blue and red spectra in Figs. 5 and 6), we also calculated the average spectra from the stations roughly perpendicular to the fault strike (black spectra in Figs. 5 and 6) for comparison. For unilateral rupture propagation, the black curve should lie between the blue and red ones, otherwise, either the rupture is propagating asymmetrically in different directions, or there is significant mismatch between the fault plane orientations of the target and EGF events. Accordingly, we reject results which show a significantly large or small value of the black curve relative to the blue and red ones. For a more convenient comparison of the averaged spectra and visual demonstration of the rupture direction, the spectral values are summed up in a band width of 0.5–1.2 Hz and plotted as bar graphs (Fig. 5 insets). The frequencies considered (marked as black boxes in Fig. 5) can be seen to be above the corner frequencies of the mainshocks (i.e., sharp change in the slope of the spectra), and below the higher frequency asymptote.

Fig. 5

Results of the S-wave spectral division to detect rupture directivity for the three mainshocks. Red is the average spectra of the stations around NW end of the fault plane, blue is the average around the SE end and black is the average for the perpendicular stations. The bars in the insets show the mean amount of the amplitude spectra in frequency range 0.5 to 1.2 Hz, using the same color scheme as the spectra curves. The target and EGF event related to each plot are specified in the bottom of the plots. Event ID’s are based on Table 1

To estimate the average spectra, a minimum of three stations were used for each direction (i.e., in the NW, SE and approximately EW directions; Fig. 1b). However, the station coverage did not provide more than four stations for averaging in each direction. For the first mainshock, MS1, using EGF1 the obtained results are considered unreliable since the average spectral amplitude in the stations perpendicular to the fault strike (gray bar) is larger than the average spectra in either NW or SE direction (Fig. 5a). For event EGF2, the S-wave data was corrupted at some stations and we could not use it for the analyses, but with EGF3 the results show clear directivity toward the NW (Fig. 5b). For the second mainshock, MS2, the relative spectral amplitude of the S-wave data indicates a clear directivity toward the NW with all the three EGF events (Fig. 5c, d, e). For the third mainshock, MS3, we get different results with different EGF events; the analysis with EGF7 suggests a uniform rupture propagation extension (Fig. 5f), with EGF8 the apparent rupture directivity is toward SE and with EGF9 the rupture apparently propagates to NW (Fig. 5g, h). According to Calderoni et al.’s investigations, if a target event does not have any directivity, then the relative spectral analysis would be affected by the directivity of the EGF event instead. Here, we consider the discrepancy of the results as an indicator of the lack of significant along-strike directivity for MS3.

4.3 Spectral analysis of significant aftershocks

To further investigate if there is a generally preferred rupture propagation direction for the studied area and the 2017 seismic sequence, we analyze directivity of the six largest aftershocks (magnitudes 5.0–5.2) which occurred from a few minutes after the first mainshock until a month after the third mainshock (listed in Table 1 as EV4 to EV9). We chose to limit our analysis to the S-wave spectral deconvolution since the process is much faster than the P-wave analyses and requires less analyst interaction. Based on the relative epicentral separation and magnitude differences, the EGF3, EGF4, EGF5 events were used as the empirical Green’s functions for events EV4 , EV5, EV8, and EGF9 was used for events EV6 , EV7, EV9 (Fig. 6). Even though the magnitude differences between the target and EGF event are smaller here (0.8–1.2), the same frequency band of 0.5–1.2 Hz as used for the mainshock investigation, was used. The frequency band is still considered appropriate since it excludes both the initial part and steep part (high frequency asymptote) of the spectra. For EV4, the first analysis (Fig. 6a) shows very small values for the black curve, violating our test criterion. The two other analyzes (Fig. 6b, c) suggest almost a bilateral rupture propagation. A significant directivity toward NW is inferred for EV5 (Fig. 6d, e, f) and events EV6, EV7, EV9 (Fig. 6j, k, l, respectively). For EV8, the gray bar is above the other two bars in analyses using EGF3 and EGF4 (Fig. 6g, h), however, the difference is not significantly large with EGF3 and we consider the test as reliable, which indicates rupture propagation toward SE. With EGF4 the difference is remarkable and the test is unreliable (Fig. 6h). With EGF5 the results exhibit a clear directivity to SE (Fig. 6i).

Fig. 6

Results of the S-wave spectral division to detect rupture directivity for the magnitude 5 aftershocks. Colors and insets as in Fig. 5

4.4 Reliability analysis

In the P-waveform deconvolution, we selected a 3-s window on the EGF waveforms for deconvolution with the mainshock. The 3 s was chosen in order to include the most similar parts of the three EGFs. Increasing this window length allows arrival of multi-paths and phase conversions which would reduce the similarities between the EGF waveforms. Furthermore, by increasing the length of the EGF window the resulting deconvolved wave will get more complicated and thus more challenging to detect the RSTF. This is shown in Fig. 7a where we extended the EGF window to 5 s for the KBAM station and the RSTF could not be retrieved successfully. However, to make sure that the choice of the window length has not introduced bias into our relative rupture duration measurements, we examined longer EGF windows at a selected station in which the waveform similarities last for a longer time period. The obtained RSTF for the first mainshock at station TKDS is shown on Fig. 7b,c,d, using EGF windows of 3, 5, 7 s, respectively. An estimate of the RSTF width was successfully recovered for all of the window lengths and the results exhibit a stable width measurement where the calculated width is equal for the 5 and 7-s windows and negligibly larger compared with the 3s window.

Fig. 7

The left column shows the selected EGF window to be deconvolved from the first mainshock and the right column is the resulting deconvolved waves. a Testing a longer EGF window (5s) for station KBAM where RSTF could not retrieved successfully. b–d Recovered RSTFs for station TKDS with EGF window lengths of 3, 5 and 7 seconds, respectively. Arrows mark the width of the RSTF

In our analysis, the P-waves were preferred for the RSTF calculations and S-waves for the spectral analysis. P-waves can be easily picked on the waveforms since they are first arrivals and by choosing the initial part of the waveforms, we can largely exclude energy from P-S conversions and multi-pathing. However, the S-waves are generally contaminated with conversed phases and coda waves, which can make it challenging even to detect the initial S arrivals. For this reason, the S-wave data (waveform sections) should be prepared very carefully if the intention is RSTF calculations. For comparison of methods and to examine our RSTF analysis and validity of the obtained results, we used S-wave data to recover RSTF for the second mainshock. The method was repeated as for the P-wave analyses with a careful data selection to pick the similar parts of the waveforms. The obtained results are shown in Fig. 8a, where a clear rupture directivity toward NW is noticeable.

Fig. 8

a Retrieved width of the RSTF for the second mainshock using S-waves. Lines and symbols as in Fig. 3bb Displacement spectra of the P-waves for the second mainshock deconvolved with the EGF5 event. Individual colors assigned for each station. c Displacement spectra averaged in strike directions. Colors and insets as in Fig. 5

The length of the P-wave window, on the other hand, is proportional to the epicentral distance of the station and thus the S-P time. That means a shorter time window would include a considerable amount of the P-wave energy at a closer station, while a longer time window is needed to include the same amount of energy in a more distant station. Since an equal data length is required for the spectral comparison, this would result in higher spectral energy for the closer stations even after the deconvolution. In this case, when averaging the spectra according to the azimuthal distribution, caution should be made to ensure a similar variety of epicentral distances at each group, to include comparable energy levels. In our dataset, the minimum epicentral distance is almost 80 km which limits the P-window to about 10 s, whereas for our maximum epicentral distance (400 km), the P window (before direct S-wave arrival) is about 47 s. The amplitude spectra plot for 8 s of the MS2 deconvolved with EGF5 (Fig 8b) shows significantly higher amplitudes for the closest station (ZRDN, the black spectra) and noticeably lower amplitudes for the furthest station (ANAR, the yellow spectra). Excluding these two stations from the spectral analysis, we calculated the average spectra according to the fault strike (Fig. 8c) and as expected rupture directivity toward the NW is detected.

5 Discussion

The EGF deconvolution method simply reduces the influence of the path and site effects on the waveform data, allowing us to better study some properties of the source. Because of the limited station distribution in our dataset and the assumptions and simplifications implicit in these types of analyses, we considered it appropriate to use “independent” data sets (P and S data separately) and different analysis methods, where getting consistent results would support both the reliability of the methods and of the obtained results. For the P-waves, our RSTF estimate involved filtered and deconvolved data, finding waveform similarities and an envelope-based algorithm for assessing pulse duration, while spectral deconvolution of the S-waves was a direct approach and just needed a proper averaging over the stations along the fault directions. However, the P waveform analysis has some advantages, including theoretical aspects such as lower sensitivity to assumptions about radiation patterns, and that data can be evaluated for every individual station allowing detailed examination of geographic variation of the pulse width.

The time and slip history of rupture is affected by a variety of parameters such as fault plane geometry (e.g., Harris and Day 2005; Poliakov et al. 2002), the material properties of the fault area (e.g., Andrews and Ben-Zion 1997; Cochard and Rice 2000) and non-uniform stress distribution on the fault (e.g., Duan and Oglesby 2007). Despite all these complexities, a rupture usually has a preferred propagation direction which affects strong ground motion and hazard estimates. In the studied sequence, the triplet mainshocks show very similar focal mechanisms and are closely located in time and space. Two have similar directionality (rupture) patterns, but the third is different. We propose two scenarios to describe the phenomenon;

• I. One assumption is that all the three mainshocks are dominated by rupture on the same fault plane. The first two mainshocks ruptured the southeastern segment of the assumed fault plane, and both propagated toward NNW, where the third mainshock occurred (Fig. 9). The third mainshock, located about 10 km NNW from the first ones, ruptured almost bilaterally on the northeastern segment of the fault. The aftershock distribution immediately after the second main shock shows migration of the aftershocks from the SE segment to the NW (Fig. 9). In this case, we could infer that preexisting barriers, motivated by the stress history, might have impeded the rupture propagation to go further NW during the first mainshock. Later, the stress distribution changes caused by the first and second mainshocks might facilitated the rupture on the NW segment, leading to the third mainshock occurrence.

• II. The second scenario assumes that the three mainshocks may involve simultaneous slip on different, close lying, faults. Depending on the local stress regime and other factors, simultaneous significant slip on more than one fault may occur, and these faults may have different orientations. To investigate this, we considered the area on a larger scale. Generally, the current tectonic state of the study area is being controlled by northward motion of the Arabian plate toward the Eurasia. The resulting deformation is accommodating by crustal shortening along the reverse faults in the mountain belts, and by the lateral displacement along the large strike-slip faults surrounding the rigid blocks (e.g., Jackson and McKenzie 1984; Vernant et al. 2004). In the area where the triplet happened, three large right-lateral strike-slip faults converge; the Gowk fault to the southeast, Kuh-Banan to the northeast, and Lakar-Kuh to the north (Fig. 1a). The overlapping arrangement of these three strike-slip faults produces a 40 km wide restraining bend where reverse faults can develop (Walker et al. 2010). The triplet is located very close to the Lakar-Kuh fault. This fault changes orientation from N-S to NNW by forming several closely spaced NNW trending splay faults in its southern end, producing a restraining bend in connection to the Gowk fault. Judging from the distribution of the aftershocks, the seismogenic area related to these events in the restraining bend is ca. 20 km wide. There is another N-S trending strike-slip fault (Nayband) which clearly cuts many structures and propagates much further south than Kuh-Banan and Lakar faults. Nayband seems to outline the eastern boundary of this restraining bend. Savidge et al. (2019) suggest that the triplet occurred on presumably young, shallow reverse faults. Such faulting is a common feature in restraining bends, which could grow with time and form new structures to accommodate displacement along the main strike-slip faults.

The main stress direction due to the Arabia-Eurasia convergence is about N10E (Vernant et al. 2004), west of the Lut-Block. The orientation of the MS3 fault is 287^∘ which is slightly different compared to the MS1 and MS2 events, 310^∘. This makes an acute angle (almost 60^∘) between the stress direction and fault strike in MS1 and MS2, but almost a right angle (almost 83^∘) with MS3. Mechanically, such orientation relative to the main stress direction might make the faults more prone to unilateral rupture propagation toward NW in the first case (MS1 and MS2), and would encourage bilateral rupture propagation in the second one (MS3).

The aftershock distribution shows that the largest aftershocks of the sequence are located almost at the two ends of the seismic sequence; three in the SE (EV4, EV5, EV8) and three in the NW (EV6, EV7, EV9; Fig. 7). All these aftershocks except two (EV4, EV8) are propagating toward the NW. These two aftershocks almost mark the southeastern limit of the ruptured area and might be affected by the changes of the fault plane orientation, leading to a different rupture propagation direction.

Besides all the similarities of the three mainshocks including magnitude, focal mechanism, and epicentral locations, the available depth estimates have both a considerable range and large uncertainties which prevents a meaningful comparison of hypocenters), there are some significant differences between the third and the first two. One is the difference in the rupture directivity that was mentioned above, and the other is the duration of the source pulse. Comparing the estimated RSTFs of the three mainshocks (Figs. 3b, and 4a, b), we clearly see that the source duration for the third mainshock is smaller. Savidge et al. (2019) also obtained a smaller rupture plane for the third mainshock. Furthermore, a 6–7-km-long surface rupture was generated during the third mainshock, while the other two are blind faults. Since the mainshocks have similar magnitudes, a relatively shorter source duration and extension of MS3 could imply higher slip values and, as suggested by Savidge et al. (2019), probably a higher stress drop.

Fig. 9

Distribution of the mainshocks (stars) and the aftershocks (circles) occurred in 5 days after each mainshock. Arrows demonstrate the detected directivity direction. Faults are from Walker et al., (2010). The red line above the third mainshock is the surface rupture detected for the third mainshock by Savidge et al. (2019). Event ID’s are based on Table 1

6 Conclusion

The basic concepts of directivity and EGF deconvolution are used in this study to investigate rupture propagation directions of a triplet in Kerman province. According to our analyses, the ruptures of first two mainshocks propagated toward NW, while for the third event propagation was almost bilaterally NW and SE. The directivity of the magnitude 5 aftershocks was also investigated and 4 out of 6 of them showed directivity toward the NW. This implies that the predominant rupture direction of the sequence is from SE to NW, suggesting that the Lakar fault is tending to bend and link to Gowk fault through a number of closely spaced splays which produced the three mainshocks and the following seismic sequence. The rupture propagation direction in the area seems to be driven by the orientation of the faults relative to the prevailing stress field where the fault orientations are affected by the local stress redistribution in a restraining bend.