Ground Deformation and Seismic Fault Model of the M6.4 Durres (Albania) Nov. 26, 2019 Earthquake, Based on GNSS / INSAR Observations

: We identify the source of the M w = 6.4 earthquake that rocked north-central Albania on November 26, 2019 02:54 UTC. We use synthetic aperture radar interferograms tied to the time series of coordinates of two permanent Global Navigation Satellite System (GNSS) stations (DUR2 and TIR2). We model the source by inverting the displacement data. Assuming in our model a half-space elastic medium and uniform slip along a rectangular fault surface, we invert the 104 picked measurements on a couple of ascending and descending interferograms to calculate the parameters of the fault. All inversions made with di ﬀ erent input parameters converge towards a stable and robust solution with root mean square (r.m.s.) residual of 5.4 mm, thus ~1 / 5 of a fringe. They reveal that the earthquake occurred deep in the crust on a low-angle fault (23 ◦ ) dipping towards east with a centroid at 16.5 km depth. The best-ﬁtting length and width of the fault are 22 and 13 km, and the reverse slip, 0.55 m. The seismic moment deduced from our model agrees with those of the published seismic moment tensors. This geometry is compatible with a blind thrust fault that may root on the main basal thrust, i.e., along the thrust front that separates Adria–Apulia from Eurasia. It is notable that there is a 123 ns yr − 1 active shortening of the crust between the GNSS stations DUR2-TIR2 (equivalent to a shortening rate of 3.6 mm yr − 1 ), and roughly in the east–west direction. Given this amount of strain the recurrence time of M6 + earthquakes along this fault should be of the order of 150 years.


Introduction
The tectonics of western and northern Albania are characterised by on-going compression due to collision between the Eurasian plate and the Adriatic block. Crustal deformation is characterised by shortening directed in a NNE-SSW to E-W orientation ( [1][2][3][4][5][6]; Figure 1a). Geological data indicate that the region between Tirana and Durrës (Figure 1b) is the site of Neogene thrusting and folding [7][8][9]. The main thrusts are west-directed [10,11] while the thrust front is dextrally offset along two NE-SW transfer zones, the "Scutari-Pec" [11,12]~60 km north of Tirana (near the town of Lezhe; Figure 1b) The November 26 M w = 6.4 event of Durrës was preceded by a M w = 5.7 that occurred on September 21, 2019 14:04 UTC, with the same type of reverse-faulting kinematics ( [21]; Figure 1b). This foreshock caused limited liquefaction effects near Durrës and damage to buildings of Durrës, Tirana and several settlements of the broader area [21]. The geometry and kinematics of both ruptures are similar. The M w = 6.4 event moment tensor solution (USGS; Table 2) shows a NNW-SSE fault plane with strike, dip and rake angles 338 • /27 • /92 • . The M w = 5.7 event moment tensor solution (USGS) shows a NNW-SSE fault plane with angles 323 • /32 • /93 • .
Here we use Sentinel-1 synthetic aperture radar interferograms, tied at the offsets measured at two permanent Global Navigation Satellite System (GNSS) stations, to infer the parameters of the seismic fault. Our best-fitting model is a low-angle reverse fault dipping at 23 • towards the northeast. The centroid of this fault is located at a depth of 16.5 km on a structure that could be a decollement at the base of the brittle crust. This finding contributes to the understanding of the present tectonic processes in the Apulia-Eurasia collision zone. Adriatic Sea bathymetry data from https://www.emodnet-bathymetry.eu/. Beachballs indicate the focal mechanisms of recent strong, shallow earthquakes [14]; USGS, Table 2; compressional quadrants are shaded in black). Shaded rectangle indicates study area shown in 1b. (B) Relief map of the Durres-Tirana area in central-western Albania showing main geological structures (thrust faults; ticks on the upthrown side), epicentre (yellow stars) determined by various agencies and Global Navigation Satellite System (GNSS) stations that recorded the earthquake. Beachball indicates the United States Geological Survey (USGS) focal mechanism of the Nov. 26, 2019 earthquake. Red star indicates the European Mediterranean Seismological Centre (EMSC) epicentre of the September 21, 2019 Mw = 5.7 earthquake (Focal mechanism from USGS). Fault lines are from [22]. Yellow line indicates the profile in Figure 6.  [22]. Yellow line indicates the profile in Figure 6.

InSAR Data Processing
We used synthetic aperture radar interferometry (InSAR) to capture the deformation produced by the Durrës earthquake ( Figure 2). In the Mediterranean, InSAR is systematically used to map the ground deformation produced by large earthquakes after removing the signal from the topography (e.g., [23][24][25][26][27][28][29][30][31]) and minimising the tropospheric noise. The red star indicates the earthquake epicentre (EMSC location), orange and green stars USGS and German Research Centre for Geosciences (GFZ) epicentres, respectively. The vector pairs represent the GPS measured (black) and modelled (white) co-seismic horizontal displacements at stations DUR2 and TIR2 respectively (see Table 3   The vector pairs represent the GPS measured (black) and modelled (white) co-seismic horizontal displacements at stations DUR2 and TIR2 respectively (see Table 3 Table 3).
We processed Sentinel-1 Interferometric Wide (IW) acquisitions from ascending and descending orbits, tracks 175 and 153, respectively. We used the open-source SNAP v7.0 ESA software [32], implemented in an offline version of Automated Interferogram Processing Station (AIPS, http: //aips.space.noa.gr) and the online facilities of the GEP-TEP (Geohazards Exploitation Platform -Thematic Exploitation Platform; http://geohazards-tep.eu). For the ascending track we use the data acquired on November 14 and 26 (16:32 UTC), and for the descending track those acquired on November 25 and December 1, 2019.
The two interferograms provided a measurement of the ground motion along two opposite line of sights (LOS; Figure 2a,b). The digital elevation model (DEM) used for the processing is the Shuttle Radar Topography Mission (SRTM) 1 Arc-Second Global ( [33]; Digital Object Identifier number:/10.5066/F7PR7TFT). We enhanced the signal to noise ratio by applying the adaptive power spectrum filter of [34] with a window size of 3, Fast Fourier Transform (FFT) size of 64 and a coherence threshold of 0.2.
The existence and amplitude of tropospheric disturbances in InSAR can be partly assessed by looking at single interferograms. This is because the tropospheric effect (signal delay) can be separated in two components: one correlated with the topography and the other not correlated. In the case of our interferograms ( Figure 2) we could find no correlation with topography. The estimated amplitude of the tropospheric noise at short spatial wavelengths (0.5-5 km) is low, below a quarter of a fringe. There are three fringes well visible in the ascending interferogram and two in the descending one, aligned in an area of~40 km elongated NW-SE north of Durrës, mostly on-shore along the coastal plain (Figures 1 and 2). The maximum ascending LOS change is~8.4 cm towards the satellite.
Four moderate magnitude aftershocks (up to M = 5.4; EMSC magnitudes https://www.emsc-csem. org/Earthquake/europe/M5/) occurred in the 36-h period following the mainshock; they did not affect the deformation signal as all those events occurred relatively deep (the GFZ moment tensor solutions for two aftershocks indicate depths comparable to that of the mainshock;~26 km) and with moderate magnitudes. We generated additional interferograms for the postseismic period (up to December 14, 2019), in order to clarify any possible impact on the estimated deformation. However, no postseismic deformation was detected (see supplementary Figure S2). We also processed Sentinel 1 interferometric pairs for the foreshock of September 21, 2019 (M w = 5.7) to investigate the source geometry and depth, but no ground displacements were detected either (see supplementary Figure S3).

GNSS Data Processing
We analysed the GNSS data of the stations DUR2 (Durrës, 19 Figure 1b for locations). The processing was made with the Precise Point Positioning (PPP) strategy [35] by means of the GIPSY/OASIS II software (ver. 6.4) developed by the Jet Propulsion Laboratory (JPL; http://gipsy-oasis.jpl.nasa.gov; [36]). We used the JPL final orbits (flinnR) and clocks, absolute antenna calibration, random walk troposphere estimation, and the FES2004 ocean loading model. We calculated the static offsets for both September 21, 2019 and November 26, 2019 events (Table 3; Table  3). Vertical lines indicate timing of earthquakes.  Table 3).

Inversion of the Geodetic Data
Vertical lines indicate timing of earthquakes.
The time series of coordinates are plotted in Figure 3 (a zoom on the week before and after the mainshock is presented in supplementary Figure S4). The long-term velocities, in the reference frame ITRF2014, are 23.6 ± 0.3 mm yr −1 in east and 16.3 ± 0.3 mm yr −1 in north at Durrës and 20.0 ± 0.3 mm yr −1 and 15.8 ± 0.3 mm yr −1 at Tirana. Therefore, the net motion of DUR2, assuming TIR2 fix, is 3.6 ± 0.5 mm yr −1 of shortening in the azimuth N80 • E. It is notable that there is a 123 ns yr −1 active shortening of the crust between the two GNSS stations located at 30 km one to the other, and roughly in the east-west direction. This 1-D estimate of tectonic strain is about 2.4 larger than the 2-D estimates from regional studies (e.g., 30-40 ns yr −1 , [5,6]). In addition, both GNSS stations exhibit long-term subsidence of −8.7 and −6.5 mm yr −1 respectively (Figure 3), presumably due mostly to anthropogenic reasons. The displacements and uncertainties retrieved from the analysis of the time series are reported in Table 3.
The LOS-projected (ascending and descending) motions measured at DUR2 for the November 26, 2019 event are used to give an absolute reference to the ascending and descending interferograms shown before (Figure 2). This absolute tie to GNSS is crucial for the inversion made in the next section. The LOS unit vector [East North Up] used is [−0.52 −0.12 0.84] and [0.63 −0.14 0.77], for the ascending and descending track, respectively.

Inversion of the Geodetic Data
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 [37]. We fed the inversion with 104 LOS measurements that were picked manually on the interferograms (Figure 2), 52 points on the ascending and 32 for descending track.
Taking various initial conditions compatible with the fault planes inferred from the seismological centres, we performed 6500 inversions. Those inversions define a stable and robust solution with root mean square (r.m.s.) residual of 5.4 mm thus~1/5 of a fringe. The parameters of the best-fitting solution are in the Table 4. Table 4 also contains the uncertainties on the parameters of our final best-fitting fault model. The inversion result for depth to top-fault vs. fault length is shown in Figure 4 which is tailored to appreciate visually the range of values (and thus uncertainties) of the parameter couple: top-fault depth-fault length. All inversion results are included in supplementary Figure S5. In Figure S6 it is shown the fit only at the picked points of the modelled and measured displacements in the LOS, along with the r.m.s. misfit, which is 3.3 mm and 8.1 mm for the ascending and descending track respectively.   Our modelling demonstrates that the earthquake occurred deep in the crust on a low-angle reverse fault (23 • ) dipping towards east with centroid at 16.5 km depth. The best-fitting length and width of the fault are 22 and 13 km and the reverse slip 0.55 m ( Figure 5). Assuming a medium rigidity of 3.3 × 10 10 Pa, the seismic moment deduced from our model, 5.19 × 10 18 N m −1 , is closer to GFZ's published seismic moment tensor ( Table 2). The modelling also indicates that the top of the fault plane is buried at a depth of~14 km, so there was no rupture in the shallow crust above that depth.
The geodetic fault-model is in agreement with published moment tensor solutions showing a NNW-SSE fault plane (see Table 2). The geodetic centroid is located~15 km to the NE of the EMSC epicentre but a few km away from either the USGS or the GFZ epicentre (Figure 1b; Figure 2). The inferred fault geometry (see Figure 5 for a surface projection) is compatible with a blind thrust fault rooting on the main basal thrust front that separates Adria-Apulia from Eurasia.

Geodetic Determination of Earthquake Parameters
Our modelling shows that the November 26, 2019 EMSC epicentre must be shifted by 15 km to match the geodetic data (see Figures 1b and 5 for locations; the same distance holds for epicentres determined by CMT and INGV; see Figure 1b). On the other hand, the USGS epicentre is located about 7 km NW of the geodetic centroid and the GFZ epicentre about 3 km to the SW, respectively. We note that such shifts between the geodetic solution and the EMSC solution (among others) has been observed for several large earthquakes that occurred along the eastern Ionian coast in the recent years, i.e., the Cephalonia doublet 2014 [26], Lefkada 2015 [30] and Zakynthos 2018 [38]. This systematic shift might be primarily due to the velocity model used for the Ionian Sea and Adriatic area by some seismological agencies.
Like in all elastic models, we cannot solve independently all parameters of the modelled fault. In particular, there is a strong trade-off between fault width and amplitude of slip, and between fault azimuth and rake. In our case fixing the rake at 90 • (based on the seismic moment tensors; Table 2) has a direct impact on the azimuth of the modelled fault. The fact that our modelled azimuth is consistent with the azimuth of the geological structures (Figure 1b; Figure 6) supports the reliability of our approach.

GNSS Magnitude of Durres Earthquake
A series of studies [40][41][42][43] produced GNSS peak-ground displacement (PGD) scaling laws and proposed algorithms for their real-time use as an unsaturated and reliable estimator of earthquake magnitude. One of them [42] modelled the magnitude scaling properties of peak ground horizontal displacements (PGD and PGD-S) for strong earthquake events using L1-norm minimisation regression. In this study we use only horizontal offsets at the two GNSS stations (DUR2 and TIR2; Figure 3, Figure S4) to calculate the GNSS magnitude of this earthquake. We used the coseismic offsets of Table 3 and the geodetic centroid location calculated by the inversion model (41.483° N, 19.604° E, 16.5 km depth; Table 4). The PGD formulas are defined as [42]: PGD = (|AN-S| + |AE-W|)/2 (1) PGD-S = (AN-S 2 + AE-W 2 ) 1/2 (2) Therefore, PGD is the mean value of the absolute horizontal displacement in two orthogonal directions (i.e., E-W, N-S; in cm), while PGD-S is the resultant horizontal displacement (in cm).
Plugging the relevant values to the PGD scaling laws (empirical equations (3) and (4) below) Our fault model also indicates that the 8 ± 2 km depth proposed by [17] for the hypocentre is not supported by the geodetic data.
We emphasize the useful role of station DUR2 and, more generally, the importance of the presence of GNSS stations with a sufficient density to give control and absolute scale to InSAR data. In the case of the Durrës earthquake, one or two more permanent stations to the north and east of the geodetic centroid would have provided great added value for the near real-time determination of the fault location. GNSS and InSAR are valuable sources of displacement data. The real challenge now is to use them jointly more and more efficiently and quickly.

Tectonic Strain and Earthquake Recurrence
With GNSS data, we were able to estimate crustal shortening between Durrës and Tirana. This shortening amounts to 3.6 ± 0.5 mm yr −1 of 1-D strain in the azimuth N80 • E. Such a large shortening rate in a small band of crust (about 30 km) and if the decollement at depth (active in that event) represents the main boundary of the blocks (i.e., the basal thrust), then we may be able to estimate earthquake recurrence along the seismic fault of November 26, 2019. Assuming a uniform slip of 0.55 ± 0.1 m that is released periodically (assuming full locking and uniform loading rates), then the recurrence time of M6+ earthquakes would be on the order of 152 ± 10 years (i.e., the amount of co-seismic slip of the event divided by the shortening rate). Further, assuming that shortening strain in the region between Durrës and Tirana is accommodated by two reverse faults (i.e., the 2019 seismic fault and its probable continuation toward south) then the mean recurrence time of a M6+ event would be in the range 70-80 years. We note that a strong M = 6.1 earthquake is reported on December 12, 1926 ([39]; 93 years ago) for the region of Durres with an elliptical isoseismal pattern (long-axis oriented NW-SE).
Using the fault model parameters (Table 4) we drew a cross section (1:1) along the direction ENE-WSW passing through Durrës and the surface projection of the east-dipping seismic fault. We then positioned the inferred seismic fault with its top near the depth of 14 km and its bottom near the depth of 18 km. We also projected the EMSC epicentre to show its position with respect to the inferred fault plane. The geometry of the inferred fault and its kinematics (reverse-slip) match the tectonic style for the region Durrës-Tirana ( [8]; see inset geological profiles in Figure 6). The geological structure is interpreted as a result of Neogene thrusting and folding with west-directed shear indicators [8]. N-S to NW-SE oriented synclines, anticlines and thrust faults have been mapped in central Albania between Durrës and Tirana [7] and published cross sections (Figure 6 inset) depict series of E-dipping late Neogene thrusts at depths down to 10 km. It is therefore reasonable to suggest that the November 26, 2019 earthquake ruptured a basement reverse fault along the Neogene thrust front.

GNSS Magnitude of Durres Earthquake
A series of studies [40][41][42][43] produced GNSS peak-ground displacement (PGD) scaling laws and proposed algorithms for their real-time use as an unsaturated and reliable estimator of earthquake magnitude. One of them [42] modelled the magnitude scaling properties of peak ground horizontal displacements (PGD and PGD-S) for strong earthquake events using L1-norm minimisation regression. In this study we use only horizontal offsets at the two GNSS stations (DUR2 and TIR2; Figure 3, Figure  S4) to calculate the GNSS magnitude of this earthquake. We used the coseismic offsets of Table 3 and the geodetic centroid location calculated by the inversion model (41.483 • N, 19.604 • E, 16.5 km depth; Table 4). The PGD formulas are defined as [42]: Therefore, PGD is the mean value of the absolute horizontal displacement in two orthogonal directions (i.e., E-W, N-S; in cm), while PGD-S is the resultant horizontal displacement (in cm).
Plugging the relevant values to the PGD scaling laws (empirical Equations (3) and (4) below) that were proposed by [42] for the Aegean region, then we obtain a moment magnitude M w (PGD) = 6.44 and M w (PGD-S) = 6.42 assuming a hypocentral distance (R) of 27.957 km to the GNSS station DUR2 (Table 5)  For station DUR2 the magnitude difference (∆M) or M w (seismology)-M w (GNSS) is less than or equal to 0.04 units using the PGD/PGD-S approach [40][41][42][43]. There is a deviation observed in the case of station TIR2, where the estimated magnitude from GNSS is 6.07 (Table 5; PGD relationship) or 6.03 (PGD-S relationship), i.e., it is underestimated by 0.33-0.37 units of magnitude. The reason for this underestimation may be (a) the azimuthal distribution of the two GNSS stations with respect to the direction of seismic slip (assuming a rake of 90 • the slip vector is directed towards N250 • E, i.e., towards the west or the city of Durrës; Figure 5) and/or (b) the uneven distribution of the thick pile of sedimentary rocks in central Albania (including several km of Triassic evaporites) that may differently attenuate seismic energy at both stations.

1.
We identify the main source of the M w = 6.4 earthquake that rocked north-central Albania on November 26, 2019 to be located within the frontal area of the basal thrust of the Dinaric-Hellenic orogen.

2.
We modelled the seismic fault by combining the ascending and descending Sentinel observations. Mixing ascending and descending orbits provides a more robust solution. We find that we can model the overall fringe pattern by reverse slip on an east-dipping fault. The fault plane is a low-angle thrust fault (22 by 13 km) that dips towards the east (23 • ).

3.
The inversion of geodetic data suggests that the upper edge of the fault is at a depth of 14 km, well constrained by the modelling of the interferograms.

4.
Geodetic data GNSS and InSAR ( Figure 2) show ground motion to the southwest and surface uplift in agreement with moment tensor solutions from seismology.

5.
The epicentre published by EMSC is located 15 km southwest of the one deduced from geodesy, this might be due to insufficient inaccuracy of the velocity model of the crust beneath the Adriatic. 6.
It is notable that there is a 123 ns yr −1 active shortening of the crust between the GNSS stations DUR2-TIR2 (equivalent to a shortening rate of 3.6 mm yr −1 ), and roughly in the east-west direction. 7.
Given this amount of strain the recurrence time of M6+ earthquakes along this fault should be of the order of 150 years. 8.
Funding: This research was supported by "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).