Groundwater Flow Modeling for Qushtapa Plain Unconfined Aquifer in Southern Erbil Basin, Kurdistan Region, Iraq

Increasing population growth and water demand for various purposes such as irrigation, domestic and industrial production in many parts of the Kurdistan Region is causing deficit in fresh water and rising groundwater dependence. Drilling many deep wells in the area unsystematically and continuously increased pumping water from groundwater reservoirs results in lowering of water table. Therefore, it is essential to assess the management of water resources. The study focuses on the groundwater modeling for the Qushtapa District plain area in particular under steady state flow conditions. The aquifer was simulated under unconfined condition and is represented by a single layer of 100 m thickness. MODPATH was used to measure contamination track lines and travel times. This approach involved the introduction of particles at sources of contaminants in the wells and the recharge area, then the identification of the path lines and the determination of the special distribution of contaminants through steady state flow conditions. The simulation of the groundwater head shows that the groundwater head starts from the northeastern part of the plain and decreases towards Lesser Zab River in the south of the plain from 420 m to 140 m above sea level. The modeled layer was calibrated under steady state conditions using hydraulic parameters obtained from observation and pumping wells. The calibrated model is effective in producing steady-state groundwater head distribution and good compliance with observed data. The standard error was estimated as 4.88 m, the normalized root mean square error is 8.3% and the residual mean is 15.79 m. The results of the forward tracking show the source of potential pollutants from the recharge area after different travel time, the particles released at the northern boundary travels to the center and the western part toward the pollution sources. The results of the backward tracking show that the particles located in the


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.

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

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

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

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

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.

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

Conceptual Model
Conceptual model construction includes the definition of the basin boundaries, aquifers recharge, and discharge sources (Anderson & Woessner, 1992

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 (Shwani, 2008).

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

S. Seeyan
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.

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.

Model Calibration
Model calibration was performed for the twenty two observation wells in the area,

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

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.

S. Seeyan
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.