The Mw = 5.6 Kanallaki Earthquake of 21 March 2020 in West Epirus, Greece: Reverse Fault Model from InSAR Data and Seismotectonic Implications for Apulia-Eurasia Collision

We identify the source of the Mw = 5.6 earthquake that hit west-central Epirus on 21 March 2020 00:49:52 UTC. We use Sentinel-1 synthetic aperture radar interferograms tied to one permanent Global Navigation Satellite System (GNSS) station (GARD). We model the source by inverting the INSAR displacement data. The inversion model suggests a shallow source on a low-angle fault (39◦) dipping towards east with a centroid depth of 8.5 km. The seismic moment deduced from our model agrees with those of the published seismic moment tensors. This geometry is compatible with reverse-slip motion along the west-verging Margariti thrust fault that accommodates part of the convergence within the collision zone between Apulia and Eurasia. We also processed new GNSS data and estimate a total convergence rate between Apulia and Eurasia of 8.9 mm yr−1, of which the shortening of the crust between the Epirus coastal GNSS stations and station PAXO in the Ionian Sea (across the Ionian Thrust) is equivalent to ~50% of it or 4.6 mm yr−1. By back-slip modelling we found that a 60-km wide deformation zone takes up nearly most of the convergence between Apulia-Eurasia, trending N318◦E. Its central axis runs along the southwest coast of Corfu, along the northeast coast of Paxoi, heading toward the northern extremity of the Lefkada island. The island of Paxoi appears kinematically as part of the Apulian plate.


Introduction
The tectonics of Epirus is characterised by on-going compression due to the active collision between the Apulian continental block [1][2][3][4] and Eurasian (Aegean) plate. The Apulian continental lithosphere subducts beneath Epirus and was imaged at 70-80 km depth by [5]. Both seismological and shows the study area in the NW corner of Greece. Panel (C) shows the relocated epicentre of the mainshock (this study; white star) and its aftershocks (red circles; period 20 March 2020-23 May 2020). Fault traces modified after [16,17]. The black line between the Paxoi and Corfu islands indicates the trace of the Ionian thrust (ticks on the overriding side).
Geological data (field mapping and borehole loggings) indicate that the region between Corfu and Ioannina ( Figure 1) is the site of Neogene thrusting and folding [2,3,12] involving mainly thinskinned thrusting onto the Apulian margin along Triassic evaporite décollements [13]. The rock formations in most of Epirus consist of Mesozoic and Cenozoic carbonate and clastic strata of the Ionian basin of the external Hellenides belt, which have been thrusted westwards onto the Apulian foreland. The inferred trace of the main basal thrust of Eurasia over Apulia, that is part of the Dinaric-Hellenic orogen, passes offshore Corfu-Paxoi islands with a general NNW-SSE orientation ( Figure  1A). In the region of Parga-Kanallaki (Figure 1), there are two onshore thrusts of the Ionian zone running N-S and are west-directed (i.e., the Margariti and Paramithia thrusts, [14]). Offshore Parga, two west-verging reverse faults are reported by [15] (mapped on an E-W oriented geological crosssection in that paper, see section B in their Figure 3) to offset Pliocene-Miocene sediments, the western one representing the contact of the overthrusted Ionian zone onto the pre-Apulian zone of the Hellenides (i.e., the Ionian Thrust, see Figure 1A). The Ionian zone is a well-known petroliferous area of Greece [3,13], yet public data on the details of the geological structure (including the detailed pattern of crustal shortening) are lacking. On the other hand, open data on the active tectonics of Epirus can be provided by Space Geodesy (InSAR and GNSS, [10,11]), which can be combined with published data from crustal seismicity surveys [6,9]. shows the relocated epicentre of the mainshock (this study; white star) and its aftershocks (red circles; period 20 March 2020-23 May 2020). Fault traces modified after [16,17]. The black line between the Paxoi and Corfu islands indicates the trace of the Ionian thrust (ticks on the overriding side).
Geological data (field mapping and borehole loggings) indicate that the region between Corfu and Ioannina ( Figure 1) is the site of Neogene thrusting and folding [2,3,12] involving mainly thin-skinned thrusting onto the Apulian margin along Triassic evaporite décollements [13]. The rock formations in most of Epirus consist of Mesozoic and Cenozoic carbonate and clastic strata of the Ionian basin of the external Hellenides belt, which have been thrusted westwards onto the Apulian foreland. The inferred trace of the main basal thrust of Eurasia over Apulia, that is part of the Dinaric-Hellenic orogen, passes offshore Corfu-Paxoi islands with a general NNW-SSE orientation ( Figure 1A). In the region of Parga-Kanallaki (Figure 1), there are two onshore thrusts of the Ionian zone running N-S and are west-directed (i.e., the Margariti and Paramithia thrusts, [14]). Offshore Parga, two west-verging reverse faults are reported by [15] (mapped on an E-W oriented geological cross-section in that paper, see section B in their Figure 3) to offset Pliocene-Miocene sediments, the western one representing the contact of the overthrusted Ionian zone onto the pre-Apulian zone of the Hellenides (i.e., the Ionian Thrust, see Figure 1A). The Ionian zone is a well-known petroliferous area of Greece [3,13], yet public data on the details of the geological structure (including the detailed pattern of crustal shortening) are lacking. On the other hand, open data on the active tectonics of Epirus can be provided by Space Geodesy (InSAR and GNSS, [10,11]), which can be combined with published data from crustal seismicity surveys [6,9]. Herein we use Sentinel-1 synthetic aperture radar interferograms, tied at the offsets measured at a nearby permanent Global Navigation Satellite System (GNSS) station (located at Gardiki, station GARD, see Figure 1A), to map the co-seismic displacement field of the 21 March 2020 Mw = 5.6 earthquake that occurred near Kanallaki, southwest Epirus ( Figure 1C). We relocated seismicity data to constrain the hypocentre of the mainshock and its moment tensor solution (Section 2). InSAR displacement data was used to invert for fault parameters (Section 3). Our inversion approach has resulted in the best-fitting fault model, providing its geometry and kinematics (strike/dip-angle/rake and concluding on its dip-direction). The geodetic observations and the seismicity data allow to discriminate the seismic fault between the two nodal fault planes of the double-couple focal mechanism of the mainshock. We then discuss our findings in terms of the geological structure (Section 4) and in particular the depth of the geodetic centroid in relation to the west-verging thrusts of the Ionian zone in that area. In this interpretation we rely on published geological data [1][2][3][4] as there are no new data (cross-sections or drilling data) investigating the subsurface of Epirus. We also analysed tectonic strain data from the GNSS velocities in Epirus aimed at estimating the rate of convergence and width of the deformation zone between Apulia and Eurasia. Our findings contribute to the understanding of the present tectonic processes in the Apulia-Eurasia collision zone.

The 21 March 2020 Earthquake-Relocation of Seismicity
On 21 March 2020 a Mw = 5.6 earthquake occurred at 00:49:52 UTC near the village of Kanallaki in Epirus ( Figure 1C), 15 km east of Parga, in north-western Greece [18]. The earthquake caused damage to old houses in Kanallaki and surrounding villages ( Figure 1) but no major injuries. The mainshock was preceded by an earthquake with a magnitude of Mw = 4.2 on 20 March 2020, 21:38 UTC (see Table 1 for moment tensor details). Until 23 May, 2020 more than 280 aftershocks (with 1.0 ≤ M L ≤ 4.0) were recorded by the Institute of Geodynamics, National Observatory of Athens (NOAGI) of which 270 events were relocated for the purposes of this study, with horizontal location uncertainties below 1 km (see Figures 1C and 2). The mainshock was followed by two small-size aftershocks, one on 23 March 2020 04:41 UTC and one on 25 March 2020 09:49 UTC with magnitudes of Mw = 3.8 and Mw = 3.9, respectively ( Table 1). The moment tensor solutions of all four events indicate NW-SE reverse faulting in agreement with the regional, compressional tectonics (Table 1). For the mainshock, two moment tensor calculations are included in Table 1. The first one is based on the rapid manual solution (Figure 3a) provided by NOA a few minutes after the occurrence of the earthquake. The second one is based on the relocated determination of its hypocentre for the needs of the present study ( Figure 3b). The calculation of the focal mechanism solutions is based on a waveform inversion methodology employing ISOLA point source inversion software ( [19], see Figure S1 for details on the quality of waveform fit). The crustal velocity model used for the relocated focal mechanism solution calculation was published by [20]. For seismicity relocation we collected catalogue and arrival-time data for the 2020 Kanallaki sequence from the online public databases of the Geodynamics Institute of the National Observatory of Athens (NOA) and the Seismological Laboratory of the National and Kapodistrian University of Athens (SL-NKUA). We merged the two databases by complementing arrival-time data for events that are common in the two catalogues, also including unique events only available on one of the two catalogues. In total, a database of 281 events was constructed for the period between 1 March and 1 May, 2020 in the study area. After several tests with different velocity models (e.g., [20]), we adopted the 1-D P-wave velocity model of [21] for the region of northwest Greece (Table 2), which produced adequate results in terms of travel-time residuals and focal depth distribution. We determined a constant ratio of Vp/Vs = 1.74 using the [22] method, considering also the root mean square (RMS) deviation curve with variable Vp/Vs, which presents a wide minimum for Vp/Vs between 1.71 and 1.75. The mean RMS error is 0.367 s (Figure 2a) while the average location errors are ERX = 0.54 km, ERY = 0.73 km, and ERZ = 2.7 km (Figure 2b). Although the azimuthal gap is relatively small for most events (mean gap = 121°, see Figure 2d), the epicentral distances from the recording stations are, in general, larger than 30 km (Figure 2e), which is why focal depths cannot be well constrained.  of 270 events were successfully relocated ( Figure 3). The relocated focal depths are mainly distributed between 3 and 11 km. The Mw = 5.6 mainshock's relocated focal depth is estimated at 6.4 km however, it should be noted that there are significant uncertainties in the absolute focal depths due to the unavailability of data from near-field stations that cannot be corrected by double-difference relocation. Soon after the Mw = 5.6 mainshock, the aftershock activity spread to a larger region further south and to the east (Figure 3). Cross-sections of seismicity are presented in supplementary Figure S2. The focal mechanism solutions produced by various agencies for the Kanallaki mainshock are reported in Table 3. The local agencies (NOA, AUTH) report centroid depths significantly lower than other agencies. Our revised moment tensor solution provided a centroid depth of 8 km. We favour the east-dipping nodal plane because it fits the regional geological data (west-verging thrusts, see Figures 1C and 4, [3,13,15]). We further justify our choice for a low-angle fault plane as the auxiliary For seismicity relocation we collected catalogue and arrival-time data for the 2020 Kanallaki sequence from the online public databases of the Geodynamics Institute of the National Observatory of Athens (NOA) and the Seismological Laboratory of the National and Kapodistrian University of Athens (SL-NKUA). We merged the two databases by complementing arrival-time data for events that are common in the two catalogues, also including unique events only available on one of the two catalogues. In total, a database of 281 events was constructed for the period between 1 March and 1 May, 2020 in the study area. After several tests with different velocity models (e.g., [20]), we adopted the 1-D P-wave velocity model of [21] for the region of northwest Greece (Table 2), which produced adequate results in terms of travel-time residuals and focal depth distribution. We determined a constant ratio of V p /V s = 1.74 using the [22] method, considering also the root mean square (RMS) deviation curve with variable V p /V s , which presents a wide minimum for V p /V s between 1.71 and 1.75. The mean RMS error is 0.367 s (Figure 2a) while the average location errors are ERX = 0.54 km, ERY = 0.73 km, and ERZ = 2.7 km (Figure 2b). Although the azimuthal gap is relatively small for most events (mean gap = 121 • , see Figure 2d), the epicentral distances from the recording stations are, in general, larger than 30 km (Figure 2e), which is why focal depths cannot be well constrained. To improve the hypocentral distribution, the double-difference relocation method was applied using the algorithm HypoDD [23]. This method works by assuming that when the distance between two neighbouring events is much smaller than their hypocentral distance from a station, then their travel-time differences can be attributed to their inter-event distance. Relative relocation between hypocentres is achieved by minimising the double-difference between observed and calculated travel-times for pairs of neighbouring events at common stations, reducing errors caused by an unmodeled velocity structure. To further enhance the relocation procedure, differential P and S travel-times from waveform cross-correlation was also incorporated. This creates links between strongly correlated events and reduces errors due to arrival-time reading inconsistencies. We acquired waveform data from stations of the Hellenic Unified Seismological Network from the National EIDA node (http://eida.gein.noa.gr/), hosted at the Institute of Geodynamics (NOA). The relocation procedure was performed using two sets of iterations, the first with stronger a priori weights on catalogue data and the second with higher weights on the cross-correlation data. A total of 270 events were successfully relocated ( Figure 3). The relocated focal depths are mainly distributed between 3 and 11 km. The M w = 5.6 mainshock's relocated focal depth is estimated at 6.4 km however, it should be noted that there are significant uncertainties in the absolute focal depths due to the unavailability of data from near-field stations that cannot be corrected by double-difference relocation. Soon after the M w = 5.6 mainshock, the aftershock activity spread to a larger region further south and to the east ( Figure 3). Cross-sections of seismicity are presented in supplementary Figure S2.
The focal mechanism solutions produced by various agencies for the Kanallaki mainshock are reported in Table 3. The local agencies (NOA, AUTH) report centroid depths significantly lower than other agencies. Our revised moment tensor solution provided a centroid depth of 8 km. We favour the east-dipping nodal plane because it fits the regional geological data (west-verging thrusts, see Figures  1C and 4, [3,13,15]). We further justify our choice for a low-angle fault plane as the auxiliary nodal plane of our revised Moment Tensor solution indicates a steep dip-angle of 59 • (see Supplementary Figure S1 and Table 3, columns NP2), which is mechanically unfavourable (based on our InSAR data). In addition, there are no mapped backthrusts (i.e., high-angle, west-dipping fault planes) in our study area (see Figure 4 for a map of active faults). Concerning the fault's strike, there are two clusters of solutions, one with NOA (initial), USGS, and GFZ around N317 • E and the others (GCMT, IPGP, AUTH, and INGV) around N336 • E. The average dip angle is 39 • towards east-northeast while the average rake angle is 103 • , i.e., almost pure reverse-slip.

Historical Seismicity and the 1895 Event
The Kanallaki region is an area of Greece where damaging earthquakes are rare during the last century and occurred with magnitudes M < 6.0 ( Figure 4, see Table S1 for the record 1915-2019) however, a disastrous earthquake took place in the 19th century. Both [24] and [25] report information on the occurrence of a major M = 6.  [24] reports, was the location of the heaviest damage, however contemporary sources (Greek newspaper "Voice of Epirus", 19 May, 1895, accessed on July 2020 through the digital archives of the Greek Parliament) and [25] mention the village Dragani (mod. Ampelia) as the meizoseismal area. According to the sources, the earthquake also destroyed the village Karvounari ( Figure 4) and damaged several others, killing Geosciences 2020, 10, 454 7 of 26 in total almost 100 people. We identified all settlements mentioned in the original sources with reported damages and plot them on a map (Figure 4), and an approximate location of the 1895 rupture can be inferred by the damage distribution (dotted red ellipse in Figure 4). Both the approximate location of the rupture and damage distribution are compatible with a M = 5.8 ± 0.3 magnitude reverse rupture along the Margariti thrust (east-dipping), as the 1895 damage zone is located on the hanging wall of this active fault. We document this finding by proper macroseismic intensity data analysis that can be found in Section 4.2 of this paper. The main shock was preceded by a strong foreshock on the previous day (see the catalogue of the Aristotle University of Thessaloniki at http://geophysics.geo.auth.gr/ss/CATALOGS/seiscat.dat). Geosciences 2020, 10, x FOR PEER REVIEW 7 of 26  [16,17]) and the epicentral area of the 1895 event (red ellipse). The instrumental seismicity data for Mw > 5.2 [26] is shown as solid circles. Dotted blue areas mark swamps and flooded fluvial areas during the early 20th century, while the used settlement names are the original ones mentioned in the 1895 earthquake reports.

Interferograms
We used synthetic aperture radar interferometry (InSAR) to capture the deformation produced by the earthquake (Figures 5 and 6). In the broader Aegean area, InSAR is systematically used to map the ground deformation produced by large earthquakes after removing the signal from the topography [27][28][29][30][31][32][33] and minimising the effect of decorrelation factors such as tropospheric noise [30,31].
We use the Sentinel-1 descending interferogram made with the 12 March and 24 March 2020 scenes acquired on track 153 ( Figure 5), and the ascending interferogram made with the 19 March and 31 March 2020 scenes acquired on track 175 ( Figure 6). The quality of the descending interferogram (track 153) is exceptional as it exhibits particularly high coherence and low tropospheric disturbances. The quality of the ascending interferogram is also good (in terms of coherence and tropospheric noise). Several interferograms of overlapping tracks were also processed  [16,17]) and the epicentral area of the 1895 event (red ellipse). The instrumental seismicity data for Mw > 5.2 [26] is shown as solid circles. Dotted blue areas mark swamps and flooded fluvial areas during the early 20th century, while the used settlement names are the original ones mentioned in the 1895 earthquake reports.
Moreover, on 16 June, 1990 02:16 UTC an offshore M = 5.5 event occurred, which was studied by [8] who determined a reverse focal mechanism. That earthquake was localised at 20.54 • E, 39.16 • N (see Figure 4), thus~15 km to the southwest of the 2020 event with a centroid depth of 8.5 ± 4 km, and strike, dip, and rake angles of 352 • ± 20 • , 37 • ± 9 • , and 110 • ± 15 • respectively, and a seismic moment of 2.49 × 10 17 N m. The focal mechanism and magnitude of the 1990 event are thus comparable with those of the 2020 earthquake.

Interferograms
We used synthetic aperture radar interferometry (InSAR) to capture the deformation produced by the earthquake (Figures 5 and 6). In the broader Aegean area, InSAR is systematically used to map the ground deformation produced by large earthquakes after removing the signal from the topography [27][28][29][30][31][32][33] and minimising the effect of decorrelation factors such as tropospheric noise [30,31].  The two interferograms capture the ground motion along two nearly-opposite lines of sight (l.o.s., see Look arrows in Figures 5 and 6). It is these measurements of l.o.s. that we are going to model to retrieve the parameters of the seismic fault. The digital elevation model (DEM) used for processing is the Shuttle Radar Topography Mission (SRTM) 1 Arc-Second Global (Digital Object Identifier number: /10.5066/F7PR7TFT). We enhanced the signal to noise ratio of the interferograms by applying the adaptive power spectrum filter of [34] with a coherence threshold of 0.3.
The two small aftershocks that occurred, one on 23 March with magnitude Mw(NOA) = 3.8, the other on 25 March with magnitude Mw(NOA) = 3.9 (see Table 1 and https://bbnet.gein.noa.gr/HL/) were not expected to affect the deformation signal of the mainshock. There was also a significant foreshock of magnitude Mw(NOA) = 4.2 on 20 March 21:38 UTC ( Table 1) that was also not expected to have produced any measurable surface deformation.  Figure 6). The quality of the descending interferogram (track 153) is exceptional as it exhibits particularly high coherence and low tropospheric disturbances. The quality of the ascending interferogram is also good (in terms of coherence and tropospheric noise). Several interferograms of overlapping tracks were also processed (see Figure S3 for ground coverage) however, the extracted information was of a lower quality (see Figure S4 for results of unused interferograms). The unit vectors from the ground to the satellites are [−0.098, 0.504, 0.858] and [−0.105, −0.569, 0.816] in the east, north, and up components for the descending and ascending views, respectively. The interferograms were produced on the Geohazards Exploitation Platform (https://geohazards-tep.eu) using the DIAPASON module developed by TRE-ALTAMIRA.
The two interferograms capture the ground motion along two nearly-opposite lines of sight (l.o.s., see Look arrows in Figures 5 and 6). It is these measurements of l.o.s. that we are going to model to retrieve the parameters of the seismic fault. The digital elevation model (DEM) used for processing is the Shuttle Radar Topography Mission (SRTM) 1 Arc-Second Global (Digital Object Identifier number: /10.5066/F7PR7TFT). We enhanced the signal to noise ratio of the interferograms by applying the adaptive power spectrum filter of [34] with a coherence threshold of 0.3.
The two small aftershocks that occurred, one on 23 March with magnitude M w(NOA) = 3.8, the other on 25 March with magnitude M w(NOA) = 3.9 (see Table 1 and https://bbnet.gein.noa.gr/HL/) were not expected to affect the deformation signal of the mainshock. There was also a significant foreshock of magnitude M w(NOA) = 4.2 on 20 March 21:38 UTC ( Table 1) that was also not expected to have produced any measurable surface deformation.

Geodetic Data Inversion: Picking Quarter of Fringes
On the two above interferograms we extracted numerically by manual picking the pixels constituting the fringe with values 32, 96, 160, and 224. On the descending interferogram (see Table 4) we selected 86 values among 247, and on the ascending interferogram (see Table 5) 86 among 289. As expected, there are less pixels for the larger displacements as this corresponds (by definition) to smaller areas. On the descending interferogram ( Figure 5), the zero value of the interferogram is very well constrained due to the high coherence of the interferogram. We assume that the zero value is the average zero counted in the triangle Parga-Igoumenitsa-Paramithia (see Figure 4 for locations). This assumption is justified by the small offset in the position time-series of the GNSS station GARD (Gardiki) which is located at the SE side of the triangle (see Supplementary Figure S5 for time series of this station and see Figures 1 and 6 for its location). We measured a −3.2 ± 0.2 mm static offset on the E-W component of station GARD (minus sign indicates displacement towards the west) which is compatible with the assumed co-seismic motion along a low-angle thrust towards west, as indicated by nodal plane 1 of the moment tensor of the mainshock ( Table 3). The InSAR uncertainty is around 1/8 of a fringe, i.e., 3.5 mm in the l.o.s. Regarding InSAR data uncertainty, there are two contributions, the first being linked to the picking of the fringes, but in our case, as we extracted only pixels with a pre-defined value this uncertainty is zero. The second is related to the tropospheric noise that is present in the scene. In our case, especially for the ascending interferogram, this noise is very low. It can be assessed visually and this was confirmed by our group. Table 6 shows how the output of our modelling changes when changing by 1/8 of a fringe the zero of our fringes which is equivalent to changing of the same amount the maximum l.o.s. displacement, here 28 to 31.5 mm. This test is made for a descending interferogram case and with strike, dip, and rake angles of 329 • , 39 • , and 103 • (i.e., the average values for NP1, see Table 3). The test shows that changing the zero of the interferogram by 1/8 of a fringe (3.5 mm) results in the change of the moment tensor by 32%, here completely absorbed by the increase of the fault length as this is the only free geometric parameter, the other free parameters being the three coordinates of the fault-top centre. As our first models predicted a geodetic moment slightly lower that the seismic moment (see Table S4), we use for the zero of the fringes the value that maximises the interferogram, thus the case of maximum line of sight change is 31.5 mm.    The ascending interferogram ( Figure 6) is noisier and therefore its zero more difficult to establish. We therefore used a different method and performed the fringe tuning by aligning the geodetic moment predicted for the descending interferogram with the one predicted for the ascending one. The best fit is obtained when using 34 mm for the displacement on the top quarter of fringe of the ascending interferogram. The fault model is now made with joint inversion of the l.o.s ground motions measured with descending and ascending Sentinel-1 interferograms. Our geodetic data provided us with a robust model and in particular on (a) the 3D location of the fault-top centre, (b) fault length and width, and (c) amount of slip.

Results of Inversion Modelling
Assuming a half-space elastic model with uniform slip along a rectangular fault surface, the source of the ground deformation was inverted using the InSAR data and the code inverse6 [35]. For the modelling of the interferograms, the same weight to the ascending and descending data was given and the same zero value for the fringes. The question of the tuning of the zero values has already been explained in the previous section. The equal weight is obtained by using the same number of ascending and descending observations (fringe pickings), here 86 for both, thus 172 observations in total (all 172 data points are compared against the model in Table S2). The mean scatter between data and model is 3.3 mm for the 86-descending data (Figure 8 top) and 3.9 mm for the 86-ascending data (Figure 8 bottom). We assume that the seismic fault plane is the one with low-angle reverse fault dipping towards ENE. This is the most consistent with the geomorphology and tectonic context of the area. Concerning the hypothesis of the antithetic plane, in our paper this hypothesis is not rejected primarily on the basis of the interferograms but on the basis of geological and seismological considerations [2,3,7,8,12,13,15]. As the fault is not very shallow, the magnitude is only M w = 5.6 and the dip angle is not very different from 45 • , the two conjugate planes produce modelled interferograms that are very similar to each other, therefore InSAR does not have the capability to discriminate. The 3D coordinates of the fault-top centre and for the fault length are inverted. We do not invert for the strike, dip, and rake angles, assuming 329°, 39°, and 103° respectively, i.e., using the average values of the various determinations of the focal mechanism parameters (see Table 3). The dip angle is, among the three angles, the one that is less constrained. For the two other angles there is a tradeoff between strike and rake. We tested models with the strike and/or the rake free but we could not constrain well those parameters. This was expected given the almost circular shape of the interferograms, especially the ascending one ( Figure 6). Indeed, the near-circular shape of the InSAR deformation pattern does not favour any angle in the azimuth and poorly constrains the dip angle as well. We may have tried to constrain the rake, but we know that in the modelling of InSAR (and also GNSS), there is a strong trade-off between rake and azimuth angles. Therefore we considered that it was more rigorous to rely on the angles provided by seismology through the focal mechanisms. The 3D coordinates of the fault-top centre and for the fault length are inverted. We do not invert for the strike, dip, and rake angles, assuming 329 • , 39 • , and 103 • respectively, i.e., using the average values of the various determinations of the focal mechanism parameters (see Table 3). The dip angle is, among the three angles, the one that is less constrained. For the two other angles there is a trade-off between strike and rake. We tested models with the strike and/or the rake free but we could not constrain well those parameters. This was expected given the almost circular shape of the interferograms, especially the ascending one ( Figure 6). Indeed, the near-circular shape of the InSAR deformation pattern does not favour any angle in the azimuth and poorly constrains the dip angle as well. We may have tried to constrain the rake, but we know that in the modelling of InSAR (and also GNSS), there is a strong trade-off between rake and azimuth angles. Therefore we considered that it was more rigorous to rely on the angles provided by seismology through the focal mechanisms. Moreover, there is a strong compatibility between the azimuth angle and orientation of the tectonic structures (see active faults in Figures 5 and 6).
In a first series of inversions we find that the width of the fault tends to become large, between 4.5 and 6 km and the length small, between 4 and 5 km, with the surface (rupture) area remaining almost constant. As the width is not well constrained, it was fixed at 4.7 km, a value the leads to an equal length of 4.7 km, thus a square fault. The amplitude of slip and rake was also released and it was found that the amplitude had the tendency to grow (i.e., from an initial 0.3 m to 0.6 m or more) and the rake to remain steady. In the second set of inversions, we constrained the slip to be 0.5 m, as the standard scaling laws (i.e., [36]) do not suggest a higher value for this earthquake, given also the size of the fault. The best fitting values for the four free parameters are reported in Table 7. Concerning the location of the fault, the horizontal coordinates are well constrained by InSAR (centre of top-fault edge) and the vertical coordinates are also well constrained by the slope of the interferograms. The three parameters that remain unknown are the fault length, width, and amount of slip. The fault length is well constrained by the spatial extension of the interferogram. This leaves two parameters remaining: The width of the fault and amount of slip. Like for azimuth and rake, there is also a strong trade-off between these two parameters, therefore it is a bit difficult to estimate a best fitting solution. We used the estimates on the centroid depth from seismology (6.4-8 km) to estimate approximately the width and it was assumed that the fault could not have a width significantly larger than its length. Under those assumptions two remaining parameters were fixed. The variable slip was not tested because the earthquake size (Mw = 5.6) is too small for that and its deformation is almost similar to the one of a point source. Non-planar fault solutions were not tested for the same reason.   Table 7). The geodetic seismic moment, assuming a rigidity of 3 × 10 10 Pa for the upper crust is 3.31 × 10 17 N m. This is within the range of values found from seismology (Table 3) and close to the NOA-revised (2.86 × 10 17 N m; Table 3) so it is consistent. The slip inferred from our model (0.5 m) is twice larger than the theoretical slip predicted by the scaling relationships of [37] (0.24 m). Those scaling laws predict a fault length of 5.5 km and a fault width of 5.8 km, respectively. In this case, the length is very well constrained by the fringes and cannot grow above a maximum of 5 km (4.7 km in our best fitting model) while the width, which is not very well constrained, could grow up to 6-6.5 km. Accepting such an unusual elongated (along-dip) shape for the fault would allow reducing the amplitude of the slip. For example a fault of 4.5 × 6.5 km and slip of 0.38 m would fit the interferogram equally well with better consistency with scaling laws. We prefer to retain a square geometry for the seismic source primarily because (a) it fits rupture-size data [37] for sources with M~5.5, (b) it comprises a continental thrust with a dip-angle of~39 • which is compatible with regional geological data. And (c) the map view of aftershocks (Figure 3b, Supplementary Figure S2) does not indicate a NE-SW distribution which would be in favour of an increased fault-width. Figure 9b presents, with the same horizontal and vertical scale, a 19-km-long cross section of the Kanallaki area running from the shore towards the interior, perpendicular to the strike of the fault (passing through the geodetic centroid). In this section we show the topography in blue and the predicted vertical deformation (surface uplift) in red, as deduced from our fault model ( Table 7). The maximum value of the vertical displacement is +4 cm at the abscissa 9.71 km, i.e., almost on top of the location of the upper tip of the fault (abscissa 8.9 km), that is on its hanging wall side. A similar surface displacement maximum (uplift) was mapped by InSAR for the 2019 Durres earthquake, as well [33]. Figure 10 shows a simplified geological map and a simplified structural section of the crust beneath Kanallaki. The 21 March 2020 earthquake has presumably occurred on one of the two main   (Table 7). The maximum value of the vertical displacement is +4 cm at the abscissa 9.71 km, i.e., almost on top of the location of the upper tip of the fault (abscissa 8.9 km), that is on its hanging wall side. A similar surface displacement maximum (uplift) was mapped by InSAR for the 2019 Durres earthquake, as well [33]. Figure 10 shows a simplified geological map and a simplified structural section of the crust beneath Kanallaki. The 21 March 2020 earthquake has presumably occurred on one of the two main reverse-slip structures of Margariti or Parga, outcropping in the area. The thrust fault of Margariti (code GR0695 on the NOAFault database, [38], https://arcg.is/04Haer) continues south towards Loutsa (Figure 9a, following [39]) and possibly offshore. It may root into the basal decollement of the Ionian zone over the crystalline basement. There are also Triassic evaporitic outcrops near Margariti (Figure 10a) and Loutsa that imply the existence of a major thrust there. As the Margariti thrust fault (see Figure 10a) is more than 12 km long its down-dip extension is comparable to its length so it could reach beneath the Kanallaki area (Figure 10b). Indeed, Kanallaki is located~10 km to the east of the trace of the Margariti thrust and such an east-dipping fault geometry also complies with a geodetic centroid depth of 8.5 km that we obtained for the 2020 seismic fault.  [38], https://arcg.is/04Haer) continues south towards Loutsa (Figure 9a, following [39]) and possibly offshore. It may root into the basal decollement of the Ionian zone over the crystalline basement. There are also Triassic evaporitic outcrops near Margariti (Figure 10a) and Loutsa that imply the existence of a major thrust there. As the Margariti thrust fault (see Figure 10a) is more than 12 km long its down-dip extension is comparable to its length so it could reach beneath the Kanallaki area (Figure 10b). Indeed, Kanallaki is located ~10 km to the east of the trace of the Margariti thrust and such an east-dipping fault geometry also complies with a geodetic centroid depth of 8.5 km that we obtained for the 2020 seismic fault.  [14,16,40]. Red diamonds mark the locations of co-seismic landslides for the 21 March 2020 earthquake, identified using Sentinel-2 optical imagery (see Table S5 for coordinates). To the northwest of Kanallaki and west of village Gliki (Figure 10a), the Lippa mountain forms a ~400 m high N-S anticline composed of hard carbonate rocks of the Ionian zone. The Kokitos and Acherontas rivers reach the sea southwest of this anticline (Figure 9a, see Figure S6 for an aerial view, other sub-parallel anticlines also exist). Whether the Lippa mountain structure is part of the active thrust system that roots at the decollement of 7-8 km depth requires further investigations (seismic profiles and drilling). In such a case, the lateral growth of the anticline would be balanced by erosion removed by the  [14,16,40]. Red diamonds mark the locations of co-seismic landslides for the 21 March 2020 earthquake, identified using Sentinel-2 optical imagery (see Table S5 for coordinates). To the northwest of Kanallaki and west of village Gliki (Figure 10a), the Lippa mountain forms a~400 m high N-S anticline composed of hard carbonate rocks of the Ionian zone. The Kokitos and Acherontas rivers reach the sea southwest of this anticline (Figure 9a, see Figure S6 for an aerial view, other sub-parallel anticlines also exist). Whether the Lippa mountain structure is part of the active thrust system that roots at the decollement of 7-8 km depth requires further investigations (seismic profiles and drilling). In such a case, the lateral growth of the anticline would be balanced by erosion removed by the two rivers. On-going tectonic activity is also evidenced by the uplift of the Epirus coastline due to the Late Tertiary convergence between the Apulian continental block and Eurasia [12]. the macroseismic epicentre from the SHEEC catalogue [43]. Then, using new intensity sites from the aforementioned sources and also revising some of the intensity values reported by [24,25], we provide a new list of macroseismic intensities in Supplementary Table S3. A detailed map of the updated macroseismic intensities is presented in Supplementary Figure S7. In this work we used the default IPEs from [44], for southern France, and [45], for metropolitan France and high-attenuation area that came along with the Quake-MD code [41]. With an intensity threshold of Ic = 6 and maximum intensity Io = 10, we determined a range of magnitudes M = 5.8 ± 0.3 for shallow focal depths H~4 km (Figure 11a,b). Other macroseismic epicentres reported in the AHEAD database were tried out, producing similar results concerning M and H, but not a better fit than with the SHEEC epicentral solution.

Revision of the Magnitude and Location of the 1895 Earthquake in West Epirus
In order to assess the magnitude and depth of the 1895 event in Dragani, we employed the QUake-MD open-source code [41]. This algorithm inverts intensity data points for a fixed macroseismic epicentre with a set of different Intensity Prediction Equations (IPE) and determines a plausible range of magnitudes and focal depths that better match the observations. We first acquired macroseismic data from the AHEAD database [42] (https://www.emidius.eu/AHEAD/event/ 18950514_0000_000, accessed on 14 September, 2020) and the macroseismic epicentre from the SHEEC catalogue [43]. Then, using new intensity sites from the aforementioned sources and also revising some of the intensity values reported by [24,25], we provide a new list of macroseismic intensities in Supplementary Table S3. A detailed map of the updated macroseismic intensities is presented in Supplementary Figure S7. In this work we used the default IPEs from [44], for southern France, and [45], for metropolitan France and high-attenuation area that came along with the Quake-MD code [41]. With an intensity threshold of Ic = 6 and maximum intensity Io = 10, we determined a range of magnitudes M = 5.8 ± 0.3 for shallow focal depths H~4 km (Figure 11a,b). Other macroseismic epicentres reported in the AHEAD database were tried out, producing similar results concerning M and H, but not a better fit than with the SHEEC epicentral solution.
The calculations with intensity observations (see Table S3) and an intensity threshold of Ic = 4, with maximum intensity Io = 11 were repeated considering several different cases for a macroseismic epicentre, including the ones reported in the AHEAD database and SHEEC, but also took into account the four sites for which a maximum intensity of Io = 11 was reported. The fit of the IPEs to the data was generally inadequate with the theoretical curves underestimating the intensity at short epicentral distances (<10 km) and overestimating larger ones. This may be due to the uncertain quality of intensity data points at intermediate and farther distances, while some high intensity observation could be attributed to local site conditions, amplifying the peak ground motion. In Figure 11c,d we used the site of Dragani, in the meizoseismal area (Io = 11) as a macroseismic epicentre. In that case, the plausible magnitude range was found at M = 5.6 ± 0.2 but the resulting focal depth, derived from the inversions, was very shallow (H~1 km, Figure 11c,d). Please note that Dragani and Karvounari, which were both completely destroyed (Io = 11) by the 1895 event, are located near the trace of a west-verging thrust, parallel to the Margariti active fault (Figure 4). To summarise this discussion, due to lack of data (such as geo-environmental effects) on this historical event we find it more reasonable to attribute the 1895 event to the Margariti fault as well.
The occurrence of two M5.6 + events (1895 and 2020) along the 38 km long Margariti fault is exceptional for Epirus and for the seismic behaviour of the active thrusts of the Ionian zone of the Hellenides. The two earthquakes likely ruptured different asperities along strike and this behaviour may indicate some sort of systematic rupture pattern of this fault and possibly of similar thrust faults in Epirus. The near 100 year earthquake record in this area (see Figure 4) is far from complete, however it is notable that no events with M > 6.1 have been recorded. Table 8. Velocities of several permanent GPS stations located in north-western Greece (Figure 1b) with longitude between 19 • -21 • and latitude between 39 • -40 • . The velocities are plotted in Figure 14 (with respect to rigid Apulia) and Figure 12

Tectonic Strain and the Determination of the Deformation Zone across the Apulia-Eurasia Plate Boundary
New GNSS data of the Greek permanent stations in Epirus and north Ionian Sea (Corfu and Paxoi islands, Table 8) were analysed. The GNSS stations belong to the NOANET geodetic network [46,47], the Hexagon SmartNet and Uranus commercial network (see Figure 1 for locations). The processing was made with the Precise Point Positioning (PPP) strategy [48] by means of the GIPSY/OASIS II software (ver. 6.4) developed by the Jet Propulsion Laboratory (JPL; http://gipsy-oasis.jpl.nasa.gov, [49]). The JPL final orbits (flinnR) and clocks, absolute antenna calibrations, random walk troposphere estimation, and the FES2004 ocean loading model were used. An in-depth description of similar GPS data processing can be found in [31][32][33]. The tectonic velocities and uncertainties retrieved from the analysis of the time series are given in Table 8 and were used to calculate plate convergence ( Figure 13).

Conclusions
1. The main source of the Mw = 5.6 earthquake that hit western-central Epirus on 21 March 2020 was identified to be located on the Margariti thrust fault, within the frontal area of the Ionian fold and thrust belt of the Hellenic orogen; 2. The seismic fault was modelled by combining the ascending and descending Sentinel-1 observations. It was found that we could model the overall fringe pattern by reverse slip on an east-dipping fault. The fault plane is a low-angle thrust fault (~5 by 5 km) that dips 39° towards east; 3. The inversion of geodetic data suggests that the upper edge of the fault is at a depth of 7 km, well constrained by the modelling of the interferograms; 4. The geodetic centroid of the slip plane is located at 8.5 km depth which is in agreement with our revised moment tensor solution that provided a centroid depth of 8 km; 5. InSAR showed ground motion towards the southwest and surface uplift in agreement with moment tensor solutions from seismology; 6. The damage distribution of the 14 May, 1895 M5.8 ± 0.3 damaging earthquake was found to be compatible with a reverse-slip rupture along the Margariti thrust; 7. Deep crustal seismicity near Ioannina is of a compressional type and is due to ongoing crustal shortening between Apulia and Eurasia; 8. GNSS data indicate that the extent of the deforming crustal area between Apulia-Epirus (from offshore Paxoi until Kanallaki), accommodating 80% of the total convergence, is 60 km, with a best fitting locking depth of 11 km; 9. The shortening rate of the crust between the Epirus coastal GNSS stations and station PAXO in the Ionian Sea (i.e., across the Ionian Thrust) is equivalent to 4.6 mm yr −1 or ~50% of the overall convergence; 10. The island of Paxoi appeared to already be on the Apulian side of the deforming zone, while it deforming area (box in Figure 12), when defined as the band that accommodates 80% of the total convergence is 60 km with a best fitting locking depth of 11 km. Unfortunately, the limited surface displacement data were not dense enough to allow for a formal inversion so as to recover the slip deficit patterns along the deformation front. Therefore, PAXO appears to be located on the Apulian side of the deforming zone, while it is exactly the opposite for the stations onshore Corfu (KASI, KERK, and KERU), with the axis of the plate boundary zone with Africa passing close to both islands and continuing along the coast of Albania. Between Corfu and Epirus (the Igoumenitsa stations IGOU-HGOU) the data and model ( Figure 12) show that there is no significant shortening within the velocity uncertainties of ± 0.2 mm   Figure 13) had a depth of 14 km, while maintaining similar kinematics. The earthquake occurred 20 km to the NW of Ioannina, Epirus, NW Greece, at an area where microseismicity had already been detected by [6]. The main shock was followed by a few hundreds of aftershocks in the period up to 25 October 2016. The hypocentre (14 km depth) is located at the bottom of the pile of deformed carbonates of the Ionian zone of the Hellenides and close to the crystalline basement as did the Durres earthquake in 2019 [33]. Due to concern for Ioannina, Sentinel differential interferograms were made immediately after that event but no deformation could be seen. This is because the moderate-size event was too deep. In addition, GPS data processing showed no displacement at the two stations of Ioannina ( Figure 1). Analysis of the moment tensor solution of the earthquake and its stronger aftershocks ( Figure 13) and map distribution of aftershocks revealed seismic slip along a shallow, east-dipping (27°) reverse fault, striking NNW-SSE [55], thus confirming the shortening regime suggested by the GNSS data. Moreover, the 2016 reverse-slip earthquake signifies the change in the stress-field near the meridian 20.9° E from ~N-S extension (Pindos mountains) to ~E-W shortening near Ioannina [11]. The deep setting of the thrust fault inside the Ionian fold-and-thrust belt of the external Hellenides denotes the ongoing accommodation of brittle deformation about 80 km to the east of the collision front with Apulia ( Figure 14). We therefore suggest that deep crustal seismicity near Ioannina is of a compressional type and it is due to ongoing crustal shortening between Apulia and Eurasia, as documented by the GNSS data ( Figure 12). The station velocities near Arta (Figure 1) also reflect the ongoing N-S extension of the Amvrakikos graben [11,50,51], so they will be left out in further analysis. Another notable exception in this kinematic field is the Hexagon-SmartNet station PAXO ( Figure 1) whose E-W velocity is 24.3 ± 0.4 mm yr −1 (N-S velocity is 16.4 ± 0.4 mm yr −1 ) and was also considered in the strain analysis made by [11].
Using the velocities of Table 8 we calculate the convergence towards Apulia in the azimuth N228 • E which is the azimuth that fit best the vectors of the stations located close to the coast (thus excluding the two stations of Ioannina). Note that the orientation of the coast of Epirus is globally perpendicular to this azimuth. The adopted Apulian plate velocity (vE:25.1 mm/yr-vN:18.4 mm/yr) is derived from [11] by averaging the velocities of 10 GPS points located, with no ambiguity, on the Apulian platform in SE Italy. The values of convergence and their uncertainties are presented in Table 9. The convergence of PAXO towards Apulia is of 1.9 ± 0.4 mm yr −1 . For the three stations of Corfu (KASI, KERK, and KERU, Figure 1) the average convergence towards Apulia is 5.6 ± 0.4 mm yr −1 (in the azimuth N228 • E). For the two stations of Igoumenitsa (IGOU and HGOU) the average convergence is 5.8 ± 0.4 mm yr −1 , close to the value of Corfu. For GARD, the rate is 6.5 ± 0.45 mm yr −1 . For the two stations of Ioannina (IOAN and IOAU) the average convergence is 7.6 ± 0.4 mm yr −1 . As the net rate of convergence between Apulia and stable Europe is only 5.2 mm yr −1 (considering stations across the Adriatic Sea), we see that the convergence is amplified in northwest Greece by a factor of~1.5. This is due to the compression and clockwise rotation imposed by the westward component of movement of the Anatolian micro-plate towards Apulia together with dextral motion along the Cephalonia Transform Fault (Figure 1), a process that can be perceived until the north of Albania [11,12,52]. It is noteworthy that the shortening rate of the crust between the Epirus coastal GNSS stations and station PAXO in the Ionian Sea (i.e., across the Ionian Thrust; see Figure 14) is equivalent to 4.6 mm yr −1 or~50% of the overall convergence. By modelling this convergence data with a simple back-slip model in elastic half-space ( Figure 13) we find the best fit for a total convergence rate of 8.9 mm yr −1 , thus an extra compression of 3.7 mm yr −1 is recovered in the regions located to the east and north of Ioannina. For back-slip modelling, we used the programme described (and distributed) in [35], simulating the locked interface like in [29], except that in this case the movement is reverse instead of strike slip. The modelled central axis of the convergence band, with azimuth N318 • E, passes along the southwest coast of Corfu, along the northeast coast of Paxoi, heading toward the northern extremity of the Lefkada island ( Figure 14). The 21 March 2020 earthquake occurred inside this band (Figures 13 and 14). The extent of the deforming area (box in Figure 13), when defined as the band that accommodates 80% of the total convergence is 60 km with a best fitting locking depth of 11 km. Unfortunately, the limited surface displacement data were not dense enough to allow for a formal inversion so as to recover the slip deficit patterns along the deformation front.
Therefore, PAXO appears to be located on the Apulian side of the deforming zone, while it is exactly the opposite for the stations onshore Corfu (KASI, KERK, and KERU), with the axis of the plate boundary zone with Africa passing close to both islands and continuing along the coast of Albania. Between Corfu and Epirus (the Igoumenitsa stations IGOU-HGOU) the data and model ( Figure 13) show that there is no significant shortening within the velocity uncertainties of ± 0.2 mm yr −1 , confirming the conclusions of [47].
In terms of strain, it is noteworthy that there is a 120 ns yr −1 active shortening of the crust between the two GNSS stations PAXO and GARD, located 38 km apart, and roughly in the ENE-WSW direction (N246 • E). This 1D estimate of tectonic strain is about 2.4 times larger than the 2D estimates from regional studies, e.g., the 50 ns yr −1 found by [11] possibly because of the strong smoothing effect of the regional strain calculations of [11]. In addition, both Ioannina GNSS stations (IOAN and IOAU) exhibit a long-term shortening of 41 and 51 ns yr −1 with respect to GARD, respectively (Figure 12), indicating ongoing compression onshore Epirus and to the east-northeast of the Kanallaki area, yet decreasing significantly when moving away from the collision front.  (Figure 14) had a depth of 14 km, while maintaining similar kinematics. The earthquake occurred 20 km to the NW of Ioannina, Epirus, NW Greece, at an area where microseismicity had already been detected by [6]. The main shock was followed by a few hundreds of aftershocks in the period up to 25 October 2016. The hypocentre (14 km depth) is located at the bottom of the pile of deformed carbonates of the Ionian zone of the Hellenides and close to the crystalline basement as did the Durres earthquake in 2019 [33]. Due to concern for Ioannina, Sentinel differential interferograms were made immediately after that event but no deformation could be seen. This is because the moderate-size event was too deep. In addition, GPS data processing showed no displacement at the two stations of Ioannina (Figure 1). Analysis of the moment tensor solution of the earthquake and its stronger aftershocks ( Figure 14) and map distribution of aftershocks revealed seismic slip along a shallow, east-dipping (27 • ) reverse fault, striking NNW-SSE [55], thus confirming the shortening regime suggested by the GNSS data. Moreover, the 2016 reverse-slip earthquake signifies the change in the stress-field near the meridian 20.9 • E from~N-S extension (Pindos mountains) to~E-W shortening near Ioannina [11]. The deep setting of the thrust fault inside the Ionian fold-and-thrust belt of the external Hellenides denotes the ongoing accommodation of brittle deformation about 80 km to the east of the collision front with Apulia ( Figure 12). We therefore suggest that deep crustal seismicity near Ioannina is of a compressional type and it is due to ongoing crustal shortening between Apulia and Eurasia, as documented by the GNSS data ( Figure 13).

1.
The main source of the Mw = 5.6 earthquake that hit western-central Epirus on 21 March 2020 was identified to be located on the Margariti thrust fault, within the frontal area of the Ionian fold and thrust belt of the Hellenic orogen; 2.
The seismic fault was modelled by combining the ascending and descending Sentinel-1 observations. It was found that we could model the overall fringe pattern by reverse slip on an east-dipping fault. The fault plane is a low-angle thrust fault (~5 by 5 km) that dips 39 • towards east; 3.
The inversion of geodetic data suggests that the upper edge of the fault is at a depth of 7 km, well constrained by the modelling of the interferograms; 4.
The geodetic centroid of the slip plane is located at 8.5 km depth which is in agreement with our revised moment tensor solution that provided a centroid depth of 8 km; 5.
InSAR showed ground motion towards the southwest and surface uplift in agreement with moment tensor solutions from seismology; 6.
The damage distribution of the 14 May, 1895 M5.8 ± 0.3 damaging earthquake was found to be compatible with a reverse-slip rupture along the Margariti thrust; 7.
Deep crustal seismicity near Ioannina is of a compressional type and is due to ongoing crustal shortening between Apulia and Eurasia; 8.
GNSS data indicate that the extent of the deforming crustal area between Apulia-Epirus (from offshore Paxoi until Kanallaki), accommodating 80% of the total convergence, is 60 km, with a best fitting locking depth of 11 km; 9.
The shortening rate of the crust between the Epirus coastal GNSS stations and station PAXO in the Ionian Sea (i.e., across the Ionian Thrust) is equivalent to 4.6 mm yr −1 or~50% of the overall convergence; 10. The island of Paxoi appeared to already be on the Apulian side of the deforming zone, while it is exactly the opposite for the GNSS stations onshore Corfu (KASI, KERK, and KERU, Figure 14), with the axis of the plate boundary zone between Africa passing close to both islands and continuing along the coast of Albania.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3263/10/11/454/s1, Figure S1: Results of the ISOLA modelling of the mainshock, Figure S2: Map and cross-sections of the relocated 2020 Kanallaki earthquake sequence, Figure S3: Map showing the Sentinel-1 frames with data used in this paper, Figure S4: Coseismic interferograms for the Kanallaki earthquake that were not used in the inversion due to low quality and noise, Figure S5: Graph of the GNSS time series of station GARD Figure S6: Aerial overview of the earthquake epicentre area, Figure S7: Detailed macroseismic intensity distribution for the 1895 earthquake, Table S1: Catalogue of strong events for the study area, 1915-2019, Table S2: Observed values (fringe picking) and modelled values, Table S3: Intensity observations for the 1895 earthquake near Dragani Table S4: Variability of the solution as a function of the focal mechanism parameters, Table S5: Location of co-seismic landslides identified using Sentinel-2 images for the Mw 5.6 21 March 2020 earthquake.
Funding: A.G. thanks "HELPOS-Hellenic Plate Observing System" (MIS 5002697) which is funded by the Operational Programme "Competitiveness, Entrepreneurship and Innovation" (NSRF 2014-2020) and co-financed by Greece and the European Union (European Regional Development Fund).