DInSAR-Based Detection of Land Subsidence and Correlation with Groundwater Depletion in Konya Plain, Turkey

Abstract

In areas where groundwater overexploitation occurs, land subsidence triggered by aquifer compaction is observed, resulting in high socio-economic impacts for the affected communities. In this paper, we focus on the Konya region, one of the leading economic centers in the agricultural and industrial sectors in Turkey. We present a multi-source data approach aimed at investigating the complex and fragile environment of this area which is heavily affected by groundwater drawdown and ground subsidence. In particular, in order to analyze the spatial and temporal pattern of the subsidence process we use the Small BAseline Subset DInSAR technique to process two datasets of ENVISAT SAR images spanning the 2002–2010 period. The produced ground deformation maps and associated time-series allow us to detect a wide land subsidence extending for about 1200 km² and measure vertical displacements reaching up to 10 cm in the observed time interval. DInSAR results, complemented with climatic, stratigraphic and piezometric data as well as with land-cover changes information, allow us to give more insights on the impact of climate changes and human activities on groundwater resources depletion and land subsidence.

1. Introduction

Protection of water resources and mitigation of the risk associated with their overexploitation represent challenging issues worldwide. In Turkey, one of the Mediterranean countries more severely affected by water scarcity [1], the Konya basin, in central Anatolia, is facing with serious water-related problems, due to climate changes [2] and lack of efficient water management policies especially in the agriculture sector [3].

Moreover, growing irrigation needs are often satisfied by uncontrolled practices of groundwater extraction: among 92,000 wells located in the Konya basin, it is estimated that 66,000 are illegal [4]. Such an aquifer overexploitation has important environmental implications for the region. Groundwater overuse and depletion not only directly impact on the hydrological cycle [5,6], but have several negative side-effects on biodiversity, ecosystem health, soil deterioration and drying up of lakes and wetlands [7,8]. Moreover, excessive groundwater withdrawal makes the region extremely vulnerable to natural and man-made hazards. As a result of piezometric level decrease, the Konya area is extensively affected by two different processes of land subsidence: (1) the progressive settlement induced by the compression of unconsolidated alluvial sediments; (2) sudden collapses of subterranean cavities formed by the combination of rock dissolution and suffosion processes. These phenomena cause damage to buildings and infrastructures, leading to high socio-economic costs for the concerned communities. To our knowledge, most of the studies conducted on the Konya area have focused on sinkholes investigations [9,10,11,12] while only a few studies on the analysis of the land subsidence caused by consolidation of alluvial fans [13,14].

In such a context, it is clear that obtaining information on the spatial and temporal evolution of land subsidence plays an important role on hazard assessment and may significantly contribute to risk prevention and mitigation. However, the continuous monitoring of surface deformations at regional/basin scale is a resource-intensive task to be accomplished by using traditional ground-based techniques. The use of satellite based information in Earth’s surface studies represents an important breakthrough, fulfilling the requirements of high spatial and temporal coverage. In particular, remote sensed data acquired by satellite Synthetic Aperture Radar (SAR) sensors and processed through Differential SAR Interferometry (DInSAR) techniques [15] greatly contribute to surface deformation analyses. DInSAR firstly relied on the use of a pair of pre- and post- event SAR acquisitions in order to retrieve ground displacements caused by single trigger events, i.e., earthquakes and volcanic unrest [16,17]. Since such early applications, relevant algorithmic advances have been achieved and several DInSAR approaches fully exploiting large SAR datasets to generate time-series of deformations have been developed. These techniques can be grouped in two main categories: the Persistent Scatterers (PS) [18,19] and the Small Baseline (SB) [20,21]. In particular, the PS techniques focus on point-like targets that are not significantly affected by decorrelation noise [22], and are suitable to monitor deformations affecting point-wise structures. Conversely, the SB methods rely on selecting SAR data pairs with short temporal (i.e., time interval between two acquisitions) and perpendicular (i.e., the distance between two satellite passes) baselines in order to reduce the noise, and allow the investigation of displacements over distributed scatterers (DS) on the ground.

Subsidence studies have greatly benefited from the application of both PS and SB DInSAR techniques which have been used for the analysis of the temporal evolution of surface deformations induced by natural and anthropogenic factors, from intense human activities in urban areas [23,24] to rock dissolution processes (i.e., karst subsidence) [25,26,27]. In particular, several works have focused on the application of such DInSAR techniques to the detection and monitoring of land subsidence induced by groundwater extraction in many regions of the world. Chaussard et al. [28] and Castellazzi et al. [29] used DInSAR deformation time-series of ALOS and RADARSAT-2 data, respectively, to study the land subsidence caused by massive groundwater over pumping in the entire central Mexico. Dehghani et al. [30] exploited both descending and ascending ENVISAT datasets to measure the magnitude of groundwater related subsidence in the Tehran basin, Iran. The effectiveness of DInSAR techniques in mapping and monitoring surface deformations has been proved in many applications from Europe to Asia [31,32,33] showing that ground subsidence triggered and/or accelerated by aquifer overexploitation is a common hazard impacting extensive areas worldwide.

In this paper, we apply the SB approach referred to as Small BAseline Subset (SBAS) [20], which is capable of detecting and measuring displacements affecting extended rural areas, thus resulting particularly suitable for the regional scale analysis of the land subsidence affecting the Konya plain. Through the generation of displacement time-series and maps of deformations rate we analyze the temporal and spatial pattern of the land subsidence caused by groundwater resource overexploitation. Furthermore, to investigate the correlation between land subsidence and its conditioning/triggering factors, we adopt a multi-technique approach, based on the use of several data, from traditional ground-based measurements (piezometric data) to remote-sensed data which recently opened new perspectives in water resources studies (GRACE data), from meteorological data to Landsat images and CORINE Land Cover maps used for extracting land cover information. Integration of such a variety of data, embedding local- and regional-scale information, allows us to get a holistic insight into the climatic, hydrogeological and human dynamics of the whole area.

2. Study Area

2.1. Geographic and Climatic Setting

The study area is located within Konya basin, central Anatolia, which is the biggest endorheic basin in Turkey extending over an area of about 54,000 km², with elevations ranging from 900 m to 3500 m above sea level (a.s.l.) (Figure 1). The basin lies on the northern side of the Central Taurus Mountains which have an important role on surface and sub-surface water supply. The Konya basin has a typical arid to semi-arid climate with a long-term average yearly precipitation of 380 mm, ranging from 250 mm in the plain areas to more than 1000 mm in the mountainous areas in the south [5], and an evapotranspiration of about 550–600 mm per year. Winters are generally cold and wet whereas summers are hot and dry. Despite having such a dry climate with low precipitation input, there are extensive and intensive agricultural activities in the region, whose main source of water is the groundwater which is extracted via more than 92,000 groundwater wells distributed around the basin [34]. In the last years the increase of irrigated land (about 3700 km²) related to some cultivations like sugar beet and corn caused an increase of groundwater consumption. The estimated groundwater consumption in 2010 was about 2.7 billion m³ [35]. Numerous smaller fresh/brackish wetlands have dried out in the last decades as consequence of groundwater depletion [5].

2.2. Geological and Geomorphological Setting

In the study area, three main geological sectors can be identified: the Konya-Karapinar plain, the volcanic field and the Obruk plateau (Figure 2). The Konya-Karapinar plain is characterized by the presence of Late Miocene to Quaternary formations, the latter mainly made up of lacustrine deposits. This formation, deriving from the fill of Pleistocenic lakes, is mainly made of unconsolidated clay, silt and sandstone. In the plain it is also possible to find aeolian and alluvial fan deposits [9,36,37]. The Late Miocene to Pliocene stratigraphy starts with a basal conglomerate and continues upward with lacustrine limestone and weakly cemented marl and limestone alternation [12]. At S-E of Karapinar, a volcanic field with age variable from Miocene to Quaternary, whose activity is related to the extensional phase that originated the Konya basin, is present; cones and other volcanic structures are visible in the study area. The Obruk plateau is characterized by Mesozoic to Oligo-Miocene formations. The Oligo-Miocene formations are composed by continental clastics, tuff, limestone and evaporites at the bottom and gypsum-shale alternation at the top [12]. Mesozoic formations are sparsely distributed in the western part of the Obruk plateau and along the southern boundary of the Konya plain, and are made up of carbonates and metamorphic rocks [38].

From the geomorphological point of view, the landscape of the study area is conditioned by several active processes: volcanism, water and wind erosion, and rock dissolution (Figure 2). Among these, water erosion and dissolution-induced subsidence are the most widespread phenomena in the study area, while aeolian processes are effective in relatively local sectors of the Konya basin.

The dissolution-induced subsidence is represented by the karst depressions called ‘sinkholes’. Sinkholes are common geomorphic features in the study area (Figure 2). These landforms are locally called “Obruks”. This term is also accepted at international level for describing the large sub-circular closed depressions observed in the Central Anatolia [39]. There are two groups of obruks according to [12]: the “old” obruks, characterized by large dimensions (up to 640 m wide), for which the time of formation is unknown (Figure 2B); the “young” obruks that have been formed during the last decades and have diameters of less than 50 m (Figure 2A). The “old” obruks are concentrated in the Obruk Plateau, a horst slightly higher than the surrounding grabens formed by the Konya plain and the Tuz Lake basin. These “old” obruks can be classified as caprock and bedrock collapse sinkholes [40]. The analysis of geological, geophysical, hydrogeological data and the groundwater chemical and isotopic compositions indicates a hypogenic origin for these large obruks [12]. There are not evidences of subsidence activity in these sinkholes so far and they may be defined as relict karst landforms. On the other hand, most of the “young” obruks can be classified as cover collapse sinkholes [40] and their formation is associated to the human-induced groundwater withdrawal [9]. In the last decade, these “young” obruks have caused social alarm in the Karapinar plain. Indeed, 13 collapse sinkholes up to 25 m in diameter occurred from 2006 to 2009. They are especially concentrated in and around the Seyithaci settlement (Figure 2A). According to the occurrence of these sinkholes, we can expect future sudden collapses if the water lowering trend will continue. Future sinkholes could damage buildings and infrastructures in the study area and represent a hazard for the local population.

Another dynamic process in the region is the sand dunes (Figure 2). The long term displacements of the sand dunes have caused the formation of sand hills and dunes [37,41]. The movement direction of the very active sand dunes, that are located at the south of Karapinar settlement (Figure 2C), changes depending on the dominant direction of the seasonal wind. However, it is observed to be mainly from south to north [37].

In the study area it is also possible to observe some active faults (Figure 2). In particular, the Seyithaci fault system is a dynamic process delimiting the Obruck plateau from the plain at north of Karapinar. The fault has a vertical component character and follows a highly linear form between the floors of plateau and plain [42]. The fault system is defined as active and is known to establish a surface rupture during Holocene [43], and has an important effect in the evolvement of Karapinar plain and in the formation of sinkholes because the karstification of the bedrock may be developed through fault structures.

3. Data and Methods

Different data and information have been integrated to carry out a comprehensive analysis of the phenomena under investigation, from satellite images (

Table 1. Satellite-based data.

Satellite	Band/Parameters	Final Product	Time Interval	Revisit Time	Spatial Resolution
ENVISAT—ASAR	C-band	Ground displacements	2002–2010	35 days	20 × 4 m2
Landsat-5/-8	Multispectral	Qualitative land cover changes/mNDVI	May 2000; May 2015	16 days	30 m
GRACE	Gravity field	Groundwater storage	2002–2010	30 days	350 km
Shuttle Radar Topography Mission (SRTM) [44]	C-band	DEM 3-arc sec	2000	Single mission	90 m

) to data acquired through in situ monitoring techniques (

Table 2. Ground-based data.

Monitoring	Product	Time Interval	Measurement Frequency
Cores	Stratigraphic data	-	-
Wells	Groundwater level	2003–2010	Monthly
Meteorological stations	Rainfall, Temperature	1935–2014; 2002–2010; 1950–2014	Monthly

).

3.1. SAR Data and DInSAR Processing

The study area has been investigated by exploiting the European Space Agency (ESA) archives of raw data collected by the ENVISAT-ASAR satellite (Table 1). These data were acquired with a 35-day revisit period at a look angle of 23° and with a C-band frequency (5.65 cm wavelength). In particular, we processed two datasets spanning overall nearly eight years: 38 SAR images taken from 5 December 2002 to 9 September 2010 on descending orbits (Track 436, Frame 2838), and 27 SAR images acquired from 14 December 2002 to 5 June 2010 along ascending orbits (Track 71, Frame 746). We applied the SBAS—DInSAR technique [20] by selecting SAR data pairs characterized by a perpendicular baseline <400 m and a temporal baseline <1500 days, and generated 111 and 73 interferograms from the ENVISAT descending (Figure 3A) and ascending (Figure 3B) datasets, respectively. As a next processing step, we computed differential interferograms (embedding only the information relevant to the deformation phase component). The phase component related to the topography of the area was estimated and removed using a Digital Elevation Model (DEM) produced from the Shuttle Radar Topography Mission (SRTM) with about 90 × 90 m² spacing (Table 1) [44,45].

Then, the differential interferograms were unwrapped by making use of the space-time Extended Minimum Cost Flow (EMCF) technique proposed in [46] which takes into account both the spatial and temporal relationships among the sequence of interferograms. After the phase unwrapping operation, differential interferograms were inverted by solving a Least Squares (LS) minimization problem [20]. SBAS allows retrieving the surface deformation evolution and is also capable to mitigate possible topographic artifacts and atmospheric phase disturbances [20]. In particular, the filtering operation for the atmospheric phase component is based on the observation that the atmospheric signal is highly correlated in space and poorly in time [18]. Moreover, the overall analysis was carried out on differential interferograms subject to complex multi-look operation with 4 looks in the range direction and 20 in the azimuth direction in order to significantly reduce the phase noise [45]. Accordingly, in the ENVISAT case the ground resolution of the DInSAR deformation rate map turns to be about 80 × 80 m² which is particularly suitable to investigate land subsidence occurring at regional level as within the Konya plain. According to the experience gathered in the Ebro Valley evaporite karst [26], this resolution could also be enough to identify progressive subsidence due to karst processes if the area affected by settlements is over ~50,000 m². Furthermore, the availability of both ascending and descending geometries allowed us to investigate the direction of the ground movement and compute the E-W and vertical components of displacements as follows [47]:

$$\begin{array}{c}{d}_{vert}=\left(d_{desc}+d_{asc}\right)/2cos\vartheta \\ d_{E-W}=\left(d_{desc}-d_{asc}\right)/2sin\vartheta \end{array}$$

(1)

where d_desc and d_asc represent the Line of Sight (LOS)—projected deformation computed from descending and ascending datasets, respectively, and ϑ the look-angle of the satellite. We remark that the retrieval of the N-S displacement component is very critical because of the near-polar orbit of the satellite (i.e., the flying direction of the satellite is nearly parallel to N-S) [48].

The applied decomposition method is performed on the LOS time-series [49], and is based on the consideration that over our study area the SAR data were acquired quite regularly from both orbits. Accordingly, the combination has been carried out on ascending and descending data that are temporally as close as possible each other. This operation is based on the assumption that no significant deformation has occurred between the two combined acquisitions, and this is reasonable for our case study.

3.2. Groundwater Level and Stratigraphic Data

Groundwater level data have been provided by the General Directorate of State Hydraulic Works of Turkey, which is in charge of monitoring the piezometric wells widely distributed along the Konya plain and the Obruk Plateau (Figure 1, Table 2). We used data recorded in nine wells (

Table 3. Piezometric wells and stratigraphic cores available for this study (for location see Figure 1).

Well/Core ID	Latitude (°)	Longitude (°)	Altitude (m a.s.l.)	Well/Core Depth (m)	Lithological Unit	Compressible Layers Thickness (m)	Groundwater Level Data
181	37.620	32.751	1014	n/a	alluvial deposits	n/a	Yes
9434	37.770	32.450	1048	n/a	alluvial deposits	n/a	Yes
10472	38.078	32.644	1000	n/a	alluvial deposits	n/a	Yes
20826	37.787	32.995	1004	n/a	alluvial deposits	n/a	Yes
42778	37.948	33.841	1075	n/a	alluvial deposits	n/a	Yes
52267	37.486	32.911	1025	n/a	alluvial deposits	n/a	Yes
52268	37.552	33.185	1013	n/a	alluvial deposits	n/a	Yes
2889	37.694	33.254	999	120	alluvial deposits	110	No
4566	37.769	33.540	996	180	alluvial deposits	130	No
4685	37.672	33.381	1001	150	alluvial deposits	40	No
35530	37.710	33.387	1006	100	alluvial deposits	0	No
38306	37.983	33.334	1046	100	Limestone/marl	0	No
42289	37.691	33.400	1006	80	alluvial deposits	0	No
43405	37.790	33.286	1062	135	Limestone/marl	0	No
43666	37.713	33.232	1007	170	alluvial deposits	80	No
47768	37.801	33.433	1040	120	Limestone/marl	0	No
47782	37.850	33.344	1070	130	Limestone/marl	0	No
48083	37.739	33.474	1047	95	Limestone/marl	0	No
49703	37.752	33.248	1025	140	alluvial deposits	80	No
51903	37.768	33.311	1055	140	Limestone/marl	0	No
47761	37.583	33.514	1002	100	alluvial deposits	50	Yes
52258	37.717	33.471	1014	175	limestone/marl	0	Yes

) which, overall, provide information on the characteristics of the whole Konya plain aquifer and on the areal extension of the groundwater depletion.

For the analysis of the correlation between groundwater level changes and surface movements we focused on the piezometric measurements carried out monthly in the period 2003–2010, thus roughly matching the period of SBAS-derived deformation time-series.

As regards the stratigraphy of the area, we used data reported in [9] distributed within the Konya plain and the Obruk Plateau (Figure 1, Table 3). The thickness of compressible layers was derived by considering the thickness of both clay and silt layers, located within the alluvial sediments of the plain, which usually influence the processes of subsidence related to groundwater extraction [31,32].

3.3. Meteorological Data

We used precipitation (P) and temperature (T) data acquired monthly by five meteorological stations of the Turkish State Meteorological Service (

Table 4. Meteorological stations used for the analysis (for location see Figure 1).

Station ID	Location	Latitude (°)	Longitude (°)	Altitude (m a.s.l.)	Measurement Time (Years)	Data
17902	Karapinar	37.7136	33.5277	997	2000–2014	P-T
17244	Konya	37.8691	32.4721	1032	1935–2014	P
					1950–2014	T
17248	Ereğli	37.5265	34.0482	1044	2000–2014	P-T
17900	Çumra	37.5662	32.7932	1016	2000–2014	P-T
17264	Karaman	37.1927	33.2212	1021	2000–2014	P-T

) [50,51]. In particular, for the analysis of climatic conditions during the 2002–2010 period covered by ENVISAT acquisitions we averaged data recorded in all five stations; for a more long-term analysis of the climate trend in the area we exploited the rainfall and temperature measurements collected by the Konya station since 1935 and 1950, respectively.

3.4. Landsat Data and Land Cover Maps

For the analysis of the land cover (LC) changes occurred in the study area, we used CORINE Land Cover (CLC) maps relevant to the years 2000 and 2012 [52,53]. For our purposes, the legend of the CLC maps was simplified to 8 classes as shown in

Table 5. Simplified CORINE land cover classes.

Simplified Classes	CORINE Classes Description
(1) Artificial surfaces	Urban fabric, industrial, transport and commercial units, mine, dump and construction sites
(2) Natural Vegetation	Natural grasslands, moors and heathland, Sclerophyllous vegetation, transitional woodland-shrub, inland marshes, peat bogs, salines, intertidal flats, broad-leaved forest, coniferous forest, mixed forest
(3) Bare rocks/soils	Bare rocks, beaches, dunes, sands, salt marshes
(4) Water bodies/Wetlands	Water courses, water bodies, coastal lagoons, estuaries, sea and ocean, inland marshes, peat bogs, salt marshes, salines, intertidal flats
(5) Non-irrigated arable land	Non-irrigated arable land
(6) Unknown or mixed irrigated/non-irrigated land	Complex cultivation patterns, land principally occupied by agricultural with natural space, vineyards and fruit tree
(7) Irrigated land	Irrigated land/rice field
(8) Pastures	Pastures

; the agricultural classes (5)–(7) were further classified on the basis of irrigated/non-irrigated lands.

Moreover, we exploited Landsat 5 TM and Landsat 8 scenes acquired, respectively, in May 2000 and May 2015 in order to analyze the irrigated land changes. Landsat images were available from the USGS Glovis Archive [54] as Standard Terrain Correction products (Level 1T—precision and terrain correction). For our analysis we used surface reflectance in the blue, near-infrared and shortwave-infrared wavelength domains (

Table 6. Landsat band combination used for the analysis.

Color	Domain	Landsat 5 TM		Landsat 8
		Band	Wavelength (μm)	Band	Wavelength (μm)
Blue	Blue	1	0.45–0.52	2	0.45–0.51
Green	NIR	4	0.76–0.90	5	0.85–0.88
Red	SWIR	5	1.55–1.75	6	1.57–1.65

). We generated false-color images by calibrating the brightness of each band to obtain similar colors for areas without any change between 2000 and 2015. Moreover, we produced classified Landsat images by applying a supervised classification based on the maximum likelihood method.

We have also calculated and calibrated a modified Normalized Difference Vegetation Index (mNDVI) for the two Landsat images in order to assess the changes of vegetation status in the main CORINE land cover types.

3.5. GRACE Water Storage Data

The GRACE satellite provides an estimate of the water storage anomaly by measuring the variation of the gravity field which is related to the variation of groundwater, soil moisture, snow coverage and water reservoirs.

In order to assess the groundwater storage variation, we used Level-2 data acquired from 2002 to 2010, available from the University of Colorado GRACE Data Analysis Website [55]. The used processed data represent the terrestrial water storage (TWS) anomaly (expressed in cm) calculated on a GRACE cell of about 100 km radius with a 25 km smooth filter. As reported by [56], GRACE is sensitive to TWS variation over areas smaller than its original footprint (400 km × 400 km) if such variation is significant. For instance, GRACE observations allow to detect 1 cm of TWS change within a basin of 200,000 km² (i.e., 2 km³ TWS volume change). Therefore, its detectability is of the order of 2 m water level change per 1000 km².

The Konya plain extends for about 6000 km² and the average groundwater variation is about 10 m, which results to be detectable by the GRACE cell. Moreover, we have corrected the GRACE data by using data available from the Global Land Data Assimilation System (GLDAS) [57] in order to reduce the effect of soil moisture and snow water equivalent (SWE) masses. The residual variation of TWS can be associated to groundwater and surface water masses.

We remark that GRACE data are used to integrate the information from the available piezometric wells, and to show that the water resources depletion is an issue to be addressed not only in our study area but at regional scale.

4. Results

4.1. SBAS-DInSAR Deformation Maps and Time-Series

Maps of the average deformation rate and, for each pixel, displacement time-series spanning the December 2002–July 2010 time interval have been produced (Figure 4 and Figure 5), allowing to detect a wide land subsidence extending for about 1200 km² (see the black dashed line in Figure 4) and measure its magnitude.

Maps showing rates and spatial pattern of surface deformations observed from both satellite orbits (i.e., descending and ascending) are reported in Figure 4A,B. Displacements are measured along the LOS of the satellite: negative values (yellow to red color) represent movements away from satellite (i.e., land subsidence) whereas positive values (cyan to blue) represent movements towards the satellite (i.e., land uplift). Green-color pixels are characterized by velocities within the range of ±0.25 cm/year and are considered stable. DInSAR measurements are related to a reference point located in an area assumed to be geologically stable, thus not susceptible to subsidence processes (for location see black star in Figure 4). The achieved results show an accuracy of 0.8 cm for the displacements and 0.15 cm/year for the deformation rate.

In order to ensure reliable results, we selected SAR pixels characterized by a temporal coherence greater than 0.8. However, despite the high coherence threshold applied, the SBAS processing allowed to achieve a density of 17 and 21 measure points/km² for descending and ascending datasets, respectively. Such values are considerably large if we consider that it is an area lacking the radar targets provided by man-made features, highlighting the capability of the SBAS-DInSAR approach to effectively investigate displacements in both rural and urban settings.

Maps of the vertical and (E-W) horizontal components of the surface displacements are shown in Figure 4C,D. The horizontal movement component can be considered negligible confirming that the regional deformation process mainly occurs along the vertical direction (Figure 4C). This is also clearly visible from the inspection of the achieved time-series revealing that significant vertical displacements reaching up to 10 cm have occurred, with an average deformation rate of about 1.2 cm/year for the period 2002–2010 (Figure 5).

4.2. Land Cover Change Analysis

The analysis of the land cover (LC) changes occurred in the area during the observation period covered by the SBAS-DInSAR measurements was carried out by using two CLC maps from 2000 and 2012 [53]. Figure 6 shows the spatial distribution of the CLC classes within the whole Konya plain, in 2000 and 2012, and the increase of irrigated land during this period (see blue polygons in Figure 6).

Table 7. Land cover changes computed by using the 2000 and 2012 CLC maps (see Figure 6).

Land Cover Classes	Area (km2)		Change (km2)
	2000	2012
Irrigated land	3723.4	3843.9	120.5
Unknown or mixed irrigated/non-irrigated land	509.1	441.0	−68.1
Non-irrigated arable land	1049.4	981.0	−68.4
Pastures	872.1	970.1	98.0
Natural vegetation	1351.2	1227.3	−123.9
Bare soil/Rock	47.5	156.3	108.8
Artificial surfaces	400.5	407.7	7.2
Water bodies/Wetland	149.7	75.5	−74.1
Total Area	8103	8103	-

reports the LC changes affecting the area during the 12-year analyzed period. It is interesting to note that irrigated lands, pastures and urban areas (i.e., artificial surfaces) increase while natural vegetation, unknown/not irrigated land and water bodies decrease. Such detected land cover changes are usually associated to an increase of demand of water supply which main source, in this area, is represented by groundwater [5]. Indeed, the increasing of crops (e.g., sugar beet and corn) translates into more irrigation needs as well as the growth of urban areas (i.e., of population and industry) and pastures (i.e., of breeding) leads to higher water consumption. Land changes that may be correlated to groundwater overexploitation are also analyzed by using a Landsat 5 TM image acquired on May 2000 and a Landsat 8 image taken on May 2015 (Figure 7A,B). Visual inspection of the two images reveals a clear pattern of crops in 2015 which is not visible in 2000, pointing out an increase of land covered by intensive agriculture and irrigated fields over the 15-year period of analysis. Such an observation is not straightforward for the area around Hotamis, probably due to a change in cultivation type that is not visible in the 2015 image. Moreover, a supervised classification of SAGA GIS based on the maximum likelihood method (Figure 7E,F) was applied. Landsat images were classified into four land cover types: (1) water bodies; (2) irrigated land; (3) not irrigated land/natural vegetation/pastures; (4) rock/bare soil/salt. Classification results, within the sample area, are in agreement with CORINE data, showing that irrigated lands increased from 7% to 14% while not-irrigated lands decreased from 26% to 11% (

Table 8. Land cover changes computed by using classified Landsat images (see Figure 7).

Land Cover Classes	Area (km2)
	2000	2015
Irrigated land	760	1634
Non- irrigated land/Natural vegetation/Pasture	3003	1271
Rock/Bare soils/Salt	7538	8389
Water bodies	20	6
Not classified	59	81
Total Area	11,380	11,380

).

We also computed the mNDVI for May 2000 and May 2015 which shows a general increase of vegetation activity in the cultivated land (Figure 7C,D).

Table 9. Average mNDVI in 2000 and 2015 for the main CLC categories.

CLC Category	mNDVI
	2000	2015	Variation
Natural vegetation	0.45	0.45	−0.001
Bare soilRock	0.28	0.22	−0.068
Pastures	0.47	0.48	0.005
Already irrigated land	0.89	0.96	0.068
New irrigated land (2012)	0.66	0.98	0.325
Non-irrigated arable land	0.52	0.74	0.216
Artificial surfaces	0.57	0.58	0.011

reports the relative variation of mNDVI for some categories of CLC. Irrigated land shows the highest mNDVI. New irrigated land shows the highest increase of mNDVI along with the non-irrigated land, this latter probably partially irrigated. The other categories of land cover do not show significant changes of mNDVI, pointing out that the vegetation activity increase is probably more related to irrigation than to climate conditions, as also reported in [5] using MODIS data.

4.3. Meteorological Analysis

The analysis of data recorded by the available meteorological stations (Figure 1, Table 4) allowed the identification of precipitation (P) and temperature (T) variations occurred in the area.

Analysis of temperature (Figure 8A), based on computing the average annual temperature data acquired by Konya station from 1950 to 2014, shows an increase of 0.75 °C during the observation period, with an higher rate starting from 1990. This increase is in agreement with the general global warming and with findings from climate-related studies [2,58,59]. Such a temperature increase may contribute to the rise of water needs for irrigation.

Analysis of precipitation, carried out by computing 1-year as well as 3-year cumulative rainfall data collected by Konya station from 1935 to 2014, did not reveal any clear trend (Figure 8B). However, it is possible to observe cycles of dry and wet periods, probably related to NAO index [60]. In particular, by analyzing the 3-year cumulative rainfall graph, it is possible to observe two minimum rainfall periods, in the early 1990s and in the 2004–2008 period, which partially overlap with the ENVISAT observation period (2002–2010).

Moreover, we performed an analysis on precipitation by averaging data acquired by all the five available stations during the 2002–2010 period (Figure 8C). The five stations show similar climate conditions from 2000 to 2010 with a mean annual rainfall ranging from 283 mm of Karapinar to 330 mm of Cumra. The analysis of 3-months cumulative precipitation data pointed out a typical seasonal trend, with less than 25 mm in July–September and up to 120 mm in November–January, while the 1-year cumulative precipitation data analysis revealed a minimum rainfall period from 2004 to 2008, followed by a wet period in 2009–2010.

5. Discussion

In this Section, we analyze the SBAS-DInSAR results in relation to the information derived from independent (remote sensed and ground-based) data, in order to understand the role of different factors (i.e., type and thickness of lithology, groundwater depletion induced by natural causes and human activities) in predisposing and/or triggering land subsidence (Figure 9).

5.1. Correlation between Land Subsidence and Lithology

As already mentioned in Section 2, alluvial sediments, mainly made up of clay and silt interbedded by sandstone and gravel, are widely present in the Konya plain. Such clay and silt units are characterized by high compressibility, playing as predisposing factor in the processes of subsidence induced by groundwater extraction.

Due to the number and distribution of stratigraphic cores within the plain, it was not possible to quantitatively investigate the correlation between subsidence and lithology. Only six of the available cores are characterized by a significant thickness of compressible layers (which reaches up to 130 m) and, therefore, the regression analysis carried out to estimate the relationship between subsidence rate and thickness of such cores is not meaningful. However, the joint analysis of Figure 9A,B allows a qualitative estimate of the relationship between the spatial distribution of DInSAR derived surface deformations and the thickness of compressible layers.

Such relationship is better pointed out in Figure 10 showing a simplified geological section against the vertical deformation rate along a profile extending from Hotamis to Karapinar (for location see Figure 9). A spatial correlation between land subsidence and alluvial deposits can be observed. Moreover, fastest subsidence affects the areas with thickest sediments as shown in

Table 10. Vertical deformation rate in relation to thickness of compressible layers.

Compressible Layer Thickness (m)	Number of Wells	Vertical Deformation Rate (cm/Year)
		Min	Max	Mean
Obruk plateau	3	0.03	0.25	0.15
<10	2	0.06	0.47	0.26
50–100	2	0.30	0.56	0.41
>100	2	0.43	0.49	0.46

, which reports the vertical deformation rate in relation to the thickness of compressible units.

5.2. Correlation between Land Subsidence and Groundwater Depletion

Groundwater level variations in the study area have been investigated by analyzing data acquired by a local network of piezometric wells. Figure 9C shows the cumulative water table depletion during the period 2003–2010, estimated from wells located in the subsidence area. Wells recorded a significant decrease of the piezometric level up to about 14 m, except for well 47,761 located at south of Karapinar, which measured a lower decrease (about 2 m).

Analysis of time-series of piezometric data reveals a seasonal variation reaching a maximum piezometric level approximately at the end of winter season (end of March) and a minimum piezometric level in October (Figure 11A–D). Such a variation, which can be related to climatic conditions and irrigation cycles, is observed in all the wells located along the plain expect for, as aforementioned, the well 47761 which shows not significant groundwater changes, probably due to the presence of a local shallow aquifer (Figure 11E).

As regards the general trend, it shows a maximum decrease of piezometric level from 2004 to 2008, which corresponds to a minimum rainfall period (Figure 8), as pointed out from the meteorological data analysis (see Section 4.3). Subsequently, despite the 2009–2010 time interval was affected by heavy precipitation events (Figure 8), the groundwater table level remains stationary or shows a little increase, due to a condition of aquifer overexploitation.

Groundwater level changes occurred at regional scale have been analyzed by using GRACE satellite data. The water storage anomaly estimated by GRACE (Figure 11F) shows a trend in general agreement with wells data. In particular, a seasonal trend (about 150 mm of storage water anomaly oscillation), probably also related to surface water (reservoir), is present. This is superimposed to a general trend showing a decrease of storage level until 2008 (about −100 mm of water storage) followed by a weak increase in the 2009–2010 period. It is important to remark that GRACE observations are characterized by a footprint larger than our study area and may be influenced by various processes beyond the groundwater level variations (see Section 3.5), thus cannot be used for a direct comparison with the piezometric measurements available for our study area. However, they can be profitably exploited for a qualitative assessment of the water resources depletion occurring at regional scale.

As regards the piezometric measurements, we noticed that a groundwater level decrease is measured in the wells located along the Konya plain and on the Obruk plateau as well. However, such a groundwater drawdown had a different impact on the surface in terms of deformation.

To investigate the susceptibility of these two geological settings (i.e., plain and plateau) to ground instability, we performed a joint analysis of DInSAR deformation time-series and piezometric data time-series for four wells, two located in the alluvial Konya plain (wells 52268 and 47761) and two on the limestone Obruk plateau (wells 42278 and 52258; Figure 9 and Figure 12). We remark that, despite the uncertainties related to the few available data which hampered us to carry out a quantitative correlation analysis between land subsidence and groundwater level changes, our analysis allowed pointing out the concomitant effect of lithological characteristics and groundwater depletion on surface displacements.

The well located on the limestone plateau (Figure 12A) and the well located on the boundary between plateau and plain (Figure 12B) recorded a significant decrease of the piezometric level reaching about 2 m and 12 m, respectively. However, the deformation times-series associated to SAR points located near these two wells show a quasi-stable trend with a subsidence rate lower than 0.2 cm/year. Therefore, such result shows no correlation between groundwater drawdown and subsidence process in the area where limestone crops out.

On the other hand, as expected, in the plain area where soft soils are present, a correlation between piezometric level changes and land subsidence is observed. This is particularly evident by comparing time-series of DInSAR deformations and piezometric data acquired by the 52268 well (Figure 12C). For this well, the groundwater level decreases of about 12 m during the 2003–2010 period and the times-series of near SAR points reveal that, in this time interval, a cumulated displacement of about 5 cm occurs in the compressible layer of alluvial sediments. The high deformability of such materials is confirmed by the DInSAR time-series of SAR points located in the proximity of well 47761 (Figure 12D). Indeed, such well, probably due to local differences of the aquifer characteristics, recorded a low piezometric level decrease (<2 m) but it is located in an area affected by moderate subsidence (0.4 cm/year).

5.3. Karst Processes and Land Subsidence

Although the groundwater drawdown has not caused detectable progressive subsidence where limestone crops out, subsidence is also present in the latter area in the form of karst collapses. As mentioned above, these processes represent a special hazard in the study area and they are also linked with the decline of the groundwater levels [9]. Unfortunately, we have not detected displacements associated to dissolution-induced subsidence, even though the time span covered by SAR measurements (2002–2010) has coincided with a period of high collapse occurrence (13 collapses between 2006 and 2009). Castañeda et al. [61] were able to detect karst subsidence in the Ebro Valley (Spain) with a similar resolution (90 × 90 m²) and using the same DInSAR processing technique. However, the karst processes described in their study area are widely different from those observed in the Konya basin. The Ebro Valley shows an evaporite karst where large depressions (>50,000 m²) are developed because of a progressive subsidence associated to the dissolution of gypsum, anhidrite, and most specially, glauberite and halite. These highly soluble minerals are responsible of a rapid dissolution phenomenon that accelerates the subsidence and influences its displacement regime (continuous, episodic or mixed). This is not the case of the Konya basin where a limestone karst has been developed and the observed landforms indicate that the processes that control the karst subsidence seem to have an episodic regime (i.e., collapse). The large sinkhole collapses (the “old” obruks) and their surroundings appeared to be inactive as our DInSAR results indicate. On the other hand, the spatial resolution of our displacement maps is a priori inadequate for detecting movements associated to the “young” obruks. Baer et al. [62] based on ERS data of similar resolution failed to detect sinkhole induced deformation along the Dead Sea shores, an area widely affected by karst subsidence processes. Galve et al. [26] show that the minimum detectable sinkhole area in their test site, using a technique that provides 369.8 points/km², was approximately 2500 m². Our analysis provides 17–21 measure points/km² and the extent of the largest “young” obruk (<500 m²) is also below the mentioned minimum detectable dimensions. Nof et al. [63] demonstrate that small sinkhole induced deformation or collapse precursory displacements can be recognized with DInSAR techniques using high resolution data such as the SAR images provided by the COSMO-SkyMed satellite (pixel size of 3 m). Moreover, a higher spatial resolution is needed for recognizing sinkhole collapse precursory displacements. Nevertheless, the displacement maps produced in this research could be used for assessing the sinkhole hazard in the study area. As the “young” obruks are related to the lowering of the water table, maps of groundwater drawdown could be derived from our displacement maps and the latter could be used for assessing the areas with the highest sinkhole hazard.

6. Conclusions

Combining remote-sensed and in situ data and information, we explored the relationship between land subsidence and its main controlling factors in one of the most vulnerable regions of Turkey, the Konya plain.

Through space-borne SBAS analyses, we detected a wide subsidence phenomenon extending for about 1200 km² in the plain, with a deformation rate up to 1.5 cm/year in the 2002–2010 period. The SBAS technique allowed achieving a significant density of measure points, demonstrating the capability of this DInSAR approach to effectively investigate displacements also in semi-rural settings. Moreover, the availability of ascending and descending SAR datasets allowed us to investigate the direction of the ground movement and retrieve vertical displacements reaching up to 10 cm in the observed time interval.

The specific role of the geological and hydrogeological conditions has been explored, revealing that the consolidation of compressible alluvial sediments associated to the groundwater level decrease drives the subsidence and controls its temporal and spatial evolution over the plain. Groundwater lowering is also one of the factors that trigger sinkholes in the Obruk plateau causing damage to infrastructures and villages.

Groundwater depletion is mainly caused by overexploitation due to uncontrolled irrigation and unsustainable agricultural practices, as pointed out by our preliminary land cover change analysis showing an increase of irrigation fields of about 120 km² during the 2000–2012 period. The analysis of data acquired by piezometric wells revealed a groundwater level decrease up to about 14 m in correspondence with the minimum precipitation period (from 2004 to 2008) detected by analyzing meteorological data. Moreover, the analysis of climatic data shows an increase of temperature of 0.75 °C in the last sixty-five years, with an higher rate starting from early 1990s; such an increase probably further contributed to the rise of water needs for irrigation.

The use of GRACE data allowed us to point out that the groundwater depletion is not limited to our study area but affects a wider region in the Anatolian Plateau, demonstrating that such satellite data can be profitably exploited for assessing water resources at regional scale.

Our study has relevant local and global implications. Locally, it provides a holistic view over the anthropogenic and environmental dynamics in the Konya plain, encouraging the development of an integrated monitoring system which may play a relevant role for the protection of the area. At a more global scale, our approach may be extended to investigate land subsidence induced by groundwater overexploitation in different climatic and geological settings worldwide.