## 3 Earthquake Data  
We considered three areas in the region of Greece (Figure 2a), which consist of distinctive seismotectonic units. The selection of the three study areas is based on criteria related to the homogeneity of the type of faulting, the comparatively intense continuous seismicity and the existence of seismic excitations during the study period. The first area, Corinth Gulf (CG) (Figure 2b), is undergoing high extensional deformation rates. The seismicity is mainly associated with eight major faults that bound the rift to the south and dip to the north [48]. The area of central Ionian Islands (CII) (Figure 2c) is characterized by the highest moment rate in the Mediterranean region. Its main seismotectonic feature is the Kefalonia  
Figure 1: (**a**) Cumulative number of background events for each algorithm, the initial ETAS catalog with black color and the mother events of the ETAS synthetic catalog with the black dotted line. The space-time evolution (**b**) of the initial ETAS catalog and of the background seismicity for the four best algorithms, (**c**) MAP-DBSCAN27 (\(J_{2}=0.647\)), (**d**) NN (\(J_{2}=0.727\)), (**e**) RB2 (\(J_{2}=0.630\)), (**f**) GK3 (\(J_{2}=0.676\)). Colored ellipses stand for large gaps and significant concentration of events.  
Transform Fault Zone (KTFZ) which extends for more than 100 km along the western coastlines of Lefkada and Kefalonia Islands, comprising two distinct main branches, the Lefkada and Kefalonia faults, respectively. Right lateral strike slip motion with a minor thrust component is the dominant faulting type [49; 50]. The third area is located in the North Aegean Sea (NAS) (Figure 2d) and is dominated by dextral strike-slip faulting, along the North Aegean Trough (NAT) and its parallel branches [51], as a consequence of the westward propagation of the North Anatolian Fault into the Aegean [52]. The driving mechanism of the active deformation in the Aegean region is the subduction of the oceanic lithosphere of the Eastern Mediterranean under the continental Aegean microplate, forming the Hellenic Subduction Zone and the extensional back arc Aegean area due to the slab rollback [53].  
For the investigation of the clustering properties in the three areas, we considered earthquake datasets from the regional catalog of the Geophysics Department of the Aristotle University of Thessaloniki [54], compiled with the recordings of the Hellenic Unified Seismological Network (HUSN). The earthquake catalogs of CG, CII and NAS, which we denote henceforth as D1, D2 and D3, include 25,595, 24,085 and 21,139 events, respectively, occurring between 2012 and 2019. For the determination of their completeness magnitude, we implemented the Goodness-of-Fit (GFT) method [55], assuming that earthquakes follow the Gutenberg-Richter (GR) law, \(logN=a-bM\). In particular, the differences between the observed and the synthetic frequency-magnitude distributions are computed for increasing magnitude bins as threshold values. The completeness magnitude is defined as the first magnitude bin at which the difference falls under the 5% residual (Figure S1). Figure S1-c show the residuals for the three datasets and Figure S1d-f present the GR law for the corresponding complete datasets. The \(b\)-value is calculated by means of the maximum likelihood method proposed by [56] and found equal to \(b=0.97,\ 0.88,