2.1. 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 identifiable 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 confirmed 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.
2.2. 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 flow 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.
2.3. 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 first step for simulation is preparation of conceptual model of aquifer based on which steady-state and transient-state models of flow would be simulated. After verification and calibration of aquifer's quantitative model, land subsidence would be predicted finally 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.
2.3.1. Preparing conceptual model of aquifer
The first 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, finite 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, flowback water from discharge of wells, runoff and rainfall infiltration, 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, coefficient of storage, and specific 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 figure 2, the overall view of aquifer's conceptual model is presented.
2.3.2. 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 infiltration of runoff and rainfall, as well as inlet boundaries to it. Under this condition, infiltration from the surface and that of groundwater flows 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 flows from the aquifer has been considered similar to inlet flow 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 deficient in terms of reservoirs; and, there have been 138 million cubic meters reservoir deficit as calculated in table 1.
Table 1. Nadjafabad Aquifer balance
Volume variation
|
Discharge
|
Recharge
|
Sum
|
Groundwater outflow
|
Groundwater evaporation
|
Discharge from well, spring, and qanat
|
Sum
|
Flow back from drinking and industrial water
|
Flow back from agriculture
|
Runoff
|
Rainfall infiltration
|
Groundwater inflow
|
-138
|
886.6
|
0
|
0.9
|
885.7
|
748.6
|
58.9
|
424.8
|
189.5
|
26
|
49.4
|
2.3.3. 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 flow 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 of the system, basic groundwater equation would be in the form of Poisson’s equation. In stable regime and in adaptation stage, correction of conveyance capability of water table or hydraulic conductivity of it as well as the natural amount of recharge and discharge are concerned. Accordingly, October 2011 has been considered as the starting point for steady-state simulation. After adaptation and calibration of groundwater steady-state model, transient state simulation of groundwater flow has been taken place. To do so, a period of five years, including 60 monthly time steps till September 2016 has been considered as transient state simulation period.
2.3.4. Calibration and verification of quantitative model
Calibration and analyzing sensitivity are two stages of operations in groundwater modelling, taking place consecutively. Accordingly, and with consideration of flow equation under steady-state and studies performed in relation to groundwater flow simulation, hydraulic conductivity parameter at steady-state have to be studied as the most important parameter for calibration. After calibration of the steady-state model of flow, transient state model of flow has been scrutinized and evaluated based on water yield parameter in free aquifer. One of the calibration methods in mathematical models of groundwater is using automatic calibration software packages. In fact, these software packages are responsible for solving progressive groundwater equations via estimation of certain parameters; and, final values of the parameter resulted from these software packages have best answer regarding the feature measured in the aquifer. To do so, the manual method and PEST model have been used in calibration. After calibration of the steady-state model and evaluating the accuracy of the results obtained, verification of the model has been taken place for the period of one year till the end of the water year of 2016-17. Verification is one of the necessities of a proper simulation and model’s behavior against discharges and information related to previous times would be measured based on observation results and model response. Duration of adaptation for verification depends on the two factors of type of aquifer and historical course of it.
2.4. 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 coefficient, 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 fine-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 fine-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 confined 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 ∆b [L]is thickness change of layers, ∆h [L]is the change of hydraulic head, Sske [L-1] is texture component of especial elastic storage, and b0 [L]is initial interlayer thickness. For free aquifer, elastic compaction component or expansion of sediments will be calculated as follows:
Where n is porosity, and nw is moisture content of soil over groundwater balance (unsaturated area) which is a fraction of total volume of porous environment.
Non-elastic storage factor: For confined aquifer, the following equation would be used to calculate non-elastic compaction of ∆ b*[L]:
Where Sskv[L-1] is the texture component of non-elastic especial storage. For free aquifer, non-elastic compaction of sediments can be stated as below:
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.