Advanced interpretation of land subsidence by validating multi-interferometric SAR data : the case study of the Anthemountas basin ( Northern Greece )

The potential of repeat-pass space borne SAR (Synthetic Aperture Radar) interferometry has been exploited to investigate spatial patterns of land subsidence in the Anthemountas basin, in the northern part of Greece. The PSI (Persistent Scatterer Interferometry) approach, based on the processing of long series of SAR acquisitions, has been applied to forty-two images acquired in 1995–2001 by ERS1/2 satellites. Interferometric results have been analysed at a basin scale as support for land motion mapping and at a local scale for the characterisation of ground motion events affecting the village of Perea in the Thermaikos municipality and the “Macedonia” international airport. PSI results revealed a moderate subsidence phenomenon along the wider coastal zone of Anthemountas basin corresponding to intense groundwater extraction. Highest values, exceeding−20 mm yr−1, were measured in the airport area where the thickest sequence of compressible Quaternary sediments occurs. Intense subsidence has been detected also in the Perea village (maximum deformation of −10 to −15 mm yr−1), where a series of fractures, causing damages to both buildings and infrastructure, occurred in 2005–2006.

These land subsidence phenomena extend gently over large areas (tens to hundreds km 2 ), present low deformation rates and take place for several decades, sometimes without being noticed at the beginning.Localized differential ground deformations can trigger damages to building structures as well as loss of functionality of linear and point infrastructures (pipeline and road network deformations, well-casing failures and protrusion etc.).
One of challenging aspects in dealing with analysis of subsidence phenomena is that their mechanism is not easily detectable.The overlapping of different sources of deformation may complicate the interpretation of the phenomenon.In most of the cases the surface lowering patterns are complex, reflecting a combination of anthropogenic and natural causes.For instance, the natural compaction of unconsolidated fine-grained deposits may be accentuated by human activities, such as over-exploitation of groundwater resource.Moreover, neotectonic activity, spatial variability of geotechnical parameters, temporal variability of water extraction rate F. Raspini: Advanced interpretation of land subsidence may contribute to the different deformation patterns throughout the affected areas.
Detecting, measuring and monitoring subsidence is fundamental for hazard zonation and risk management.Analysis and research on deformation mechanisms are usually limited by the lack of data concerning the rate, spatial extent and temporal evolution of subsidence.Distribution of deformation can be assessed by using conventional, groundbased geodetic instruments, such as levelling (e.g.Phien-wej et al., 2006) and GPS (e.g.Ikehara, 1994) or exploiting new remote-sensing methods (e.g.Galloway and Hoffman, 2007;Galloway and Burbey, 2011).
The current study focuses on the land subsidence phenomena occurring at the Anthemountas basin, located East of Thessaloniki (Fig. 1).
This area draw the attention of the geo-scientists in 2005 when a series of fractures, causing damages to both buildings and roads, occurred at the Perea village, on the southern section of the basin's coastal zone.These fractures were attributed to the overexploitation of the aquifers, although they occurred along an active fault (Anthemountas Fault) (Koumantakis et al., 2008).The application of satellite SAR interferometry (InSAR) for the detection of land motion phenomena revealed that along the coast other large area besides the Perea village are affected by subsidence.The main objectives of the current work are to identify the main causes of the observed ground deformations and to assess the contribution of the remote-sensing data on the study of the phenomena.

Geological background
The broader Thessaloniki area belongs to the NNW-SSE trending alpine Circum Rhodope Belt Thrust System (CRBTS), overthrusting from the East at the Axios geotectonic zone (Peonia subzone).The Circum Rhodope Belt is part of the Inner Hellenic orogen, which has been characterised by repeated SW-directed thrust sheets (Tranos et al., 1999(Tranos et al., , 2003)).The pre-alpine and alpine basement, since the Miocene, is subjected to a brittle extensional deformation, forming NW-SE and E-W directed continental-type basins (Pavlides and Kilias, 1987;Tranos, 1998;Tranos et al., 1999Tranos et al., , 2003)).These basins, among which the Anthemountas basin, are bordered by high-angle normal faults and have been filled with Neogene and Quaternary sediments.
The geological units occurring in the area of the Anthemountas basin are the Mesozoic bedrock formations occupying the bordering mountains, the Neogene deposits outcropping at the hilly areas and the foot of the mountains and the Quaternary deposits occupying the plain area (IGME, 1966(IGME, , 1978a;;Rozos et al., 1998;Anastasiadis et al., 2001).The Mesozoic formations consist of metamorphic (phyllite and gneiss with marble intercalations) and igneous (granite, gabor and peridotite) rocks.The Neogene deposits consist of two sequences, the upper sand and gravel sequence and the lower sandy marls -red clays sequence, outcropping along the borders of the basin.The Quaternary formations occupy the central part of the plain, with increasing thickness towards the coastal area.
According to data coming from geotechnical and geophysical studies as well as from geotechnical and hydrogeological drills (Thanasoulas, 1983;Rozos et al, 1998;Zervopoulou, 2010), the Quaternary formations consist of alternating layers of clastic and fine grained sediments (Fig. 1), generating favourable conditions for the development of an upper phreatic and several successive semi-confined aquifers.Close to the coastline sediments extend to depths of 100 to 140 m, with thickness decreasing gradually to the East (Fig. 1).
As already mentioned above, the wider coastal zone is affected by land subsidence phenomena.Focusing at this particular area, the Quaternary formations can be divided in three horizons (Fig. 1).The top horizon, from the surface to a depth of 30 m, consists of coarse to fine sands with gravel intercalations.The second horizon, with a thickness up to 30 m, is made of impermeable clay to silty clay layers intercalated by fine sand layers.The third horizon, extending down to the Neogene formations consists of coarse to fine sands with a few clay intercalations.Beneath the Quaternary formations, the Neogene sand and gravel sequence occurs.
The Anthemountas basin is bordered by two main normal faults (Fig. 1).The Anthemountas fault (F-An), with an E-W orientation and a length of 32 km, borders the basin to the south, extending from the Galarino village to the Agelochori peninsula.It is the longest active seismic fault close to the city of Thessaloniki, characterised as active because of the clear morphotectonic evidences and the continuous microseismic activity (Zervopoulou, 2010).It constitutes the border line between the rigid Neogene formations and the compressible Quaternary deposits of the plain.Its track is clearly related with the surface ruptures recorded at the Perea village, leading some researchers to the conclusion that they occurred due to the faults activity (Zervopoulou et al., 2007;Zervopoulou 2010).
The second normal fault borders the basin to the north, extending from the Vasilika village to the Mikro Emvolo peninsula.It is a NW-SE oriented fault, dipping to the SW, with a length of approximately 21 km and a possible extension of 8 km into the Thermaikos gulf.According to former studies (Tranos et al., 1999;Karamitrou et. al., 2008;Zervopoulou, 2010) and because of the lack of activity indications this fault is characterised as inactive.
Along the same side of the basin three more active or probably active faults are recorded (Tranos et al., 1999;Karamitrou et. al., 2008;Zervopoulou, 2010).These faults, with a varying length of 5 to 8 km, intersect the above mentioned normal fault presenting an E-W orientation and a south dip direction.

Hydrogeological setting
Evaluating the logs from more than 100 water drills (Nagoulis and Loupasakis, 2001) three aquifers systems were distinguished at the plain of the Anthemountas basin.These systems are: 1.The shallow phreatic aquifers system.The shallow aquifers occupy the upper coarse-grained Quaternary deposits and they are mainly replenished by the stream network percolation.Along the coastal zone they extend down to a maximum depth of 30 m and they are sufficiently separated by the underlying semi-confined aquifers by a thick fine-grained layer (Fig. 1).
2. The semi-confined alternating aquifers.The semiconfined aquifers extend down to depths from 50 to more than 200 m, close to the coastline, occupying the lower Quaternary and the upper sand and gravel sequence of the Neogene deposits.Deeper in the Anthemountas plain (Fig. 1 -Drill Profile 2), where the thickness of the Neogene and Quaternary sediments decreases, these aquifers also extend down to the fractured Mesozoic formations.They are mainly replenished by the bordering mountains' aquifers inflow and also by the shallow aquifer infiltration.
3. The deep confined (artesian) aquifers.Two very deep, isolated aquifers can be distinguished in this category.The first aquifer system extends South of the Anthemountas fault, occupies the lower Neogene deposits and the Mesozoic limestone and contains a subacidic sparkling water rich in calcium and magnesium (known as Souroti Natural Mineral Water).The other deep confined aquifer is the low enthalpy thermal aquifer gushing out along the normal fault bordering the basin to the north and close to the Thermi village.These systems do not seem to be affected by the variations of the ground water level of the two shallower aquifer systems.
The rapid urban growth, the industrial development and the intensification of agricultural activity led to the overexploitation of the aquifers.All the villages along the coastline of Anthemountas basin and many others located on the bordering mountains derive drinking water from local low-lying alluvial aquifers.Some of these urbanized areas constitute the most rapidly developing suburbs of Thessaloniki, experiencing both an increasing urbanization trend and a significant population growth in the last few decades.Between 1991 and 2001 an increase of population of about 108 % and 45 % has been recorded in the municipality of Thermaikos and in the whole basin, respectively (WATERinCORE, 2011).Moreover, several industries, economic activities and infrastructures have developed around the Thessaloniki International Airport.Without a proper management, increase of water demand led to an overexploitation of the groundwater resource.
The overexploitation started affecting the ground water conditions since the 90s.As shown by the isopiezometric curves referring to the conditions of the deep semi-confined aquifers in 1993 (Fig. 2b) and 1998 (Fig. 2a), the ground water level close to the sea was lowered by more than 5 m, during that period.This condition was intensified during the following years, leading to a degradation of the ground water quality and to the occurrence of land subsidence phenomena.The continuation of groundwater overexploitation is expected to spread land subsidence over the entire coastal zone of the basin and to extend deeper at the Anthemountas plain, wherever the thickness of the Quaternary formations is significant, causing great damages both to urban areas and infrastructure.

The InSAR technique applied to subsidence analysis
Several methods apply for measuring, mapping and monitoring spatial extent and temporal evolution of regional and local subsidence.Assessment of ground deformation is historically based on conventional geodetic methods (GPS and leveling network above everything else).
Conventional methods for the production of subsidence maps rely on materialization of a network of geodetic benchmarks, designed to cover the extension of the likely subsiding area.Repeat surveys of benchmarks, referenced to a stable control point, allow the estimation of the deformation extent and rates.These techniques, providing a static picture of the investigated area at each measurements acquisition campaign, are time consuming and resources intensive, since a great deal of time and economic resources are required for timely updates.In addition to conventional geodetic monitoring systems, Earth Observation (EO) techniques have successfully demonstrated in the last decades to be highly valuable in measuring land motion in a wide range of application fields (Tralli et al., 2005), including geohazard-related land motions.

The PSI approach
EO technologies widened their range of applications with the development of Differential Interferometric Synthetic Aperture Radar (DInSAR) techniques.DInSAR constitutes a remote-sensing technique relying on the analysis of phase variations, or interference, between two different radar images (Zebker and Goldstein, 1986;Gabriel et al., 1989;Massonnet and Rabaute, 1993;Massonnet and Feigl, 1998;Rosen et al., 2000), gathered over the same area at different times using the same acquisition mode and properties (beam, orbit, off-nadir angle, etc).The main purpose of DIn-SAR techniques is to retrieve measurements of the ground displacement that occurred between the two different acquisitions.
The application of conventional satellite DInSAR is limited by temporal and geometrical decorrelation (Zebker and Villasenor, 1992).Also, by using individual interferograms the accuracy on the estimated displacement values can be degraded by atmospheric artifacts (Massonnet and Feigl, 1995).
To overcome the main limitations of single-pair interferogram analysis, different approaches (under the "family name" Persistent Scatterer Interferometry), relying on the processing of multi-temporal stacks of satellite SAR images of the same target area, have been developed.Typically, at least 15 images should be available for carrying out a proper PSI analysis.Of course, the larger the number of images the more precise and robust the results.Two classes of multiinterferometric processing techniques, the Permanent Scatterer, PS (e.g.Ferretti et al., 2000;Hooper et al., 2004;van der Kooij et al., 2006) and the Small BAseline Subset (SBAS) technique, are used for processing long series of SAR imagery.
Within the first family, PSInSAR (Ferretti et al., 2000(Ferretti et al., , 2001) ) was the first technique specifically implemented for the processing of multi-temporal radar imagery.Signal analysis of a network of coherent radar targets (persistent point targets), exhibiting high phase stability over the entire observation time period, allows estimating occurred displacement, acquisitions by acquisition.
SBAS technique is an alternative method, originally developed and presented by Berardino et al. (2002) and further implemented by Lanari et al. (2004).This approach relies on the use of a large number of SAR acquisitions distributed in small baseline subsets.Small baseline methods are based on combining a set of unwrapped interferograms, computed in order to minimise perpendicular, temporal and Doppler baseline and to reduce phenomena of phase decorrelation between different SAR acquisitions.The technique allows linking independent SAR acquisition datasets, separated by large baselines, using a combination of differential interferograms.This interferograms are produced by data pairs characterised by a restricted to small orbital separation, leading to the generation of deformation velocity maps and displacement time series.
Mixed approaches have also been proposed (e.g.Hooper, 2008).These hybrid methods identify and exploit the phase of mixed scatterers, selected from Small Baseline interferograms, with the phase of pixels selected from PS interferograms.
The main idea underpinning PSI techniques is to discriminate phase contributions related to displacement from those due to atmosphere, topography and noise through the analysis of the so-called PS (Persistent Scatterers) (Ferretti et al., 2000(Ferretti et al., , 2001;;Werner et al., 2003).Thanks to the Tandem pairs (i.e.ERS-1/2 data with a temporal baseline of 1 day), a conventional InSAR DEM can be reconstructed and subtraction of topographic information from each interferogram can be performed.As an alternative, pre-existing DEM (e.g.such as the Shuttle Radar Topography Mission, SRTM or the ASTER Global DEM) can be used as reference.Once the topographic phase is removed from interferograms, the remaining signal is composed of two contributions: the deformation and the atmosphere-related signal.Spurious atmospheric effects are estimated by using image stacking and filtered out through a statistical analysis of the signals.Atmospheric artifacts are strongly correlated in space within each SAR scene, but are uncorrelated in time.Conversely, target motion usually shows strong correlation in time and can exhibit different degrees of spatial correlation.Once the atmospheric effects are estimated and removed, the remaining contribution maps changes in the satellite-to-target path between the acquisition times of the two SAR scenes.
The phase difference of two SAR scenes, called interferometric phase, is ambiguous, because it is continuously sampled in a discrete wrapped phase (i.e. it can only take discrete values in the interval −π, +π ).Any method seeking for the correct integer number of phase cycles that needs to be added to each phase measurement to retrieve the correct full range (unwrapped) phase can be defined as a phase unwrapping technique.The resulting unambiguous data is called absolute interferometric phase.During phase unwrapping, ambiguity related to the discrete interval sampling of wrapped phase can remain unsolved: low PS density, low temporal sampling and phase aliasing may affect the correct estimation of displacement.
Over urban fabric, where many stable (i.e. with temporally coherent phase) reflectors can be identified, LOS (Line of Sight) deformation rate can be estimated with an accuracy theoretically better than 0.1 mm yr −1 (Colesanti et al., 2003).Unlike the DInSAR approach, PSI analysis is designed to generate time-series of ground deformations for individual elementary reflectors (i.e. the Persistent Scatterers, PS), assuming a linear model of deformation (Ferretti et al., 2001;Werner et al., 2003), or exploiting algorithms estimating the linear and nonlinear components of the displacement (e.g.Mora et al., 2003).The accuracy of the single measurement in correspondence of each SAR acquisition ranges from 1 to 3 mm (Colesanti et al., 2003).Each measurement is referred temporally and spatially to a unique reference image and to a stable reference point.
Nowadays interferometric-based approaches are mature and widely employed by scientific community and practitioners dealing with geo-hazard managements.Since the pioneering studies of Galloway et al. (1998), who introduced In-SAR for the detection of aquifer system compaction effects, it was clear that the availability of SAR imagery archives (like those provided by the European Space Agency, ESA) could support the compilation of subsidence maps, reducing the time for their production and saving economic resources required for timely update.Exhaustive review on the application of remote sensing techniques for mapping subsidence accompanying groundwater overexploitation is given by Galloway and Hoffman (2007) and by Galloway and Burbey (2011).
One of the main priorities to guarantee further advances and continuity of the EO-based geohazards services in the next years and to make them more consolidated and accepted by the user communities, consists in broader use of wide-area processing strategies of satellite SAR imagery.

SAR dataset and processing details
Aiming to map the ground deformations on countries and continents by using the PSI technique, in the framework of the ESA Terrafirma project (http://www.terrafirma.eu.com/), a Wide Area Product (WAP) has been developed by the German Space Agency (DLR) (Adam et al., 2011).
Nine satellite image frames have been processed to produce a PSI ground motion map covering a 65 000 km 2 wide area of Greece, approximately half of the country's territory.This WAP product over Greece has been based on ERS1/2 images.Among these 9 satellite image frames, a dataset of 42 SAR imagery in C-band (5.6 cm wavelength; frequency 5.3 GHz) has been employed for the reconstruction of the history and spatial patterns of land subsidence in the Anthemountas basin.These images have been acquired along descending orbits between 10 April 1995 and 1 January 2001, with a temporal resolution of 35 days, which corresponds to the revisiting time of ERS1/2 satellites.
The reference point of the stack -whose selection during the processing activity is of crucial importance -was chosen in the southeast urban fabric of downtown Thessaloniki.Apart from the phase stability throughout the dataset, the reference point is chosen in area assumed unaffected by ground motions.The reference point is located on rigid over consolidated formations of Upper Miocene-Lower Pliocene age, mainly compact sandstones, conglomerates and marls horizons.These rock to hard soil formations can be safely con- sidered as stable, unable to be affected by local scale deformation processes related to hydrogeological factors.
PSI technique uses large stacks of SAR images to generate differential interferograms with respect to one master image.As shown in Fig. 3 the master image for this dataset has been chosen centrally on 31 December 1996.In the processing chain the master image is selected such as the dispersion of values of the geometrical, temporal and Doppler baseline is as low as possible.The master image is selected to maximize the coherence of the computed interferograms.The stack coherence is larger when the master is selected centrally in time.

Analysis of PSI dataset
The main purpose of the WAP is to implement wide area processing strategies devoted to map countries with the PSI technique.To provide accurately georeferenced, country-scale information on ground motion by using medium resolution SAR sensors (i.e.ERS1/2 with a ground-range and azimuth resolution of 20 m and 5 m, respectively) it is often a challenge.
The validation of the EO measurements accuracy is vitally important.The georeferencing accuracy of the PSI datasets is often difficult to assess, especially when dealing with PSI datasets resulting from ERS processing, due to its low spatial resolution.This fact makes difficult to identify exactly what kind of object is actually acting as "persistent scatterer".The first step for the validation of the PSI datasets was to identify prominent features or isolated objects and to study the PS point distribution around them.Examples of the georeferencing accuracy are given in Fig. 4. As presented, although the theoretical accuracy for the WAP processing is ±25 m in both X and Y direction, much smaller offset values were actually assessed.For instance, the PS points along the Thessaloniki breakwater (Fig. 4a for approximately 6-8 m, while PSs along Thessaloniki pier (Fig. 4b) fit well with the structure of the pier itself.
The potential of PSI techniques is higher in urban areas, because of the wide availability of bright man-made objects.As expected, in the study area, the density of measurement points is higher for the urban area of Thessaloniki (maximum density up to about 300 PS km −2 ), as well as over the urban fabric of the Perea village (up to hundred of point km −2 ).On the contrary over the agricultural and vegetated terrains the density of measurement points is much lower (few point km −2 ), representing a limiting factor of this technique.Although a relatively limited number of reflectors were retrieved, spatial distribution of ground information is sufficient to identify the spatial patterns of deformation.
A frequency histogram has been generated showing the distribution of the 10.061 PS velocity values in the studied region (Fig. 5a).The deformation distribution presents two main components for the time span [1995][1996][1997][1998][1999][2000][2001].First of all, the most common deformation rates are between 50-100 m a.s.l.).This is also consistent with the hypothesis that subsidence largely affects the plain sectors of the Anthemountas basin, as a consequence of intense overexploitation of the aquifers located in the low-lying alluvial basins.The spatial distribution of the deformations the PSs points was plotted on a Visual Earth image (Fig. 6).PSs with LOS velocities between +1.5 mm yr −1 and −1.5 mm yr −1 have been considered as stable points.PSs with LOS velocities < −1.5 mm yr −1 and > +1.5 mm yr −1 have been considered as subsiding and uplifting points, respectively.The colour scale indicates green points as stable.The gradation from yellow to dark red represents increasing subsidence rates.The interpretation of the observed spatial distribution is exhaustively analysed in the following paragraphs referring to specific parts of the basin.

Case studies on specific areas
The city of Thessaloniki and its closest conurbation show very low LOS deformation rates, indicating relatively stable ground conditions during the investigated time period (Fig. 6).Indeed, large parts of the urban fabric are built over Neogene rigid formations, mainly compact sandstones, marls and over consolidated red silty clays, which outcrop extensively in the urban areas.
At the nearby Anthemountas basin the conditions regarding the deformations are different.Unstable areas clearly manifest in the central-lower part of the basin, presenting a mean deformation rate of −10 mm yr −1 .The measured subsidence rates increase towards the coastline (moving from southeast to northwest), where the thickest sequences of compressible Quaternary layers occur.Considering the decreasing drawdown of the semi-confined aquifers piezomet-ric level as well as the decreasing thickness of the loose Quaternary formations towards Vasilika (Fig. 2), the deformations are expected to nullify accordingly.The lack of PSI data concerning the plain area between the subsidence bowl, close to the coastline, and Vasilika village does not allow any further estimation about the extension of the attenuating subsidence phenomenon to the southeast.
Issues of particular interest are presented by the airport area, with several points exceeding −20 mm yr −1 , and by the Perea village in the Thermaikos municipality, with LOS deformation rates of −10 to −15 mm yr −1 .The latter subsidence phenomena were first noticed in [2005][2006], when a series of fractures, causing damages to both buildings and infrastructures, occurred at the Perea village.Since then, no further indications have been recorded at the rest of the area, and the above described extension of the deformations distribution was never recorded or studied by means of any other technique.From the first moment the phenomena occurred at the Perea village were attributed to the overexploitation of the aquifers (Koumantakis et al., 2008) and, as mentioned in Sect.2.2, the dominant process triggering the observed ground deformations along the central part of the basin's coastal zone seems to be the same.

Thessaloniki airport area
Located at the south of Anthemountas River outlet (Fig. 1), about 15 km SE of the city centre of Thessaloniki, the International Airport "Macedonia" is the largest state owned airport in Greece.With 3.9 million passengers during 2010 it is the main airport of Northern Greece and serves the city of Thessaloniki, the surrounding cities of the region and many popular tourist destinations.Considering the geological setting, the sensitive anthropogenic context and the logistic importance of this infrastructure, the analysis of surface deformation is a priority.As highlighted by Ge et al. ( 2009), potential damages of such infrastructures are not only loss of money but also loss of significant landmark in the world.
In order to deeply understand the deformation pattern along infrastructures such as airports, characterised by large linear extent, a monitoring technique coupling wide coverage and high precision is required.During the last years, interest and attention on these problems increased, and analysis of linear infrastructures stability has been successfully carried out with PSI technique (Pigorini et al., 2010).
Subsidence of about −5 to −15 mm yr −1 has been observed in the airport area, with several points exceeding −20 mm yr −1 along the two runway areas.Observed subsidence is without any doubt related to the terrain compaction accompanying the excessive groundwater withdrawal occurred in the Anthemountas plain.Indeed, the airport area is located within a wider scale subsidence bowl which affects a large extension of the basin (Fig. 6).
In order to better investigate the subsidence zone along runways areas, two buffer zones of 500 m have been considered.Two sets of PSs have been extracted within the buffer zones (Fig. 7) to create interpolated subsidence rate maps for both runways areas (about 2400 m in length).
The two sets of scattered points have been used to extend the information (i.e. to interpolate) to areas with lower PS density, assigning values to unmeasured location through the IDW (Inverse Distance Weighted) method (Shepard, 1968).
To create an interpolated surface, the IDW method uses a weighted average of the available neighborhood PS.The weight of input points decreases as the distance between the known point and the interpolation point increases.
Due to the intrinsic characteristics of the satellite radar instruments, the PSI technique can estimate only the sensorto-target displacement, i.e. the projection along the LOS of the real motion.Assuming the occurrence of predominant vertical motion in the Anthemountas plain (as confirmed by Raspini et al., 2012 for groundwater exploitation-related ground subsidence), the vertical deformations (Vv) can be estimated solving the following formula: Vv = LOS/cosθ.Considering the ERS1/2 satellites, the look angle θ is about 23 • .The relation above corresponds to an increment of the LOS of about 8.6 % By interpolating the mean velocity values retrieved by PSI technique within the buffer zone, two subsidence profiles can be generated (Fig. 8a), showing the distribution and magnitude of subsidence along the runway areas of the Thessaloniki airport.Profiles correspond to the longitudinal lines along the centre of the runways.The major subsidence zone could be found in the central part of both profiles (distance range between 1100 and 1800 m), close to the section where the two runways cross each other.The maximum subsidence rate values at this particular section reach the −23 mm yr −1 .
Linear infrastructures, like airports, are mainly affected by differential displacements, and less affected by the absolute values of subsidence.Differential displacements are mainly caused by lateral variations of the thickness of the compressible soil layers (imposed by the stratigraphy or the tectonics

Perea village
PSI results of Fig. 6 reveal that the coastal area of the Perea village (Thermaikos municipality) shows very low LOS deformation rates, ranging between −1.5 and 1.5 mm yr −1 , indicating relatively stable ground conditions since 1995.Nevertheless, in the southern part of the urban area (upper Perea) subsidence can be observed, with maximum LOS deformation rates of about −10 to −15 mm yr −1 (Fig. 9).
The potential of repeat-pass space-borne SAR interferometry can be exploited not only to map the extension of af-fected areas but also to evaluate their deformation history.Displacement time series available for each PS in the area of interest are ideally suited for monitoring temporally continuous geohazard-related ground motions.Example of time series, referring to three Permanent Scatterers located in the Perea's section affected by ground motion, are reported in Fig. 10.
The extraction of groundwater is considered as the main factor causing the land deformation measured, from April 1995 to January 2001, in the village of Perea.The most pernicious consequence of water level drawdown and relative terrain compaction has been the development, in [2005][2006], of a series of fractures affecting the urban fabric.Surface ruptures deformed streets pavement and caused serious damages on buildings, making some of them uninhabitable and forcing their demolishment (Fig. 11).
The abovementioned ground ruptures extend parallel to the coastline, approximately along an E-W direction, dipping towards North.They extend for a length of about one kilometre along the scarp of the Anthemountas fault, defining the boundary between the thick Quaternary alluvial deposits and the Neogene formations, delimiting the subsiding basin.Actually, the manifestation of the ruptures at this particular location, along the scarp of the fault, has been arranged by the abovementioned boundary.
The groundwater withdrawal affecting the area took place over the last twenty years.Analysis of the hydrogeological conditions of the village of Perea is included in Koumantakis et al. (2008), who first identified the overexploitation of confined aquifer as the main cause of the observed ground deformation.As mentioned above, the Quaternary alluvial deposits in the Perea wide area are characterised by both phreatic and semi-confined aquifers (Nagoulis and Loupasakis, 2001).The continuous increase on water demand of the Thermaikos municipally is exclusively covered by the ground water exploitation.A great number of public and private wells are in operation in the wider area.The public wells are deep, reaching down to depths from 120 to 380 m, providing large quantities of good quality drinking water.Despite the public network for water supply and distribution, several private and usually uncontrolled water wells have been drilled without careful management.
Besides the ground water drawdown taking place along the coastal zone of the Anthemountas basin, the narrow area affected by the surface ruptures is located inside the intersecting depression cones of three public and probably several unknown private wells.The locations of the three public wells (G1, G2 and G3) are clearly indicated in Fig. 9 and, as presented, they are only a few tens of metres away from the ruptures.So, it is clear that the narrow area affected by the ruptures is subjected to an extra ground water drawdown, amplifying the subsidence mechanism, active in large parts of the Anthemountas basin.
The intensity of the ground water drawdown can be seen in the graphs of Fig. 12  Tensile ruptures in Perea village are the ultimate effects of terrain compaction related to water level decline.Other visible traces of subsidence, like well-casing protrusion, have been reported northern to surface ruptures.Besides that, no other effects related to the groundwater withdrawal can be clearly detected without analysing PSI data.

Discussion
The main purpose of the WAP approach is to offer a synoptic view of the deformations caused by extensive natural disaster phenomena such as active tectonics, hydro-geological hazards, landslides, throughout the Greek mainland.This huge amount of information, including over a million persistent scatterers over half of Greece's territory, have been exploited scanning wide areas to identify "hotspots", corresponding to sites affected by geo-hazards.Among the numerous sites around Greece verified to suffer from geohazards, the Anthemountas basin case study was selected due to the primitive stage of the occurring land subsidence phenomena (hydrogeological hazard) and furthermore due to the unawareness of the local society and scientists about their actual occurrence and extent.
The analysis of the deformations affecting the Anthemountas basin confirms the importance of early detection and measuring of ground movements related to subtle phenomena like subsidence.Indeed, the mitigation of land subsidence and related effects is very difficult at a late stage.Models explaining and simulating the expected ground subsidence induced by aquifers exploitation are summarised by Galloway and Burbey (2011) and applied by Loupasakis and Rozos (2009) and Herrera et al. (2009).
The compaction of the formations is irreversible and furthermore the progress of the deformations cannot be stopped rapidly by just stopping the overexploitation of the aquifers.The artificial recharge of the aquifers does not seem to give a solution to the problem as the consolidated formations  present reduced transmissivity values delaying the restoration of the ground water pore pressure and, as a result, the deceleration of the deformations.Certainly, the best way to prevent the consequences of the land subsidence is to block their progress before they start causing damages.The integrated water resource management is the best way to cope early stage subsidence.The only problem, as already mentioned above, is that the deformations can take place for several years before being noticed.In that case, this study proved that SAR interferometry and in particular PSI techniques are valuable tools for the early stage detection of the vertical deformations caused by the overexploitation of the aquifers.
Furthermore, main outcomes from the deformation analysis at the Thessaloniki airport give the opportunity to make some considerations about the effectiveness of PSI techniques for mapping and monitoring deformation pattern of linear transportation networks (i.e.railways, highways, main roads), aerial utilities infrastructures (aqueducts, oleoduct and pipe networks for transportation of goods in general) and underground excavations (metro and tunnel constructions).First of all it is worth noticing the density of the two sets of measurement points created for the analysis: (a) the buffer dataset created for the NNW-SSE oriented runway area includes 122 PSs, with a density of 81 PS km −2 and (b) the buffer dataset for the WNW-ESE includes 149 PSs, with a density of 96 PS km −2 .Conventional geodetic networks, despite their robustness and reliability, are not able to provide such a large density of measurement points, as PSI technique does.So, the numerous measurement points, coupled with the mm yr −1 accuracy, enhance the overall understanding and confidence on ground motion occurring in the investigated area, improving the spatial and temporal characterisation of the subsidence.
Full coverage of the runway areas cannot be expected with medium resolution satellites like ERS1/2.Azimuth and range resolutions (along and perpendicular to the flight direction, respectively) of the SAR sensor installed on the ERS1/2 platform are quite poor.The resolution mainly depends on the bandwidth and on the length of radar antenna.The ERS SAR sensor has a bandwidth of 15.6 MHz (wavelength of 5.6 cm), and an antenna length of 10 m.Consequently, the ground range resolution is about 20 m and the azimuth resolution is 5 m.Taking into account the flight path almost parallel to Earth's meridians, east-west oriented linear structures are imaged poor due to the lower range resolution.On the contrary linear structures with a north-south oriented component are expected to be more coherent in the SAR images, due to the finer azimuth resolution.
Due to the high revisiting time (24-35 days) PSI data cannot be used as a real-time monitoring tool neither during construction works, nor for timely updated stability analysis of existing structures.This consideration stands for both previously in orbit satellites (the European ERS1/2 and Envisat) and still operating satellites (the Canadian Radarsat1/2).Nevertheless, since 2008, high resolution imagery acquired by two X-band (3.1 cm of wavelength) radar satellite systems are available: (a) the Italian COSMO-SkyMed SAR constel-lation of satellites, and (b) the German TerraSAR-X satellite.With an increased bandwidth of 300 MHz more coherent targets can be retrieved, providing a higher density of PS points.A better ground resolution (up to 1 m in both azimuth and range direction) and a reduced revisiting time (from 4 to 11 days), are ideally designed for local scale analysis of deformation patterns and for investigating their temporal evolution.Furthermore, given their intrinsic characteristics, this generation of radar sensors allows monitoring faster movements.These enhanced characteristics of the new generation of SAR sensors have improved the capability of PSI technique for land motion mapping (Crosetto et al., 2010).Better spatial resolution, wider applicability to more recent, complex ground movements gave the unique opportunity for a more effective use of SAR data also as monitoring tool in emergency situations (Covello et al., 2010).
Moreover, in the next few years, the launch of new satellite systems (Sentinel 1 and Radarsat Constellation Mission) will further facilitate the research and analysis of land surface phenomena, ensuring the continuity of C-band SAR data (Snoeij et al., 2008) and coupling wide coverage, high precision and short revisiting time.

Conclusion
The potential of SAR interferometry has been exploited through the innovative WAP approach, recently implemented by DLR and aimed at measuring land deformation over large areas.
Interferometric results have been analysed at a basin scale in order to investigate spatial patterns of land motion in the wider Anthemountas plain from 1995 to 2001.The WAP results turned out to be a valuable tool for the characterisation of the land subsidence in the wider plain as, up to now, the only indications of land subsidence phenomena were identified at the village of Perea, affected in [2005][2006] by a series of tensile ground ruptures due to excessive groundwater withdrawal.The PSI data revealed that not only the village of Perea but also large parts of the lower Anthemountas plain, where the International airport of Thessaloniki is located, are affected by land subsidence phenomena.
The detection of subsidence phenomena at an initial stage is extremely important, as further extension of the affected area and damages on settlements and infrastructure can be prevented.

Figure 1 .Fig. 1 .
Figure 1.Location and geological map of the Anthemountas basin.Geology modified from 802 I.G.M.E., 1978b.Three characteristic drill profiles present the vertical distribution of the 803 geological formations from the coastline to Vasilika village.804 805

Figure 3 .Fig. 3 .
Figure 3. Perpendicular baseline for each SAR scene of the PSI dataset, labeled with the 812 temporal baseline.Parameters are referred to the master image which is indicated with the red 813 diamond.814 815

Figure 4 .Fig. 4 .
Figure 4. Examples of the georeferencing accuracy of the WAP dataset: Thessaloniki 817 breakwater (A) and Thessaloniki commercial harbour pier (B).818 819 Fig. 4. Examples of the georeferencing accuracy of the WAP dataset: Thessaloniki breakwater (A) and Thessaloniki commercial harbour pier (B).

Figure 5 .Fig. 5 .
Figure 5. Histogram of LOS deformation velocities distribution in the 821 (A).Comparison between LOS deformation velocities with relative ele 822 of PS points (B).823 824

825Figure 6 .Fig. 6 .
Figure 6.LOS deformation rates in the Anthemountas basin.PSI map is overlaid on Visual 826 Earth imagery.Adopted classification reflects the classes of the histogram of Fig.5.In the 827 inset: Line Of Sight (LOS; red vector) for ERS1/2 in descending mode; green vector, azimuth 828 direction; yellow vector, vertical component of deformation.829 830 831 832 Figure 7. LOS velocity map along the runway areas of "Macedonia" Airport.PSs inside the 834 buffer zones have been extracted and analyzed to evaluate differential settlement.835 836

Fig. 8 .
Fig. 8. Vertical deformation (Vv) of the LOS velocities measured along the runway areas of the Thessaloniki airport (A).Differential displacement along the runways area (B).Origin is set on the coastline.