**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–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

**Citation:** Bountzis, P.; Papadimitriou, E.; Tsaklidis, G. Identification and Temporal Characteristics of Earthquake Clusters in Selected Areas in Greece. *Appl. Sci.* **2022**, *12*, 1908. https://doi.org/ 10.3390/app12041908

Academic Editor: Amadeo Benavent-Climent

Received: 31 December 2021 Accepted: 1 February 2022 Published: 11 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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 nonparametric definition of a cluster considering the space-time-magnitude proximity among earthquakes. Zaliapin and Ben-Zion [23] introduced a binary threshold, *η*, 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, *Nt*, with a piece-wise constant intensity rate determined by the underlying Markov process, *Jt*. 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 sugges<sup>t</sup> 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 sugges<sup>t</sup> 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–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 (*<sup>a</sup>*, *K*, *p*, *c* and *μ*) 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 (*Nt*, *Jt*)*<sup>t</sup>*∈*R*<sup>+</sup> , where *Nt* counts the number of earthquakes up to time *t* that occur with a rate *λ<sup>t</sup>*. Its value is associated with the unobservable states *i* = 1, . . . , *K*, of the Markov process *Jt*. In particular, when the process *Jt* is in state *i*, earthquakes occur according to a Poisson process with rate *λi* and, therefore, the sojourn time in this state follows an exponential distribution with expected value 1/*λi*. When an earthquake occurs, the MAP can transit with probability *pij* to another state *j*, so now, earthquakes occur according to a Poisson process with rate *<sup>λ</sup>j*, or remain in the current state *i* with probability *pii*.

To represent the MAP model, we need the *K* × *K* rate matrices **D**0 and **D**1, where **D**0 is a diagonal matrix whose non-negative elements we denote as *λ*1,..., *λK*, and which correspond to the *K* Poissonian rates, each one assigned to a hidden state of process *Jt*, and **D**1 consists of the transition rates among the states, along with the occurrence of an earthquake, which we denote as *qij*. 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, *λthr*, according to which a potential sequence starts when the rate of the counting process *Nt* achieves *λt* > *λthr* and ends as soon as the process *Jt* 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].
