Hindawi Journal of Chemistry Volume 2017, Article ID 2041648, 16 pages https://doi.org/10.1155/2017/2041648 Research Article Regional Groundwater Quality Management through Hydrogeological Modeling in LCC, West Faisalabad, Pakistan AamirShakoor,1ZahidMahmoodKhan,1MuhammadArshad,2HafizUmarFarid,1 MuhammadSultan,1MuhammadAzmat,3MuhammadAdnanShahid,4andZafarHussain5 1DepartmentofAgriculturalEngineering,BahauddinZakariyaUniversity,Multan,Pakistan 2DepartmentofIrrigationandDrainage,UniversityofAgriculture,Faisalabad,Pakistan 3Geo-InformaticsEngineering,NationalUniversityofSciencesandTechnology(NUST),Islamabad,Pakistan 4WaterManagementResearchCenter,UniversityofAgriculture,Faisalabad,Pakistan 5DepartmentofForestryandRangeManagement,BahauddinZakariyaUniversity,Multan,Pakistan CorrespondenceshouldbeaddressedtoHafizUmarFarid;farid [email protected] Received 8 May 2017; Revised 11 August 2017; Accepted 23 August 2017; Published 19 October 2017 AcademicEditor:PedroA´vila-Pe´rez Copyright©2017AamirShakooretal.ThisisanopenaccessarticledistributedundertheCreativeCommonsAttributionLicense, whichpermitsunrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperlycited. Theintensiveabstractionofgroundwateriscausinganumberofproblemssuchasgroundwaterdepletionandqualitydeterioration. Tomanagesuchproblems,thedataof256piezometersregardinggroundwaterlevelsandqualitywereacquiredfortheperiodof 2003to2012incommandareaofLowerChenabCanal(LCC),WestFaisalabad,Pakistan.MODFLOWandMT3Dmodelswere calibratedfortheperiodof2003–2007andvalidatedforyears2008–2012withrespecttoobservedgroundwaterlevelsandquality data,respectively.Afterthesuccessfulcalibrationandvalidation,twopumpingscenariosweredevelopeduptoyear2030:Scenario I(increaseinpumpingrateaccordingtothehistoricaltrend)andScenarioII(adjustedcanalwatersuppliesandgroundwater patterns).ThepredictedresultsofScenarioIrevealedthat,uptoyear2030,theareaundergoodqualitygroundwaterreduced significantlyfrom50.35to28.95%,whilemarginalandhazardousgroundwaterqualityareaincreasedfrom49.65to71.06%.Under ScenarioII,thegoodqualitygroundwaterareaincreasedto6.32%and12.48%areapossesseslesshazardousqualityofgroundwater. Itwasconcludedthatthecanalwatersupplyshouldshiftfromgoodqualityaquiferzonetopoorqualityaquiferzoneforproficient managementofgroundwateratthestudyarea. 1.Introduction wells of capacity ranging from 0.015 to 0.056m3/sec have beeninstalledinPakistanforirrigationatdepthof30–85m In Pakistan, groundwater is the second largest source of having diameter of 15–30cm, out of which 86% tube wells irrigatedagriculturebecause ofthearidclimaticconditions are installed in Punjab [6]. The farmers of tail ends of the [1,2].Therelianceongroundwaterhassignificantlyincreased distributariesandwatercoursesareforcedtorelyheavilyon during the past two decades to meet the food and fiber groundwater,particularlywherecanalwatersuppliesarecon- requirements of growing population [3]. Furthermore, the strained.Withoutgroundwateravailability,notonlyPunjab surfaceirrigationsystemhashighdegreeofconveyanceand butthewholecountrywouldfaceaseverewatershortagethat application water losses. Thus, the system operates at an leads to food shortages as Punjab produces more than 90% efficiencyoflessthan40%thatmeetsonlylessthan40%ofthe of total grains [7]. However, the unchecked abstraction of cropwaterrequirements[4].Theprovisionofirrigationwater groundwaterhascreatedseriousnegativeconcernsinterms from other sources such as groundwater is indispensable ofloweringwatertableandsaltwaterintrusionthatmaylead for potential productivity [5]. The growth rate of private totheissuesofsustainabilityofusablegroundwaterresources. tube wells has increased significantly at the rate of 60% in Asaresultofsalinegroundwaterintrusion,about200public Punjabprovincefrom1991to2000.About1.20milliontube tubewellsinitiallyinstalledinthefreshgroundwaterzoneof 2 JournalofChemistry PunjabandSindhprovinceshadbeenabandoned[8].Shah conceived. Kori et al. [17] linked groundwater flow model [9]reportedsimilarproblemsexistedinmostoftheirrigated (MODFLOW)withsolutetransportmodel(MT3D),andsev- regions of the world; those were further increasing rapidly eralsimulationrunswerecarriedoutaftersuccessfulcalibra- andnegativelyaffectingagriculturalproductivity. tionofmodelfortwosampledsiteslocatedatJRS-57andJRS- Therearewidespatialandtemporalvariationsinground- 60 tube wells at Nawabshah-Pakistan. Carretero et al. [24] waterqualityintheIndusbasin,whichisduetothepatternof appliedgroundwatermodelingtechniquetostudytheimpact groundwatermovementintheaquifer[10].Thebeltoffresh ofapossibleriseof1minsealevelagainstthelow-lyingcoast groundwaterisgenerallyavailablenearthemainriversand ofPartidodeLaCosta,Argentina.Similarly,Abu-el-Sharand canalonaccountofhighrechargeoffreshseepagewater.But, Hatamleh[25]simulatedthegroundwatermodel(PMWIM) then,thequalityofgroundwaterchangestounfitaslaterally for the Azraq Basin, Jordan, to manage groundwater. The awayfromtherivers[5,11].Thecontinuousandintensiveuse largenumberoflatestresearchpapers—Asokaetal.[26],Gal- ofgroundwaterforirrigationisaddingplentyofsaltscausing itskayaetal.[27],Kambaleetal.[28],AbdullahandMorteza secondary salinization because groundwater generally has [29], and Durand et al. [30]—dealing with groundwater moresaltsthancanalwater.BakhshandAwan[12]reported managementusingmodelingtechniqueswerepublishedisa thatapplicationofgroundwaterhavingtotaldissolvedsolids testimonytotheimportantroleplayedbythemodels. −1 (TDS)of1000mgL uptoasoildepth370mmchangedthe Therefore,theassessmentofgroundwaterqualitythrough top300mmdepthofnonsalineintoasalinesoilthatimpaired modeling is important to ensure the safe use of these crop productivity. In many agricultural areas of Pakistan, resources on a sustainable basis. If unchecked pumping of the usage of poor quality tube well water for irrigation groundwater continues to remain in the study area, the is considered as one of the major causes of salinity and, irrigation tube wells would not be able to lift water at their consequently,lowerfoodproductivity[13].Accordingtoan presentquality.Hence,thereisadireneedtoinvestigatethe estimate in Pakistan, the secondary salinization degraded impactofgroundwaterflowconditionsandoverexploitation the crop land which reduced the production potential of on groundwater quality. Thus, the followings were the two majorcropsby25%,valuedatanestimatedlossofUS$250 mainobjectivesofthecurrentstudy:(1)todeveloparegional million/year[14].Lowqualitywaterandsoilsalinitycanaffect hydrogeological groundwater flow model and observe its plantgrowthandsoilstructureinseveralways,directlyand futuretrendand(2)toformulategroundwatermanagement indirectly[15]. strategyforitsproficientutilization. Thus,theassessmentofgroundwaterqualityisimportant toensurethesafeuseoftheseresourcesonasustainablebasis. 2.MaterialandMethods Groundwatermodelsaremostwidelyusedtoolsforefficient managementofpreciousgroundwaterresourcesandtopre- 2.1. Study Area. The research study was carried out in the dict different future scenarios [16]. Different groundwater commandareaofLowerChenabCanal(LCC),WestFaisal- modelingcodesareavailable,eachwiththeirowncapabilities, abad, Pakistan, having longitude of 73.85∘ to 72.18∘ (E) and operationalcharacteristics,andlimitationssuchasPMWIN, latitudeof32.32∘ to30.85∘ (N)(Figure1).Thesitehasgross FEFLOW, SVFlux, and GWVistas. But the most extensively commandareaof1.16Mhaandculturablecommandareaof usedthree-dimensionalgroundwaterflowmodelamongthe 0.98Mha, respectively. It falls in rice-wheat agroecological available models is PMWIN (Processing MODFLOW for zone of province Punjab, Pakistan. It is comprised of vast Window) [17, 18]. Its popularity has continued, in part due canalnetworksuchasmaincanals,branchcanals,andminor to the modularity of the program, resulting ability, and distributaries.TheRiverChenabandGugeraBranchCanals user friendly interface [19]. It uses a block-centered finite- arelocatedatthenorthwestandsoutheastside,respectively. difference scheme for saturated zone. The advantages of Similarly, Qadirabad-Baluki Link Canal is located in the PMWIN include numerous facilities for data preparation, northeastsideandTrimmu-SidhnaiLinkCanalispresentin easyexchangeofdatainstandardform,extendedworldwide thesouthwestsideofthestudyarea(Figure2). experience, continuous development, availability of source code, and relatively low price or being freely available 2.2.Climatic. ThesummerseasonstartsfromAprilandcon- [20]. A number of research studies have been conducted tinuestillOctober.Thetemperaturevariesintherangeof21to regardingapplicationofmodelingapproachforgroundwater ∘ 51 Cduringthesummerseason.Similarly,thewinterseason managementindifferentpartoftheworld.Moecketal.[21] starts from October and lasts till April with temperature developed a 3D groundwater model for simulating existing ∘ ranges of 1 to 27 C. The average annual precipitation was andproposedwatermanagementstrategiesasatooltoensure theutmostsecurityfordrinkingwaterinBasel,Switzerland. estimated to be 439mm. The ten-year average value of ET𝑜 was1413mm/year. Gebreyohannes et al. [22] developed regional groundwater flowmodelforGebaBasin,northernEthiopia,andreported thatnoneofthehydrogeologicalformationscanbeexploited 2.3.HydrogeologicalConditions. Lithologyofaquifersystem for large-scale groundwater exploitation. Rahmawati et al. in the study area has different classification according to [23] conducted research to study saltwater intrusion from differenttexturalcharacteristics.Thesurfacesoiltexturesare 1995 to 2108 in Semarang city based on well log data and largelyfineandmoderatelymediumwithgoodpermeability MODFLOW numerical model was used. For salt intrusion properties. The aquifers of study area were formed as a projectioninthefuture,thesealevelriseprojectionalsowas result of sediment deposition due to flat topography. These JournalofChemistry 3 N 35∘00N 30∘00N 25∘00N 0 75150 300 450 600 (km) 65∘00E 70∘00E 75∘00E Study area Punjab Pakistan N Hafizabad 32∘00N Pindi Batiyan Chiniot Safdarabad Bhowana Shahkot 31∘300N Faisalabad Jhang Gujra 31∘00N TT Singh 0510 20 30 40 (km) 72∘300E 73∘00E 73∘300E 74∘00E Main city Piezometers Study area Figure1:Geographicallocationofstudyarea. 4 JournalofChemistry N Qadirabad headworks Q B L 32∘00N in k 31∘300N h CChaenanlab River Branch Canal lanaC .rB ilA naiMGugera CaBnarl. Canal Branc Rakh pper g U n 31∘00N TrimhmeauTSdworks Jha L ow er Gu g era Br. C a n al L in k C 0510 20 30 40 a na (km) l 72∘00E 72∘300E 73∘00E 73∘300E 74∘00E Irrigation network Chenab River Canal Figure2:Canalnetworkofthestudyarea. sediments were carried by the river waters from the vast Shakoor et al. [6]. The northern and northwestern part of alluvialbasinofriversasaresultofmaterialswasheddown thestudyarea(alongChenabRiver)containedgoodquality from Himalayan Mountain. In general, like other aquifers groundwater.Themarginalandhazardousqualitywaslaidin in Punjab, the aquifers of study area are unconfined. The thesouthernpartofthestudyarea. groundwaterlevelsintheyear2012werepresentintherange The annual water balance of the whole study area from of145to205m(Figure3).Thegroundwaterlevelwashigher theyear2003totheyear2012isshowninFigure5.Thewater intheupperpartandgraduallydecreasedtowardsthelower balancedescribesthevolumeofwaterentering,subtraction, side,whichmeansthatthegroundwaterflowwasfromupper and net storage in the aquifer system. The water entering tolowerside. parametersarerechargeandriverleakage,whilesubtraction Thegoodqualitygroundwaterispresentinthenorthern parameters are wells and ET. The positive storage term partwhereaspoorqualitygroundwaterisinthesouthernpart showedthatvolumeofwaterextractedfromtheaquiferwas ofthestudyarea.Theaveragespatialvariationofmeasured morethantherechargewhichmeanscropwaterdemandwas groundwaterTDSintheentirestudyareafortheyear2012is fulfilled from aquifer storage and vice versa. The storage in giveninFigure4.AcalibratedTDSmeter(TDS-98302)was the aquifer in 2003 and 2012 was −444.5 and 326.6MCM, usedtomeasuretheTDSofgroundwaterwithanaccuracyof respectively, clearly indicating the replenishment of aquifer ±2%.ThecalibrationofTDSmeterwasaccomplishedinthe in 2012. Similarly, the volume of water withdrawn through standardsolutionintheWaterQualityandEnvironmentLab- tube wells increased from 3298 to 3759MCM from years oratory,DepartmentofIrrigationandDrainage,Universityof 2003to2012,respectively,alsoshowingincreaseddemandof Agriculture,Faisalabad.ThemeasuredTDSvaluesindicated groundwater. thatalmosthalfoftheaquiferofthestudyareahasmarginalto hazardousqualityofgroundwater.A30kmwidestripalong 2.4. Model Selection. PMWIN 5.3 (Processing MODFLOW theChenabRiverhadgroundwaterofgoodquality,whereas for Window) was used in this research. It is a simulation hazardous quality groundwater existed in 20km wide strip system based on the modular three-dimensional finite- at a distance of 40km away from the Chenab River. The difference technique for modeling groundwater flow and groundwaterqualityzonesoffreshandsalinewaterwerefor- contaminationingroundwaterwithawiderangeofnatural mulatedaccordingtocriteriadevelopedbyWAPDA[31]for systems. PMWIN is used widely throughout the world and irrigation water quality (good: TDS < 1000mg/L; marginal: it was applied in many groundwater modeling applications TDSof1000–1700mg/L;andhazardous:TDS>1700mg/L), [17,23,25]. JournalofChemistry 5 N 2 0 5 32∘00N 5 0 2 200 185 190 195 0 8 1 31∘300N 5 6 0 5 1 7 7 0 1 1 6 1 5 5 1 0 5 1 31∘00N 5 14 0 510 20 30 40 (km) 72∘300E 73∘00E 73∘300E 74∘00E Groundwater head Contour 2012 Figure3:Spatialvariationinmeasuredgroundwaterlevelsfortheyear2012. 2 2.5. Model Input and Calibration. The MODFLOW and 6800km irrigatedagriculturalregionlocatedalongtheSea MT3D models are packages of PMWIN 5.3. They were ofCortezinSonora,Mexico. calibrated with respect to observed groundwater level and Lithology of aquifer system in the study area was quality data, respectively, using inverse modeling method. obtainedfromtheWaterandPowerDevelopmentAuthority Thedataof256piezometersregardinggroundwaterleveland (WAPDA) [31]. The soils have different classification quality were acquired from the Department of Land and accordingtodifferenttexturalcharacteristics.Thesurfacesoil Reclamation,Faisalabad,fortheperiodfrom2003to2012to texturesarelargelyfineandmoderatelymedium,withgood achievethemodelcalibrationandvalidation(Figure1).The permeabilityproperties.Theareasof4451.3(38.48%),4987.3 piezometerwellshavemaximumdepthsof55m.Thedepth (43.08%),1621(14%),464.1(4%),and52.3km2 (0.45%)have ofpiezometersvariedatsomelocationdependinguponthe soiltextureoffine,moderatelymedium,medium,moderately watertableconditions.Thepiezometerwellshavediameters course,andcourse,respectively.Theaquiferofthestudyarea of5cmandstrainerlengthsof3m. wasdefinedwithfourdifferentlayersdependingupontheir The model area had a square geometry and the whole lithologicaldata[31].Thespatialdomainrepresentedinthe area was divided into 76 columns and 76 rows. The total model consisted of four layers (0–7, 7–30, 30–90, and 90m number of cells was of 5776 cells. Each square cell has a tobedrock).Thehorizontal(𝐾ℎ)andvertical(𝐾V)hydraulic dimension of 2.5km × 2.5km. The model area had 3653 conductivitieshavelargevariationfromonetotheotherside inactivecells,whichwereoutsidetheboundaryofthestudy of the study area (Table 1). The minimum and maximum area,whilethemodelareahad2123activecellslocatedwithin values of 1–265m/day and 1–15m/day, respectively, for theboundaryofstudyarea.Thecellsizeof2.5km× 2.5km horizontal and vertical hydraulic conductivities were used 2 (6.25km )wasalsousedforgroundwatermodelingstudies in the model, as majority of the tube wells were installed byAl-Fatlawi[32]inUmmErRadhuma,theWesternDesert, in second layer because 80% part of study area has 𝐾ℎ in Iraq, and by Khan et al. [33] in Rechna Doab. Abu-el-Shar the range of 70–100m/day. The similar range of values for and Hatamleh [25] developed groundwater model for the hydraulic conductivity within Punjab province domain was AzraqBasinandthebiggestcellsizeof8.69km2wasselected. used by Jehangir et al. [35]; Ahmad [36]; Arshad [37]; and Similarly, Schoups et al. [34] used cell size of 2km × 2km Khanetal.[33].Thespecificstoragevaluesforlayer1ranged −1 to calibrate groundwater model of the Yaqui Valley, having as0.0001–0.001m andfortheremaininglayersthevalues 6 JournalofChemistry Table1:Lithologicaldataofalllayers. Subsurface 𝐾ℎ % 𝐾V % Specific−s1torage % Specificyield % layers (m/day) area (m/day) area (m ) area area 1–15 20 >0-1 80 1𝐸−4–5𝐸−4 40 0.05–0.013 15 (Layer1) 15–70 60 1–5 10 5𝐸−4–6𝐸−4 25 0.013–0.17 60 70–160 20 5–15 10 6𝐸−4–1𝐸−3 35 0.17–0.25 25 1–70 10 >0–2 50 1𝐸−5–2𝐸−4 5 0.05–0.15 10 (Layer2) 70–100 80 4–10 25 3𝐸−4–1𝐸−3 15 0.2–0.25 80 100–165 10 10–15 25 2𝐸−4–3𝐸−4 80 0.15–0.2 10 20–80 20 >0–2 50 1𝐸−5–2𝐸−4 5 — — (Layer3) 80–120 30 4–10 25 3𝐸−4–1𝐸−3 15 — — 120–265 20 10–15 25 2𝐸−4–3𝐸−4 80 — — 20–80 20 >0–2 50 1𝐸−5–2𝐸−4 5 — — (Layer4) 80–120 30 4–10 25 3𝐸−4–1𝐸−3 15 — — 120–265 20 10–15 25 2𝐸−4–3𝐸−4 80 — — 𝐾ℎ:horizontalhydraulicconductivity.𝐾V:verticalhydraulicconductivity. N 32∘00N 31∘300N 31∘00N 0 510 20 30 40 (km) 72∘300E 73∘00E 73∘300E 74∘00E Groundwater TDS 2012 Good Marginal Hazardous Figure4:SpatialvariationsofmeasuredgroundwaterTDSintheyear2012. −1 of0.00001–0.0003m wereused.Thevaluesofspecificyield Pakistan, Kharif and Rabi. Kharif starts from June and July forlayer1rangedas00.05–0.25andfortheremaininglayers and goes to October and November, while the Rabi season the values of 0.05–0.20 were used. The effective porosity of starts from September and October and continues to April 0.25 was given to all layers [31]. The simulation time unit andMay. “days”andsimulationflowtype“transient”wasselected.Two The cumulative evapotranspiration during a period stressperiodsineachyearwereconsideredtorepresentthe (Kharif or Rabi) was divided by its duration in days and KharifandRabiseasonshaving183and182days,respectively, thus the evapotranspiration rate per day was calculated with six time steps in each stress period, as there are two usingCROPWAT.Theevapotranspirationratesof0.006and maincroppingseasonsbasedonagroclimaticconditionsin 0.003m/day were used for odd (Kharif) and even (Rabi) JournalofChemistry 7 4000 volume of aquifer representing sources (positive) and sinks 3000 (negative)[T−1],𝐶𝑠 istheconcentrationofsourcesorsinks 2000 [ML−3], ∑𝑁𝐾=1𝑅𝑘 is the chemical reaction term [ML−3T−1], M) 1000 and𝑡istime[T]. C M 0 The calibration of MT3D was achieved by adjusting me ( −1000 the values of advection and dispersion values. Advection olu −2000 depictsmasstransportbasicallyduetothebulkflowofwater V −3000 in which the mass is dissolved or movement of solute as −4000 a result of groundwater flow. Advection is the dominant process of solute transport due to groundwater flow. For a −5000 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 given time step, the movement of solutes due to advection Time (year) was simulated using the corresponding field velocity which was computed by the calibrated MODFLOW model and Storage ET theeffectiveporosityoftheaquifer.MXPARTofmaximum Wells Recharge 1300000particleswasallowedinasimulationwhilethevalue River leakage PERCEL(Courantnumberornumberofcellsinanyparticle) Figure5:Waterbalancefortheyears2003to2012. wasselectedas0.79.Itwasallowedtomoveinanydirection within one transport step and generally ranged from 0.5 to 1[39].Thedefaultvaluesoftheremainingparameterswere stress periods, respectively, throughout the model stress usedforMT3Dmodel. periods for model calibration. The recharge package was Thelongitudinaldispersivity(𝛼𝐿)andtransversedisper- used to simulate spatial distribution of recharge water at sivity (𝛼𝑇) are required to estimate dispersive transport of fieldfromrainfallandirrigationtothegroundwatersystem. thesolute.Thedispersivityisacharacteristicpropertyofthe Theminimumandmaximumrechargevaluesfrom0.00065 porous medium that describes the spreading of the solute to 0.0013m/day were used for odd (Kharif) stress periods, in a medium. Longitudinal dispersivity is used when the whereas values from 0.00026 to 0.0005m/day were used spreadingofsoluteisinthedirectionofbulkflow,whereas for even (Rabi) stress periods. The river package was used thetransversedispersivityistheperpendicular(verticaland to simulate the flow between an aquifer and surface-water horizontal) spreading of solute to the direction of bulk features such as rivers, canals, lakes, and reservoir. The flow. Thelongitudinaldispersivity(𝛼𝐿)wascomputedfrom moleculardiffusionandlocaldispersioncoefficientbasedon hydraulic features, such as canal length, bed width, full- heterogeneityofthemedium[40,41]asgivenin supply level, hydraulic conductance, and discharge of the canals of the study area, were used. The recharge flux was 𝐷 𝐷 computed by the model through the canal system. The net 𝛼𝐿 =( V𝑑 + V𝐿)+𝐴𝐿, (2) groundwater demand was calculated by subtracting the net canalwatersuppliesfromnetcropwaterrequirement.After where𝛼𝐿islongitudinaldispersivity(cm),𝐷𝑑iscoefficientof providing all the input data, the calibration process was moleculardiffusion(16×10−6cm2/s),𝑉isvelocityofflow(6.8 started.Inthiscalibrationprocess,inputparameterssuchas × 10−6cm/s), 𝐷𝐿 is dispersion coefficient (30 × 10−6cm2/s), recharge, pumping rate, and hydraulic conductivities were and 𝐴𝐿 is asymptotic longitudinal dispersivity (500cm). adjustedinMODFLOW. Asymptoticlongitudinaldispersivityisduetoheterogeneity of the medium, mainly dependent on the variance of the 2.6. Solute Transport Model. The MT3D numerical model log transformed conductivity and correlation length in the withPMWIN5.3userinterfacewasusedtosimulateadvec- mean direction of flow. As no field data were available tion and dispersion contaminants in the three-dimensional for solute transport in aquifer media, the value of 𝐴𝐿 was groundwater flow system of the study area. The MT3D adopted as 500cm after [42]. The Peclet number, 𝑃𝑒 = simulated groundwater flow by numerically solving the Δ𝑥/𝛼𝐿, was calculated to measure the relative significance groundwater flow and solute transport equations. The par- of advection and dispersion in the study area, where Δ𝑥 tial differential equation describing the three-dimensional is the size of a cell equal to 2500m. The unitless Peclet transportofdissolvedsolutesinthegroundwater[38]canbe number is usually used to decide the dominant factor in writtenin solutetransportamongadvectionanddispersion[43].After thesuccessfulcalibrationofbothmodels(MODFLOWand 𝜕 𝜕𝐶 𝜕 𝑞 𝑁 𝜕𝐶 𝜕𝑥 (𝐷𝑖𝜕𝑋 )− 𝜕𝑥 (V𝑖𝐶)± 𝜃𝑠𝑐𝑠+∑𝑅𝑘 = 𝜕𝑡, (1) MyeaDr3sD20)0f8o–r2y0e1a2rasn2d0t0w3–o2f0u0tu7,rethsecemnaordioeslswvearleidsaimteudlaftoerd.the 𝑖 𝑖 𝑖 𝑘=1 where 𝐶 is the concentration of contaminants dissolved in 2.7. Future Scenarios. The model prediction was accom- groundwater[ML−3],𝑥𝑖 isthedistancealongtherespective plishedinordertoinvestigatetheresponseofthemodelfor Cartesiancoordinateaxis[L],Diisthehydrodynamicdisper- twofuturescenariosregardingthepumpingrateuptoyear sioncoefficient[L2T−1],V𝑖istheseepageorlinearporewater 2030.Itwasassumedthattherewillbenouncertainchange velocity [LT−1], 𝑞𝑠 is the volumetric flux of water per unit or tragedy in climate and irrigation system. In Scenario I, 8 JournalofChemistry Table2:Statisticalanalysisofthefieldandmodeleddata. 220 210 Results m) 200 Parameters Formula MOD(mFL)OW M(mTg3/DL) ads ( 118900 1 𝑛 he Meanerror ME= 𝑛∑𝑖=1(ℎ𝑜−ℎ𝑠) −1.10 19.5 deled 116700 Mabesaonluteerror MAE= 𝑛1∑𝑛 󵄨󵄨󵄨󵄨ℎ𝑜−ℎ𝑠󵄨󵄨󵄨󵄨 1.72 139 Mo 114500 𝑖=1 1 𝑛 130 Rsqouoatrmeeerarnor RMSE=√𝑛∑(ℎ𝑜−ℎ𝑠)2𝑖 2.24 184 130 140 150 160 170 180 190 200 210 220 𝑖=1 ∑𝑛 (ℎ −ℎ)2 Field observed heads (m) Meffiocdieelncy MEF= ∑𝑖𝑛=1(ℎ𝑜 −ℎ𝑠)2 0.98 0.91 Figure 6: Scattergram of measured versus modeled heads in 𝑖=1 𝑜 MODLOW. (𝑅2) 0.89 0.87 3600 thepumpingwillincreaseaccordingtothehistoricaltrend. g/L) 3100 As in Punjab province of Pakistan, about 1.20 million tube m 2600 wellshavebeeninstalledandincreasingat5.5%annually[44]. DS ( 2100 T This indicated that the amount of withdrawal from aquifer ed 1600 may increase from 3738MCM in year 2012 to 6008MCM el d in year 2030 whereas recharge would be about 2664MCM Mo 1100 according to historic trend. In Scenario II, one of the water 600 management options for the study area was proposed. In 100 theupperpartofthestudyarea(PindiBhatiyan-Safdrabad) 100 600 1100 1600 2100 2600 3100 3600 wheregroundwaterhasgoodquality,therateofgroundwater Field observed TDS (mg/L) abstractionwasincreasedandrechargethroughanirrigation Figure7:ScattergramofmeasuredversusmodeledTDSinMT3D. systemwasdecreasedby35%,whileinthelowerpart,where groundwater has poor quality, the groundwater abstraction wasdecreasedandrechargethroughanirrigationsystemwas increasedbythesamerate. and simulated values. The average coefficient of determina- tion (𝑅2) for the selected piezometers for calibration and validation period was calculated as 0.89 and 0.87 (Figures 3.ResultsandDiscussion 6 and 7) for MODFLOW and MT3D models, respectively, which indicated a close agreement between the calculated 3.1. Calibration and Validation of Models. The degree of and measured groundwater head. Hagos [49] calibrated fit between model simulations and field measurements was PMWIN model for Raya valley, Ethiopia, with respect to quantified by statistical means (Table 2) and all parameters groundwater level. The values of ME, MAE, and RMSE of were found in acceptable range. The minus sign of mean calibratedresultswere−1.4m,7.8m,and10.7m,respectively, error(ME)representedthatthemodelsimulatedvalueswere withcoefficientofdeterminationof0.97andreportedtobe higherthanthemeasuredhead.Thecalibrationcriterionfor satisfactory. hydraulic head is root mean squire error (RMSE) which is lessthanorequalto10%ofheadvariationwithintheaquifer 3.2.SensitivityAnalysis. Thesensitivityanalysisshowedthat being modeled [45] and the same criterion is also followed rechargeandtransmissivityoftheaquiferweremostsensitive in groundwater TDS. The head in the aquifer within the parameters. The factors of 0.5, 0.8, 0.9, 1.1, 1.2, 1.3, and studyareavariesfromapproximately145to205m,resulting 1.5 were multiplied with the calibrated values of recharge in an acceptable RMSE of 6m or less. Similarly, the TDS and transmissivity.The resultinghydraulic heads were then in the aquifer at selected points are varied from 384 to comparedwiththeobservedheadsandRMSEwascalculated 3768mg/L,resultinginanacceptableRMSEof339mg/Lor foreachparameter.Itwasobservedthattheminorvariation less. Anderson and Woesner [46] and Moriasi et al. [47] intransmissivityorrechargeratevaluesaffectedthehydraulic reported that the RMSE is generally considered the best head impressively. The resulting plots of sensitivity showed calibrationindicator.Asgharetal.[48]saidthatthenegative the nonlinear response to recharge and transmissivity (Fig- valueofModelEfficiency(MEF)indicatesthehighvariability ures8and9). between the observed and simulated values. A zero value ofMEFindicatesapoorsimulation.Ifthemodelsimulated values exactly match the observed value then the MEF = 3.3. Predicted Groundwater Level. The contour lines of pre- 1. So, all the values showed a good fit between measured dictedgroundwaterlevelfortheyear2030underScenarioI JournalofChemistry 9 140 Table3:DispersivityvaluesusedinMT3D. 120 Possiblerangeand Parameter Usedvalue 100 source m) MS error ( 6800 Ldiosnpgeirtsuivdiitnya(lm) 5.06 31–.18–1951.–0254[5[[5443]2]] R 40 Horizontaltransverse 0.45 0.01–10[42,55] 20 dispersivity(m) Verticaltransverse 0 0.015 0.01–10[42,55] dispersivity(m) 0.4 0.6 0.8 1 1.2 1.4 1.6 Multiplying factor Figure8:Sensitivityanalysiswithrespecttorecharge. Table4:Valueofparametersfordispersivetransport. Layernumber TRPT(—) TRPV DMCOEF 35 1 0.09 0.003 0.1 30 2 0.09 0.003 0.1 m) 25 3 0.09 0.003 0.1 ror ( 20 TdiRspPeTrsiisvtihtye;raTtRioPoVftihsethhoerirzaotniotalotfratnhsevevresretidciaslpetrrasinvsivtyertsoethdeislpoenrgsiitvuidtyintaol MS er 15 tchoeeffilocnigenittu[dLin2Tal−d1]is.persivity;DMCOEFistheeffectivemoleculardiffusion R 10 5 0 costofprivatetubewellinPakistanwasUS$530fortheareas 0.4 0.6 0.8 1 1.2 1.4 1.6 wherethewatertabledepthwaslessthan6manditwasUS Multiplying factor $3206fortheareaswithmorethan24mdepth. Figure9:Sensitivityanalysiswithrespecttotransmissivity. 3.4.TransportofSalts. Thevalueoflongitudinaldispersivity was calculated as 5.06m (Table 3). Ahmad [53] reported that the values of longitudinal dispersivity were between andScenarioIIarepresentedinFigures10and11,respectively. 1.89 and 5m for the aquifer of the Indus Basin of Pakistan. Under Scenario I, the maximum depletion in groundwater Gelhar et al. [42] reviewed various researches and reported level would be up to 18m near Bhawana, where irrigation that the value of longitudinal dispersivity was between 3 is mainly dependent upon groundwater. Qureshi et al. [7] and 15.24m, whereas Engesgaard et al. [54] reported the reported that the depletion of groundwater was more pro- range of longitudinal dispersivity from 1 to 10m. The value nouncedinuncommandareasofthePunjab,wheresurface- of horizontal and vertical transverse dispersivity was found watersupplieswereconstrainedandagriculturewasheavily to be 0.45 and 0.015m, respectively, based on soil type dependent on groundwater. In the upper part of the study andhydraulicconductivity[42].Shiehetal.[55]concluded area,comparativelylessdepletion(2montheaverage)inthe thatthetransverse(horizontalandvertical)dispersivitywas groundwaterlevelwasobserved.UnderScenarioII,decline between0.01and10m.Narayanetal.[56]developedSUTRA, ofgroundwaterlevelby3-4mwasobserved.Incomparison asolutetransportmodelfortheLowerBurdekinDelta,North withScenarioIthedifferenceisofabout1-2m,whichisnot Queensland,andfoundlongitudinaldispersivityof2.5mand too high. In the lower part near Jhang and TT Singh, there transversedispersivityof0.5m. wouldbeariseofgroundwaterlevelbyupto2m(Figure11). ThedispersivityvaluesusedintheMT3Dmodelaregiven Shifting of canal water supplies would have good effect for inTable4.ThevalueofTRPT(ratioofthehorizontaltrans- replenishinggroundwater. versedispersivitytothelongitudinaldispersivity)andTRPV Thedepletionofgroundwaterlevelhasdirectandindirect (ratiooftheverticaltransversedispersivitytothelongitudinal impact on pumping cost. Kori et al. [17] discussed that the dispersivity) was 0.09 and 0.003, respectively. The effective declineingroundwaterlevelwouldincreasetheabstraction molecular diffusion coefficient (DMCOEF) describes the cost. The construction cost of a deep electric tube well diffusive flux of a solute in water from an area of greater (>20m)wasreportedasUS$5000ascomparedtoUS$1000 concentrationtowardsanareawhereitislessconcentrated. forashallow(<6m)tubewell.Obrienetal.[50]reportedthat Themoleculardiffusioncoefficientisgenerallyverysmalland 3 cost of pumping 103m water was US $8.61 and US $18.78 negligiblecomparedtothemechanicaldispersion,soitwas for31mand91mlift,respectively.Basharat[51]reportedthat 0.1[17,39].ThevalueofPecletnumberwasfoundtobe494. cost of pumping per cubic meter of groundwater increased ForsuchhighervalueofPecletnumber,thesolutetransport about3.5timesbecausethedepthofwatertabledroppedfrom isdominatedbytheadvection,thatis,thetransportofsolute 6to21m.Qureshietal.[52]calculatedthattheinstallation duetobulkflowofwater[43,57]. 10 JournalofChemistry N 5 0 2 32∘00N 0 0 2 18 190 195 0 145 155 170 185 5 31∘300N 16 5 7 1 0 4 1 135 0 6 1 0 31∘00N 5 15 4 1 0 510 20 30 40 (km) 72∘300E 73∘00E 73∘300E 74∘00E Groundwater head Contour 2030 Figure10:Predictedgroundwaterlevelin2030(ScenarioI). N 5 0 2 5 9 1 2 0 0 32∘00N 185 1 5 9 8 0 1 1 8 5 5 31∘300N 17 1 8 0 0 0 7 6 1 1 5 6 1 5 5 1 31∘00N 0 0 510 20 30 40 5 1 (km) 72∘300E 73∘00E 73∘300E 74∘00E Groundwater head Contour 2030 Figure11:Predictedgroundwaterlevelin2030(ScenarioII).