Modeling of land subsidence caused by groundwater withdrawal in Konya Closed Basin, Turkey

Land subsidence is a threat that occurs worldwide as a result of the withdrawal of fluid and also underground mining. The subsidence is mainly due to excessive groundwater withdrawal from certain types of rocks, such as fine-grained sediments. Mitigating the effects of land subsidence generally requires careful observations of the temporal change in groundwater level and ideally modeling of groundwater flow and subsidence. In Turkey, land subsidence is a crucial issue in the Konya Closed Basin. When simulating the effect of long-term groundwater withdrawal on the spatial variation of subsidence rates, various coupled numerical groundwater-flow and subsidence models have been used. Also, GPS, InSAR and ENVISAT SAR images have been used for verification of the models' parameters. In the work reported here, a novel numerical solution based on consolidation theory was developed in MATLAB to predict the land subsidence of the Konya Closed Basin. In order to adjust the model to the local conditions, historical data from the study area for the years 2011–2014 were used. The presented solution allowed for subsidence model development which can support the prediction of the ground movement for the Konya Closed Basin in Turkey.

1 Introduction

Land subsidence occurs when groundwater has been withdrawn from certain types of rocks such as fine-grained sediments. The rock compacts due to the pressure decrease in the aquifer strata. Eliminating the adverse effects of land subsidence requires observations of the temporal change in elevation coupled with groundwater flow modeling (Özyurt et al., 2018). In the areas where excessive groundwater is withdrawn from susceptible aquifer systems, land subsidence induced by aquifer-system compaction is typically observed, resulting in severe socio-economic damage for the affected communities. Generally, it is the lowered groundwater levels caused by groundwater pumping that cause land subsidence in compressible aquifer systems. Such systems typically are made up of basin-fill and unconsolidated alluvial layers that include both aquifers and aquitards (Galloway and Burbey, 2011). Either slow or sudden, the motion of the ground due to subsidence is a hazard that threatens civil and environmental resources. Present and potential future hazards have been assessed by models that are based on basic relations between groundwater levels or aquifer hydraulic head, effective or intergranular stress, the compressibility of groundwater and the aquifer skeleton, and groundwater flow (Özyurt et al., 2018). These models use two different approaches: the first is based on groundwater flow theory (Jacob, 1940) and, secondly, the theory of linear poroelasticity (Biot, 1941). Various models implementing, e.g., the De Glee method and software programs like IMOD have been used to simulate land subsidence (Guzy et al., 2018). Land subsidence is a consequence of the over-exploitation of groundwater resources from susceptible aquifer systems. According to an assessment, measured subsidence rates in different locations in the world range from 6 mm yr⁻¹ in Kolkata, India, between 1992 and 1998 to 300 mm yr⁻¹ in Mexico City, Mexico, between 2004 and 2006 (Osmanoğlu et al., 2011). And among 18 selected sites distributed globally the mean and median subsidence rate is 100 and 55 mm yr⁻¹, respectively (Özyurt et al., 2018).

In Turkey, land subsidence is prominent in the Konya Plain. The MODFLOW numerical groundwater flow model coupled with the subsidence (SUB) package (Hoffman et al., 2003) has been used to simulate the effect of long-term groundwater abstraction on the spatial variation of subsidence rates (Özyurt et al., 2018). Also, Caló et al. (2017) used the Small BAseline Subset DInSAR technique to process two datasets of ENVISAT SAR images spanning from the 2002–2010 period. The produced ground deformation maps and associated time series allowed the detection of widespread land subsidence extending for about 1200 km² and measured vertical displacements of as much as 10 cm in the observed time interval. This area, which is under strong anthropogenic pressures and is confronted with serious water-related problems, is widely affected by land subsidence (Caló et al., 2017). For that reason, approaches have been taken to investigate this fragile environment of the Konya Plain, located in central Turkey to help understand and predict the occurrence of land subsidence in order to quickly remediate it (Özyurt et al., 2018). Groundwater levels in the Konya Closed Basin have been declining at a rate of about 1 m yr⁻¹ since the early 1980s. Recent GPS observations reveal subsidence rates of 22 mm yr⁻¹ over the southern part of the Konya Closed Basin (KCB). The MODFLOW numerical groundwater flow model coupled with the SUB package was used to simulate the effect of long-term groundwater abstraction on the spatial variation of subsidence rates (Özyurt et al., 2018). This research reported here seeks to use MATLAB software to model and predict the land subsidence of the KCB. MATLAB predictive analytics uses historical data to predict future events. Typically, historical data are used to build a mathematical model that captures important trends. That predictive model is then used on current data to predict what will happen next or to suggest actions to take for optimal outcomes. Predictive analytics has received a lot of attention in recent years due to advances in supporting technology, particularly in the areas of big data and machine learning.

2 Study area

KCB is located in the central part of Turkey between 36^∘51^′ N and 39^∘29^′ N and 31^∘36^′ E and 34^∘52^′ E. It is recognized for its agricultural and industrial productivity as well as its geomorphological and hydrogeological characteristics (Üstün et al., 2013). The KCB, which encompasses a 55 000 km² area of the Anatolian peninsula, is the biggest closed basin of Turkey involving Konya, Karaman, Nigde and Aksaray provinces (Canaslan Comut et al., 2015). The Taurus Mountains enclose the basin from the south and southwest (Üstün et al., 2015). Because the mountains around the KCB prevent water from reaching the Mediterranean Sea from the basin, it is a closed basin with unique characteristics (Üstün et al., 2014). The elevations vary from 850 to 1000 m over the flat plains and rise to 3900 m on the Taurus Mountains. From a morphological point of view, the KCB is an orogenic plateau located to the north of the Taurus Mountains Belt extending along the eastern Mediterranean Sea. The plains of the southern and northern subbasins are located roughly at 1100 and 950 m, respectively. The elevation along the water divide between the two subbasins declines from about 2200 m in the Taurus Mountains to about 1100 m in the most northern part (Özyurt et al., 2018). The amount of water used for agricultural irrigation can be as much as 70 % of the water use for the whole country. The basin has a significant groundwater-resource potential, but because of small amounts of annual rainfall in normal years, susceptibility to drought, and overuse of the groundwater resource, groundwater levels are decreasing (Canaslan Comut et al., 2015). According to the records of the Turkish State Meteorological Service, the average annual rainfall in the KCB is 287 mm (Yavuz, 2010). Due to the overuse of groundwater, subsidence occurs in the area surrounding the city of Konya (a current problem) and also in the southwestern part of the basin (Canaslan Comut et al., 2015). A simplified scheme of the hydrogeological system includes a confined Neogene aquifer and the low-permeability Plio-Quaternary lake sediments on top. The Neogene aquifer includes lacustrine carbonates in which extreme secondary porosity and permeability develop due to karstification and low-permeability interbeds rich in clay and silt. Groundwater obtained from the Neogene aquifer is the principal source of water in the Konya subbasin, where extensive agricultural production has been carried out since the 1970s. Because annual groundwater abstraction is about 2 times the annual natural recharge, regional groundwater levels have declined at a rate of about 1 m yr⁻¹ since the 1980s. The land subsidence accompanying groundwater over-exploitation is mostly subtle, but sudden collapses have also occurred with an increasing annual frequency during the last decade(s). The meter/decameter-sized collapse structures threaten life and property (Fig. 1) (Canaslan Comut et al., 2015).

Figure 1 An earth fissure caused by land subsidence in the southwestern part of the KCB.

3 Methodology

The research methodology consisted of eight stages (Fig. 2).

The aim of the first three stages was to calculate land subsidence for years 2011, 2012, 2013 and 2014 at particular test sites. For that purpose, Terzaghi's law of soil consolidation was used as shown in Eq. (1) (Terzaghi et al., 1967):

$$\begin{array}{}\text{(1)}& S={\displaystyle \frac{\mathrm{\Delta }e}{e_{\mathrm{i}}+\mathrm{1}}},\end{array}$$

where S is land subsidence; Δe is change in initial void ratio; e_i is initial void ratio; and D is thickness of the aquifer.

The changes in void ratio and initial void ratio, Δe and e_i, respectively, were determined during laboratory tests for 42 wells distributed in the study area and obtained from Ali Ferat Bayram of Konya Technical University, Faculty of Engineering and Natural Sciences. However, due to the lack of full recognition of geological conditions in the research area, it was assumed that the entire aquifer is compressed. Therefore, D was assumed to be 650 m (Özyurt et al., 2018).

In the fourth stage of the research, previously calculated values of land subsidence for the years 2011–2014 were used in predicting its future occurrence. Therefore, a model for the known period under review was developed using Eq. (2) (Koopejan, 1948):

$$\begin{array}{}\text{(2)}& S=\left(\left({\displaystyle \frac{\mathrm{1}}{\mathrm{Cp}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{Cs}}}\times \log t\right)\ln \left({\displaystyle \frac{p_{\mathrm{i}}+\mathrm{\Delta }{p}_{\mathrm{i}}}{p_i}}\right)\right)D,\end{array}$$

where Cp and Cs are consolidation constants; t is time; p_i is initial pressure; and Δp_i is change in pore pressure.

Because values for Δe and e_i were available in the data acquired for the study area, a relationship between p_i and Δp_i (Eq. 2) and e_i and Δe was used to formulate Eq. (3):

$$\begin{array}{}\text{(3)}& S=\left(\left({\displaystyle \frac{\mathrm{1}}{\mathrm{Cp}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{Cs}}}\times \log t\right)\times \left({\displaystyle \frac{\mathrm{\Delta }e\times \mathrm{Cc}}{{\left(e_{\mathrm{i}}+\mathrm{1}\right)}^{\mathrm{2}}}}\right)\right)D,\end{array}$$

where Cc is the compression index.

Figure 2 Flow chart of the research methodology.

Having developed the equation for subsidence as a function of time in the fifth stage of the research, the minimum and maximum subsidence values calculated for 42 wells for the period 2011–2014 were used in developing two simultaneous equations. This was done to determine Cp and Cs. However, Cc could be determined using Eq. (4) or (5) (Habibbeygi et al., 2017):

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathrm{Cc}=\mathrm{0.29}\times \left(e_{\mathrm{i}}-\mathrm{0.27}\right),\text{(5)}& {\displaystyle }\mathrm{Cc}=\mathrm{0.43}\times \left(e_{\mathrm{i}}-\mathrm{0.25}\right),\end{array}$$

For this research purpose, e_i values of the corresponding minimum and maximum subsidence values were used in both Eqs. (4) and (5). Therefore, four values of Cc were computed. After Cc was determined, Cp and Cs were then determined based on Eq. (3). These calculations were carried out using sets of equations with geotechnical properties associated with minimum subsidence for one equation and those properties associated with maximum subsidence for the other equation as determined based on Eqs. (4) and (5). Finally, with t expressed in days, a subsidence prediction for a period of 30 years was carried out based on Eq. (3). In the prediction, values obtained for the forecasted land subsidence were not realistic. For that reason, further research highlighted calculations on consolidation constants. In research by De Glopper and Ritzema (2006), the direct effect stands for the subsidence that occurs the first day after the load increase and is based on the assumption that Cp and Cs are independent of the size of the load, but where Cp depends on the selected time period (1 d). The secular effect stands for that part of the subsidence that occurs as the excess water is drained out of the aquifer-system material. This is a very slow process, especially in clays, because of their low hydraulic conductivity. The secular effect will cause the subsidence to continue indefinitely. Therefore, using Cp =1 and applying both equations of Cc (Eqs. 4 and 5), as well as in the simultaneous equations of minimum and maximum subsidence as a function of time, Cs =0.292 and Cs =0.471 were computed based on the minimum subsidence obtained for both equations of Cc. Also, Cs =0.3347 and Cs =0.5301 were computed based on the maximum subsidence obtained for both equations of Cc. After all parameters in the subsidence as a function of the time equation were known, subsidence values were then predicted using the individual geotechnical parameters over a period of 30 years with reference to the years under review.

4 Discussion

From the MATLAB model for estimating the vertical consolidation which is the predicted subsidence, both Eqs. (4) and (5) for the compression index Cc were applied, and the results obtained from both equations as shown in Fig. 3 had no significant difference between the two compression indexes used. This was confirmed from the dependent t-test on the minimum and maximum subsidences using both equations of Cc. The p values of the minimum and maximum subsidences in all cases were 0.4 (> 0.05) and as such indicated no significant difference between the two equations of Cc used. The minimum and maximum subsidences to be observed were plotted and a general regression analysis was done. Linear regression, exponential regression, logarithmic regression, polynomial regression, and power regression were used. Logarithmic regressions resulted in R² values of 1 in both cases of Cc for minimum and maximum subsidences, whereas the other regression analyses had R² values less than 1. Therefore, in calculating the minimum level of subsidence, the model generated for a logarithmic regression was either Eq. (6) using Eq. (4) as shown in Fig. 3a or Eq. (7) using Eq. (5) as shown in Fig. 3b.

$$\begin{array}{}\text{(6)}& {\displaystyle }S\left(t\right)=\mathrm{0.4787}\ln \left(t\right)+\mathrm{3.1460}\text{(7)}& {\displaystyle }S\left(t\right)=\mathrm{0.4550}\ln \left(t\right)+\mathrm{3.1781}\end{array}$$

In calculating the maximum level of subsidence, the model generated for a logarithmic regression was either Eq. (8) using Eq. (4) as shown in Fig. 3c or Eq. (9) using Eq. (5) as shown in Fig. 3d.

$$\begin{array}{}\text{(8)}& {\displaystyle }S\left(t\right)=\mathrm{2.0145}\ln \left(t\right)+\mathrm{13.4380}\text{(9)}& {\displaystyle }S\left(t\right)=\mathrm{1.9079}\ln \left(t\right)+\mathrm{13.5850}\end{array}$$

Figure 3 Graphical representation of vertical consolidation against time in years.

The subsidence results obtained were compared with related land subsidence from previous research (Üstün et al., 2013, 2015) for the period 2011–2014. In that regard, land subsidence results were in the range of 0.2–4.0 cm yr⁻¹ using GPS and 1.0–4.0 cm yr⁻¹ using InSAR data analysis. Thus, the maximum subsidence can be assumed to be the maximum range for the 4 years (2011–2014), which is 16.0 cm. Moreover, data from ENVISAT ASAR images with ENVISAT ASAR data between 2003 and 2010 above the city center, compared with GPS measurements, estimated the deformation rate of KCKY GPS station to be approximately 2.0 cm yr⁻¹ (Canaslan Comut et al., 2015). Also, considering the 8-year period of that study gives a maximum vertical deformation or consolidation of 16.0 cm. The above values are confirmed in this research with the use of Terzaghi's 1D consolidation equation and the aid of the MATLAB model which predict land subsidence ranging from 3.8 to 16.2 cm. This was analyzed for the same period of 4 years, which results in a subsidence rate in the range of 1.0 to 4.1 cm yr⁻¹ approximately in conformity with past research of 1.0 to 4.0 cm yr⁻¹.

5 Conclusions

The research focused on Terzaghi soil consolidation theory and subsidence as a function of time with MATLAB analysis. With the subsidence, compression index and consolidation constants predicted, for the study area, the value for minimum subsidence or vertical deformation 30 years after the year 2011 is predicted to be approximately 5.0 cm, whereas the maximum land subsidence is predicted to be 20.0 cm assuming the groundwater withdrawal is the same as the rate of water recharge into the aquifer. However, it is highly recommended that further accurate data concerning porosity, moisture content, compression index and effective stress of the field be obtained for more accurate subsidence results for a 3D consideration.

Finally, yet importantly, consolidation constants should also be measured and identified from the laboratory. Results can then be used to predict the trend of land subsidence with a given lifespan of the study area.