Detection of Tectonic and Crustal Deformation using GNSS Data Processing: The Case of PPGnet

Aitolo-Akarnania prefecture, western Greece, is an area with strong earthquakes and large active fault systems. The most prominent are the Katouna sinistral strike slip fault and the Trichonis Lake normal fault system. Their proximity to large cities, and the lack of detailed information on their seismogenic potential, calls for multiparametric research. Since 2013, the area’s crustal deformation has been monitored by a dense GNSS Network (PPGNet), consisting of five stations, equipped with Leica and Septentrio receivers. The objective of this network is to define the rate of deformation across these two main fault systems. Data is recorded using two sampling frequencies, 1 Hz and 10Hz, producing hourly and daily files. Daily data is processed using Bernese GNSS Processing Software using final orbits of International GNSS Service. Double-difference solution is computed using phase measurements from the PPGNet network complemented by four stations from Athens’ National Observatory GNSS network and six stations from METRICA network. First results show a NNE movement at PVOG station of 12 mm/y and a similar movement at RETS station of about 9 mm/y. This means that the Trichonis Lake normal fault system, located between these two stations, depicts a slip rate of 3 mm/y. KTCH and RGNI stations move eastwards at a velocity of about 5 mm/y due to the Katouna-Stamna fault system. Data from PPGNet has provided important results on crustal deformation in the area, i.e. slip rates have been attributed to specific fault systems. The comparison and links of these data with broader geodynamic models is now possible and we expect, in a later phase that will provide a more detailed image of the associated seismic hazard for Aitolo-Akarnania.


Introduction
Greece is an earthquake-prone country located at the convergence boundary between Nubian and Eurasian plate ( Figure 1). The fast convergence rate, ~30mm/y at southern Greece and Aegean Sea, is transformed to crustal collision in northwestern Greece and Adriatic Sea. This transition is accomplished through the Cephalonia dextral strike slip fault located in the Ionian Sea [1,2]. To the east of the Ionian Islands and in the central Greece area, crustal extension in a NS direction dominates. A major tectonic element in this area is the Corinth Gulf rift, which is a fast-extending continental rift (~15 mm/y) [3,4] (Figure 1). 15 The area west of Corinth Gulf rift appears to be a "triple junction" with the E-W normal faults of Corinth Rift to the east, the NW-SE striking Katouna-Stamna fault system (KSF) [5] and the Trichonis Graben to the north at the Aitolo-Akarnania area [6] and the Achaia-Elia (AE) strike-slip fault, of NE-SW direction, to the south [7]. These fault systems are connected to the Kefalonia-Lefkada fault system and the subduction at the southwest, bounding a crustal block that seems internally undeformed.
Recently, Pérouse et al. (2016) [8] proposed the existence of a microplate in the area, i.e. the Ionian Islands-Akarnania Block (IAB), based on GPS and tectonic data. Based on the above it is clear that Aitolo-Akarnania is a key area for understanding the deformation pattern in Western Greece. Chousianitis et al. (2015) [9], using a high-quality dataset of ~100 continuous GPS stations in Greece, identified a shear zone in Aitolo-Akarnania in a NW-SE direction compatible with the KSF fault system. Moreover, they compared the geodetic and seismic moment rates and found a large deficit of the seismic moment rate. This can be attributed either to lack of events from the seismic catalog or to aseismic behavior of the fault zone. These results are compatible with the findings of Pérouse et al. (2012) [10]. Similar behavior was also found in Yuasa et al. (2020) [11].
Although the rough characteristics and the main tectonic elements in the area have been recognized by traditional geological mapping, there is a lack of continuous and dense monitoring of the crustal deformation in the area, i.e. about the activity of main faults, which is an important input in seismic hazard calculation. In this way, by the addition of seismicity information the construction of a detailed fault model will be possible, see e.g. the work by Astupina et al. (2019) [12] for the subduction interface in Peru. Seismic hazard is important for the area, since major cities and significant infrastructures are either close or even cross-cutted by faults, e.g. the highway No. 5 with many tunnels and bridges goes near the Katouna-Stamna fault.
We decided to use regional GNSS networks which, according to Murray et al. (2019) [13], are also suitable for crustal deformation monitoring, in addition to seismicity monitoring. Since 2013 the area is monitored by a network of five GNSS stations (PPGNet) and complementary to permanent seismic network PSLNet [14]. The main purpose of PPGNet is to densify the regional GPS network and monitor the crustal deformation in order to further investigate the existence of the IAB microplate and elucidate the role of the active faults in the area. The main active tectonic elements in Aitolo-Akarnania are the Katouna sinistral strike slip fault and the Trichonis Lake normal fault system. In Figure 2 a simplified tectonic map of the area is presented. Major events, for the time period 1900-2010, with M>5.5, are taken from Papazachos et al. (2010) [15]. Seismicity is of course denser in the Ionian Islands (Lefkada and Cephalonia), where catastrophic events have occurred in the past.
Nevertheless, strong events have also occurred in Aitolo-Akarnania. These events are clearly connected with the Trichonis Lake faults or with the northern part of the Katouna-Stamna fault. Their magnitude is close to 6Mw or slightly larger, thus they pose a significant threat for the major cities in the area ( Figure 2).
What is surprising, and important for the seismic hazard of the region, is the absence of strong events in the central-south part of the Katouna-Stamna fault. This could be attributed to the large return period of events in this part of the fault system or in the creeping behavior of the fault as suggested by Pérouse et al. (2016) [8]. In this article we start by presenting a detailed description of the recently established dense GNSS network and its daily operation. The methodology and the analysis of the data follows and finally the first results of the GNSS network are presented as well as the aims of our research. Stations are equipped by Leica, Septentrio and Trimble instruments. The Septentrio and Trimble receivers produce GPS NAVSTAR data only. Data are stored in RINEX format using two sampling frequencies, 1 Hz for stations RETS, RGNI, VALY and 10Hz for KTCH, LEPE and PVOG. Hourly and daily files are produced, and the daily files with 30 second sampling intervals are freely available. Data with 10Hz sampling are available upon request in the frame of the CzechGeo/EPOS project. Basic network information is listed in Table 1. It describes the state of the PPGNet network up to December 2017.

Processing of the Network
The procedure of data processing and analysis of results is shown in Figure 4  Data quality plots are produced automatically, on a daily basis, using TEQC software [16] https://www.unavco.org/software/data-processing/teqc/teqc.html ( Figure 5). These plots provide a first indication of the data characteristics and assist in data quality management.
Two approaches are used to process phase and code data for determination of crustal deformation: processing data from one station using PPP (precise-point positioning) such as in Avallone et al. (2011) [17]; or processing measurements based on double differences from a pair of stations (double differences) such as in Henrion et al. (2020) [18]. We use the second approach -processing based on double differences. The daily data is processed using the Bernese GNSS Processing Software version 5.0 [19].
For processing automation, the Bernese Processing Engine was used with complementary scripts from EPN (EUREF Permanent Network www.epncb.oma.be) and the local analytic centre Geodetic Observatory Pecny (GOP).
The processing strategy is based on good practice defined during long term processing of EPN subnetworks.  As fiducial sites, we used 9 Class A stations of EUREF Permanent Network (AUT1, COST, DUB2, DUTH, GOPE, MATE, ORID, SRJV, USAL). Data from the GPS NAVSTAR system at elevations above 3 degrees were used because we do not have GLONASS data from all processed stations. Input data in RINEX format with sampling interval of 30s in time period from 2013-11-20 till 2017-10-01 were processed. We obtained data from all stations for more than two years. The IGS05 system with minimum constraint (only translations) to fiducial sites was used for reference system definition, final coordinates are in the IGS05 system in epoch of measurement. Final orbits of the International GNSS Service in the ITRS coordinate system were used for satellite positions. The absolute antenna phase centres model was used at all stations, individual models were used on the stations when available. a) The processing method was the double differenced observation with minimum baseline counts. Quasi-Iono Free (QIF) method was used for ambiguity strategy, which also eliminates the ionospheric refraction error by the Iono-Free L3 linear combination of the L1 and L2 frequencies of observation. The tropospheric refraction was modelled in two steps. The dry part was determined by a priori Niell model, the wet part was computed by the Niell model with estimated parameterszenith troposphere delay one per hour and troposphere gradients one per day. Average repeatability of daily solutions was around 1.5 mm in horizontal coordinate components and up to 6 mm in vertical component. Each coordinate time series (example at Figure 6) was checked for searching of jumps (usually as a result of near earthquake), co-seismic effect (which introduce exponential behavior of time seriesas in Ansari and Bae (2020) [20]) and yearly periodicities (usually as result of temperature deformation of the building where the station is located - Figure 5a). Due to using double differenced processing with good daily repeatability, we don't use the commonmode component filtering technique [21] for improving signal-to-noise ratio of the velocities. The results, i.e. coordinates from the double-differenced processing, are always relative. The velocities (changes of coordinates in time) derived from coordinate time series are primary referenced to the average motion of set of used fiducial sites. The velocities were primary fixed to Patras/Rio GNSS stations PATR/PAT0 (to the average of velocities of both stations) for the first solution shown in Figure 7a. For a better description of fault movements in the area of KFZ fault zone, the second solution (Figure 7b) of velocities was fixed to Paravola GNSS station PVOG (Table 2).

Triangulation Methodology -Strain Regime in the Area
Following the raw data processing we attempt an initial interpretation of the results based on the area's tectonics and on known regional deformation patterns from previously published papers e.g. Chousianitis et al. (2015) [9], Konstantinou et al. (2017) [24]. We have applied the triangulation method, as proposed by UNAVCO (https://www.unavco.org/), to calculate the first order deformation pattern in Aitolo-Akarnania. This methodology is using raw geodetic data on three GPS stations that form a triangle. It takes the velocity at each of the three GPS stations and determines what types of transformations the region between them is undergoing. Finally, it breaks the total measured GPS velocities into components of the different types of transformations -translation, rotation, extension, and strain. Application of the methodology on many triangles results in the deformation pattern of the entire area. For calculating the crustal strain, the GPS Triangle Strain Calculator software developed by UNAVCO was used. Various strain related parameters were computed but we'll focus mainly on maximum/minimum horizontal extension, area strain and rotation which are more tightly connected with the tectonics of the region. Due to the limited extent of the network we will concentrate on the two main fault systems that are the aim of PPGNet.
The maximum horizontal extension (MHE, e1H) is defined as MHE = ( -0)/ 0, where is the final length along the longest axis of the strain ellipse and 0 is the original length, is expressed in nano-strain and it's an indicator of tectonic activity in an area; its direction is perpendicular to the strike of the active fault zones [25]. The minimum horizontal extension (e2H) is defined in a similar way, but along the shortest axis of the strain ellipse. The relative values of the maximum/minimum horizontal extensions as well as their signs, relate to the prevailing tectonic regime i.e. extensional or compressional type. For example, if the maximum horizontal extension is a positive number and the minimum is zero, then we have a purely extensional tectonic regime with the direction of the extension coinciding with the azimuth of the maximum horizontal extension.
The results of the maximum horizontal extension in the area are presented in Table 3 for two triangle layouts, i.e. RETS-PVOG-AGRI and RGNI-KTCH-RETS. For the first triangle layout results indicate pure extension in an almost north-south direction and are compatible with the area's tectonics, i.e. the Trichonis normal fault system is developed in an EW direction perpendicular to maximum extension. Similar results have been proposed for the area by Chousianitis et al. (2015) [9]. The style of deformation is changing in the second triangle layout which is located further to the west, in a continuation of the Katouna-Stamna strike slip fault system to the south. The minimum horizontal extension axis is now a bit larger indicating an area affected by shortening also.
To further investigate the style of faulting in the area, the area strain was also calculated. The area strain reflects the change in area (if any) during distortion and is equal to the sum e1H + e2H measured in nano-strain. A positive value of area strain indicates that the area has increased (extension), and a negative indicates that the area has decreased (shortening). Consequently, the extension is mainly associated with normal faults, or with transtensional strike-slip zones, while the shortening is mainly related with reverse faults, or with transpressional strike-slip zones. The results of the area strain are presented in Table 4 for two triangle layouts i.e. RETS-PVOG-AGRI and RGNI-KTCH-RETS. The first triangle layout results indicate an extensional regime, in accordance with the Trichonis lake fault zone, while the second indicate a reverse or transpressional regime, compatible with the Katouna-Stamna fault system. This is an important result as it supports the continuation of this fault system to the south. This information is important for the area's seismic hazard since the fault trace cannot be mapped directly on the surface. Furthermore, the area strain in the second triangle is a lot smaller, suggesting lower tectonic activity in comparison to Trichonis lake area.
In Figure 7 we attempt to display these results in a graphic form and compare our results with seismological data. The velocities (changes of station's coordinates with time) were calculated with respect to fixed PATR/PAT0 (to the average of velocities of both stations) ( Figure 7a) and with respect to fixed PVOG station (Figure 7b). In Figure 7a we observe a NNE movement of PVOG and AGRI stations at ~12 mm/y and an almost parallel movement of RETS station at a lower velocity ~9 mm/y, this suggests that there is a 3mm/yr slip distributed in the area between them i.e. in the Trichonis lake fault zone. On the other hand, the RGNI and KTCH stations are moving to the NE, with respect to fixed PATR/PAR0 suggesting shortening. In Figure 7b it is also clear that the RGNI and KTCH stations, located at the western part of the Katouna-Stamna fault system, are moving to the south with almost the same velocity, following the fault's kinematic character. The RETS station on the other hand is also moving to the south but with significantly lower velocity, indicating that this station is affected by other factors also, e.g. the Trichonis lake fault zone extension or the Patraikos gulf faults. In Figure 7b the available focal mechanisms of major events in the area are plotted for comparison with the results from GNSS data analysis. It is clear that the focal mechanisms also suggest an extension to the Trichonis lake fault system, roughly in a NS direction, while further west at the Katouna-Stamna fault system the focal mechanisms turn to a strike -slip or oblique thrust type, in accordance with the suggestion of GNSS data analysis of a transpressional seismotectonic regime.

Conclusions
The PPGNet's scope is the study of the tectonic deformation in Aitolo-Akarnania Prefecture, in Western Greece, using GNSS data. Unlike previous field studies and regional observations we had so far, the installation of the GNSS network and the recorded data provides us, for the first time, with detailed and accurate measurements of the area's crustal deformation. The first results are based on the processing of coordinates of GNSS stations from years 2013 -2017. Time series of coordinates were analyzed and the velocities (changes of coordinates in time) and periodicities were determined.
These first results (Figure 7a, b) show a NNE movement of PVOG station at 12 mm/y and a similar movement of RETS station at about 9 mm/y (with respect to the PATR/PAT0 stations in Partas/Rio). This means that the Trichonis Lake normal fault system, that is located between these two stations, depicts a slip rate of 3mm/y. The KTCH and RGNI stations move eastwards at a velocity of about 5 mm/y. Keeping PVOG fixed (Figure 7b) stations RGNI and KTCH depict a SSE movement, while station RETS moves in the same direction but with significant smaller velocity. These motions are in accordance with the Katouna sinistral strike slip fault, while similar motions of RGNI and KTCH further support the existence of a rigid block in the area. Thus, the existence of a microplate border in the area, i.e. the Ionian Islands-Akarnania Block (IAB) proposed by Pérouse et al. (2016) [8], can be confirmed. Although a direct comparison of our results with previous studies in the area is not easy, due to the results analyzed here being new and not used in previous studies, we attempt a comparison of our results with the most recent GNSS based study of deformation in Greece [9]. In this study the authors use ~100 continuous GPS stations in Greece (only two are located in Aitolo-Akarnania) and derive results on the strain deformation. According to Chousianitis et al. (2015) study [9] in Aitolo-Akarnania a shear strain boundary exists, in accordance with the findings in the present study.
Accuracy of determined velocities from time series longer than 2 years will be better than 1 mm/y. Based on results in Henrion et al. (2020) [18] or Ansari and Bae (2020) [20], we can assume that the longer time series will allow more accurate determination of horizontal velocities as well as determination of the vertical component of velocity at individual stations. Therefore, it is believed that data from PPGNet will provide valuable information on the Aitolo-Akarnania area internal deformation and relative motions at the main faults in the area and eventually will help us understand how this deformation is linked to the major active structures in the broader area.

Data Availability Statement
The data presented in this study are available in article.

Funding
Financial support from the CzechGeo/EPOS (LM2015079) project is acknowledged by V. Plicka, V. Filler and J. Kostelecky. E. Sokos and E. Lyros acknowledge financial support by the HELPOS Project, "Hellenic Plate Observing System" (MIS 5002697).

Conflicts of Interest
The authors declare no conflict of interest.