## 4 Results
### 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(t)=c^{p-1}(p-1)(t-t_{i}+c)^{-p}\), and the productivity law that is expressed by \(N=k(M_{i})=Ke^{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\geq 30\) of each area. In particular, \(LL_{i}=\sum_{j=1}^{n_{i}}\log\lambda(t_{j})-\int_{l_{0}}^{t_{end}}\lambda(t)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,\ldots,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 \(LL_{i}\), and define as the common log-likelihood  
\[LL=\sum_{i=1}^{N^{*}}LL_{i}, \tag{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. There are 27 sequences in CG (Figures 4 and 5), 9 in CII (Figure 6) and 17 in NAS (Figure 7) from 2012 until 2019 with \(N\geq 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,\ldots,5\), then,  
Figure 7: (**a**) Spatial distribution of the centroids of the identified clusters for the area of North Aegean Sea along with the trace of North Aegean Trough (yellow lines). The size of the circles is proportional to the earthquake number in each cluster, whereas the duration is represented by the color scale. (**b**) Spatial distribution of the clusters with \(N\geq 30\) events. The index of each cluster is given in the inset box. we compute the corresponding \(LL_{i}^{new}\), \(i=1,\ldots,N^{*}\), log-likelihood values and store the new parameters under the condition \(LL_{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.