Depletion zones of groundwater resources in the Southwest Desert of Iraq

: Aquifers o ﬀ er a reliable supply of high-quality water, making groundwater signi ﬁ cant in arid and semi-arid regions. Climate change is predicted to result in a decrease in rainfall and an increase in droughts. The prolonged drought severely devastated Iraq and is the main reason for the ongoing increase in groundwater consumption over the last decade. In this study, the stresses and depletion potentials of the Dammam con ﬁ ned aquifer, which extended along the Najaf and Muthanna governorates, are identi ﬁ ed and analyzed. By using the Groundwater Modelling System (GMS v10.4) software, a numerical simulation of groundwater ﬂ ow was used to study the Dammam aquifer system. The upper layer was modeled as uncon ﬁ ned, while the bottom layer was con ﬁ ned. The ﬁ ndings from the steady-state calibration indicate that the hydraulic conductivity (HK) within the study area varies between 1.47 to 20.0 m/day. Additionally, the recharging rate (RH) was estimated to be approximately 1.66 × 10 − 6 m/day. These parameters were utilized as the initial condition for conducting the transient analysis. Two operating scenarios were employed to perform unsteady simulations. The initial scenario involved the utilization of 89 production wells, while the second scenario included all 139 stand-by production wells, resulting in a total of 228 wells for the operational period from July 2021 to June 2022. The results of the ﬁ rst operation scenario showed that the drawdowns ranged from 0.4 to 5.8 m, whereas the second operation scenario showed that the drawdown increased from 1.0 to 22 m. The depletion rate in the groundwater static heads was measured by the percentage of relative di ﬀ erence. Hence, the depletion rate for the ﬁ rst scenario varied from 4.32 to 33.34%. On the other hand, the second scenario ranged from 7.45 to 33.34%.


Introduction
Groundwater is considered one of the most important freshwater resources, especially in arid and semi-arid regions with low annual rainfall and frequent droughts that are lately emerged as a serious global concern [1]. The prolonged neglect of regional water diplomacy has left Iraq's transboundary water resources vulnerable. Currently, Iraq suffers from reduced surface water in the Euphrates and Tigris rivers and their tributaries, which cannot meet the water demand. Generally, surface water will become more susceptible to contamination due to declining levels and rising concentrations of pollutants. Therefore, attention and emphasis must be directed to groundwater as an alternative water supply [2]. To address the increasing stresses on water resources, it is important to create a water-efficient society that focuses on water conservation and optimal use of the available alternative water resources. Numerous researchers studied groundwater modeling, evaluating management scenarios, and allowing decision makers to compare alternative measures and make management decisions to achieve efficiency goals without violating certain constraints [3]. Specifically, Al-Asadiy and Atiaa [4] studied the hydrological budget of the shallow Dibdibba aquifer in the Safwan-Zubair districts in Basra, Iraqi, for the period 1980-2000 to provide effective management of groundwater resources in the studied area. The amount of groundwater recharge was assessed according to the water surplus concept, corresponding to an average value of 0.252106 m 3 /day. AL-Fatlawi and Jawad [5] developed a mathematical model using MODFLOW software to evaluate groundwater availability and the changes in groundwater levels due to the pumping from the Umm Er Radhuma aquifer, Western Desert, Iraq. The calculation shows the aquifer receives 0.0578 mm/day as average recharge due to rainfall. Al-Mussawy [6] used MODFLOW 2000 package under the GMS 7.1 program to simulate the flow regime of the Dammam confined aquifer in the Karbala desert region, Iraq. The study presented that the vertical recharge rate (RH) varied from 2.74 × 10 −9 to 8.49 × 10 −8 mm/day. Al-Sudani [7] estimated the groundwater recharge in the Khan Al-Baghdadi district northwest of the Anbar governorate, Iraq, which equals 0.018 mm/day. Furthermore, Bayat et al. [8] investigated the effects of effective parameters on groundwater recharge, such as precipitation, surface recharge, and well water harvesting. GMS software was used to calibrate and validate the Karvan aquifer in Iran for 86 months to study climate change, water consumption control, and artificial recharge to improve water storage in the aquifer. In addition, the groundwater potential in Concepcion, Tarlac, Philippines, was explored by Quitaneg [9] using a groundwater flow model designed by GMS-MODFLOW to explain the aquifer's behavior in response to changing circumstances. The simulation reveals a reduction in hydraulic heads due to increased pumping induced by a 10-year groundwater demand projection equivalent to a 38.5% rise in the discharge rate. Mohammadi Arasteh and Shoaei [10] examined the hydrogeological effects of size and sediment granulation of the Zanjan Plain aquifer in western Iran on groundwater resources. Thus, between 2009 and 2012, hydrodynamic parameters were generated from a calibrated MODFLOW model, showing that hydraulic conductivity (HK) values are higher in foothill regions due to coarse-grained alluvial materials and that runoff inside Zanjan has a larger RH. This study aims to build a numerical model using the MODFLOW package integrated with the GMS program for a confined aquifer in the Dammam geological formation in the southwest desert of Iraq, which extends along the Najaf and Muthanna governorates, to estimate the hydraulic properties of the aquifer such as the HK and the RH at the steady state and develop a contour map for the piezometric heads. Also, the calibrated estimated parameters from the steady state are used as initial conditions for the transient state to investigate the effect of the pumping process from the wells that are distributed all over the study area by applying two realistic operating scenarios based on the water demand requirements in this region of Iraq to indicate the aquifer's severe depletion and stress zones by mapping the groundwater drawdown for the investigated area.

Research methodology
The study's tools and technology are discussed in the following sections. Numerical models can predict how the hydrodynamic system will react to human activities (groundwater exploitation) or natural processes (droughts) by visualizing the flow field in 3D, calculating areal water budgets, predicting exploitable groundwater reserves, etc. Thus, the Groundwater Modeling System (GMS 10.4) was used to calibrate a steadystate model and analyze HK and recharge in the Dammam confined aquifer. The transient condition runs using steadystate simulations to develop the groundwater stress map.

Study area
The study area is located in Najaf and Muthanna governorates in the southwest desert of Iraq. The geographic coordinates for the study area are between longitudes 44°19′ 07.68″ to 46°33′26.55″ east and between latitudes 29°06′ 15.74″ to 32°19′55.97″ north. The study area is of about 76,556 km 2 , and the base map and layout of the study area are shown in Figure 1.
The ground surface in the southern desert generally rises from the Euphrates River in the east toward the southwest in the direction of Saudi Arabia by about 50 m every 10-15 km [11]. The middle part is almost flat with many deep and shallow depressions; the depth of the Al-Salman depression in the Al-Muthana governorate reaches 40.0 m. The eastern part of the study area is mostly flat, partly covered by a thin layer of soil [12]. Topographic characteristics significantly impact the direction of groundwater movement and its recharged and discharged areas [13].
The topographic map that shows the major elevations of the study area was produced using ArcGIS 10.7 and a digital elevation model (DEM) of the SRTM type, with a 30 m accuracy, as shown in Figure 2.

Modeling of groundwater in the study area
The mathematical representation of governing equation for three-dimensional groundwater flow with a constant density nonequilibrium in a heterogeneous and anisotropic porous medium is expressed as follows [14]: where K x , K y , and K z are the HK values along the x, y, and z directions in (L/T); h is the water head in (L); W is a volumetric flux per unit volume and represents sources/ sinks of water in (1/T); S s is the specific storage for the porous media in (1/L), and t is the time scale in (T). For the case of (∂h/∂t = 0), which represents the steady state of groundwater flow, equation (1) is reduced to the following form: Equations (1) and (2) are solved numerically using the finite difference method. A three-dimensional, numerical time-dependent flow model of finite differences was performed using a GMS to simulate and solve the transient flow condition. In the study of the impacts of the actual withdrawal rates of the Dammam confined aquifer, however, a conceptual model is developed by creating two layers of flow domain, setting up the boundary conditions, initial values for the aquifer parameters, and conducting the estimation for parameters and the final calibration of the model.
The geospatial analysis, visualization, and DEM processing are performed by ArcGIS v10.7 software. The methodology used in the numerical modeling of the groundwater flow for the selected study area is presented in Figure 3.

Model conceptualization
Model conceptualization characterizes a description of field conditions by systematically identifying coverage  used to define sources/sinks and boundary conditions. Areal properties were used to define the top and bottom elevation of the aquifer layers and the aquifer properties such as RH and HK for steady-state simulation and storage coefficient for unsteady-state simulation. Observation points are used to define the observation well point, represented as (Head) in the case of static/steady state and the term of (Trans. Head) for the dynamic/transient state. The required coverages were created to describe the groundwater flow processes for the modeled area. After completing all needed coverages, the 3D grid is generated by converting the conceptual model into the 3D grid model.

Grid design
A finite-difference cell-centered 3D grid of 200 rows and 200 columns was developed using the MODFLOW package integrated with the GMS program. The locations of cells are labeled in terms of rows (I), columns (J), and numbers of layers (K). The modeled domain consisted of 80,000 cells, 41,530 were active cells, and 38,470 were inactive; each cell is compared to each of the coverage's polygons; if a cell does not fall inside any of the polygons, it is considered to be outside the model's domain and is inactivated while cells inside the domain adopted as active cells [15]. Each cell had a dimension of 1953.5 m × 1879 m with an area of 3.64 km 2 , for each cell, as shown in Figure 4.
The model consisted of two layers, confined and unconfined (K = 2). The efficient groundwater models depend on the accuracy of the aquifer parameters [16]. Consequently, the spatial variation of aquifer thickness was considered, and the spatial distribution maps were generated using 3D analysis tools powered by the ArcMap program, as presented in Figures 5 and 6.
The top of the first layer equals the static water level, as shown in Figure 7. The Kriging interpolation method powered by the Arc-GIS toolbox was used to form the groundwater levels contour map for the Dammam aquifer. Generally, in the groundwater contour map for the studied area, the zonal groundwater elevation ended with 20.78-69.08 m in the discharge area (NE), while it was 208.423-257.65 m in the recharge area (SW), which indicate that the groundwater flows from the southwest to the northeast. The second layer is the confined Dammam aquifer, which had a thickness ranging between 85 and 165 m. The top of this layer was  equal to the bottom of the unconfined aquifer, while the bottom layer was 125 m below sea level; this two-layers system was built by assigning the grid cells' HK, transmissivity, and storativity.

Boundary and initial conditions
The boundary conditions can be either physical or hydraulic. For example, impermeable rocks, lakes, or waterways are common types of physical boundaries. On the other hand, water divides or flow lines are common examples of hydraulic boundaries [17]. Boundary conditions refer to hydraulic conditions along the perimeter of the problem domain and can be mathematically classified as follows [14]: The specified Head Boundary is also known as the Dirichlet condition. The head along this boundary is set to a known value. Depending on the space, the heads of this type can vary. The constant head boundary is a special case of this boundary, where the heads along the boundary  have a constant value. The specified flow or a Neuman boundary condition gives a specified derivative of the head along with the boundary, as shown in Figure 8.
The flux boundary condition imposed at the Δx Δz face can be written as follows: where (∂h/∂y) is the y-component of the head gradient and q y is the y-component of the specific discharge vector. The no-flow boundary is considered a special case of this boundary, where the flow across the boundary is set to zero. The head-dependent or Cauchy condition, in this case, the flow throughout the boundary, is determined by the gradient between the head outside the boundary and the head estimated by the model at the node that is on the boundary or very close to it, which can be expressed as follows: where h i j k , , is the computed head at the ell center and h b is the specified head along the face pf (ΔxΔz) located at a distance of (Δy/2) from the cell center.
Based on the regional groundwater flow pattern of the Dammam aquifer, the constant head boundary (CHD) was set away from the groundwater wells field to limit influences on simulated heads within the model domain [18]. The Euphrates river was modeled using an arc with defined endpoints (nodes); each node required a river head stage and bottom elevation provided by the General Committee for Groundwater/Ministry of Water Resources, Iraq, as shown in Figure 9 and presented in Table 1.
For the initial condition, the values of input parameters such as RH and HK are required for the parameter estimation tool (PEST) process in the steady state.
The HK of the groundwater in layered aquifer systems is based on the equivalent magnitude of the whole system's hydraulic conductivities [15]. GMS provides more than one way to assign HK to the groundwater model. In this study, the HK was assigned to polygons on a single coverage due to the use of the conceptual model approach and then mapped to the MODFLOW model. Consequently, the study area was divided into five zones to input the initial values of HK, as shown in Figure 10. After inputting the initial HK values, the estimation process is carried out using an inverse model tool (PEST) under GMS software in which the user-defined set of initial input parameters (HK and RH) is systematically adjusted until the difference between the computed and observed values is minimized. These initial values were adopted from the General Committee for Groundwater/Ministry of Water Resources, Iraq, as presented in Table 2.
In contrast, the initial recharge value was calculated from the metrological datasets from 1976 to 2020 [19]. The   annual rainfall was 103.3 mm for the Muthanna station and 98.8 mm for the Najaf station. The suggested initial value for the RH was 5% of the daily average rainfall [15]. As a result, the initial RHs for the studied governorates were 0.000014 and 0.0000135 m/day. Since the difference between the calculated recharge values for the studied governorates is insignificant, the study area was considered a single zone with an initial recharge of 0.000014 m/day.
The initial and final estimates for HK and RH are shown in Table 2.

Evaluation of steady state
The field observed values would be compared automatically with the calculated values obtained by the model during the calibration process. The residual values were computed for a (95%) confidence interval and an allowable head interval of 1.5 m. The calibration targets illustrated the validity of the calibration process. If an observed value has been assigned to a feature object, the calibration error at each object can be plotted using a calibration target presented in Figure 11. A set of calibration targets provides useful feedback on the calibration error's magnitude. The observed value is represented by the center of the target. The upper limit of the target is equivalent to the observed value added to the interval, while the lower limit corresponds to the observed value reduced by the interval. The colored bar denotes the presence of an error. When the bar is completely contained within the target, it is represented by the color green in the visualization. If the error in the observation target exceeds the permissible limits by a small margin, it will be indicated by a yellow highlight on the error bar. On the other hand, when the error exceeds 200%, it is denoted by the color red. The calculated and observed groundwater level measurements exhibited a high degree of agreement in the modeled region. The contour map showing the simulated piezometric heads of the aquifer has been presented in Figure 12. Out of the 25 observation targets in the steady state simulation, 21 of them have been denoted with green    error bars. Additionally, one of the observation target wells has been found to be in perfect agreement with the observed values. Merely three error bars were colored in yellow; however, they were found to be within the permissible limits of the target, thereby signifying a substantial enhancement over the initial solution.
The study conducted an error analysis, as reported by Anderson and Woessner (1992) [14], in which the Mean Error (ME), Mean Absolute Error (MAE), and Root Mean Square Error (RMSE) were employed as metrics to assess the model's reliability.
Mean error (ME) is the difference between the residual errors and can be expressed as follows: where n is the target number, i = 1,…,n, h m is the measured head in (m), and h s is the simulated head in (m). Mean absolute error (MAE) is the mean of the residual's absolute value, which was often greater than and stronger predictor of model fit compared with ME; this measure can be expressed as follows: Root mean square error (RMSE) is the average of the squared residuals. Typically, the RMSE is larger than the MAE. It can be expressed as follows: The ME, MAE, and RMSE were: 0.084 m, 1.166 m, and 1.314 m, respectively. The model is calibrated if the RMSE is less than a certain percentage of the calibration target range of values; There is no allowable magnitude of the ME, MAE, or RMSE; other than that, these values should be minimized [14].
The calculated values of these statistical measures indicate that the conceptual model, boundary conditions, and final hydrological parameters used in the model were reasonable and could be further used since the head calibration target is set to (1.5 m), and the calculated RMSE was found to be 1.314 m. Statistical error measures for the steady-state calibration are shown in Figure 13.

Calibration process for the transient state
The groundwater levels have fluctuated from one-time step to another in the transient state of the flow of groundwater. In such a case, the hydraulic heads obtained from the steady state were adopted as the initial heads to conduct the transient analysis; the analysis required specifying storage parameters for the aquifer. A realistic operating scenario was established based on the water demand requirements in this region of Iraq, and a 1-year operation period was studied using a 30-day time step to examine the stresses on the storage of the Dammam aquifer system throughout the study area. The calculated distribution of hydraulic head by the analysis shows different measured findings of drawdown due to the pumping process, which are compared with the observed hydraulic head for three wells in different locations due to the availability of transient head data for the selected period, as presented in Figure 14. Therefore, the transient state's calibration process was started from July 1, 2021, until June 1, 2022. The data used in the calibration and validation included the observed transient head for the selected stress period. All the hydrogeological and hydraulic parameters were extracted from the calibrated steady state. The calibration process for the transient state has been performed for a specific storage parameter that changed during the transient period by using a trial-and-error procedure.

Results and discussion
The results of groundwater flow in both steady and unsteady states will be discussed in the following sections:

Groundwater flow in steady state
The results of the calibrated steady-state model were revealed by a contour map of piezometric head distribution for the confined aquifer, as shown in Figure 15, which indicated that for 25 calibration targets performed by the model, practically all errors were within the specified interval, although several piezometric heads matched actual values. The statistical measures of error, such as MAE and RMSE, confirm this inference, that is, MAE = 1.166 m and RMSE = 1.314 m. Figure 16 shows the relation between observed and computed head values for the wells in the study area. It was found that the coefficient of determination (R 2 ) is equal to 0.9997, which reflects a perfect correlation since the simulated head values are close enough to the targeted head values constrained by an allowable error interval equal to (1.5 m).   The residual head values, representing the variance between the observed and estimated head, were discovered to be below the allowed error limit (1.5 m), as presented in Table 3. Furthermore, for the selected wells for the validation process, their calibration targets were displayed in green color, which indicated a good validity of the model of R 2 equaled 0.9997, i.e., excellent fit.

Scenarios of transient state of groundwater flow
After a good match between the observed and calculated groundwater heads at the transient state was achieved, the specific storage varied from 0.00005 to 0.00065. Then the two realistic operational scenarios were introduced to help researchers better understand where groundwater stress might occur in the study area. The first scenario was adopted based on the fact that the study area is facing a decline in surface water inflow and a significant increase in water demand. As a result, 89 wells were pumped at their full capacity simultaneously. Wells have an actual pumping rate ranging from 259.2 to 1,728 m 3 /day, depending on the conditions. The wells averaged a pumping rate of 565.3 m 3 /day adopted for each well of the 89 wells. For a whole year, the model was subjected to these conditions. It was calculated that the total daily discharge from all 89 wells is (50310.72 m 3 /day).   Figure 17 shows the results of the transient state of operation of the groundwater flow model, which indicated that the drawdown in the groundwater level ranged from 0.4 m to its maximum value of 5.8 m during the operation period from July 1, 2021, until June 1, 2022. The second operation scenario includes using 228 wells in the study area. Due to the increased water demand, 139 were added beyond the existing wells to cover the water demand requirements. The wells operated by their actual pumping rate of an average discharge rate equal to 581.3 m 3 /day. The discharging rates were found to be 132546.24 m 3 /day from 228 pumping wells.
Drawdown calculations indicate that the drawdowns range from 1 to 22 m. The simulation results of such a scenario are shown in Figure 18.
The high drawdown values were recorded in the populated areas as a direct result of the increase in the numbers and the rate of discharge of water wells in such areas due to the increase in the water demand, particularly the expansion in agricultural activities and the high growth of population in the surrounding regions and its extended to the study area. In contrast, low drawdown values are due to the scattered well distribution. Consequently, many wells can be drilled in those areas of low drawdowns.
The depletion rate in the groundwater static heads was measured by the percentage of relative difference (RD%), a measure used to compare the predicted values by model by a reference value expressed as a percentage. The reference values were represented by the minimum and maximum values of the steady-state head, which ranged from 30 to 255 m, while the head distribution ranged from 20 to 244 m and from 20 to 236 m for the first and second operational scenarios, respectively; therefore, RD% is written as follows:  where h p and h r are the predicted and reference values of static heads in (m), respectively, in the specified location within the groundwater flow domain. Consequently, for the first operation scenario, the values of RD% ranged from 4.32 to 33.34%. On the other hand, in the second operational scenario, RD% values ranged from 7.45 to 33.34%. Groundwater potentiality is significantly affected by land use and land cover (LULC), which in turn affects the control of surface runoff and soil infiltration [20]. Land use/land cover affects the presence and the change of groundwater in the area, either by lowering runoff and encouraging infiltration or increasing surface runoff and impeding infiltration. The LULC increase is expected to increase the evaporation rate and depletion in the groundwater level significantly and quality; consequently, groundwater depletion is the source of major economic and environmental consequences [21,22].
As a result, the groundwater potential is quite low in urban areas due to the prevalence of impervious surfaces such as roofs and roads that reduce the amount of water that naturally infiltrates into the ground and eventually reaches the water table [23]. The LULC map was prepared by processing Sentinel-2 images for 2021, as shown in Figure 19. LULC is an important factor in the existence of groundwater and how it changes [24]. Consequently, the region under examination was divided into six classes: water, herbs wetland, herbaceous vegetation, cropland, and bare/sparse vegetation; each class represents the priority effect on groundwater accumulation. Changes in particular land cover classes, such as the amount of water bodies and forest vegetation present, may impact groundwater RH and other hydrological components. These changes affect the interception and infiltration processes [25].

Conclusions
The evaluation of groundwater resources in the Dammam confined aquifer system in the southwest desert of Iraq has been carried out depending on numerical simulation by three-dimensional groundwater conceptual model under GMS software for both steady and transient states. The calibration processes and error analysis revealed that the model was effectively calibrated due to the calibration targets and statistical measures. Two realistic operating scenarios were adopted.
Based on the simulation findings of the GMS and GIS program for the study area, the following may be concluded: • According to the groundwater contour map that was developed for the examined region by the ArcMap program, the groundwater flows from southwest to northeast with a groundwater elevation of 208.423-257.65 m in the recharging area (SW)), while it is 20.78-69.08 m in the discharge area (NE). • The calibrated steady-state simulation using GMS revealed that the study region's HK varied from 1.47 to 20.0 m/day. In contrast, the RH was equal to × − 1.66 10 6 m/day. • Based on findings of the simulation of static conditions, the contour map of piezometric heads showed that the higher heads occur to the southwest of the study area at an elevation of about 260 m, which corresponds to a recharge zone. In contrast, lower heads are located in the discharge zone close to the Euphrates river at an elevation of 10 m, showing a good agreement with the groundwater flow map prepared by the ArcMap software using field data. • The transient-state calibration showed a satisfactory agreement between the measured and predicted groundwater heads from 1 Jul 2021 to 1 Jun 2022. This indicates that the specific storage coefficient in the study region ranged from (5 × 10 −5 to 6.5 × 10 −4 ). • The results of the first scenario indicated that the depletion rates in groundwater heads ranged between 4.32 and 33.34%. In comparison, the second scenario ranged from 7.45 to 33.34% • A drawdown of 0.4-5.8 m in the groundwater level is indicated by the operation of 89 wells at their full capacity. While in the case of operating 228 wells, the drawdown is 1-22 m due to an increase in total pumping rates from 50310.72 to 132546.24 m 3 /day. Consequently, low drawdown areas identify a good potentiality zone for drilling additional water wells since the simulated model reveals that the pumping process has moderately influenced such areas.