Geodetic Upper Crust Deformation Based on Primary GNSS and INSAR Data in the Strymon Basin, Northern Greece—Correlation with Active Faults

: The Strymon basin (Northern Greece) belongs to the geodynamically active regime of the Aegean and, as expected, it hosts active faults. Nevertheless, the study area exhibits a low instrumentally and historically recorded seismicity. In order to comprehend the crustal deformation, we implemented GNSS- and InSAR-based techniques. Global Navigation Satellite System (GNSS) primary geodetic data recorded by 32 permanent stations over 7 years were analyzed and input in the triangulation methodology so as to calculate a series of deformational parameters. Moreover, a geostatistical methodology indicated the spatial distribution of each parameter, showing strain delimited up to 2750 × 10 − 9 . These results are in broad agreement with palaeoseismological surveys and active fault mapping. Moreover, InSAR analysis, based on a 6-year data recording, concluded that no horizontal rates have been traced in the E–W direction; if they do exist, they would be below resolution (less than 2 mm/yr). Peak vertical subsidence values of a few mm/yr are traced towards the hanging wall of the Serres fault zone within the Quaternary sediments at the eastern margin of Strymon basin but are attributed mainly to groundwater extraction. However, it is noteworthy that geodetic strain analysis implies: (a) that a couple of areas need further study to trace potentially active faults by palaeoseismological means; (b) the fault trace of the Serres fault zone might be further prolonged 8–10 km eastwards, where Quaternary sediments cover the fault.


Introduction
The advance of Global Navigation Satellite System (GNSS) technology is a valuable asset for the study of crustal deformation, providing insights to the strain and velocity fields and, furthermore, to the geodynamic evolution of a region (e.g., [1,2]). Crustal deformation is the result of elastic strain accumulation and, when a threshold is exceeded, almost all of the elastic deformation turns permanent by fault rupturing [3] and, consequently, into earthquake genesis. Mapping and measuring strain rates by combining GNSS and Interferometric SAR provides not only a significant temporal resolution, but, most importantly, a high spatial resolution and coverage, which are essential to fault-specific seismic hazard assessment [4]. The strain rate can be calculated by applying geological (e.g., [5]), seismological (e.g., [6]), and geodetic methods (e.g., [7,8]), or a combination of these [9]. During the last few decades, the latter has mostly been based on the continuous monitoring of a permanent GNSS station network and/or Interferometric Synthetic Aperture Radar (InSAR), providing time series [10] which are well suited to study interseismic strain accumulation (e.g., [11]).
Satellite geodetic analysis is a useful tool for estimating the deformation of the upper crust in a regional scale. Although it covers a relatively short-term period of the seismic cycle (a few years up to a few decades) that maybe inadequate for capturing the full length (which can range from a few hundred up to several thousands of years), it provides a valuable extended snapshot of the present-day deformation field. As a result, it must be emphasized that the estimated short-term geodetic rate may not necessarily reflect the long-term deformation rate, as the latter may not be stable during the seismic cycle [12,13]. Moreover, no conclusion can be drawn about the maturity stage of the fault, and, hence, the time of the next rupture. Nevertheless, satellite geodesy can depict the strain accumulation rate along known faults or along faults which were never mapped before with conventional techniques or provide important strain rate thresholds.
Our study area consists of the Strymon basin in North Greece, which develops in a NW-SE direction from the Greek-Bulgarian border down to the estuary in the Orphanos Gulf (North Aegean Sea). The basin is bounded by active faults controlling the marginal topography, causing subsidence, while low strain rates (30-70 × 10 −9 /yr) are documented by broad-coverage geodetic studies [7,[14][15][16][17]. Low strain rates infer long seismic recurrence intervals; as a result, many faults have not been activated during the historically and instrumentally recorded periods in earthquake catalogues. Moreover, the extracted geological rates are elusive or subtle due to low strain rate per fault, and seismological data are insufficient due to the long recurrence intervals and the incompleteness of the historical earthquake record. Hence, geodetic data are valuable in such a tectonic regime.
The geodetic datasets' processing is based on the triangulation methodology (its basic principle is based on combinations of three different GNSS stations), while the extracted results are further geostatistically processed, implementing interpolation (kriging) methodology, in order to determine the results' spatial distribution in a grid pattern.
Taking all the above into account, the scope of this study is to: (i) estimate the upper crust deformation, based on primary geodetic data, (ii) offer a comparison between GNNS and InSAR results, and (iii) correlate geodetically extracted strain with the surrounding active faults and define the potential seismic hazard.

Geodynamic Setting and Tectonic Pattern
Our study area, located in Central-Eastern Macedonia (mainland Greece), belongs to the broader extensional regime which prevails after the southward migration and roll-back of the Hellenic Arc since the Early Miocene (e.g., [18][19][20][21][22]). The variously proposed (micro-) plate models for the broader Aegean and Anatolian regions (e.g., [23] and references therein) primarily suggest that Central-Eastern Macedonia lies either on an undivided Eurasian plate overthrusting the African plate and undergoing internal deformation, or on a separated Eurasian plate into smaller micro-plates, such as the stable Eurasian to the north and the Aegean and the Anatolian to the south, all overriding the African plate.
The dominant structure north of the Hellenic Arc is the North Anatolian Fault (NAF) that accommodates 24-30 mm/yr of dextral motion [14,24] and its western prolongation in the Aegean Sea, i.e., the North Aegean Trough (NAT) and, further to the west, the North Aegean Basin (NAB). Whereas NAT is a quasi-pure shear fault zone, the NAB is characterized by transtensional deformation and represents a transitional region, which links pure strike-slip faulting along the NAF and pure extension in central Greece (e.g., [25,26]). This lithospheric-scale transcurrent tectonic structure is believed by many researchers to represent the northern margin of the Aegean-Anatolian microplate(s), along which the latter slides westwards, pushing the former to the SW and affecting the surrounding area, including Macedonia and Thrace (North Greece), by forming a transtensional field. The same stress field is responsible for the formation of large, tectonically controlled basins usually delimited by antithetic normal, dip-slip, tectonically active faults [18,20,[27][28][29][30][31][32][33][34][35][36][37].
Although this region is characterized by low seismicity, compared to other regions of the Aegean, a strong earthquake has been recorded. In particular, this seismic event occurred in Amphipolis (597 AD, ~M6.8), while no seismic events exceeding Μ5.0 have been documented for the last 120 years [47,48].
Low seismic activity in Central-Eastern Macedonia is in agreement with the low values of the geodetic strain rate. As seen in Figure 1 [16], the strain rate in this area is contrastingly lower than the neighbouring NAT or NAB. The same pattern has been suggested by many other researchers as well [7,8,16,17,[49][50][51]. The strain tensor analysis from previous works may show some variances in the main axes' direction, but they all agree with the relatively much lower values compared to the southern parts of the Aegean [7,8,[14][15][16][17]24,[51][52][53][54][55][56][57]. In particular, Mouslopoulou et al. [58], who focused on the Strymon basin fault system, suggested significantly higher strain values compared to the previous works, while no strain axes direction was proposed; however, even these values are lower compared to South Aegean region.  [20,25,26,[59][60][61][62][63][64][65][66][67][68]). The North Aegean Basin (NAB) and the North Aegean Trough (NAT) comprise the northern strand of the North Anatolian Fault (NAF), while the Central Aegean Trough represents the corresponding southern strand. The Corinth Gulf (CG) is the most extensionally growing region in the Aegean. In the Ionian Sea, the Apulian-Aegean Collision Zone (AACZ) is shown to the north and the Hellenic Arc to the south, separated and displaced by the right-lateral Cephalonia Transfer Fault Zone (CTFZ) in the middle. Hatched areas represent extension. Moment tensor solutions since 1997 from the RCMT catalogue [69] are also shown. The white rectangle indicates the location of the inset map (c). (c) Hillshade map showing the active/Neotectonic faults of the study area (F = Fault, FZ = Fault Zone). Recent earthquake spatial distribution (M L ≥ 0.5) since 2008 by the Institute of Geodynamics, National Observatory of Athens (IG-NOA), is also shown.
Tectonic strain is partitioned among the active faults of the study area; therefore, a holistic approach studying strain should incorporate the active structures. Overall, major fault zones in the study area are striking both E-W (Sochos, Drama, etc.) and NW-SE (e.g., Symvoli-Fotolivos, Tholos-Nea Zichni), which is in agreement with the palaeoseismic results [45]. The most important known fault zones in our study area are: 1.
The Symvoli-Fotolivos fault zone The Symvoli-Fotolivos fault zone (SFFZ) is a NW-SE-striking normal fault zone, dipping to the NE, forming the Drama-Phillipi basin western margin. It is approx. 53 km long and comprises several fault segments. Palaeoseismic trenching, supported with radiometric dating, revealed that: (i) the fault plane dips at 55 • towards NE (052 • ), in agreement with the strike of the fault zone extracted from the surficial geological and geomorphological mapping, and (ii) it is Holocene active, but with low slip-rate (0.14 mm/yr) and, hence, long earthquake recurrence intervals in the order of several thousands of years [45]. No major historically or instrumentally recorded earthquakes can be correlated with this fault, most probably due to its long recurrence interval.

2.
The Tholos-Nea Zichni fault zone The Tholos-Nea Zichni (TNFZ) is a NW-SE-trending major normal fault zone that dips and downthrows to the northwest (225 • ), forming the eastern boundary of the Strymon basin [45]. The area bounded by the two fault zones (TNFZ and SFFZ) comprises of Alpine Mesozoic metasediments, which are unconformably overlain by Neogene fluvial and lacustrine deposits.
The TNFZ displays common tectonic geomorphology features such as a relatively linear topographic boundary, inverted wine-glass valleys, and highly incised catchments that result from footwall uplift (e.g., Figure 2). Moreover, the direction differentiation of the ongoing extension, as well as the fault strike, may indicate that a limited horizontal component may occur. According to recent geodetic studies [58], a sinistral lateral component of motion is proposed. No significant events have been recorded instrumentally near the fault (e.g., all events were below M = 5.0). Only one strong historical event (M~6.8), which occurred in AD 597, has been recorded a couple of tens of km from the fault into its hanging wall, where the maximum subsidence is expected. This event changed the flow of the Strymon River near Amphipoli [70]. Therefore, this fault is one of the few candidates (possibly the most proximal one) that could have hosted this event. Palaeoseismic trenching, supported with radiometric dating, revealed that this is a Holocene active complex fault zone, with its deformation distributed across several synthetic and antithetic planes, and is interpreted as having a low slip-rate fault (0.20 ± 0.1 mm/yr), with recurrence intervals in the order of 2500 ± 1000 yrs [45]. Overall, the fault trace separates uplift from subsidence, controlling the present-day sedimentation processes, and also demonstrating that the fault is active.
According to the initial morphotectonic observations [75,76], the Mt Belles-Kerkini was considered as an E-W-trending horst (Pliocene-Quaternary), delimited by two adjacent grabens, expanding northerly and southerly, respectively. The Quaternary fault activity affected the Pliocene alluvial deposits, lying along the alpine bedrock-composed foothills. In particular, geomorphological features such as well-developed fault-scarp zones, the deposition of new generations of fans upon the older ones close to the fault zone, post-depositional gully erosion, etc., resulted in the conclusion that the western part of Mt Kerkini is more active than the eastern part.
Later on, Mountrakis et al. [41] separated the fault zone into two fault segments, i.e., an 18 km-long segment to the west and a 24 km-long segment to the east. The former, consisting of two subparallel strands, dips to the south in a steep angle, and, westerly, it merges with an ENE-WSW-trending fault, extending up to Doirani Lake. The latter separates the sediments of Strymon basin from the alpine metamorphic rocks (Mt Belles-Kerkini). The authors also observed that the Petritsi fault intersects the NE-SW (to ENE-WSW)-trending faults north of Petritsi village.

The Sochos fault zone
The Sochos fault zone forms a smaller intermontane basin on the Mt Vertiskos crystalline mass [77,78]. Mountrakis et al. [41,79] reassessed its location and length and characterized it as an active structure. It consists of at least two segments, named the Sochos and Mavrouda faults, in a right-stepping geometry [79,80]. Mountrakis et al. [41] assessed
According to the initial morphotectonic observations [75,76], the Mt Belles-Kerkini was considered as an E-W-trending horst (Pliocene-Quaternary), delimited by two adjacent grabens, expanding northerly and southerly, respectively. The Quaternary fault activity affected the Pliocene alluvial deposits, lying along the alpine bedrock-composed foothills. In particular, geomorphological features such as well-developed fault-scarp zones, the deposition of new generations of fans upon the older ones close to the fault zone, postdepositional gully erosion, etc., resulted in the conclusion that the western part of Mt Kerkini is more active than the eastern part.
Later on, Mountrakis et al. [41] separated the fault zone into two fault segments, i.e., an 18 km-long segment to the west and a 24 km-long segment to the east. The former, consisting of two subparallel strands, dips to the south in a steep angle, and, westerly, it merges with an ENE-WSW-trending fault, extending up to Doirani Lake. The latter separates the sediments of Strymon basin from the alpine metamorphic rocks (Mt Belles-Kerkini). The authors also observed that the Petritsi fault intersects the NE-SW (to ENE-WSW)-trending faults north of Petritsi village.

4.
The Sochos fault zone The Sochos fault zone forms a smaller intermontane basin on the Mt Vertiskos crystalline mass [77,78]. Mountrakis et al. [41,79] reassessed its location and length and characterized it as an active structure. It consists of at least two segments, named the Sochos and Mavrouda faults, in a right-stepping geometry [79,80]. Mountrakis et al. [41] assessed that this fault zone was activated on 29 September 1932 (M6.2, according to Papazachos and Papazachou [47]), and its isoseismals were parallel to the fault zone [81].

5.
The Serres fault zone The Serres fault zone (SFZ) lies along the margin of the corresponding basin, i.e., the Strymon basin, while it consists of discontinuous fault traces of a general E-W direction downthrowing to the south.
The SFS is located along the central-east margin of the Strymon basin, disrupting its general NW-SE trend. The eastern part of the fault zone shows a slight change of direction to the ESE, toward the Tholos-Nea Zichni fault zone. This part is shown on the respective 1:50,000-scale geological map of the Institute of Geology and Mineral Exploration (IGME) [82] as affecting the Quaternary sediments that lie along the foothills of Mt Menikion. The western, E-W-striking part of the fault zone, passing by and through the city of Serres, was later mapped by Tranos [83], Tranos and Mountrakis [84], and Mountrakis et al. [41], who characterized it as active. In fact, the part which passes through the city of Serres is the westmost part of the zone showing an ENE-WSW orientation. Normal dip-slip slicken lines overprinted on older strike-slip ones were found on the main boundary faults separating the crystalline basement of Mt Menikion from the post-alpine sediments of the basin [84].

6.
The Drama fault zone The Drama fault zone (DFZ) shows common properties with SFZ, as they both trend towards E-W, downthrowing to the south, and are comprised of discontinuous fault traces.
However, a difference between them is that the DFZ bounds entirely the northern part of the Drama basin against the Mt Phalakron. It is characterized as an active normal fault since it bounds Quaternary deposits and the travertine sedimentation and displays dip slip striations on the fault plane [41]. The DFZ might have been activated during the 5 May 1829 catastrophic earthquake in Drama (M7.3, X intensity), whose isoseismals were elongated approx. E-W to WSW-ENE, and thus parallel to the fault zone [81]. A more recent event M = 5.5 that occurred on 9 November 1985 in Drama (VII intensity) confirmed again the activity of this zone and exhibited a similar fault parallel-elongated isoseismal pattern [81]. The main active faults and fault zones have been described above and concern either E-W-or NW-SE-trending faults.

7.
The Nigrita fault zone It forms the western boundary of the Strymon basin and is traced towards the southern part of the basin. It is divided into (i) a NW-SE-trending segment that is aligned to the older Pliocene basin and (ii) a younger WNW-ESE-striking fault. Both structures merge towards the town of Nigrita.
The younger structure is a 20 km-long, WNW-ESE-trending normal fault that dips and downthrows to the NNE, shaping the southwestern boundary of the Strymon basin. It has a weak topographic expression within the Quaternary sediments, implying that it is active but with low slip-rate and long recurrence intervals. It could be an antithetic fault to the Serres fault since the latter displays a more profound imprint on the topography.
Finally, according to Tranos [73], both strike-slip and normal, and sometimes oblique, dip-slip faulting occur in the broader Strymon region, with the former preceding the latter, respectively. In particular, strike-slip faulting is characterized by NW-SE and NE-SW to ENE-WSW striking, while the normal and oblique dip-slip faults exhibit E-W striking, which corresponds to the N-S extension. In conclusion, it should be noted that secondary faults also occur ( Figure 3) with poor geological evidence, though with distinctive, although subtle, geomorphological expression.  [85], as well as the active-neotectonic faults as shown in Figure 1; (b) 16 permanent GNSS stations (black triangles), monitoring the broader Strymon region.

Seismic History of the Study Area
No strong events (Μ > 5) have been instrumentally recorded in the area. However, there is microseismic activity with several events up to M5.0; since 2008, when the Institute of Geodynamics, National Observatory of Athens (IG-NOA)'s catalogue became more accurate and complete, seismic activity has been recorded mostly at depths down to 40 km [48], i.e., within the upper crust. The hypocentral spatial distribution is occasionally concentrated in specific areas, in clusters, highlighting the occurrence of faults. As seen in Figure 1, most of the seismicity is correlated with geomorphological depressions, such as the Strymon and Drama-Philippi basins, and particularly close to or along their margins. Mountainous areas, such as Mt Vertiskos, show gaps in spatial distribution. The maximum earthquake magnitude recorded over the past 13 years did not exceed M4.9.
Regarding the historical seismicity, only one major and destructive earthquake is documented, the ~M6.8, which occurred in AD 597 towards the southern part of the basin near Amphipolis. At this point, the earthquake changed the flow direction of the Strymon River [47]. However, there are no isoseismals from this very old event to help us constrain the locality and trend of the activated fault.

Seismic History of the Study Area
No strong events (M > 5) have been instrumentally recorded in the area. However, there is microseismic activity with several events up to M5.0; since 2008, when the Institute of Geodynamics, National Observatory of Athens (IG-NOA)'s catalogue became more accurate and complete, seismic activity has been recorded mostly at depths down to 40 km [48], i.e., within the upper crust. The hypocentral spatial distribution is occasionally concentrated in specific areas, in clusters, highlighting the occurrence of faults. As seen in Figure 1, most of the seismicity is correlated with geomorphological depressions, such as the Strymon and Drama-Philippi basins, and particularly close to or along their margins. Mountainous areas, such as Mt Vertiskos, show gaps in spatial distribution. The maximum earthquake magnitude recorded over the past 13 years did not exceed M4.9.
Regarding the historical seismicity, only one major and destructive earthquake is documented, the~M6.8, which occurred in AD 597 towards the southern part of the basin near Amphipolis. At this point, the earthquake changed the flow direction of the Strymon River [47]. However, there are no isoseismals from this very old event to help us constrain the locality and trend of the activated fault.

Primary Geodetic Raw Data Analysis
The extraction of geodetic raw data, collected by permanently installed GNSS stations, is a cutting-edge, instrumental approach for estimating the upper crust strain rate. These primary geodetic raw datasets are partially related to the general, geodetic velocity pattern, calculated throughout the Greek region [93]. Particularly, in the wider Strymon basin, 16 permanent GNSS stations have been installed, recording the aforementioned geodetic datasets; the denser GNSS network within this area is characterized by a uniform distribution. It is mentioned that these stations are distributed within five different permanent networks (HxGN/SmartNet Greece, Hepos, NoaNet, EUREF Permanent Network-EPN and HermesNet; [94,95]), providing additional value and reliability to the recorded data.
Regarding the geodetic velocity calculation, a 7-year (2008-2014) interval was analyzed, while the data recording was performed in a 30 s rate during 24-h operation. In addition, the proposed geodetic velocities reference frame is the European Terrestrial Reference Frame 2000 (ETRF2000), according to which the Eurasian plate is considered fixed. The GAMIT/GLOBK software suite (release 10.5) was implemented [96,97]. Firstly, the raw data were processed for the entire time period in order to estimate daily loosely constrained solutions of sites, while the orbital and Earth Orientation parameters remained fixed to the IGS final and IERS Bulletin A values with weights, respectively. In our strategy, we estimated Zenith Total Delay (ZTD), adjusted every 2 h interval, and the Vienna Mapping Function1 (VMF1) was selected as a mapping function [96]. Concerning the Solid Earth Tides, the IERS2003 conventions were followed, in addition to Ocean Loading, while we applied FES2004 model, as it was considered appropriate to apply Atmospheric Pressure Loading corrections. In the second step, we combined all the individual sessions of the local subnetworks, using a seven-parameter Helmert transformation, with WRMS 2.0, 2.2, 4.8 mm on the east, north, and vertical components, respectively, estimating site coordinates by using the individual daily solution and geodetic velocities produced by the time series analysis with the linear and extended model. The differences between the two models are <1 mm/yr for the data span, exceeding the 2.5 years, while the extended model is described by the following equation: where α and β are the position and velocity of monitoring station, c k and d k describe periodic motions, and the equation is used to describe the discontinuities caused by the earthquakes and equipment changes, with magnitudes g k at epochs T gk and H is the Heaviside step function. On next step, after the outliers' removal, we estimated station coordinates and velocities using the Kalman filtering sequential approach. The realization of the frame was estimated by minimizing the deviations between positions and velocities, relative to the current terrestrial frame with 7-parameter Helmert transformation and Post RMS 2 mm and 0.6 mm/yr for the coordinates and velocities, respectively. Table 1 shows the station velocities of the processed permanent GNSS sites relative to ETRF2000.

Implementation of Triangulation Methodology-Geostatistical Analysis (Interpolation) of Triangulation Methodology Results
The primary geodetic data processing, analyzed above, resulted in the estimation of the east and north velocity component of each GNSS station, as well as in the corresponding errors of these components. The triangulation methodology, implemented widely by various authors [98][99][100][101][102][103], was considered as the most suitable between different methodologies (least square methodology, grid-based methodologies, etc.) in order to process these datasets. The principle of the triangulation methodology is the velocity components' combination of three GNSS stations, and their errors, leading to the total velocity calculation. In particular, these stations form a triangle, while each of them is considered as a triangle vertex. Then, the triangle is examined before and after the deformation process [57,104,105]. The undeformed triangle is represented by an inscribed inner circle, while the corresponding one is deformed by an ellipse, respectively [104,105]. The inner circle and ellipse axes' combination lead to the strain parameters' extraction. It should be noted that the procedure mentioned above is based on mechanics principals (more details about the triangulation methodology are provided in the Supplementary Materials).
As has already been mentioned, the wider Strymon basin area is surrounded by 16 permanent GNSS stations-a fact that leads to a large number of potential stations combinations. However, the total number of triangles that were eventually constructed and considered for the deformation estimation meet the following criteria: (i) the vast majority of triangles are equilateral in order to be affected equally by the deformation pattern of the three combined GNSS stations and (ii) the smallest possible distance between the GNSS stations is considered in order to efficiently observe the crustal deformation related to specific tectonic structures while significant interaction between the combined stations is documented at the same time [103]. Therefore, 216 triangles were constructed, represented by the corresponding triangle centroids. For all these centroids, a series of strain parameters was estimated, including ('GPS triangular calculator' software; Beyond the point spatial distribution of these parameters, a geostatistical process was performed in order to determine the grid pattern deformation. The kriging interpolation methodology (linear approach) was considered to be the most reliable geostatistical methodology, as it smoothly distributes the parameter values within the study area while emphasizing the extreme value concentrations as well as the neutral-transition areas, which are vital for the tectonic interpretation, as they are associated with tectonic regime changes.

Maximum (MaHE) and Minimum (MiHE) Horizontal Extension
The parameters of Maximum (MaHE) and Minimum (MiHE) Horizontal Extension are directly related to the identification of active tectonic structures. Particularly, the MaHE parameter indicates the extensional tectonics occurrences (extension in strain terminology), while, on the other hand, the MiHE, which is perpendicular to the MaHE, represents compressional tectonics occurrences (compression in strain terminology). The MaHE and MiHE parameters are interconnected, and, therefore, the three following cases are observed: It is noted that the MaHE develops along the major axis of the strain ellipse, as was explained above, while the MiHE develops along the minor axis of the strain ellipse, respectively. The following equations (Equations (3) and (4)) define mathematically the MaHE and MiHE parameters: where l f : the final length along the strain ellipse major axis, l 0 : the original length along the strain ellipse major axis, l f : the final length along the strain ellipse minor axis, and l 0 : the original length along the strain ellipse minor axis. It is mentioned that l 0 = l 0 in case of circle occurrences.

Velocity (V)
Velocity (V) is associated with the motion of an area, which indicates partially the geotectonic evolution of a larger scale region. Regarding the Aegean region, the highest values are documented in the region (and adjacent regions) of the Eurasia-Nubia subduction zone [7,14,16,107]; therefore, low V values are expected within the study area, as it is not located close to the subduction zone. The V is calculated by implementing the east and north velocity components, which are perpendicular to each other.

Maximum Shear Strain (MaxSS)
The Maximum Shear Strain (MaxSS) is a parameter indicating the extensional tectonic activity, occurred by shearing activity, developing along active faults or fault zones. These shearing effects result in the brittle, crustal deformation [108]. It is noted that shearing, and therefore the MaxSS parameter, is related to two kinematic types: (i) normal dip-slip and (ii) strike-slip. On the contrary, reverse faults or fault zones lack shear, as they are mainly associated with compressional tectonic activity. The following equation (Equation (5)) represents the mathematical expression of MaxSS:

Area Strain (AS)
The Area Strain (AS) is a parameter which decisively contributes to the determination of the crustal deformation type; two major types are identified: (i) dilatation and (ii) compaction. Dilatation receives the positive AS values, while it is directly related to extensional active tectonics; therefore, dilatation is identified within areas, where normal dip-slip or transtensional strike-slip faulting occurs. On the other hand, compaction is represented by the negative AS values, while it is associated with compressional active tectonics; hence, compaction occurrences are observed within areas of reverse dip-slip or transpressional strike-slip faulting. The mathematical equation of the AS parameter is the following (Equation (6)):

SAR Data and Interferometric Processing
Imaging geodetic techniques such as Syntheric Aperture Radar (SAR) Interferometry have been proven to be an irreplaceable tool for measuring surface motion in different contexts and environments, including land subsidence [109], slope instabilities [110,111], volcanoes [112], and earthquakes [113,114]. Interferometric SAR (InSAR) has evolved over the past years into an operational technique for monitoring surface motion, allowing wide area coverage (including remote areas) at relatively low cost. Towards this direction, the systematic availability of open and free SAR data from the Copernicus Sentinel-1 mission [115] played an important role, ensuring the applicability of InSAR techniques on a global scale without the requisite for planning or ordering SAR acquisitions for specific sites. However, the requirement for big data processing both in terms of computational power and storage capacity led to the devolvement of cloud-based solutions [116]. Among such solutions, the Geohazards Exploitation Platform (GEP) is an initiative setup by the European Space Agency (ESA), aiming to support further exploitation of satellite Earth Observation (EO) data for geohazards assessment [117]. GEP is a cloud-based environment that provides a set of processing tools and services that allow for the mapping of hazardprone land surfaces and the monitoring of terrain motion. It aims to support scientists and EO practitioners to better understand geohazards and their impact while introducing innovative concepts, such as the exploitation of the e-collaboration environment. The service has been already applied in scientific and operational frameworks [118,119].
For mapping the ground displacements of the Strymon basin, we utilized the Surface motioN mAPPING (SNAPPING) service on the GEP. SNAPPING enables the Persistent Scatterers Interferometric (PSI) analysis of Sentinel-1 SLC data for measuring average displacements as well as the investigation of their temporal evolution (e.g., displacement time series). The service is implemented to facilitate on-demand processing and the easier monitoring options of a user-defined area of interest.
The SNAPPING medium resolution service is considered sufficient for wide area processing, followed by more detailed investigation when a relevant signal (in our case, fault-related motion) is detected. The wide overview of ground displacements at medium resolution, reducing relevant production costs, is a proper practice that it well demonstrated in relevant research and development activities [120]. Medium resolution refers to the merging of Persistent Scatterers (PS) point measurements within a 100 m radius. However, PS targets exhibiting different phase properties are not affected, thus maintaining local displacement patterns.
In practice, InSAR time series measurements reach an accuracy of 1-2 mm/yr for average displacement rates and 3-5 mm for each displacement estimate in the time series [121]. This should be critically assessed when a region of relatively low interseismic motion is being investigated. However, the abovementioned accuracy level might be affected by the time span of the observation period and the amount of SAR imagery used in the processing. To ensure the robustness of our analysis, we expanded the InSAR observation period from April 2015 to December 2020. A period of approx. 6 years is acceptable for obtaining a millimeter level accuracy. In turns, Sentinel-1A and -1B satellite units were exploited, repeated every 6 days, to better address non-linear temporal signals while minimizing error.
A total number of 164 acquisition dates from ascending track 102 (10 April 2015-27 December 2020) and 159 dates from descending track 7 were considered (4 April 2015-9 December 2020) following a fully unsupervised processing on the GEP platform. An area of about 10,200 km 2 (120 × 85 km) was examined. To restrict the analysis over high-quality measurements points, interferometric processing was performed, setting the amplitude dispersion for PS detection to a value of 0.35. We avoided the definition of a local reference area for the PSI processing. This is a preferable solution to minimize potential influence by unknown deformation signal, especially when no relevant information on the stability of the areas is available a priori. Thus, each solution was referenced to the average motion of the entire area (mean velocity is set to zero). The above decision does not affect the relative motion between neighboring points or the observed temporal patterns in the time series; rather, it only affects the absolute velocity of the point measurements. SNAPPING PSI measurements are provided in the satellite Line-of-Sight (LoS) direction. Positive values correspond to uplift or motion towards the satellite, whereas negative ones correspond to subsidence or motion away from the satellite platform. A detailed description of the SNAPPING service and the processing steps involved are described elsewhere [122].
Finally, a decomposition scheme was applied [123] to retrieve vertical and horizontal E-W motion components combining ascending and descending SNAPPING measurements. Initial steps involved the rasterization of point measurements into a common 100 m grid, assigning the median displacement rate value when multiple points are identified within the raster resolution cell.

Maximum (MaHE) and Minimum (MiHE) Horizontal Extension
The results of Maximum and Minimum Horizontal Extension parameters (MaHE and MiHE, respectively) are characterized by a wide range of strain field values. In particular, the MaHE values vary between −55 and 2649 × 10 −9 , while the corresponding MiHE values vary between −631 and 146 × 10 −9 , respectively (Figure 4b,c). In addition, the e xx , e xy , e yy (maximum, minimum, and average) values, and the corresponding uncertainties (Table 2)     The MaHE and MiHE results were also geostatistically processed, implementing interpolation (kriging) methodology in order to facilitate visualization through spatial distribution. Initially, the MaHE geostatistical analysis (Figure 4b) shows three major high-value locations (1700-2649 × 10 −9 ) located at the eastern margin (two sites) and the northwestern margin (third site). All three locations correlate well with the occurrence of faults and fault zones (e.g., Serres and Tholos-Nea Zichni fault zones), while the mediumand low-value locations (800-1700 × 10 −9 ) indicate fewer deforming faults or fault zones (e.g., Belles-Kerkini fault zone). Regarding the respective MiHE analysis (Figure 4c), the absolute high-value locations (−631-−350 × 10 −9 ) almost correspond to the MaHE ones; however, a significant difference between the magnitude of the MiHE and MaHE values is documented. Moreover, the lack of the high values of MiHE at the northwestern part of the study area confirms the pure, dip-slip, normal faulting regime of this part, while, on the other hand, the significant MiHE high values along the eastern Strymon basin boundary indicates increased shear, related to the oblique-slip normal faulting. Moreover, the similar MaHE and MiHE values at the southern part of the study area, which are approaching the North Aegean Sea, indicate the potential occurrence of strike-slip tectonic structures, probably related to the major tectonically active strike-slip structures located within the North Aegean Trough (NAT).

Velocity (V)
The calculated Velocity (V) values range between 1.

Maximum Shear Strain (MaxSS)
The Maximum Shear Strain (MaxSS) values within the broader Strymon area range between 0.7 and 2743 × 10 −9 (Figure 6b). In particular, the maximum MaxSS values are documented along the eastern boundary of the Strymon basin, while similar values are also observed at the northern part of the corresponding western boundary (Figure 6a). In addition, much lower but considerable MaxSS values also occur at the southern part of the western boundary Strymon basin. The geostatistical analysis (Figure 6b) shows a clearer view of the MaxSS distribution, indicating the locations of the higher values along a WNW-ESE-trending imaginary axis crossing the Strymon basin from the southern tip of Lake Kerkini to the west to Nea Zichni to the east. It is noteworthy that extremely low, almost negligible values occur in all of the Drama basin. Therefore, the high shear values indicate the significant active tectonic structures along the Strymon basin boundaries (especially the eastern one).

Area Strain (AS)
The Area Strain (AS) parameter is characterized by two different types: (i) dilatation, associated with extensional tectonics, and ii) compaction, related to compressional tectonics. The AS results indicate a clear dilatation dominance over compaction (Figure 7a), with values ranging between −596 × 10 −9 (compaction) and 2715 × 10 −9 (dilatation), suggesting an almost 4.5 times higher dilatation than compaction (Figure 7b,c).
Although the dilatation type prevails within the study area (Figure 7b), it unevenly coexists in several cases with compaction, implying the involvement of shear. In particular, the maximum values' concentrations are documented strictly within the Strymon basin region, along the eastern boundary, as well as south of Kerkini Lake (northwestern part of the Strymon basin boundary). Moreover, it is worth mentioning the significant dilatation values observed at the central and eastern part of Mt Belles at the northwestern part of the study area. Regarding the high-value compaction concentrations (Figure 7c), they are mainly observed along the eastern margin of Strymon basin, southwesterly of the western Strymon basin margin within the broader Mt Vertiskos area, near the coastal part

Area Strain (AS)
The Area Strain (AS) parameter is characterized by two different types: (i) dilatation, associated with extensional tectonics, and ii) compaction, related to compressional tectonics. The AS results indicate a clear dilatation dominance over compaction (Figure 7a), with values ranging between −596 × 10 −9 (compaction) and 2715 × 10 −9 (dilatation), suggesting an almost 4.5 times higher dilatation than compaction (Figure 7b,c). A significant difference between dilatation and compaction values indicates extensional (dilatation > compaction) or compressional (compaction > dilatation) active tectonics occurrences, while similar dilatation and compaction values show shear tectonics activity [57]. For instance, in the broader Kerkini Lake and Mt Belles area, the compaction values are absent or extremely low, indicating that this part of the study area is dominated by extensional active tectonics, manifested by dip-slip, normal faults (or fault zones). Regarding the sensu stricto area of Strymon basin, a clear dilatation (qualitative and quantitative) prevalence over compaction is documented; however, the compaction occurrence Although the dilatation type prevails within the study area (Figure 7b), it unevenly coexists in several cases with compaction, implying the involvement of shear. In particular, the maximum values' concentrations are documented strictly within the Strymon basin region, along the eastern boundary, as well as south of Kerkini Lake (northwestern part of the Strymon basin boundary). Moreover, it is worth mentioning the significant dilatation values observed at the central and eastern part of Mt Belles at the northwestern part of the study area. Regarding the high-value compaction concentrations (Figure 7c), they are mainly observed along the eastern margin of Strymon basin, southwesterly of the western Strymon basin margin within the broader Mt Vertiskos area, near the coastal part of the Orfanos Gulf region and in the north part of the Drama basin (close to Kato Nevrokopi village).
A significant difference between dilatation and compaction values indicates extensional (dilatation > compaction) or compressional (compaction > dilatation) active tectonics occurrences, while similar dilatation and compaction values show shear tectonics activity [57]. For instance, in the broader Kerkini Lake and Mt Belles area, the compaction values are absent or extremely low, indicating that this part of the study area is dominated by extensional active tectonics, manifested by dip-slip, normal faults (or fault zones). Regarding the sensu stricto area of Strymon basin, a clear dilatation (qualitative and quantitative) prevalence over compaction is documented; however, the compaction occurrence is noteworthy, showing that, despite the dilatation dominance, a significant shear component is observed, leading to active oblique-slip, normal faults' (or fault zones') occurrence. Finally, the similar dilatation and compaction values at the southern part of the study area (Amphipolis-Paleokomi-Galipsos-Ofrynion area), where the Strymon River discharges into the North Aegean Sea, indicates potentially active strike-slip tectonic structures, extending within the offshore North Aegean region, which are associated with the tectonically and seismologically active strike-slip North Aegean Trough (NAT) fault zones.

EO-Based Terrain Motion Measurements
Interferometric PSI processing revealed multi-temporal ground displacements over the entire Strymon basin, including the surrounding mountain ranges (Figures 8 and 9). Following SNAPPING PSI processing, approx. 50 k and 70 k point targets were detected for the ascending and descending tracks, respectively. This can be translated to an average density of 6-8 points per square kilometer, which is acceptable given the processing parameters used (low amplitude dispersion), the significant vegetation cover (i.e., agricultural lands and forested regions), and the medium-resolution PSI service used. PS points are located, as expected for the PSI technique, mainly over buildup areas and non-vegetated natural terrain.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 19 of 29 Figure 8. Average LoS displacement rates for the period 2015-2020, derived by SNAPPING service on GEP using Sentinel-1 data from ascending track 102. White lines (a-e and f-k) correspond to spatial profiles shown in Figure 8. Please note that, for the selected color scale, rates higher than ±4 mm/yr appear saturated to allow for the proper visualization of spatial variability of motion. Figure 8. Average LoS displacement rates for the period 2015-2020, derived by SNAPPING service on GEP using Sentinel-1 data from ascending track 102. White lines (a-e and f-k) correspond to spatial profiles shown in Figure 8. Please note that, for the selected color scale, rates higher than ±4 mm/yr appear saturated to allow for the proper visualization of spatial variability of motion.
The overall agreement between the independent geodetic techniques (i.e., GNSS and InSAR) for the calculated annual motion rates is clearly shown in Figure S1. Correspondingly, observations from different SAR geometries show consistency in terms of recognized patterns and motion sign, which serves as an inter-verification of results by independent EO datasets. For the majority of point measurements (approx. 95%), displacement rates do not exceed ±2 mm/yr, and are within the uncertainties of the PSI measurements (an order of a couple of millimeters per year; see Supplementary Figures S2 and S3). This signifies for the overall stability of the area or for the presence of rather low motion gradients. Figure 8. Average LoS displacement rates for the period 2015-2020, derived by SNAPPING service on GEP using Sentinel-1 data from ascending track 102. White lines (a-e and f-k) correspond to spatial profiles shown in Figure 8. Please note that, for the selected color scale, rates higher than ±4 mm/yr appear saturated to allow for the proper visualization of spatial variability of motion. Figure 9. Average LoS displacement rates for the period 2015-2020, derived by SNAPPING service on GEP using Sentinel-1 data from descending track 7. White lines (a-e and f-k) correspond to spatial profiles shown in Figure 8. Please note that, for the selected color scale, rates higher than ±4 mm/yr appear saturated to allow for the proper visualization of spatial variability of motion. Figure 9. Average LoS displacement rates for the period 2015-2020, derived by SNAPPING service on GEP using Sentinel-1 data from descending track 7. White lines (a-e and f-k) correspond to spatial profiles shown in Figure 8. Please note that, for the selected color scale, rates higher than ±4 mm/yr appear saturated to allow for the proper visualization of spatial variability of motion.
Regions of relatively higher displacement gradients, reaching −16 mm/yr, correspond to local patterns (i.e., subsidence) that are mostly recognized on flat areas within Strymon basin (Figures 8 and 9). The circular to elliptical shape of these local patterns most probably relates to groundwater extraction phenomena. The temporal variability of ground displacements for these areas, as derived by the analysis of the corresponding time series (Figure 10), also highlights the relationship of these temporal deformations to water table variations in underground aquifers [124,125]. Positive motion rates, negligible in terms of magnitude, and corresponding to uplift at the highest parts of the mountains, can be attributed to a residual topo-dependent atmospheric signal that has not been compensated for. Focusing on more tectonic-related signals, the lack of linear patterns along the fault zones at the boundary of the basin does not allow for the correlation of any InSAR signal to the interseismic long-term slip-rates of these tectonic structures. This is more evident given that the PSI displacement rates do not exceed the obtained uncertainties nearby these fault zones and in the vicinity of regions of high strain accumulation, as indicated Focusing on more tectonic-related signals, the lack of linear patterns along the fault zones at the boundary of the basin does not allow for the correlation of any InSAR signal to the interseismic long-term slip-rates of these tectonic structures. This is more evident given that the PSI displacement rates do not exceed the obtained uncertainties nearby these fault zones and in the vicinity of regions of high strain accumulation, as indicated by the geodetic GNSS analysis. Further investigation by means of spatial profiles across fault zones clearly shows maximum motion occurring, not at the edge of the mountain ranges, but relatively shifted towards the center of the sedimentary basin ( Figure 11). Deviations between measured PSI displacement rates from different viewing geometries should be, at least partially, related to the contribution of horizontal components, especially at the margins of regions showing relatively high subsidence rates.
2, x FOR PEER REVIEW 21 of 29 Figure 11. Spatial profiles (a-e and f-k) showing LoS displacement rates across major fault zones from both ascending and descending satellite tracks (see Figures 6 and 7). For the calculations, average motion grids of 100 m spacing were used. Bars indicate the variability of motion within a buffer of 1-2 km, whereas dots the average calculated rates along the profiles.

Discussion
The combination of different techniques and methods is crucial for better assessing the tectonic regime of a region, as different aspects are enlightened. Initially, the GNSS data analysis indicates a low strain pattern for the broader Strymon region, except for a few locations, mostly lying at the eastern boundary of the Strymon basin. These data compare well with existing geological data derived from palaeoseismic trenching, which similarly highlight the most significant tectonic activity at the eastern part of the study area. Regarding the InSAR analysis, no horizontal rates (e.g., below resolution) have been traced in the E-W direction, whereas limited vertical motions have been recorded at the eastern margin. Although GNSS, INSAR, and geological data may partially differ, since they cover different time periods, they all indicate that the area is characterized by rela- Figure 11. Spatial profiles (a-e and f-k) showing LoS displacement rates across major fault zones from both ascending and descending satellite tracks (see Figures 6 and 7). For the calculations, average motion grids of 100 m spacing were used. Bars indicate the variability of motion within a buffer of 1-2 km, whereas dots the average calculated rates along the profiles.
Hence, apart from the co-interpretation of independent results from different SAR geometries, the investigation of the actual motion components (vertical and E-W) provides valuable insights on the nature of motion and the corresponding deformation sources. It is worth noting that most of the deformation patterns recognized in the different SAR geometries are maintained in the vertical motion component, whereas motion in the E-W direction is practically insignificant, with rates ranging within 1-2 mm/yr ( Supplementary Figures S4 and S5).
The absence of a pronounced tectonic-related signal in the InSAR measurements poses no requirement for processing complementary satellite tracks (i.e., different satellite orbits) or PSI reprocessing at a higher spatial resolution. This is also valid since a majority of the surface discontinuities (i.e., basin margins and terrain slopes) are well sampled by the already acquired results.

Discussion
The combination of different techniques and methods is crucial for better assessing the tectonic regime of a region, as different aspects are enlightened. Initially, the GNSS data analysis indicates a low strain pattern for the broader Strymon region, except for a few locations, mostly lying at the eastern boundary of the Strymon basin. These data compare well with existing geological data derived from palaeoseismic trenching, which similarly highlight the most significant tectonic activity at the eastern part of the study area. Regarding the InSAR analysis, no horizontal rates (e.g., below resolution) have been traced in the E-W direction, whereas limited vertical motions have been recorded at the eastern margin. Although GNSS, INSAR, and geological data may partially differ, since they cover different time periods, they all indicate that the area is characterized by relatively low extensional deformation that is predominantly localized in the eastern part of Strymon basin.
Spatially detailed information on the vertical motion component by InSAR indicates that relatively higher subsidence rates are limited within the quaternary deposits of the Strymon basin and predominantly towards the immediate hanging wall of the Serres fault zone. These correspond to localized patterns not easily attributed to tectonics, but mostly to groundwater extraction, while their temporal behavior that reflects the seasonal variability of the groundwater table indirectly confirms the above statement ( Figure 10).
Regarding the average horizontal motion of the study area, the V results are in accordance with the geodynamic regime of the broader Aegean region. In fact, the northmost Aegean region, where our study area belongs to, shows a slower ca. southward motion compared to the central and southern parts of the Aegean, where the velocity gradually, but also rapidly, increases. The V results show consistency with other geodetic investigations conducted in the wider Aegean region [7,14,107].
After correlating the geodetic strain parameters with the major active tectonic structures, it becomes clear that the highest extension occurs where the E-W-trending Serres fault and the NW-SE-trending Nea Zichni fault are located. Both structures are recognized as normal faults (Serres as a pure normal fault and Nea Zichni as a predominant normal fault with a possible minor strike-slip component). The Serres fault has a clear imprint on the topography, as well as in the hydrographic network. On the other hand, the Nea Zichni fault has also been studied through palaeoseismic trenching, confirming its Holocene activity and showing that deformation is widely distributed in several parallel strands.
To a lesser but certain degree, geodetic analysis marks some tectonic activity at the western Strymon basin's margin. In particular, moderate values of geodetic parameters are documented in the northwestern part of the study area, south of Kerkini Lake and NW of the village of Zevgolatio. Within this area, a dense concentration of epicenters is also shown in Figure 1, completing the geodetic results. Until now, no active tectonic structures had been identified or mapped in this part of the Strymon basinal margin. Along the same margin, further to the southeast, NE-dipping marginal faults do occur. All evidence and indications suggest the northwestward continuation of the marginal faults up to Lake Kerkini with the missing 'Zevgolatio fault'. However, one may wonder whether these results and observations belong to an antithetic fault along the opposite side of the basin. The possibility of the occurrence of such a fault along the northeastern margin of the Strymon basin cannot be excluded; on the contrary, it is quite probable that the entire eastern margin of the basin, from Nea Zichni to Sidirokastro and the Belles fault zone, is tectonically controlled. Therefore, considering the geodetic parameters' spatial distribution, as well as the local-scale tectonic setting (i.e., the E-W, south-dipping Belles-Kerkini fault zone, located north of Kerkini Lake), it is assumed that a potential, non-identified, antithetic, normal dip-slip fault occurs in this part, forming a graben with the (antithetic) Belles-Kerkini fault zone; however, this significant geodetically derived strain may also correlate with the NW-SE-trending boundary of the Pliocene basin, indicating potential tectonic activity. On the other hand, limited tectonic activity is imprinted in the geodetic results of the western-southwestern part of the study area. Whatever the answer, we are afraid that it cannot be resolved with current data. Future studies might examine further the possible causes of the moderate strain we have measured.
Given that the active faults in our area of interest mostly occupy the basinal margins and sometimes cross the inner parts of the basins, it is almost impossible to detect the full extent of these structures within the basins with simple geological observations. The main reasons for this are three-fold: (i) the sediments that fill the basins can quickly cover any recent fault trace, (ii) the geomorphic markers that faults leave are prone to erosion since the affected formations are recent unconsolidated sediments, and (iii) the agricultural activities that wipe out any fault trace. Fault detection in such conditions requires costly surveys, usually geophysical and/or palaeoseismological ones, which also need some indications to applying them at the proper location. The high deformational values around the Serres fault and their western, basinward extension strongly suggest that the fault trace is further prolongated westwards for approx. 9 km. This is a nice example of how geodetic data can provide key information on active fault lengths. In terms of seismic hazard assessment (SHA), such information is crucial for the earthquake potential of a seismic source. Especially for the deterministic assessment of a seismic hazard (DSHA), earthquake magnitude depends on the fault length and ground motion simulation of a site, which can also be affected by the parameters of the seismic source.
Comparing our strain results with the previous studies', they show adequate consistency both in the extension values and in the documented direction of extension, which varies between N-S [16,52] and NE-SW [7,8,17,51,126,127]. In addition, our results show common characteristics with the local work of Mouslopoulou et al. [58] despite differences in the input data such as the GNSS stations' density and networks and the recording interval. To be more specific, we improved the spatial resolution using 16 permanent GNSS stations of five different networks in a smaller area (covering approx. 47,000 km 2 ), in contrast with the 20 permanent GNSS stations of one network scattered within an area that is double the size (approx. 93,000 km 2 ) used by Mouslopoulou et al. [58]. Furthermore, we extended the observation time interval, as our geodetic datasets correspond to a 7-year (84 months) time period, implementing a 30 s-rate (24-h operation) recording, while Mouslopoulou et al. [58] implemented datasets from a 5.5-year (66-month) time period; the solution of each month was based on 48 h of continuous recording. Overall, this study offers a higher spatial resolution in the area surrounding the Strymon basin and a longer but different temporal window.
Despite the consistency, a major difference is spotted at the western boundary of the Strymon basin. According to Mouslopoulou et al. [58], this area accommodates 3.3 mm/yr deformation. Our GNSS, InSAR, and active fault data agree that this part is either inactive or most probably of negligible or low deformation (less than 1 mm/yr), and that strain is mostly accommodated towards the eastern part of the basin. Moreover, Mouslopoulou et al. [58] also referred to the disarray between their geodetic high slip-rate that would require 35-40 m of Holocene displacement with the lack of evidence of recent faulting along the western boundary of the basin. Indeed, such high rates should have exerted a major influence on the stratigraphy, geomorphology, and probably (micro-)seismic activity of the study area. Finally, we believe that, if more active faults were incorporated in their modelling (Sochos, Serres, Nea Zichni, etc.), then the strain would not have localized into one boundary fault, but rather be distributed along several structures.

Conclusions
Space geodetic techniques bring to light how much and in which manner the upper crust is deformed, and, consequently, how tectonic processes are occurring at present. The application of this analysis in the broader Strymon basin, which is based on the qualitative and quantitative calculation of specific geodetic strain parameters (MaHE, MiHE, MaxSS, and AS), strongly suggests that the study area is characterized by low predominantly extensional deformation, an outcome that is also supported by the poor seismic activity and the low slip-rate faults. Hence, considering the GNSS results, combined with the InSAR analysis, the following concluding remarks are drawn:

1.
A dominant N-S (up to NE-SW) direction of extension is documented by the geodetic parameters, which is responsible for the major E-W (up to NW-SE) fault zones observed throughout the Strymon basin; 2.
The eastern margin of Strymon basin includes the most tectonically active structures (e.g., Serres and Tholos-Nea Zichni dip-slip, normal fault zones) of the study area; the Belles-Kerkini dip-slip, normal fault zone, located at the northwestern part, may be considered potentially active. It is noted that the Serres fault zone accommodates the greatest strain within the study area; 3.
The southwestern margin of the Strymon basin shows much lower, but not negligible, deformation. This deformation could be associated with the low and scattered seismicity in this area. According to the deformation indices, shear faulting is the expected kinematics; 4.
GNSS data indicate some moderate strain rates at the northwest boundary of the Strymon basin (southeast of Lake Kerkini and near Zevgolatio), implying that some are highlighted by both deformational parameters and seismic activity; 5.
The strain parameters suggest that the Serres fault zone should be further extended to the west, indicating that its length has been underestimated; 6.
The geotectonically related GNSS results (V parameter) highlight the SSE motion of the study area, being consistent with the Aegean microplate motion, which moves towards the Hellenic Subduction Zone; 7.
No horizontal rates have been traced by InSAR in the E-W direction along the study area; if they do exist, they would be below resolution (less than 2 mm/yr); 8.
Vertical InSAR subsidence of a few mm/yr is traced in the immediate hanging wall of the Serres fault zone, within the Quaternary basin; however, it is mainly attributed to non-tectonic effects such as groundwater extraction; 9.
The detection of potential active faults in areas where geodetic strain is recorded, but active faults have not been mapped, or the prolongation of fault traces in areas where local conditions prevent direct observations, highlight the essential role of satellite geodesy and InSAR in the process of seismic hazard assessment, as it may decisively contribute to the determination of fault properties.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/app12189391/s1, Figure S1: C (vertical component), for GNSS sites located over the InSAR processing extend. The closest PS target to the GNSS sites, within a search radius of 200 m, was considered. The corresponding measurements' uncertainties for each technique are also shown (error bars).; Figure S2: Uncertainties of PSI displacement rates for Sentinel-1 ascending track 102 (2015-2020); Figure S3: Uncertainties of PSI displacement rates for Sentinel-1 descending track 7 (2015-2020); Figure S4: Vertical displacement rates for the period 2015-2020 at 200 m spatial resolution, as derived by combining ascending and descending PSI measurements. Positive values correspond to uplift, whereas negative values to subsidence; Figure S5: E-W displacement rates for the period 2015-2020 at 200 m spatial resolution, as derived by combining ascending and descending PSI measurements. Positive values correspond to motion towards East, whereas negative ones correspond to westward motion.