Optimizing Reduction of Groundwater Discharge to Control Land Subsidence (Case study: Nadjafabad Aquifer)

This study evaluates and predicts the ground subsidence that happens due to the haphazard operation of groundwater resources. Also, several strategies have been developed to control this unpleasant phenomenon. For this purpose, groundwater ow simulation has been conducted using MODFLOW numerical model, and subsidence simulation in Najafabad plain has been done using SUB package under three climatic scenarios for future periods. Examination of the simulation results shows that the amount of land subsidence will increase with the aquifer operation's continuation. The maximum amount of subsidence for 6 years in drought conditions will be 23 cm at the aquifer's outlet. According to the land subsidence results at the aquifer, risk zoning of the aquifer operation was done to develop a solution to reduce the withdrawal of groundwater resources to control subsidence. Therefore, risk zoning was performed using land use and the extent of operation of groundwater resources. The results showed that the north-eastern part of the aquifer has the maximum risk of subsidence. According to the obtained results from subsidence risk zoning, scenarios of reduced water withdrawal from the aquifer in its outlet were developed. The treatment strategies results showed that the maximum amount of subsidence in wet, normal and dry conditions will be 10, 14 and 18 cm, respectively. These results indicate a 14% improvement in the quantitative condition of the aquifer in wet conditions, 10% in normal conditions and 7% in dry conditions in the total aquifer of Najafabad. Improvement of conditions by simulation shows the impact of the importance of optimal utilization of groundwater resources.


Introduction
Unplanned usage and lack of correct management of groundwater resources have led to qualitative and quantitative problems in addition to land subsidence as an irreparable damage imposed to natural resources (Li et al, 2019). In fact, land subsidence will happen due to compaction of loose sediments resulted from the reduction of water table. Non-uniform decline of groundwater in different parts of aquifer will lead to heterogeneous compaction of alluvium and non-uniform land subsidence as well as creation of cracks on the land surface (Ranjbar and Ehteshami, 2019). This gradually will lead to imposing of various damages to structures, buildings, and pipelines. The most crucial danger in this respect is the cut of pipelines such as oil, gas, fuel and water pipelines (Guzy and Malinowska, 2020).
On the other hand, sensitive structures would be prone to deformation, which will probably increase irreparable damages. In some cases, land subsidence will create a big sinkholes in the region and residents to feel insecure. According to the de nition provided by geological institute, land subsidence as a phenomenon includes collapse, or downward subsidence of the land surface can be accompanied by some minor horizontal displacements. This movement is not limited in terms of severity, and expansion of regions engaged; and, subsidence can be created based on natural geological phenomena such as liquidation, ice melting, slow movements of lithosphere, as well as exit of lava from solid crust and/or human activities such as mining, groundwater discharge or extraction of petroleum. This is mainly irreparable, costly, and destructive.
Due to groundwater discharge and water coming out of porosities, material can be compacted up to 300m in depth. The more water discharge takes place; the more compact would become the material. Land subsidence will lead to creation of deep cracks on the land surface, leaning of well casing, destruction of buildings, and siphon gamy of wells. Land subsidence is resulted from nature of the land and groundwater exit (Pardo et al., 2013). Various important land subsidence factors include height, land slope, lithological direction, distance from fault, distance from river, normalized vegetation difference index, and type of soil, amperage indicator, topographic wetness index, and land use (Pradhan et al., 2014). Groundwater resource management and provision of remedial mechanism creates possibility of controlling the phenomenon. One of the appropriate tools to study this phenomenon is simulation process, using numerical model. Mathematical model of groundwater provides simulation of a hydrological system that follows laws of physics and mathematics. Local land subsidence resulted from groundwater discharge has been numerically modeled by Mohammadpoor et al. (2016) in south-western plain of Tehran and the procedure has been predicted for up to the year 2018, based on data related to the period of 1984-2012. The modeling method used in the research has been PMWIN, being well matched with measured values. Based on numerical results of the research and assuming xed pumping rate in future; till 2018 land subsidence in south-western plain of Tehran will reach 33cm. land subsidence resulted from groundwater pumping is considered as a serious threat for south-western plain of Tehran. Land subsidence in south Spain has been simulated by Tessitore et al. (2016); and subsidence of Saveh aquifer has been modeled by Jafari et al. (2016). Simulation performed by Shen et al. (2017) has been performed via elastic model; and, subsidence resulted from gas reservoir exploitation has been simulated by Giani et al. (2018). Subsidence simulation through CLM5 model has been studied by Ekici et al. (2019). Reviewing studies performed showed that the phenomenon is a serious alarm and the role played by behavior of an aquifer is of high importance in its formation. In general, land subsidence potential all over the aquifer can affect status of groundwater discharge. So, applying numerical simulation of aquifer to predict land subsidence amount through numerical models can be suggested as one of the components of comprehensive and sustainable management of water resources (Wang et al, 2017). Present research has been performed to the aim of predicting land subsidence amount resulted from immethodical usage of groundwater in one of the arid regions is looking for remedial mechanism. The research results are considered as a warning for authorities to control and manage groundwater exploitation.

Land Subsidence
Land subsidence is a morphological phenomenon affected by the downward movement of the land. The occurrence of this phenomenon is highly commonplace; however, due to the slow movement of the land, it cannot be understood and measured most of the times. So, this event is mostly identi able in those regions where the morphology of land surface and especially facilities and equipment has been affected and resulted in destruction and some damages to happen. During the past decade and upon the prevalence of new tools such as GPS, monitoring those areas prone to land subsidence has been taken into consideration to specify land subsidence rate. This technique has been con rmed in terms of accuracy and correctness; however, problems like relatively high installation or operation cost, monitoring of permanent stations and lack of simplicity in determining subsidence domain have resulted to failure. Moreover, the unorganized monitoring time due to annual budget changes and credit amounts is another reason for the unsuccessful result. Therefore, just a few cases have been studied for years.

Study Area
The area under study is Nadjafabad, one of the main limits in Gavkhooni Basin in central part of Iran with an area equal to 1711.6km2. This area is located between 50 53' to 51 42' longitude from the east and 32 18' to 32 50' latitude from the north. The main volume of water supply resources in the region is groundwater used via wells available in the region. Groundwater depletion volume in the region is equal to 850 million cubic meters/year, which considering the aquifer replenishment sources, about 80cm groundwater decline is observed annually. The general direction of groundwater ow is from north to the south and south-eastern parts of the aquifer. Considering high amounts of unmanaged water discharge in the aquifer, it has a very critical status; and, some parts of this aquifer have cracked on the land surface. From geological point of view, the main part of the aquifer is formed by alluvium and southern and northern parts of the aquifer limited to heights are formed of limestone formations that play an important role in supply of groundwater fronts entering to the aquifer. Figure (1) shows the location of Nadjafabad Aquifer.

Aquifer Simulation
Considering the research objective, i.e., simulating land subsidence phenomenon and evaluating groundwater discharge amount in this respect; numerical ModFlow model has been used for quantitative simulation; then, SUB package has been used to determine land subsidence amount.
Simulation has been taken place using SUB package in GMS 10 interphase to use land subsidence model. The rst step for simulation is preparation of conceptual model of aquifer based on which steady-state and transient-state models of ow would be simulated. After veri cation and calibration of aquifer's quantitative model, land subsidence would be predicted nally and based on prediction of its future status via quantitative model. Accordingly, modelling is aimed at studying the following options: method of stabilizing groundwater balance, prediction of the future status of the aquifer under different hydrological scenarios, verifying conceptual model's data, and following standard mathematical modelling methods along with provision of calibration.

Preparing conceptual model of aquifer
The rst step in simulating groundwater model is to prepare aquifer's conceptual model. In fact, the model is indicative of parameters related to the water resources in the region and considered as the main input to MODFLOW. To solve groundwater equation, nite difference method has been applied. Accordingly, the aquifer has been spatially divided into a series of square or rectangular networks or cells (500x500m).
Observation wells as control points, exploitation wells, underground fronts of aquifer outlet and groundwater evaporation have been entered into conceptual model as parameters of groundwater discharge; and, owback water from discharge of wells, runoff and rainfall in ltration, as well as inlet groundwater front as parameters related to aquifer recharge have been entered into the conceptual model. After preparation of the aquifer's conceptual model, the structure of aquifer, including the land topography, initial groundwater balance, and bedrock have been entered into the model. Hydrodynamic parameter of the aquifer including hydraulic conductivity, coe cient of storage, and speci c yield of the aquifer also have been entered into the model as initial values based on pumping test results in extraction region (Kardan Moghaddam et al., 2019). In gure 2, the overall view of aquifer's conceptual model is presented.

Groundwater balance parameters of the aquifer
Considering parameters related to the groundwater balance are being introduced in the conceptual model, aquifer recharge is resulted from in ltration of runoff and rainfall, as well as inlet boundaries to it. Under this condition, in ltration from the surface and that of groundwater ows have been entered into the model respectively in expanded and GHB (general head boundary) forms. Groundwater discharge and outlet fronts include the main volume of aquifer discharge. Outlet groundwater ows from the aquifer has been considered similar to inlet ow in the form of GHB based on groundwater levels and put in the model as a known boundary. Groundwater discharge within the modelling limit only takes place from wells. The most important control points to evaluate the model are observation wells in the plain. Considering recorded statistics, 17 observation wells have been analyzed. Considering the calculations performed, the aquifer balance is highly de cient in terms of reservoirs; and, there have been 138 million cubic meters reservoir de cit as calculated in table 1.

Preparing a quantitative model of aquifer
After preparing a conceptual model to simulate the aquifer, primarily steady-state model has been used for simulation. The water ow regime in the aquifer has been assumed to be stable, and changes in water volume in the underground reservoir have been ignored. Under stable condition and despite discharge and recharge

Land Subsidence Simulation
After quantitative simulation of the aquifer, SUB package has been used to simulate the amount of land subsidence due to water discharge. SUB package simulates land subsidence modelling for those aquifers with compressible layers. This computer program has been designed to simulate vertical compression in a groundwater table; and, it simulates changes of storage coe cient, as well as compression in discontinuous or continuous interbed layers under pressure with the capability of compressibility resulted from changes of pressure in the aquifer. In this modelling, geostatic pressure is a function of groundwater balance and value of compression is a function of pressure changes effective on aquifer bottom (USGS Chapter 23 Book 6, 2997).
Groundwater exits aquifer due to reduction of hydraulic head. The volume of water discharged from aquifer is proportionate to the compressibility of soil and water; because the reduction of hydraulic head will lead to increase of effective stress in solid soil texture or skeleton and reduction of available aquifer's water pressure. Increase of effective stress on solid soil texture will lead to compaction of soil texture. SUB package calculates the volume of discharged water from aquifer storage and simulates elastic and non-elastic compaction of compressible ne-grained bed in aquifer due to discharge of groundwater.
Preconsolidation stress or head, elastic storage, non-elastic storage and initial compression have to be applied in the model as the main simulation parameters of land subsidence.
Preconsolidation head: it is the amount of minimum head with which aquifer would not be prone to subsidence. For each of the model cells that determined pre-consolidation head would be more than the initial hydraulic head, these two would be considered to have equal values. When compressible ne-grained sediments would be subjected to higher stress than maximum previous stress in the aquifer (pre-consolidation stress); compaction would be non-elastic.
Elastic storage factor: For con ned aquifer, elastic compaction or swelling of sediments is proportionate or relatively proportionate to values of the hydraulic head of aquifer. SUB package uses the following equation to calculate thickness changes of layers (positive for compaction and negative for swelling of layers).
Where Where n is porosity and nw is moisture content of soil over groundwater balance which is a fraction of total volume of porous environment.
Initial compaction: Compaction values calculated via SUB package (subsidence and aquifer system compaction) will be added to initial compaction value so that initial values would be included in compaction and land subsidence recorded values. Initial compaction will not affect the calculation of recorded or compaction changes.

Evaluating quantitative model of the aquifer
To evaluate the quantitative model of the aquifer for simulation of land subsidence amount mainly aimed at in the research; sensitivity analysis of groundwater modelling parameters has been performed, primarily. Accordingly, upon change of input parameters to the steady-state model of groundwater ow with consideration of 5 to 20% increase or decrease of the parameter; three parameters of hydraulic conductivity, the height of bedrock, and recharge of the aquifer have been identi ed as sensitive parameters. The results extracted compared to those of other numerical simulation studies are indicative of this being correct; and, the hydraulic conductivity parameter has been selected as the most sensitive parameter for calibration. Accordingly, model calibration has been performed based on RMSE (root mean square deviation); and, this value being lower than one meter has been set as calibration goal. After calibration of steady-state groundwater ow, transient state simulation has been performed for the simulation period, and speci c water yield parameter has been calibrated. Finally, veri cation of the model for the whole last year has been performed; and, the results have been indicative of proper accuracy of the results. Quantitative model's error analyses in relation to the steady-state, transient state, and veri cation have been performed, and the results have been provided in Table 2. Accordingly, preliminary parts of the aquifer with more coarse grains have higher amounts of hydraulic conductivity and middle as well as end parts of it are ne-grained with lower hydraulic conductivity. Speci c yield of the aquifer compared to hydraulic conductivity has been high and as a general procedure, the amount of speci c yield has been reduced from aquifer's inlet towards the outlet. Accordingly, in Fig. 4, the difference between the simulated and observed groundwater table is shown. According to the results, the maximum difference is observed in observation wells No. 16, 13, and 15.
Accordingly, the nal analysis of the steady-state model and transient state model of groundwater ow in aquifer after calibration is presented in Fig. 5.
Considering the complexity of groundwater system and the role played by natural factors in its behaviour measurement, veri cation of the model has to be taken place to identify an appropriate prediction model. So, considering type of aquifer and changing procedure of balance parameters of water resources in the region, veri cation has been performed for the last year of the simulation period. For 12 months, error bars for the nal year has been drawn as shown in Fig. 6; and, the results are indicative of proper accuracy of the veri cation stage. As it is clear, the difference between observation and simulated groundwater balance has a normal distribution with a mean value close to zero.

Predicting the future status of the aquifer After veri cation of aquifer's quantitative model and evaluation of the proper correlation between observation and simulated groundwater
balance, the future status of the aquifer has been predicted with consideration of current status to be continued for 6 years. Under such condition, climatic situation under three drought, normal, and wet periods have been considered in addition to the amount of groundwater exploitation based on current procedure. Accordingly, the highest stress and decline level of groundwater has been related to the drought period of the region. That is, the maximum decline has been at outlet part of the aquifer and equal to 40.6m. Under the same conditions, decline amount under wet period in the same observation well has been predicted to be 13.8m. Analyzing results related to aquifer local decline status under these three climatic conditions are presented in Fig. 7.
Also, in Fig. 8, the uncertainty band related to the prediction results of aquifer's hydrograph under three climatic conditions is shown. Hydrograph results of the aquifer show that drought scenario will result in the intensi cation of volume reduction of groundwater resources.

Simulating Land subsidence amount
Land subsidence resulted from groundwater exploitation would be formed due to the increase of interparticle pressure in alluvium. This pressure increase on free water tables results from solid particles losing their otation in a part of the table in which water level is low. In the harmonious decline of groundwater level and inhomogeneous density of alluvium will lead to inharmonious land subsidence and creation of some cracks in earth's crust. To model land subsidence in the region studied, after preparation of the quantitative model of aquifer, SUB package has been used so that data related to the structure of saturated and unsaturated environment would be entered. Information required for land subsidence simulation includes elastic storage, non-elastic storage, initial compaction, and pre-consolidation head. The non-elastic coe cient has been considered equal to 0.001 to 0.1, according to the range provided by the United States Geological Survey (USGS2010, 6-A23). One of the most important parts of modelling is de ning compressible layers within aquifer's limit that are usually obtained based on drilling logs. Under such condition, outlet limits of the aquifer and those limits with the maximum number of exploitation wells will have maximum land subsidence amount due to water discharge. Compressible layers are usually made of clay and absorbing water. They nd higher speci c gravity which will lead to aquifer decline and reduction of groundwater balance. In the model, two sections of compressible layers called No-delay interbeds and delayed interbeds have been de ned. No-delay interbeds were related to drainage rate, and it did not show any delay in water conveyance. Water ow rate under this condition will be calculated via conceptual model under the non-steady state of the model. Preconsolidation head as well is of high importance within the aquifer's limit. In the present study, interpolation and raw data as well as estimation of pre-consolidation head coe cients, have been used to divide the aquifer and values entered into land subsidence model. After the entrance of this information and using those data related to the quantitative model of the aquifer (hydraulic conductivity and initial groundwater balance), land subsidence has been simulated. After simulation of land subsidence model, subsidence amount for a period of six years upon the continuation of the current condition has been predicted for three droughts, normal, and wet periods. In Fig. 9, land subsidence amounts at the end of the prediction period under three climatic conditions are shown. As shown by the results, outlets of aquifer, especially in the north-eastern part of it have the highest land subsidence amount. Studying the aquifer model and behaviour shows that this part of aquifer has the highest numbers of exploitation wells from which the highest volume of water would be extracted. So, continuous unplanned usage would be indicative of aquifer's water decline and increase of land subsidence up to maximum 23cm, under normal climatic condition. Studying the relationship between the concentration of water wells and discharge of them in the region showed that outlets of the aquifer with several exploitation wells have been prone to maximum damages with consideration of aquifer decline. Also, studying land subsidence amount in other sections of aquifer shows that central and south-eastern parts of it are affected by the current procedure which in turn will lead to increase of damage.

Evaluating risk level of land subsidence
After land subsidence simulation and prediction of it for the future period and to study the role played by exploitation level on the amount of land subsidence; the correlational relationship between the two parameters according to the Fig. 10 has been extracted via Spearman method.
The results show that the average correlation between exploitation and amount of land subsidence is 0.62 and indicative of an effective relationship between exploitation level and amount of land subsidence based on statistical analysis. Maximum correlation is related to the south-eastern part of the aquifer and equal to 0.86; and, minimum correlation is 0.43 and related to north-western parts. Correlation results show that more than 70% of aquifer area is of low correlation which matches the results from land subsidence and predictions.
After specifying signi cant relationship between exploitation level and land subsidence, to formulate remedial mechanisms; exploitation level of the aquifer has to be decreased. To do so and to specify the exploitation level, exploitation risk indicator has been applied. This indicator has been extracted with consideration of correlation level between land subsidence and exploitation, as well as land use. Land use in the region has been extracted through Landsat 8 satellite images. In this indicator, the value of land subsidence probability has been estimated through correlation value and vulnerability value based on the amount of land subsidence and type of land use. Combining probability and vulnerability levels has been introduced as the exploitation risk map.
Risk: vulnerability X probability In Fig. 11, exploitation risk map of the aquifer is provided due to land subsidence in different climates.

Restoration scenarios for groundwater
Considering risk estimation related to land subsidence and prediction of land subsidence amount mechanism for reduction of water discharge has been simulated for different levels. Accordingly, the reduction of discharge between 5 to 30% has been simulated. To apply reduction level of discharge in the aquifer area and considering risk map; the aquifer has been classi ed, and the increase of discharge reduction has been applied to high-risk classes. Also, the lower reduction of discharge has been applied to lower risk classes. A linear relationship between subsidence risk and the discharge amount of groundwater resources has been used. In Fig. 13, local clustering of the aquifer from the perspective of land subsidence risk is shown in addition to the reduction level of discharge. Highest aquifer area is related to region 1 with an area about 74% of the region; however, the maximum number of exploitation wells and discharge amount is related to the region 2 with about 1700 wells. More than 46% of aquifer discharge from groundwater resources is related to this area. In region ve as the most sensitive region in terms of land subsidence, there are 50 wells with annual 18 million cubic meters of discharge.
Considering three climatic conditions, optimizing reduction level of discharge has been performed via LP algorithm with consideration of minimum land subsidence and improvement of the quantitative status of the aquifer. Reduction level of discharge has been applied in the conceptual model of the aquifer; and, the results have been provided in Table 3 as analysis of the reduction level of discharge in each region under three scenarios.  Maximum discharge reduction was related to cluster 5 with 30% reduction under all of the three climates. However, the maximum volume of discharge reduction was related to region 4 under drought period and equal to 22 million cubic meters. The results show that under the wet period, the best situation is observed in terms of reduction percentage of discharge and under this condition; the highest optimization level of quantitative status and land subsidence is observed. That is, in case of continuity of current situation, land subsidence has reached maximum 23cm in the outlet; and, under wet period and applying reduction of discharge, this has been reduced to 10cm. As shown by the result analysis, despite the effectiveness of the reduction level of discharge on quantitative status and maximum land subsidence in Nadjafabad Aquifer; average land subsidence of the aquifer had no signi cant change. This has been created due to the high expansion of aquifer and conditions created in previous periods.

Conclusion
Unplanned exploitation and no attention paid to exploitability potential of groundwater resources have led to the creation and expansion of land subsidence. The research studies land subsidence phenomenon and predicts the future status of the aquifer under this condition in three arid, normal and wet climatic conditions. Considering the effectiveness of immethodical discharge of aquifer on land subsidence, Nadjafabad Aquifer quantitative simulation has been performed; and, model veri cation and calibration has been performed based on initial simulation. showed that maximum land subsidence of Nadjafabad aquifer during the next six years would be 23cm at the outlet of the aquifer. Also, average land subsidence in the aquifer as a whole is 3cm. After prediction of land subsidence status due to groundwater discharge and to formulate remedial mechanism, aquifer risk map has been prepared based on which Nadjafabad Aquifer has been divided into ve areas, applying land use and correlation level between land subsidence and water discharge. This division has been done to apply the reduction condition of groundwater discharge. Then using LP algorithm, reduction of water discharge from the aquifer for each area under three climatic conditions has been optimized. Final results showed that under wet period, maximum land subsidence is 10cm which provides 14% improvement in the quantitative status of the aquifer. Under normal and drought condition as well, maximum land subsidence will reach 13 and 18cm, respectively. The results showed that average land subsidence in the aquifer as a whole has had no signi cant change with consideration of such phenomenon and expansion of the aquifer.
With consideration of the results obtained regarding current and future status of Nadjafabad aquifer; mechanisms provided in the region have to be implemented so that irrecoverable damages would not be imposed to the region's substructures. According to the research results, there is a direct relationship between land subsidence level of the aquifer and immethodical discharge and maximum subsidence affects those parts of the aquifer that have a high concentration of exploitation wells and high amounts of water discharge. Figure 1 The region studied. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.  Correlation between observation and simulated groundwater balance at the veri cation stage Predicting Nadjafabad aquifer decline status under three climatic conditions. a) Under drought period: highest amount of decline at the end of the six-year period at the outlet of the aquifer equal to 40.6m has been predicted. At north and rst part of the aquifer and with consideration of lack of concentration of exploitation wells, decline level is minor. b) Under normal climatic condition: maximum decline has been predicted at the end of the six-year period in the outlet of aquifer and equal to 18.3m. The minimum decline is related to the north part and equal to 32cm. c) Under wet period: Minimum decline compared to drought and normal periods will have no change; however, the maximum decline would be related to the outlet and it equals to 13.8m during six years. The results from land subsidence modelling at the end of the prediction period under normal climatic conditions Figure 10 Correlation between land subsidence and exploitation level. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

Figure 12
Aquifer clustering from perspective of land subsidence and reduction level of discharge. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.