MATHEMATICAL MODEL OF GROUNDWATER FLOW IN TEEB AREA , MISSAN PROVINCE , SOUTH OF IRAQ

A two-dimensional mathematical model is constructed to simulate the flow regime of the upper part of Quaternary Deposits in Teeb area, south of Iraq. The upper aquifer is mainly consist of gravel, sand, silt, and clay beds with good porosity. Their constituent grade from relatively coarse material near the apices, to the fine grained, nearly impermeable silts and clays in the Delta plain. The suggested conceptual model which is advocated to simulate the flow regime of aquifer is fixed for one layer, i.e. the activity of the deeper aquifer is negligible. The model is calibrated using trial and error procedure in two stages, steady state followed by transient state. According to the calibration process, the hydraulic characteristics of the upper aquifer has been identified, the hydraulic conductivity in the study area ranged (1,10) m/day, while the specific yield between (0.1, 0.4).


1.Introduction
A groundwater model is a computer-based representation of the essential features of a natural hydrogeological system that uses the laws of science and mathematics.It's two key components are a conceptual model and a mathematical model.The conceptual model is an idealized representation (ie. a picture) of our hydrogeological understanding of the key flow processes of the system.A mathematical model is a set of equations, that are subjected to certain assumptions and quantifies the active physical processes in the aquifer system(s) being modeled.
The equations that describe the groundwater flow may be solved using different types of models include analytical and numerical models.Because of the simplifications inherent with analytical models, it is not possible to account for field conditions that change with time or space, such as groundwater flow rate , groundwater direction, and other hydraulic properties.Also many analytical models require that the medium should be homogeneous and isotropic, for this reasons the use of numerical models is a very realistic situation (Wang and Anderson, 1982).Numerical models are capable of solving the more complex equations that describe groundwater flow.This model has been extensively used for groundwater analysis since the mid-1960's as high speed digital computers become widely available (Mercer and Faust, 1981).
In the present research a finite difference two dimensional model is used for modeling the groundwater flow for the upper aquifer in Teeb area by using MODFLOW program.Teeb area is located in north and north east of Missan province as shown in Figure 1.It occurs along the foot of mountains of the Iraqi-Iranian frontier in south of Iraq, between longitudinal-line ( 47°06'-47°36' ) and latitude-line ( 32°06'-32°30' ).The considered area is about 1860 km 2 .It extended from Teeb area close to the Iraqi-Iranian border to Shikh Fars area.This area represents the eastern edge of the Mesopotamian plain which is part of it.The shallow aquifer is composed of clay, silt, sand, and gravel which grade from relatively coarse material near the apices, to the fine grained, nearly impermeable silts and clays of the Delta plain.There is two rivers in this region, Teeb and Duiraige river.Teeb river crosses the area from north to south and ends in a local marsh.Duiraige river lies in south-east of this area.Marshes receive their water from the distributaries of Teeb and Duiraige river, and other small streams that generally flow toward the west and south ward from the foothills of Himreen along the Iraqi-Iranian frontier in south of Iraq.

Mathematical Model
The governing partial differential equation used in MODFLOW is (McDonald and Harbaugh, 1988).
, , are values of hydraulic conductivity along the x, y, and z coordinates axes, which are assumed to be parallel to the major axes of hydraulic conductivity (Lt -1 ).h is potentiometric head (L).W is a volumetric flux per unit volume and represents sources and/or sinks of water (t -1 ). is the specific storage of the porous material (L -1 ), and t is time (t).
In general , , , and may be functions of space ( , = , , etc.) and W may be a function of space and time ( ; equation (1) describes groundwater flow under nonequilibrium condition in a heterogeneous and anisotropic medium, provided the principal axes of hydraulic conductivity are aligned with the coordinate directions.
Analytical solutions of equation ( 1) are rarely possible for very simple systems, so various numerical methods must be employed to obtain approximate solutions.One such approach is the finite difference method, wherein the continuous system described by equation ( 1) is replaced by a finite set of discrete points in space and time, and the partial derivatives are replaced by terms calculated from the differences in head values at these points.The process leads to systems of simultaneous linear algebraic difference equations; their solution yields of head at specific points and times.

3.Conceptual Model
One specific but very important component of the modeling process is the conceptualization.The proper definition of the nature of a given hydrologic system is the corner stone in the selection of the applicable mathematical model.Fig. 2 shows the suggested conceptual model of the shallow aquifer of the Quaternary deposits in the study area.Two dimensions mathematical model is presented to simulate the flow regime in the shallow aquifer, in other word this simulation is fixed for upper layer only.It is assumed that the base of the upper aquifer (aquitrad) is an impermeable boundary, i.e. the activity of the deeper aquifer is negligible.There is no sufficient data for the static water level of the deeper aquifer is available, this is the main reason for selecting the present conceptual model.Also constructing the peizometeric system for monitoring the static water level of deeper aquifers is more costly.

4.Boundary Conditions and Model Grid
Fig. 3 shows a spatial distribution of the aquifer in the study area, the present model consists of (55) column and (50) row.Where the area of cells is equal to (1000 × 1000m).All boundaries in the present model was modeled as head-dependent boundary to allow inflow to the modeled region at a rate proportional to the head difference between the aquifer outside the simulated area and the model boundary.The top of the model was represented as unconfined aquifer.The water table elevation changes as part of the model solution.The bottom of the model was represented as a no flow condition.The vertical location of this boundary was selected to correspond with the base of the aquifer (aquitrad).
Shallow aquifer is composed of clay, silt, sand, and gravel Aquitrad Water table

Fig.3 Configuration of Nodal Network and Boundary
Condition of the Study Area.

Initial and Historical Water Table Elevations
The initial conditions refer to initial values of elements that may increase or decrease in the course of time inside the model domain and they cover largely the same phenomena as the boundary condition do.
Ten monitoring wells distributed over the study area are selected to measure the initial and historical groundwater level for one year.Fig. 1 shows the local distribution for observed wells.In the monitoring process used a number of different devices, such as level, tap, sounder and GPS to measure the groundwater level per month for the period (April, 2010;March, 2011).

Extracting Rate and Distribution of Pumped Wells
The spatial distribution of existing wells in the study area is shown in Fig. 4. The average discharge of the existing wells in the western parts of the study area ranged from 4 l/sec to 8 l/sec.While most existing wells in eastern part of the study are along the foothill of Hemrin abstracted groundwater from the Bai Hassan and Mukdadiya Formations by artesian wells.All hand-dug and tube wells in the study area withdraw groundwater for domestic and stock use only.Tube wells are commonly used for apstracting groundwater compared with large diameter hand dug wells.Hand dug wells are conducted randomly with non-uniform shapes.

Initial Assessment of Hydraulic Characteristics
Twenty soil samples are selected to obtain the soil texture class.The average value of effective porosity is then interpolated using Kriging techniques in Geostatistical analysis extension of Arc GIS 9.3 to produce the soil hydrologic group layer in the study area as shown in Fig. 5. Based on the covered area by soil hydrologic groups, and the data which is presented in Tables 1 and 2, the initial assessment of hydraulic conductivity and effective porosity are supplied to the numerical program.

Direct Recharge (Refined CMB Method)
Groundwater recharge may be described as the process whereby the amount of water present in, or flowing through, the interstices of the sub-soil increases by natural or artificial means.(Huisman and Olsthoorn, 1981).In the present model a refined CMB (chloride mass balance) equation is used, which is both mathematically and statistically sound; no arbitrary averages, such as the weighted average, are incorporated, as is the case in the classical CMB approach (Subyani and Sen, 2006).The perturbation methodology dictates that additional statistical parameters, such as the standard deviations and correlation coefficient of the basic variables, appear in the final formulation of the average recharge estimation by the CMB method.The equation of refined CMB method is Where: : recharge rate : average annual rainfall is the correlation coefficient between the rainfall and its Cl-concentration, and are respectively the standard deviation of rainfall and its Cl-concentration.And and are the weighted average Cl-concentration in rainfall and groundwater respectively.In order to apply the refined CMB method as presented in equation ( 2), mean monthly rainfall records (October, 2010-April, 2011) and samples of rainfall to measure Cl -concentration are collected.Fifteen groundwater samples for average Cl -concentration are also taken from the wells in the study area, see Table 3. Chloride concentrations in these samples are measured in the laboratory of Chemical Engineering Department/ Basrah University.There are no irrigation activities in this area, nor is there active dry deposition of Cl -.The annual rainfall average for the period (1980 -2007) is 178.2 mm.All the data needed for the application of the refined CMB method are given in Table 4.

Calibration of the Model
Calibration of present model was carried out in two sequential stages, a steady state calibration followed by transient calibration.Steady state calibration permits the adjustment of hydraulic conductivity, where aquifer storage changes are not significant.Thus, dynamic stress and storage effects are excluded.Comparison between the observed and calculated heads for this process is shown in Table 9.A trial and error calibration based on mean absolute error and root mean squared error as the following formulas:  After the steady state calibration was achieved, transient calibration was undertaken to calibrate the aquifer storage and direct recharge.The unsteady state calibration results were evaluated by comparing the temporal variation in simulated heads with those of observed ones at ten observed wells, see Figs. (6 to 15) and Table 10 shows the unweighted groundwater levels calibration statistic.Figs.16 and 17 show the finally distribution of calibrated hydraulic conductivity and specific yield respectively.Also the calibration process shows the best value of direct recharge as percentage from rainfall is equal to 10%.

Fig. 1
Fig.1 Location of Study Area and Observed Wells.

Fig. 2
Fig.2 Conceptual Model in the Study Area.

Fig. 4
Fig.4 Spatial Distribution of Existing Wells in the Study Area.

Fig. 5
Fig.5 Hydrological Soil Groups in the Study Area.

Fig
Fig.6 Comparison between Observed and Simulated Heads for well No. 1.

Fig. 13
Fig.13 Comparison between Observed and Simulated Heads for well No. 8.