Identification and Temporal Characteristics of Earthquake Clusters in Selected Areas in Greece

Abstract

The efficiency of earthquake clustering investigation is improved as we gain access to larger datasets due to the increase of earthquake detectability. We aim to demonstrate the robustness of a new clustering method, MAP-DBSCAN, and to present a comprehensive analysis of the clustering properties in three major seismic zones of Greece during 2012–2019. A time-dependent stochastic point model, the Markovian Arrival Process (MAP), is implemented for the detection of change-points in the seismicity rate and subsequently, a density-based clustering algorithm, DBSCAN, is used for grouping the events into spatiotemporal clusters. The two-step clustering procedure, MAP-DBSCAN, is compared with other existing methods (Gardner-Knopoff, Reasenberg, Nearest-Neighbor) on a simulated earthquake catalog and is proven highly competitive as in most cases outperforms the tested algorithms. Next, the earthquake clusters in the three areas are detected and the regional variability of their productivity rates is investigated based on the generic estimates of the Epidemic Type Aftershock Sequence (ETAS) model. The seismicity in the seismic zone of Corinth Gulf is characterized by low aftershock productivity and high background rates, indicating the dominance of swarm activity, whereas in Central Ionian Islands seismic zone where main shock-aftershock sequences dominate, the aftershock productivity rates are higher. The productivity in the seismic zone of North Aegean Sea vary significantly among clusters probably due to the co-existence of swarm activity and aftershock sequences. We believe that incorporating regional variations of the productivity into forecasting models, such as the ETAS model, it might improve operational earthquake forecasting.

1. Introduction

Earthquake clustering is an essential property of seismicity and is manifested as the concentration of earthquakes in space and time. Due to the improvement of seismic monitoring worldwide and the development of new powerful algorithms for earthquake detectability [1] additional information is available, which is crucial for reliable regional estimates of aftershock forecasting probabilities [2,3] and the determination of faulting geometry [4,5], among others. In addition to the necessity of cluster identification, for many studies it is important to reliably separate the background seismicity from clustered events for the development of long-term seismic hazard maps [6,7,8] or the regional optimization of background rates [9].

Among the methods that are available for the detection of seismic clusters, one of the most widely used is the window-based approach [10] with known drawback the large gaps after the occurrence of strong earthquakes [11]. The link-based model [12] considers stress redistribution and the Omori law for the determination of the spatiotemporal interactions among earthquakes. The stochastic declustering method of Zhuang et al. [13] is based on the modeling of earthquake occurrences by the ETAS model [14,15], where events are separated into background and clustered ones according to the estimated probabilities. Important clustering features can be inferred using the stochastic algorithm [16,17], whereas it can be also used for declustering to optimize the background seismicity rate estimates [9,18]. Marsan et al. [19] introduced an ETAS model with a time-dependent background component for the detection of aseismic transients. The modified ETAS model is used efficiently to reveal both main shock-aftershocks and earthquake swarms [20,21]. Finally, the Nearest-Neighbor metric proposed by Baiesi and Paczuski [22] adopts a non-parametric definition of a cluster considering the space-time-magnitude proximity among earthquakes. Zaliapin and Ben-Zion [23] introduced a binary threshold, $\eta$, according to which the earthquakes are classified into clusters and background seismicity. They applied the algorithm on a global scale, revealing a link between the clustering properties of seismicity and the heat flow level [24]. The method is proven to be efficient in detecting main shock–aftershock sequences in Northeastern Italy [25], as well as earthquake repeaters in the Sea of Marmaras [26]. Bayliss et al. [27] introduced a new approach to create probabilistic cluster networks based on the intersection between the background and clustered component of the nearest-neighbor distances.

Another approach is based on the assumption of a common physical trigger during a seismic sequence, expressed by fluctuations in the occurrence rate [28,29]. In our method, change-points of the intensity rate in the earthquake occurrences are detected with the use of the MAP model [30]. The temporal distribution of the events is approximated essentially by a non-homogeneous Poisson process, $N_t$, with a piece-wise constant intensity rate determined by the underlying Markov process, $J_t$. Simulation studies and applications on real datasets showed that the model efficiently identifies the changes in the earthquake occurrence rates [31]. Recent works by Lu [32] and Benali et al. [33] are based on non-stationary Poisson models whose rate is modulated by a hidden Markov process to determine a set of change-points for seismicity rate. Concerning earthquake clustering, Bountzis et al. [34] proposed the combination of the MAP model with a density-based clustering algorithm, DBSCAN [35], for the detection of earthquake clusters in the spatiotemporal domain. They used the method on a micro-seismicity earthquake catalog in central Ionian Islands, Greece, efficiently revealing the clustered seismicity.

Several studies suggest that the clustering properties of seismicity (spatiotemporal distribution, productivity rates) might be controlled by the tectonic regime. Llenos and Michael [36] showed that the adoption of region-specific aftershock parameters can improve forecast estimates, as the information from the tectonic region is particularly useful, and suggest the determination of clustering features in smaller regions where high-quality earthquake data are available. More recently, Hardebeck et al. [37] updated the generic parameters of sequences in California incorporating the regionalization of the former work for their determination. In this way, there was an improvement of the aftershock forecasts’ accuracy. The temporal ETAS model assumes that seismicity is evolving in the form of independent events (background seismicity) who generate their aftershocks with each one producing their own. It incorporates two empirical laws, the Omori–Utsu law [38] and the productivity law used to explain the distribution of aftershocks, as well as a stable Poissonian rate for the background seismicity. Our work utilizes the estimated parameters of the ETAS model to investigate regional variabilities in the productivity of the sequences and gain insights into the involved triggering mechanisms [39,40,41].

The aim of this study is firstly to demonstrate the efficiency of the proposed clustering method to separate triggered from background seismicity and subsequently to investigate and compare the clustering properties among three major seismic zones of Greece. In particular, we focus on the statistical analysis of the detected clusters based on the ETAS model producing generic and sequence specific parameters for each area.

The paper is organized as follows. In Section 2, the MAP-DBSCAN method for the identification of the clusters is described and simulation results for the evaluation of the method are deployed. MAP is used as a tool for the detection of changes in the seismicity rate and DBSCAN is implemented to reveal spatially high-density areas. In Section 3, the study area along with the datasets used in the application are presented. Finally, in Section 4, details on the detected clusters for each area are given and the ETAS regional parameters based on the identified clusters are derived. The regional variability of the clustering properties among the three areas is investigated based on the generic estimates of ETAS parameters (a, K, p, c and $\mu$) and sequence-specific parameters. A brief discussion of the results is given in Section 5 and the main conclusions are presented in Section 6.

2. MAP-DBSCAN Method

2.1. MAP as a Tool for the Detection of Seismicity Rate Changes

In the first step, the temporal distribution of seismicity is approximated by a stochastic point model, the Markovian Arrival Process. The MAP is a two-dimensional Markov process ${(N_t,J_t)}_{t\in R^{+}}$, where $N_t$ counts the number of earthquakes up to time t that occur with a rate ${\lambda }_t$. Its value is associated with the unobservable states $i=1,\dots ,K$, of the Markov process $J_t$. In particular, when the process $J_t$ is in state i, earthquakes occur according to a Poisson process with rate ${\lambda }_i$ and, therefore, the sojourn time in this state follows an exponential distribution with expected value $1/{\lambda }_i$. When an earthquake occurs, the MAP can transit with probability $p_{ij}$ to another state j, so now, earthquakes occur according to a Poisson process with rate ${\lambda }_j$, or remain in the current state i with probability $p_{ii}$.

To represent the MAP model, we need the $K\times K$ rate matrices ${\mathbf{D}}_0$ and ${\mathbf{D}}_1$, where ${\mathbf{D}}_0$ is a diagonal matrix whose non-negative elements we denote as ${\lambda }_1,\dots ,{\lambda }_K$, and which correspond to the K Poissonian rates, each one assigned to a hidden state of process $J_t$, and ${\mathbf{D}}_1$ consists of the transition rates among the states, along with the occurrence of an earthquake, which we denote as $q_{ij}$. Additional details on the estimation procedure of the MAP model and its properties are given in Appendix A.

Concerning our clustering algorithm, we are interested in the evaluation of the transitions among the hidden states, namely, the detection of changes in the seismicity rate. In particular, we define a rate threshold, ${\lambda }_{thr}$, according to which a potential sequence starts when the rate of the counting process $N_t$ achieves ${\lambda }_t>{\lambda }_{thr}$ and ends as soon as the process $J_t$ moves for the first time to a state with a Poisson rate below that threshold. Our main assumption is that each state corresponds to a distinct evolution phase of a seismic sequence, independently of its underlying mechanism. In this way, the model has the ability to approximate the temporal evolution of earthquake catalogs that incorporate both aftershock sequences and earthquake swarms [31], as well as datasets with non-stationary characteristics [34].

2.2. Temporal Constraints

The earthquakes above the defined rate threshold comprise the potential clusters. However, results on methods that are based solely on changes in the seismicity rate can sometimes be misleading. One such case is when the rate at the tail of aftershock sequences has reached the level of the background seismicity, so it becomes difficult to discriminate these events from background ones. One similar case is related to the sparse foreshock activity, which, as it is shown in Lippiello et al. [42], exhibits significantly smaller frequency than the aftershock activity. Therefore, a day rule, $dt$, is assigned in the sense that events in $\pm dt$ from the potential cluster are included within. Another case that we observed is related to the existence of fluctuations during a seismic excitation, when the seismic activity that is triggered by the same underlying mechanism is divided into smaller clusters. For this reason, we assign a time window, T, so that clusters in temporal distance smaller than or equal to T are merged into one.

2.3. DBSCAN Algorithm

The merged clusters comprise seismicity concentrated in time. However, events with temporal proximity can be spatially sparse and are falsely assigned into the same cluster. To overcome this ambiguity, a density-based clustering algorithm, DBSCAN, is applied to separate events in space based on a distance metric on the earthquakes’ epicentral distribution. Depending on the adopted distance metric, the algorithm can be used for grouping events with waveform similarities [4] as well as earthquakes with related rupture styles and orientations (focal mechanism clustering) [43]. Density-based algorithms search for areas where the event density exceeds a threshold, $ϵ$. The boundaries of these areas are set where the spatial density falls below that threshold. The DBSCAN algorithm in particular requires as input two parameters, the upper threshold, $ϵ$, and the minimum number of neighboring events, $N_{pts}$. A cluster is defined if an earthquake i exists along with at least $N_{pts}$ events within distance $d\le ϵ$, including itself. Earthquake i is then considered a core point of the cluster and the algorithm moves to the investigation of the other events. If $N_{pts}$ neighbors are identified, they are also considered core events; otherwise, they consist of the boundary points of the cluster and the algorithm stops. Events that have not been assigned to any cluster at the end of the procedure are included to the background seismicity and are merged with events that occurred during periods with estimated rate under the rate threshold, ${\lambda }_{thr}$. In this way, the algorithm can remove events that are sparsely distributed in space. It has been efficiently applied for detecting similarities among earthquake locations, origin times and focal mechanisms [34,44]. An advantage of the algorithm is that it does not require as input a predefined number of earthquake clusters, such as the k-means algorithm, where further optimization techniques for the determination of the clusters number are necessary [45].

2.4. Performance Evaluation

2.4.1. ETAS Framework

The efficiency of the method to correctly identify spatio-temporal correlated seismicity is evaluated on a simulated ETAS catalog, where the underlying structure of the clusters is known a priori. Additionally, we aim to demonstrate the performance of the method compared to widely used clustering algorithms. In particular, our approach is compared with the Nearest-Neighbor (NN), Gardner and Knopoff (GK) window-based and Reasenberg (RB) link-based algorithms. A detailed review on each one of them is given in Appendix B.

The ETAS model belongs to a wide class of branching processes where the occurrence rate of earthquakes, known as intensity function, depends on the history of all previous seismicity and consists of two parts given by

$$\lambda (t,\mathbf{x}=(x,y)/H_t)=\mu \left(\mathbf{x}\right)+K\sum _{j:t_j<t,x_j\in {\Sigma }_0}^{}{e}^{a(m_j-m_c)}g\left(t\right)f\left(\mathbf{x}\right),$$

(1)

where j runs over all past earthquakes with magnitude larger than or equal to $m_c$. The intensity function is evaluated at each event that occurred during the time interval $[t_{in},T]$ and inside the target region $\Sigma \subseteq {\Sigma }_0$. However, events in the broader region ${\Sigma }_0$ and time interval $[t_0,T]$, with $t_0<t_{in}$, can be considered triggering events, therefore, they should be included in the evaluation of $\lambda (t,\mathbf{x})$. The first term of the right-hand side of Equation (1) expresses the background seismicity $\mu \left(\mathbf{x}\right)$ which is assumed stationary in time (mother events) and clustered in space due to the fault network geometry. The latter term represents the space-time-dependent seismicity (daughter events) expressed through the following empirical laws:

• the productivity law, $K{e}^{a(m_j-m_c)}$, which gives the number of aftershocks triggered by a main shock with magnitude $m_j$;
• the modified Omori law, $g\left(t\right)=(p-1)c^{(p-1)}{(t-t_j+c)}^{-p}$, with $p>1$, which describes the temporal decay of aftershocks;
• the spatial distribution of aftershocks, $f(\mathbf{x}/M)=\frac{q-1}{\pi d{\left(m_j\right)}^{q-1}}[\parallel \mathbf{x}-{\mathbf{x}}_j{{\parallel }_2^2+d\left(m_j\right)]}^{-q}$ with $q>1$, and $d\left(m_j\right)=d_0{10}^{\gamma (m_j-m_c)}$, which assumes an isotropic distribution of aftershocks around the main shock.

Each mother event generates its daughters (first generation), the daughters generate their own descendants (second generation), and so on. In this way, a cluster is defined as a sequence of events with a common mother event (first event in the cluster). There are also single events, i.e., mother events without any subsequent triggered earthquake.

2.4.2. Simulation Procedure

For the simulation of an ETAS earthquake catalog, first, we need to generate the mother events. The heterogeneity in space can be preserved from the spatial coordinates of a declustered earthquake catalog. In particular, we implement a declustering procedure to the earthquake catalog of Greece (${\Sigma }_0=[{19}^{\circ }E-{29}^{\circ }E]\times [33.7^{\circ }N-{42}^{\circ }N]$) for earthquakes with $m_c=2.5$ during the period of 2011–2019, using the NN method, and then we produce $N_{main}$ mother events according to a Poisson distribution with mean value equal to the number of the identified background events. Their coordinates are sampled with replacement from the declustered catalog by adding a random factor. The occurrence times are simulated from a uniform distribution $U(t_0,T)$, where $t_0=0$ and $T=20$ years.

The magnitudes are independent from the earthquakes’ spatial and temporal distribution and follow the Gutenberg–Richter (GR) law truncated from the left at the completeness magnitude, $m_c=2.5$, and from the right at the maximum observed magnitude of the instrumental earthquake records in Greece plus a small factor, $m_{max}=7.8$. The functional form of their distribution is the following, $s\left(m\right)=\left(\beta e^{-\beta m}\right)/(e^{-\beta m_c}-e^{-\beta m_{max}})$, where $\beta$ relates to the b-value of the GR law with $\beta =blog\left(10\right)$, and we chose $b=1.0$. After the generation of the background events, we simulate their aftershock number, following Poisson distribution with expected rate equal to the productivity of the model, $k\left(m_j\right)=K{e}^{a(m_j-m_c)}$. Their occurrence times are sampled from the modified Omori law, $g\left(t\right)$, and the locations from the isotropic spatial distribution function, $f\left(\mathbf{x}\right)$. For next-generation daughters, the triggering step is repeated until there are no more generated events. In

Table 1. ETAS parameters used for the synthetic earthquake catalog with $m_c=2.5$. The target area is Corinth Gulf, Greece, with $\Sigma x[t_{in},T]=[21.3^{\circ }E-23.2^{\circ }E]\times [37.9^{\circ }N-38.6^{\circ }N]\times [2,20]$. The number of clustered events is $N=4253$ and the number of mother events is $N_{bg}=1595$.

Parameter		Parameter
K	0.1	d	2.41 $\times {10}^{-5}$
a	2.19	q	1.805
p	1.13	$\gamma$	0.59
c	0.024 (days)	$\mu$ (events/day)	4.50

, we give the parameter set that produced the ETAS catalog.

2.4.3. Evaluation

The ETAS earthquake catalog is divided into either single events or mother events with their descendants. We define by $\mathbf{X}={\left\{{\mathbf{X}}_k\right\}}_{k=1,\dots ,N_c}$ the true partition of the catalog, where $N_c$ corresponds to the number of clusters, and with ${\mathbf{Y}}_i={\left\{{\mathbf{Y}}_n\right\}}_{n=1,\dots ,N_{y_i}},i=1,\dots ,K,$ the partition after the implementation of the MAP-DBSCAN and the other $K-1$ methods, including different tuning of the algorithm’s parameters. $N_{y_i}$ is the number of clusters after the implementation of method i.

Next, we define the Jaccard index [46], which is a measure to quantify the overlap between two partitions, in our case, the true one of the ETAS catalog, $\mathbf{X}$, and the one of the $i$-th implemented algorithm, ${\mathbf{Y}}_i$. The Jaccard index is expressed by $J_1(\mathbf{X},\mathbf{Y})=a_{11}/(a_{11}+a_{10}+a_{01})$, where $a_{11}$ indicates the number of pairs of elements which are correctly assigned into the same cluster (true links), $a_{01}$—the number of pairs of elements which are in the same cluster in the ETAS catalog and in different clusters in the estimated one (missed links) and $a_{10}$—the number of pairs of elements which are wrongly identified as clustered events (false links). If all the initial clusters are correctly identified by the implemented method, then $a_{10}=0=a_{01}$ and $J_1(\mathbf{X},\mathbf{Y})=1$. Conversely, if all pairs are wrongly identified as clustered or independent, then $a_{11}=0$ and, as a consequence, $J_1(\mathbf{X},\mathbf{Y})=0$.

In addition, we introduce a generalization of the Jaccard index, $J_2(\mathbf{X},\mathbf{Y})=b_{11}/(b_{11}+b_{10}+b_{01})$, to identify the partition $\mathbf{Y}$ with the best discrimination between the background seismicity and clustered elements, following the definition in Lippiello and Bountzis [47]. We consider as background seismicity single events and the mother events of each cluster, i.e., the one that initiated a cascade of events. Here, $b_{11}$ represents the number of common background events in the two partitions, $b_{10}$ is the number of elements wrongly identified as mother events in the partition $\mathbf{Y}$, whereas $b_{01}$ corresponds to the number of true mother events identified as clustered elements in the partition $\mathbf{Y}$. In

Table 2. $J_i,i=1,2,$ values for 3 parameter sets (PS) of RB and GK algorithms, respectively, and the corresponding values of the MAP-DBSCAN and NN methods.

PS	RB1	RB2	RB3	GK1	GK2	GK3	MAP-DBSCAN (PS27)	NN
$J_1$	0.530	0.593	0.648	0.382	0.397	0.585	0.627	0.756
$J_2$	0.612	0.630	0.617	0.418	0.192	0.676	0.647	0.727

, we show the Jaccard index values $(J_i,i=1,2)$ after the implementation of the different clustering algorithms. In particular, for the MAP-DBSCAN algorithm we show the one with the best results in terms of the Jaccard index (MAP-DBSCAN27). In Appendix B, the results for all input parameters are given.

The window-based method removes all the events within $d\left(M\right)$ kilometers and $t\left(M\right)$ days after a main shock with magnitude M. We used three different temporal and spatial intervals, given by Equations (A2) (GK1), (A3) (GK2) and (A4) (GK3). The Reasenberg algorithm combines a deterministic spatial window and a probabilistic temporal one, determined by the Omori law. We used three sets of parameters in the ZMAP tool. RB1 (

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

PS	${\mathit{\tau }}_{\mathit{min}}$	${\mathit{\tau }}_{\mathit{max}}$	${\mathit{p}}_1$	${\mathit{x}}_{\mathit{k}}$	${\mathit{x}}_{\mathit{meff}}$	${\mathit{r}}_{\mathit{fact}}$
RB1	1	10	0.95	0.5	2.5	10
RB2	1	10	0.95	0.5	2.5	20
RB3	0.5	20	0.95	0.5	2.5	20

) corresponds to the original parameters proposed in Reasenberg [12]. In the second set, RB2, we extend the spatial zone by increasing the factor $r_{fact}$ from 10 to 20 km, whereas, in the third set, RB3, we also extend the temporal window modifying the parameters ${\tau }_{min}$ and ${\tau }_{max}$ (see Table A1). The NN method separates seismicity into background and clustered events according to the bimodal distribution of the rescaled time and distance metrics. The algorithm needs as input two parameters, the b-value and the fractal dimension of the earthquake locations, which are considered equal to $b=1.0$ and $d_f=1.51$, respectively. For the MAP-DBSCAN method, first, the optimal MAP model and the corresponding rate threshold (${\lambda }_{thr}={\lambda }_1$) are determined. Then, different temporal constraints are tested for the merging of the potential clusters (consecutive events above the rate threshold) and, subsequently, the DBSCAN is implemented. We set the minimum number of neighbors equal to $N_{pts}=2$, similar to the minimum size of a cluster that can be given as output from the other algorithms. For the distance threshold, $ϵ$, we tested a wide range of possible values (see

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

$\mathit{ϵ}$	${\mathit{N}}_{\mathit{pts}}$	PS	T	$\mathit{dt}$	PS	T	$\mathit{dt}$
[2.5 5 7.5 10 12.5]	2	1–5	0	0	16–20	0	7
		6–10	7	0	21–25	7	7
		11–15	14	0	26–30	14	7

).

The NN method shows the best performance in the construction of the clusters ($J_1=0.756$) and in the detection of the mother events ($J_2=0.727$). This is also evident by its cumulative number of background seismicity (purple line in Figure 1a), which is the closest one to the initial catalog (dotted black line). The temporal evolution of background seismicity is shown in Figure 1b–f across the longitude for ease of reading as west–east normal faults dominate the area. For the NN method, no large gaps are evident in the space-time evolution of the declustered seismicity, although there is a significant concentration of events between the 7th and 8th year of the catalog, which is also persistent in both RB2 and GK3 methods (orange ellipses in Figure 1d–f) and less apparent on MAP-DBSCAN method (orange ellipse in Figure 1c). The high efficiency of the NN method is probably related to the metric it uses, which is similar to the ETAS one with ${\lambda }_j(t_i,x_i)={(t_i-t_j)}^{-1}{r}_{ij}^{-d_f}{10}^{b{m}_j}$ and $c=0,p=1,d=0,q=d_f$ and $a=b$. The windowing technique seems to overestimate the temporal and spatial windows, since it removes large amounts of seismicity (blue ellipse in Figure 1f), in accordance with previous results [11]. The same gap between the 14th and 15th year of the catalog is also evident in the background seismicity from the MAP-DBSCAN method, however, it is smaller and some sparse seismicity is left (blue ellipse in Figure 1c). On the other hand, Reasenberg’s declustered catalog has more events than any other method (pink, magenta and green line, Figure 1a) and significant concentrations of events are visible in the space-time evolution of the background seismicity (orange and purple ellipses in Figure 1e).

Best overlapping among the true, $\mathbf{X}$, and the estimated partition, $\mathbf{Y}$, does not mean necessarily the best detection for the declustered seismicity. For instance, the GK3 partition is characterized by a lower index, $J_1=0.585$, than the MAP-DBSCAN27 partition, $J_1=0.627$, however, its declustering catalog is more accurate (GK3-$J_2=0.676>J_2=0.647$-MAP-DBSCAN27). Nevertheless, both indexes combined, the MAP-DBSCAN partition shows a higher efficiency than the rest of the algorithms, except for the NN. The Jaccard index values for the rest of the MAP-DBSCAN input parameters are quite stable with small fluctuations from the best parameter set (MAP-DBSCAN27), apart from the smallest distance cutoff, $ϵ=2.5$ km, which seems inadequate for capturing the spatial correlations among the events (Appendix B.4).

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 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 S1a–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,0.89$, for the D1, D2, D3 datasets, respectively (

Table 3. The magnitude of completeness, $M_c$, for the datasets of the three areas CG, CII and NAS, along with the productivity, a, and the b-value of the GR law. N and $N_c$ denote the initial number of events and the ones with $M\ge M_c$, respectively.

Region	Notation	N	${\mathit{M}}_{\mathit{c}}$	${\mathit{N}}_{\mathit{c}}$	a	b
CG	D1	25,595	1.5	13,043	5.57	0.97
CII	D2	24,085	2.2	6981	5.80	0.88
NAS	D3	21,139	2.1	8328	5.79	0.89

). The resulting magnitude threshold for the three datasets is equal to $M_c=1.5,2.2,2.1$, with 13,043, 6981, 8328 events (Table 3), respectively, the epicenters of which are shown in the maps of each study area (Figure 2).

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 $p_i\left(t\right)=p(J_t=i)$ for $i=1,\dots ,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 $J_t$, the one with the highest probability, i.e., ${argmax}_{0\le i\le K}{p}_i\left(t\right)$, with $p_i\left(t\right)=p_i\left(t_k\right)$ for $t_k\le t<t_{k+1}$. Each state i corresponds to an occurrence rate, ${\lambda }_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 $t_k\le t<t_{k+1}$ indicates the state with the maximum probability at the current time and the legend contains its corresponding occurrence rate.

The temporal patterns of dataset D1 indicate the dominance of state 2 (yellow color, Figure 3a) with occurrence rate ${\lambda }_2=3.01$ events/day for almost the entire period. Nevertheless, there is a slight decrease in the occurrence of earthquakes (${\lambda }_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 ${\lambda }_{thr}={\lambda }_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 ($M_{w}6.1$ and $M_{w}6.0$) [58] and the $M_{w}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 ${\lambda }_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 $M_{w}5.8$ [60], the 2014 $M_{w}6.9$ Samothraki [61] and the 2017 $M_{w}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 ${\lambda }_{thr}={\lambda }_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 ${\lambda }_{thr}$ are classified into groups which we call potential clusters, and then, we test four different sets of temporal constraints, $(T,dt)$, to the three datasets. Potential clusters within a temporal interval T are merged into one and events that occurred in $\pm 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 ($N_{pts}=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.

4.2. Cluster Analysis

Table 4. Cluster statistics and the parameter set of the clustering algorithm for the three datasets. $N_{clust}$ corresponds to the number of clustered events and $N_{bg}$ to the background seismicity frequency. $\overline{\tau }$ and $\overline{n}$ are the mean duration in days and size of the clusters, respectively.

Dataset	$(\mathit{T},\mathit{dt},\mathit{ϵ},{\mathit{N}}_{\mathit{pts}})$	${\mathit{N}}_{\mathit{clust}}$	${\mathit{N}}_{\mathit{bg}}$	# Clusters	$\overline{\mathit{\tau }}$	$\overline{\mathit{n}}$
D1	$(0,5,2.5,4)$	4662 (36%)	8381 (64%)	255	12.50	18.28
D2	$(5,5,2.5,4)$	5221 (75%)	1770 (25%)	45	54.60	118.43
D3	$(5,5,5,4)$	3688 (44%)	4640 (56%)	187	15.08	19.72

gives the chosen parameter set of the clustering algorithm for each dataset based on the analysis in Appendix C and a summary on the statistics of the detected clusters. In the CG area, we identified the largest number of seismic clusters (255) due to the increased detectability of micro-seismicity (low completeness magnitude threshold), however, they are short in size ($\overline{n}=18.28$) and duration ($\overline{\tau }=12.50$) Conversely, the CII area is characterized by a small number of seismic clusters (45) but with large mean size ($\overline{n}=118.43$) and duration ($\overline{\tau }=54.60$). The clustered seismicity is prevalent ($75\%$), whereas in CG and NAS, the background component is more dominant than clustered seismicity with $64\%$ and $56\%$, respectively (Table 4). In CG, this is explained by the lack of large main shocks during the study period and the occurrence of few moderate events, the largest number with $M=$5.2.

Concerning the CG area, the majority of the clusters are located on the western subarea where 22 out of 27 clusters with $N\ge 30$ occurred. The main activity is located offshore between Aigion and Trizonia Island, but also north of the Psathopyrgos fault (Figure 4a). The eastern subarea comprises smaller clusters that are mainly concentrated offshore Xylokastro and Perachora faults, as well as near Itea Gulf (Figure 5a). The seismicity of CII is dominated by the two major main shock–aftershock sequences, each sequence comprising 2829 and 1396 events, respectively. Essentially, 4225 out of the 5221 clustered events belong to these sequences (Table 4). Furthermore, 45 clusters are detected in total with the main activity concentrated along the KTFZ (Figure 6a). The NAS area comprises 187 clusters, including both main shock–aftershock sequences and earthquake swarms (Table 4). Figure 7a shows that the main clustered activity is concentrated along the NAT and the sub-parallel branches, as well as in the southeastern subarea.

4.2.1. Corinth Gulf Area

The western subarea of the Corinth Gulf is characterized by rich seismic activity, especially in 2013–2014, when 13 out of the 22 clusters with $N\ge 30$ occurred. One of the major detected sequences is the 2013 Aigion swarm ($C6$ in Figure 4b) which initiated on 21 May 2013 with a bulk of small events and several bursts associated to earthquakes with magnitudes ranging between 3.3–3.7 (Figure S3) [63,64]. Two distinct excitations followed ($C8$ and $C10$ in Figure 4b) in accordance with the ones observed by Michas et al. [65]. The first cluster began on 7 July with some activity prior to the $M=3.7$ event on 15 July 2013, and lasted until 27 August, 2013 (Figure S3). The second half of 2014 is also a well-studied period with intense seismic activity. Five clusters with $N\ge 30$ are detected ($C15,C16,C18,C19$ and $C20$) in the western subarea, including the offshore $M4.8$ earthquake on 7 November 2014, associated with $C19$ (Figure 4b), and the $M4.6$ event on 21 September 2014, associated with the earthquake swarm located between Nafpaktos and Psathopyrgos [66] ($C15$ in Figure 4b). Persistent activity since 22 July 2014 is also observed offshore Aigion ($C16$), close to the earthquake swarm, $C15$, which began on 7 November 2014 (Figure S5). In 2012, fewer clusters are observed, mostly during the first semester, with three clusters comprising $N\ge 30$ events, $C1,C2$ and $C3$, and a plethora of smaller clusters (Figure 4b and Figure S2). Between November 2013 and July 2014, the activity is sparse with three relatively large clusters, $C11,C12$ and $C14$ (Figure 4b and Figure S4). Six more clusters with $N\ge 30$ are observed until the end of 2017 ($C21,C22,C23,C24,C25,C27$, Figure 4b).

The eastern subarea is characterized by more sparse activity. A major seismic sequence, Offshore Perachora ($C4$ in Figure 5b), is detected, including two sub-sequences, the first initiated on 22 September and the second on 30 September 2012 (Figure S6). Two relatively large clusters, $C13$ and $C17$, are observed near Itea; the former lasted almost two weeks at the end of March, 2014, and the latter—almost three months between August and October 2014 (Figure 5 and Figure S7).

4.2.2. Central Ionian Islands Area

The two main shocks of sequence $I1$ (Figure 6b) with $M=6.1$ and $M=6.0$ occupy the southern and the central part of the onshore area of Kefalonia Island. The 2014 Kefalonia earthquake sequence ($I1$ Figure 6) started on 19 January with the first main shock occurring on 26 January ($M=6.1$), and aftershock activity extending over 35 km [58], part of which hosted the second main shock ($M=6.0$) that occurred on 3 February and the compound aftershock activity. A sub-cluster is also detected offshore to the southwest of Kefalonia Island ($I2$ in Figure 6b) that is deployed concurrently with the main sequence (Figure S8). In addition, two distinct clusters, $I3$ and $I4$ (Figure 6b), are revealed, which occurred between November and December 2014 (Figure S9), across the edges of the double rupture. They might be triggered by the stress transfer of the main ruptures, indicating activation of adjacent fault segments. The seismic activity of cluster $I5$ (Figure 6b) retains the most interest because it is essentially two seismic excitations evolving at the same time. The first initiated in the Myrtos Gulf and the second offshore the south part of Kefalonia Island. It comprises 164 earthquakes in about 100 days (Figure S9). The activity of the $I7$ cluster (Figure 6b and Figure S10) spreads along the western coastline of Lefkada and Kefalonia Islands, far beyond both sides of the 2015 Lefkada main rupture. To the south the aftershock activity is sparse, probably due to the large amount of stress released in the main rupture, revealing that the main slip is associated with a fault of about 17 km in length [59]. Apart from cluster $I4$, two additional clusters ($I6$ and $I9$ in Figure S9 and Figure S11, respectively) are detected in the area between Lefkada and Kefalonia, extending to about 15 km, which is considered as a transition zone encompassing step-over structures [58]. All of them relate to the E–W-oriented, parallel step-over faults, similar to the ones detected in the microseismicity cluster analysis between September 2016 and December 2019 in the study area [34].

4.2.3. North Aegean Sea Area

The first seismic excitation with $N\ge 30$ events ($N1$ in Figure 7b) is a sequence of interest since two moderate events ($M=5.2$ and $M=5.3$) occurred in 3 weeks, both producing their own aftershocks (Figure S12). The 2013, January 8 $M=5.8$ North Aegean earthquake [60] along with its aftershock activity (cluster $N3$ in Figure 7b) is also detected. The aftershock activity is temporally divided into two clusters (Figure S13). The 24 May 2014 $M=6.9$ Samothraki main shock was followed by aftershock activity confined to three major clusters ($N4,N5,N6$ in Figure 7b) and some secondary clusters with $N\ge 10$ events (Figure S14), which are in accordance with the ones observed by Saltogianni et al. [61]. The seismic activity that took place near the Aegean coast of NW Turkey during January–October 2017 [67] is divided into three clusters with $N\ge 30$ ($N10,N11$ and $N14$ in Figure 7b) and two minor clusters with 22 and 23 events, respectively (Figure S15). The strong main shock ($M=6.4$) that occurred on the 12th of June 2017 offshore, south of the SE coast of Lesvos Island, along with its intense aftershock activity, is revealed and illustrated in Figure 7b ($N12$). Two major ($N\ge 30$) secondary outbursts of clustered activity occurred concurrently on the west ($N17$) and east ($N16$) side of the sequence (Figure S16). A thorough analysis revealing multiple spatial clusters of the sequences is conducted by Papadimitriou et al. [62].

4.3. Regional Variability of Clustering Properties

In this section, we investigate regional variations in the clustering behavior of the detected seismic sequences, in particular, on their productivity rates and on their temporal evolution that can differ among areas with distinct seismotectonic characteristics. Therefore, we adopt the temporal ETAS model that expresses two empirical relationships that characterize the temporal and size distribution of earthquakes, the normalized Omori–Utsu law, given by $g\left(t\right)=c^{p-1}(p-1){(t-t_i+c)}^{-p}$, and the productivity law that is expressed by $N=k\left(M_i\right)=K{e}^{a(M_i-m_c)}$, where N is the number of triggered events by an earthquake of magnitude, $M_i$, K is a constant of proportionality, which depends on the number of triggered events per mainshock above the catalog cutoff, and a describes the impact of magnitude on the number of triggered events.

We compute the generic parameters for the 3 areas, CG, CII and NAS, by jointly inverting the ETAS parameter set $\theta =(p,c,a,K,\mu )$ from the identified sequences with $N\ge 30$ of each area. In particular, $L{L}_i={\sum }_{j=1}^{n_i}log\lambda \left(t_j\right)-{\int }_{t_0}^{t_{end}}\lambda \left(t\right)dt$ denotes the log-likelihood of the i-th sequence for each area, namely, the logarithmic probability of observing $n_i$ events with occurrence times $t_j,j=1,\dots ,n_i$, during the period of the sequence $(t_0,t_{end})$, and no other events between them. The intensity function, ${\lambda }_t$, of the model is given by Equation (1), neglecting the spatial component. We then stack all the sequences of each area, compute their corresponding logarithmic probabilities $L{L}_i$, and define as the common log-likelihood

$$LL=\sum _{i=1}^{N^{*}}L{L}_i,$$

(2)

where $N^{*}$ is the number of sequences. The optimal inverted parameters are the ones that maximize Equation (2). The results of the ETAS parameter estimation for the three regions are shown in

Table 5. Generic ETAS parameter values for the three study areas. $N^{*}$ denotes the number of sequences with $N\ge 30$.

Area	p	c	a	K	$\mathit{\mu }$	$\mathit{\beta }$	${\mathit{N}}^{*}$
CG	1.23	0.0171	0.82	0.74	0.43	2.13	26
CII	1.31	0.11	1.29	0.44	0.15	2.21	9
NAS	1.26	0.0324	1.04	0.51	0.28	2.03	17

. There are 27 sequences in CG (Figure 4 and Figure 5), 9 in CII (Figure 6) and 17 in NAS (Figure 7) from 2012 until 2019 with $N\ge 30$ events, however, we removed cluster C26 from the computations, since it is located at the boundaries of the study area (Figure 5) with part of the aftershock data being omitted. For the maximization of the common log-likelihood $LL$, we implement an iterative procedure where at each step we update the model parameters by a random factor so ${\theta }_k^{new}={\theta }_k+u$, for $k=1,\dots ,5$, then, we compute the corresponding $L{L}_i^{new},i=1,\dots ,N^{*}$, log-likelihood values and store the new parameters under the condition $L{L}_{new}>LL$ moving to the next iteration. After some iterations, the logarithm converges and the algorithm stops. Essentially, this is a grid-based procedure, since we use a large number of iterations.

The parameter a for CG ($a=0.82$) is the lowest among the three areas, indicating the dominance of swarm activity presumably due to fluid flow in accordance with many relevant studies [65,68]. Low a values characterize areas with high heat flow [39], even though the estimated value can be underestimated due to magnitude incompleteness after the occurrence of the main shock or due to the existence of time-dependent background seismicity [69]. Conversely, in CII, the estimated value ($a=1.29$) is relatively larger compared to the former region ($a=0.82$), indicating the dominance of typical main shock–aftershock sequences. In the NAS area, a moderate value is acquired ($a=1.04$), probably due to the co-existence of swarm activity and aftershock sequences. Another indicator for the existence of swarm activity in CG is the large value of the background seismicity ($\mu =0.43$) compared to NAS and CII. High values of the background rate can indicate the existence of aseismic loading transients [40]. LLenos et al. [70] observed increased values of the background component of the fitted ETAS model when it was applied to pre-swarm and swarm activity, respectively.

For the comparison of the productivity among the three areas, since they have different completeness magnitudes, we use the following relation,

$$N=k\left(M_i\right)P(M\ge m_c^{*})=K{e}^{a(M_i-m_c)}{e}^{-\beta (m_c^{*}-m_c)},$$

(3)

which yields the number of earthquakes above magnitude $m_c^{*}$, generated by a main shock of magnitude $M_i$. Figure 8 shows the number of direct triggered events, $K\left(M_i\right)$, from an earthquake of magnitude $M_i$. We consider $m_c^{*}=2.2$, which is the maximum completeness magnitude among the three datasets. The exponent of the exponential magnitude distribution is expressed by $\beta$ and is defined as $\beta ={\sum }_{i=1}^N{\beta }_i/N^{*}$, where $N^{*}$ is the number of clusters for each dataset and ${\beta }_i$ their corresponding exponent values. Concerning the distribution of aftershocks in time, the normalized Omori law distribution is used, given by $g\left(t\right)$.

In CII, the seismic sequences seem to be more productive, as shown in Figure 8, with NAS and CG to exhibit smaller values. Combined with the higher background rate for the area of CG ($\mu =0.43$), we may say that a significant part of Corinth Gulf’s sequences cannot be contributed to the triggering effect of mainshocks but different underlying mechanisms seem to play an important role. Conversely, in CII area, mainshock–aftershock sequences seem to dominate, generating a rich number of aftershocks (very low background rate, $\mu =0.15$, and high productivity of mother events).

4.4. Sequence-Specific Clustering Properties

Next, we estimate the a-values for the individual sequences of each area by maximizing $LL$ as a function of a, K and the background rate $\mu$, while keeping the rest of the parameters fixed for clusters with $N<80$. In this way, we increase the robustness of the inversion procedure since there are sequences with few events. A similar procedure was followed by Page et al. [3] and Llenos and Michael [36] who demonstrated that fitting multiple parameters for a single sequence can be unstable and Hardebeck et al. [37] who implemented this method for the estimation of California aftershock parameters. We intend to investigate potential differences in the productivity ($a,K$) and the background rate, $\mu$, among sequences of each area and their relation to different underlying triggering mechanisms. Productivity parameters a and K are correlated, so we enabled both to run during the iterative procedure. We also examine the value of the background rate among sequences since it can be also an indicator of aseismic transients in a region. Both parameters, a and K, are not influenced by $\mu$, as we verified it by implementing the inversion procedure, also keeping parameters a and K fixed.

4.4.1. Corinth Gulf

In

Table 6. Details on the 26 clusters with $N\ge 30$ events in CG area and the inverted ETAS parameters. The generic values of the Omori law, p and c, are adopted for clusters with $N<80$.

ID	${\mathit{T}}_{\mathit{in}}$	${\mathit{T}}_{\mathit{end}}$	N	p	c	b	a	K	$\mathit{\mu }$	${\mathit{M}}_{\mathit{max}}$
C1	12/1/12	23/1/12	33	1.23	0.017	1.20	0.49	0.79	1.03	3.1
C2	13/1/12	27/1/12	33	1.23	0.017	0.83	1.69	0.23	1.25	3.1
C3	4/3/12	6/4/12	65	1.23	0.017	1.03	1.53	0.26	1.05	3.0
C4	22/9/12	3/10/12	69	1.23	0.017	0.99	0.36	0.94	1.44	5.0
C5	27/12/12	1/1/13	34	1.23	0.017	0.82	1.84	0.21	1.32	3.8
C6	22/5/13	28/6/13	310	1.45	0.012	0.96	0.20	0.90	0.47	3.7
C7	8/6/13	28/6/13	144	1.11	0.007	1.22	0.60	1.30	1.00	3.0
C8	7/7/13	27/7/13	128	1.04	0.001	0.77	0.34	2.48	0.59	3.7
C9	8/9/13	13/9/13	65	1.23	0.017	1.19	1.28	0.79	2.74	2.8
C10	29/10/13	6/11/13	68	1.23	0.017	1.27	0.10	0.91	2.87	3.1
C11	19/1/14	16/1/14	33	1.23	0.017	0.92	1.26	0.50	1.37	3.8
C12	29/1/14	10/2/14	70	1.23	0.017	0.81	1.39	0.29	1.92	3.9
C13	21/3/14	1/4/14	52	1.23	0.017	0.83	2.97	0.009	3.41	4.0
C14	8/6/14	11/6/14	74	1.23	0.017	0.81	0.92	0.64	4.86	4.3
C15	21/7/14	31/10/14	506	1.37	0.051	1.04	1.38	0.34	1.32	4.6
C16	22/7/14	1/11/14	95	1.26	0.014	1.15	0.72	0.45	0.44	2.8
C17	24/7/14	26/10/14	61	1.23	0.017	0.94	1.72	0.16	0.35	3.4
C18	23/7/14	31/10/14	121	1.25	0.131	0.95	1.77	0.24	0.05	4.7
C19	7/11/14	18/12/14	228	1.07	0.071	0.92	1.80	0.55	0.76	4.8
C20	7/11/14	14/12/14	36	1.23	0.017	1.05	1.27	0.41	0.42	3.1
C21	1/10/15	6/10/15	44	1.23	0.017	1.16	1.97	0.49	1.61	2.8
C22	27/7/16	5/8/16	32	1.23	0.017	0.75	3.50	0.09	0.45	2.7
C23	1/8/16	8/8/16	147	2.79	0.160	0.98	0.10	0.85	2.98	3.4
C24	9/1/17	23/1/17	104	2.79	0.702	0.82	1.70	0.15	1.05	4.5
C25	14/7/17	17/7/17	39	1.23	0.017	0.43	0.73	0.40	5.95	4.2
C27	30/10/17	2/11/17	31	1.23	0.017	0.50	1.68	0.10	6.19	3.5

, the inverted parameters for the 26 clusters of dataset D1 with $N\ge 30$ are given. We adopt the generic values of Omori law (p and c) for clusters with $N<80$ to increase the stability of the inverted parameters. We observe relatively high background rates for most of the sequences and low a values, in particular, $a<1$ for 10 out of the 26 clusters.

Concerning the 2013 Aigion earthquake swarm and its subsequent swarms (clusters $C6,C9$ and $C10$), we observe relatively low productivity values of the ETAS model ($a=1.31,1.29,1.13$, Table 6) in accordance with studies suggesting pore-fluid pressure as the main triggering mechanism during the excitation [63]. Clusters $C11$ and $C12$ are part of the same swarm (Figure S4) that occurred offshore Psathopyrgos. Their relatively high background rates ($\mu =1.34,1.92$) show that a significant part of the clustered seismicity cannot be explained by the empirical laws of the triggering part of the ETAS model. Cluster 14 is part of a major swarm that began on 8 June 2014 (Figure S4). Michas et al. [65] did not find high diffusion rates that are related to fluid pore pressure. However, the large background rate found in our study ($\mu =4.86$) and the low a value ($a=0.92$) suggest the existence of a non-typical mainshock–aftershock sequence, with more complex triggering mechanisms being responsible, such as aseismic creep. Similarly, the largest cluster in the dataset, the $C15$, located offshore Nafpaktos, is characterized by relatively high background rate ($\mu =1.32$) and low productivity ($a=1.38$), more typical values for swarm activity. In contrast, clusters $C18$ and $C19$ that are more typical mainshock–aftershock sequences with a distinct in magnitude event in the initiation of the sequence (Figure S5), have low background rates ($\mu =0.05,0.76$) and relatively high productivity rates ($a=1.77,1.80$). The two clusters near Itea Gulf show contradictory results, in particular, the first one, $C13$, is characterized by a high background rate ($\mu =3.41$), whereas the second, $C17$, which occurred four months later, exhibits a much smaller background value ($\mu =0.35$) more typical for mainshock–aftershock sequences. However, biases can exist in the inversion of the parameters for clusters with a small number of events, so we should be cautious with the inference.

4.4.2. Central Ionian Islands

In

Table 7. Details on the 9 clusters with $N\ge 30$ events in CII area and the inverted ETAS parameters. The generic values of the Omori law, p and c, are adopted for clusters with $N<80$.

ID	${\mathit{T}}_{\mathit{in}}$	${\mathit{T}}_{\mathit{end}}$	N	p	c	b	a	K	$\mathit{\mu }$	${\mathit{M}}_{\mathit{max}}$
I1	19/1/14	16/9/14	2829	1.42	0.24	0.79	1.31	0.40	0.17	6.1
I2	23/1/14	14/9/14	55	1.31	0.11	1.23	1.38	0.30	0.12	3.7
I3	5/11/14	11/12/14	134	1.36	0.06	0.99	1.44	0.29	0.99	5.1
I4	13/11/14	12/12/14	66	1.31	0.11	0.93	1.43	0.38	0.37	4.9
I5	5/1/15	27/4/15	164	1.05	0.01	0.93	2.82	0.10	0.76	4.4
I6	18/1/15	24/4/15	71	1.31	0.11	1.08	1.91	0.36	0.15	3.8
I7	13/11/15	26/6/16	1396	1.45	0.30	0.86	1.51	0.29	0.45	6.5
I8	20/11/15	25/6/16	65	1.31	0.11	0.84	0.94	0.53	0.07	4.3
I9	4/4/17	4/5/17	67	1.31	0.11	0.95	2.26	0.18	0.70	3.9

, the inverted parameters for the nine clusters identified in the area of CII with $N\ge 30$ events are given. We maintain fixed the Omori law parameters p and c (generic values) for clusters with $N<80$. The estimated ETAS parameters of the sequence $I1$ are in accordance with the existence of a main shock–aftershock sequence described in Section 4.2. In particular, the background rate is relatively low ($\mu =0.17$), indicating that the seismicity is adequately described by the triggering part of the ETAS intensity function. The seismic activity of clusters $I3$ and $I5$ (shown by green and blue color in Figure S9, respectively) are characterized by relatively high background rates ($\mu =0.99,0.76$, Table 7). The space-time evolution of the former indicates a rapid migration in the beginning of the sequence (Figure S9), whereas, for the latter, it is characterized by the smallest K value ($K=0.04$) in the area although the a value is rather large. Taking into account the lack of distinct main shocks at the initiations of the sequences, they can be characterized as earthquake swarms, one of the few observed in an area which comprises mostly main shock–aftershock sequences. Concerning cluster $I9$, located in the transition zone between Lefkada and Kefalonia Islands, there is evidence for swarm activity due to the relatively high background seismicity rate ($\mu =0.70$). Ultimately, the major main shock–aftershock sequences in the area, $I1,I7$, have the highest p values ($p=1.42,1.45$), meaning that they are characterized by high aftershock decay in time.

4.4.3. North Aegean Sea

In NAS area cluster $N1$, which consists of two moderate events ($M=5.2$ and $M=5.3$) within a time period of 3 weeks, exhibits the lowest a value ($a=1.10$) among the main detected clusters, which could be an indicator of fluid diffusion in the area (

Table 8. Details on the 17 clusters with $N\ge 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$.

ID	${\mathit{T}}_{\mathit{in}}$	${\mathit{T}}_{\mathit{end}}$	N	p	c	b	a	K	$\mathit{\mu }$	${\mathit{M}}_{\mathit{max}}$
N1	14/2/12	4/4/12	136	1.41	0.03	0.96	1.10	0.40	0.54	5.3
N2	27/4/12	3/5/12	30	1.26	0.03	0.53	1.64	0.16	1.07	4.8
N3	8/1/13	6/3/13	285	1.07	0.06	0.88	2.39	0.06	0.55	5.8
N4	24/5/14	9/7/14	94	1.41	0.78	0.74	1.82	0.01	0.16	6.9
N5	24/5/14	11/7/14	153	1.60	0.16	0.69	1.76	0.16	0.60	4.5
N6	24/5/14	22/6/14	83	1.49	0.04	0.64	1.25	0.30	0.29	4.4
N7	6/12/14	29/12/14	41	1.26	0.03	0.67	1.60	0.15	0.31	4.9
N8	26/3/15	2/4/15	30	1.26	0.03	0.97	1.45	0.36	1.76	4.1
N9	29/10/16	31/10/16	49	1.26	0.03	0.89	2.44	0.28	2.88	3.4
N10	26/1/17	28/3/17	568	1.29	0.04	0.73	1.31	0.36	1.00	5.1
N11	7/4/17	12/5/17	38	1.26	0.03	1.05	1.29	0.11	0.91	3.4
N12	12/6/17	8/8/17	614	1.48	0.12	0.79	1.46	0.25	0.86	6.4
N13	13/6/17	29/7/17	48	1.26	0.03	1.03	2.42	0.17	0.35	3.7
N14	15/8/17	23/10/17	38	1.26	0.03	1.06	1.13	0.39	0.26	3.5
N15	16/8/17	11/11/17	34	1.26	0.03	1.08	1.46	0.36	0.15	3.5
N16	17/8/17	8/11/17	39	1.26	0.03	1.24	2.39	0.14	0.31	3.2
N17	24/8/17	11/11/17	35	1.26	0.03	1.01	2.23	0.13	0.27	3.6

). Another case worth mentioning is the 24 May 2014, $M=6.9$, Samothraki seismic sequence which is divided into three major clusters ($N4,N5,N6$, in Figure 7). The estimated background rates of the three major clusters are relatively small ($\mu =0.16,0.60,0.29$), whereas the opposite holds for the scaling parameter, a, for the first two clusters ($a=1.82,1.76$). Concerning the seismic excitation that consists of clusters $N10,N11$ and $N14$, the relatively low productivity rates of the ETAS model ($a=1.31,1.29,1.13$) and, conversely, the relatively high background rates for the first two, $N10$ and $N11$, clusters ($\mu =1.00,0.91$) may indicate fluid intrusion. This observation is in accordance with the study of Mesimeri et al. [67] who derived high background rates after the estimation of the ETAS model to the empirically divided 5 sub-clusters of the primary seismic activity (January–March 2017). A fast-diminishing aftershock activity is observed for the main shock ($M=6.4$) that is located SE of Lesvos Island ($N12$), which is translated into a high Omori exponent, $p=1.48$. Additionally, low background rates characterize the three main clusters, $N12$, $N16$ and $N17$, indicating that they are probably related to tectonic and coseismic stress transfer from previous seismicity [62]. Worth mentioning are the remarkable high background rates for clusters $N8$ ($\mu =1.76$) and $N9$ ($\mu =2.78$), which could be an indicator for seismic activity driven by transient forces, however, the number of events is rather small and could have led to significant biases in the inversion of the parameters.

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 $J_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 (Figure 4, Figures 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 (Figure 5 and Figure S7) and offshore Perachora and Xylokastro (Figure 4 and Figure 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 ($\mu$) 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 $C6$, $C9$ and $C10$, Table 6, Figure 4 and Figure S3) in Corinth Gulf [63] and the 2017 Tuzla earthquake swarm (clusters $N10$, $N11$ and $N14$, Table 8, Figure 7 and Figure 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 $\mu >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 $\mu <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.

6. Conclusions

In this study, we present the efficiency of our clustering method, MAP-DBSCAN, on a simulated earthquake catalog where the structure of the clusters is known a priori and its competitiveness against well-known clustering algorithms, as in most cases, shows better results. The main seismic clusters in the Corinth Gulf, Central Ionian Islands and North Aegean Sea during 2012–2019 are detected by our method and their clustering properties are investigated. The results show the existence of regional variability in aftershock productivity and background rates. In particular, the Corinth Gulf is characterized by low productivity values and high background rates related to the dominance of earthquake swarms, whereas seismicity in the Central Ionian Islands is comprised by main shock–aftershock sequences with high productivity. Sequence-specific parameters verify the dependence between low productivity values and high background rates with pore-pressure due to fluids migration. We believe that future studies on Operational Earthquake Forecasting should incorporate localized parameters into the models to improve the forecasting accuracy.