## 1. Seismicity parameters (Mc, a, b)

By applying the temporal-spatial window, the aftershocks related to each mainshock were determined. The effect of aftershocks smaller than Mc (magnitude of completeness) should be removed from the catalog. This was done by drawing the Gutenberg-Richter diagram using the maximum likelihood estimation method in ZMAP software (Wyss et al., 2001). The values of a, b (a-value and b-value) and Mc are listed in Table 2. The range of Mc is from 1.5 to 3.4, while the b-value is from 0.65 to 1.13.

Table 2

a-value, b-value and Mc for each sequence

Number | Event Name | Mag | Mc | Mc+- | a-value | a-value (annual) | b-value | b+- |

1 | Kaki | 6.3 | 2.9 | 0.44 | 4.57 | 4.33 | 0.65 | 0.02 |

2 | Khanaqin | 5.7 | 2.3 | 0.12 | 4.1 | 4.18 | 0.75 | 0.06 |

3 | Khanehzeniyan | 5.4 | 2.6 | 0.23 | 4.96 | 5.18 | 1.13 | 0.2 |

4 | Arad | 5.7 | 2.9 | 0.01 | 5.12 | 5.19 | 1.06 | 0.14 |

5 | Bastak | 5.5 | 3.4 | 0.17 | 5.16 | 5.39 | 1.04 | 0.33 |

6 | Borazjan | 5.6 | 2.7 | 0.13 | 4.18 | 4.27 | 0.85 | 0.1 |

7 | Jam | 5.6 | 2.7 | 0.09 | 4.63 | 4.8 | 1.05 | 0.16 |

8 | Sisakht | 5.3 | 2.3 | 0.18 | 3.33 | 3.58 | 0.7 | 0.22 |

9 | Mormori | 6.2 | 2.5 | 0.1 | 4.82 | 4.58 | 0.7 | 0.04 |

10 | Sarpolezahab | 7.3 | 2.2 | 0.1 | 4.98 | 5.09 | 0.75 | 0.04 |

11 | Somar1 | 5.3 | 2.6 | 0.21 | 4.25 | 4.51 | 0.86 | 0.13 |

12 | Somar2 | 5.6 | 2.4 | 0.11 | 4.37 | 4.45 | 0.71 | 0.04 |

13 | Sangan | 5.8 | 2.3 | 0.22 | 4.47 | 4.45 | 0.91 | 0.15 |

14 | Masjedabolfazl | 5.7 | 2.9 | 0.06 | 4.51 | 4.61 | 1.02 | 0.19 |

15 | Hojedk | 6.2 | 2.3 | 0.08 | 4.89 | 4.7 | 0.73 | 0.03 |

16 | Mohamadabad | 6.5 | 3.3 | 0.11 | 5.42 | 5.13 | 0.94 | 0.09 |

17 | Negar | 5.8 | 2.7 | 0.3 | 4.22 | 4.28 | 1.01 | 0.36 |

18 | Zahan | 5.6 | 2 | 0.4 | 3.98 | 4.08 | 1.07 | 0.36 |

19 | Ahar-Varzaghan | 6.5 | 1.7 | 0.04 | 4.72 | 4.4 | 0.71 | 0.02 |

20 | Qotor | 5.7 | 2.3 | 0.25 | 4.54 | 4.58 | 0.82 | 0.11 |

21 | Torud | 5.9 | 1.8 | 0.05 | 4.43 | 4.35 | 0.85 | 0.03 |

22 | Torkmanchay | 5.9 | 1.8 | 0.07 | 4.43 | 4.35 | 0.85 | 0.04 |

23 | Neyshabor | 5.4 | 1.5 | 0.16 | 3.44 | 3.66 | 0.85 | 0.1 |

24 | Pishqale | 5.7 | 1.7 | 0.13 | 3.45 | 3.48 | 0.72 | 0.06 |

25 | Sefidsang | 6 | 1.6 | 0.1 | 4.29 | 4.14 | 0.66 | 0.03 |

26 | Goharan | 6.1 | 3.2 | 0.22 | 4.46 | 4.61 | 0.66 | 0.08 |

The background b-values were also estimated by a similar method in order to compare these values with the values obtained from the sequence of aftershocks. For this purpose, earthquakes data were collected for each event with a radius of 1 degree from each side and from 2006 to before the occurrence of the investigated events. Figure 2 shows the bar graph of both b-values side by side. Although in most cases, there is not much difference between the values of b, there is a significant difference in cases such as Ahar-Varzghan, Goharan and Khanehzeniyan. The biggest difference is observed in Ahar-Varzaghan; so that the value of background b-value is estimated to be 1.43, while this value is determined to be 0.71 using aftershock.

### 2. Aftershock Parameters (P, C, K)

After determining the spatial and temporal ranges of aftershock occurrence for each main shock and extracting the aftershocks, their cumulative graph was drawn and the best fitted model was obtained. Of course, aftershocks smaller than Mc were removed from the catalog so that incomplete data does not affect the calculation process. As mentioned before, 4 models can be considered. As an example, 2 cases of fitted curves on the cumulative aftershocks diagram are shown in Fig. 2. The results related to each seismic sequence are listed in Table 3 based on these 4 models (based on above-mentioned four models). 14 seismic sequences (more than half of them) followed model 1 and 8 of them were fitted with model 4. Only 4 sequences can be fitted with models 2 and 3. The range of p and c parameters for Iranian plateau is calculated from 0.75 to 2.7 (average 1.32) and from 0.01 to 3.92 (average 0.57), respectively. The range of changes of parameter k is more than the previous two parameters and it is from 10.0 to 1014.2 with an average of 86.7.

Table 3

Values of p, c and k for each sequence

Number | Event Name | Mag | Mc | model | P | C | k |

1 | Kaki1 | 6.3 | 2.9 | 4 | p1 = 1.15, p2 = 2.70 | c1 = 0.39, c2 = 0.64 | k1 = 87.8, k2 = 10.0 |

2 | Khanaqin | 5.7 | 2.3 | 4 | p1 = 0.99, p2 = 2.70 | c1 = 0.08, c2 = 0.41 | k1 = 22.1, k2 = 10.1 |

3 | Khanehzeniyan | 5.4 | 2.6 | 1 | 1 | 0.07 | 12.1 |

4 | Arad | 5.7 | 2.9 | 1 | 0.89 | 0.20 | 10 |

5 | Bastak | 5.5 | 3.4 | 1 | 1.19 | 1.52 | 10 |

6 | Borazjan | 5.6 | 2.7 | 1 | 0.83 | 0.03 | 10.5 |

7 | Jam | 5.6 | 2.7 | 1 | 0.98 | 0.12 | 10 |

8 | Sisakht | 5.3 | 2.3 | 1 | 1.15 | 3.92 | 22.8 |

9 | Mormori2 | 6.2 | 2.5 | 4 | p1 = 1.65, p2 = 0.78 | c1 = 2.39, c2 = 0.02 | k1 = 1014.0, k2 = 17.2 |

10 | Sarpolezahab | 7.3 | 2.2 | 4 | p1 = 0.82, p2 = 1.23 | c1 = 0.66, c2 = 0.30 | k1 = 164.4, k2 = 92.0 |

11 | Somar1 | 5.3 | 2.6 | 2 | 1.14 | 0.66 | k1 = 13.6, k2 = 10.0 |

12 | Somar2 | 5.6 | 2.4 | 3 | p1 = 1.05, p2 = 2.70 | 0.44 | k1 = 115.0, k2 = 10.0 |

13 | Sangan | 5.8 | 2.3 | 1 | 0.75 | 0.01 | 14.0 |

14 | Masjedabolfazl | 5.7 | 2.9 | 1 | 1.32 | 0.30 | 10.0 |

15 | Hojedk | 6.2 | 2.3 | 4 | p1 = 0.94, p2 = 2.70 | c1 = 0.42, c2 = 0.87 | k1 = 136.8, k2 = 17.3 |

16 | Mohamadabad | 6.5 | 3.3 | 4 | p1 = 0.92, p2 = 1.82 | c1 = 0.04, c2 = 0.58 | k1 = 11.5, k2 = 39.4 |

17 | Negar | 5.8 | 2.7 | 1 | 1.11 | 0.48 | 10.0 |

18 | Zahan | 5.6 | 2 | 1 | 0.89 | 0.01 | 10.1 |

19 | Ahar-Varzaghan | 6.5 | 1.7 | 4 | p1 = 1.04, p2 = 1.71 | c1 = 1.74, c2 = 0.83 | k1 = 582.7, k2 = 171.1 |

20 | Qotor | 5.7 | 2.3 | 1 | 0.99 | 0.11 | 45.3 |

21 | Torud | 5.9 | 1.8 | 1 | 1.19 | 0.14 | 38.8 |

22 | Torkmanchy | 5.9 | 1.8 | 3 | p1 = 1.09, p2 = 2.70 | 0.41 | k1 = 134.1, k2 = 10.0 |

23 | Neyshabor | 5.4 | 1.5 | 1 | 1.04 | 0.032 | 15.9 |

24 | Pishqale | 5.7 | 1.7 | 1 | 0.85 | 0.015 | 10.3 |

25 | Sefidsang | 6 | 1.6 | 4 | p1 = 1.09, p2 = 1.07 | c1 = 0.67, c2 = 0.04 | k1 = 304.2, k2 = 11.1 |

26 | Goharan | 6.1 | 3.2 | 2 | 1.33 | 0.86 | k1 = 66.2, k2 = 13.2 |

Earthquakes that have a secondary sequence have a magnitude greater than or equal to 5.9; Except for the Somar 2 and Khanaqin earthquakes, which registered magnitudes of 5.6 and 5.7, respectively. Mainly, in the secondary sequence, the value of P is greater than the primary part, that is, after the secondary sequence, the decay rate of aftershocks increases and it will reach the background seismicity of that area more quickly.

K and c parameters depend on the frequency of events in most cases. In order to better display and compare the changes of these parameters, we normalized them in the range of 0 to 1 and showed them in a bar graph (Fig. 3). Normalization was done according to the relation (x-min(x))/(max(x)-min(x)). As can be seen, both k and c parameters have a direct relationship with the frequency of aftershocks. Of course, there are exceptions to this. The earthquakes of Sarpolezahab and Mormori in the diagram related to k and the earthquakes of Bestak, Sisakht and Mormori in the diagram related to c are very different from other events.

### 3. Distribution Of Waiting Time

In this part, the waiting time distribution of the aftershock sequence as one of the temporal characteristics of aftershocks has been investigated. In some studies, a function in the form of gamma distribution has been used (e.g. Corral (2004); Hainzl et al. (2006)):

$$f\left(\tau \right)=C{\left(\frac{\tau }{{\tau }_{0}}\right)}^{\gamma -1} exp(-\frac{\tau }{{\tau }_{0}})$$

2

Where C and τ0 are constant of normalization and scaling parameter, respectively. Also γ characterizes the power-low decay. Michas and Vallianatos (2018) proposed the q-generalized gamma distribution for nonstationary earthquake time series:

$$f\left(\tau \right)=C{\left(\frac{\tau }{{\tau }_{0}}\right)}^{\gamma -1} . {\left[1+(1-q)\left(\frac{\tau }{{\tau }_{0}}\right)\right]}^{1/(1-q)}$$

3

Where the last term of this equation is the function of q-exponential (expq(\(\frac{\tau }{{\tau }_{0}}\))). This distribution has been widely used in temporal studies of seismic and aftershock sequences (e.g. Vallianatos et al. 2012, 2016, 2018; Vallianatous and Michas, 2020; Michas et al. 2022). To investigate the probability distribution p(τ), the histogram of waiting times τ is constructed, preferably in logarithmically spaced bins. Then, p(τ) is estimated by counting the number of τ that falls in each bin, further normalized with the bin width, and divided by the total number of counts (e.g., Corral (2004); Michas and Vallianatos (2018); Michas et al (2022)). The diagram of this distribution for all aftershock sequences is shown in Fig. 4. As can be seen, there is no significant difference in the diagrams of the 5 seismotectonic provinces. The average of all these lines was estimated and the best line that can be fitted was determined. The best fitted line has these values: C = 0.039 ± 0.004, τ0 = 9.2 ± 4.1, γ = 0.87 ± 0.08, q = 1.79 ± 0.23.