Depletion zones of groundwater resources in the Southwest Desert of Iraq

2.1 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 k m 2 , and the base map and layout of the study area are shown in Figure 1.

Figure 1

Layout map of the study area.

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.

Figure 2

Study area topography.

2.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.

Figure 3

Groundwater modeling flowchart for the study area.

2.2.2 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.

Figure 4

3D grid of the modeled study area (Z magnification = 50).

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.

Figure 5

Thickness of the first layer of the Dammam aquifer.

Figure 6

Thickness of the second layer of the Dammam aquifer.

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.

Figure 7

Map for groundwater flow in the Dammam formation.

2.2.3 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.

Figure 8

Control volume of groundwater flow in the porous medium [14].

The flux boundary condition imposed at the ∆x ∆z face can be written as follows:

(3) ∂ h ∂ y = − q y K y ,

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:

(4) q y = − K y h b − h i , j , k ( ∆ y / 2 ) ,

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.

Figure 9

Boundary conditions of the study area.

Table 1

Euphrates river head stage and bottom elevation at each endpoints node

Governorate name	Head stage (m)		Bottom elevation (m)
	Start point	Endpoint	Start point	Endpoint
Najaf	21	21	15	14.5
Muthanna	14	13.5	8	7

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.

Figure 10

Initial HK in (m/day) for the study area zones.

Table 2

Parameter estimation outputs

Parameter key	Initial value (m/day)	Estimate value (m/day)
HK_20	2.30	1.46569
HK_40	5.06	2.67806
HK_60	7.97	4.12694
HK_80	9.56	20.0
HK_100	7.97	13.3331
RH_120	1.4 × 10−5	1.66 × 10−6

^*HK is the HK of the aquifer.

^*RH is the recharge rate for the aquifer.

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.

2.3 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.

Figure 11

Calibration target.

Figure 12

Model calibration target for the present study.

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.

Figure 13

Statistical measures for the model calibration process.

2.4 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.

Figure 14

Location of the transient-state observation wells.

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.

3.1 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 15

Contour map of piezometric heads in meters.

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).

Figure 16

Observed vs computed heads correlation for the study area.

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.

3.2 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.

Figure 17

Groundwater drawdown contour map for the first scenario.

Drawdown calculations indicate that the drawdowns range from 1 to 22 m. The simulation results of such a scenario are shown in Figure 18.

Figure 18

Groundwater drawdown contour map for the second scenario.

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].

Figure 19

The LULC map.