*4.2. Cluster Analysis*

Table 4 gives the chosen parameter set of the clustering algorithm for each dataset based on the analysis in Appendix C and a summary on the statistics of the detected clusters. In the CG area, we identified the largest number of seismic clusters (255) due to the increased detectability of micro-seismicity (low completeness magnitude threshold), however, they are short in size (*n*¯ = 18.28) and duration (*τ*¯ = 12.50) Conversely, the CII area is characterized by a small number of seismic clusters (45) but with large mean size (*n*¯ = 118.43) and duration (*τ*¯ = 54.60). The clustered seismicity is prevalent (75%), whereas in CG and NAS, the background component is more dominant than clustered seismicity with 64% and 56%, respectively (Table 4). In CG, this is explained by the lack of large main shocks during the study period and the occurrence of few moderate events, the largest number with *M* = 5.2 .

Concerning the CG area, the majority of the clusters are located on the western subarea where 22 out of 27 clusters with *N* ≥ 30 occurred. The main activity is located offshore between Aigion and Trizonia Island, but also north of the Psathopyrgos fault (Figure 4a). The eastern subarea comprises smaller clusters that are mainly concentrated offshore Xylokastro and Perachora faults, as well as near Itea Gulf (Figure 5a). The seismicity of CII is dominated by the two major main shock–aftershock sequences, each sequence comprising 2829 and 1396 events, respectively. Essentially, 4225 out of the 5221 clustered events belong to these sequences (Table 4). Furthermore, 45 clusters are detected in total with the main activity concentrated along the KTFZ (Figure 6a). The NAS area comprises 187 clusters, including both main shock–aftershock sequences and earthquake swarms (Table 4). Figure 7a shows that the main clustered activity is concentrated along the NAT and the sub-parallel branches, as well as in the southeastern subarea.


**Table 4.** Cluster statistics and the parameter set of the clustering algorithm for the three datasets. *Nclust* corresponds to the number of clustered events and *Nbg* to the background seismicity frequency. *τ*¯and *n*¯ are the mean duration in days and size of the clusters, respectively.

#### 4.2.1. Corinth Gulf Area

The western subarea of the Corinth Gulf is characterized by rich seismic activity, especially in 2013–2014, when 13 out of the 22 clusters with *N* ≥ 30 occurred. One of the major detected sequences is the 2013 Aigion swarm ( *C*6 in Figure 4b) which initiated on 21 May 2013 with a bulk of small events and several bursts associated to earthquakes with magnitudes ranging between 3.3–3.7 (Figure S3) [63,64]. Two distinct excitations followed ( *C*8 and *C*10 in Figure 4b) in accordance with the ones observed by Michas et al. [65]. The first cluster began on 7 July with some activity prior to the *M* = 3.7 event on 15 July 2013, and lasted until 27 August, 2013 (Figure S3). The second half of 2014 is also a well-studied period with intense seismic activity. Five clusters with *N* ≥ 30 are detected ( *C*15, *C*16, *C*18, *C*19 and *C*20) in the western subarea, including the offshore *M*4.8 earthquake on 7 November 2014, associated with *C*19 (Figure 4b), and the *M*4.6 event on 21 September 2014, associated with the earthquake swarm located between Nafpaktos and Psathopyrgos [66] ( *C*15 in Figure 4b). Persistent activity since 22 July 2014 is also observed offshore Aigion ( *C*16), close to the earthquake swarm, *C*15, which began on 7 November 2014 (Figure S5). In 2012, fewer clusters are observed, mostly during the first semester, with three clusters comprising *N* ≥ 30 events, *C*1, *C*2 and *C*3, and a plethora of smaller clusters (Figures 4b and S2). Between November 2013 and July 2014, the activity is sparse with three relatively large clusters, *C*11, *C*12 and *C*14 (Figures 4b and S4). Six more clusters with *N* ≥ 30 are observed until the end of 2017 ( *C*21, *C*22, *C*23, *C*24, *C*25, *C*27, Figure 4b).

**Figure 4.** (**a**) Spatial distribution of the centroids of the identified clusters for the western subarea of Corinth Gulf along with major faults (yellow lines). The size of the circles is proportional to the earthquake number in each cluster, whereas the duration is represented by the color scale. (**b**) Spatial distribution of the clusters with *N* ≥ 30 events. The index of each cluster is provided in the inset box.

The eastern subarea is characterized by more sparse activity. A major seismic sequence, Offshore Perachora ( *C*4 in Figure 5b), is detected, including two sub-sequences, the first initiated on 22 September and the second on 30 September 2012 (Figure S6). Two relatively large clusters, *C*13 and *C*17, are observed near Itea; the former lasted almost two weeks at the end of March, 2014, and the latter—almost three months between August and October 2014 (Figures 5 and S7).

**Figure 5.** (**a**) Spatial distribution of the centroids of the identified clusters for the eastern subarea of the Corinth Gulf along with the major faults (yellow lines). The size of the circles is proportional to the earthquake number in each cluster whereas the duration is represented by the color scale. (**b**) Spatial distribution of the clusters with *N* ≥ 30 events. The index of each cluster is provided in the inset box.

#### 4.2.2. Central Ionian Islands Area

The two main shocks of sequence *I*1 (Figure 6b) with *M* = 6.1 and *M* = 6.0 occupy the southern and the central part of the onshore area of Kefalonia Island. The 2014 Kefalonia earthquake sequence (*I*1 Figure 6) started on 19 January with the first main shock occurring on 26 January (*M* = 6.1), and aftershock activity extending over 35 km [58], part of which hosted the second main shock (*M* = 6.0) that occurred on 3 February and the compound aftershock activity. A sub-cluster is also detected offshore to the southwest of Kefalonia Island (*I*2 in Figure 6b) that is deployed concurrently with the main sequence (Figure S8). In addition, two distinct clusters, *I*3 and *I*4 (Figure 6b), are revealed, which occurred between November and December 2014 (Figure S9), across the edges of the double rupture. They might be triggered by the stress transfer of the main ruptures, indicating activation of adjacent fault segments. The seismic activity of cluster *I*5 (Figure 6b) retains the most interest because it is essentially two seismic excitations evolving at the same time. The first initiated in the Myrtos Gulf and the second offshore the south part of Kefalonia Island. It comprises 164 earthquakes in about 100 days (Figure S9). The activity of the *I*7 cluster (Figures 6b and S10) spreads along the western coastline of Lefkada and Kefalonia Islands, far beyond both sides of the 2015 Lefkada main rupture. To the south the aftershock activity is sparse, probably due to the large amount of stress released in the main rupture, revealing that the main slip is associated with a fault of about 17 km in length [59]. Apart from cluster *I*4, two additional clusters (*I*6 and *I*9 in Figure S9 and Figure S11, respectively) are detected in the area between Lefkada and Kefalonia, extending to about 15 km, which is considered as a transition zone encompassing step-over structures [58]. All of them relate to the E–W-oriented, parallel step-over faults, similar to the ones detected in the microseismicity cluster analysis between September 2016 and December 2019 in the study area [34].

**Figure 6.** (**a**) Spatial distribution of the centroids of the identified clusters for the area of Central Ionian Islands along with the trace of the Kefalonia Transform Fault Zone (yellow lines). The size of the circles is proportional to the earthquake number in each cluster, whereas the duration is represented by the color scale. (**b**) Spatial distribution of the clusters with *N* ≥ 30 events. The index of each cluster is given in the inset box.

#### 4.2.3. North Aegean Sea Area

The first seismic excitation with *N* ≥ 30 events (*N*1 in Figure 7b) is a sequence of interest since two moderate events (*M* = 5.2 and *M* = 5.3) occurred in 3 weeks, both producing their own aftershocks (Figure S12). The 2013, January 8 *M* = 5.8 North Aegean earthquake [60] along with its aftershock activity (cluster *N*3 in Figure 7b) is also detected. The aftershock activity is temporally divided into two clusters (Figure S13). The 24 May 2014 *M* = 6.9 Samothraki main shock was followed by aftershock activity confined to three major clusters (*N*4, *N*5, *N*6 in Figure 7b) and some secondary clusters with *N* ≥ 10 events (Figure S14), which are in accordance with the ones observed by Saltogianni et al. [61]. The seismic activity that took place near the Aegean coast of NW Turkey during January– October 2017 [67] is divided into three clusters with *N* ≥ 30 (*N*10, *N*11 and *N*14 in Figure 7b) and two minor clusters with 22 and 23 events, respectively (Figure S15). The strong main shock (*M* = 6.4) that occurred on the 12th of June 2017 offshore, south of the SE coast of Lesvos Island, along with its intense aftershock activity, is revealed and illustrated in Figure 7b (*N*12). Two major (*N* ≥ 30) secondary outbursts of clustered activity occurred concurrently on the west (*N*17) and east (*N*16) side of the sequence (Figure S16). A thorough analysis revealing multiple spatial clusters of the sequences is conducted by Papadimitriou et al. [62].

**Figure 7.** (**a**) Spatial distribution of the centroids of the identified clusters for the area of North Aegean Sea along with the trace of North Aegean Trough (yellow lines). The size of the circles is proportional to the earthquake number in each cluster, whereas the duration is represented by the color scale. (**b**) Spatial distribution of the clusters with *N* ≥ 30 events. The index of each cluster is given in the inset box.

#### *4.3. Regional Variability of Clustering Properties*

In this section, we investigate regional variations in the clustering behavior of the detected seismic sequences, in particular, on their productivity rates and on their temporal evolution that can differ among areas with distinct seismotectonic characteristics. Therefore, we adopt the temporal ETAS model that expresses two empirical relationships that characterize the temporal and size distribution of earthquakes, the normalized Omori–Utsu law, given by *g*(*t*) = *cp*−<sup>1</sup>(*p* − 1)(*t* − *ti* + *<sup>c</sup>*)−*p*, and the productivity law that is expressed by *N* = *k*(*Mi*) = *Kea*(*Mi*−*mc* ), where *N* is the number of triggered events by an earthquake of magnitude, *Mi*, *K* is a constant of proportionality, which depends on the number of triggered events per mainshock above the catalog cutoff, and *a* describes the impact of magnitude on the number of triggered events.

We compute the generic parameters for the 3 areas, CG, CII and NAS, by jointly inverting the ETAS parameter set *θ* = (*p*, *c*, *a*, *K*, *μ*) from the identified sequences with *N* ≥ 30 of each area. In particular, *LLi* = ∑*ni j*=1 log *λ*(*tj*) − *tend t*0 *λ*(*t*)*dt* denotes the log-likelihood of the *i*-th sequence for each area, namely, the logarithmic probability of observing *ni* events with occurrence times *tj*, *j* = 1, ... , *ni*, during the period of the sequence (*<sup>t</sup>*0, *tend*), and no other events between them. The intensity function, *λ<sup>t</sup>*, of the model is given by Equation (1), neglecting the spatial component. We then stack all the sequences of each area, compute their corresponding logarithmic probabilities *LLi*, and define as the common log-likelihood

$$LL = \sum\_{i=1}^{N^\*} LL\_{i\prime} \tag{2}$$

where *N*∗ is the number of sequences. The optimal inverted parameters are the ones that maximize Equation (2). The results of the ETAS parameter estimation for the three regions are shown in Table 5. There are 27 sequences in CG (Figures 4 and 5), 9 in CII (Figure 6) and 17 in NAS (Figure 7) from 2012 until 2019 with *N* ≥ 30 events, however, we removed cluster C26 from the computations, since it is located at the boundaries of the study area (Figure 5) with part of the aftershock data being omitted. For the maximization of the common log-likelihood *LL*, we implement an iterative procedure where at each step we update the model parameters by a random factor so *θnew k* = *θk* + *u*, for *k* = 1, ... , 5, then, we compute the corresponding *LLnew i* , *i* = 1, ... , *N*<sup>∗</sup>, log-likelihood values and store the new parameters under the condition *LLnew* > *LL* moving to the next iteration. After some iterations, the logarithm converges and the algorithm stops. Essentially, this is a grid-based procedure, since we use a large number of iterations.

**Table 5.** Generic ETAS parameter values for the three study areas. *N*∗ denotes the number of sequences with *N* ≥ 30.


The parameter *a* for CG (*a* = 0.82) is the lowest among the three areas, indicating the dominance of swarm activity presumably due to fluid flow in accordance with many relevant studies [65,68]. Low *a* values characterize areas with high heat flow [39], even though the estimated value can be underestimated due to magnitude incompleteness after the occurrence of the main shock or due to the existence of time-dependent background seismicity [69]. Conversely, in CII, the estimated value (*a* = 1.29) is relatively larger compared to the former region (*a* = 0.82), indicating the dominance of typical main shock– aftershock sequences. In the NAS area, a moderate value is acquired (*a* = 1.04), probably due to the co-existence of swarm activity and aftershock sequences. Another indicator for the existence of swarm activity in CG is the large value of the background seismicity (*μ* = 0.43) compared to NAS and CII. High values of the background rate can indicate the existence of aseismic loading transients [40]. LLenos et al. [70] observed increased values of the background component of the fitted ETAS model when it was applied to pre-swarm and swarm activity, respectively.

For the comparison of the productivity among the three areas, since they have different completeness magnitudes, we use the following relation,

$$N = k(M\_i)P(M \ge m\_c^\*) = K e^{a(M\_i - m\_c)} e^{-\beta(m\_c^\* - m\_c)},\tag{3}$$

which yields the number of earthquakes above magnitude *<sup>m</sup>*<sup>∗</sup>*c* , generated by a main shock of magnitude *Mi*. Figure 8 shows the number of direct triggered events, *<sup>K</sup>*(*Mi*), from an earthquake of magnitude *Mi*. We consider *<sup>m</sup>*<sup>∗</sup>*c* = 2.2, which is the maximum completeness magnitude among the three datasets. The exponent of the exponential magnitude distribution is expressed by *β* and is defined as *β* = ∑*Ni*=<sup>1</sup> *βi*/*N*<sup>∗</sup>, where *N*∗ is the number of clusters for each dataset and *βi* their corresponding exponent values. Concerning the distribution of aftershocks in time, the normalized Omori law distribution is used, given by *g*(*t*).

In CII, the seismic sequences seem to be more productive, as shown in Figure 8, with NAS and CG to exhibit smaller values. Combined with the higher background rate for the area of CG (*μ* = 0.43), we may say that a significant part of Corinth Gulf's sequences cannot be contributed to the triggering effect of mainshocks but different underlying mechanisms seem to play an important role. Conversely, in CII area, mainshock–aftershock sequences seem to dominate, generating a rich number of aftershocks (very low background rate, *μ* = 0.15, and high productivity of mother events).

**Figure 8.** (**a**) The number of events triggered by an earthquake of magnitude *Mi* at CG (blue), NAS (orange) and CII (brown), respectively. (**b**) The temporal distribution of triggered aftershocks.

#### *4.4. Sequence-Specific Clustering Properties*

Next, we estimate the *a*-values for the individual sequences of each area by maximizing *LL* as a function of *a*, *K* and the background rate *μ*, while keeping the rest of the parameters fixed for clusters with *N* < 80. In this way, we increase the robustness of the inversion procedure since there are sequences with few events. A similar procedure was followed by Page et al. [3] and Llenos and Michael [36] who demonstrated that fitting multiple parameters for a single sequence can be unstable and Hardebeck et al. [37] who implemented this method for the estimation of California aftershock parameters. We intend to investigate potential differences in the productivity (*<sup>a</sup>*, *K*) and the background rate, *μ*, among sequences of each area and their relation to different underlying triggering mechanisms. Productivity parameters *a* and *K* are correlated, so we enabled both to run during the iterative procedure. We also examine the value of the background rate among sequences since it can be also an indicator of aseismic transients in a region. Both parameters, *a* and *K*, are not influenced by *μ*, as we verified it by implementing the inversion procedure, also keeping parameters *a* and *K* fixed.

## 4.4.1. Corinth Gulf

In Table 6, the inverted parameters for the 26 clusters of dataset D1 with *N* ≥ 30 are given. We adopt the generic values of Omori law (*p* and *c*) for clusters with *N* < 80 to increase the stability of the inverted parameters. We observe relatively high background rates for most of the sequences and low *a* values, in particular, *a* < 1 for 10 out of the 26 clusters.

Concerning the 2013 Aigion earthquake swarm and its subsequent swarms (clusters *C*6, *C*9 and *C*10), we observe relatively low productivity values of the ETAS model (*a* = 1.31, 1.29, 1.13, Table 6) in accordance with studies suggesting pore-fluid pressure as the main triggering mechanism during the excitation [63]. Clusters *C*11 and *C*12 are part of the same swarm (Figure S4) that occurred offshore Psathopyrgos. Their relatively high background rates (*μ* = 1.34, 1.92) show that a significant part of the clustered seismicity cannot be explained by the empirical laws of the triggering part of the ETAS model. Cluster 14 is part of a major swarm that began on 8 June 2014 (Figure S4). Michas et al. [65] did not find high diffusion rates that are related to fluid pore pressure. However, the large background rate found in our study (*μ* = 4.86) and the low *a* value (*a* = 0.92) sugges<sup>t</sup> the existence of a non-typical mainshock–aftershock sequence, with more complex triggering mechanisms being responsible, such as aseismic creep. Similarly, the largest cluster in the dataset, the *C*15, located offshore Nafpaktos, is characterized by relatively high background rate (*μ* = 1.32) and low productivity (*a* = 1.38), more typical values for swarm activity. In contrast, clusters *C*18 and *C*19 that are more typical mainshock–aftershock sequences with a distinct in magnitude event in the initiation of the sequence (Figure S5), have low background rates (*μ* = 0.05, 0.76) and relatively high productivity rates (*a* = 1.77, 1.80). The two clusters near Itea Gulf show contradictory results, in particular, the first one, *C*13, is characterized by a high background rate (*μ* = 3.41), whereas the second, *C*17, which occurred four months later, exhibits a much smaller background value (*μ* = 0.35) more typical for mainshock–aftershock sequences. However, biases can exist in the inversion of the parameters for clusters with a small number of events, so we should be cautious with the inference.

**Table 6.** Details on the 26 clusters with *N* ≥ 30 events in CG area and the inverted ETAS parameters. The generic values of the Omori law, *p* and *c*, are adopted for clusters with *N* < 80.


#### 4.4.2. Central Ionian Islands

In Table 7, the inverted parameters for the nine clusters identified in the area of CII with *N* ≥ 30 events are given. We maintain fixed the Omori law parameters *p* and *c* (generic values) for clusters with *N* < 80. The estimated ETAS parameters of the sequence *I*1 are in accordance with the existence of a main shock–aftershock sequence described in Section 4.2. In particular, the background rate is relatively low (*μ* = 0.17), indicating that the seismicity is adequately described by the triggering part of the ETAS intensity function. The seismic activity of clusters *I*3 and *I*5 (shown by green and blue color in Figure S9, respectively) are characterized by relatively high background rates (*μ* = 0.99, 0.76, Table 7). The space-time evolution of the former indicates a rapid migration in the beginning of the sequence (Figure S9), whereas, for the latter, it is characterized by the smallest *K* value (*K* = 0.04) in the area although the *a* value is rather large. Taking into account the lack of distinct main shocks at the initiations of the sequences, they can be characterized as earthquake swarms, one of the few observed in an area which comprises mostly main shock–aftershock sequences. Concerning cluster *I*9, located in the transition zone between Lefkada and Kefalonia Islands, there is evidence for swarm activity due to the relatively high background seismicity rate (*μ* = 0.70). Ultimately, the major main shock–aftershock sequences in the area, *I*1, *I*7, have the highest *p* values (*p* = 1.42, 1.45), meaning that they are characterized by high aftershock decay in time.

**Table 7.** Details on the 9 clusters with *N* ≥ 30 events in CII area and the inverted ETAS parameters. The generic values of the Omori law, *p* and *c*, are adopted for clusters with *N* < 80.


#### 4.4.3. North Aegean Sea

In NAS area cluster *N*1, which consists of two moderate events ( *M* = 5.2 and *M* = 5.3) within a time period of 3 weeks, exhibits the lowest *a* value (*a* = 1.10) among the main detected clusters, which could be an indicator of fluid diffusion in the area (Table 8). Another case worth mentioning is the 24 May 2014, *M* = 6.9, Samothraki seismic sequence which is divided into three major clusters ( *N*4, *N*5, *N*6, in Figure 7). The estimated background rates of the three major clusters are relatively small (*μ* = 0.16, 0.60, 0.29), whereas the opposite holds for the scaling parameter, *a*, for the first two clusters (*a* = 1.82, 1.76). Concerning the seismic excitation that consists of clusters *N*10, *N*11 and *N*14, the relatively low productivity rates of the ETAS model (*a* = 1.31, 1.29, 1.13) and, conversely, the relatively high background rates for the first two, *N*10 and *N*11, clusters (*μ* = 1.00, 0.91) may indicate fluid intrusion. This observation is in accordance with the study of Mesimeri et al. [67] who derived high background rates after the estimation of the ETAS model to the empirically divided 5 sub-clusters of the primary seismic activity (January–March 2017). A fast-diminishing aftershock activity is observed for the main shock ( *M* = 6.4) that is located SE of Lesvos Island ( *N*12), which is translated into a high Omori exponent, *p* = 1.48. Additionally, low background rates characterize the three main clusters, *N*12, *N*16 and *N*17, indicating that they are probably related to tectonic and coseismic stress transfer from previous seismicity [62]. Worth mentioning are the remarkable high background rates for clusters *N*8 (*μ* = 1.76) and *N*9 (*μ* = 2.78), which could be an indicator for seismic activity driven by transient forces, however, the number of events is rather small and could have led to significant biases in the inversion of the parameters.
