**4. Results**

#### *4.1. Triggered and Background Seismicity Separation*

We fit MAP models with between two and seven states for each earthquake subcatalog, computationally a very demanding process as the number of states increases, especially for large datasets such as D3 with 13043 events. According to the BIC values, six, seven and again six states are sufficient to approximate the temporal distribution of earthquakes for the D1, D2 and D3 datasets, respectively. Next, we evaluate the transitions among the hidden states of the models. Firstly, the state probabilities *pi*(*t*) = *p*(*Jt* = *i*) for *i* = 1, ... , *K* are estimated with the use of the forward and backward vectors given in Equation (A1), and then we assign as state of the hidden process *Jt*, the one with the highest probability, i.e., argmax0≤*i*≤*K pi*(*t*), with *pi*(*t*) = *pi*(*tk*) for *tk* ≤ *t* < *tk*+1. Each state *i* corresponds to an occurrence rate, *λi*, therefore, by evaluating the transitions among the states of the model, we detect change-points in the seismicity rate. Figure 3a–c illustrate the transitions among the states for the datasets D1, D2 and D3, respectively. The colored box at each temporal

interval *tk* ≤ *t* < *tk*+<sup>1</sup> indicates the state with the maximum probability at the current time and the legend contains its corresponding occurrence rate.

**Figure 3.** Most probable path of the hidden states of the model along with the daily frequency of events (gray vertical lines) with (**a**) *M* ≥ 1.5 for D1, (**b**) *M* ≥ 2.2 for D2 and (**c**) *M* ≥ 2.1 for D3 datasets, respectively. Each color is assigned to a different state *i* with seismicity rate *λi*. Inset magnifies the transitions among the states, which are otherwise difficult to visualize due to the short sojourn times compared to the study period. The rate threshold, *λthr*, is set equal to *λ*2 = 3.01, *λ*1 = 0.58, *λ*1 = 1.53 for the D1, D2 and D3 datasets, respectively.

The temporal patterns of dataset D1 indicate the dominance of state 2 (yellow color, Figure 3a) with occurrence rate *λ*2 = 3.01 events/day for almost the entire period. Nevertheless, there is a slight decrease in the occurrence of earthquakes (*λ*1 = 1.23 events/day) in the second part of the catalog, starting from 02/2016 with transitions to state 1 (red color, Figure 3a) until almost the end of the catalog in 12/2019. This is probably related to the lack of seismic sequences during the last part of the study period compared to the previous intense seismic activity especially during the period 2013–2014 in the western subarea of the CG [57]. The rate threshold is set equal to *λthr* = *λ*2, which we consider as the background rate during the study period.

The seismicity of the CII area is dominated by the two major sequences during the study period, the 2014 Kefalonia doublet (*Mw*6.1 and *Mw*6.0) [58] and the *Mw*6.5 2015 Lefkada earthquake sequence [59]. States 7 (brown), 6 (dark cyan), 4 (dark blue), 3 (orange) and 2 (yellow) in Figure 3b are clearly associated with the aftershock evolution of the two sequences—essentially, they approximate the Omori temporal distribution. Background

seismicity is described by state 1 (red) with occurrence rate *λ*1 = 0.58 events/day, which we set as rate threshold for the primary classification of the clusters.

Finally, dataset D3 also contains some major sequences, the 2013 *Mw*5.8 [60], the 2014 *Mw*6.9 Samothraki [61] and the 2017 *Mw*6.4 Lesvos earthquake sequences [62], whose aftershock temporal distribution is approximated by states 6 (dark cyan), 4 (dark blue), 3 (orange) and 2 (yellow) of the model (Figure 3c). The rate threshold value is set equal to *λthr* = *λ*1.

In general, we observe significant variations in the temporal evolution of the seismic excitations between the CG and the CII and NAS areas. Figure 3 illustrates that the daily frequency of events during seismic sequences in CII and NAS is decreasing in time, typical of mainshock–aftershock sequences, whereas in CG we observe large fluctuations in the daily frequency, common for earthquake swarms, as in 2014 when multiple seismic excitations occurred in the western subarea of CG.

Consecutive events above the rate threshold *λthr* are classified into groups which we call potential clusters, and then, we test four different sets of temporal constraints, (*<sup>T</sup>*, *dt*), to the three datasets. Potential clusters within a temporal interval *T* are merged into one and events that occurred in ±*dt* time from the potential cluster are also included. Next, the DBSCAN algorithm is implemented to the merged clusters in order to separate them based on their spatial density. The minimum number of neighbors for the determination of a cluster is set equal to 4 ( *Npts* = 4) for avoiding insignificant cases with fewer events. This is an appropriate choice for two-dimensional data according to Ester et al. [35]. For the determination of the distance threshold, , we computed the *k*-distances which is a procedure proposed by Ester et al. [35], which is commonly used to constrain the distance threshold [4]. In Appendix C we provide more details on the choice of the parameters and how they affect the spatio-temporal evolution of background seismicity.
