Groundwater Flow Modeling for Qushtapa Plain Unconfined Aquifer in Southern Erbil Basin, Kurdistan Region, Iraq ()
1. Introduction
The continuous exploitation of groundwater from aquifers for all purposes leads to the depletion of groundwater in many parts of the world (Ramesh & Fritz, 2016).
Groundwater models contribute additional insight into the complex system practice and can assist in advance special management strategies.
Groundwater modeling is a crucial tool to estimate and evaluate the groundwater flow and quantifying it’s potential. A simulation of groundwater flow using mathematical model indirectly by means of dominant equations represents the physical process in the system, with equations that describe the head or flow along the boundaries of the model (Anderson & Woessner, 1992). The same simulation applies to groundwater flow and contaminant transport modeling (Shwani, 2008) indicate that there is groundwater pollution in northwestern part of the study area as consequence of the human contamination, Agricultural crops and using insecticide in the soils to prevent the bacterial effect in the area of investigation.
MODFLOW-2000 is the most widely used ground-water flow simulation program which is the only ground-water flow model capable of calculating accurate sensitivities for simple and complex flow systems (McDonald & Harbaugh, 1988; Harbaugh & McDonald, 1996; Harbaugh et al., 2000).
Designing MODFLOW finite difference model direct approach is less intuitive, precisely for complex boundary conditions. Therefore, a model can be developed either using a grid or conceptual model approach (Sohrabi et al., 2013).
The main objective of this research is to investigate the groundwater flow system in the Qushtapa District area to develop conceptual models and accordingly to construct numeric groundwater modeling which can be used to simulate the groundwater flow under steady state flow conditions, also to generate hydraulic head in the area, calculate velocity distributions, and determine travel time and travel paths of contaminants in the area to selecting the wellhead protection area.
The model can be used to predict future groundwater flow conditions, estimate an aquifer’s hydraulic response, and predict the pumping rate needed to monitor the well discharge; thus, it can be used to predict the area’s water availability and sustainability.
In Kurdistan Region groundwater represents a main source of freshwater and water resource, and therefore it is important to study the groundwater systems in order to maintain this vibrant source and provide necessary information about this source in the area.
2. Modeled Area Description
Study area is located on the southern Erbil city about 33 km between latitude (36˚02'15"E to 35˚45'43"E) and longitude (44˚07'45"N to 43˚52'19"N). The area of interest is 426 km2, the elevation ranges between 150 - 700 m above sea level (Figure 1). The area bounded by lesser Zab River in the southern part, and the Avana anticline in the western part. The modeled aquifer in the area is unconfined aquifer.
![]()
Figure 1. Location map and digital terrain model of the studied area with observation wells locations.
3. Tectonic Setting and Stratigraphy
Tectonically; the area lies within the Unstable Shelf Zone that was affected by the Alpine orogeny in Mesozoic in Chamchamal-Butma sub-zone of the Foothill Zone which is structurally highest part of the zone (Buday & Jassim, 1987; Jassim & Goff, 2006).
The Foothill Zone region is characterized by widespread development of the fluviatile sediments filling broad and relatively shallow synclines between the narrow anticlinal ranges.
Stratigraphically, Bai Hassan formation and Quaternary deposits are the main outcrops in the area (Figure 2). The formation consists of molasse sediments represented by alternation of claystones, conglomerates, sandstones and siltstones. The claystones are reddish and pale brown, green, yellowish green colors, fairly hard, occasionally silty and sandy with small lenses of conglomerates.
The conglomerates form the main constituents the upper parts of the formation. The gravels of the conglomerates are mainly of limestone, chert and greenstone (Buday & Jassim, 1987). The depositional environments of Bai Hassan Formation is fluviatile which interchanges with alluvial fans or lacustrine environments (Mahmood, 1986).
![]()
Figure 2. Geological map of the study area (after Stevanovic, 2003).
Quaternary Deposits cover rock units in the study area, the deposition and stratigraphic sequence of the Quaternary sediments depend on the climatic oscillations, causing periodically repeated phases of accumulation and erosion. Quaternary Deposits are divided into four types according to (Youkhanna & Sissakian, 1986); they are Slope Deposit, River Terraces, Polygenic Deposit and Flood Plain Deposit. Most part of the study area covered by Slope Deposits and River Terraces.
Geomorphologically, the studied area is considered as a plain area, the main drainage patterns are dendritic. The elevation ranges between 650 to 250 m above sea level, the highest area in the north eastern part while the lowest area in the southern part near the Lesser Zab River. The dominant erosion types are eolian and water erosion (Al-Mahdi, 1979).
4. Materials and Method
Parameters are defined in the groundwater flow process input files and can be used to calculate most model inputs including model layers, horizontal and vertical hydraulic conductivity, horizontal and vertical anisotropy, specific storage, specific yield, head boundaries, areal recharge rate, evapotranspiration, hydraulic and constant head boundaries (Figure 3).
Twenty one Observation wells and 4 pumping wells was used for creating the model with coordination (UTM “WGS-84”) and water tables information to determine the head equipotential and drawdowns in the study area and measuring gains and losses along head-dependent boundaries like river and measuring flows through constant head boundaries (Table 1).
![]()
Figure 3. Flow chart shows the steps of the modeling using visual MODFLOW.
![]()
Table 1. Observation wells used in the model, UTM (WGS-84) coordination system.
4.1. Data Collection
The hydrogeological and geological information and data collected for the groundwater model include information on surface and subsurface geology water table, precipitation, evapotranspiration, pumped abstraction, stream flows, boundary condition, hydraulic properties, geological formations, lithological descriptions and a topographic map with observation wells and boreholes location.
4.2. Model Design and Geometry
Study area defined into the model as a one layer unconfined aquifer. The study was horizontally discretized 100 (x) × 100 (y) (10,000 horizontal grid cells). Maximum and minimum thickness of the model was 650 m, 50 m, respectively. The geometry of the area was defined by using a map with UTM coordinates (Figure 4).
![]()
Figure 4. The modeled area design with locations of the observation wells (cross points) and Lesser Zab River (blue line).
The surface layer was imported from the digital elevation data of SRTM 30. Although the entire area is built up by sediments it is assumed that a continuum approach can be applied. Thus, the Darcy equation can be used to describe and model the groundwater flow. A finite difference approach was chosen with equal grid spacing. The elevation of the bottom of the aquifer is unknown due to lack of geophysical data; therefore the constant value for all cells was used to simplest assumption for the bottom of aquifer.
4.3. Conceptual Model
Conceptual model construction includes the definition of the basin boundaries, aquifers recharge, and discharge sources (Anderson & Woessner, 1992). It is the first step for building groundwater model which is simplified representation of reality with a focus on the geological and hydrogeological conditions. Groundwater system in Qushtapa Plain conceptualized as a single unconfined (Inter-granular) aquifer consisting of two areas having different hydrogeological properties; Quaternary deposits represented by alluvial deposits and river terraces; and Tertiary formations represented by Bai Hassan Formation.
4.4. Boundary Conditions
Groundwater flows from north eastern to the south part of the study area. Different types of boundary conditions (BC) were used for the modeled area representing by; no flow boundary conditions, River boundary condition represented by Lesser Zab River (the river setting is shown in Table 2), Constant head boundary conditions, Recharge and Evapotranspiration boundary conditions (Figure 5).
Constant head values were selected according to the observation wells close to the boundary by using the line-option in the model. Constant head values were selected according to the observation wells close to the boundary by using the line-option in the model. Recharge of the area set to 38.17 mm/year according to (Shwani, 2008).
![]()
Table 2. Rivers boundary conditions values for model area.
![]()
Figure 5. Boundary condition map (Red line is the constant head B.C, Blue line is River B.C and Green area is no flow B.C).
4.5. Hydrodynamic Characterization
Using the Aquifer Test program, hydraulic properties of the model areas are collected from the well test data for four pumping wells. It is possible to estimate the physical properties of the aquifer from single well test (Kruseman & de Ridder, 1994; Dellur, 1999; Schaaf, 2004). The parameters include transmissivity, hydraulic conductivity, specific yield and storage coefficient. Jacob and Theis methods were used to determine hydraulic conductivity and transmissivity.
Hydraulic conductivity and specific storage of the unconfined aquifer system of Qushtapa plain were obtained from the analysis of pumping tests that were carried out in the boreholes in the plain (Table 3).
Theis’s recovery method (Theis, 1935) and Jacob method (Jacob, 1948) has been used for the calculation of Transmissivity and Hydraulic conductivity. These methods was used because they are the commonly applied methods for pumping test analysis and recommended methods for analysis of single-well analysis also takes into account the well-bore storage. (Al-Sawaf, 1977; Todd, 1980) also used these two methods.
Specific yield calculated in this research according to Johnson (Johnson, 1955). Specific capacity calculated for this study according to Fetter (Fetter, 1994) by the formula below (Table 4).
In Jacob method the relationship between correction drawdown and time drawing on semi-logarithmic paper, while in Theis’s method the relationship between residual drawdown and time drawing on semi-logarithmic paper (Appendix).
![]()
Table 3. Hydraulic parameters according to the is and Jacob method by using aquifer test program, and specific yield according to Johnson for study area.
![]()
Table 4. Specific capacity of the wells in the studied area according to by Fetter.
The effective porosity was assumed to be 0.15% and total porosity was assumed to be 0.3%. The storage coefficient is simply expressed by the product of the volume of aquifer lying between the water table at the beginning and at the end of period of time and the average specific yield of the formation (Todd, 2005). The storage coefficient was set to 1e-5 and specific yield to 9.65%. The Initial head set to 231.7 m.
5. Result and Discussion
5.1. Model Result
The simulation of the groundwater head in the area shows that the groundwater head starts from the northeastern part of the plain and decreases towards the south of the plain from 420 m to 140 m above sea level (Figure 6). The equipotential lines show that the Lesser Zab River drains the shallow groundwater. A minor relationship between groundwater and the corresponding river can be recognizing due to either drainage system or because of water well pumping. The flow direction is from northeastern part to the southern part toward Lesser Zab River.
5.2. Model Calibration
Model calibration was performed for the twenty two observation wells in the area, the average groundwater head observations were compared to the simulated groundwater heads. The calibration procedure involves manual trial and error adjustments of the most sensitive parameter (hydraulic conductivity, recharge, and conductance) within reasonable ranges until reaching the minimum residuals between computed and observed heads. The standard error of the estimate is 4.881 m, the root mean square is about 27.381 m, Normalized RMS is 8.373%, the maximum and minimum residuals are 51.41 m in well number 21, and 1.166 m in well number 2 respectively, and the residual mean and the absolute residual mean are15.793 m and 24.293 m respectively. Correlation coefficient of the comparison results is 0.967 (Figure 7).
5.3. Particle-Tracking Simulations
MODPATH is a post-processing particle tracking program designed to measure three-dimensional flow paths using steady-state or transient simulations of groundwater flow. It uses a semi-analytic particle tracking scheme that allows each finite-difference grid cell to obtain an analytical expression of the particle flow. Particle tracking simulations were performed on steady state conditions in both forward and backward mode.
The particle paths are determined by monitoring particles from one cell to the next until the particle hits a boundary, an internal source, or some other termination requirements are met. MODPATH offers the possibility of monitoring particles forward in the direction of groundwater flow and backward in the direction
![]()
Figure 6. Groundwater head equipotential contour line in the study area (a), Groundwater Flow direction (b), contour interval is 50 m.
![]()
Figure 7. Comparison of measured and simulated groundwater heads (m a.s.l) under steady state flow condition.
of recharge points. In the forward tracked; the recharged water at the mountain front to estimate the travel time from recharge to discharge area, which the particles terminate at points of recharge, rather than points of discharge. The particles were located in the potential recharge areas in the north and north eastern part; the particles were tracked forward along the flow paths were 20 particles assigned in the recharge area (Figure 8).
The backward simulation used to estimate contaminant’s travel time to the extraction wells were 10 particles assigned in the center of each abstraction well for best path line visualization and the particles were tracked backward.
The results of the forward tracking show the source of potential pollutants from the recharge area after 5, 10, 25, 50, and 100 years travel time, the particles released at the northern boundary travels to the center and the western part toward the pollution source near the Qurshaghlu village (Figure 9). The results of the backward tracking shows the particles tracked toward the recharge area in the north eastern part of the study area. The particle distance from the extraction wells after 5, 10, 25, 50 and 100 years are shown in Figure 9.
The coordination of the extraction wells is shown in Table 5. The particles located in the extraction wells moved toward the north and north eastern part to the recharge area. The recession in the mountain front particle track can result from the low permeability in this region of the compacted rocks.
![]()
Figure 8. Particle tracking simulation for potential recharge areas cells for forward tracking (travel time of 5, 10, 25, 50, and 100 years).
![]()
Table 5. Extraction wells coordination (UTM-WGS84).
![]()
Figure 9. Particle tracking simulation for potential recharge areas cells for forward tracking (travel time of 5, 10, 25, 50, and 100 years).
6. Conclusion
The steady state flow model of the Qushtapa plain was applied to determine the flow direction and estimate the contaminant sources at different travel times. The aquifer was simulated under unconfined condition and is represented by a single layer of 100 m constant thickness. The model shows good agreement between observed and calculated water levels. The simulation of the groundwater head equipotential in the area shows that the head starts from the northeastern part of the plain and decreases towards the south of the plain from 420 m to 140 m above sea level.
The particles in the forward tracking released at the northern boundary travel to the center and the western part toward the pollution source near the Qurshaghlu village, while the results of the backward tracking show that the particles located in the extraction wells moved toward the recharge area in the north and northeastern part of the study area.
Acknowledgements
This work was supported by the Kurdistan Ministry of Higher Education and Scientific Research, Salahaddin University—Erbil, Agricultural Engineering Science College, Soil and Water Department.
Appendix
![]()
Figure A1. Drawdown-time (S-t) and residual drawdown-time (S'-t/t') Relationship for the pumping wells; A—Jacob method (Drawdown); B—Theis method (Recovery).
![]()
Figure A2. Lithological profile of the pumping wells.