Temporal properties of aftershock sequences of large earthquakes in Iran - Analysis of primary and secondary aftershocks of the Ezgeleh sequence

In this study, the decay of earthquake aftershock sequences of some major earthquakes in different tectonic regimes in the Iranian plateau is discussed. The studied earthquakes are Rigan [2010], Ahar-Varzaghan [2012], Goharan [2013], Sefidsang [2017] and Ezgeleh [2017]. The spatial and temporal windows are considered based on the method proposed by Gardner and Knopoff [1974] to compute decay parameters for each sequence. The decay rates of sequences were compared to well-known models to find the best fit for each sequence. The results showed that the modified Omori is the best fit for Ahar-Varzaghan and Ezgeleh sequences, for Rigan and Sefidsang sequences the modified Omori and the Kisslinger ones found as the best fits. The values of the p parameter of the Reasenberg and Shcherbakov models were larger compared to the Omori model, but the parameter of the Kisslinger model was slightly smaller compared to the Omori one. The c parameter showed an inverse relation to the threshold magnitude. The correlation between the p and c parameters and also the and the Gutenberg and Richter (G-R) parameters were investigated. In addition, we made use of a graphical method to analyze the seismic sequence of the Ezgeleh earthquake during 13 months after the main event. The graphical method was successful to estimate the occurrence of an event with an approximate magnitude of M=6.4 in the sequence.


Introduction
Seismic sequences are referred to as a series of earthquakes in a region that occur over a certain period of time.
Typically, any seismic sequence consists of a main event and its associated events; foreshocks and aftershocks. The important parameters of any seismic sequence are the number of aftershocks, spatial distribution and the reduction of the number over time. The design of procedures for identifying and differentiating these stages in the seismic regime can help a better understanding of the earthquake phenomenon and consequently can reduce the impacts of earthquakes. Aftershock events are associated with smaller magnitudes and with a time interval of between several minutes to several years after the main seismic earthquake, and their abundance decreases over time [Omori, 1894].
Aftershocks can be affected by various factors. The simplest one could be considered as when the stress in a region increases, not all of the stress is released by the main event, but some part is released by aftershocks [Das and Scholz, 1981a, b]. The concentrations of stress resulting from asperities and barriers can create aftershocks [Scholz, 2002]. In some cases, the adjacent faults to the main events are broken down and can lead to aftershocks or more energetic events than the previous events. Of course, in some cases, seismic waves caused by the main event can lead to change the stress of the fluids into adjacent fault holes, which reduces the strength (static friction) on the adjacent faults and creates aftershocks [Shcherbakov et al., 2005]. Also, as stated by Riga and Balocchi [2017] the aftershocks are classified by Kisslinger [1996] into three sets: aftershocks occurring in the rupture area of the main fault plane, aftershocks occurring in the main fault plate, but outside the rupture area and aftershocks outside the main fault.
Aftershock studies are considered as a major step towards knowing better the physical process of earthquakes [Kisslinger, 1996]. Many studies are carried on the temporal and spatial distribution patterns of aftershocks. One of the oldest and most important researches in this regard is carried by Utsu [1969]. Also, Reasenberg [1989] studied the California earthquakes and presented a relation to combine the two experimental relationships of Gutenberg and Richter [1944] and Omori. Later, Kisslinger [1993] studied the dependence of the aftershock distribution on physical properties of fault zones and environmental conditions, especially resistance, stress and temperature.
Then, Shcherbakov et al., [2004] presented a generalized formula for aftershock decay rates, by combining the three relations between the Omori law, the modified G-R and the modified Bath's law. Gasperini and Lolli [2009] performed a comparison between decay rate models in Southern California and Italy based on the modified Omori model, the modified Kisslinger model and the band limited power law. They resulted that the modified Kisslinger model and the band limited power law generally explain decay rate aftershocks better than the modified Omori model. More recent studies on the spatial and temporal analyses of the distribution of aftershock sequences carried on in different parts of the world; Wiemer and Katsumata [1999], Bayrak and Öztürk [2004], Öztürk et al. [2008], Nuannin et al. [2012] and Öztürk and Şahin [2019]. Ommi et al., [2016] determined the mean values of coefficients of the Omori law for 15 mainshock-aftershock sequences in different seismotectonic provinces of Iran. Also, Ansari [2017] studied the spatial and temporal seismicity parameters, stress state and seismic energy of the Shonbe earthquake sequence.
In this study, decay rates of the earthquake aftershock sequences of some major earthquakes in different tectonic regimes in the Iranian plateau were discussed. The studied earthquakes are introduced in Table 1 The primary purpose of this study is analyzing the decay rates for different sequences of some large earthquakes in Iran applying some well-known decay models. For the second purpose, we studied the dependence of the parameters of the Omori law on threshold magnitudes and temporal windows. Moreover, the first primary and secondary aftershocks are proposed for the most recent sequence of Ezgeleh 2017, applying a graphical method proposed by Riga and Balocchi [2016].

Tectonic setting
The Iranian plateau is located in the Alps-Himalayas seismic belt, one of the most active and seismic regions on the Earth, in the collision zone of the Arabian, Indian and Eurasian plates. The collision of these plates and the pressure they have caused deformations, cracks, mountains and earthquakes. The seismiotectonics of Iran were studied by individuals, including Berberian [1981] and Nowroozi [1976]. All have come to terms with the fact that the Iranian plateau contains different seismiotectonic regimes. In this study, we used seismiotectonic provinces proposed by Mirzaei et al., [1998]. They subdivided the Iranian plateau into five major seismotectonic provinces:

Decay models
There are different methods for separating foreshocks, main shocks and aftershocks. Van Stiphout et al., [2012] made a complete review of different declustering algorithms. We applied the Gardner and Knopoff [1974] spatial and temporal windows. According to studies of Rezapour and Mohsenpur [2013] and Walker et al. [2013] the Rigan earthquake included two seismic sequences occurred on two cross-faults, we considered the first event in a period of 39 days before the second sequence was started. The Ahar-Varzaghan earthquake has a secondary sequence that occurred approximately 90 days after the main shock. The aftershocks migrated in both along-strike and up-dip directions [Rezapour, 2016]. We studied the Ahar-Varzaghan earthquake in a 90 days period, this is why the time windows of the two earthquakes mentioned are not accordance with the time windows in Gardner and Knopoff [1974], as presented in Table 2.
One of the most important factors in measuring seismic parameters is the threshold magnitude. The threshold magnitude is defined as the lowest magnitude at which 100% of the earthquakes in a space-time volume are detected [e.g. Mignan and Woessner, 2012]. This parameter was determined for different seismotectonic provinces of Iran by Mousavi et al. [2014] and Ommi et al. [2016]. In order to determine the spatial and temporal distributions of earthquakes, the seismic pattern, first the G-R law was investigated for each event and the a, b and parameters were calculated applying the ZMAP package Wiemer [2001] as shown in Figure 2. The ZMAP software derives the best fitted model using the Akaike Information Criteria [Akaike, 1974], denoted as AIC, the best model corresponds to the model with lower AIC value. The parameters of the model are derived using the maximum likelihood using the likelihood over the time span of the sequence. We examined the impact of different values of threshold magnitudes and the behavior of their aftershock sequences is compared with different decay models; Omori [1894] which later modified by Utsu [1961], Reasenberg and Jones [1989], Kisslinger [1993] and Shcherbakov et al., [2004]. Then, the parameters and decay rates of the abovementioned sequences were discussed and analyzed. There are two empirical laws that can be applied to the aftershock sequences; 1) The frequency-magnitude distribution of the aftershocks following a main event can be modeled by the G-R law as in relation (1); (1) where is the number of earthquakes with magnitudes greater than or equal to , and is the slope of the frequency-magnitude distribution.
2) The magnitude difference between a main event and its largest recorded aftershock according to the Bath's law [Bath, 1965] in relation (2) is: ( 2) where is the magnitude of the main event and is the magnitude of its largest recorded aftershock.

Iran large earthquakes aftershock sequences
The aftershock sequences are discussed based on different decay models (relations 3 to 7), and the relevant parameters are extracted. Relation (3) is the rate of aftershock occurrence in time , follows the modified Omori law proposed by Utsu [1961] as where the parameter indicates the speed of the decay rate of aftershocks. Also, the slope of the graph decreases more rapidly with increasing the value, until the decay rate of the aftershock sequence reaches the background seismicity. It has been shown by Mogi [1967a, b] that the decreasing rate of aftershock activity ( ≥1.3) occurs in the regions where the crust is warmer. Additionally, changes in the value can be related to the degree of heterogeneity of the fault or the liquid fluid in the fault, which weakens and reduces the shear strength of the fault [Nur and Booker, 1971]. Areas that experience greater slip during the main event specify larger values. The c parameter represents the behavior of missing aftershocks in the early stage of an aftershock sequence. Mostly as in the initial stages of sequences of aftershocks the smaller events overlap on the seismogram or because the major earthquake faults have a large seismic motion and activity during the first hours, then the apparent reduction of small events near to epicenter locations, reduce in the very beginning times of the catalog of earthquakes. At this part, aftershocks seem not to follow a decreasing rule. The Omori , c and parameters can be estimated by maximizing the loglikelihood function [Ogata, 1983;Utsu, 1995]. Reasenberg and Jones [1989] presented a parametric model for describing earthquake sequences based on the California earthquakes and the probability of occurring any larger earthquake or strong aftershocks. This model also presents a method for estimating the probability of occurrence of strong aftershocks or larger earthquakes at any time interval. The model is based on a nonhomogeneous Poisson process in time, that follows the modified Omori law, as well as the distribution of magnitudes that follow the G-R relation. Then, the decay rate, , of aftershocks with magnitude or larger, at time following the main event with magnitude is according to relation (4) as where , , and are constants.
On the other hand, Kisslinger [1993] applied the stretched exponential relaxation function, presented by Williams and Watts [1970], to the aftershocks. Kisslinger showed the decay rate of aftershock sequences follows to relation (5) as (5) where N*(0) is a finite number of potential events at the beginning of the sequence, 0 is the relaxation time for the overall process, which has not yet been related to a specific physical mechanism and is a coefficient that has a value of 0 < ≤ 1. Relation (5) is accepted as a valid universal function which it has been compatible with many aftershock sequences in the world. Then Kisslinger [1993] compared the relation (5) with the modified Omori law and by fixing the parameters to obey the modified Omori law, he revised it as in relation (6), Later Shcherbakov et al. [2004] modified the Bath's law, so that the largest aftershock was obtained by replacing (N ≥ m) = 1 in the G-R relation. They proposed relation (7) which is a combination of the G-R relation, the Bath law and the Omori law as (7) where ( , ) is the decay rate of occurrence of aftershocks with magnitudes greater than per day, * is the modified Bath's law, is the time after the main event, is the main event and and ( ) are characteristic times.

Decay rates of large earthquakes in Iran
The parameters of each decay rate are obtained based on the least square method. All seismological data were taken from the Iranian Seismological Center (IRSC), the magnitude reported by IRSC is , a modified Nuttli magnitude. The decay rates of the studied earthquakes are shown in Figure 3 (a to e). The mean value of the parameter for the studied earthquakes for the primary aftershock, in the Omori law is 1.02, and in the Reasenberg, Kisslinger and Shcherbakov models are 6.48, 0.72, and 16.15, respectively. Due to the incomplete data for Rigan earthquake, the Shcherbakov parameter is not properly estimated. A summary of the results of decay rate models for the studied earthquakes is presented in Table 3, including the deduced R-square values.
7 Iran large earthquakes aftershock sequences The speed of the decay rate has a direct relation with the energy drop of the main event. Therefore, the slope of decay rates of the aftershock sequences could be a proper criteria representative of evacuation of the stored energy.
The dominant decay rate in this research is the Omori relation. The parameter, depends on the seismicity of the region. It is obvious that by increasing the threshold magnitude and shortening the temporal window, the parameter is decreased. For discussing about changes in the and parameters further examinations are required.
According to some findings [e.g. Utsu, 1962;Ranalli, 1969;Yamakawa et al., 1969;Hamaguchi and Hasegawa, 1970], the and parameters are independent of the threshold magnitude . There are also few researches [e.g. Utsu et al., 1995;Motoya and Kitagamae, 1971] showed the dependency of on . We investigated the dependency of the and parameters for Rigan, Ahar-Varzaghan, Goharan, Sefidsang and Ezgeleh earthquakes. Figure 4 depicts the and parameters of the Omori law versus different threshold magnitudes, for the studied earthquakes. According to Figure 4, changes in the and parameters versus different magnitude thresholds are small for the mentioned earthquakes.
In order to investigate the impact of different time intervals of the seismic sequence on the parameters of models, different temporal windows are tested according to Table 4. There are slight changes in the values of the and parameters, particularly differences are obvious for 90 and 180 time windows in the Ezgeleh earthquake.
8    Table 4. The percentage of relative differences between the parameters , and vin different temporal windows.

Mahdieh Lavasani and Elham Shabani
Also, the correlation between each couple of , and parameters were studied for Goharan, Sefidsang and Ezgeleh earthquakes, the results are presented in Figure 5. Where the stars indicate the time period of 180 days and the plus signs are assigned to the time period of 90 days. According to Figure 5, the correlation between and parameters also between and the G-R parameters for Goharan and Ezgeleh earthquakes are not constant over the time. The correlations between with , , parameters for different threshold magnitudes are not constant at different threshold magnitudes.

The primary and secondary aftershocks of the Ezgeleh earthquake
The 2017 Ezgeleh earthquake of magnitude 7.3 occurred in a region where large earthquakes have not been documented for several centuries. Understanding fault locations, geometries, and seismic behaviors are highly noticeable in this region [Gombert et. al., 2019]. In this regard, we studied the sequence of the Ezgeleh. The epicenter and aftershocks of the Ezgeleh earthquake are depicted in Figure 6. Also the seismicity map of the region is presented in Figure 7. The type of aftershocks of the Ezgeleh earthquake is analyzed through the graphical method introduced by Riga and Balocchi [2017] and the primary and secondary aftershocks were distinguished. Figure 8 presents the branch structure for the Ezgeleh earthquake during 15 months; starting from three months before the occurrence of the main event to 12 months after that. The main event = 7.3 is marked with a red star and the colored lines indicate the developmental stages, which include the black, pink, green, light and dark blue lines representing the first to the fifth order of energy accumulation stages for the aftershock sequence of Ezgeleh. The pink stars with the letter F represent foreshocks for the next large earthquake with = 6.4. According to Figure   8, there is a gradual decrease initially in the magnitude values, after the pink F named points the magnitude of events are increasing up to the occurrence of the next large earthquake with = 6.4. If a straight line is drawn from the midpoint of the fifth order branch, an estimate of the minimum magnitude in the energy release stage is made. This predictive method can confirm the accuracy of our words to express that the F points could be considered as foreshocks of the subsequent large earthquake = 6.4. Figure 9 presents a simple, convenient method according to Riga and Balocchi [2016], to identify the first relative minimum value after the main event = 7.3, the Ezgeleh earthquake. This minimum could be considered as the beginning of the secondary aftershocks in the sequence. The point is the first relative minimum event = 4.1 relative to the foreshock = 3.7.

Conclusion
Our study shows that the Kisslinger or the modified Omori models are the best-fitted decay rates to the aftershock sequences of selected large earthquakes in Iran. Although the Kisslinger exponential model shows poorly fitting with the data at the beginning of the sequences, but over time the Kisslinger decay rate shows better fitting with data compared to the Omori decay rate. This could be due to the parameter which controls the number of calculating events at the times immediately after the main event. Generally, in most of the studied earthquakes, the decay rate of the Omori in the short-term and the Kisslinger at the long-term show the best fitting with the sequence data. The behaviors of the two models in the long time are easily recognizable. The results show that the and parameters of the Omori law depend on the threshold magnitude and the chosen time window. And the parameter decreases with increasing the threshold magnitude, because the larger threshold magnitude value makes the smaller number of aftershocks to be eliminated at the beginning of the sequence in the catalog, therefore, the parameter decreases over time. However, we found no systematic pattern of changes for the parameter. The branch structure method [Riga and Balocchi, 2017] on the Ezgeleh earthquake sequence could successfully estimate the occurrence of an event with magnitude about 6.4, which its evidence was the occurrence of the earthquake = 6.4 on 25/11/2018.