InSAR Campaign Reveals Ongoing Displacement Trends at High Impact Sites of Thessaloniki and Chalkidiki, Greece

Abstract

We studied the broader area of Thessaloniki in northern Greece and Chalkidiki and performed an InSAR campaign to study the surface deformation phenomena that have been known to exist for at least two decades. Sentinel-1 data (2015–2019) together with drill measurements were exploited to focus on specific sites of interest. Our results indicate an ongoing displacement field. At the region of Kalochori and Sindos—where intense subsidence in the 1990s was previously found to have had a natural surface rebound in the 2000s—a new period of subsidence, caused by the enlivenment of the groundwater overexploitation, was reported. The uplifting trend of Oreokastro is still active and subsidence in Anthemountas graben is ongoing; special focus was set on the Makedonia Airport, where significant displacement is occurring. The study also reveals a new area at Nea Moudania, that was not known previously to deform; another case corresponding to anthropogenic-induced surface displacement. Thessaloniki is surrounded by different persistent displacement phenomena, whose main driving mechanisms are anthropogenic. The sensitivity of the surface displacements to the water trends is highlighted in parts of the study area. Results highlight the plan of a water resources management as a high priority for the area.

1. Introduction

Thessaloniki is located in Northern Greece and is the second largest city of the country (Figure 1). Past studies of the broader area have highlighted several subsiding deformations occurring in the wider area (e.g., [1,2]). Synthetic Aperture Radar Interferometry (InSAR) techniques have been successfully applied using data from radar satellites for the monitoring of surface displacement in several cases globally (e.g., [3,4]) and also in Greece (e.g., [5,6,7,8]). The advent of technological advancements of radar satellites has provided significant input data and has given rise to a series of InSAR studies at the wider region of Thessaloniki. Some of these studies have detected and measured the ground surface displacement (e.g., [9,10,11]). Other InSAR studies used in-situ data and ground campaigns to interpret and find the cause of the reported displacement signals [12,13,14,15,16,17]. It is widely accepted that underground water activity can cause surface displacement [18]. This paper aims to study the surface deformation occurring in the vicinity of Thessaloniki (Figure 1) and to consider the role of anthropogenic activities in the evolution of these deformation phenomena. Aquifer overpumping had been identified as the main cause of deformation for the areas of Kalochori–Sindos, for the broader area of Oreokastro, as well as for the Anthemountas basin (Figure 1) [12,13,14,15,16,17].

In this study, we exploit the Sentinel-1 satellite, operated by the European Space Agency (ESA), to answer the question of whether there is a considerable continuing, persisting hazard at the sites of interest or whether it has diminished or been replaced by other types of geo-hazards. The geo-hazard interpretation is strongly supported by using data from in-situ conducted campaigns in the specific sites of interest. Focus is set on the areas of Kalochori–Sindos, Oreokastro, and Anthemountas, which were highly impacted by surface displacement, the latter resulting occasionally in numerous damages of houses, roads, and industrial buildings. Moreover, new deformation evidence has emerged at the area of Nea Moudania (Figure 1). No previous study has reported any type of deformation occurring in this area. Apart from the scientific curiosity that is the dominant driver for the study, significant motivation has also been the fact that the broader Thessaloniki metropolitan area consists a hub for the Balkans and south-eastern Europe; here lay important economic activities and critical infrastructure, and many actors are placed belonging to the primary and touristic sectors as well.

2. Methods and Data Used

We used 91 Sentinel-1 satellite images, descending pass, orbit no. 7, and frame 457, spanning the period 12 December 2015–15 December 2015 (

Table 1. Radar data used.

Date	Normal Baseline	Date	Normal Baseline	Date	Normal Baseline
12 December 15	−88	28 February 17	37	26 November 18	15
24 December 15	6	12 March 17	−15	08 December 18	−38
05 January 16	48	24 March 17	84	20 December 18	47
17 January 16	94	05 April 17	77	01 January 19	45
29 January 16	89	17 April 17	75	13 January 19	75
10 February 16	7	29 April 17	49	06 February 19	19
05 March 16	33	11 May 17	−26	02 March 19	−19
17 March 16	112	20 October17	0	14 March 19	117
29 March 16	65	01 November 17	−14	26 March 19	58
10 April 16	35	13 November 17	53	07 April 19	−5
22 April 16	5	25 November 17	45	19 April 19	−40
16 May 16	73	07 December 17	37	01 May 19	−40
28 May 16	78	19 December 17	−51	13 May 19	98
09-Jun-16	73	31 December 17	−31	06-Jun-19	61
15-Jul-16	95	05 February 18	30	18-Jun-19	58
27-Jul-16	103	01 March 18	−15	30-Jun-19	16
08-Aug-16	59	11-Jun-18	18	12-Jul-19	86
20-Aug-16	−12	23-Jun-18	62	24-Jul-19	98
01-Sep-16	81	05-Jul-18	163	05-Aug-19	82
13-Sep-16	61	17-Jul-18	85	17-Aug-19	−25
25-Sep-16	23	29-Jul-18	77	29-Aug-19	27
07 October16	33	10-Aug-18	24	10-Sep-19	81
19 October16	66	22-Aug-18	116	22-Sep-19	48
31 October16	−36	03-Sep-18	44	04 October19	17
12 November 16	36	15-Sep-18	71	16 October19	46
24 November 16	40	27-Sep-18	23	28 October19	31
06 December 16	22	09 October18	19	09 November 19	36
11 January 17	−30	21 October18	34	21 November 19	68
23 January 17	68	02 November 18	49	03 December 19	44
04 February 17	68	14 November 18	15	15 December 19	−29
16 February 17	86

). The detection and measurement of surface displacement was achieved by performing Multi Temporal InSAR (MT-InSAR). Our processing chain was Parallelized Persistent Scatterer Interferometry (P-PSI; Papoutsis et al. [19]) which is a novel, distributed processor for the fully automated and computationally efficient assessment of line-of-sight ground velocities through PSI, tailored to scale to the large multi-temporal archive of Sentinel-1 data.

Our P-PSI implementation invokes two main software packages that have been modified to parallelize the execution of the different processing tasks. We use the InSAR Scientific Computing Environment (ISCE) with the topsStack processor [20] for creating a stack of coregistered Single Look Complex (SLC) datasets. Several processing steps with ISCE software were parallelized, including the generation of overlap interferograms and the computation of geometrical offsets among all slave SLCs and the stack master, for precise coregistration. The SLC stack was then used for PSI analysis with a parallelized implementation of the Stanford Method for Persistent Scatterers (StaMPS; Hooper et al. [21,22]), a process that produces ground velocities and deformation histories for the Persistent Scatterers (PS) that have been detected. The full Sentinel-1 frame was divided in sub-regions, i.e., “patches”. Our processing chain exploited multi-core programming techniques to execute simultaneously the processing steps of the Stanford method for Persistent Scatterers in multiple patches. StaMPS relies on the spatial correlation of the deformation, rather than on a temporal linear deformation model. It is therefore suitable to also detect non-linear deformation.

To estimate the topography effect, the Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM) [23] with a resolution of 1 arc-second, was employed. In order to have a reasonable trade-off between a larger pixel spacing and pixels with random phase, we opted for a threshold of 20 m for the maximum topographic error. To exclude pixels considered as “noisy” or pixels selected as PSI, due to signal contribution from neighboring ground resolution elements, a threshold of 1 was selected for the phase noise standard deviation for all pixel pairs. Therefore, pixels exceeding this threshold were excluded from further steps of PS processing.

The open-source Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) [24] was employed, by implementing a phase-based linear correction on interferograms produced for PS processing. Since the areas of interest investigated are characterized by divergent topography, the linear correction applied, which compensates for the tropospheric delay relative to the topography, provides a satisfying accuracy in the tropospheric delay correction.

The outcome of P-PSI is a collection of pixels (PS) with stable phase over time. These points provide estimates of annual velocities and a quality measure named velocity standard deviation (Figure S1). High values for the latter indicate either a noisy estimate or provide an indication of non-linear deformation. In addition, for each PS, a deformation time-series in the period 2015–2019 was provided that highlights various deformation trends.

3. Results

The interferometric pairs that were formed are schematically shown in Figure 2 with respect to their temporal and perpendicular baselines. We processed the entire set of available Sentinel-1 images for the area and period of interest and identified a dense and spatially well distributed network of 1,551,817 stable-phase pixels, for which we estimated ground displacements (Figure 3).

Since PSI estimates are a relative measurement on both the spatial and temporal domain, we used in our PSI the SAR acquisition of 20 October 2017 as the master scene for temporal reference, and the VRAS continuous GPS station (Figure 3) as a spatial reference point. VRAS has data for four years in the period 2013–2016, is 80 km away from Thessaloniki, and is located in an area for which there is no prior evidence that it has been deforming. VRAS data shows that movement in the Up direction is negligible (~0.1 mm/yr) with some seasonal movement, as expected. The stability of the VRAS station ensures us that we do not need to take account of a constant bias in our velocities estimations.

The results indicate that there is surface deformation at the areas of interest, enclosed by the blue rectangles shown in Figure 3. More specifically in the area of Kalochori–Sindos there is a subsidence signal, Oreokastro has an uplifting trend while a subsidence trend appears to exist at the basin of Anthemountas. A new area not previously known to deform is the Nea Moudania plain (Figure 1 and Figure 3).

The following sections focus on each of these areas, analyze the findings from InSAR processing and ground truth data when available, and attempt to provide interpretations on the reported signals.

3.1. Kalochori–Sindos Region

The villages of Kalochori and Sindos are located at the west of Thessaloniki (Figure 1). The area is known for the intense industrial activity. It is located on Quaternary deposits (Figure 4) with the sediments consisting of sands, clays, and silty clays [25,26]. The upper 5 to 15 m of the Quaternary deposits are occupied by fine-to-medium grained sand intercalated by silty sand horizons. Underneath, a 25 to 35 m thick black-grey silty sand horizon with organic materials extends, forming an impermeable boundary, determinant for the hydrogeological conditions of the site. The lower strata, down to the Neogene bedrock, consist of alternating brown sand and black-grey silty sand horizons [26,27].

From a hydrogeological point of view, two aquifers can be identified: a shallow aquifer unconfined within the sandy layer of the upper strata, and a deeper confined aquifer system within the brown sand of the lower strata. The quality of the water of the shallow aquifer is poor, in contrast to the high quality water of the second confined aquifer system, extending below the impermeable organic silty sand horizon. The overexploitation of the confined aquifer has been proven as the main cause of the land subsidence phenomenon occurring at the site.

Ground truth techniques had been used to study the subsiding pattern that was known since the 1960s [25,29,30,31,32]. The 1990s InSAR displacement was first identified by Raucoules et al. [10] using acquisitions from the European Remote sensing Satellite (ERS) satellites. A detailed investigation by Raspini et al. [13] used the ERS satellite results (1990s dataset), together with in-situ data for the interpretation of the detected subsidence. However, in a later study, an uplifting pattern was identified in Svigkas et al. [11] corresponding to the post-2000 period, covering the Environmental Satellite (ENVISAT) era. This change in the deformation trend was attributed to the natural recharge of the aquifers, due to the cessation of groundwater pumping [14]. A Global Navigation Satellite System (GNSS) study covering the period 2013–2015 showed a mixture of subsidence and uplift displacement pattern [33].

The Sentinel-1 interferometric results presented in Figure 5 indicate that the uplifting trend of the 2000s has been reversed again and now a new period of subsidence is occurring in the area. Maximum displacement values in the vicinity of Kalochori are of the order of −10 mm/yr, while in the broader area, maximum values reach a rate of −14 mm/yr.

Hydrogeological data of the region enabled the validation and interpretation of the PSI Sentinel-1 based results. The drill data of the same period indicate a decrease in the underground water levels (Figure 6), which in most drills appears to have been initiated between 2016 and 2017. The amount of dropdown from November 2015 to November 2019 agrees spatially with the areas of significant displacement (Figure 7) depicted by the PSI processing.

In order to have a holistic and temporally complete view of the phenomenon that has been occurring at Kalochori over the last decades, a graph presenting the groundwater level of a drill, together with the surface deformation estimated at the exact neighborhood of that drill, is presented in Figure 8. The plot shows that the subsidence taking place during the 1990s was followed by a rise of the hydraulic head—due to natural rebound—which in turn resulted in a surface uplift the following years [14]. The current analysis of satellite and in-situ hydrogeology data, indicates a potential onset of a new cycle of subsidence, derived from the recent downward trend of the aquifer water level (Figure 8). The subsidence trend of Sentinel-1 SAR results appeared with a time-lag of ~2 years with respect to the initiation of the lowering of the hydraulic head. Similarly a ~2-year time-lag had been also identified at all former changes of the trend of the curves; namely when the aquifer started to recover and when the recovery trend was reduced during 2006 (Figure 8).

3.2. Oreokastro Industrial Zone

Oreokastro is located at the exact proximity of Thessaloniki (Figure 1) and is an area strongly related with the industrial activities of the city. The town was founded on a Pliocene lake and fluvial deposits. The industrial zone of Oreokastro is on Pliocene clays intercalated by sand, lying above the Pre-Alpine bedrock (Figure 9 and Figure 10). The main tectonic features are the Efkarpia fault and the Asvestochori fault, both considered as potentially active [34]. No major seismic event is known to exist and the area is mainly related with microseismicity [35,36,37].

InSAR techniques had previously revealed that the area has been deforming [10]. The subsidence of the 1990s was followed by an uplifting trend after 2000. In Svigkas et al. [15], by combining remote sensing analysis with seismological, geological, and hydrogeological data of the area, it was shown that the observed displacement—which also caused building failures—was related with the activity of aquifers. Also, the displacement was found to be a fault-controlled pattern, in a sense that the faults at the area were most likely acting like water-flow barriers.

Sentinel-1 PSI results indicate a continuation of the uplifting pattern (Figure 11), previously detected to occur during the 2002–2010 ENVISAT based monitoring period. Specifically, the maximum uplifting rate was found to be as high as 8 mm/yr, which is essentially the same order of value as the one observed during the ENVISAT era (9 mm/yr).

To analyze with higher detail the displacement pattern, velocity profiles were created (Figure 12). It can be seen in all the profiles under study that the Asvestochori fault (Figure 9) is affecting the displacement trend.

Based on the ENVISAT results, Svigkas et al. [15] proposed that at the area there is an unknown feature highlighted by abrupt changes in the deformation trends. The interferometric analysis conducted showed that between the uplifting pattern and the Asvestochori fault, there was a thin subsiding zone just north of the fault and before the expression of the uplift. The same pattern appears also in the present study (Figure 11) of the Sentinel-1 based interefrometric processing. However this time the deformation pattern presents a slightly different geometry than the one of the period 2003–2010, a phenomenon that might be related to the difference of underground water level activity, during the two monitoring periods. An explanation for this could be the different amounts of the underground water level uprising that might cause a different surface displacement pattern, depending also on the various permeabilities of the earth strata that the water flow meets.

3.3. Anthemountas Basin

Anthemountas basin is located west of the city of Thessaloniki (Figure 1). In the area the development trends had already started before the 1990s [40]. Here lay the longest tectonic structures close to the city of Thessaloniki; the northern branch of Anthemountas fault (F1 in Figure 13) that has a length of 21 km, considered inactive [34], and the south branch of Anthemountas fault (F-2). Morphotectonic studies indicate it as being an active fault structure [34,41].

Regarding the geological structure of the site, the top Quaternary formations consist of alterations of clastic and fine sediment [34,42,43]. At the Quaternary deposits (Figure 13 and Figure 14), three main horizons can be identified: a) coarse to fine sand b) clay to silty clay horizon, and c) Neogene deposits characterized by coarse to fine sands. The Neogene deposits consist of a top sand and gravel layer and a bottom clay-marl layer [42,44]. The area has a Mesozoic bedrock where phillites, gneiss, and granites lie. The thickness of the sediments decreases when moving towards the east.

In the upper strata of the Quaternary deposits, there is a shallow aquifer system interacting with a deeper semi-confined aquifer extending down to the Neogene formations [45]. The water demands of the area have caused a decrease of the groundwater level.

Raucoules et al. [10] used acquisitions of ERS satellites and by applying InSAR time-series detected subsidence during the 1990s. It was found that the displacements were extending from the coastal zone of the graben throughout the entire structure, up to the far East of the basin [16]. Raspini et al. [12] attributed the subsidence phenomena of the 1990s–that appeared to have also caused building failures–to aquifer overpumping. The subsidence in Anthemountas continued after 2000 and had been occurring continuously and increasingly for at least 20 years [17].

The observations based on PSI processing of Sentinel-1 images are presented in this study and indicate that the subsidence phenomenon is still active (Figure 15). The maximum subsidence value at the site appears to be −18 mm/yr, the same as the maximum values previously detected during the ENVISAT imagery era [17]. Detailed velocity profiles are shown in Figure 16.

An important infrastructure located at the Anthemountas graben is the International Airport of Thessaloniki (cyan rectangle Figure 15). It is shown here that during the monitoring period, the airport has a subsidence signal reaching −12 mm/yr (Figure 17). This is located at the intersection of the two airport runways. The current terminal of the airport was found to have been subsiding with a rate of −8 mm/yr, during 2015–2019.

3.4. Nea Moudania

Nea Moudania plain drew the attention of the researchers for the first time, due to recently occurring deformation pattern identified at the site. The wider region belongs to the Paeonian subzone, part of the Axios geotectonic zone. The Litho-stratigraphy of the study area is divided into the following formations [47,48].

Quaternary formations, located mainly near the coastal areas, including alluvial deposits (sands, clay), coastal deposits (beach ridges etc.), lagoon deposits (sands, clayey sand), and alluvial fans (Figure 18). Furthermore, along the streams, two terrace systems can be identified. The lower consists mainly of sands, schists, and limestone pebbles and the upper of metamorphic rock pebbles. Neogene formations cover the pre-Neogene bedrock formations. They can be distinguished into the lower red clays series, occupying most of the study area and the upper marls with sandstone series, outcropping at the west of the study area. The red clays series consists mainly of red silty clays with fine-grained quartz and mica, intercalated by sand, marl, and conglomerate lances. The marl with sandstone series consists of altering beds of sand, clay, and clay-marl. Regarding the bedrock, it consists of party recrystallized bluish limestone, greenish-brown two mica gneiss, and muscovite gneiss.

Two aquifer systems can be identified at the study area: a shallow unconfined aquifer and a deep confined aquifer [49]. The shallow phreatic aquifer mainly occupies the alluvial Quaternary deposits and the deep confined aquifer the permeable strata intercalated into the Neogene formations, namely sand, conglomerate, and sandstone layers. The confined aquifer system’s thickness exceeds 200 m in depth [50]. The water provided laterally by the karstic aquifers of the bedrock formations, as well as the water infiltrated from the streams and the precipitation, are the main sources of water for the aquifers.

The Sentinel-1 based PSI results indicate a significant surface displacement (Figure 19). The subsidence rates have a maximum value of -23 mm/yr, during the time span of 2015–2019.

The small extent of the alluvial deposits, in combination with the overexploitation of the shallow aquifer over the past decades has led to the draining of the phreatic aquifer [51]. Since the 1990s, numerous deep wells have been constructed in the area to meet the water needs [52]. The confined aquifer system of the Neogene formations is considered to be the main source of underground water.

According to data provided by Veranis and Xatzikirkou [50] during the late 2000s, the distribution of the piezometric lines have been clearly indicating the initiation of a sea water intrusion, along the coastal area close to Sozopoly, as well as a depression cone at the plane area NW of Nea Moudania town, both due to the overexploitation of the aquifers (Figure 20). The depression cone had been covering the plane area extending up to the north close to Portaria and Zographou villages with a maximum depression of approximately -15 m. According to the 2016 data provided by Gavrihlidou [53], since 2008 the conditions radically changed. The zero level contour line (blue line in Figure 20) moved half way towards Nea Triglia village and the depression cone extended much more to the NW including clearly Portaria and Zographou villages (Figure 18).

Unluckily, the available groundwater level data referring to 2016, do not allow the production of piezometric curves covering the same area with the one covered by the 2008 data, but nevertheless an expansion of the depression cone is clearly evident (Figure 20). Also, the changes observed in the flow direction of the ground water prove the expansion of the depression cone. As shown from the comparisons of the two upper inset maps in Figure 20 (Contours May 2008 and May 2016), the flow direction heading to the sea during 2008 changed by turning to a SE flow, towards the center of the plane area. As the extent of the depression cone coincides with the spatial distribution of the PS indicating subsidence, it can be inferred that the deformations are connected with the consolidation of the formations, triggered by the overexploitation of the aquifer.

4. Discussion

In all cases presented, the detected surface displacements are mainly related to the trends of the hydraulic head. During the recharge of the underground water level, there is an increase of the pore fluid pressure and a decrease of the effective stresses. This is the period where we can have a surface uplift, like for example the one identified in Oreokastro during the 2000s and also from ~2015 to 2019 or like in Kalochori–Sindos during the 2000s. On the contrary, a decrease of the groundwater level causes an increase of the effective stresses and as a result, subsidence. As identified by the InSAR studies this was occurring at Kalochori–Sindos during the 1990s and from ~2015 to 2019, in Oreokastro during the 1990s and in Athemountas during 1990s, 2000s, and during the period of 2015–2019.

In several areas around the world artificial aquifer recharge strategies are followed aiming to overcome the subsidence phenomena. It has been shown that at least for the areas of Oreokastro and Kalochori–Sindos, the surface rebound and the return to a non-deforming state, can be achieved naturally. This could potentially be the result of the generation of rational water management. However, if the groundwater decline continues, the subsidence will continue and the failures will be ongoing in the future. Even though significant hazards have occurred at those areas, the detected natural recharge is a positive factor to be taken into account; at those areas, there is still a regime of rebound deformation. If at some point the underground water level decrease becomes severe enough to cause stresses that exceed the preconsolidation stresses, a rearrangement of the granural material of the aquitards could occur. Then, the area will start experiencing inelastic deformation (e.g., [54,55,56,57,58]). In this case, the results would be irreversible or only partially reversible, even in the case of an aquifer recharge [59].

5. Conclusions

Based on Sentinel-1 and new in-situ underground water level data, we analyzed several areas in the broader area of the city of Thessaloniki, which were found to be deforming during the years of 2015–2019. For the area of Kalochori–Sindos, even though there was a natural rebound after 2000—pausing the subsidence rates of the past—this study showed that a clear subsidence trend has been initiated again. It also became clear that the surface of the area is very sensitive to the changes of the aquifer trends. Moreover, the study shows that the broader area of Oreokastro is still uplifting since the pause of its subsidence trend that occurred after the 1990s. The groundwater level activity enabled us to see discontinuities whose nature and cause could be related to tectonic or other types of geological structures. At this stage, we consider that a geophysical survey at Oreokastro might provide additional insights. The displacement of Anthemountas basin is still ongoing, with the same magnitude as the values detected during the 2000s. Another important fact presented here is that the displacement trend still continues throughout the entire basin and not only on the coastal part of Anthemountas. Significant subsidence occurs also at the area of the International Airport of Thessaloniki where as before, a continuous monitoring seems necessary. A new area that was found for the first time to be affected by subsidence is the Nea Moudania plain. Our analysis highlights that it is also attributed to an anthropogenically derived hazard, linked with water overpumping.

Since most of the study areas presented here have been deforming, occasionally with different trends, for at least ~25 years, it becomes clear that a water management plan should be a forefront priority in the agenda of the decision makers.