**5. Discussion**

The consistency and efficiency of the MAP-DBSCAN method is examined on a simulated earthquake catalog of 18 years that produces the main features of seismicity in the region of Greece. In particular, we showed that our method is able to identify the connections among the events generated by a spatiotemporal ETAS model, as well as the mother events that initiated each cluster. The knowledge of the links among the events enabled the comparison of the method with some well known clustering algorithms, like the Gardner and Knopoff, the Reasenberg and the Nearest-Neighbor, by the use of the Jaccard index. This is a tool for measuring the overlap between the original partition of events into clusters and background seismicity, and the estimated one after the implementation of each clustering method. The results show that MAP-DBSCAN method is very competitive and in most cases outperforms the tested algorithms. The NN achieves the best reconstruction

of the clusters (Table 2), which is probably related to the similarity of its metric with the ETAS metric that is used for the generation of the seismicity. The window-based method overestimates the clustered seismicity in accordance with work by Peresan and Gentili [11], whereas the Reasenberg link-based method seems to overestimate the background events (Figure 1).

The advantage of using the MAP model lies in its efficiency in capturing the changes in seismicity rate, independently of the mechanisms responsible for each seismic sequence. Furthermore, in case of non-stationary background seismicity, the MAP model can approximate the different phases by embedding multiple states into the Markov process *Jt*, i.e., distinct occurrence rates, and adopting a multiple rate threshold alternating according to the phase of the process each time. In this way, although it is more complicated, we can model both the non-stationary background seismicity and the triggered events without declustering the earthquake catalog [33]. The DBSCAN algorithm does not assume any specific spatial distribution of earthquakes and settles them into groups based solely on their spatial density.

We applied our method to three seismic zones in Greece during 2012–2019, identifying the major seismic sequences and a plethora of smaller ones. The rich seismic activity during 2013–2014 in the western subarea of the Corinth Gulf is detected in detail, a nontrivial issue, especially for the area between Nafpaktos-Psathopyrgos and offshore Aigion, where multiple excitations occurred in close proximity and within short periods (Figures 4, S4 and S5). Seismicity in the eastern subarea of the Corinth Gulf is found to be more sparse with few major clusters located near Itea Gulf (Figures 5 and S7) and offshore Perachora and Xylokastro (Figure 4 and S6). On the contrary, seismicity in the Central Ionian Islands is dominated by the two major main shock–aftershock sequences associated with the 2014 Kefalonia and the 2015 Lefkada seismic sequences (Figure 6). Together they comprise the 81% of the clustered seismicity in this area. Many large clusters are identified in the North Aegean Sea area that includes both main shock–aftershock sequences and earthquake swarms.

We investigated the properties of clustering seismicity among the three study areas with the use of the ETAS model. The results indicate that there are differences in aftershock productivity rates between Corinth Gulf, Central Ionian Islands and North Aegean Sea, showing that productivity can vary regionally. As showed by Page et al. [3] and LLenos and Michael [36] adopting the regional variations of productivity can produce a significant gain on aftershock forecasts. In the Central Ionian Islands, main shock–aftershock sequences seem to be more productive with the North Aegean Sea and the Corinth Gulf to follow (Figure 8). The sequences in the Corinth Gulf in particular are characterized by the highest background rate among the three areas (Table 5), meaning that a significant portion of clustered seismicity is not caused by the triggering of a main shock coseismic slip, but by the contribution of different triggering mechanisms. Many studies have focused on this area, suggesting pore-pressure changes due to fluid migration and aseismic creep as possible triggering mechanisms for the clustered seismicity [57,71]. In the North Aegean Sea, the swarm activity coexists with aftershock sequences, implying that for forecasting purposes, a finer regionalization might be more appropriate.

We also investigated potential differences in the productivity and the background rates among sequences of each region and their relation to different underlying triggering mechanisms. Results show that the high background seismicity (*μ*) and low productivity (*a*) values of the ETAS model are related to earthquake swarm activity triggered by fluid pore-pressure changes, such as the 2013 Aigion swarm (clusters *C*6, *C*9 and *C*10, Table 6, Figures 4 and S3) in Corinth Gulf [63] and the 2017 Tuzla earthquake swarm (clusters *N*10, *N*11 and *N*14, Table 8, Figures 7 and S15) in North Aegean Sea [67]. This is in accordance with studies suggesting the dependence of low productivity values to the existence of fluids [39,69]. In general, 18 out of 26 clusters in Corinth Gulf have background rates *μ* > 1 and low productivity values (11 out of 26 with *a* < 1), whereas in the Central Ionian Islands, where main shock–aftershock sequences dominate, we observe very low background rates of the ETAS model (all with *μ* < 1) and relatively high productivity values. In the North Aegean Sea area, we cannot observe a clear pattern, however, the majority of the detected clusters are characterized by low background rates and relatively high productivity, suggesting the dominance of typical main shock–aftershock sequences.


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