Modeling Contamination Transport (Nitrate) in Central Basin Erbil, Kurdistan Region, Iraq with Support of MODFLOW Software

Abstract


Introduction
Since the water in rivers and lakes frequently derives from groundwater, a significant amount of the freshwater on the earth's surface is found underground.It is not exaggerated to say that all the water we use for drinking, agriculture, and industry is regarded as groundwater, or at least it was once groundwater in the natural water cycle (Huseen and Abed, 2020).Groundwater is a significant and influential factor in both the environmental and economic spheres.In dry months with little rain, they provide water to lakes, rivers, and wetlands (Huseen and Abed, 2020).Groundwater supplies in many areas are significantly threatened by excessive groundwater usage and a lack of knowledge of the hydrogeological context (Dhungel and Fiedler, 2016).Numerical simulation models have developed into potent tools in science and engineering during the last several decades for describing and researching phenomena and physical systems.Numerical simulation models describe hydrological processes, water movement in the subsurface, and contamination transport in groundwater hydrology.Simulators are also used to analyze or forecast the long-term effects of water withdrawals and pollutant migration and to investigate different groundwater management options (Karatzas, 2017).To model groundwater flow in an area, simulations require a variety of variables, including hydrogeological, hydrological, topographical, and climatic data (Rapantova et al., 2017).To estimate the derivatives of the algebraic equations, numerical groundwater flow models employing a finite difference method were initially introduced by Anderson (1984), Mercer andFaust (1984), andRemson et al. (1971), and the groundwater flow equation stated that the integrated finite difference (Narasimhan and Witherspoon, 1978).The USGS Modular Three-Dimensional Ground-Water Flow Model (MODFLOW), the U.S. Geological Survey, was created in the early 1980s, and it is now the world's most well-known 3D finite-difference groundwater flow model.The MODFLOW was the initial name of the program (Harbaugh et al., 2000, McDonald andHarbaugh, 2003).MT3DMS code uses to model solute transport and investigate the temporal and spatial distribution of contaminants (Qi et al., 2022).One of the most well-known 3D models in the world is the GMS Software (Aquaveo, 2007).This program provides a valuable tool to explore groundwater production and its behavior in the different soil layers and has high accuracy for analyzing groundwater movement (Ahmed, 2009;Panagopoulos, 2012;Huseen and Abed, 2020).Some papers studied the simulation of groundwater flow can summarize as follows: Zhu et al. (2000) simulated transport of petroleum in limestone karst fissures in the Dawu water source area in Zibo Area, Shandong Province, China.The result showed that the permeability coefficient and effective porosity significantly impact the solute transport in the aquifer (Yang et al., 2011).Visual MODFLOW was utilized for modeling groundwater flow in Tongliao, China.The hydrological data were used for building a conceptual model.The result shows a good correlation between the observed and calculated head.The result can use to model the contamination transport in the area.Mustafa et al. (2017) simulated groundwater flows to study the impact of changes in Tigris water surface elevation on layers of the aquifer for the Nuclear Research Center in the Al-Tuwaitha area, as well as to evaluate the ability of the proposed pumping well to collect groundwater and change the direction of flow by using processing MODFLOW V. 5.3.The result shows that surface water elevation increased the head of water in observation wells and the flow velocity 9.77 x10 -6 , 3.6 x 10 -4 , and 0.0015 m/d for three-layer, respectively.The study suggested methods to change the groundwater direction by using a pumping well-Simulated nitrate concentration in Qazvin plain using GMS.The result showed that high nitrate concentration due to sewage might lead to more contamination in the central part of the aquifer.The purpose of this study is to model groundwater flow and nitrate transport in Central Basin-Erbil City using GMS version 10.6.2 and predict the future effect of nitrate after 50 years from now.

Description of the Study Area
The study area was located in the central basin of Erbil toward the Greater Zab River, which covers about 579.72 km 2 , as shown in Fig. 1.Erbil dumpsite finds in the main basin.Many contaminants leach from it to the soil and can reach the groundwater.The digital elevation map for the study area is presented in Fig. 2. The geological formations of the study area are Bai Hassan (Upper Pliocene) and Mukdadiya (Lower Pliocene) sandstone, gravely sandstone, and red mudstone are found in fining upward cycles in the Mukdadiya formation (Jassim and Goff, 2006).Most of the research region is covered by the Bai Hassan Formation, made up of molasses sediments characterized by an alternating of claystone and conglomerates that are included with some sandstone and siltstones (Hassan, 1998).The groundwater aquifer in the study area is an unconfined aquifer that mainly consists of gravel, sand, silt, and clay.The geological formation map of the soil is shown in Fig. 3.

Mathematical background
The equation of groundwater flow is based on the conservation equation and Darcy's law (Bear and Cheng, 2010;Dellerur, 1999 andJavandel et al., 1984): (1) Where: : is specific storage of the porous material ( ), t: is the time (T), W: is the flux of volumetric per unit volume describing sources/sinks of water ( ), h: is the Hydraulic head (L), : are the hydraulic conductivity in x, y and z direction (L/T).
The equation covers solute transport in aquifer (Bear and Cheng, 2010;Dellerur, 1999, andJavandel et al., 1984) is: (2) Where: -C: groundwater pollutants concentration, t: time (t), vi: is the seepage velocity, qx: is the volumetric water flux per unit volume of the aquifer, Cs: is the sources or sinks concentration, Dij: is the coefficient of hydrodynamic dispersion, n: is the porous medium porosity, Rc: is the term of a chemical reaction.

Methodology
The study area covers about 579.72 km 2 .The data about the groundwater observation head was obtained from the general director of groundwater -in Erbil city; the data includes the observation head for 23 wells for the years 2021-2022.Groundwater Modeling System (GMS) software version 10.6.2 was used.MODFLOW is used to simulate the groundwater flow in the study under steady-state conditions; MODPATH is a particle tracking code used to study the direction of groundwater flow and tracking of particles.The model domain is divided into 500 columns and 156 rows with cell dimensions 100 m x100 m and three layers.The total cells of model 78000 with 57972 cells active and 20028 cells inactive.The initial model head is assumed as the default value of GMS software for the (first, second and third) layers are 434 m a.s.l., 134 m a.s.l., and 109 m a.s.l., respectively.The annual groundwater recharge for the central basin from 1980 to 2016 was 0.08782 m/year.For this study, the groundwater recharge was considered a constant and was applied to the study area's top layer (Al-Kubaisi et al., 2019).The boundary condition is constant from three sides (Top, bottom, and left) which vary between 210 and 342 m a.s.l, and the fourth side is the Greater Zab River (right side).The river conductance is calculated for this side.
River conductance in MODFLOW is defined as the river area (width times length) in the cell multiplied by the riverbed hydraulic conductivity materials divided by the vertical thickness of the riverbed materials (Groundwater Modeling System, 2019).

Soil layer
The aquifer elevation is 400 m, and the deepest well in the study area is selected for this study.The study area was divided into three layers.The first layer is silty sandy gravel with a 300 m depth.Thesecond layer is clay at 25 m and the third layer is silty sand at 75 m depth.The aquifer is unconfined.The lithology description of the deepest well is presented in Fig. 4.

Aquifer parameters
Hydraulic conductivity is defined as the ability of fluid to flow through pores.Each soil type has a range of hydraulic conductivity values.The hydraulic conductivity of silty sandy gravel soil, clay, and silty sand can be presented in Fig. 5.The porosity of soil layers is showed Table 1.The wells coordination and observation head are presented in Table 2. Longitudinal dispersive ranges from 10 -2 to 10 4 in this study are assumed to be 250 m (Gelhar et al., 1992), and transverse dispersive to longitudinal dispersity is 1/10.Diffusion is the process the solute can transport from an area of high concentration to an area of low concentration, the effective molecular diffusion coefficient results from a multiplying of the diffusion coefficient and tortuosity, which is a measurement of the impact of the shape of the water molecules' flow path in a porous medium, the effective molecular diffusion is very small and neglected if compare with the effect of chemical dispersion, the effective molecular diffusion coefficient for this study is considered to be zero (Chiag, 2012;Saatsaz et al., 2013;Khayyun, 2018).

Model Calibration
The model was calibrated by changing parameters such as hydraulic conductivity, porosity and recharge rate.Trial and error method was used to reduce the difference between the predicted and calculated groundwater head.The model repeatedly runs until the calculated head closely resembles the observed head with an acceptable degree of accuracy.Calibration was achieved when the interval of estimated error ±0.05 of the measured value.The best fit curve between the observed and calculated value is shown in Fig. 6. the result is shown the model was more susceptible to hydraulic conductivity than porosity, the calibration value of hydraulic conductivity for layers 1, 2 and 3 are 1.4 m/d, 0.0025 m/d, and 1.35 m/d, respectively, and the porosity for layers 1, 2 and 3 are 0.23, 0.095, and 0.37, respectively.The value of mean absolute error (MAE), mean error (MR), root mean square error, and R 2 equal 0.398044, 0.398044, 3.60265, and 0.9917, respectively.The groundwater contour map for layers 1, 2, and 3 are presented in Figs.7a, b, and c.The flow budget is essential to calculate the flow for the model as subregions and the budget the entire model area.The flow budget for the first layer is presented in Table 3, and the percent of discrepancy is 0.0000069, which is less than one; this means the equation of groundwater flow is correctly solved by MODFLOW.

Model Validation
The model was validated with data from the observation groundwater head for year 2021, which was obtained from the General Director of Groundwater.This data was applied to the model with the same hydraulic properties used in the calibration part.The correlation between the observed head and predicted head R 2 = 0.8387, as shown in Fig. 9.The objective of model validation is to check if MODFLOW Model work accurately with information about groundwater table for the same wells for different years and achieve possible significant trust MODFLOW Model.

MT3DMS model: Solute Transport
The result of the MODFLOW model used to simulate the contamination transport from landfill sites depends on the data for 2022.Nitrate contamination plume is predicted for the next (10, 30 and 50) years, as shown in Figs.11a, b, and c.MT3DMS model run with packages advection and dispersion, the results showed that dispersion package was affected more in the spread of contamination than advection because of the high permeability of the first layer, which contains silt, sand, and gravel soil, the nitrate spreads horizontally and vertically downstream.Landfill leachate can contaminate a 1.5 km upstream landfill and 1.5 km downstream landfill; the nitrate contamination can reach the second and the third layer after 50 years with a high concentration.Figs.12a and 12b are presented the nitrate contamination plume in the second and third layer after 50 years.

Conclusions
Groundwater Modeling System is a powerful software for modeling groundwater flow and contamination transport.The groundwater flow model is highly affected by hydraulic conductivity in the calibration step; a good correlation between the observed head and calculated head is obtained as R 2 = 0.9917.The processing effect on nitrate transport in the aquifer is dispersion, and the advection process has less impact.The landfill leachate can contaminate about 1.5km upstream and downstream of the landfill; the nitrate contamination was transported vertically and reached the third layer.The effect of landfill leachate was not reached the Greater Zab River only affected wells in the surrounding area.

Fig. 3 .
Fig.3.Geological formation map of soil for study area, modified after(Sissakian and Fouad, 2014) conductance per unit length [(L 2 /T)/L] or [L/T], t: the thickness of the river bed, the material [L], w: the width of the material along the length of the river [L], L: the length of along river reach, K: hydraulic conductivity of river bed sediment.The hydraulic conductivity of river bed sediment is equal to 1*10 -6 m/s = 0.0864 m/d.The thickness of river bed sediment is equal to 4 m(Seeyan and Merkel, 2015), and the width of Greater Zab River within the study area is a minimum width of 56.07 m and a maximum width of 349.82 m.The average width is 136.1969m, and the river length within the study area is 19.73 km.In this study, the river conductance= 2.94 (m 2 /d)/m according to Equation 3. The river water surface elevation varies between 219.95 to 242.92 m a.s.l., which was obtained from the Directorate of Irrigation Erbil.MT3DMS code is used to model groundwater contamination.The nitrate contamination leached from the Erbil landfill site as a point source contamination is considered a constant concentration; nitrate concentration in landfill leachate is 5264 mg/l.The properties of landfill leachate were tested in Erbil Environmental Office.The initial nitrate concentration for the model was selected to be like the World Health Organization (WHO) standard 50 mg/l (WHO, 2004).

Fig. 4 .
Fig. 4. Lithology description of deepest well in the study area (General Director of Groundwater)

Fig. 6 .
Fig. 6.Observed groundwater head vs. calculated groundwater head at steady state condition

Fig. 7a .
Fig. 7a.Groundwater flow contour map at steady state condition in study area layer 1, b. Groundwater flow contour map at steady state condition in study area layer 2, c.Groundwater flow contour map at steady state condition in study area layer 3.

Fig. 9 .
Fig. 9. Model validation best fitting curve between observed and calculated head

Fig. 10 .
Fig. 10.Partials Flow Path from the source of contamination

Fig. 11a .
Fig. 11a.Nitrate Plume after 10 Years in Study area, b.Nitrate Plume after 30 Years in Study area, c.Nitrate Plume after 50 Years in Study area.

Fig. 12a .
Fig. 12a.Nitrate contamination Plume in the second layer after 50 Years, b.Nitrate contamination Plume in the third layer after 50 Years

Table 1 .
The porosity of soil layers.

Table 2 .
Groundwater observation head for year 2022

Table 3 .
Model flow budget for the first layer