Ground deformation detection of the greater area of Thessaloniki (Northern Greece) using radar interferometry techniques

. In the present study SAR interferometric techniques (stacking of conventional interferograms and Permanent Scatterers), using images from satellites ERS-1 and 2, have been applied to the region of Thessaloniki (north-ern Greece). The period covered by the images is 1992– 2000. Both techniques gave good quantitative and qualitative results. The interferometric products were used to study ground surface deformation phenomena that could be related to the local tectonic context, the exploitation of underground water and sediments compaction. detected and assessed.


Introduction
Deformation measurements for greater area of Thessaloniki were carried out using methods of Stacked Differential Interferometry and Permanent Scatterers SAR Interferometry (PSInSAR). Processing of Synthetic Aperture Radar (SAR) interferometric pairs from ERS1 and ERS2 satellites and their correlation with the potential sources of deformation are presented in this article for the period between 1992 and 2000. The observed regular ground subsidence being produced by different causes (such as abstraction of earth material by mining or tunneling or fluid pumping of groundwater, oil or gas) is not the same type of hazard as sudden and catastrophic natural events like earthquakes. However, high vulnerability to the effects of subsidence is exhibited by cities. In these regions and especially in those that are highly populated and heavily industrialized, such as greater Thessaloniki, where important infrastructure is present or are under development, the exposure to the hazard of subsidence increases.
The first section of this article describes the context of ground deformation and the theoretical background of advanced interferometric techniques. The following sections deal with data processing, detection of deformation and the interpretation of the interferometric results.

Geotectonic and geodynamic accounts
The City of Thessaloniki is established as the second most important urban centre of Greece. It is located in Central Macedonia and is situated in the inner part of the Thermaikos Gulf ( Fig. 1). It has an extended industrial zone in its suburbs and a major international port that constitutes the centre of merchant shipping for the Balkan countries. The main morphological units in the area are the N-S trending Axios river basin, the smaller with the same trend Gallikos river basin and the Mygdonian basin in the northern part with a NW-SE to E-W trend. The boundaries of the deltaic deposits of the Axios and Gallikos rivers cannot clearly be distinguished and therefore the term "deltaic complex" is generally used hereafter. The deltaic plains, with a bird-foot shape, at the river mouths indicate the dominance of fluvial over marine process in forming the most active part of the Thermaikos shoreline. It has been estimated on the basis of hydrological charts that the area of the deltaic plain of the Axios river grew seawards by 175 km 2 between 1850 and 1987 (Poulos et al., 1994).
The broader area is geologically (Fig. 2) composed of Neogene and Quaternary sediments. Only the hilly and mountainous areas are formed from pre-Alpine and Alpine formations. The pre-Alpine rocks are mostly two-mica schist and gneiss, biotitic gneiss and amphibolites. Deepsea meta-sediments, meta-volcanoclastic rocks, phylites and quarzites and rare exposures of ophiolites compose the Alpine rocks. Over the pre-Alpine and Alpine formations Neogene and Quaternary deposits have unconformably been laid down. These Neogene deposits of the Late Miocene-Pliocene, which are mainly filling the basins of Axios, Gallikos and Mygdonian, are extensive at the north and east of Thessaloniki, but in the Mygdonian Graben are met to a lesser extent. The Thessaloniki area consists of consolidated red beds with abundant intercalations of sands and gravels. East of Thessaloniki, some red beds are overlain by claylymarly sediments of Pontian age. In the Mygdonian Graben the Neogene deposits are localized along the mountain frontiers and consist at their base of conglomerate beds overlain conformably by sandstones, and continue with rhythmic alternations of silts, sands and red beds of Pleistocene age. The Quaternary deposits presented in the area comprise of undivided deposits consisting of scree and fan deposits, Pleistocene deposits, concerning undivided deposits of brown-red sediments and Holocene deposits that cover a large part over the area under investigation.
The development of basins started during the Neogene. The existent Miocene planation surfaces were broken as a result of faulting. Some of the blocks have been uplifted by about 300-400 m during the Neogene and Quaternary forming horsts, while other blocks have been subsided forming depressions that subsequently were filled by sediments. The noticed subsidence in the Gallikos and Axios basins was Nat. Hazards Earth Syst. Sci., 8, 779-788, 2008 www.nat-hazards-earth-syst-sci.net/8/779/2008/ about 400 to 600 m during the Quaternary. The geodynamic regime of the broader area is characterized by continuous extensional deformation associated with active normal faulting, trending mainly E-W, WNW-ESE and NE-SW. A number of inactive faults of NW-SE direction in the south-eastern part of Thessaloniki, and of NE-SW direction in the northeastern part are also present (Fig. 3). An almost continuous seismic activity is taking place in the area of the Mygdonian Basin. On 20 June 1978, a destructive earthquake of magnitude Ms=6.4 took place in this basin (Papazachos and Papazachou, 1997). In the broader area of Thessaloniki, though, including the Axios Basin smaller magnitude earthquakes with focal depths ranging between 5 and 15 km and magnitudes ranging from 1.5 to 4, are often recorded. According to the above concepts the expected deformations should be attributed both to active tectonic and to sediment compaction.

Conventional SAR interferometry (InSAR)
Repeat-pass space-borne radar interferometry (Gabriel et al., 1989;Hanssen, 2001;Massonnet and Feigl, 1998) is based on the use of phase information that is obtained by Synthetic Aperture Radar (SAR) instruments onboard satellites. The phase-difference image called the "interferogram" depends on both the topography and possible deformations of the ground, which may occur between the acquisitions of the two satellite images. Once the topography contribution has been modeled and removed from an interferogram, the obtained differential interferogram gives access to the ground deformations that have occurred between the two acquisitions. One fringe of an interferogram corresponds to a deformation of half a wavelength (28 mm in the case of ERS satellites) in the Line Of Sight of the satellite (LOS).The application of conventional Differential interferometry (DIn-SAR) regardless of the applied method (three passes or two passes plus DEM) is limited by the following constraints: (i) loss of coherence, which occurs when the physical/geometrical nature (e.g. vegetation, water or ploughed field) of the ground changes and the stability of the phase signal is lost (Zebker and Villasenor, 1992); (ii) atmospheric artifacts: the fluctuations of the atmospheric (tropospheric and ionospheric) layers between the two acquisitions induce subtle phase variations, which can be misinterpreted in terms of ground deformation signatures (Massonnet and Feigl, 1995); (iii) uncompensated topography; (iv) instrumental limitations, orbital repeat-cycle maximum deformation gradient detectable, pixel size, archive data availability, 1-D measurement in the line of sight of the sensor are other factors to be considered.

Stacking interferometry
Stacking of differential interferograms aims to combine the information from several differential interferograms in order to extract common information. The most basic procedure is to compute linear combinations (generally sums or averages) of interferograms. More elaborated methods have been proposed. For example, Sandwell and Sichoix (2000), stacked interferograms by combining phase gradients. For this study, we will sum multiple differential interferograms to obtain single interferograms. Interferogram stacking is useful in overcoming the following two shortcomings of conventional InSAR.
(i) Low coherency over long temporal separations. If reasonable coherency levels can only be obtained over short time periods (for example, in the case of rural settings in temperate climates) then several shorttime-period temporally-contiguous interferograms can be summed (subject to data availability) to produce a pseudo-interferogram over a longer period. This enables low magnitude displacements to be monitored over longer periods, where no single coherent interferogram exists.
(ii) Atmospheric influence. When multiple differential interferograms exist that brackets an instantaneous event (such as an earthquake or other sudden ground displacement) they can be summed to increase the (displacement) signal to (atmospheric) noise ratio. This is possible because the displacement signal is constant in each interferogram, whereas the atmospheric signal varies randomly.
The computational effort required to stack interferograms is limited (once the interferograms themselves have been produced). However, such a strategy should only be applied in the selected circumstances described above. Also an understanding of the temporal evolution of the deformation is desirable for determining if interferogram stacking is relevant. This generally means that a previous interpretation of the interferogram series is required before stacking is considered (Raucoules et al., 2003b;Strozzi et al., 2001).

Permanents Scatterers (PS) interferometry
The PSInSAR technique (a detailed description can be found in Ferretti et al., 2000 andColesanti et al., 2003) is an advanced processing tool allowing the joint exploitation of a series of interferometric SAR images all referred to a unique master acquisition.
The processing is carried out on a network of points (PS) of the SAR images having stable radiometric characteristics. The PS selection is based on statistics on the SAR amplitude values. The processing is thereafter performed on the selection of points providing results.
The PS approach to overcome limitations due to atmospheric conditions is based on a few basic observations. Atmospheric effects show a strong spatial correlation within every single SAR acquisition, but they have low temporal correlation . Conversely, target motion is usually strongly correlated in time and can exhibit different degrees of spatial correlation depending on the phenomenon at hand (e.g. subsidence due to water pumping, fault displacements, localized sliding areas and collapsing buildings). Based on these observations, atmospheric effects can be estimated and removed by combining data from long time series of SAR images, such as those available in the ESA ERS archive, which includes data since late 1991. As in all differential interferometry applications, results are not absolute both in time and space. Deformation data are referred to the master image (in time) and results are computed with respect to a reference point of known elevation and motion (in space). Despite this remark and the fact that it provides just one component of the deformation, PS is a sort of natural geodetic network allowing the analyses of surface deformation phenomena over hundreds or thousands of km 2 . It provides a complement to traditional monitoring methods like GPS and optical levelling and an alternative for sites that were not instrumented before an event (Colesanti et al., 2001). Finally, we have to note that the statistical approach (fitting of the deformation and height models to the data, Ferretti et al., 2001) used in the technique for separating the different components (atmosphere, deformation and height) is reliable only if the amount of data is sufficient. A set of at least 30 images is required and the precision of the final assessment increases with the number available . In addition, the standard PSInSAR approach used in this study assumes that the deformation is linear in time (constant velocity).
The complementary use of stacked conventional DInSAR and standard PSInSAR could provide information for a better identification of the deformation. Standard PSInSAR measurement is more precise but can be limited where the deformation is too fast or non-linear (Crosetto et al., 2007) whereas conventional DInSAR (although less precise) based on well selected SAR images could provide an assessment of the highest deformations values or when the PSInSAR is lost due to characteristics (value/linearity) of the movement. This will be observed in this study ( Fig. 4 and 5).

The SAR Data Used
A set of 47 SAR images from the satellites ERS1 & ERS2 was obtained from the European Space Agency (ESA), covering the period 1992-2000 (Table 1). These radar scenes were used for the production not only of conventional differential interferograms (DInSAR), but also in the analysis of the advanced techniques of Stacked and PS Interferometry.
Conventional interferograms were computed (with a 2×5 multi-looking) using the Gamma software (Wegmuller et al., 1998) by re-sampling each SAR scene to a common geometric reference from an image (in our case orbit 12 176) whose perpendicular component of the orbital baseline is close to the mean of the set in order to avoid excessive distortions between images. Using such a common geometry simplifies and allows the automation of most of the processing. We have to note that this geometric reference image will not be the interferogram reference (master). To produce a given interferogram from two images, their phases are subtracted to produce the interferogram. The flat earth contribution has been removed using a Fourier transform to determine the corresponding fringe rate. This step provides an estimation of the perpendicular baseline, which is then used to simulate and remove the contribution of the topography using a Digital Elevation Model having a resolution of 50×50 m. The final resolution of the geocoded products is 25 m×25 m. All the interferometric pairs, with perpendicular baselines smaller than 200m, were automatically generated. By limiting the baseline, the most incoherent interferograms were rejected. Subsequently, the most relevant interferograms were visually selected by rejecting the pairs affected by phase noise and atmospheric effects. The following criteria were finally considered for this selection.  4. Interferometric stacking image produced after processing of three relevant interferograms (dates 1996.09.17/1997.06.24; 1999.10.11/2000.06.13; 1999.11.16/2000.10.31. The green rectangle shows a coherent area of the interferogram (West of Kalochori) with potential deformation larger than 50 mm.yr −1 .
1. The noise on the interferogam (in fact the possibility of a clear separation of the fringes and not only a threshold on the coherence).
2. The fact that the presumed deformation signatures appeared on several independent interferograms (rejecting the possibility that they could be atmospheric artifacts).
Conversely interferograms with an excessive number of atmospheric effects (observed artifacts) were rejected.
3. If non-linear deformation signatures were exhibited (e.g. if a deformation signature was identified in a short time-interval with respect to the full studied period or if the sign of deformation changes during the studied period), these interferograms would be particularly interesting to interpret. This is because a standard PSInSAR product is not able to detect non-linear deformations. Nevertheless, in this study, strongly non-linear deformations were not observed.
Such a visual selection on large interferogram sets was used in previous studies and provided reliable results (Le Mouelic et al., 2005;Parcharidis et al., 2006;Raucoules et al., 2007). In a first selection 100 interferograms were chosen (Table 1). After this interpretation we considered the pairs of: 1996.09.17/1997.06.24; 1999.10.11/2000.06.13 and 1999.11.16/2000.10.31 as the most significant in terms of the description of the deformation and therefore these were chosen to be stacked.
For the application of the PSInSAR technique, the selected master image for our data set corresponds to the orbit 7376 (i.e. image acquired on 1996.09.17) that reduces the average baseline magnitude. All other scenes are coregistered and resampled to the reference image. In the frame of this study a "standard PS analysis" velocity map was produced (Fig. 5).

Correlation between the observed deformation and potential sources of deformation
Interesting generic observations may be made by examining the two final inerferometric images. The interferometric stacking image shows clear fringe patterns of deformation in the area of Sindos-Kalochori, Oreokastro and less clearly in Langadhas Basin, and the airport area. A large part of the image is covered by noise, specifically in the Axios and Mygdonian Basin areas. The PS interferometric image shows a very high-density concentration of PSs in the area of the city of Thessaloniki and in urban centers in greater Thessaloniki. Outside the urban centers, as in the agricultural fields and the mountains, Nat. Hazards Earth Syst. Sci., 8, 779-788, 2008 www.nat-hazards-earth-syst-sci.net/8/779/2008/ there is a very low density of PSs, which is inadequate for the sampling of small-scale deformation phenomena, representing thus a limiting factor of this technique. Deformation areas are noted in the Sindos-Kalochori area, Oreokastro, Langadhas area and in the area of Thessaloniki Airport. An issue of particular interest is that there are no observed deformation patterns in both images within the metropolitan area of Thessaloniki. A detailed description of deformation for the areas affected by subsidence phenomena and their potential correlation to deformation sources is outline in the following section.
We have to remark that the interferometric products provide the Line Of Sight component of the deformation. Because of the rather small (23 • ) incidence angle of ERS, the measurement is much more sensitive to vertical deformation than to horizontal deformation. The comparison with levelling, carried out in this study, would be fully relevant if no horizontal deformation occurs. In the following sections (and based on the previous knowledge of the deformation on this area), we assume that the deformation is mainly vertical.

Sindos-Kalochori
This area (Fig. 6) constitutes the industrial zone of Thessaloniki (especially Sindos) while the coastal zone (Kalochori) is characterized as an agricultural zone. The observed subsidence rate in Kalochori (more than 40 mm.yr −1 ) is extended to another subsidence bowl to the west of Sindos. Even though the coherence is poor in this area due to its agricultural nature, some sparse PS's in that zone confirm that the area is subsiding. However, we can observe that the maximum deformation assessed by PSInSAR is lower (less than 30 mm.yr −1 ) than assessed by conventional DInSAR. This is probably due to a loss of information on the PSs undergoing higher deformations. A more thorough examination on the stacked interferograms shows a small coherent area at the border of the deformation with a subsidence rate of 50 mm.yr −1 . This zone is therefore characterized by subsidence of 50 mm.yr −1 over an area of more than 10 km in diameter.
According to Andronopoulos et al. (1991) the Quaternary deposits in the area are classified in three horizons: sandy, silty and black silty clays. Since 1955, overpumping of the ground water in the Kalochori area to supply the city of Thessaloniki has enhanced the subsidence phenomena in the area. The pumping was drastically reduced in the early 1980s and only small amounts of water are extracted, mainly from private wells. The observed subsidence phenomena should, therefore, rather be attributed to the intense and extended water pumping from the 1960s (Andronopoulos et al., 1991). Subsidence is estimated at 2.5 m for the period 1955-1985(IGME, 1989. This caused a gradual and significant fall in the level of the water table. This fall caused the drainage of saturated sediments and their consolidation to be visible as subsidence of the ground surface. Therefore, the area between Kalochori and Sindos seems to have been affected by uniform subsidence reaching up to 3 m from 1955 to 1980 (Hatzinakos et al., 1990). Between 1980Between to 1985Between and 1985Between to 1999, rates of relative height changes between 8 to 10 cm.yr −1 respectively, are observed based on triangulation re-measurement results (Stiros, 2001).
It is, hence, concluded that the lack of spatial and time correlation between fluctuations of piezometric levels, topographic changes and pumping indicates that the observed subsidence should be regarded as the cumulative effect of several factors, such as consolidation of near-surface sediments due to the decline of the piezometric level and the partial abandonment of the delta, oxidation of peat soils in the vadose area, syn-sedimentary deformation, as well as loading-induced consolidation of deeper sediments. The general local (over a few kilometres) subsidence has increased near shore water depths and thus provoked an increase in wave activity. Due to this the sea barriers that protect the deltaic plain were destroyed and catastrophic floods have occurred several times.  (Chatzipetros et al., 2005;Gounaris et al., 2007;Martinod et al., 1997;Stiros, 2001;Stiros and Drakos, 2000).

Mygdonian Graben
A subsidence rate of 20 mm.yr −1 is observed in the western part of the Mygdonian Graben and specifically in the Langadhas area (Fig. 7). The former is filled by Holocenic deposits such as pebbly gravels, pebbly sands, sandy clays and lake deposits. The basin is demarked from the Assiros-Analipsi fault zone to the north and Lagina-Agios Vassilios fault zone to the south (Langadhas sheet map). Both faulting zones are characterized as active but not as seismic because there is no information connecting the two faults with historical seismic events. Between 1958 and 1978, an observed subsidence of about 25 cm (annual average of 12 mm) was observed along a distance of 30 km in the basin from highprecision spirit-levelling results (Stiros and Drakos, 2000). Papazachos et al. (2001) studied active crustal deformation for the most active sections of the Mygdonia Basin using the combined stress pattern and corresponding moment-rate tensors. The results show a N-S extension for the central part at an average rate of 3.5 mm.yr −1 , which is consistent with available differential GPS re-measurements (Martinod et al., 1997). According to the former study, comparison of results between 1979 and 1994 shows that this part of the Mygdonian Graben was subjected to a N-S horizontal extension of about 80 mm. This horizontal extension was interpreted either as a long-term post-seismic relaxation processes within the graben or as a continuous aseismic slip motion. Based on paleoseismological trenching in different sites within the basin Chatzipetros et al. (2005) evidence aseismic creeping with a maximum slip-rate of 2.5 mm.yr and a minimum of 0.7 mm.yr. According these authors aseismic creeping plays a very important role in this area.
According to estimates by the Directorate of Water resources and Agricultural Engineering of the Thessaloniki Prefecture in the greater area of Langadhas there exist 2.100 legal and illegal water boreholes. In relation with the drought period of the last two of decades (especially during the first half of the 1990s) a considerable decline in the piezometric level of about 80 m has occurred mainly during the summer seasons.
Based on the above considerations the observed deformations could be attributed either to the dominant active normal type of faulting of the Mygdonian Basin, or to intense water pumping in the basin for irrigation purposes. Finally the combined action of the above two factors could be the most fair interpretation of the detected subsidence.

Oreokastro and airport areas
The deformation pattern clearly evident in both images for this zone (Fig. 8) is located in the NNW (Oreokastro) and SE (area of airport) suburbs of Thessaloniki and presents a particular interest due to a subsidence rate 10 to 20 mm.yr −1 . For these areas there are no bibliographic references or technical reports concerning subsidence phenomena. The Oreokastro area is covered by Neogene and Quaternary formations and a gentle landscape with the presence of small Nat. Hazards Earth Syst. Sci., 8, 779-788, 2008 www.nat-hazards-earth-syst-sci.net/8/779/2008/ gullies. The western extremity of the Asvestochori-Polichni active normal fault crosses the area with a WNW-ESE direction and is dipping to NNE. The subsided area belongs to the northern hanging wall. A series of micro-earthquakes could be related to this fault. The geological basement of the area around the airport area is identical to that of Sindos-Kalochori, namely pebbly gravel, sands, clays and coastal deposits. It appears that there is no active fault in the area, as indicated by the neotectonic map (Thessaloniki sheet). In order to fully interpret these last two areas additional investigations and ground-based information is needed.

Conclusions
The application of space-borne interferometry to the area of Thessaloniki (Greece) to detect and assess possible ground deformation phenomena is presented in this article. There is good agreement between the interferometric results and those available from other sources and summarised in Table 2.
In the Sindos-Kalochori area, InSAR and PSInSAR seem to have underestimated the deformation rates. Similar observations were outlined by Crosetto et al. (2007), showing that interferometric measurements (using standard procedures without prior information on the observed deformation pattern) present limitations for high rates of deformation. This is because an insufficient image acquisition sampling does not allow correcting for the phase ambiguities if the deformation rate is too high (generally, several cm.y −1 ) in the PSInSAR processing and to a possible excessive filtering of conventional interferograms.
Even though the two advanced interferometric techniques (both stacking and PS) did not allow the identification of the causes that provoked the subsidence phenomena, they both gave a clear synoptic view of quantitative and qualitative ground deformation, the areal extent of the subsidence. In addition, they helped to locate areas of previously unknown subsidence, as in the northern and south-eastern suburbs of Thessaloniki. Finally, they suggest where future control and detailed field studies are necessary.