Investigation of the Factors Controlling the Duration and Productivity of Aftershocks Following Strong Earthquakes in Greece

: Strong crustal earthquakes in Greece are typically followed by aftershocks, the properties of which are important factors in seismic hazard assessment. In order to examine the properties of earthquake sequences, we prepared an earthquake catalog comprising aftershock sequences with mainshocks of M w ≥ 5.5 from 1995 to 2021. Regional aftershock parameters were estimated to highlight variations in aftershock decay and productivity among regions with similar seismotectonic characteristics. A statistically based method of estimating aftershock duration and a metric of relative aftershock productivity to examine the variations among the different cases were employed. From the detailed analysis of the selected seismic sequences, we attempt to unravel the physical mechanisms behind deviations in aftershock duration and productivity and resolve the relative contribution of background seismicity, the Omori–Utsu law parameters and the mainshock faulting properties. From our analysis, the duration of aftershock sequences depends upon the rupture process of the mainshock, independently of its magnitude. The same applies to aftershock productivity, however, other tectonic setting (e.g., seismic coupling) or source-related (e.g., focal depth, stress drop) parameters also contribute. The estimated regional parameters of the aftershock rate models could be utilized as initial ones to forecast the aftershock occurrence rates at the early stage following a mainshock.


Introduction
The generation of aftershocks following a strong (e.g., M w ≥ 5.0) mainshock remains a critical key to understanding the fundamental characteristics of seismogenesis.It is widely accepted that they originate from the relaxation of stress concentrations due to the dynamic rupture caused by the mainshock or to the stress redistribution after its occurrence [1].The first empirical model describing the nature of aftershocks was introduced by Omori et al. [2], who claimed that the aftershock occurrence rate decays over time following a hyperbolic function, analogous to 1/t.Decades later, the observation that several earthquake sequences presented higher aftershock decay than predicted by the 1/t function led Utsu [3] to formulate a generalized version of the Omori law, most commonly known as the Modified Omori Law (MOL): where R(t) is the rate of aftershocks at time t, t is the time elapsed from the mainshock occurrence and K, c and p are parameters constrained in each seismic sequence.Among them, the c-value is the most ambiguous, indicating the temporal delay due to the shortterm aftershock incompleteness that is observed at the early stage of an aftershock sequence (e.g., [4]).The parameter p expresses the power law decay of the aftershock rate following a mainshock, thus, for sequences with higher estimated p-values, their aftershock rates decay faster with time than those obtaining a lower value.Unlike c and K, the p-value is independent of the cutoff magnitude, and is widely used to draw conclusions about the differences in the behavior of different seismic sequences.Among the MOL parameters, the K-value is considered to be more closely related to aftershock productivity.However, its value appears to depend on the total number of aftershocks and, consequently, on the completeness magnitude (M c ) [5,6].According to Kisslinger [4], the K-value is related to the earliest part of the sequence and originates from the fact that early aftershocks do not follow a constant decay rate, instead their rate is increased in the first hours after the mainshock occurrence and then starts to decrease.The frequency-magnitude distribution of earthquakes has been described and extensively analyzed through the empirical Gutenberg-Richter (G-R) law [7] which states that earthquake magnitudes are exponentially distributed as: where N is the total number of earthquakes with magnitudes greater than or equal to M, and a GR , b are constants with the latter being the slope, denoting the ratio of small-to large-magnitude earthquakes.The b-value is considered to vary among different regions reflecting the different stress regimes coming in play (e.g., [8,9]).Consequently, spatiotemporal variations in b-values have been and still remain a subject of discussion in numerous studies attempting to interpret its physical properties (e.g., [10]).Regarding Greece, studies have shown that the spatial variation in b is characterized by a general trend of gradually decreasing values from the outer Hellenic Arc towards NE Greece [11][12][13].
Another widely used empirical relationship in seismology is Båth's law, which states that within a typical mainshock-aftershock sequence, the difference between the magnitude of the mainshock (M m ) and the one of the largest aftershock (M a ) is approximately equal to 1.2 [14], independently of the mainshock magnitude.Since then, numerous studies have attempted to validate and extend this empirical law (e.g., [15][16][17].Drakatos and Latoussakis [18] estimated the relationship between M m and M a for earthquakes in Greece (up to 1997), finding an average value close to 1.0 but with high variance.
In the present study, we attempt to systematically review the aftershock properties of strong (M w ≥ 5.5) earthquakes in Greece by exploiting the basic models that describe aftershock processes with the ultimate goal being to reveal systematic patterns and properties, enabling us to highlight any anomalies and draw conclusions concerning the interconnection between the examined parameters.The aftershock productivity and the total duration of aftershock sequences are examined taking into account, among others, the magnitude and the faulting style of the mainshock, as well as region-specific characteristics.In Section 2, the dataset of aftershock sequences used and the seismotectonic frame of the study area are described.An overview of the applied methodologies to model aftershock rates and the total duration of seismic sequences follows, along with an objective measure of relative aftershock productivity.Results of the implemented methods are presented in Section 3. In Section 4, the outcomes of the analysis are discussed, focusing on the reasons controlling both regional and intersequence variability in aftershock productivity, as well as the effect of tectonic setting properties on aftershock duration.Finally, the main conclusions of the present study are summarized in Section 5.

Selection of Aftershock Sequences
In order to examine the aftershock properties of earthquake sequences in Greece, we used an earthquake catalog compiled from the bulletins of the Central Seismological Station of the Geophysics Department of the Aristotle University of Thessaloniki (geophysics.geo.auth.gr/ss/).Aiming to use accurate data with a satisfactory number of aftershocks, a mainshock magnitude threshold of M w ≥ 5.5 was set, covering the time period from 1995 to 2021.Spatiotemporal constraints were applied to select aftershocks associated with each mainshock.An earthquake with M w ≥ 5.5 is characterized as mainshock in cases where a higher magnitude one has not occurred within 3 fault lengths (3*L) from its epicenter, estimated using the well-established empirical relation of Wells and Coppersmith [19], 60 days before and 10 days after its occurrence.Then, aftershocks were associated with each mainshock if they occur within 3*L from its epicenter, a value already adopted for global aftershock forecasts [20], also taking into account that aftershock zone scaling is considered independent from the faulting type of the mainshock [21].From the initial catalog, a list of 73 aftershock sequences was obtained (Figure 1; Table 1) for the period 1995-2021.Isolated mainshocks (as they were identified with the above criteria, but without aftershocks) were excluded from the analysis.Additionally, the rupture type of each mainshock was classified according to their focal mechanism as provided by the Global Centroid Moment Tensor (GCMT) database.The classification was performed using the FMC program developed by Alvarez-Gomez [22] which is based on a hierarchical clustering analysis using the parameters provided by the focal mechanism solutions.nia parameters.Since then, Hardebeck et al. [23] updated these parameters by intro updated data and more detailed regionalization.The RJ89 model has also been s fully used on a global scale.A prime example is the work of Page et al. [20] who large global dataset and grouped tectonically uniform areas to determine the RJ89 eters for each tectonic regime.The RJ89 model has significant advantages in esti generic parameters but it comes with several assumptions, with the most importan that all aftershocks are triggered by the mainshock, which is not always true.As a in cases where large aftershocks trigger their subsequent aftershocks and in swarms, where similar-sized events occur, the model may perform inadequately.
Following the approach of Page et al. [20], we estimated mean regional afte parameters to optimally detect spatial variations in productivity and aftershock de may signify variations in the predominant stress field and tectonic loading process study area was divided into distinct tectonic regions (Figure 1), based on the gen nation scheme suggested by Vamvakaris et al. [13].Ultimately, six seismic zone defined based on common seismotectonic characteristics.Zones I and III are ge related to reverse or thrust faulting, Zones II and VI to strike-slip faulting and I normal faulting (Figure 2)

Modeling of the Aftershock Rates
Aftershock rate models constitute a powerful tool to describe aftershock sequences and most importantly to evaluate the predicted aftershock decay rate which can be of crucial significance in disaster management following a strong earthquake.The statistical framework proposed by Reasenberg and Jones [5] (RJ89), which has been utilized for decades by USGS as a forecast tool to estimate the probability of occurrence of strong aftershocks following M ≥ 5.0 mainshocks in California, was employed to model aftershock rates.The RJ89 model describes aftershock occurrence as a non-stationary Poisson process combining the G-R relationship and MOL.The rate of aftershocks at time t after a mainshock (M m ) is estimated as: where b is the slope of the frequency-magnitude distribution (G-R), a is the parameter representing aftershock productivity and p, c are the parameters of MOL denoting the temporal aftershock decay and the time shift before its beginning, respectively.M min is the lower magnitude threshold of the sequence.Reasenberg and Jones [5] fitted Equation (3) using a dataset of seismic sequences in California and provided the median parameters of a = −1.67,b = 0.91, p = 1.08 and c = 0.05 which were established as the "generic" California parameters.Since then, Hardebeck et al. [23] updated these parameters by introducing updated data and more detailed regionalization.The RJ89 model has also been successfully used on a global scale.A prime example is the work of Page et al. [20] who used a large global dataset and grouped tectonically uniform areas to determine the RJ89 parameters for each tectonic regime.The RJ89 model has significant advantages in estimating generic parameters but it comes with several assumptions, with the most important being that all aftershocks are triggered by the mainshock, which is not always true.As a result, in cases where large aftershocks trigger their subsequent aftershocks and in seismic swarms, where similar-sized events occur, the model may perform inadequately.Following the approach of Page et al. [20], we estimated mean regional aftershock parameters to optimally detect spatial variations in productivity and aftershock decay that may signify variations in the predominant stress field and tectonic loading processes.The study area was divided into distinct tectonic regions (Figure 1), based on the general zonation scheme suggested by Vamvakaris et al. [13].Ultimately, six seismic zones were defined based on common seismotectonic characteristics.Zones I and III are generally related to reverse or thrust faulting, Zones II and VI to strike-slip faulting and IV, V to normal faulting (Figure 2).
Zone I is characterized by the active compressional stress field due to the continental collision between the Apulian and the Aegean microplate [24,25].Active deformation is related to reverse faulting, with a general NW-SE strike following the coastline of NW Greece.Moving to the south (Zone II), this stress field is terminated by a major transform fault zone, namely, the Kefalonia Transform Fault Zone (KTFZ) which is characterized by dextral strike-slip sense of motion, often accompanied with thrust components [26][27][28].Mostly NE-SW to NNE-SSW dextral strike-slip faults dominate the area, being responsible for the highest seismicity among the other regions [29].The southernmost part of the KTFZ is connected to the NW edge of the Hellenic subduction front, the western part of which forms Zone III.The limits of this zone illustrate the amphitheatrical shape of the Hellenic Subduction zone, starting from a general NW-SE trend to the west, then bending (E-W) and finally obtaining a SW-NE direction reaching the SW coasts of Turkey (Figure 2).Thrust faulting dominates this zone, producing large-magnitude earthquakes, however, significant strike-slip earthquakes have also occurred in its easternmost part (e.g., [30]).Zone IV includes the extension stress field dominating the Corinth Gulf forming faults with a general E-W strike, formed under a N-S extensional deformation stress field.Normal faulting is predominant in this zone, although seismicity is often manifested as seismic swarms (e.g., [31,32]) with distinct triggering and driving processes, different than in the other normal faulting-related areas such as the ones enclosed in Zone V.This broad seismic zone includes (i) the area north of the subduction zone, with the presence of N-S to NNW-SSE striking normal faults, (ii) the relatively low-seismicity area of the central Aegean Sea, (iii) the eastern Aegean area with normal faults, generally striking E-W, which have hosted several strong events in the past [33], (iv) the Central Greece area, excluding the Corinth Gulf, and, finally, (v) the area of northern Greece within which the levels of seismicity vary spatially but contains several strong earthquakes hosted around specific normal fault zones [34], the latter being parts of the back arc area.Lastly, Zone VI encompasses the North Aegean region, where dextral strike-slip motion dominates at the western prolongation of the northern branch of the North Anatolian Fault along the NE-SW trending North Aegean Trough (Figure 2).

Aftershock Duration
Aftershock duration ( ) is defined as the time interval in which the overall seismicity is dominated by aftershock activity.It signifies, therefore, the time elapsed since the mainshock occurrence up to when the seismicity rate returns to the average background seismicity rate.Dieterich [35] suggested that under the rate-and-state friction law [36], an aftershock sequence decays hyperbolically as described in the Omori law, and returns to the background seismicity rate after an elapsed time Ta, which is independent of the mainshock magnitude and varies inversely with the fault stressing rate (Equation ( 4)): where A is a constitutive parameter of the time-dependent friction, σ is the normal stress on the fault and  is the stressing rate.Other researchers, however, consider this definition to be the "apparent aftershock duration" and the "true" aftershock duration to be the triggering time T which is the duration of the physical triggering process of a single event and can be significantly longer than the first one [37].
There are several approaches to assess  which are mainly differentiated by the selection of spatiotemporal criteria to estimate the seismicity rates.In this study, the aftershock areas were designed to include aftershocks both on the fault plane and in the surrounding area (3*L), which may be triggered by stress transfer.We followed the method proposed by Sebastiani et al. [38], which is based on the mean values difference significance test.Prior to the calculations, aftershock catalogs were compiled for each seismic sequence, starting from each mainshock and their corresponding  was estimated by applying the Goodness-of-Fit Method (GFT) proposed by Wiemer and Wyss [39] considering the 95% confidence interval level of residuals.In cases where this threshold was not reached, the 90% level was kept instead.

Aftershock Duration
Aftershock duration (T a ) is defined as the time interval in which the overall seismicity is dominated by aftershock activity.It signifies, therefore, the time elapsed since the mainshock occurrence up to when the seismicity rate returns to the average background seismicity rate.Dieterich [35] suggested that under the rate-and-state friction law [36], an aftershock sequence decays hyperbolically as described in the Omori law, and returns to the background seismicity rate after an elapsed time T a , which is independent of the mainshock magnitude and varies inversely with the fault stressing rate (Equation ( 4)): where A is a constitutive parameter of the time-dependent friction, σ is the normal stress on the fault and .
τ is the stressing rate.Other researchers, however, consider this definition to be the "apparent aftershock duration" and the "true" aftershock duration to be the triggering time T which is the duration of the physical triggering process of a single event and can be significantly longer than the first one [37].
There are several approaches to assess T a which are mainly differentiated by the selection of spatiotemporal criteria to estimate the seismicity rates.In this study, the aftershock areas were designed to include aftershocks both on the fault plane and in the surrounding area (3*L), which may be triggered by stress transfer.We followed the method proposed by Sebastiani et al. [38], which is based on the mean values difference significance test.Prior to the calculations, aftershock catalogs were compiled for each seismic sequence, starting from each mainshock and their corresponding M c was estimated by applying the Goodness-of-Fit Method (GFT) proposed by Wiemer and Wyss [39] considering the 95% confidence interval level of residuals.In cases where this threshold was not reached, the 90% level was kept instead.
For each seismic sequence, the following procedure was executed.A temporal sequence of daily rates of events with magnitudes above the M c threshold was computed for all catalog data (T all ), with the occurrence time of the mainshock considered as the origin time.T all was then split into two subsets of equal length.The first one contained the last T-values of daily rates, whereas the other was not stable and was formed each time by incorporating successive T-values from the first T all -T-values of the full catalog.Then, the mean value (µ) and variance (S 2 ) of each pair of subsets (that never intersect) were estimated and formed the sequence w t : The lengths of T s and T were calibrated for each seismic sequence separately, ensuring that the sample size T is adequate and T s is wide enough to include the end of the aftershock sequence.T a was then considered as the first time that w t is less than 3 (99.9%significance level) assuming that w t can be approximated by a Gaussian random variable [40].

Aftershock Productivity
It is well documented that the number N of aftershocks following a mainshock increases with the mainshock magnitude M m .This pattern is expressed via the productivity law (e.g., [18,[41][42][43]): where k depends on the number of aftershocks over the magnitude of completeness threshold and α controls the relative number of aftershocks activated corresponding to mainshock magnitude (M m ).Variations in α have been reported among different regions, for example, Felzer et al. [42] suggested a value of 1.0 for California, whereas Helmstetter [41] indicated a value of 0.8 for southern California.On the other hand, Drakatos and Latoussakis [18] calculated a value of 0.9 for the Greek territory.In this study, the productivity law was fitted through linear least squares regression after dividing the data into 0.1 M m bins and using the median number of aftershocks for each bin (Figure 3).The spatial selection of aftershocks of seismic sequences from 1995 to 2021 was the same as described in Section 2.1, but this time M c was estimated using the full catalog in order to obtain a uniform magnitude threshold (M c = 3.5).Three temporal windows (10, 60 and 120 days) were used to count aftershocks and test the sensitivity of the calculations.Variations between each temporal window were negligible, resulting in a value of α = 0.7.Based on this, each aftershock sequence was associated with a Relative Productivity (RP) value following the definition of Dascher-Cousineau et al. [43]: where N is the detected number of aftershocks of the seismic sequence and Ñ(M) is the expected number of aftershocks predicted from Equation ( 6) for the mainshock magnitude M m .The abovementioned approach was utilized to effectively compare the aftershock productivity of seismic sequences of different mainshock magnitudes and investigate the reasons behind their fluctuations.It is, therefore, crucial that the used measure of productivity is independent of mainshock magnitude to avoid linking variations in productivity with deviations in the average size of earthquakes [43].Then, having assigned a value of relative aftershock productivity to each seismic sequence, we attempt to analyze certain properties such as their faulting style and other source-related parameters to seek answers about the reason behind their increased or decreased aftershock production.

Aftershock Rates
Stacked fits of the mean aftershock rates were performed using different time intervals after the mainshock occurrence (10, 60, 120 days) to test the sensitivity of the results.Furthermore, the spatial constraint that was used (3*L) was deemed as the optimum solution to secure the entirety of selected aftershocks and reduce the background events as much as possible.The b-value of the G-R law for the catalog was estimated at 1.06 through the maximum likelihood method [44], with a  = 3.5 considering the 95% confidence interval level of residuals.The aftershock rate analysis was then performed through two approaches.First, stacked fits were performed assuming a uniform b-value (b = 1.06) for all zones (Figure 4) and then region-specific ones were estimated (Table S1).Throughout the analysis, the c-value was set to a uniform mean value (c = 0.1) that was obtained from sequence-specific MOL fits.
Figure 4 shows the variations in aftershock productivity between the selected seismic zones in the first 60 days after each mainshock occurrence.The 60-day stacked fit of all regions provide the "generic" mean parameters for Greece (a = −1.66,p = 0.86, c = 0.1) and can be used to estimate interval probabilities of strong aftershocks as a forecasting tool.Zones I and II are revealed to be the most productive ones, whereas Zones III, IV and VI are the most deficient in aftershocks.In particular, sequences belonging to Zones I and II appear to be 2-4 times more productive.As far as the estimated p-values are concerned, aftershock sequences located within Zones I and VI decay faster than the rest, whereas Zones III and IV exhibit the slowest decay rates.However, it is worth noting that Zone I encloses the smallest amount of data, making the computations more sensitive compared to other zones.Using mainshocks with  ≥ 5.0 for this zone leads to a drop in p-value (below 1.0) whereas the results regarding the other zones remain stable.In the Supplementary Materials, the results of 10-and 120-day stacked fits are provided (Table S2).In general, with a longer considered time window, the aftershock decay rate slightly increases, whereas the a-values show minimal adjustments.In all cases, the relative ranking of values regarding each zone remains the same.On the other hand, when using the region-specific b-values (Table S1), the productivity is scaled depending on their deviation from the generic b-value of the uniform catalog.This makes Zone VI more productive than previously, and the opposite applies to Zone I.The stacked fits corresponding to the

Aftershock Rates
Stacked fits of the mean aftershock rates were performed using different time intervals after the mainshock occurrence (10, 60, 120 days) to test the sensitivity of the results.Furthermore, the spatial constraint that was used (3*L) was deemed as the optimum solution to secure the entirety of selected aftershocks and reduce the background events as much as possible.The b-value of the G-R law for the catalog was estimated at 1.06 through the maximum likelihood method [44], with a M c = 3.5 considering the 95% confidence interval level of residuals.The aftershock rate analysis was then performed through two approaches.First, stacked fits were performed assuming a uniform b-value (b = 1.06) for all zones (Figure 4) and then region-specific ones were estimated (Table S1).Throughout the analysis, the c-value was set to a uniform mean value (c = 0.1) that was obtained from sequence-specific MOL fits.
Figure 4 shows the variations in aftershock productivity between the selected seismic zones in the first 60 days after each mainshock occurrence.The 60-day stacked fit of all regions provide the "generic" mean parameters for Greece (a = −1.66,p = 0.86, c = 0.1) and can be used to estimate interval probabilities of strong aftershocks as a forecasting tool.Zones I and II are revealed to be the most productive ones, whereas Zones III, IV and VI are the most deficient in aftershocks.In particular, sequences belonging to Zones I and II appear to be 2-4 times more productive.As far as the estimated p-values are concerned, aftershock sequences located within Zones I and VI decay faster than the rest, whereas Zones III and IV exhibit the slowest decay rates.However, it is worth noting that Zone I encloses the smallest amount of data, making the computations more sensitive compared to other zones.Using mainshocks with M w ≥ 5.0 for this zone leads to a drop in p-value (below 1.0) whereas the results regarding the other zones remain stable.In the Supplementary Materials, the results of 10-and 120-day stacked fits are provided (Table S2).In general, with a longer considered time window, the aftershock decay rate slightly increases, whereas the a-values show minimal adjustments.In all cases, the relative ranking of values regarding each zone remains the same.On the other hand, when using the region-specific b-values (Table S1), the productivity is scaled depending on their deviation from the generic b-value of the uniform catalog.This makes Zone VI more productive than previously, and the opposite applies to Zone I.The stacked fits corresponding to the "generic" b-value are herein presented as a more robust measure of highlighting their differences.
"generic" b-value are herein presented as a more robust measure of highlighting their differences.

Båth's Law
We further verified that the quantity ΔΜ ( - ), the difference in magnitude between the mainshock ( ) and the strongest aftershock ( ), is independent of  and the same also applies for ΔΤ, the time difference between the mainshock occurrence and the origin time of the largest aftershock.Average and median ΔΜ for a magnitude range of 5.5 to 7.0 are in the order of 1.2 (Figure 5), in full agreement with Båth's empirical law.Among the different faulting styles, strike-slip-related sequences tend to produce stronger, larger aftershocks (ΔΜ~0.9),whereas mainshocks associated with reverse faulting tend to have a bigger ΔΜ (~1.4).

Båth's Law
We further verified that the quantity ∆M (M m -M a ), the difference in magnitude between the mainshock (M m ) and the strongest aftershock (M a ), is independent of M m and the same also applies for ∆T, the time difference between the mainshock occurrence and the origin time of the largest aftershock.Average and median ∆M for a magnitude range of 5.5 to 7.0 are in the order of 1.2 (Figure 5), in full agreement with Båth's empirical law.Among the different faulting styles, strike-slip-related sequences tend to produce stronger, larger aftershocks (∆M~0.9),whereas mainshocks associated with reverse faulting tend to have a bigger ∆M (~1.4).
Regarding ∆T, the results indicate that in most cases (~56%), the largest aftershock occurred during the first 24 hours following the mainshock, which is in very good agreement with the findings of Drakopoulos [45] who reported that 60% of the examined cases of mainshocks in Greece are followed by their largest aftershock during the first day after their occurrence.Other studies using global data also acknowledge the fact that the largest aftershock is expected to occur rather early during a certain aftershock sequence than later [46].Regarding ΔΤ, the results indicate that in most cases (~56%), the largest aftersho occurred during the first 24 hours following the mainshock, which is in very good agr ment with the findings of Drakopoulos [45] who reported that 60% of the examined ca of mainshocks in Greece are followed by their largest aftershock during the first day af their occurrence.Other studies using global data also acknowledge the fact that the larg aftershock is expected to occur rather early during a certain aftershock sequence than la [46].

Aftershock Duration (𝑇 )
Aftershock duration ( ) is an indispensable component for the comprehensive u derstanding of aftershock patterns and crucial to seismic hazard assessment, given tha  is not properly constrained, aftershocks may be misinterpreted as background even leading to the overestimation of seismic hazard [47].We attempted to estimate  for selected seismic sequences and investigate the factors that influence their fluctuations.T final subset of reliably estimated  out of the entire dataset is limited to the cases th have an adequate number of aftershocks with magnitudes above the  , which var from case to case, and to those who pass the significance level test (Table S3).From o analysis, no correlation between  and mainshock size (for comparable magnitude bi was identified, as several researchers have also pointed out (e.g., [47][48][49]).A main goal this study was to determine whether the different tectonic settings influence  .Figu 6a,b reveal that seismic sequences associated with pure Reverse Faulting (RF) pres lower  -values compared to pure Strike-Slip (SS) and Normal (NF) ones.In particul the median value of  for reverse faulting sequences was found equal to 117 da whereas the other ones (395 days for strike-slip and 407 days for normal) exhibit larg and relatively similar durations. for mainshocks with an oblique rupture process oblique, RF oblique) are also shown separately in Figure 6b.Among them, two can distinguished from their counterparts.First, mainshock 7 (Table 1) exhibits a relativ short  compared to both NF and SS sequences.The area hosting this sequence belon to Zone VI (Figure 1) and more specifically to the area west of Chios Island, which is ch acterized by the transition from strike-slip to normal faulting displaying oblique faulti [50].Evidence of swarm-like activity is explicit in this area [50,51] and given that the co trolling factors in the duration of swarms (e.g., fluid diffusivity) are different from

Aftershock Duration (T a )
Aftershock duration (T a ) is an indispensable component for the comprehensive understanding of aftershock patterns and crucial to seismic hazard assessment, given that if T a is not properly constrained, aftershocks may be misinterpreted as background events, leading to the overestimation of seismic hazard [47].We attempted to estimate T a for the selected seismic sequences and investigate the factors that influence their fluctuations.The final subset of reliably estimated T a out of the entire dataset is limited to the cases that have an adequate number of aftershocks with magnitudes above the M c , which varies from case to case, and to those who pass the significance level test (Table S3).From our analysis, no correlation between T a and mainshock size (for comparable magnitude bins) was identified, as several researchers have also pointed out (e.g., [47][48][49]).A main goal of this study was to determine whether the different tectonic settings influence T a .Figure 6a,b reveal that seismic sequences associated with pure Reverse Faulting (RF) present lower T a -values compared to pure Strike-Slip (SS) and Normal (NF) ones.In particular, the median value of T a for reverse faulting sequences was found equal to 117 days, whereas the other ones (395 days for strike-slip and 407 days for normal) exhibit larger and relatively similar durations.T a for mainshocks with an oblique rupture process (SS oblique, RF oblique) are also shown separately in Figure 6b.Among them, two can be distinguished from their counterparts.First, mainshock 7 (Table 1) exhibits a relatively short T a compared to both NF and SS sequences.The area hosting this sequence belongs to Zone VI (Figure 1) and more specifically to the area west of Chios Island, which is characterized by the transition from strike-slip to normal faulting displaying oblique faulting [50].Evidence of swarm-like activity is explicit in this area [50,51] and given that the controlling factors in the duration of swarms (e.g., fluid diffusivity) are different from the mainshock-aftershock type sequences, that could be the reason behind its short T a .On the other hand, mainshock 63 (Table 1), namely, the Zakynthos 2018 sequence (e.g., [45]), exhibits a prolonged T a , especially compared to pure RF sequences.According to Karakostas et al. [52], the mainshock rupture (R-SS) led to the synchronous activation of multiple nearby fault segments through either static or dynamic stress transfer.As a result, a broad area was activated with an intense and long-lasting aftershock activity.exhibits a prolonged  , especially compared to pure RF sequences.According to Karakostas et al. [52], the mainshock rupture (R-SS) led to the synchronous activation of multiple nearby fault segments through either static or dynamic stress transfer.As a result, a broad area was activated with an intense and long-lasting aftershock activity.We then examined the dependence of  on the background seismicity rate and fault slip rate, both being proxies of the stressing rate.Background seismicity rates ( ) were calculated using seismicity up to 4 years prior to the mainshock occurrence and geodetic slip rate data were acquired from the literature (e.g., [53]).Fault slip rates were estimated according to the coupling ratio to account solely for its coseismic portion [54].Our analysis suggests that both  and slip rates exhibit weak to moderate inverse correlation with  (Figure 7).Therefore,  appears to be more closely related to stressing rate than to the mainshock magnitude, a view that is consistent with the rate-and-state friction law [35].We then examined the dependence of T a on the background seismicity rate and fault slip rate, both being proxies of the stressing rate.Background seismicity rates (R bg ) were calculated using seismicity up to 4 years prior to the mainshock occurrence and geodetic slip rate data were acquired from the literature (e.g., [53]).Fault slip rates were estimated according to the coupling ratio to account solely for its coseismic portion [54].Our analysis suggests that both R bg and slip rates exhibit weak to moderate inverse correlation with T a (Figure 7).Therefore, T a appears to be more closely related to stressing rate than to the mainshock magnitude, a view that is consistent with the rate-and-state friction law [35].
exhibits a prolonged  , especially compared to pure RF sequences.According to Karakostas et al. [52], the mainshock rupture (R-SS) led to the synchronous activation of multiple nearby fault segments through either static or dynamic stress transfer.As a result, a broad area was activated with an intense and long-lasting aftershock activity.We then examined the dependence of  on the background seismicity rate and fault slip rate, both being proxies of the stressing rate.Background seismicity rates ( ) were calculated using seismicity up to 4 years prior to the mainshock occurrence and geodetic slip rate data were acquired from the literature (e.g., [53]).Fault slip rates were estimated according to the coupling ratio to account solely for its coseismic portion [54].Our analysis suggests that both  and slip rates exhibit weak to moderate inverse correlation with  (Figure 7).Therefore,  appears to be more closely related to stressing rate than to the mainshock magnitude, a view that is consistent with the rate-and-state friction law [35].

Relative Productivity (RP)
Each seismic sequence was assigned with its corresponding relative productivity (RP) value (Table S3).Then, the dataset was subdivided according to the faulting style of the mainshock, as detailed in Table 1.The focal mechanism of the mainshock appears to have an influence on aftershock productivity, as demonstrated in Figure 8. Sequences originated from reverse faulting mainshocks appear to have 2-2.5 times fewer aftershocks than sequences related to normal or strike-slip faulting and the same applies to their T a (Figure 6b).Normal faulting mainshocks produce the greater number of aftershocks on average with the strike-slip faulting ones following closely.Considering the oblique focal mechanisms, it is observed that strike-slip motion accompanied with either thrust or a normal component enhances aftershock production and the same applies to reverse faulting with pure thrusting being at the bottom of average RP-values.On the contrary, more aftershocks are promoted by pure normal faulting than in instances where combined relative motion comes in play.To assess whether the tectonic setting is the only contributing factor to the relative productivity, we also examined the resulting RP-values according to the zonation scheme described in Section 2.2 to define certain geographic areas and draw conclusions.A synthesis of the parameters tested is shown in ascending order of median RP-values (Figure 8).The results indicate that the low RP-values regarding reverse faulting are mainly driven by mainshocks that occurred along the Hellenic Subduction zone (Zone III), especially its central and eastern part, taking into account that reverse faulting ones outside this zone (e.g., II and V) yield higher values.The primacy of the central Ionian Islands area (Zone II) is also evinced, being the most productive among all areas.If the tectonic setting was the only controlling factor on productivity, then strike-slip mainshocks outside this region should not have fewer aftershocks.However, the results indicate that strike-slip sequences within other regions (e.g., Zones V, VI) produce fewer aftershocks, implying that other parameters, such as the stressing rate, may decisively contribute.
(RP) value (Table S3).Then, the dataset was subdivided according to the faulting style of the mainshock, as detailed in Table 1.The focal mechanism of the mainshock appears to have an influence on aftershock productivity, as demonstrated in Figure 8. Sequences originated from reverse faulting mainshocks appear to have 2-2.5 times fewer aftershocks than sequences related to normal or strike-slip faulting and the same applies to their  (Figure 6b).Normal faulting mainshocks produce the greater number of aftershocks on average with the strike-slip faulting ones following closely.Considering the oblique focal mechanisms, it is observed that strike-slip motion accompanied with either thrust or a normal component enhances aftershock production and the same applies to reverse faulting with pure thrusting being at the bottom of average RP-values.On the contrary, more aftershocks are promoted by pure normal faulting than in instances where combined relative motion comes in play.To assess whether the tectonic setting is the only contributing factor to the relative productivity, we also examined the resulting RP-values according to the zonation scheme described in Section 2.2 to define certain geographic areas and draw conclusions.A synthesis of the parameters tested is shown in ascending order of median RP-values (Figure 8).The results indicate that the low RP-values regarding reverse faulting are mainly driven by mainshocks that occurred along the Hellenic Subduction zone (Zone III), especially its central and eastern part, taking into account that reverse faulting ones outside this zone (e.g., II and V) yield higher values.The primacy of the central Ionian Islands area (Zone II) is also evinced, being the most productive among all areas.If the tectonic setting was the only controlling factor on productivity, then strike-slip mainshocks outside this region should not have fewer aftershocks.However, the results indicate that strike-slip sequences within other regions (e.g., Zones V, VI) produce fewer aftershocks, implying that other parameters, such as the stressing rate, may decisively contribute.along with 11 cases of isolated mainshocks which were not assigned with an RP-value.From the spatial distribution of RP-values, regional variations are once again highlighted Figure 9 illustrates the RP ranking of the selected seismic sequences (Tables 1 and S3) along with 11 cases of isolated mainshocks which were not assigned with an RP-value.From the spatial distribution of RP-values, regional variations are once again highlighted as well as presumable anomalies in adjoining regions.First, the western part of the Hellenic Arc is more productive than its central and eastern parts, hosting intense seismic sequences such as the M w 6.8 Zakynthos 2018 (e.g., [52]), the M w 6.8 Methoni 2008 (e.g., [55]) or the remarkably prolific M w 5.7 Zakynthos 2006 seismic swarm (e.g., [56]) (mainshocks 63, 36 and 29, respectively; Table 1).Moving towards the eastern part, the vast majority of seismic sequences exhibit negative RP-values apart from the cases of M w 6.4 S. Crete 2009 [57] and the M w 6.1 Kasos 2015 [58] (mainshocks 42 and 55; Table 1).
Kozani 1995 sequence (e.g., [64]) (mainshock 1; Table 1) that occurred in an area not often visited by strong events, causing severe damage, which altered the seismic hazard of the area [65] and was among the events that led to the update of the Greek Building Code.A multi-segment rupture process led to the generation of numerous aftershocks, however, their almost synchronous activation served to limit  to customary levels (Table S3).The highest average RP is obtained in the central Ionian Islands area, irrespectively of mainshock magnitude.Aftershock sequences originating from strike-slip faulting outside this area seem to vary from case to case.As a case in point, the two sequences belonging to the N. Aegean region (zone VI), namely, the M w 6.1 Skyros 2001 (e.g., [59]) and the M w 6.9 N. Aegean 2014 (e.g., [60]) (events 26 and 53; Table 1), demonstrate a contrasting behavior.Despite their difference in moment release (by orders of magnitude), the latter is highly unproductive compared to the former, having almost half its T a (Table S3).Regarding normal faulting sequences, they appear to have a smaller spread in RP-values compared to the other faulting styles.For example, the recent strong events in the Eastern Aegean, namely, the M w 6.4 Lesvos 2017 (e.g., [61]), the M w 6.6 Kos 2017 (e.g., [62]) and the M w 7.0 Samos 2020 (e.g., [63]) sequences (mainshocks 59, 60 and 70, respectively; Table 1), have RP-values assigned close to zero, implying that their aftershock production was typical, given their respective order of magnitude.A rare antithetical example is the M w 6.5 Kozani 1995 sequence (e.g., [64]) (mainshock 1; Table 1) that occurred in an area not often visited by strong events, causing severe damage, which altered the seismic hazard of the area [65] and was among the events that led to the update of the Greek Building Code.A multi-segment rupture process led to the generation of numerous aftershocks, however, their almost synchronous activation served to limit T a to customary levels (Table S3).

Discussion
The application of the Reasenberg and Jones [5] (RJ89) model on a regional level uncovered several properties of aftershock sequences following strong (M w ≥ 5.5) mainshocks in Greece.The central Ionian Islands area (Zone II) appears to be almost four times more productive compared to the Corinth Gulf (Zone IV) which presents lower productivity values, with the Hellenic Subduction (Zone III) and the N. Aegean (Zone VI) areas following closely.These three areas exhibit similar a-values, however, the underlying mechanisms behind their low aftershock production are different.The Corinth Gulf area is characterized by the dominance of swarm-like activity [51,66] associated with the presence of fluid flow [67] which in turn has been linked to low aftershock productivity (e.g., [68]), whereas in the N. Aegean area, swarm-like and mainshock-aftershock activity coexist [51].The differences between the generation process of these two distinct manifestations of seismicity affect their expansion both in space and time.Swarms are primarily a result of fault weakening processes due to fluid intrusions onto a fault zone, whereas mainshock-aftershock activity is an outcome of stress changes caused by the fault slip (or possible afterslip) of the mainshock.As a result, the productivity and duration of swarms are highly variable, being affected by properties such as diffusivity and crustal permeability (e.g., [69]).The decreased aftershock productivity of the Hellenic Subduction zone area may be attributed to the focal depths of the mainshocks located within this zone, which are the largest among the others.Focal depth has been attested to be inversely related to productivity [43,70].This could be explained, to a certain extent, with the variation in the stress drop among relatively shallow and deep mainshocks and the stress homogeneity of the rupture area due to elevated pressure and temperature values at these focal depths.Deeper mainshocks tend to have higher stress drop than shallower mainshocks of comparable magnitude [71] which in turn is associated with lower aftershock production.Moreover, higher stress drop mainshocks are generally connected with smaller rupture dimensions, as a larger amount of energy is released per unit fault area.As a result, a more restricted area around the mainshock is exposed to stress transfer and triggering in the adjacent area.
Another property examined in this study is the aftershock duration (T a ) of the selected sequences, which, contrary to productivity, does not scale with mainshock magnitude [47,48].The most important outcome of our analysis is the dependence of T a on the style of faulting of the mainshock.Sequences whose mainshocks are originated from extensional tectonic settings (normal faulting) tend to have longer T a than compressional ones (reverse faulting).Transform fault zones (strike-slip) also exhibit long durations on average, however, they are more scattered, revealing region-specific dependencies.Valerio et al. [72] reached a similar conclusion using two different techniques, stating that sequences related to extensional setting are more abundant in aftershocks and last longer than those related to contractional regimes.A possible explanation for this can be ascribed to the different type of energy release dominating in each type of faulting.Thrust fault blocks move against gravity and quickly consume the energy to elevate the fault hanging wall upwards (fewer aftershocks, smaller T a ), whereas normal fault blocks move in favor of gravity, thus gravitational forces prevail over elastic stress release [72,73].As a result, when normal faults are activated, inertia is perpetuated by gravity (more aftershocks, higher T a ) until a gravitational equilibrium is reached [74].
Our data suggest a possible negative correlation between T a and background seismicity rate as well as the fault slip rate, both being proxies of the stressing rate.However, our dataset is relatively too small to draw definitive conclusions due to the fact that the reliable estimation of T a was possible only for a subset of the total data.According to the rate-andstate framework [35], as the stress step increases, the subsequent aftershock productivity becomes higher.In other words, an area with a low background seismicity rate R bg (related to a background stressing rate τ bg ) will need more time to respond in an abrupt stress step τ than a region with higher τ bg , where less time is needed for seismicity rate to return to its background level.However, we often find typical examples of highly productive sequences in areas of increased R bg .This has been addressed by Page and van der Elst [75] who argued that over a limited spatial scale, variations in stressing rate and background seismicity may reflect the local fault population in a given area.It is therefore projected that both background events and aftershocks are more likely to occur in areas with a higher concentration of faults.
Deviations in intersequence productivity among the zones are more intricate to interpret.For this reason, a measure of Relative Productivity (RP) was employed.RP is closely related to the a-value of the RJ89 model but is more efficient to assess variations in productivity on an earthquake-by-earthquake basis [43], provided that it is independent of the mainshock magnitude.The map illustrating RP-values in Greece for the study period (Figure 9) reveals once again the influence of the tectonic regime.The different behavior of these regimes, apart from the faulting style that comes in play, may be accredited to the level of seismic coupling characterizing these regions.The convergence between the African and the Eurasian plate along the Hellenic Arc occurs mostly aseismically, indicating only a fractional coupling across its interface (e.g., [54,76]).On the contrary, the Ionian Islands are considered a fully coupled area, with almost all seismic moment released seismically [77].Seismic coupling and aftershock productivity could well be interconnected, as the accumulated elastic energy (strain rate) depends on the locking degree of fault systems [78].It is therefore expected that areas with highly locked fault zones store a larger extent of stored elastic energy which in turn leads to more aftershocks.
Another source-related parameter that may explain variations in productivity and T a between aftershock sequences is the stress drop of the mainshock.Studies on the stress drop parameter in Greece (e.g., [79,80]) have revealed significant disparity between reverse faulting mainshocks compared to normal and strike-slip faulting ones, with the latter obtaining notably lower values.Our results imply an inverse relationship between productivity and stress drop, providing space for physical interpretation.For instance, two highlighted cases are the M w 6.5 Kozani 1995 and M6.5 Aigio 1995 sequences (e.g., [81]) (mainshocks 1 and 2, respectively; Table 1) which exhibit low stress drops [79] whilst having elevated RP-values (Figure 9) and a profuse duration compared to the rest.Our interpretation of the connection between stress drop and aftershock productivity lies in the available extent that is prone to failure after the mainshock occurrence.For a given mainshock magnitude, higher stress drop is akin to a smaller rupture area that subsequently leads to diminished aftershock production in contrast with lower stress drop values [71].

Conclusions
A comparative analysis of aftershock sequences following strong (M w ≥ 5.5) mainshocks in Greece was performed using two methods to assess aftershock productivity on a regional as well as on a case-by-case basis.We tested the efficacy of Båth's law and implemented a statistical method to estimate the total duration of the sequences.The central Ionian Islands area emerges as the one with the highest aftershock productivity whereas the Hellenic Trench exhibits the most inhibited aftershock activity.We find that aftershock productivity depends on the rupture process of the mainshock, however, other source-related parameters such as the focal depth and stress heterogeneity are contributing factors, leading to regional variability.The herein obtained results indicate that, irrespectively of the mainshock magnitude, the total duration of aftershock sequences is controlled by the rupture process of the mainshock.Extensional stress regimes exhibit longer durations compared to the other ones, presumably due to the different form of energy release coming in play.The examination of more source-related (i.e., by exploiting finite-fault solutions) or other parameters may constitute the scope of future studies in an effort to highlight variations amidst the aftershock sequences following strong (M w ≥ 5.5) mainshocks in Greece.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/geosciences12090328/s1.Table S1 lists the b-values calculated for each zone through the maximum likelihood estimate along with the confidence interval level of residuals.In Table S2, the results of the stacked fits of the RJ89 model using variable time windows (10, 60, 120 days) and mean (b = 1.06) or region-specific b-values are presented.Table S3 demonstrates

Figure 1 .
Figure 1.Map of the study area along with the selected seismic zones (Zone I: Continental Collision; Zone II: Central Ionian Islands; Zone III: Oceanic Subduction; Zone IV: Corinth Gulf; Zone V: Back Arc Area; Zone VI: North Aegean).Stars denote the identified mainshocks (M w ≥ 5.5) and circles the epicenters of earthquakes that occurred with M ≥ 2.5 from 1995 to 2021.Both are color-coded according to each seismic zone.

Geosciences 2022 , 20 Figure 2 .
Figure 2. The main geodynamic features of Greece.Red lines denote the active boundaries and the red arrows the relative motions (NAT: North Aegean Trough; KTFZ: Kefalonia Transform Fault Zone; RTF: Rhodes Transform Fault).

Figure 2 .
Figure 2. The main geodynamic features of Greece.Red lines denote the active boundaries and the red arrows the relative motions (NAT: North Aegean Trough; KTFZ: Kefalonia Transform Fault Zone; RTF: Rhodes Transform Fault).

Figure 3 .
Figure 3. Logarithm of the number of aftershocks of mainshocks with Mw ≥ 5.5 for 60 days after the mainshock occurrence for seismic sequences that occurred between 1995 and 2021.The productivity law is fitted through linear least squares regression using the median log number of aftershocks for each 0.1 magnitude bin.

Figure 3 .
Figure 3. Logarithm of the number of aftershocks of mainshocks with M w ≥ 5.5 for 60 days after the mainshock occurrence for seismic sequences that occurred between 1995 and 2021.The productivity law is fitted through linear least squares regression using the median log number of aftershocks for each 0.1 magnitude bin.

Figure 4 .
Figure 4. Stacked fits of the RJ89 model for seismic sequences located within each seismic zone of Figure 1.Red lines denote the maximum likelihood fit of the aftershock rates for M ≥ Mc.

Figure 4 .
Figure 4. Stacked fits of the RJ89 model for seismic sequences located within each seismic zone of Figure 1.Red lines denote the maximum likelihood fit of the aftershock rates for M ≥ M c .

Figure 5 .
Figure 5. Magnitude difference (ΔΜ) between the mainshock and its largest aftershock.Dash green and red lines correspond to the median and average ΔΜ, respectively.

Figure 5 .
Figure 5. Magnitude difference (∆M) between the mainshock and its largest aftershock.Dashed green and red lines correspond to the median and average ∆M, respectively.

Figure 6 .
Figure 6.(a) Examples of the methodology applied to estimate Ta for a strike-slip (green), normal (blue) and reverse (orange) faulting mainshock.The dashed horizontal line denotes the threshold value of wt = 3 corresponding to 99.9% significance level.(b) Results for all seismic sequences are presented classified by faulting style (NF: Normal Faulting; SS: Strike-Slip; RF: Reverse Faulting; SS Oblique, RF Oblique).

Figure 7 .
Figure 7. Aftershock duration (Ta, in days) as a function of (a) fault slip rate (blue squares) and (b) background seismicity rate (orange crosses).Data points are binned by mainshock magnitude.

Figure 6 .
Figure 6.(a) Examples of the methodology applied to estimate T a for a strike-slip (green), normal (blue) and reverse (orange) faulting mainshock.The dashed horizontal line denotes the threshold value of w t = 3 corresponding to 99.9% significance level.(b) Results for all seismic sequences are presented classified by faulting style (NF: Normal Faulting; SS: Strike-Slip; RF: Reverse Faulting; SS Oblique, RF Oblique).

Figure 6 .
Figure 6.(a) Examples of the methodology applied to estimate Ta for a strike-slip (green), normal (blue) and reverse (orange) faulting mainshock.The dashed horizontal line denotes the threshold value of wt = 3 corresponding to 99.9% significance level.(b) Results for all seismic sequences are presented classified by faulting style (NF: Normal Faulting; SS: Strike-Slip; RF: Reverse Faulting; SS Oblique, RF Oblique).

Figure 7 .
Figure 7. Aftershock duration (Ta, in days) as a function of (a) fault slip rate (blue squares) and (b) background seismicity rate (orange crosses).Data points are binned by mainshock magnitude.

Figure 7 .
Figure 7. Aftershock duration (Ta, in days) as a function of (a) fault slip rate (blue squares) and (b) background seismicity rate (orange crosses).Data points are binned by mainshock magnitude.

Figure 8 .
Figure 8. Relative aftershock productivity (RP) in relation to the style of faulting classification (RF, SS, NF, RF oblique, SS oblique, NF oblique) and according to the zonation scheme of this study (Zones I-VI).Red diamonds represent the median RP-value, and the caps denote the extent of calculated RP-values.

Figure 9
Figure 9 illustrates the RP ranking of the selected seismic sequences (Tables 1and S3)along with 11 cases of isolated mainshocks which were not assigned with an RP-value.From the spatial distribution of RP-values, regional variations are once again highlighted

Figure 8 .
Figure 8. Relative aftershock productivity (RP) in relation to the style of faulting classification (RF, SS, NF, RF oblique, SS oblique, NF oblique) and according to the zonation scheme of this study (Zones I-VI).Red diamonds represent the median RP-value, and the caps denote the extent of calculated RP-values.

Figure 9 .
Figure 9. Map of mainshocks (M w ≥ 5.5), color-coded according to the relative productivity value of their respective sequences.Gray circles indicate isolated mainshocks (without aftershocks).The serial numbers of each mainshock correspond to those of Table1.