Earthquake Source Investigation of the Kanallaki, March 2020 Sequence (North-Western Greece) Based on Seismic and Geodetic Data

: The active collision of the Apulian continental lithosphere with the Eurasian plate charac-terizes the tectonics of the Epirus region in northwestern Greece, invoking crustal shortening. Epirus has not experienced any strong earthquakes during the instrumental era and thus there is no detailed knowledge of the way the active deformation is being expressed. In March 2020, a moderate size (Mw 5.8) earthquake sequence occurred close to the Kanallaki village in Epirus. The mainshock and major aftershock focal mechanisms are compatible with reverse faulting, on NNW-ESE trending nodal planes. We measure the coseismic surface deformation using radar interferometry and investigate the possible fault geometries based on seismic waveforms and InSAR data. Slip distribution models provide good ﬁts to both nodal planes and cannot resolve the fault plane ambiguity. The results indicate two slip episodes for a 337 ◦ N plane dipping 37 ◦ to the east and a single slip patch for a 137 ◦ N plane dipping 43 ◦ to 55 ◦ to the west. Even though the area of the sequence is very close to the triple junction of western Greece, the Kanallaki 2020 activity itself seems to be distinct from it, in terms of the acting stresses.


Introduction
On March 21, 2020 (00:49:51.8 UTC), a moderate-size earthquake sequence, with the main event Mw 5.8, occurred near the village of Kanallaki in Epirus, western Greece ( Figure 1). The seismic activity caused no fatalities but extended house damage in several villages, forced the local inhabitants to abandon their homes [1,2].
Epirus undergoes intense tectonic activity, which is evidenced in the topography, with a close alteration of high massifs and low valleys. Diapiric movements of Triassic evaporites contribute to the activation of reverse faults and to the development of a compressional regime [3]. Micro-earthquake studies in the Epirus area indicate a mixture of focal mechanisms [4][5][6]; this is also reflected by the presence of various types of tectonic structures: anticlines and synclines (NNW-SSE) [7], multiple thrusts, and overthrusts acting parallel (NNW-SSE) [7,8]. Reverse faults have a NW-SE strike, normal and oblique-normal faults have a NE-SW strike, and strike-slip faults have an E-W trend [3].
The 2020 seismic sequence occurred in a region where the Apulian platform converges with western Greece, and this active ENE-WSW compression is expressed mainly by thrust belts that trend NNW-SSE [9,10]. From a tectonic context, its location is significant because it occurred at the vicinity of a triple junction, which connects the compressional regime occurring to the north with an extensional and a strike-slip regime to the SE and SW, respectively ( Figure 1). More specifically, the Kanallaki thrust zone bifurcates to connect to the SE with the~E-W graben of the Amvrakikos Gulf, originated by extensional tectonics [11,12]. South of the Amvrakikos Gulf is the Katouna Fault [12,13], that operates as a sinistral strike-slip fault ( Figure 1). The other branch of the triple junction, SW of Kanallaki, is the Cephalonia-Lefkada Transform Fault Zone (CTF) [14,15], which shows strike-slip and thrust movement components and is the most seismically active area in Eastern Mediterranean. At the northernmost tip of the CTF, compressional and shear motions are acting simultaneously, creating a limited zone of transpressional tectonics [16] ( Figure 1). showing the triple junction of western Greece (light green notations). K-A fault notation is the Katouna-Amfilochia fault. The green dashed-line sketch of the Cephalonia-Lefkada fault Zone is based on [16]. Focal mechanisms (beachballs) indicate the representative kinematics at each area. Faults are from the GreDaSS fault database [17]. The inset shows the study area in a broader tectonic context. The Kanallaki 2020 sequence was extensively studied by [1]. They specifically reviewed the historical seismicity and the most important destructive earthquake in the region, the May 14, 1895, M~6 event, located to the north of Kanallaki village ( Figure 2). The dip direction of the causative fault of the mainshock is not known a priori. As mentioned in their study, geodetic inversion results cannot resolve the fault plane ambiguity; however, based on the geology of the area, they propose an east dipping geometry for the structure that caused the Kanallaki 2020 mainshock. Based on this assumption and their analysis, they proposed that the ruptured fault could be the Margariti thrust fault [18] ( Figure 2).

Figure 2.
Spatial distribution of the Kanallaki 2020 sequence (green dots) for the period March-November 2020. Beach balls are focal mechanisms of the sequence (green) alongside previous available solutions (grey) for earthquakes with Mw > 5.0. The location of the Mw~6 event of 1895 is also shown (asterisk). The Margariti thrust fault sketch is based on [1], and GRCS330 and GRCS377 faults are from the GreDaSS Fault database [17]. The inset shows the magnitude histogram of the sequence.
In [1], they present a source model assuming a uniform slip; we assume here a different model in order to investigate the complexity of the rupture process, exploiting geodetic but also seismic waveform data. The scope of this study is to focus on the rupture process expressed during the Kanallaki sequence in view of the broader seismotectonic context. To this end, we first use seismic data to invert for moment tensor solutions and examine the faulting type of the major earthquakes of the sequence. We also measure the coseismic surface deformation using the satellite InSAR (Interferometric Synthetic Aperture Radar) technique. We then investigate possible slip distribution source models, using InSAR surface deformation results and seismic waveform data. We compare these two individual datasets and present a series of possible slip distribution models that we consider all to be similarly plausible. Our goal is to have an insight about the possible style of rupture processes that could be expected in this poorly studied area. Finally, we estimate the principal stress orientations that represent the Kanallaki 2020 sequence to investigate the potential relation with the broader geotectonic framework of the plate collision occurring at Epirus.

Seismic Analysis
Digital three-component, full-broadband, and strong motion seismic waveform data from regional stations ( Figure S1) were exploited for the estimation of focal mechanisms and slip distributions. The data were downloaded from the ORFEUS EIDA node (see Data Availability section). The original waveforms were baseline corrected, tapered, corrected for the instrument response, converted to displacement, band pass filtered between 0.05 to 0.08 Hz for stronger events and 0.05 to 0.10 Hz for weaker events, and re-sampled to 1 Hz. Green's functions were computed by the frequency-wave number method and the 1D velocity model of Novotný et al. [19]. This model has been found to be an appropriate representation of average crustal properties for the Aegean area, as it effectively describes regional wave propagation in the Aegean and accounts for the characteristics of the waveforms in the low frequencies used. The synthetic waveforms (Green's functions) were also filtered similarly to the real data.
To retrieve the moment tensor solutions for Mw ≥ 3.3 events, we adopted the Time-Domain Moment Tensor inversion algorithm from the Berkeley Seismological Laboratory [20][21][22]. We constrain the inversion to solve for the best fitting deviatoric solution only, with no isotropic component. The double couple (DC) of the moment tensor is represented by the strike, dip, and rake of the two nodal planes. Source depth is found iteratively by finding the solution that yields the largest variance reduction (VR). In this approach, it is assumed that the location of the event is well represented by the high frequency hypocentral location; thus, a low frequency centroid location is not determined. Furthermore, the source time history is assumed to be synchronous for all the elements of the moment tensor, and it is approximated by a delta function.
To calculate the slip models from seismic data, we adopted a non-negative least squares solver [23,24]. This inversion scheme has proved to produce finite source parameters, which compare quite well with those obtained using local strong motion records. In the source model, the rupture propagates with constant speed over a grid of point sources, each with constant dislocation rise time, and the Green's functions are shifted in time to account for relative hypocentre-subfault-station distances and the time for the passage of a circular rupture front. The regional nature of the inversion procedure allows the simplification of the problem, considering only constant values of dislocation rise time and rupture velocity.
We also estimate the local stress filed expressed by the Kanallaki 2020 sequence, exploiting the here calculated and available moment tensor solutions. By using the STRESS-INVERSE code of Vavryčuk [25], we jointly inverted for the stress and fault orientation. It is assumed that there is no interaction between the earthquakes, they do not affect the background tectonic stress, the stress is uniform in the region, and the slip vector points are in the shear traction direction [26][27][28][29]. The algorithm is based on [30,31] and the instability criterion of [32]. Details and applications of the method can be found in [25,33,34].

Geodetic Analysis
For the measurement of the surface deformation caused by the seismic activity, we used InSAR [35,36]. Sentinel-1 satellite radar data (Table 1), from both the ascending and descending modes, were acquired from the European Space Agency (ESA) to apply the Differential InSAR (DInSAR) technique, using the SARscape software (sarmap, CH). We acquired images before and after the mainshock to calculate differential interfero-grams, showing the Line-of-Sight (LoS) ground displacement that occurred between the two dates. We used Goldstein filtering [37] and the minimum cost flow algorithm [38] for the phase unwrapping. The topographic correction was performed with the SRTM (Shuttle Radar Topography Mission) Digital Elevation Model (DEM) [39], while GACOS (Generic Atmospheric Correction Online Service for InSAR) [40][41][42] was used for the atmospheric corrections. For the InSAR data modelling, at first uniform slip non-linear inversion [43,44] was performed to estimate the best-fit parameters of a finite source dislocation embedded in an elastic half-space Okada [45] to predict the displacement over a large set of points sampled from all the available unwrapped interferograms; external constraints based on the focal mechanism parameters were also introduced to find the solution. Our goal was to get the slip distribution across the fault plane. Thus, after obtaining a uniform slip solution, we performed a linear inversion [46][47][48]. For the latter, we used the rupture mechanism and position derived from non-linear inversion, extending the uniform-slip source dimension to taper to zero the distributed slip. The linear inversion was carried out after subdividing the fault plane into 1 × 1 km subfaults and applying a further non-negativity constraint; to avoid back-slip, we used a trial-and-error approach to define the system damping, which is the empirical parameter balancing the data fit and the slip distribution roughness [49].

Results
The spatial distribution of the sequence and the focal mechanisms of the mainshock and largest events (Tables 2 and 3

Ref
Origin

Geodetic Source Modelling
During geodetic inversions, both the east and west dipping geometries gave the same rms values (Table S1); InSAR cannot discriminate between the two dipping geometries, as already mentioned by [1]. Thus, we report the InSAR slip distribution models for both cases.
The best fit geometric parameters derived from the non-linear inversion are: (a) West dipping plane: Strike/Dip/Rake = 137 • /43 • /109 • , and (b) East dipping plane: Strike/Dip/ Rake = 337 • /37 • /102 • . To estimate the uncertainties of the uniform-slip models, an empirical approach was followed; by adding realistic noise to the interferograms: a nonlinear inversion was applied again beginning from the configuration of the optimum parameters. We repeated this procedure 50 times and gathered results that represent the uncertainties and trade-offs of the parameters (Figures 4 and 5). Finally, two slip distributions were estimated: we carried out a linear inversion using the parameters presented above. The east dipping plane showed two main slip patches with a maximum value of 50 cm at~7 km and also some shallower slip ( Figure 6). The east dipping plane expressed a geodetic moment of 9.44·10 17 N m. For the case of the west-dipping plane, InSAR results show that there is only one prevalent slip patch with a maximum value of 40 cm, together with small portions of slip at shallow depths. The geodetic moment in this case is 6.96·10 17 N m. The maximum slip values are close to the uniform slip value adopted by [1] and the moment magnitudes are in accordance with the Global Centroid Moment Tensor Project (GCMT) estimated scalar seismic moment (4.78·10 17 N m). Both geometries can predict the original data well. The fitting of the InSAR predicted displacements with the observed ones is shown in Figures 7 and 8. The root mean square (rms) values of the linear inversion are almost the same for the two oppositely dipping planes (Table 4). InSAR source results of this study are reported in Table 5.

Seismic Source Modelling
We used 33 seismic waveform components to calculate slip models that best reproduce the data for both nodal planesm likewise with the InSAR modelling. For the west dipping case, we adopted the fault strike 137 • , dip 55 • and rake 78 • , while for the east dipping case, the parameters strike 337 • , dip 37 • , and rake 106 • (see Table 2, this study). The initial fault model had dimensions larger than the expected from empirical scaling relations to allow for bilateral rupture. Thus, a model of 25 × 25 km was used, discretized into 1 × 1 km subfaults. The dislocation rise time is assumed to be 10% of the approximate total duration, and in this parameterization, it was set equal to 0.5 s. The rupture speed was taken to be equal to 2.8 km/s (0.8 of Vs at the source depth).  The two slip distribution models derived from seismic waveforms are compared to the InSAR ones in Figure 6. Like with InSAR, when using the seismic data, both the east and west dipping geometries can predict the observed data equally well, with a variance reduction of~70% in both cases. Interestingly, the east-dipping geometry predicts a two-lobe slip pattern. For the west-dipping fault scenario, the slip is confined to a single patch. The resolved seismic moment for both fault parameterizations is of the order of 4.3·10 17 N m (Mw 5.7) in accordance with the MT modelling. The maximum slip amplitudes differ between the geodetic and seismic data: these values are shaped by the smoothing parameterization adopted within the inversions. In all cases though, the average slip from the seismic data is of the order of~5 cm, compatible with the values expected from the region empirical relations [50].

Estimation of the Stress Orientation
We exploited the Kanallaki seismic focal mechanisms to estimate the local stress field. The estimated principal stress axes (Figure 9) indicate that the maximum principal stress axis σ1 has an azimuth of 225 • with small plunge and the minimum principal stress axis σ3 is nearly vertical; the stress setting expresses the thrust faulting regime (Table 7; Figure 9). The σ1 orientation results are in a general accordance with the World Stress Map and with previous stress studies [51,52] and should be also considered to reflect the convergence direction between the Apulian platform and western Greece.   (Table 6). (d) Shape ratio histogram. Using as input the seismic focal mechanisms, the STRESSINVERSE algorithm proposes two optimally oriented faults (Table 6). It can be shown in the Mohr cycle diagram (Figure 9c) that most events of the sequence (blue crosses), except one, are situated in the lower hemisphere. The lower hemisphere represents an optimally oriented fault that is striking 340 • (optimal fault 2 in Table 6), which is indicating an east dipping fault.

Discussion and Conclusions
The moderate-size instrumental seismicity and the scarcity of historical major events hinder the detailed study of compressional tectonics of Epirus in northwestern Greece. We present here a source study of the Kanallaki March 2020 moderate-size sequence. Focal mechanisms denote reverse or low angle thrust faulting along NNW-ESE trending nodal planes.
The InSAR and seismic waveform inversion model results are compatible with both east-or west-dipping sources. Whatever is the case, the region is fabricated by closely spaced individual structures with variable geomorphic expressions. Even though most thrust sheets dip eastwards, an activation of a west dipping structure in the same stress field cannot be ruled out.
Regarding the slip distribution pattern, for both fault geometries there are a number of stable features: − There is a spatial correlation of the estimated main slip patches with the locations of the observed landslides. Both fault dip directions explain their relative location. − The Kanallaki village, where the major damage was observed, is close or exactly above the main slip patch in all the presented models. − The InSAR and seismic waveform slip distributions agree as to the slip pattern: A west dipping fault plane would have slipped in a single main patch, whereas an east-dipping fault would have expressed two main slip patches.
The Cephalonia-Lefkada fault (which is one branch of the triple junction of western Greece- Figure 1) is known to have originated major earthquakes that consisted of double or even triple sub-events [53]. The source time function of the Kanallaki 2020 mainshock, as reported by IPGP (http://scardec.projects.sismo.ipgp.fr/ (accessed on 9 October 2020)) [54], depicts two moment release episodes. Thus, it is possible that the Kanallaki mainshock was a double-source earthquake, which could mean the existence of two main asperities like those presented in our east-dipping models.
Moreover, this means that multiple-source events are plausible even at the northern branch of the triple junction of western Greece. This possible double-source indication is an interesting insight for the area of Epirus, also because the event was only of moderate size.
Even though the stress analysis was limited to the Kanallaki sequence, the orientation of the principal axes, that is, a horizontal σ1 and a nearly vertical σ3, are compatible with the activity of the lithospheric plate movements. The thrust character of the sources and the stress orientation could indicate that the Kanallaki 2020 sequence is the result of the plate collision. Moreover, it seems that the stresses expressed by the specific sequence are separate from the activity of the other parts of the triple junction, even though the area is very close to it.
Thus, our results indicate the possibility that the Kanallaki 2020 sequence has an indicative tectonic significance. At the same time, we suggest that a moderate seismic sequences should be each time considered with skepticism with regard to what extent they can be considered a reliable indicator of the regional stress pattern or of the regional tectonics. Future major seismic events are likely to shed more light regarding the complex tectonics of North-Western Greece.