Abstract
Land subsidence in the coastal zone of the Neo Faliro, Moschato, and Kallithea municipalities, along the Faliro bay, has been recorded since the mid 1960’s. This phenomenon has caused damage to buildings, pavements, and roads. Aiming to identify the main causes of the observed ground deformations, data referring to the geological, geotechnical, and hydrogeological settings of the study area has been evaluated. Subsidence has been quantified by the use of space-born Synthetic Aperture Radar interferometry (InSAR) techniques. SVD (Singular Value Decomposition) and IPTA (Interferometric Point Target Analysis) techniques have been applied for the production of deformation maps, referring to the time period between 2002 and 2010. Furthermore, aiming to extend the study of the phenomenon further to the past, Persistent Scatterer Interferometry (PSI) data for the time period from 1992 to 2001 were also evaluated. Finally, the results of the InSAR analysis have been crosschecked with measurements acquired by a vertical geodetic control network as well as by ground truth data, referring to damage inventory of the site. The current research presents an interesting case study of an urban site affected for a long-lasting period by the activities of a neighboring industrial zone. The development of an extensive depression cone, mainly due to the overexploitation of the aquifers for industrial use, is the main cause of the land subsidence phenomenon, without excluding a component of motion due to the natural compaction of the compressible soil in the area of interest. The complexity of the geological, hydrogeological, and geotechnical conditions and the interaction of the numerous land use activities make this study far more interesting.
1. Introduction
Land subsidence is considered as one of the most frequent geological hazard that usually occurs as a consequence of a number of causes, namely, groundwater overexploitation [1–5], organic soils oxidation [6, 7], and collapse of underground cavities [8–10].
The effects of this geohazard include significant damage to buildings as well as loss of functionality of linear and point infrastructures (pipeline and road network deformations, well-casing failures, and protrusion, etc.) [11, 12].
Detecting, measuring, and monitoring land subsidence is important for urban and infrastructure planning as well as for risk management. Repeat-pass Synthetic Aperture Radar interferometry (InSAR) is a unique tool for large-scale monitoring of surface deformation, at a low cost and great precision [13, 14], and particularly for land subsidence regardless of its cause [15–30].
Numerous studies using InSAR techniques have been conducted in Greece concerning areas affected by land subsidence. Characteristic examples of such areas are Thessaly plain in central Greece [31–34], the wider area close to Kalochori and Sindos villages west of the Thessaloniki plain [35–38], the Anthemountas basin at the east of Thessaloniki [39–42].
The current study focuses on the land subsidence phenomenon occurring at the coastal zone of the N. Faliro, Moschato, and Kallithea municipalities, along the Faliro bay. In this urban area, subsidence has been observed since the mid-60s. Although numerous studies have been conducted by using ground-based techniques, like geodetic surveys [43–47]; up until now, there is not any publication referring to any InSAR analysis results as well as no any publication investigating their causal factors.
The study area is a highly populated urban environment, where the human factor is exposed through this geohazard. The socioeconomic importance of the region is excessive. Several major infrastructure projects have taken place in the study area during the last few decades. Furthermore, the plans for the future development of the region are huge, as the entire sea frond along Faliro bay is considered to be the Greek Riviera. Major residential and commercial constructions combined with marinas are planned to be constructed. The urban expansion and the remaining industrial activities affect the natural evolution of the landscape and the physical processes.
The application of the state of the art technology of Synthetic Aperture Radar (SAR) Interferometry successfully verified the subsidence phenomenon adding new information on its spatial distribution and deformation rates. Land motion mapping data, produced by Singular Value Decomposition (SVD) and Interferometric Point Target Analysis (IPTA), combined with geotechnical characteristics and hydrogeological data are the stepping stone for the development of an efficient management and mitigation of this geohazard.
The objective of the present study is to identify the main causes of the observed ground deformations and to upgrade the existing knowledge for their spatial and temporal patterns by means of advanced DinSAR techniques. To achieve this objective, all parameters affecting the evolution of the land subsidence phenomenon have been evaluated, including geological, geotechnical, hydrogeological, and remote sensing data.
2. The Study Area
The site under investigation is located along the Faliro bay coastal front at the southern part of Athens. In particular, it includes the coastal zone of the Neo Faliro, Moschato, and Kallithea municipalities. As indicated in Figure 1, the area is crossed by two major rivers Kiffisos and Ilissos flowing into the bay.
The study area has experienced both an increasing urbanization trend and significant population growth in the past century. Besides the high-value residential buildings, numerous high impact infrastructures, such as hospitals, stadiums, as well as the National Library and the National Opera House of Greece, are located there. Also, several urban rail lines, sensitive to displacements, are crossing the site. As already mentioned, the future foresees an even more intensive development, supported by the location of the site along the coastline.
Since the 50’s, an extensive industrial zone has been established north of Moschato municipality, the Tavros Industrial zone. It was a mixed industrial park including companies dealing with food processing, leather and fabric processing, footwear and cloth production, and construction materials production; most of them highly water consuming. The industrial activity was severely affected by the financial crisis period, 2009-2018, when the use of the park gradually changed from industrial to commercial, as most of the remaining companies are dealing with logistic activities.
2.1. Geological Setting
Regarding the geological setting of the study area (Figure 2), the top layers are occupied by loose Quaternary deposits; including slightly unconsolidated fine-grained fluvial sediments of Kiffisos and Ilissos rivers and coastal deposits. The deeper strata consist of Pleistocene fine-grained deposits, while the underlay Neogene formations are considered to be the bedrock formation [48, 49].
The loose Quaternary coastal deposits mainly consist of soft clay horizons intercalated with sandy silt to silty sand and silt. These deposits cover a wide zone along the coastline extending to the north for more than 1.5-2.0 Km, reaching down to a maximum depth of approximately 20 m (Figure 2).
The Quaternary fine-grained fluvial deposits consist mainly of sandy clay of low to medium plasticity; intercalate with clayey-silts, clayey-sands, silty sands, and clayey-silty gravels. According to data coming from geotechnical boreholes, their maximum thickness is approximately 30 m, and their lower strata are occupied by the coarser sandy and gravel horizons.
The Pleistocene fine-grained deposits consist of a red-brownish sandy clay (sandy loam) horizon with fine gravels and pisolites, intercalated with loose conglomerates. The Neogene basement is represented by alternate beds of marls, marlstones, marly limestones, dolomitic limestones, and locally intercalations of sandstones, conglomerates, and occasionally lignite. The thickness of the Neogene sediments is approximately 180 m, outcropping at the southwest of the Athens basin (Figure 2).
A number of probable inactive faults of NNW- SSE and NE-SW direction have been identified intersecting the Pleistocene and the Neogene formations of the wider study area [50] (Figure 2).
2.2. Geotechnical Setting
Considering the geotechnical properties of the Quaternary soil formations, it is clear that, according to the oedometer tests, the silty horizons of the coastal deposits are highly compressible presenting compression index (Cc) values ranging from 0.06 to 0.23 (Table 1). Besides that, the Standard Penetration Tests (SPT) reveled that they are loose, as the Nspt values are low, varying from 2 to 14. These horizons certainly contribute to the high compressibility potential of the coastal deposits.
Furthermore, evaluating the overall properties of the coastal deposits, the oedometer tests reviled that the preconsolidation pressure values of several samples have been systematically recorded slightly below the effective geostatic stresses, proving that the loose coastal deposits occur slightly underconsolidated, and as a result susceptible for the manifestation of deformations due to natural compaction. So, a wide zone extending from the coastline up to the southeast boundaries of the industrial zone is susceptible for the manifestation of land subsidence also due to the natural compaction of the relatively unconsolidated formations.
Regarding the Quaternary fluvial deposits, the sandy clay horizons, constituting 70% of the formation, appear to be stiff with low to medium plasticity, normally consolidated and with low to medium compression index (Cc) values, ranging from 0.09 to 0.32 (Table 2). Respectively, the silty sands and clayey-silty gravels are stiff, with even higher compression index (Cc) values. So, the fluvial deposits are less susceptible to the manifestation of subsidence phenomena than the coastal deposits.
The Pleistocene fine-grained deposits include mainly fine-grained horizons up to a percentage of 80%. As indicated by the high NSPT values, they are stiff (Table 3), and as a result, they are not expected to have any particular contribution on the land subsidence mechanism.
2.3. Hydrogeological Setting
Two interacting aquifer systems can be identified, a shallow—unconfined and a deeper semi-confined. The shallow aquifer occupies the top coarse-grain horizons of the Quaternary deposits, while the semi-confined the deeper horizons extending deep down to the Pleistocene deposits.
Furthermore, a confined aquifer system is formed at the semipermeable sandstone and conglomerate horizons of the Neogene formations. Although these aquifers appear to be protected from the seawater intrusion by the impermeable marl horizons, they are not exploited due to their poor recharge rates.
Previous hydrochemical studies [59, 60] have identified the intrusion of seawater along the coastline, during the late 90s. The overpumping of the groundwater, at the industrial area of Tavros, had been cutting off the underground feed of the aquifers located along to the coastline causing their depression and as a result the intrusion of the seawater. The depression cone extending from the Tavros industrial zone down to the coastline can be clearly identified at the piezometric map of July 1997 [59] (Figure 3(a)).
(a)
(b)
The annual water demand in the Tavros industrial area was continuously increasing until the mid-2000s. In the last years, the water demand was reduced and according to piezometric measurements conducted during 2015, the groundwater level gradually recovered. So, by comparing the July 1997 piezometric lines [59] with those based on the data collected during the May 2015 campaign, it is clear that the aquifers have recovered throughout the entire study area (Figures 3(a) and 3(b)). As presented in Figure 3(b), the zero level contour line moved back close to the coastline. The aforementioned recovery is related to the reduction of the water consumption for industrial purposes, triggered by the economic crises of the late 2000s and the resulting changes on the industrial activities.
A clear view of the changes taking place in the last two decades at the groundwater piezometry can be seen in Figure 4, presenting the equal drawdown contour lines produced by comparing the July 1997 and May 2015 piezometric surfaces. As clearly presented, at the vicinity of the Tavros industrial area, the groundwater level has risen up 3 to 7 m since 1997. Along the coastline, the recovery uplifted the groundwater level for 1 to 2 m, moving the zero level piezometry contours back along the coastline, regenerating the barrier protecting the inland from the seawater intrusion.
3. Land Subsidence Historical Background and Ground Truth Data
Since the mid 60’s, the coastal zone of the N. Faliro, Moschato, and Kallithea municipalities has experienced land subsidence phenomena associated mainly with groundwater exploitation [43, 59]. A leveling network has been established since 1989 and geodetic surveys have been carried out in the study area [44]. These measurements indicate significant vertical displacements at the area extending between the outlets of the Kiffisos and Ilissos rivers (Figure 5).
(a)
(b)
(c)
(d)
The maximum displacement rates identified by the vertical geodetic control network can be distinguished into three periods [44]. The first period extends from 1989 to 1995 (Figure 5(a)). At the time, the deformations were affecting a zone extending approximately 500 m north of the Avenue running along the coastline. The maximum total subsidence values were reaching up to -16 mm (approx. -3 mm/yr) [43]. The second period extending from 1995 to 2001 (Figures 5(b) and 5(c)) was a quiet period with maximum total subsidence displacements less than -4 mm [45, 46]. During the period from 2001 to 2008, the accumulated deformation reached a maximum value of -23 mm (approx. -3.5 to -4 mm/yr) [47]. The deformation affected the entire area between the two riverbeds, up to the Tavros industrial zone. Curves with positive deformation values, indicating uplift, are probably related with control points attached to rigid constructions presenting tilting or rotational movements, uplifting them.
At this point, it should be noted that at the early 2000’s, due to the pre-Olympic game financial growth of the country, the activities at the Tavros industrial area were expanding. The construction activities, resulting to extensive draining works, were also growing at the entire coastal zone. Furthermore, Kiffisos River was fully channelized cutting off the outflows to the aquifers. So, although no groundwater level measurements exist for that time period, it is fully acceptable that the groundwater table was susceptible to drawdown, justifying the amplified subsiding deformations observed between 2001 and 2008.
Following the spatial distribution of the deformations, numerous surface ruptures and structural failures can be identified at the area between the riverbeds of Kiffisos and Ilissos (Figure 6).
4. Multitemporal InSAR Methods and Processed Data
In the present study, twenty Single Look Complex (SLC) Advanced Synthetic Aperture Radar (ASAR) C-band images of ESA (European Space Agency) Envisat satellite have been processed. The images refer to the time period between 2002 and 2010, and they are acquired along the descending orbit. Due to the poor dataset of available scenes covering a period of 8 years, two multitemporal InSAR techniques have been applied namely (a) SBAS-type processing and (b) the Persistent Scatterer Interferometry using the IPTA (Interferometric Point Target Analysis) algorithm for the production of time series deformation maps. The first and common step after the introduction of the images in the system is a Digital Elevation Model- (DEM-) assisted coregistration procedure, implemented at the SLC images. The SBAS-type method is characterized by a first multireference interferometric processing to a single reference image by applying the SVD algorithm [61], obtaining the time series surface deformation. The second one is the classical time series interferometric processing. Both techniques at the last step of processing transform the deformation phase into line-of-sight (LOS) displacement and geocoded maps (Figure 7).
Aiming to extend the study of the phenomenon further to the past, already processed Persistent Scatterer Interferometry (PSI) data covering the time period from 1992 to 2001 were also evaluated. The PSI analysis was conducted within the framework of the ESA GMES Terrafirma project (GMES—Global Monitoring for Environment and Security) funded by the European Space Agency (ESA) [40]. During this project, 9 satellite image frames were processed using a special semiautomated processor to produce a PSI ground motion map covering a 65,000 Km2 area of Greece. This Wide Area Product (WAP) over Greece is based on strip map ERS 1/2 images, obtained from the ESA.
The SBAS-type processing using the SVD algorithm method allows connecting disjoint subsets of interferometric pairs separated by large baselines [61, 62]. SVD relies on the minimum-norm criterion of the deformation rate to acquire time series [63, 64]. The benefit of this technique is that the nonlinear deformation can be estimated without the need of a priori assumption of a deformation model or past knowledge and as any other multi-temporal approach, the mitigation of decorrelation phenomena and topographic inaccuracies are limited [62]. More precisely, this approach is characterized by a first multireference interferometric processing. So, a large number of interferometric pairs are produced to a single reference image interferometric processing, obtaining the time series surface deformation by applying the SVD algorithm [59]. In this way, any potential deformation values that may occur between two dates are added to the subsequent one, for each image acquisition until the last one. Therefore, creating a cumulative deformation history for each pixel that exists in the available time interval and thus enabling the creation of time-series graphs.
At the IPTA technique, based on the coregistered SLC scenes, two lists of candidate point scatterers were initially generated using different selection approaches. Next, two point lists were generated and combined into a single one, which was used in additional processing after removing duplicate records. The single reference method has been used in order to form multiple pairs of interferograms. The criteria by which the reference scene was selected are the following: (i) forming interferometric pairs with the minimum baseline backprojection (Bp), (ii) acquisition date near the central date of the time period for which there are available SAR acquisitions, and (iii) reference scene to present low atmospheric distortions.
Consequently, the initial point differential interferograms were generated using the coregistered SLCs, the DEM heights, and the predefined Persistent Scatterers (PS) list. The topographic component was subtracted using the DEM. The generated point differential interferograms were further analyzed in the temporal and spatial domains in order to obtain information on the atmospheric phase term, deformation phase, and baseline errors.
The following steps involved least-square regression using the differential phases in an iteration process to estimate scattering heights and the deformation rate relative to the selected reference area.
Based on the regression analysis, the quality of the PS candidates was further evaluated through the estimated standard deviation of the phase difference.
Furthermore, since residual phases contain the atmosphere, nonlinear deformation components, and other error terms, they were further processed in order to compensate for these effects. Thus, residual phases were filtered assuming that atmospheric contribution is high-pass in time and low-pass in space. Further iteration was done in order to consider these corrections in the final regression model. It is worth mentioning that phase unwrapping is performed exclusively in the temporal domain, avoiding issues that could arise from error propagation when the spatial unwrapping of individual layers is considered. The generated results consist of height corrections, linear deformation rates, atmospheric phase, temporal coherence levels, and nonlinear deformation histories for each PS. Finally, the deformation phases were transformed into displacements and geocoded to a selected map.
5. Results and Discussion
5.1. Available PSI Land Motion Mapping Data (1992-2001)
The available land motion mapping data, produced by PSI analysis, revealed vertical displacements at the wider study area during the time period from 1992 to 2001 (Figure 8(a)).
(a)
(b)
As presented in Figure 8(a), the Persistent Scatterers (PS) indicating subsiding movements (yellow to red dots) are concentrated along the coastal zone, with a higher density between the riverbeds of Kiffisos and Ilissos. Also, a stack of PSs can be identified at the northern part of the Tavros Industrial zone. The maximum displacement rates at the coastal zone reach up to values of -1.5 to -4 mm/yr, while even larger displacement values, up to -6.9 mm/yr, can be identified at the industrial zone. Also, the analysis of the time series clearly shows a steady downward trend in the time intervals covered by ERS 1/2 dataset (Figure 8(b)).
The spatial distribution of the deformations fits with those identified by the geodetic control network, although the maximum displacement values appear to be larger at the InSAR dataset.
5.2. IPTA and SVD Land Motion Mapping Data (2002-2010)
Both DinSAR techniques applied at the Envisat data reveal an identical pattern of deformations for the time period extending from 2002 to 2010. As presented in Figures 9(a) and 10., both datasets indicate that the maximum displacements can be identified at the area included between the riverbeds of Kiffisos and Ilissos, at the south of the Tavros industrial zone. Nevertheless, slight differences can be identified at the maximum displacement rates. The Interferometric Point Target Analysis (IPTA) indicated deformation rates of -3 to -5 mm/yr, along the satellite’s LoS (Figure 9(a)), whereas the Singular Value Decomposition (SVD) interferometric technique indicates lower LoS deformation rates, ranging between -1.5 and -3 mm/yr (Figure 10).
(a)
(b)
Displacement time series are also available for each PS target, ideally suited for monitoring the temporal evolution of displacement for the investigated period of time. Characteristic graphs referring to two PSs affected by the ground motion are reported in Figure 9(b).
5.3. Correlation of Space Borne, Geodetic, Hydrogeological, and Geotechnical Data
The space-born interferometric data indicate clearly that the area affected the most by the land subsidence is the one extending between the riverbeds of Kiffisos and Ilissos. This conclusion fits perfectly with the results coming from the vertical geodetic control network. Evaluating the displacement rate values, the PSI dataset coming from the processing of the ERS 1/2 acquisitions indicate slower displacement rates during 1992 to 2001 than those identified by the Envisat data for the period 2002-2010. So, as in the case of the geodetic measurements, the deformations were accelerated during the second decade. Nevertheless, the displacement rate values of the first decade (1992-2001) do not fit perfectly with the geodetic data. Particularly, the geodetic measurements indicate displacements with velocities lower than -1 mm/yr while the PSI dataset velocities from -1.5 up to -4 mm/yr. On the contrary, the IPTA dataset of the second decade (2002-2010) fit perfectly with the geodetic data as both indicate displacement rates from -3.5 to -5 mm/yr.
Comparing the July 1997 groundwater contour lines (Figure 3(a)) with the distribution of the displacements acquired by both space-born and geodetic data, it is clear that the maximum displacements are located at the area affected by the depression cone, extending from the industrial zone down to the coastline. Unluckily, there are not any groundwater level measurements referring to the period between 2002 and 2010, in order to correlate them with the acceleration of the displacements. Nevertheless, as already mentioned, the pre-Olympic games financial growth of the country, during the early 2000’s, expanded the activities of the Tavros industrial zone, as well as the construction activities along the coastal front, tensioning the depression of the aquifers. The full channelization of the Kiffisos River outlet, completed before the Olympic games of 2004, supported the evolution of the depression cone by cutting off the outflows of the river’s water to the aquifers. The effect of the rebound of the aquifer on the deceleration of the displacements cannot be witnessed as space-born or geodetic data are not available for this time period. As mentioned, the rebound of the aquifer, identified at the contour lines of 2015 (Figure 3(b)), can be attributed to the financial crises pausing all financial activities (industrial and construction) at the site.
Regarding the geotechnical properties of the geological formations, the loose silty horizons of the Quaternary coastal deposits present low compression index (Cc) values making them highly susceptible to the manifestation of land subsidence due to consolidation. The deformations attributed to the remaining natural compaction of those formations combined with the occasional consolidation events caused by the newly build constructions triggers additional deformations at the areas occupied by the coastal deposits.
6. Conclusions
The study area is a high-interest urban site along the Faliro bay coastal zone of the city of Athens, including districts of the Neo Faliro, Moschato, and Kallithea municipalities. The site, planed to be the Greek Riviera, is subjected to continuous deformations due to the land subsidence phenomena, occurring since the mid 60’s.
The present study provides a methodological approach for the identification of the subsidence mechanism by combining the geotechnical characteristics with the hydrogeological data and by verifying the displacements with the use of space-born SAR interferometry (InSAR) techniques and ground truth geodetic measurements.
As proved by the groundwater leveling data, the industrial zone extending at the north of the study area, due to the overexploitation of the aquifers, causes the development of an extensive depression cone, extending down to the coastline. Furthermore, the costal deposits covering a big part of the study area are proved to be compressible due to the low compression index (Cc) values identified in some of their silty horizons. So, the development of the depression cone is expected to be the main cause of the land subsidence affecting the entire area extending between the riverbeds of Kiffisos and Ilissos. Also, an additional deformation component is expected at the areas occupied by the coastal deposits as they have been proved to be slightly unconsolidated.
The manifestation of the expected deformations has been verified with excellent accuracy by the use of space-born SAR interferometry (InSAR) techniques. Furthermore, the results of the InSAR techniques have been crosschecked with measurements acquired by a vertical geodetic control network as well as by ground truth data, referring to damage traced on building and rigid constructions.
The current research demonstrates an interesting case study of an urban site affected for a long-lasting time period by the activities of a neighboring industrial zone. The complexity of the geological, hydrogeological, and geotechnical conditions and the interaction of the numerous land use activities make this study far more interesting.
Furthermore, this study demonstrates the potential of space-borne InSAR techniques as a suitable tool for the validation of the subsidence movements; a cost-efficiency method, complimentary to ground-based measurements.
Eventually, the current study indicates that by applying the proper techniques and methodology the early detection of surface deformations can be established, allowing taking measures before severe subsidence occurs and therefore allows for timely protection of the affected areas.
Data Availability
The ERS 1/2 data as well as the Envisat raw data have been downloaded form: https://vertex.daac.asf.alaska.edu. The hydrogeological and geotechnical data are coming from national archives and they are only assessable by visiting the libraries.
Conflicts of Interest
The authors declare no conflict of interest regarding this publication.
Acknowledgments
The Terrafirma Extension project has funded the SAR imagery processing of the ERS 1/2 acquisitions covering the time period from 1992 to 2001. The project is one of ten services being supported by the Global Monitoring for Environment and Security (GMES) Service Element Programme, promoted and financed by the European Space Agency (ESA). The authors gratefully acknowledge the German Aerospace Centre (DLR) for having processed the SAR data of the ERS 1/2 acquisitions.