Space-based Deformation Monitoring of Coastal Urban Areas: The Case of Limassol’s Coastal Front

: In the last five years, the urban development of Limassol City has rapidly increased in the sectors of industry, trade, real estate, and many others. This exponentially increased urban development introduces several concerns about the aggravation of the land subsidence in the Limassol coastal front. Fifty Copernicus Sentinel-1 data from 2017-2021 have been processed and analyzed using the Sentinel Application Platform (SNAP) and the Stanford Method for Persistent Scatters (StaMPS). A case study for the identification and analysis of the elements (PS) in pixels in a series of interferograms, and then, the quantity of the land displacements in the Line of Sight, in the Limassol coastal front, is presented in this research, with the subsidence rates up to about (-5 to 4 mm / year). For the validation of the detected deformation, accurate ground-based geodetic measurements along the coastal area were used. Concordantly, taking into account that there are a significant number of skyscrapers planned to be built, this study attempts a preliminary assessment of the impact these structures will pose on the coastal front of the area of Limassol.


Introduction
The notable urban and infrastructure growth of the last five years in Limassol, Cyprus, has determined the necessity to provide a detailed deformation monitoring analysis of land subsidence of the coastal zone. Following the crisis events of 2013 [1]. The Cyprus government promoted incentives for land development to restart the Cyprus economy. Consequently, Limassol arose to be the fastest growing city in Cyprus in construction, with skyscrapers and tall buildings built one after the other along the coastal front, in almost 20km. This massive development attracted tourists and real estate investments generating increased concern that a combination of factors, such as overexploitation of groundwater, the structures' load, the sea level rise driven by the global climate change, and the intense earthquake activity, may holistically trigger land subsidence phenomena with severe impacts [2]. Land subsidence can be defined as the differential sinking of the ground surface with respect to surrounding terrain or sea level. The later poses an imminent threat for the socio-economic equilibrium of the country, as well as an incremental factor for the risk of possible floods in the specific area [3].
The integrated use of multiple space-borne InSAR (Interferometric Synthetic Aperture Radar) techniques are, undoubtedly, among the most effective and terminated methods to monitor land subsidence [4] and, therefore, assess the impact of urban infrastructures on the coastal zones [5]. Techniques such as the Persistent Scatterer Interferometry (PSI) presented by Feretti et.al [6] and Differential Interferometry (D-InSAR) and Small Baseline Subset Interferometry SBAS [6] became indispensable parts in ground deformation monitoring analysis of urban areas as it may provide cm-to mm-level accuracy products [8]. Investigating the land subsidence in an urban environment using D-InSAR, may be inadequate due to the lack of coherence in multitemporal images in cities, caused by temporal/geometrical decorrelations [9,10] and atmospheric inhomogeneities [11], as well. Thus, the use of PSI could be characterized as inevitable. PSI technique detects and measures specific points (buildings, stable rocks, roads etc) on the surface of the Earth that are phase-coherent and stable over a period of time [12,13]. The displacements are only measured in the direction of the axis connecting the target and sensor (LoS); this axis has a different orientation in space for ascending and descending orbits. [14]. Only descending orbits in the current research are considered, since ascending ones were affected by various illumination effects arising from SAR images' side-looking scene illumination. Especially in densely built-up areas with huge buildings, large portions of the data can be interfered with by the shadowing, layovering, and foreshortening effects due to the slope orientation, the geometry, and height of the buildings in an urban environment, as well as the double bouncing of the radar signal [15].
A plethora of studies has proved the ability of land subsidence monitoring in coastal urban areas in many nations worldwide, combining different space-based techniques, such as PSI and SBAS, using various satellite constellations. As shown in Error! Reference source not found., Sentinel-1 provided by ESA satellite mission is most usable in the past few years. Specifically, almost all cities are faced with land subsidence worldwide, consisting of alluvial soil, and the subsidence rates are ranged between mm to cm/year. Utilizing the technology and analyzing the findings so far on these landslides, it can be relatively concluded that many coastal cities are highly affected and may even be at high risk in the future. Using space-based deformation monitoring techniques to investigate the land subsidence occurring along the Limassol Coastal Front, the purpose of this work is twofold. As the first objective, this research focuses on quantifying the land subsidence of the area using remote sensing techniques and validating the detected deformation, if any, using conventional in-situ ground-based data. As a secondary objective, this research investigates the possibility of automatically monitoring land subsidence in urban coastal areas, consisting of skyscrapers and tall buildings through space-based techniques. Concordantly, taking into account that there is a significant number of skyscrapers that is planned to be built, this study attempts a preliminary assessment, characteristics and discussion of the impact these structures will pose on the coastal front of the area of Limassol through the Cyprus Continuously Operating Natural Hazards Monitoring and Prevention System (CyCLOPS) strategic research infrastructure unit [36].

Case Study-Limassol Coastal Front
Limassol city belongs to the south part of the island and is the second-largest urban city of Cyprus after Nicosia, covering an area of almost 53km². The Port of Limassol is one of the busiest ports in the Mediterranean transit trade and the largest port in Cyprus. Furthermore, Limassol is the base of Cyprus University of Technology, one of three state universities established in 2004. The coastal front of Limassol City encompasses an area of about 20 km 2 . The area of interest is relatively flat, with an altitude that reaches up to almost 1.5m above Mean Sea Level. From a geological point of view, as the majority of the coastal cities, the geology of the Limassol coastal zone consisted of alluvial soils. As of late 2013, Cyprus has been experiencing a construction regeneration, and Limassol is the focal point, upgrading the image of the city as a cosmopolitan and attractive destination. Since 2013, when legislation was passed encouraging construction, a significant number of skyscrapers and tall buildings are introduced to the Cypriot community each year. The construction industry in Cyprus counts for 7% of the country's GDP. Up to date, there are about 70 buildings that are proposed and/or under construction that will stand taller than 50m when completed, of which 29 are skyscrapers. Specifically, the Error! Reference source not found. presents, in detail, the name of the projects, the height, and the floors, as well as the year of their construction. A recently-completed residential highrise building, One, is the tallest tower in Cyprus, the tallest seafront residential building in Europe, and the 48th-tallest building in the European Union. Among others, the underconstruction project is the 'City of Dreams Mediterranean,' Europe's largest casino resort. Since Cyprus island is located in the Mediterranean fault zone, interacts between the Eurasian and the African plates, a unique site for geodynamic analysis is exhibited. Taking into consideration the geology of the area, the seismicity of Cyprus in the last five years [37], the sea level rise [38] conducted by the thermal expansion of seawater due to ocean warming and water input from land ice melt, the significant vulnerability to climate change and the rapid construction of skyscrapers and tall buildings in the area land subsidence monitoring is an increasingly eminent component of risk-planning and decisionmaking regarding infrastructure. No cases of land subsidence have already been identified over the study area and quantified through various InSAR techniques such as PSI.

Data Processing
Monitoring the land subsidence in the coastal front of Limassol city has been operated using several techniques ranging from traditional leveling to PSI, as shown in Error! Reference source not found..

Satellite Data Processing
Considering the novelty of the research due to the lack of previous investigations in urban areas in Cyprus and the fact that a significant number of buildings is planned to be built, Open Data policies, through multiple hubs, that involve the use of non-commercial satellite data are chosen, in order to guarantee the reusability of the current study approach on larger spatial and time scales. More specifically, Copernicus satellites imagery data, particularly focusing on Sentinel-1 imagery data, are used. This satellite system provides free data with global coverage at a weekly revisit time [39][40][41]. In conjunction with short temporal baselines, the possibility of coherent phase information between the primary and secondary images increases, thus contributing to a potential greater density of Persistent Scatterers (PS) for estimating deformation over time [42]. Processing and analyzing of Copernicus Sentinel-1 data, covering a time interval ranging from March 2016 to March 2021, implemented by using the Sentinel Application Platform (SNAP), the Stanford Method for Persistent Scatters (StaMPS/MTI) [43], and the 'snap2stamps' python workflow. StaMPS/MTI is a Sentinel-1 TOPS data PSI processing open-source software package that is demonstrated in automated mode. At the same time, 'snap2stamps' contains python scripts to automate the pre-processing of Sentinel-1 SLC data and their preparation for ingestion to StaMPS [44][45][46].
Forty-six Sentinel-1B Advanced Synthetic Aperture Radar images acquired by the European Space Agency from the satellite on descending acquisition mode from 2017-2021 were used to identify and calculate the deformation information, as presented in Error! Reference source not found., functioning in C-band (5.4GHz). It is noteworthy that the dataset follows a repeat sequence interval of almost 30 days between every image in order to avoid the decorrelation in the temporal analysis of the baseline. This period was selected for two reasons: Firstly, most of the buildings' development started in late 2016, so the earlier study will be insufficient. Secondly, the benchmarks' network for leveling reasons was designed and implemented by the Laboratory of Geodesy and Hydrographic Surveying of the Cyprus University of Technology through the 2017 year. The sentinel-1 dataset is composed of Single Look Complex (SLC) images with Interferometric Wide (IW) mode, with similar characteristics, creating a time-series sequence over the passing five years. The Interferometric Wide swath (IW) mode, applied in the present research, acquires data with a 250 km swath, with a high spatial resolution of 5 × 20 m, and by using Terrain Observation with Progressive Scans SAR (TOPS), captures one sub-swath in range-direction [47]. Specifically, in the Sentinel-1 scenes observed by the 167 relative orbit, the orbit passed over Limassol in descending acquisition mode. Moreover, according to Kopel et al. [48], the favorable polarization mode for monitoring land subsidence in an urban environment is the VV polarization mode, as is used in the dataset for the specific study since it strongly returns part of the signal to the sensor due to the artificial scatterers of the surface (e.g. buildings).

SNAP-STAMPS PSI Processing
The PSI processing is divided into two independent workflows. The first workflow concerns the preprocessing of the images and the single-master DInSAR processing using ESA SNAP and 'snap2stamps', and the second workflow contains the PSI processing using StaMPS/MTI.
Initially, the image acquired on the 17th of July 2019 was selected and used as the 'master' image in the pre-processing procedure. The master image minimizes the sum decorrelation of all interferograms by eliminating the perpendicular baseline values and the Doppler phenomenon among all images and concordantly maximizing the expected coherence in the interferometric stack [49]. This results in the geometric and interferometric co-registration of the dataset with respect to a common frame (master image). In each of the sub-swaths in the pre-processing workflow, three bursts were separated in the azimuth direction. The orbital information demonstrated accurate satellite position and velocity information following the SNAP image processing. Apart from that, the geocoding of all images was applied, using the SRTM-1m [50]. Equally important was the enhancement of the spectral diversity, as well as the extraction of the topo-phase, and the generation of the interferograms, as shown in Error! Reference source not found.. Subsequently, after preceding the phase noise estimation(SNR) and the phase unwrapping, the Persistent Scatterers (PS) are created. Specifically, the number of the final persistent scatterers was computed in three main stages. Firstly, all the candidate scatterers are identified in the study area. Persistent scatterers are pixels or points characterized by stability over time of the backscatter electromagnetic signal. Following this step, extracting the possible scatterers with a minimum value of coherence threshold equal to 0.80 was held. Due to the density of the buildings in a small distance, the coherence value should be to the maximum extent possible to avoid the matching of the same PSs to more than one buildings [50]. Last but not least, the final scatterers that correspond to tall buildings or skyscrapers in every pixel were isolated. The amount of persistent scatterers is affected by various factors, including the size of the building, the material, and structures of the roof, the orientation of the building relative to radar look direction, and nearby vegetation [51]. The    In the current study, the orthometric height was computed by leveling surveys using the industrial grade digital level Leica LS15 and two fiberglass staffs[Error! Reference source not found.] . The measurements were provided with an invar bar code that allows a resolution of 0.01 mm in electronic height measurements, offering a height accuracy of ±1.0 mm (standard deviation) per 1 km double-run measurements. The procedures conducted before (instrument alignment), during (alternated back and forward sightings), and after ( the orthometric correction) the measurements eliminate the main sources of errors such as rod alignment, the curvature and refraction errors as well as other several systematic errors. The tolerance error for the first order-class II survey of a section of 2 km (maximum) line length is 4 mm√K, and for an entire leveling line is 5 mm √K, where K is the length in kilometer units [52]. The heights of leveling benchmarks of the vertical geodetic network were referred to the Cyprus Geodetic Reference System of 1993 (CGRS93).

Figure 7. Invar and Leveling Benchmarks installed by Laboratory of Geodesy and Hydrographic
Surveying, Cyprus University of Technology.

Combination of the methodologies in GIS Environment
After the data processing is complete, all information mentioned above associated with the building characteristics, PS velocity and displacements, and leveling measurements were imported in ArcGIS Pro software, provided by ESRI, for further analysis.
Firstly, concerning the buildings as a unit, a database was created in order to assemble all possible geographical and descriptive information about each one. Specifically, the database includes the geolocation of the building exported by Google Earth Pro. In addition to this, important information about the height and the floors of the structure was added to the database. Equally important was the starting date of the construction, as well as the ending date, in order to understand the result values of the PS. Apart from that, recent photos of each structure were included in the database. As mentioned before, the results of this study are preliminary; hence, a well-structured and detailed database will provide future studies over time.
As is well known, the PSI technique provides spatially scattered measurements of mm to cm accuracy displacements with spatial density depending on the resolution of each sensor. Taking into consideration the high density of PS on the infrastructures and, therefore, the validity of the PSI results, a more detailed analysis, isolating and analyzing PS separately, was carried out on the structural components. The total number of PSs extracted by STaMPS in Limassol equals 57 916. Ergo, after investigating PS around/on each construction, a buffer zone of 200m that includes PS and structure information around every building was created. The number of 50m is used as an assumption of demonstrating the land subsidence in each building and its surrounding area through PS analysis.
No interpolation method (e.g. Kriging, IDW etc.) was used because the study focuses on the structures and not in the region generally.

PSI Classification and Trends of PSI Results
Following the above methodology, the geographic information (location) of the corresponding PSs and the displacements in LoS in millimeters/year and their statistics were extracted. Concerning the classification of the PSs values, it was held in ten (10) classes for visualization purposes. Each class equally distances in 1mm/year, ranging from -5mm/year to 4 mm/year. For further analysis, the PSs included in a 200m buffer zone were analyzed for each structure and determined the mean value for the displacement in the vertical component in LoS in each case.
As shown in Error! Reference source not found. and Error! Reference source not found., where the location distribution of the PSs in the buffer zone is presented on an excerpt of a map, there is a significant variance in the point density in each case. This is ascribable to the starting date of each construction and the presence of the buildings in all Sentinel-1 scenes. Thereby the buffer zone of 200m around each building was created, including PSs, as presented in Error! Reference source not found.. For each construction, statistics data describing the mean, minimum, and value of the subsidence in LoS, are computed. For the presentation of the PSI results, a categorizing of the buildings into skyscrapers (<100m ) and tall buildings(<100m) according to height component [53] is created.

PSI Results Analysis of Skyscrapers in Limassol Coastal Front
The skyscrapers category compasses Trilogy, The One, The Dream Tower, Limassol Del Mar, The ICON, Sky Tower, MARR Tower, and DTA Tower. As shown in Error! Reference source not found., the land subsidence of the area was affected the most by the constructions in the region of Trilogy. Trilogy construction includes three similar buildings( two facing the sea and one facing the Limassol city) with a mean value of -1.80 mm/year. Comparily with other skyscrapers except the building named ONE, the mean value is over-doubled. The minimum value of land subsidence in LoS equals -4.51 mm/year. Noteworthy is that the third building facing the north of the complex is still under construction. The structures complete the list of the skyscrapers in Limassol Coastal Front, showing that from the starting day of their building up to the present, the land subsidence of the built area reaches a steady-state that is leaned to zero. . PSI geospatial analysis of the tall buildings reveals that the Olympic Residence building, which follows a two-building complex construction, has an average deformation rate up to -1.7 mm/year in the vertical component. Equally notable is that the specific buildings have the past date of construction, among others, the year 2012. The lowering observed minimum value of -4.17 mm/year along all tall buildings corresponds to Olympic Residence. Concerning the Oval and Arcship building, the average displacements on altitude are equal to -0.41 and -0.42 mm/year , respectively.

Final Land Subsidence's Results in Limassol Coastal Front through PSI and leveling techniques
The final results of the study are presented in Error! Reference source not found.. Studying the distribution and the number of the PSs in every buffer zone can be characterized as crucial since it provides information about the sufficiency of the final results. The average number of distribution is 182, ranging between 86 (Dream Tower) and 395 (ICON). Concerning the distribution of PSs in the buffer zone of 200m around each building, this depends on the surrounding area and the construction phase in each case. A similar measurement unit in millimeters is indispensable for the compatibility correctness of the results. Ergo, the multiplication by 4 (years of study: 2017-2021) of the value corresponds to PSI mean vertical deformation (mm/year) computes the final land subsidence rates for PSI analysis from 2017-2021(mm).
The validation of PSI results is accomplished by leveling measurements. The leveling procedure has a completed accuracy of 0.01cm. For the application of the comparison and validation of two techniques, the nearest benchmark to the according building was selected. This leans to a new assumption; the rate of the vertical displacement of leveling benchmark is accordingly compared with the mean deformation in LoS of the PSs in their buffer zone.
As it is shown in Error! Reference source not found., only three cases exceed 2mm/year; The One, Trilogy, and Olympic Residence, with vertical land deformation values to be -2.393 mm, -5.008mm, and -4. 906mm, respectively. As it can be seen in the description of the case study Error! Reference source not found., those structures are built one next to the other at a distance of one kilometer, with two of them being skyscrapers and one tall building. Buildings like the Dream Tower, Sky Tower, iHome, DTA Tower, and MARR Tower have an average value of vertical deformation below 0.5mm, with the minimum value owned to MARR tower. Due to the medium spatial resolution of Sentinel-1 images, for the estimation of the displacements in the vertical component, the PSs created by the satellite data processing were isolated in a 200m buffer zone around each building for the current study. The building most affected by land subsidence is a skyscraper named Trilogy, with an average value of -5.008 mm, while the less affected is MARR Tower, with a slight uplift to 0.245 mm. Near the region of Trilogy, the One and Olympic Residence buildings are also built, posing the specific area as a high-risk. Similar values of land subsidence were presented in Olympic Residence two-buildings complex with a -4.906 mm over five years. Beyond that, affected buildings with subsidence are the Dream Tower, Limassol Del Mar, The Oval, iHome, DTA Tower, and Arc Ship Tower. Contrarily, in some cases, it can be noticed that uplift happens. Those areas are The ICON, Sky Tower, and MARR tower. Although uplift values rated to almost zero milimeters in five years, the cause of the specific change shall be identified. Concerning being in an area reclaimed from the sea, the deformations may be due to ground settlement taking place in the filler materials under the surface. Positive accumulated ground motion trends can also be shown as relative movements towards the satellite.
The notable importance of studying the behavior of displacements in the environment surrounding artificial objects is that it is advisable to explore the whole area. Suppose land subsidence Limassol Coastal Front could be looked at regionally instead of buildings as a unit. In that case, it can be clearly seen that the area where those buildings are constructed is affected the most. Furthermore, especially in the Molos area of Limassol, in the west part of the studied area, there are places that the land subsidence is equal to -5 mm/ year. Generally analyzing each location in the coastal plain, it can be deduced that except for the new massive development, the erosion of the coastal plain should be studied caused by various factors. Various types of land subsidence may occur in the Limassol Coastal front, namely subsidence due to groundwater extraction, the Mean Sea Level rise, subsidence induced by the load of the constructions, and geotechnique subsidence caused by the alluvium soil. Additionally, it is possible to highlight the presence of new construction sites created in the area of interest.
Finally, due to the recent construction dates, further research in the study area could be directed to applying new combinations of technologies and methodologies. Although leveling is a time-consuming procedure, additional leveling benchmarks can be installed to dense the current altitude network, and every benchmark can correspond to one building. Apart from that, GNSS measurements can be used as the comparison and validation ground-based method for InSAR. A more accurate displacement assessment can show more sufficient results in the following years since the number of satellite images will be extended and Electronic Transponders can be installed, determining Limassol as a perfect city for studying the urban environment through SAR techniques. All the procedures can be automatized, and the possibility of real-time monitoring of land subsidence can be carried out. Through CyCLOPS (INFRASTRUCTURES/1216/0050) program, corner reflectors have been installed all around Cyprus with quite good results.

Conclusions
The combination of using space and ground-based data for detecting and monitoring changes in the coastal zone of Limassol after the exponentially rapid growth, can be characterized as innovative since no case of detecting changes in urban areas in Cyprus has been studied in the past. This first assessment was accomplished using a combination of techniques such as PSI and leveling. Fifty (50) Sentinel-1 images were used in the PSI process, and thousands of PSs were created. Those PSs were isolated in a 200m buffer zone and the average/mean value of the land subsidence in the LoS was computed. For the validation of the PSI results, a leveling procedure was performed, measuring altitude differences in seventeen (17) benchmarks from 2017 to 2021. The detected vertical deformations using the PSI technique have shown land deformation in the majority of the buildings ranging from almost -6mm/year to 4mm/year in a five-year study. These displacements could be related to the coastal incidence and retried caused by the rapid development growth in the area combined with various atmospheric components. Leveling results have shown similar results with a land deformation range of -3.012mm to 0.092 mm. Combining those methodologies, the land deformation in the Limassol Coastal Area was computed between -5mm to 4 mm from 2017 to 2021. The current research has determined the beginning of monitoring the land deformation locally (Limassol) and regionally (Cyprus).
Author Contributions: Conceptualization, CD; methodology, KF and CD; validation, KF, DK, MP, RB, ME and CD; formal analysis, KF, DK and CD; investigation, KF, DK, GM and CD; resources, CD; writing-original draft preparation, KF; writing-review and editing, KF and CD; visualization, KF; supervision, CD; project administration, CD; funding acquisition, CD. All authors have read and agreed to the published version of the manuscript.