## 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  
\begin{table}
\begin{tabular}{l c c c c c c c c c} \hline
**ID** & \(T_{in}\) & \(T_{end}\) & \(N\) & \(p\) & \(c\) & \(b\) & \(a\) & \(K\) & \(\mu\) & \(M_{max}\) \\ \hline I1 & 19/1/14 & 16/9/14 & 2829 & 1.42 & 0.  
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 [PERSON] and [PERSON] [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 \(I_{t}\), 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 [PERSON] et al. [3] and [PERSON] and [PERSON] [36] adopting the regional variations of productivity can produce a significant gain on aftershock forecasts.