Spatio-Temporal Assessment of Land Deformation as a Factor Contributing to Relative Sea Level Rise in Coastal Urban and Natural Protected Areas Using Multi-Source Earth Observation Data

The rise in sea level is expected to considerably aggravate the impact of coastal hazards in the coming years. Low-lying coastal urban centers, populated deltas, and coastal protected areas are key societal hotspots of coastal vulnerability in terms of relative sea level change. Land deformation on a local scale can significantly affect estimations, so it is necessary to understand the rhythm and spatial distribution of potential land subsidence/uplift in coastal areas. The present study deals with the determination of the relative vertical rates of the land deformation and the sea-surface height by using multi-source Earth observation—synthetic aperture radar (SAR), global navigation satellite system (GNSS), tide gauge, and altimetry data. To this end, the multi-temporal SAR interferometry (MT-InSAR) technique was used in order to exploit the most recent Copernicus Sentinel-1 data. The products were set to a reference frame by using GNSS measurements and were combined with a re-analysis model assimilating satellite altimetry data, obtained by the Copernicus Marine Service. Additional GNSS and tide gauge observations have been used for validation purposes. The proposed methodological approach has been implemented in three pilot cases: the city of Alexandroupolis in the Evros Delta region, the coastal zone of Thermaic Gulf, and the coastal area of Killini, Araxos (Patras Gulf) in the northwestern Peloponnese, which are Greek coastal areas with special characteristics. The present research provides localized relative sea-level estimations for the three case studies. Their variation is high, ranging from values close to zero, i.e., from 5–10 cm and 30 cm in 50 years for urban areas to values of 50–60 cm in 50 years for rural areas, close to the coast. The results of this research work can contribute to the effective management of coastal areas in the framework of adaptation and mitigation strategies attributed to climate change. Scaling up the proposed methodology to a continental level is required in order to overcome the existing lack of proper assessment of the relevant hazard in Europe.


Introduction
Sea-level rise (SLR) is one of the most significant effects of climate change that has been the focus of international attention due to its high future projected rates that would impact coastal areas around the world [1]. Particularly, coastal low-lying regions will be affected, and their land will be decreased due to coastal erosion and inundation [1,2]. In the early 1990s, 33.5% of the global population lived within In this framework, MT-InSAR approaches were extensively and successfully exploited to investigate coastal ground deformation [32][33][34]. These techniques, such as permanent scatterers (PSs) [35], SBAS [36], and interferometric point target analysis (IPTA) [37] use a medium-to-large dataset of SAR images acquired at different times over the observed area to follow the temporal evolution of a deformation phenomenon, retrieving the mean deformation rate and the time series for each point target. Moreover, these techniques can benefit from the availability of free and open access SAR data archives provided by European and international space agencies.
Although several studies have been conducted for studying land deformation using MT-InSAR [38][39][40], PSI TERRAFIRMA products [41], and GNSS [42][43][44][45] in Greece, only a few have been targeting the coastal hazard [46][47][48]. Moreover, none of them have used both the aforementioned space geodesy methods in combination with altimetry measurements to define the relative sea-level rise. In this respect, the aim of this study is to improve the understanding of and quantify the coastal relative sea level rise in three coastal regions of Greece by using a combination of multi-source Earth observation-SAR, GNSS, tide gauge, and altimetry data. To that end, MT-InSAR technique will be used by taking advantage of the most recent Copernicus Sentinel-1 data, calibrated by GNSS data and combined with a reanalysis model assimilating satellite altimetry data. Additional GNSS and tide gauge observations have been used for validation purposes. The proposed methodological approach has been implemented in three pilot areas with special characteristics, contributing to a better understanding of the coastal hazards, constraining ground deformation phenomena occurring along urban and natural protected coastal areas.

Study Areas
This paper deals with the spatio-temporal assessment of land deformation, as a factor contributing to relative sea level rise, in three pilot cases with different combinations of characteristics ( Figure 1). The first pilot case is the city of Alexandroupolis extending up to the Evros Delta region (A), which is a combination of coastal urban area (in the West side) and the unique characteristics of the trans-boundary Evros Delta (on the east side) in terms of the importance of the wetlands. The second pilot case is the coastal zone of the Thermaic Gulf (B) which is a combination of a big urban center (Thessaloniki city on the east side) with its port and airport infrastructures and the unique formation of Delta Axios (in the west side) and the Axios-Loudias-Aliakmonas National Park. The third one is the coastal area of Killini (west side)-Araxos (northeast side) (C) in which there are two important infrastructures (port and airport) and the special Lagoon of Kalogria area (northeast side).

Coastal Area of Alexandroupolis City-Evros Delta Region
The coastal area from the Alexandroupolis city extending up to the Evros Delta (Gulf of Alexandroupolis) covers an area of about 350 km 2 [49]. Alexandroupolis, the capital of the prefecture of Evros, is the largest city of Thrace and the region of Eastern Macedonia in terms of size and population. This very fast-growing city has an important port and it is a summer resort centre in SE Balkans due to the great location at the center of sea and land routes connecting Greece with Turkey [50]. The coastal zone area under investigation extends to more than 50 km and can be distinguished in the western unity and the eastern unity in terms of geomorphology. The western unity from the Mesimvria coast up to Alexandroupolis is characterized as hilly to mountainous, whereas the eastern one covering the coastal zone from Alexandroupolis to the Evros river mouth is a plain area relief [50].
The Evros Delta region, shared by Greece and Turkey, has been characterized as one of the most important wetlands on a national and European level [51]. A major part of the delta on the Greek side has been characterized as a wetland in special protection area (SPA), as a site of community importance (SCI) in the Natura 2000 Network, and as internationally important under the Ramsar Convention (1971), due to its significant and rare species of plants, fauna and birds [52]. The Turkish part of the delta is also included in the list of wetlands of international importance, while Lake Gala, in close proximity to the delta, has been declared a National Park area [51]. The delta is formed by

Coastal Area of Alexandroupolis City-Evros Delta Region
The coastal area from the Alexandroupolis city extending up to the Evros Delta (Gulf of Alexandroupolis) covers an area of about 350 km 2 [49]. Alexandroupolis, the capital of the prefecture of Evros, is the largest city of Thrace and the region of Eastern Macedonia in terms of size and population. This very fast-growing city has an important port and it is a summer resort centre in SE Balkans due to the great location at the center of sea and land routes connecting Greece with Turkey [50]. The coastal zone area under investigation extends to more than 50 km and can be distinguished in the western unity and the eastern unity in terms of geomorphology. The western unity from the Mesimvria coast up to Alexandroupolis is characterized as hilly to mountainous, whereas the eastern one covering the coastal zone from Alexandroupolis to the Evros river mouth is a plain area relief [50].
The Evros Delta region, shared by Greece and Turkey, has been characterized as one of the most important wetlands on a national and European level [51]. A major part of the delta on the Greek side has been characterized as a wetland in special protection area (SPA), as a site of community importance (SCI) in the Natura 2000 Network, and as internationally important under the Ramsar Convention (1971), due to its significant and rare species of plants, fauna and birds [52]. The Turkish part of the delta is also included in the list of wetlands of international importance, while Lake Gala, in close proximity to the delta, has been declared a National Park area [51]. The delta is formed by the Transboundary Evros River alluvial deposits, filled rapidly over millennia and affected by the interaction with the sea [53].

Coastal Zone of Thermaic Gulf, City of Thessaloniki-Axios Delta
The City of Thessaloniki, located in Central Macedonia and situated in the inner part of the Thermaic Gulf, is established as the second most important urban center of Greece. It has an extended industrial zone in its suburbs and an international port that constitutes the major center of merchant shipping for the surrounding Balkan countries [54]. The catchment area of the Thermaic Gulf, located in the southern Balkan Peninsula, is approximately 40,000 km 2 , and the main rivers are Axios, Aliakmon, Loudias, and Gallikos [55,56]. The Plain of Thessaloniki is formed by the sub-aerial deltaic plains of these rivers, together with River Loudias and the artificially drained lake of Giannitsa [55].

Coastal Zone of Thermaic Gulf, City of Thessaloniki-Axios Delta
The City of Thessaloniki, located in Central Macedonia and situated in the inner part of the Thermaic Gulf, is established as the second most important urban center of Greece. It has an extended industrial zone in its suburbs and an international port that constitutes the major center of merchant shipping for the surrounding Balkan countries [54]. The catchment area of the Thermaic Gulf, located in the southern Balkan Peninsula, is approximately 40,000 km 2 , and the main rivers are Axios, Aliakmon, Loudias, and Gallikos [55,56]. The Plain of Thessaloniki is formed by the sub-aerial deltaic plains of these rivers, together with River Loudias and the artificially drained lake of Giannitsa [55]. The wetlands of the Axios Delta at the Thessaloniki plain provide a typical example of wetland destruction in Greece. In 1917, 36% of the plain was wetland, but this area now amounts only to 5.5% [56]. This is also confirmed by Psimoulis et al. [57], who stated that the areas of Kalochori and Sindos used to be a delta some thousand years ago. The main environmental pressures such as water discharge decrease, drainage works, urbanization, and pollution negatively affected the ecological character of the deltaic area, leading to the destruction of 70% of the original wetlands during the 20th century [56]. Moreover, the decrease in rainfall in combination with the overuse of water for irrigation, has resulted in severe salinization of the delta area, which has impacted on the flora and fauna of the wetlands [58]. Currently, some of these activities have ceased, and their impacts have already been mitigated [56].

Coastal Area of Killini-Araxos (Patras Gulf), NW Pelloponese
The coastal area from the port of Killini to the Araxos airport is part of the Gulf of Patras and is located in the northwestern Peloponnese. The Peninsula of Killini, located at the westernmost end of Peloponnese, contains an isolated hilly area of about 130 km 2 . It is formed by a morphological rise (Kastron at a highest point of 244 m) connected to mainland Peloponnese in the east by the plain of the Pineios River (Elis plain) [59]. Araxos is a village in the northwestern part of Achaea that is located in the coastal plains near Cape Araxos, which separates the Gulf of Patras from the Ionian Sea. From Araxos to Killini, there are three lagoons-Prokopos, Kalogria, and Kotychi. Part of the area was recognized as a Wetland of International Importance in 1975, when it was included in the 10 wetlands of Greece protected under the Ramsar Convention. A few decades later, parts of the area were recognized as SPA for bird species, in the framework of Directive 2009/147/EE, as well as SCI in line with Directive 92/43/EEC, which led to the formation of the European NATURA 2000 Network of protected areas [60]. Most of the coastal area lies on sand dune formations, while small areas near Kounoupeli, in the eastern area of the Kalogria lagoon and the Mavra Vouna hill, are all composed of hard limestone. The eastern part of the area behind the dunes is covered with clay deposits with depths of a few centimeters to more than 2 m. The seashore consists of unconnected single-grained medium and fine-sized sand with a very small amount of silt [61].

Materials and Methods
In this section, the data and methods that were used for this study are presented. More specifically, data specifications are given for all the satellite SAR, altimetry, GNSS, and tide gauge measurements. The mature technology of the MT-InSAR was used for this research work to provide land deformation rates which were set to a reference frame with the use of a GNSS data. The deformation rates were combined with the altimetric SSH rate to calculate the relative sea level change. Additional GNSS and tide gauge observations were used for validation purposes. The workflow of the proposed methodological approach is illustrated in the flowchart of Figure 2.

Data
This study, apart from coastal urban centers, focuses on natural protected areas. For this reason, the appropriate selection of sensors is important in order to achieve high spatial coverage of surface

Data
This study, apart from coastal urban centers, focuses on natural protected areas. For this reason, the appropriate selection of sensors is important in order to achieve high spatial coverage of surface deformation information. In this case, high temporal resolution is a crucial parameter. Copernicus satellites Sentinel-1 A and B show high spatial and temporal resolution operating at C band. Today, the latest and most advanced SAR mission is the Sentinel-1 constellation. The Sentinel-1A, the first of the twin satellites, was launched in April 2014 by ESA. The Interferometric Wide (IW) swath acquisition mode, has a swath of 250 km and spatial resolution of 5 m in range (dimension perpendicular to the satellite track) and 20 m in azimuth (dimension along the satellite track), with no multilook and a repeat cycle of 12 days. In April 2016, the Sentinel-1B was also successfully placed into orbit, decreasing the temporal resolution for the constellation to six days [62].
In the TOPSAR operation mode (the most common of Sentinel-1), in addition to steering the beam in range as in ScanSAR (the most common mode of ESA's ASAR/ENVISAT and ERS missions), the beam is also electronically steered from backward to forward in the azimuth direction for each burst, resulting in an homogeneous image quality throughout the swath [62]. With this mode, the same coverage and resolution as in ScanSAR is achieved, but with an improved signal-to-noise ratio (SNR). Sentinel acquisitions are available systematically and free of charge through the Copernicus Open Access Hub (https://scihub.copernicus.eu/dhus/#/home).
For the case of Alexandroupolis, Evros Delta, a set of 91 ascending Sentinel-1B, IW, Level-1, single look complex (SLC) acquisitions between 27 September 2016 and 30 October 2019 (Table 1) was used. For the case of the wide area of the Thermaic Gulf, Thessaloniki, Axios Delta, a set of 130 ascending Sentinel-1 A and B, IW, Level-1 SLC acquisitions was collected. The set of Sentinel-1 images was acquired along ascending passes covering the period 12 October 2014 until 24 June 2019 with 70 Sentinel-1A images and 60 Sentinel-1B images (Table 1). Finally, for the case of the coastal area of Killini, Araxos in the northwestern Peloponnese, a set of 98 ascending Sentinel-1A and 60 Sentinel-1 B, Level-1, SLC, IW acquisitions was used. The set of 158 images spans the period from 11 November 2015 to 23 July 2019 (Table 1). For this research work, the monthly mean "sea-surface height above sea level" of "MEDSEA_REANALYSIS_PHYS_006_004" product at 0.042 degree spatial resolution from Copernicus Marine Environment Monitoring Services (http://marine.copernicus.eu/) was used. This product is the result of a physical re-analysis component from the Nucleus for European Modelling of the Ocean (NEMO) assimilating satellite data. The assimilated dataset includes mono altimeter satellite along-track SSH computed with respect to a seven-year mean [63].

Methodology
The aim of the persistent scatterers technique is to overcome interferometric coherence degradation over time using a set of phase-stable pixels (the persistent scatterers). This technique is performing better in urban areas, where temporal decorrelation is minimized due to the high density of stable structures acting as scatterers (buildings, bridges, etc.). The flowchart of the methodology is illustrated in Figure 2.
MT-InSAR deformation rates from Sentinel-1 data alone suffer from lack of anchor or reference point. A reference PS could be selected in an area located on limestone bedrock with high coherence, being away from deltaic and fluvial deposits or recent seabed outcropped deposits and generally from brittle material. All three GNSS solutions were calculated from NOA GNSS processing station (NGProS) (http://aips.space.noa.gr) using version 6.3 of GIPSY/OASIS II software (http://gipsy-oasis.jpl.nasa.gov), developed by the Jet Propulsion Laboratory of the National Aeronautics and Space Administration (NASA), Pasadena, California, CA, USA [66]. The orbits of the GNSS satellites used were the precise ones in the reference frame system of ITRF2014. The ocean load model that was used is the Goddard/Grenoble Ocean Tide (GOT 4.3) empirical model maintained by NASA-Goddard Space Flight Center. The ocean load values for the three GNSS was calculated from the "free ocean tide loading provider" (http://holt.oso.chalmers.se/loading/). The vertical rate for each GNSS station was calculated using the least squares method ( Table 2). The daily time series of the three GNSS stations were calculated and the vertical linear rate values were estimated as shown in Table 2 and Figure S1. In the case of both Alexandroupolis-Evros Delta and Kyllini-Araxos, the Parallel SBAS (P-SBAS) [67,68] service under the Geohazards Exploitation Platform (GEP) of ESA was exploited. GEP is a cloud platform that provides a rich set of ready to use EO Data processing services for geohazards analysis and monitoring (https://geohazards-tep.eu/). A threshold of temporal coherence of 0.8 was used. As reference points, the GNSS stations ALEX and RLSO were used accordingly. The P-SBAS service is a parallel version of SBAS algorithm. It was selected for these two cases due to its higher effectiveness in rural areas in comparison with PSI which is more powerful in urban areas such as the Thermaic Gulf surrounded by urban centers.
In the case of the Thermaic Gulf, SAR PROcessing tool by periZ software (SARPROZ, Razer Limited, Hongkong, China) [69] was used as the MT-InSAR processor, extending the standard linear PS technique to estimate non-linear motion having no prior information. Shuttle Radar Topography Mission (SRTM) 1 arc-sec DEM [70] was used for the PSI processing. A master date was chosen for this period based on the minimization of the perpendicular baseline, the temporal baseline and the weather data. After selecting a master image (acquired on 24 March 2017), the proposed algorithm was applied only on pixels showing an inverse amplitude stability (i.e., the ratio between the mean calibrated image amplitudes and its standard deviation) higher than 0.7. The GNSS station AUT1 was chosen as a reference point (Table 2). Once the data was unwrapped, the low pass component was estimated using a three-sample temporal base-line weighted moving average and assumed this residual phase term is an estimation of the non-linear motion contribution. These steps were performed considering only neighboring pixels. The same processing steps were then applied, considering only differences between the reference point and all persistent scatterers including the estimated atmospheric phase screens (APS). Finally, for each PS the linear rate was calculated by using least squares approximation.
All three MT-InSAR datasets are originally in slant range geometry, i.e., along the LOS direction. We assumed that the deformation is totally vertical, a hypothesis that is not valid if an earthquake or slow aseismic slip of tectonic origin or landslide occurs. For the study period and the three study areas according to global Centroid Moment Catalog (https://www.globalcmt.org/) there were no earthquakes greater than Mw = 5.0 and shallower than 15 km within a distance of 40 km. Nevertheless, the earthquake is a sudden event and this kind of deformation has little impact or is filtered out by the MT-InSAR process. The aseismic processes take place in large fault zones (e.g., North Anatolian fault as well as its continuation in the Aegean Sea) or smaller ones (e.g., Gulf of Corinth) and the vast majority of the deformation is concentrated within a narrow zone along the feature. The landslides affect the deformation field only very locally. For the Alexandroupolis city-Evros Delta region the east branch of the north Anatolian fault (a right lateral structure) sits almost 15 km south of the southern coast of Evros Delta. Potentially, spatially horizontal (east-west oriented) deformation due to the continuous deformation of this fault zone is considered negligible due to the long distance involved and the limited spatial coverage of each case. Currently, the distribution of the available GNSS stations, with adequate temporal baselines in the three areas, is not dense enough to differentiate for horizontal deformation velocities. Accordingly, considering the vertical component of the Line of Sight vector for each PS, the MT-InSAR rates were transformed to the vertical ones. In order to anchor them with the GNSS ones the following procedure was used. An average value was calculated, at each of the three GNSS locations, for the PSs (of the corresponding MT-InSAR dataset) closer than 300 m from each location (except for AUT1 that is located in a semi-urban area, where the distance is 600 m). The offset of this value to the GNSS one was subtracted from each one of the three MT-InSAR datasets. Thus finally, all the MT-InSAR products were transformed to vertical direction and anchored to a single reference frame. These products are labelled as "calibrated" in Figure 2. From the SSH, using least squares, we calculated the linear SSH rates as in, e.g., Figure 4. Further, the relative sea-level rise rate was calculated simply by subtracting the vertical deformation rate from the sea-surface height rate. Finally, the additional units of cm per 50 years were used to facilitate the discussion on a basis of longer time periods. For validation purposes, we used the vertical component of GNSS solutions from a Nevada institute as well as tide gauges from PSMSL where they were available. Results of localized deformation rates and thus relative sea-level rise are following in the next section. The GNSS uncertainties are used along with the standard deviation values of the spatially distributed MT-InSAR deformation rates in order to provide a level of uncertainty of each detected localized deformation.

Results
In this section, we present the localized deformations detected by MT-InSAR processing and sea-surface height, following the methodological approach detailed in the previous section (Figure 2) in the three Greek coastal areas under investigation effectively covered by Sentinel-1 SAR acquisitions. The absolute sea level trend across Greece for the period 1993-2019 ranges between 1 mm/year and 3 mm/year and, in some specific areas such as in the southern part of the Peloponnese and in the northeastern part of Crete, to more than 3 mm/year.

Alexandroupolis, Evros Delta
In the case of the Alexandroupolis, Evros Delta region, the velocity of deformation has been derived using 91 Sentinel-1 acquisitions (Table 1) and the persistent scatterers technique over an area of 212 km 2 . The urban areas are adequately covered with natural scatterers, while the lack of them in the vegetated and wetland areas present low coherence (Figure 3). For the period 2016-2019, the characteristic vertical deformation ranges between −1.5 mm/year and 5 mm/year. The city of Alexandroupolis shows relatively stable patterns with a mean rate of 0.4 mm/year, while the Turkish part of the Evros Delta seems to uplift with high rates. In order to investigate this further, the MT-InSAR time series of 199 PSs vertical deformation rate lying in the coast of the Delta (bounded by red dashed polygon in Figure 3), along with its mean values and the linear trend, are plotted in Figure 4. We noted that the seasonal trend is highly visible and that the loading periods of the acquirer (before summer) have a higher gradient than the unloading. Moreover, during the maximum acquirer unloading, the vertical deformation is not reverted to its previous situation but is higher. The mean deformation rate of the 199 PS points is 4.5 mm/year. It is obvious that, for the study period, the rate of the loading/drainage is positive; thus, the superficial land over the acquirer is uplifting. In order to have more concrete results, a wider time range covering drought periods is necessary.
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 21 dashed polygon in Figure 3), along with its mean values and the linear trend, are plotted in Figure 4. We noted that the seasonal trend is highly visible and that the loading periods of the acquirer (before summer) have a higher gradient than the unloading. Moreover, during the maximum acquirer unloading, the vertical deformation is not reverted to its previous situation but is higher. The mean deformation rate of the 199 PS points is 4.5 mm/year. It is obvious that, for the study period, the rate of the loading/drainage is positive; thus, the superficial land over the acquirer is uplifting. In order to have more concrete results, a wider time range covering drought periods is necessary.     The SSH rate for the coastal area of Alexandroupolis, Evros Delta is~2.1 mm/year ( Figure 3). For this case, we considered six sub-areas-two are bounded by polygons, and four are labelled by their place names. For each one, in Table 3, the statistics of the PSs (i.e., number, mean, minimum, maximum, standard deviation, with confidence level of 85%) were calculated. Additionally, the ALEX GNSS station (set as reference) uncertainty of 1.0 mm/year was added to the standard deviation values. The highest sea level rise rates are located at the airport with a value of 18.5 ± 18 cm in 50 years and at WaterLand Park (green dashed polygon in Figure 3) with 16.5 ± 21 cm in 50 years. Alexandroupolis presents a uniform behavior, with 8.5 ± 11 cm in 50 years. Finally, the areas south of the Delta (i.e., Yayla Sahili and Erickli Plaji) present low values of sea level rise (1.5 ± 10 cm in 50 years).
According to the Nevada Geodetic Laboratory, the vertical deformation rate of the GNSS station, located at the west end of the city center ALE3 (Figure 3), near the ring road, has a value of 1.02 ± 8.74 mm/year. Uncertainty is very high due to the missing data of almost four years, between 2016 and 2020. The mean vertical deformation value of the PSs with a distance less than 200 m from the ALE3 location is 0.4 mm/year, a value close to the GNSS one.
The PSMSL tide gauge (Table 3) located in the harbor (red circle in Figure 3) measures a relative sea level trend of 1.84 ± 0.67 mm/year (with 95% confidence) based on the monthly mean sea level data between 1969 and 2017. The mean vertical deformation rate of selected PSs located in the same deck as the tide gauge is −0.1 mm/year, which corresponds to a relative sea level rate of 2.2 mm/year (or 11 cm in 50 years), a value close to the tide gauge measurement (9 cm in 50 years). Table 3. Statistics of the sub-areas of Alexandroupolis, Evros Delta. Amongst others, the characteristic (mean) values of vertical deformation rate, the relative sea-level rise, and the uncertainty are shown. The sub-area Coastal Delta zone presents negative sea-level rise, thus sea-level falling. PS stands for permanent scatterers; GNSS for Global Navigation Satellite System; RSL for relative sea level.

Location
Nr

Thermaic Gulf; Thessaloniki City, Axios Delta
In the case of the Thermaic Gulf, the velocity of deformation was derived using 130 Sentinel-1 images (Table 1) and MT-InSAR technique over an area of 142 km 2 . The land deformation map ( Figure 5) shows a very high-density concentration of PSs in the area of Thessaloniki and in the surrounding urban centers. Outside, as in the agricultural fields and the mountains, there is a low density of PSs. A factor indicating the limits of this technique [71]. According to the Nevada Geodetic Laboratory, the vertical deformation rate of the GNSS THS1 station, located in the Aristotle University of Thessaloniki (Figure 5), has a value of 0.2 ± 0.58 mm/year. The mean vertical deformation value of the PSs with a distance less than 200 m from the THS1 location is −0.4 mm/year, a value close to the GNSS value. The PSMSL tide gauge located in the harbor (red circle in Figure 5) measures a relative sea level trend of 3.83 ± 0.66 mm/year (with a 95% confidence) based on a monthly mean sea level data from 1969 to 2017. The mean vertical deformation rate of selected PSs located in the same deck as the tide gauge is −2.4 mm/year which corresponds to a relative sea level rate of 4.3 mm/year (or 21 cm in 50 years), a value very close to the tide gauge measurement.
Raucoules et al. [54] using permanent scatterers and stacking interferometry and exploiting 47 ERS acquisitions for the period 1992-2000 produced a deformation map in the slant range (LOS) direction, including a large part of the map of this study. There is an overall consistency with this study except that the areas of maximum subsidence of 1992-2000, i.e., the airport and Kalochori areas. currently present smaller subsidence. For the former area (in the period of 1992-2000), the subsidence was 10-20mm/year in LOS, and currently, it is 4 mm/year in the vertical direction. The latter area (between 1992 and 2000) was characterized by a subsidence of 50 mm/year, and currently, the maximum values reach 10 mm/year in the vertical direction.
Constantini et al. [72], using combined PSI and SBAS approaches and exploiting a number of 20 ASAR/ENVISAT acquisitions for the period of 2004-2010, produced a deformation map in the slant range (LOS) direction, including a large part of the map of this study. There is also an overall consistency with this study. The area of Kalochori (for the period of 2004-2010) presented a subsidence rate of 15 mm/year in the LOS that is also larger than the current rate, estimated at 4 The PSMSL tide gauge located in the harbor (red circle in Figure 5) measures a relative sea level trend of 3.83 ± 0.66 mm/year (with a 95% confidence) based on a monthly mean sea level data from 1969 to 2017. The mean vertical deformation rate of selected PSs located in the same deck as the tide gauge is −2.4 mm/year which corresponds to a relative sea level rate of 4.3 mm/year (or 21 cm in 50 years), a value very close to the tide gauge measurement.
Raucoules et al. [54] using permanent scatterers and stacking interferometry and exploiting 47 ERS acquisitions for the period 1992-2000 produced a deformation map in the slant range (LOS) direction, including a large part of the map of this study. There is an overall consistency with this study except that the areas of maximum subsidence of 1992-2000, i.e., the airport and Kalochori areas. currently present smaller subsidence. For the former area (in the period of 1992-2000), the subsidence was 10-20mm/year in LOS, and currently, it is 4 mm/year in the vertical direction. The latter area (between 1992 and 2000) was characterized by a subsidence of 50 mm/year, and currently, the maximum values reach 10 mm/year in the vertical direction.
Constantini et al. [72], using combined PSI and SBAS approaches and exploiting a number of 20 ASAR/ENVISAT acquisitions for the period of 2004-2010, produced a deformation map in the slant range (LOS) direction, including a large part of the map of this study. There is also an overall consistency with this study. The area of Kalochori (for the period of 2004-2010) presented a subsidence rate of 15 mm/year in the LOS that is also larger than the current rate, estimated at 4 mm/year in vertical direction.
The main reason for the subsidence in the area of Kalochori to Kimina is the over-exploitation of the underground water. It is obvious that the subsiding of these areas has been progressively slowed down since 1992. The most likely explanation is that regulations for water drilling were successfully applied. For the detailed study, the whole area was divided into 11 sub-areas, as shown in Figure 5. These sub-areas were selected as separate locations with a similar deformation rate near the coast or in low altitudes (in case of PSs absence in the coastal zone). For each sub-area, in Table 4, the statistics of the PSs (i.e., number, mean, minimum, maximum, and standard deviation) were calculated. Additionally, the GNSS uncertainty of 0.5 mm/year was added to the standard deviation values. Finally, considering the sea-surface height rate, the relative sea level rise was calculated, together with their uncertainty values. We observe that the maximum values are located at the villages of Chalastra and Kimina (both at an altitude of~5 m) with a relative sea level rate of 67 ± 23 and 39 ± 14 cm in 50 years. Kalochori. Airport and Perea sub-areas present a rate of 31 ± 17 cm in 50 years. All the others have a rate less than 30 ± 9 cm in 50 years. The lowest values are located in Kalamaria. being 17 ± 7 cm in 50 years.

Killini, Araxos
For the case of Killini, Araxos, in the northwestern Peloponnese, the velocity of deformation ( Figure 6) was derived using 158 Sentinel-1 images (Table 1) and the SBAS technique over an area of 129 km 2 . For the time period of 2015-2019, the mean vertical deformation rate for the coastal areas ranges between −3.6 mm/year and 3.8 mm/year. The areas between P. Kalamakiou and Brinia are subsiding, and the rest are uplifting.
The area was divided into five sub-areas. For each one, in Table 5, the statistics of the PSs (i.e., number, mean, minimum, maximum, standard deviation, with a confidence level of 85%) were calculated. Additionally, the GNSS uncertainty of 1.1 mm/year was added to the standard deviation values. Finally, considering the sea-surface height rate, the relative sea level rise was calculated, with their uncertainty values. Specifically, in sub-areas from K. Achaea to P. Kalamakiou and Ioniko, there is a small relative sea level rise. The sub-areas presenting the highest sea level rise rates are P. Kalamakiou to L. Kalogria and Strofilia Forest to Briana, with values of 20 ± 14 cm and 28 ± 21 cm in 50 years. Finally, at Killin, the relative sea level rise is 6 ± 14 cm in 50 years. Table 4. Statistics of the sub-areas of the Thermaic Gulf; Thessaloniki, Axios Delta. Amongst others, the characteristic (mean) values of vertical deformation rate, the relative sea level rise, and the uncertainty are shown. PS stands for permanent scatterers; GNSS for Global Navigation Satellite System; RSL for relative sea level.  The area was divided into five sub-areas. For each one, in Table 5, the statistics of the PSs (i.e., number, mean, minimum, maximum, standard deviation, with a confidence level of 85%) were calculated. Additionally, the GNSS uncertainty of 1.1 mm/year was added to the standard deviation values. Finally, considering the sea-surface height rate, the relative sea level rise was calculated, with their uncertainty values. Specifically, in sub-areas from K. Achaea to P. Kalamakiou and Ioniko, there is a small relative sea level rise. The sub-areas presenting the highest sea level rise rates are P. Kalamakiou to L. Kalogria and Strofilia Forest to Briana, with values of 20 ± 14 cm and 28 ± 21 cm in 50 years. Finally, at Killin, the relative sea level rise is 6 ± 14 cm in 50 years.

Discussion
In the present work, the relative sea level rise rate was estimated for several coastal areas within the three selected case studies. The MT-InSAR vertical deformation rate data (exploiting Copernicus Sentinel-1 acquisitions) is referenced to the vertical component of GNSS measurements (calculated in ITRF2014 reference frame system) in order for their rate to be "calibrated". Additionally, GNSS measurements not calculated by this study were used for validation purposes. A reanalysis model (Copernicus Marine Environment Monitoring Services) assimilating satellite altimetry dataset that provided sea-surface height time series was used for extracting the respective rate. The "calibrated" MT-InSAR vertical ground displacement rates were combined with the SSH rates to estimate the relative rates. Combined maps of land deformation and SSH rates show the spatial distribution of these measurements. The coastal areas are segmented into sub-areas of similar behavior and their statistics, characteristic velocities, and finally the relative sea level rise rate (expressed in cm over a period of 50 years) are estimated, along with uncertainty values. Tide gauge measurements, which were available from PSMSL network, are in accordance with our findings.
In the three case studies the sea-surface height rate is varied between 1.9 mm/year off-shore the Thermaic Gulf and 2.1 mm/year off-shore the Evros Delta. The characteristic relative sea-level rise rates vary between 1.5 ± 9 and 67 ± 23 cm in 50 years (Figure 7). The case study presenting the smallest characteristic rates is Alexandroupolis and Evros Delta with its values varying between 1 ± 10 and 18.5 ± 18 cm in 50 years (Figure 7). The coastal area of Evros Delta in the southwest seems to be currently in a phase of positive aquifer loading due to high precipitation, presenting high uplifting rates and a relative sea level fall of 16 ± 14 cm in 50 years. Longer temporal study periods are needed to ensure more accurate results in the Delta (on-shore). The highest sea level rise rates are located at the airport and at the Waterland Park northwest of the Delta. In Alexandroupolis, the rate is lower at 9 ± 11 cm in 50 years. In the southeast, the two areas of Evros Delta reach the lowest values. In Killini, Araxos, the relative sea rise rate varies between 6 ± 14 and 28 ± 21 cm in 50 years (Figure 7) in the sub-areas of Killini and Strofilia Forest-Briania, respectively. Sea level fall is observed in two sub-areas. Finally, in the Thermaic Gulf, the highest rates observed reach 67 ± 23 and 39 ± 14 cm in 50 years in the villages of Chalastra and Kimina, respectively. The Kalochori. Airport and Perea sub-areas have a rate of~30 ± 14 cm in 50 years. The city center presents a rate between 14 ± 7 (Kalamaria) and 22 ± 5 cm in 50 years (P. Paralia) characteristic rates is Alexandroupolis and Evros Delta with its values varying between 1 ± 10 and 18.5 ± 18 cm in 50 years (Figure 7). The coastal area of Evros Delta in the southwest seems to be currently in a phase of positive aquifer loading due to high precipitation, presenting high uplifting rates and a relative sea level fall of 16 ± 14 cm in 50 years. Longer temporal study periods are needed to ensure more accurate results in the Delta (on-shore). The highest sea level rise rates are located at the airport and at the Waterland Park northwest of the Delta. In Alexandroupolis, the rate is lower at 9 ± 11 cm in 50 years. In the southeast, the two areas of Evros Delta reach the lowest values. In Killini, Araxos, the relative sea rise rate varies between 6 ± 14 and 28 ± 21 cm in 50 years (Figure 7) in the subareas of Killini and Strofilia Forest-Briania, respectively. Sea level fall is observed in two sub-areas. Finally, in the Thermaic Gulf, the highest rates observed reach 67 ± 23 and 39 ± 14 cm in 50 years in the villages of Chalastra and Kimina, respectively. The Kalochori. Airport and Perea sub-areas have a rate of ~30 ± 14 cm in 50 years. The city center presents a rate between 14 ± 7 (Kalamaria) and 22 ± 5 cm in 50 years (P. Paralia) Some of the land deformations in this study are due to the changes in aquifer, i.e., loading and unloading procedures as we have seen in the areas of Evros Delta (uplift) and Kalochori to Kimina (subsidence), respectively, resulting from the equilibrium between rainfall and underground water exploitation. Additionally, since the deltaic areas are prone to subsidence due to their alluvial sediments' compaction, we expect a component similar to what has been observed in many deltaic areas [46,48], which increases toward the shore with a rate of~2 mm/year per 1 km [73]. Shallow faulting may be another deformation source especially when the delta is located in the hanging-wall of the fault. In this case, the subsidence is the accumulated effect of both deformation sources [46]. Macroscopically, this type of faulting has not been observed, but a more focused analysis, not being within the scope of the current study, may reveal such concurrent effects.
Future efforts will focus on coupling the land deformation rates with a series of rainfall timelines; more accurate relative sea level rise data; a very high-resolution DEM; land use/land cover maps; and wind, storm surges, and tidal levels through modeling in order to identify the most vulnerable areas and thus to contribute to response and mitigation actions. Moreover, instead of the hypothesis of only vertical deformation that we have assumed in this study, a combination of ascending and descending pass acquisitions can be exploited in order to estimate the horizontal deformation rates [46] and accordingly extract the true vertical component. Therefore, a scenario of future relative sea-level and spatio-temporal assessment of vulnerability of coastal areas could be performed. In a very recent study, Vousdoukas et al. [5] presented economically efficient protection scenarios from rising extreme sea levels in coastal areas of the European continent. This very interesting study in Nature Communications highlighted the fact that there is a lack of proper assessment of this hazard on a continental scale, which reinforces the need for such research work to be implemented on a larger scale.

Conclusions
The combination of measurements. Sensors, platforms, and processing techniques is fundamental for improving the understanding of the ongoing processes along coastal urban and natural protected areas, such as land subsidence and sea-level rise. Satellite data availability and processing methodologies are mature enough to be exploited for relative sea level rise studies. The spatio-temporal assessment of such phenomena could be measured accurately with the available free and openly accessible Copernicus satellite data and the available tools implemented in InSAR processing platforms such as GEP and SARPROZ. The combination of interferometric data with relative sea-level-rise data can be exploited to identify the potential coastal areas at risk and contribute to the decision-making process.
The present research has shown that the variation in the relative sea level in the three case studies is high. ranging from close to zero values i.e. from 5-10 and 30 cm in 50 years for urban areas to characteristic values of 30 and even very locally to 60 cm in 50 years for rural areas close to the coast ( Figure 7). This demonstrates the corresponding high variability in the eastern Mediterranean region, implying the need of inclusion of its coastal mapping amongst other hazard properties (like coastal erosion) as an additional feature in risk analysis processes. The new knowledge generated can contribute to the effective management of coastal areas in the framework of adaptation and mitigation strategies attributed to climate change.
The potential upscaling through the extensive exploitation of MT-InSAR (such as Copernicus European Ground Motion Service, https://land.copernicus.eu/user-corner/technical-library/europeanground-motion-service) and GNSS data will assist the study of more coastal areas in Greece, in the Mediterranean Sea, and also in lower latitude areas where the sea level rise is expected to be higher, as in the case of west Africa. The proposed methodological approach is scalable and can be implemented at local, regional, and international scales. Author Contributions: Conceptualization: I.P. and P.E.; Methodology: I.P., P.E., and G.B.; Software: P.E.; Validation: P.E. and G.B.; Formal analysis: P.E., G.B., and T.P.; Investigation: P.E. and T.P.; Resources: P.E.; Data curation: P.E., G.B., and T.P.; Writing-original draft preparation: T.P., P.E., and G.B.; Writing-review and editing: T.P., P.E., and I.P.; Visualization: P.E., G.B., and T.P.; Supervision: I.P. and P.E.; Project administration: T.P. and I.P.; Funding acquisition: I.P. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by Greece and the European Union (European Social Fund-ESF) through the Operational Programme "Human Resources Development, Education and Lifelong Learning 2014-2020" in the context of the project "Spatio-temporal assessment of land deformation as a factor contributing to relative sea level rise in coastal urban and natural protected areas by using space-based Earth observation data" (MIS 5006912). The APC was funded by Greece and the European Union (European Social Fund-ESF) through the Operational Programme "Human Resources Development, Education and Lifelong Learning 2014-2020" in the context of the project "Spatio-temporal assessment of land deformation as a factor contributing to relative sea level rise in coastal urban and natural protected areas by using space-based Earth observation data" (MIS 5006912). Funding: This research was funded by Greece and the European Union (European Social Fund-ESF) through the Operational Programme "Human Resources Development, Education and Lifelong Learning 2014-2020" in the context of the project "Spatio-temporal assessment of land deformation as a factor contributing to relative sea level rise in coastal urban and natural protected areas by using space-based Earth observation data" (MIS 5006912). The APC was funded by Greece and the European Union (European Social Fund-ESF) through the Operational Programme "Human Resources Development, Education and Lifelong Learning 2014-2020" in the context of the project "Spatio-temporal assessment of land deformation as a factor contributing to relative sea level rise in coastal urban and natural protected areas by using space-based Earth observation data" (MIS 5006912).
Aristotle University of Thessaloniki for providing GNSS station AUT1 data, the National Observatory of Athens for station RLSO, and the Tree Company Greece for station ALEX, used as reference stations for the MTInSAR data, as well as the Nevada Geodetic Laboratory for providing solutions for validation purposes of stations AUT1 and ALE3.