**Appendix B**

#### *Appendix B.1. Gardner and Knopoff Window-Based Method*

The procedure introduced by Gardner and Knopoff [10] for the detection of aftershocks is based on specific magnitude dependent space-time windows. It is known as the windowbased method, and it is one of the simplest forms of aftershock identification. For each earthquake with magnitude *M*, the subsequent events are assigned as aftershocks if they occur within a temporal window *t*(*M*) and a spatial interval *d*(*M*), respectively. Foreshocks are treated as aftershocks when a larger earthquake occurs later in the sequence. The event is considered as an aftershock and the algorithm is repeated based on the largest magnitude.

We give in Equation (A2) the functional form of the spatial and temporal windows suggested in Gardner and Knopoff [10], which are denoted as GK1. Additionally, in Equations (A3) and (A4) we present alternative window parameter settings that can be found in van Stiphout et al. [75]. We denote them as GK2 and GK3, respectively.

$$d = 10^{0.1238 \ast M + 0.983} \text{ (km)} \text{ and } t = \begin{cases} 10^{0.032 \ast M + 2.7389} & M \ge 6.5\\ 10^{0.5409 \ast M - 0.547} & M < 6.5 \end{cases} \text{ days} \right\} \tag{A2}$$

$$d = e^{1.77 + \sqrt{0.037 + 1.02 \ast M}} \text{ (km) and } t = \begin{cases} 10^{2.8 + 0.024 \ast M} & M \ge 6.5\\ e^{-3.95 + \sqrt{0.62 + 17.32 \ast M}} & M < 6.5 \end{cases} \text{ days} \right\} \tag{A3}$$

$$d = e^{-1.024 + 0.804 \ast M} \text{ (km)} \text{ and } t = e^{-2.87 + 1.235 \ast M} \text{ days} \tag{A4}$$

#### *Appendix B.2. Reasenberg Linked-Based Method*

In this method, an interaction zone among earthquakes is assumed that is modeled based on estimates of the stress redistribution for the spatial extent and on a probabilistic model, the Omori law, for the temporal extent, respectively. Any earthquake that occurs within the interaction zone of a prior earthquake is considered an aftershock and is included in the cluster. For the Reasenberg algorithm, we used the ZMAP tool [76] and we adopted 3 different sets of parameters given in Table A1. The parameters *τmin* and *τmax* correspond to the minimum and maximum elapsed time since the last event, in order to observe the next correlated earthquake at a certain probability, *p*1. Additionally, *xmef f* denotes the minimum magnitude threshold for the earthquake catalog, whose value in the clusters is raised by a factor *xk* of the largest earthquake within. Finally, the parameter *rf act* corresponds to the radii we adopt to consider linking a new event with the cluster.

**Table A1.** Input parameters for the Reasenberg clustering algorithm. The first row corresponds to the standard parameter set [12].


#### *Appendix B.3. Nearest-Neighbor Method*

The approach is based on the space-time-magnitude distance metric among two earthquakes given by Baiesi and Paczuski [22]:

$$
\eta\_{\rm ij} = (t\_j - t\_i) r\_{\rm ij}^{d\_f} \mathbf{1} \mathbf{0}^{-bm\_i} \,\tag{A5}
$$

where *rij* is the epicentral distance between events *i* and *j*, *df* is the spatial fractal dimension and *b* is the component of the Gutenberg–Richter distribution. Each event *j* is connected to its nearest neighbor *i* = *argmini*:*tj*>*tiηij* if their distance, *ηj*, is lower than a predefined threshold *η*0. The earthquake catalog is then partitioned on distinct clusters, each containing at least one event. For the selection of the threshold value, *η*0, the logarithm of the nearest neighbor distance *η*∗ = {*ηj*}*j*=1,...,*<sup>N</sup>* is considered, where *N* the number of events. It follows an 1D Gaussian distribution with two components, which is essentially a mixture model of two Gaussian densities with parameters *<sup>N</sup>*(*μ*1, *<sup>σ</sup>*1), *<sup>N</sup>*(*μ*2, *<sup>σ</sup>*2) and *a*1, *a*2 weights, respectively. Then, the intersection of the two functional forms gives the threshold value.

There are only two free parameters, the fractal dimension *df* and *b* value, which are considered equal to *df* = 1.51 and *b* = 1.0, respectively. The logarithm of the separation distance is equal to log *η*0 = −5.04, based on the intersection of the two modes in the 1D density distribution of distances (Figure A1).

**Figure A1.** Distribution of the NN distances among all pairs of earthquakes of the ETAS synthetic catalog. (**a**) 1D density distribution of log *η*, with estimated Gaussian densities for clustered (yellow) and background (orange) components. (**b**) 2D joint distribution of rescaled space and time distances.

#### *Appendix B.4. MAP-DBSCAN Method*

A MAP with 7 states is chosen based on BIC, the rate threshold is set to *λthr* = *λ*1, and different temporal windows are tested for merging the potential clusters. Finally, the DBSCAN algorithm is implemented for 5 different distance thresholds (). The minimum number of events is set to *Npts* = 2 for a better comparison with the other methods where clusters with at least 2 events can be defined. In Table A2 we present details on the parameter tuning.


**Table A2.** The 30 different parameter sets used for the detection of the clusters.

The method seems rather insensitive to the parameter selection. In particular, Figure A2 presents the Jaccard index values that describe the efficiency of the method to correctly reconstruct the initial clusters (*J*1) as well as to identify the single events (*J*2). We observe that the Jaccard index values are quite stable with small fluctuations, apart from the smallest upper-distance cutoff, = 2.5 km, which seems inadequate to capture the spatial correlations among the events. Furthermore, the contribution of the temporal constraints to the clustering procedure seems negligible, with the exception of the two peaks for PS12 and PS27. This is an indicator that the MAP model has already achieved a sufficient separation between background and triggered seismicity based on the embedded multiple rates of the model.

**Figure A2.** The Jaccard index values, *J*1 with blue and *J*2 with orange color, respectively, for all the input parameters of the MAP-DBSCAN method.
