Regional Groundwater Quality Management through Hydrogeological Modeling in LCC, West ...

17 Pages·2017·2.19 MB·English

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