Estimation of land subsidence using coupled particle swarm optimization and genetic algorithm: The case of Damghan aquifer

Land subsidence refers to the vertical and downward movement of the land's surface. Various factors such as human activity, natural processes, or both can cause subsidence phenomena. The subsidence caused by natural processes often happens gradually and in the long term and does not leave many harmful effects. However, the subsidence caused by human activities happens suddenly and causes harmful effects (Goorabi et al. 2020; Li et al. 2022).

Over the past few years, the world's population has increased the demand for water in agriculture, domestic, and industrial sectors. This has increased groundwater consumption, especially in areas with arid or semi-arid climates. However, to prevent the overuse of this resource, it is important to explore other sources of water, such as wastewater (Emamgholizadeh et al. 2014; Saraiva et al. 2020; Ingrao et al. 2023; Levintal et al. 2023). Therefore, it is necessary to use other water resources such as wastewater. It is worth noting that the discharge of wastewater (i.e., brine) degrades groundwater quality and thus water cannot be directly used for potable water (via desalination) and industrial applications (Panagopoulos 2020a, 2020b, 2022).

Over-pumping groundwater due to rapid population growth, urbanization, and industrialization has increased global concern. Over the past 50 years, a causal relationship existed between the over-pumping of groundwater and many hydrogeological hazards (Gorelick & Zheng 2015). Between 1950 and 1970, in parallel with industrialization and population growth in the world, the risk of land subsidence due to groundwater depletion was reported (Karemi et al. 2013).

In Iran and some countries, the arid and semi-arid climatic conditions have led to the overexploitation of groundwater resources and, thus, subsidence occurrence (Sharifikia 2011; Konikow 2015). The first subsidence report in Iran dates back to 1967 in the Rafsanjan Plain. In 30 years (1969–1999), the groundwater table in this plain has dropped about 25 m and caused subsidence equivalent to 15 cm (Tourani et al. 2018). The study result of Goorabi et al. (2020) shows that in Isfahan, the amount of land subsidence is estimated at a rate of −5 to −100 mm/year.

Since land subsidence is very complex and there are different intensities and domains of land subsidence at various times and places (Mohammady et al. 2019), various aspects (e.g., monitoring of ground motions) are addressed in research on subsidence (Amelung et al. 1999; Schmidt & Bürgmann 2003; Galloway & Hoffmann 2007; Higgins et al. 2013; Lo et al. 2022). The increasing risk of land subsidence induced by groundwater extraction demands the development of efficient optimization tools to manage groundwater usage (Gorelick & Zheng 2015; Bagheri-Gavkosh et al. 2021). Generally, subsidence is divided into three major types: The first one includes natural and local subsidence induced by the empty spaces underground, karst or pseudo-karst regions, and mining (Jalini et al. 2018). The second one occurs elastically in coarse-grained soils, and its extent depends on the amount of applied load and the soil type (Moarefvand & Shamsadin Saeid 2013). Finally, the third type occurs plastically in fine-grained soils due to groundwater depletion as a gradual and irreversible process (Tomás et al. 2005; Conway 2016). Subsidence can lead to uneven changes in elevation and slope of streams, canals, water distribution structures, well wall pipes failure and protrusion, compressive stress caused by compression of aquifers, disruption in groundwater exploitation, and irreversible depletion of all or part of groundwater reservoirs. This phenomenon can also destroy the vital arteries of essential structures and infrastructures. In this regard, identifying the plains with the potential for subsidence and estimating their amount can certainly play a significant role in managing and controlling this phenomenon.

Monitoring and mapping the subsidence vulnerability of the region is the first step to studying land subsidence. In the present study, the subsidence vulnerability index (SVI) of the Damghan Plain is evaluated using the ALPRIFT model (Nadiri et al. 2018, 2020, 2022). This model was developed by Nadiri et al. (2018) to evaluate the SVI based on seven factors affecting the land subsidence, i.e., the types of subsoil, land use, pumping, recharge, the thickness of the plain aquifer, distance to the fault, and groundwater depletion. In the next step, the measured field data in the study area are used to assess the ALPRIFT model. In this research, seven factors affecting land subsidence, including the aquifer media (A), land use (L), pumping (P), recharge (R), aquifer thickness impacts (I), distance from the fault (F), and depletion in the water table (T) were used as input data to determine the amount of vulnerability. These factors were proposed by Nadriri et al. (2018, 2020, 2022) to evaluate land subsidence. Hence, combining the input data can play a significant role in obtaining the exact value of subsidence vulnerability. To the best of our knowledge, no research has been conducted using the particle swarm optimization-genetic algorithm (PSO-GA) to estimate land subsidence. Thus, this model was used in this study to determine the optimal input coefficient of the ALPRIFT model. The capability of the proposed PSO-GA was also examined by comparing its results with those of the GA and gene expression programming (GEP) models.

According to the data obtained from drilling logs, exploratory wells, observation wells, and exploitation wells, alluvial deposits in the entire study area are formed from a wide range of fine and coarse grains. Also, evaporites (i.e., gypsum and salt) are observed in some parts of the plain. There are wide fluvial deposits on the southern foothills of the Alborz Mountains. The alluvial deposits of the Damghan Plain contain the following sediments:

• Alluvial fans sediments of foothill and plain,

• Fine-grained sediments of the plain, and

• Desert sediments (i.e., clay and salt marsh, saline and moist clay, clay, and salt).

Damghan Plain aquifer consists of the following five geological units:

• Unit A (the northern highlands region including limestone formations, limestone marl, lime-dolomite, sandstone, and flysch),

• Unit B (the southern part of the aquifer, including lime-marl formations, clay, and siltstone-marl, along with a significant amount of evaporite, salt, and chalk),

• Unit C (the southwest part of the aquifer containing impenetrable volcanic formations),

• Unit D (the eastern part of the aquifer containing dolomite, limestone, marl, sand, and mudstone formations), and

• Unit E (the western part of the aquifer including the alluvial terrace, young alluvial deposits, limestone, and sand).

In the study area, population growth has increased in recent years, followed by increasing agricultural and industrial developments, leading to the formation of economic, social, and trade activities. Furthermore, the number of wells has increased in recent years, and traditional irrigation with low efficiency is used for agricultural purposes. Then, the extraction and exploitation of groundwater resources exceeded the potential capacity of the exploited plain. Consequently, the amount of water level decline and land subsidence has constantly increased in this area, and land subsidence in the Damghan Plain is now growing widely.

Groundwater balance calculations for the Damghan Plain in the 17-year period (2006–2022) indicate that the amount of net recharge is equal to 11,000,000 m³, and the reservoir deficit in 2006 and 2022 was 21,000,000 and 31,000,000 m³, respectively.

The transmissivity of the plain's alluvial aquifer is estimated between 40 and 1,200 m²/day, and changes in the storage coefficient of the aquifer are estimated between 2 and 9%. Based on the inverse distance weighted (IDW) method, the average transmissivity and storage coefficient of the aquifer for the whole plain was estimated as 185.1 m²/day and 5.1%, respectively.

Based on the logs of wells drilled in the study area, the alluvium thickness of the Damghan Plain in the northern and northeastern parts of the plain is more than 63 m and about 74 m, respectively. Moreover, due to the bedrock uplift, this thickness reaches 27, 71, and 73 m in some of the northeast, southwest, and southeast parts of the plain, respectively.

In this study, to estimate the land subsidence in the study area, seven factors that affect land subsidence are used as input data for artificial intelligence (AI) models. These parameters are the aquifer media (A), land use (L), aquifer pumping (P), aquifer recharge (R), aquifer thickness impact (I), distance from the fault (F), and groundwater depletion in the water table (T) (see Table 1). These seven parameters were introduced by Nadiri et al. (2018) and used as input data in the ALPRIFT model to obtain the potential SVI of the aquifer. Also, these parameters were used by different researchers to calculate SVI and land subsidence (Budhu & Adiyaman 2013; Manafiazar et al. 2019; Nadiri et al. 2020; Sadeghfam et al. 2020). The main reason why these parameters have been used as an influencing factor in the amount of land subsidence is described in Table 1. As described by Nadiri et al. (2018), each of these seven parameters is divided into sub-categories, which are allocated a rank from 1 to 10 based on the degree of influence on the occurrence of the subsidence phenomenon (Table 2). A score of 10 and a score of 1 mean the most and the least effective parameter in the amount of land subsidence, respectively.

To model land subsidence, GIS software was used to prepare the input data layer for the AI models according to the explanations given in Table 2. In the following, each of these seven parameters is explained as follows.

According to the field observations and evaluations performed in the Damghan alluvial plain, three types of land use were defined: (1) agricultural areas, (2) residential areas, and (3) wastelands. The most subsidence is observed in the agricultural zones because of the excessive groundwater withdrawal in these areas. On the other hand, the lowest amount of subsidence has occurred in the wastelands. The study area was classified and scored based on the type of land use (see Table 2). Also, the study area's land use raster layer map is represented in Figure 5(b).

The annual withdrawals of exploitation wells, measured by the Semnan Regional Water Authority, were used to prepare the pumping layer. Based on the wells' coordinates and each well's discharge, the interpolation was conducted for the entire area using the GIS software using the IDW method. According to Table 2, the discharge of the wells was classified and scored based on the amount of withdrawal. The pumping layer raster map is represented in Figure 5(c).

The Piscopo method (Piscopo 2001) was used to create a recharge layer. In this method, a recharge map is obtained by combining three layers of slope, rainfall, and soil permeability. In the present study, the recharge layer was created using the three above maps, and finally, each data class was scored and classified based on the description in Table 2. A raster map of the recharge is shown in Figure 5(d).

To prepare this layer, the measured water level data of piezometric wells in 2021 were used. Also, the depth of bedrock in the study area was determined using the data measured by the Geological Survey and Mineral Exploration of Iran (2006). After determining the water table and the depth of the bedrock, the thickness of the aquifer saturation zone in the studied plain was estimated. After determining the aquifer thickness at the location of the piezometric and exploratory wells, interpolation was done for the entire area using the IDW method. The raster layer map of the aquifer thickness, which was created in the GIS software, was classified using the information provided in Table 2 (Figure 5(e)).

The distance from the fault is one of the factors affecting the subsidence and consequently determining the SVI of the plain. The greater the distance between the point and the fault the lower the risk of subsidence vulnerability. In this study, the fault map was acquired from the Geological Survey and Mineral Exploration of Iran to obtain the distance between the point of the plain and the faults in the study area. Examining the fault map in the study area indicates that the main faults are located in the north and northeast of this area. A part of the fault is also located southeast of the area. By determining the faults around the study area, the map of the distance between the plain points and the fault was calculated using Euclidean distance. According to Table 2, each point of the plain was scored based on the classification. The raster layer of the distance between the points and the fault for Damghan Plain is given in Figure 5(f).

The data of piezometric wells in the Damghan Plain for October 2021–October 2022 was used to calculate this layer. The IDW method was used to calculate the groundwater depletion in other points of the study area. Then, the raster layer of drawdown was classified and scored based on Table 2. The raster layer of groundwater depletion is illustrated in Figure 5(g).

where O_i and P_i represent the observed and predicted land subsidence, respectively, and N is the number of datasets. According to the RMSE, MAE, and R² indices, the model has high performance and accuracy in predicting land subsidence when the value of R² is close to 1 and the values of RMSE and MAE are low.

GEP is an evolutionary algorithm, and it was first introduced by Ferreira (2006). It uses the advantage of both GA and genetic programming for modeling (Ferreira 2006). It can be used to solve complex problems and for estimating the parameter under study (land subsidence). Many researchers used this model in different studies such as the prediction of groundwater table and land subsidence (Parhizkar et al. 2015; Nadiri et al. 2020) and soil science (Emamgholizadeh et al. 2017; Emamgholizadeh & Mohammadi 2021; Bazoobandi et al. 2022).

The procedure of the estimation of land subsistence using the hybrid PSO-GA can be summarized in the following steps:

In this equation, (the estimated land subsistence) is the response variable and it is the function of the predictor variables (A, L, P, R, I, F, T), and W₁, … , W₇ are the coefficients (or weights).

Step 3: The initial value was considered for W_i (i = 1–7) to estimate land subsistence (⁠⁠).

Step 4: Using Equation (8), and the initial value of W_i (step 3), the land subsidence value (⁠⁠) is estimated.

Step 5. The estimated land subsidence value (⁠⁠) is compared with the measured land subsidence value (⁠⁠), and the RMSE error is calculated using Equation (3).

Step 6. If the RMSE error is not acceptable, the hybrid PSO-GA algorithm (see the flowchart of this algorithm presented in Figure 6) was used to calculate W_i, i = 1–7, with the goal of minimizing the error, and finally, the best value of W_i is achieved.

As shown in Figure 7, the value of the SVI index is between 57 and 171, which includes three bands: Band 1 (low), Band 2 (medium), and Band 3 (high). Also, according to this figure, the highest value of the SVI is related to the northeast, central, and some parts of the south and southwest regions. This result is consistent with the subsidence values measured by the field method in the study area. In Figure 7, zones with the highest SVI are marked as regions A, B, C, and D, and the area of these regions is 19.78, 36.92, 125.92, and 30.4 km², respectively. In these areas, average subsidence of over 10 cm/year was measured. The explanation is that since there are 384 exploitation wells in these areas and 966 wells in the whole study area, 40% of the exploitation wells are concentrated in these four areas. A total of 62.9% (equivalent to 69.23 million cubic meters) of water is extracted annually from these three areas. Generally, subsidence of more than 4 mm/year is considered critical worldwide. Moreover, every 10 cm of subsidence, depending on the type of soil layers, causes about 1–2 m of the plain aquifer to be lost.

According to the evaluations carried out on the drilling logs of the well sections, which represent the condition of soil layers, the soil type in the mentioned areas (A, B, C, and D) is mainly silt and clay, which does not allow proper recharge to the aquifer due to low permeability. In addition, owing to the lack of elasticity and a high coefficient of consolidation, it undergoes irreversible consolidation after groundwater withdrawal and causes land subsidence.

The position of each model in the plot shows how closely the estimated land subsidence values were with the observation values. As the result of this figure, predictions of the hybrid PSO-GA algorithm are in agreement with the observations (R > 0.95).

According to the results of this study, out of the ALPRIFT, GA, GEP, and PSO-GA models, the PSO-GA model performed the best. In the PSO-GA model, the optimal weights for the input parameters (A, L, P, R, I, F, and T) were 1.21, 3.85, 4.82, 3.55, 1.49, 1.06, and 4.95, respectively. A comparison of the optimal weights illustrates that the most significant parameter affecting land subsidence was the combination of two factors, groundwater depletion in the water table (T) and aquifer pumping (P). Therefore, in the study area, controlling the amount of land subsidence requires managing the groundwater table and reducing pumping from the aquifer.

In the studied area, pistachios are the main cultivated crop of the plain (approximately 65.8% of the total). Our research has revealed that about 40% of this crop is irrigated traditionally using flood irrigation systems, which are found to be less than 50% efficient. Considering the relatively high water consumption of this crop, it is necessary to switch to a drip irrigation system, which has an efficiency of around 90%. This will reduce water consumption and increase the efficiency of the irrigation system. In addition, since most of the water consumption is related to the agricultural sector, it is necessary to pay serious attention to changing the cultivation pattern or improving pistachio cultivation.

The Damghan Plain is one of the most important plains of the Semnan province in Iran for agricultural use. In this plain, due to excessive extraction of groundwater to meet water needs (agriculture, drinking, and industry), the studied area is at risk of subsidence. Therefore, it is necessary to investigate and evaluate the risk of subsidence in this plain to manage and control it. In this research, the ALPRIFT model was used to investigate the vulnerability of subsidence in the Damghan Plain. The ALPRIFT model uses seven factors affecting land subsidence, i.e., the types of subsoil, land use, pumping, recharge, the thickness of the plain aquifer, distance to the fault, and groundwater depletion, to estimate the potential of land subsidence. Three AI models, including GA, GEP, and hybrid PSO-GA, were also used to optimize the weights of the seven affecting factors in the ALPRIFT model. According to the results obtained from the comparison of the mentioned models, the hybrid PSO-GA model with the lowest RMSE (1.11 cm) and the highest values of R² (0.91) outperformed the GA (with RMSE = 0.73 cm and R² = 0.73), GEP (with RMSE = 1.75 cm and R² = 0.78), and ALPRIFT models (with RMSE = 0.3.41 cm and R² = 0.64). Thus, this model can be used to optimize the weights of the proposed ALPRIFT model.

Comparing the results of the GA, GEP, ALPRIFT, and hybrid PSO-GA models showed that the hybrid PSO-GA model had better performance than others. Therefore, according to the satisfactory results of the PSO-GA model in optimizing the weights of input layers, this model can be used in the studies of other plains and obtain a correct estimation of the land subsidence value in the studied area.

Also, since the plain's central, northern, and northwestern regions are more prone to subsidence, the cultivation pattern (replacing low-consumption plants) should be changed while controlling the excessive groundwater withdrawal. Furthermore, the surface water resulting from rainfall in the upstream basin of the study area can be used to inject the groundwater table.

The authors gratefully acknowledge the support and funding provided by the Semnan Regional Water Company Research Center to conduct this research.

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.