The intensive abstraction of groundwater is causing a number of problems such as groundwater depletion and quality deterioration. To manage such problems, the data of 256 piezometers regarding groundwater levels and quality were acquired for the period of 2003 to 2012 in command area of Lower Chenab Canal (LCC), West Faisalabad, Pakistan. MODFLOW and MT3D models were calibrated for the period of 2003–2007 and validated for years 2008–2012 with respect to observed groundwater levels and quality data, respectively. After the successful calibration and validation, two pumping scenarios were developed up to year 2030: Scenario I (increase in pumping rate according to the historical trend) and Scenario II (adjusted canal water supplies and groundwater patterns). The predicted results of Scenario I revealed that, up to year 2030, the area under good quality groundwater reduced significantly from 50.35 to 28.95%, while marginal and hazardous groundwater quality area increased from 49.65 to 71.06%. Under Scenario II, the good quality groundwater area increased to 6.32% and 12.48% area possesses less hazardous quality of groundwater. It was concluded that the canal water supply should shift from good quality aquifer zone to poor quality aquifer zone for proficient management of groundwater at the study area.

1. Introduction

In Pakistan, groundwater is the second largest source of irrigated agriculture because of the arid climatic conditions [1, 2]. The reliance on groundwater has significantly increased during the past two decades to meet the food and fiber requirements of growing population [3]. Furthermore, the surface irrigation system has high degree of conveyance and application water losses. Thus, the system operates at an efficiency of less than 40% that meets only less than 40% of the crop water requirements [4]. The provision of irrigation water from other sources such as groundwater is indispensable for potential productivity [5]. The growth rate of private tube wells has increased significantly at the rate of 60% in Punjab province from 1991 to 2000. About 1.20 million tube wells of capacity ranging from 0.015 to 0.056 m3/sec have been installed in Pakistan for irrigation at depth of 30–85 m having diameter of 15–30 cm, out of which 86% tube wells are installed in Punjab [6]. The farmers of tail ends of the distributaries and watercourses are forced to rely heavily on groundwater, particularly where canal water supplies are constrained. Without groundwater availability, not only Punjab but the whole country would face a severe water shortage that leads to food shortages as Punjab produces more than 90% of total grains [7]. However, the unchecked abstraction of groundwater has created serious negative concerns in terms of lowering water table and saltwater intrusion that may lead to the issues of sustainability of usable groundwater resources. As a result of saline groundwater intrusion, about 200 public tube wells initially installed in the fresh groundwater zone of Punjab and Sindh provinces had been abandoned [8]. Shah [9] reported similar problems existed in most of the irrigated regions of the world; those were further increasing rapidly and negatively affecting agricultural productivity.

There are wide spatial and temporal variations in groundwater quality in the Indus basin, which is due to the pattern of groundwater movement in the aquifer [10]. The belt of fresh groundwater is generally available near the main rivers and canal on account of high recharge of fresh seepage water. But, then, the quality of groundwater changes to unfit as laterally away from the rivers [5, 11]. The continuous and intensive use of groundwater for irrigation is adding plenty of salts causing secondary salinization because groundwater generally has more salts than canal water. Bakhsh and Awan [12] reported that application of groundwater having total dissolved solids (TDS) of 1000 mgL−1 up to a soil depth 370 mm changed the top 300 mm depth of nonsaline into a saline soil that impaired crop productivity. In many agricultural areas of Pakistan, the usage of poor quality tube well water for irrigation is considered as one of the major causes of salinity and, consequently, lower food productivity [13]. According to an estimate in Pakistan, the secondary salinization degraded the crop land which reduced the production potential of major crops by 25%, valued at an estimated loss of US $250 million/year [14]. Low quality water and soil salinity can affect plant growth and soil structure in several ways, directly and indirectly [15].

Thus, the assessment of groundwater quality is important to ensure the safe use of these resources on a sustainable basis. Groundwater models are most widely used tools for efficient management of precious groundwater resources and to predict different future scenarios [16]. Different groundwater modeling codes are available, each with their own capabilities, operational characteristics, and limitations such as PMWIN, FEFLOW, SVFlux, and GWVistas. But the most extensively used three-dimensional groundwater flow model among the available models is PMWIN (Processing MODFLOW for Window) [17, 18]. Its popularity has continued, in part due to the modularity of the program, resulting ability, and user friendly interface [19]. It uses a block-centered finite-difference scheme for saturated zone. The advantages of PMWIN include numerous facilities for data preparation, easy exchange of data in standard form, extended worldwide experience, continuous development, availability of source code, and relatively low price or being freely available [20]. A number of research studies have been conducted regarding application of modeling approach for groundwater management in different part of the world. Moeck et al. [21] developed a 3D groundwater model for simulating existing and proposed water management strategies as a tool to ensure the utmost security for drinking water in Basel, Switzerland. Gebreyohannes et al. [22] developed regional groundwater flow model for Geba Basin, northern Ethiopia, and reported that none of the hydrogeological formations can be exploited for large-scale groundwater exploitation. Rahmawati et al. [23] conducted research to study saltwater intrusion from 1995 to 2108 in Semarang city based on well log data and MODFLOW numerical model was used. For salt intrusion projection in the future, the sea level rise projection also was conceived. Kori et al. [17] linked groundwater flow model (MODFLOW) with solute transport model (MT3D), and several simulation runs were carried out after successful calibration of model for two sampled sites located at JRS-57 and JRS-60 tube wells at Nawabshah-Pakistan. Carretero et al. [24] applied groundwater modeling technique to study the impact of a possible rise of 1 m in sea level against the low-lying coast of Partido de La Costa, Argentina. Similarly, Abu-el-Shar and Hatamleh [25] simulated the groundwater model (PMWIM) for the Azraq Basin, Jordan, to manage groundwater. The large number of latest research papers—Asoka et al. [26], Galitskaya et al. [27], Kambale et al. [28], Abdullah and Morteza [29], and Durand et al. [30]—dealing with groundwater management using modeling techniques were published is a testimony to the important role played by the models.

Therefore, the assessment of groundwater quality through modeling is important to ensure the safe use of these resources on a sustainable basis. If unchecked pumping of groundwater continues to remain in the study area, the irrigation tube wells would not be able to lift water at their present quality. Hence, there is a dire need to investigate the impact of groundwater flow conditions and overexploitation on groundwater quality. Thus, the followings were the two main objectives of the current study: to develop a regional hydrogeological groundwater flow model and observe its future trend and to formulate groundwater management strategy for its proficient utilization.

2. Material and Methods

2.1. Study Area

The research study was carried out in the command area of Lower Chenab Canal (LCC), West Faisalabad, Pakistan, having longitude of 73.85° to 72.18° (E) and latitude of 32.32° to 30.85° (N) (Figure 1). The site has gross command area of 1.16 Mha and culturable command area of 0.98 Mha, respectively. It falls in rice-wheat agroecological zone of province Punjab, Pakistan. It is comprised of vast canal network such as main canals, branch canals, and minor distributaries. The River Chenab and Gugera Branch Canals are located at the northwest and southeast side, respectively. Similarly, Qadirabad-Baluki Link Canal is located in the northeast side and Trimmu-Sidhnai Link Canal is present in the southwest side of the study area (Figure 2).

2.2. Climatic

The summer season starts from April and continues till October. The temperature varies in the range of 21 to 51°C during the summer season. Similarly, the winter season starts from October and lasts till April with temperature ranges of 1 to 27°C. The average annual precipitation was estimated to be 439 mm. The ten-year average value of was 1413 mm/year.

2.3. Hydrogeological Conditions

Lithology of aquifer system in the study area has different classification according to different textural characteristics. The surface soil textures are largely fine and moderately medium with good permeability properties. The aquifers of study area were formed as a result of sediment deposition due to flat topography. These sediments were carried by the river waters from the vast alluvial basin of rivers as a result of materials washed down from Himalayan Mountain. In general, like other aquifers in Punjab, the aquifers of study area are unconfined. The groundwater levels in the year 2012 were present in the range of 145 to 205 m (Figure 3). The groundwater level was higher in the upper part and gradually decreased towards the lower side, which means that the groundwater flow was from upper to lower side.

The good quality groundwater is present in the northern part whereas poor quality groundwater is in the southern part of the study area. The average spatial variation of measured groundwater TDS in the entire study area for the year 2012 is given in Figure 4. A calibrated TDS meter (TDS-98302) was used to measure the TDS of groundwater with an accuracy of ±2%. The calibration of TDS meter was accomplished in the standard solution in the Water Quality and Environment Laboratory, Department of Irrigation and Drainage, University of Agriculture, Faisalabad. The measured TDS values indicated that almost half of the aquifer of the study area has marginal to hazardous quality of groundwater. A 30 km wide strip along the Chenab River had groundwater of good quality, whereas hazardous quality groundwater existed in 20 km wide strip at a distance of 40 km away from the Chenab River. The groundwater quality zones of fresh and saline water were formulated according to criteria developed by WAPDA [31] for irrigation water quality (good: TDS < 1000 mg/L; marginal: TDS of 1000–1700 mg/L; and hazardous: TDS > 1700 mg/L), Shakoor et al. [6]. The northern and northwestern part of the study area (along Chenab River) contained good quality groundwater. The marginal and hazardous quality was laid in the southern part of the study area.

The annual water balance of the whole study area from the year 2003 to the year 2012 is shown in Figure 5. The water balance describes the volume of water entering, subtraction, and net storage in the aquifer system. The water entering parameters are recharge and river leakage, while subtraction parameters are wells and ET. The positive storage term showed that volume of water extracted from the aquifer was more than the recharge which means crop water demand was fulfilled from aquifer storage and vice versa. The storage in the aquifer in 2003 and 2012 was −444.5 and 326.6 MCM, respectively, clearly indicating the replenishment of aquifer in 2012. Similarly, the volume of water withdrawn through tube wells increased from 3298 to 3759 MCM from years 2003 to 2012, respectively, also showing increased demand of groundwater.

2.4. Model Selection

PMWIN 5.3 (Processing MODFLOW for Window) was used in this research. It is a simulation system based on the modular three-dimensional finite-difference technique for modeling groundwater flow and contamination in groundwater with a wide range of natural systems. PMWIN is used widely throughout the world and it was applied in many groundwater modeling applications [17, 23, 25].

2.5. Model Input and Calibration

The MODFLOW and MT3D models are packages of PMWIN 5.3. They were calibrated with respect to observed groundwater level and quality data, respectively, using inverse modeling method. The data of 256 piezometers regarding groundwater level and quality were acquired from the Department of Land and Reclamation, Faisalabad, for the period from 2003 to 2012 to achieve the model calibration and validation (Figure 1). The piezometer wells have maximum depths of 55 m. The depth of piezometers varied at some location depending upon the water table conditions. The piezometer wells have diameters of 5 cm and strainer lengths of 3 m.

The model area had a square geometry and the whole area was divided into 76 columns and 76 rows. The total number of cells was of 5776 cells. Each square cell has a dimension of 2.5 km × 2.5 km. The model area had 3653 inactive cells, which were outside the boundary of the study area, while the model area had 2123 active cells located within the boundary of study area. The cell size of 2.5 km × 2.5 km (6.25 km2) was also used for groundwater modeling studies by Al-Fatlawi [32] in Umm Er Radhuma, the Western Desert, Iraq, and by Khan et al. [33] in Rechna Doab. Abu-el-Shar and Hatamleh [25] developed groundwater model for the Azraq Basin and the biggest cell size of 8.69 km2 was selected. Similarly, Schoups et al. [34] used cell size of 2 km × 2 km to calibrate groundwater model of the Yaqui Valley, having 6800 km2 irrigated agricultural region located along the Sea of Cortez in Sonora, Mexico.

Lithology of aquifer system in the study area was obtained from the Water and Power Development Authority (WAPDA) [31]. The soils have different classification according to different textural characteristics. The surface soil textures are largely fine and moderately medium, with good permeability properties. The areas of 4451.3 (38.48%), 4987.3 (43.08%), 1621 (14%), 464.1 (4%), and 52.3 km2 (0.45%) have soil texture of fine, moderately medium, medium, moderately course, and course, respectively. The aquifer of the study area was defined with four different layers depending upon their lithological data [31]. The spatial domain represented in the model consisted of four layers (0–7, 7–30, 30–90, and 90 m to bedrock). The horizontal and vertical hydraulic conductivities have large variation from one to the other side of the study area (Table 1). The minimum and maximum values of 1–265 m/day and 1–15 m/day, respectively, for horizontal and vertical hydraulic conductivities were used in the model, as majority of the tube wells were installed in second layer because 80% part of study area has in the range of 70–100 m/day. The similar range of values for hydraulic conductivity within Punjab province domain was used by Jehangir et al. [35]; Ahmad [36]; Arshad [37]; and Khan et al. [33]. The specific storage values for layer 1 ranged as 0.0001–0.001 m−1 and for the remaining layers the values of 0.00001–0.0003 m−1 were used. The values of specific yield for layer 1 ranged as 00.05–0.25 and for the remaining layers the values of 0.05–0.20 were used. The effective porosity of 0.25 was given to all layers [31]. The simulation time unit “days” and simulation flow type “transient” was selected. Two stress periods in each year were considered to represent the Kharif and Rabi seasons having 183 and 182 days, respectively, with six time steps in each stress period, as there are two main cropping seasons based on agroclimatic conditions in Pakistan, Kharif and Rabi. Kharif starts from June and July and goes to October and November, while the Rabi season starts from September and October and continues to April and May.

The cumulative evapotranspiration during a period (Kharif or Rabi) was divided by its duration in days and thus the evapotranspiration rate per day was calculated using CROPWAT. The evapotranspiration rates of 0.006 and 0.003 m/day were used for odd (Kharif) and even (Rabi) stress periods, respectively, throughout the model stress periods for model calibration. The recharge package was used to simulate spatial distribution of recharge water at field from rainfall and irrigation to the groundwater system. The minimum and maximum recharge values from 0.00065 to 0.0013 m/day were used for odd (Kharif) stress periods, whereas values from 0.00026 to 0.0005 m/day were used for even (Rabi) stress periods. The river package was used to simulate the flow between an aquifer and surface-water features such as rivers, canals, lakes, and reservoir. The hydraulic features, such as canal length, bed width, full-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 groundwater demand was calculated by subtracting the net canal water supplies from net crop water requirement. After providing all the input data, the calibration process was started. In this calibration process, input parameters such as recharge, pumping rate, and hydraulic conductivities were adjusted in MODFLOW.

2.6. Solute Transport Model

The MT3D numerical model with PMWIN 5.3 user interface was used to simulate advection and dispersion contaminants in the three-dimensional groundwater flow system of the study area. The MT3D simulated groundwater flow by numerically solving the groundwater flow and solute transport equations. The partial differential equation describing the three-dimensional transport of dissolved solutes in the groundwater [38] can be written in where is the concentration of contaminants dissolved in groundwater [ML−3], is the distance along the respective Cartesian coordinate axis [L], Di is the hydrodynamic dispersion coefficient [L2T−1], is the seepage or linear pore water velocity [LT−1], is the volumetric flux of water per unit volume of aquifer representing sources (positive) and sinks (negative) [T−1], is the concentration of sources or sinks [ML−3], is the chemical reaction term [ML−3T−1], and is time [].

The calibration of MT3D was achieved by adjusting the values of advection and dispersion values. Advection depicts mass transport basically due to the bulk flow of water in which the mass is dissolved or movement of solute as a result of groundwater flow. Advection is the dominant process of solute transport due to groundwater flow. For a given time step, the movement of solutes due to advection was simulated using the corresponding field velocity which was computed by the calibrated MODFLOW model and the effective porosity of the aquifer. MXPART of maximum 1300000 particles was allowed in a simulation while the value PERCEL (Courant number or number of cells in any particle) was selected as 0.79. It was allowed to move in any direction within one transport step and generally ranged from 0.5 to 1 [39]. The default values of the remaining parameters were used for MT3D model.

The longitudinal dispersivity and transverse dispersivity are required to estimate dispersive transport of the solute. The dispersivity is a characteristic property of the porous medium that describes the spreading of the solute in a medium. Longitudinal dispersivity is used when the spreading of solute is in the direction of bulk flow, whereas the transverse dispersivity is the perpendicular (vertical and horizontal) spreading of solute to the direction of bulk flow. The longitudinal dispersivity was computed from molecular diffusion and local dispersion coefficient based on heterogeneity of the medium [40, 41] as given in where is longitudinal dispersivity (cm), is coefficient of molecular diffusion (16 × 10−6 cm2/s), is velocity of flow (6.8 × 10−6 cm/s), is dispersion coefficient (30 × 10−6 cm2/s), and is asymptotic longitudinal dispersivity (500 cm). Asymptotic longitudinal dispersivity is due to heterogeneity of the medium, mainly dependent on the variance of the log transformed conductivity and correlation length in the mean direction of flow. As no field data were available for solute transport in aquifer media, the value of was adopted as 500 cm after [42]. The Peclet number, , was calculated to measure the relative significance of advection and dispersion in the study area, where is the size of a cell equal to 2500 m. The unitless Peclet number is usually used to decide the dominant factor in solute transport among advection and dispersion [43]. After the successful calibration of both models (MODFLOW and MD3D) for years 2003–2007, the models validated for the years 2008–2012 and two future scenarios were simulated.

2.7. Future Scenarios

The model prediction was accomplished in order to investigate the response of the model for two future scenarios regarding the pumping rate up to year 2030. It was assumed that there will be no uncertain change or tragedy in climate and irrigation system. In Scenario I, the pumping will increase according to the historical trend. As in Punjab province of Pakistan, about 1.20 million tube wells have been installed and increasing at 5.5% annually [44]. This indicated that the amount of withdrawal from aquifer may increase from 3738 MCM in year 2012 to 6008 MCM in year 2030 whereas recharge would be about 2664 MCM according to historic trend. In Scenario II, one of the water management options for the study area was proposed. In the upper part of the study area (Pindi Bhatiyan-Safdrabad) where groundwater has good quality, the rate of groundwater abstraction was increased and recharge through an irrigation system was decreased by 35%, while in the lower part, where groundwater has poor quality, the groundwater abstraction was decreased and recharge through an irrigation system was increased by the same rate.

3. Results and Discussion

3.1. Calibration and Validation of Models

The degree of fit between model simulations and field measurements was quantified by statistical means (Table 2) and all parameters were found in acceptable range. The minus sign of mean error (ME) represented that the model simulated values were higher than the measured head. The calibration criterion for hydraulic head is root mean squire error (RMSE) which is less than or equal to 10% of head variation within the aquifer being modeled [45] and the same criterion is also followed in groundwater TDS. The head in the aquifer within the study area varies from approximately 145 to 205 m, resulting in an acceptable RMSE of 6 m or less. Similarly, the TDS in the aquifer at selected points are varied from 384 to 3768 mg/L, resulting in an acceptable RMSE of 339 mg/L or less. Anderson and Woesner [46] and Moriasi et al. [47] reported that the RMSE is generally considered the best calibration indicator. Asghar et al. [48] said that the negative value of Model Efficiency (MEF) indicates the high variability between the observed and simulated values. A zero value of MEF indicates a poor simulation. If the model simulated values exactly match the observed value then the MEF = 1. So, all the values showed a good fit between measured and simulated values. The average coefficient of determination for the selected piezometers for calibration and validation period was calculated as 0.89 and 0.87 (Figures 6 and 7) for MODFLOW and MT3D models, respectively, which indicated a close agreement between the calculated and measured groundwater head. Hagos [49] calibrated PMWIN model for Raya valley, Ethiopia, with respect to groundwater level. The values of ME, MAE, and RMSE of calibrated results were −1.4 m, 7.8 m, and 10.7 m, respectively, with coefficient of determination of 0.97 and reported to be satisfactory.

3.2. Sensitivity Analysis

The sensitivity analysis showed that recharge and transmissivity of the aquifer were most sensitive parameters. The factors of 0.5, 0.8, 0.9, 1.1, 1.2, 1.3, and 1.5 were multiplied with the calibrated values of recharge and transmissivity. The resulting hydraulic heads were then compared with the observed heads and RMSE was calculated for each parameter. It was observed that the minor variation in transmissivity or recharge rate values affected the hydraulic head impressively. The resulting plots of sensitivity showed the nonlinear response to recharge and transmissivity (Figures 8 and 9).

3.3. Predicted Groundwater Level

The contour lines of predicted groundwater level for the year 2030 under Scenario I and Scenario II are presented in Figures 10 and 11, respectively. Under Scenario I, the maximum depletion in groundwater level would be up to 18 m near Bhawana, where irrigation is mainly dependent upon groundwater. Qureshi et al. [7] reported that the depletion of groundwater was more pronounced in uncommand areas of the Punjab, where surface-water supplies were constrained and agriculture was heavily dependent on groundwater. In the upper part of the study area, comparatively less depletion (2 m on the average) in the groundwater level was observed. Under Scenario II, decline of groundwater level by 3-4 m was observed. In comparison with Scenario I the difference is of about 1-2 m, which is not too high. In the lower part near Jhang and TT Singh, there would be a rise of groundwater level by up to 2 m (Figure 11). Shifting of canal water supplies would have good effect for replenishing groundwater.

The depletion of groundwater level has direct and indirect impact on pumping cost. Kori et al. [17] discussed that the decline in groundwater level would increase the abstraction cost. The construction cost of a deep electric tube well (>20 m) was reported as US $5000 as compared to US $1000 for a shallow (<6 m) tube well. Obrien et al. [50] reported that cost of pumping 103 m3 water was US $8.61 and US $18.78 for 31 m and 91 m lift, respectively. Basharat [51] reported that cost of pumping per cubic meter of groundwater increased about 3.5 times because the depth of water table dropped from 6 to 21 m. Qureshi et al. [52] calculated that the installation cost of private tube well in Pakistan was US $530 for the areas where the water table depth was less than 6 m and it was US $3206 for the areas with more than 24 m depth.

3.4. Transport of Salts

The value of longitudinal dispersivity was calculated as 5.06 m (Table 3). Ahmad [53] reported that the values of longitudinal dispersivity were between 1.89 and 5 m for the aquifer of the Indus Basin of Pakistan. Gelhar et al. [42] reviewed various researches and reported that the value of longitudinal dispersivity was between 3 and 15.24 m, whereas Engesgaard et al. [54] reported the range of longitudinal dispersivity from 1 to 10 m. The value of horizontal and vertical transverse dispersivity was found to be 0.45 and 0.015 m, respectively, based on soil type and hydraulic conductivity [42]. Shieh et al. [55] concluded that the transverse (horizontal and vertical) dispersivity was between 0.01 and 10 m. Narayan et al. [56] developed SUTRA, a solute transport model for the Lower Burdekin Delta, North Queensland, and found longitudinal dispersivity of 2.5 m and transverse dispersivity of 0.5 m.

The dispersivity values used in the MT3D model are given in Table 4. The value of TRPT (ratio of the horizontal transverse dispersivity to the longitudinal dispersivity) and TRPV (ratio of the vertical transverse dispersivity to the longitudinal dispersivity) was 0.09 and 0.003, respectively. The effective molecular diffusion coefficient (DMCOEF) describes the diffusive flux of a solute in water from an area of greater concentration towards an area where it is less concentrated. The molecular diffusion coefficient is generally very small and negligible compared to the mechanical dispersion, so it was 0.1 [17, 39]. The value of Peclet number was found to be 494. For such higher value of Peclet number, the solute transport is dominated by the advection, that is, the transport of solute due to bulk flow of water [43, 57].

3.5. Predicted Groundwater Quality

Predicted variations in groundwater TDS under Scenario I up to year 2030 are presented in Figure 12. For the year 2030, the 28.95% area will have good quality groundwater, 38.35% area will have marginal quality, and the remaining 32.71% area will have hazardous quality groundwater of the aquifer with TDS concentration of greater than 1000 mg/L. The results showed that 21.40% good quality area will be converted from marginal into poor quality after 18 years. So, marginal and hazardous quality groundwater areas are expected to increase by 19.28 and 2.13%, respectively. The area of hazardous quality water almost remained the same but its severance increased due to less recharge of fresh water than discharge. The change in groundwater TDS concentration under Scenario I up to year 2030 is shown in Figure 13.

The negative sign indicated the increase in groundwater TDS concentration and vice versa. The groundwater TDS of about 85% of the aquifer would increase up to 1500 mg/L and the remaining area up to 2500 mg/L located at southern part of the study area near Jhang and TT Singh. The groundwater quality of <1% area would increase up to 140 mg/L near Hafizabad because of having relatively high potential of recharge irrigation network and rainfall. The predicted results revealed that groundwater quality of the northern part of the study area will not be affected with the increase in the pumping due to balance between discharge and recharge. The irrigation networks in that upstream area are well established and generally received more canal water than the allocated share [50, 51].

Simultaneously, the groundwater quality in middle and lower part will be affected severely by salinity due to imbalance between discharge and recharge of fresh water. The farmers of downstream areas received about 21% less canal water than designed discharge imposed farmers to use groundwater as a primary source [51, 52]. This abundant abstraction of groundwater is deteriorating the groundwater quality due to both horizontal and vertical saline intrusion [58, 59]. Prinos [60] described the multiple pathways for saltwater intrusion into freshwater zone in Florida near the well field. Khan et al. [33] said that overpumping of groundwater from aquifer induced the lateral saltwater intrusion into the fresh groundwater area and the vertical upconing of the saline interface is resulting in degradation of aquifers. The farmers in the study area abstracted groundwater without any quality check. The continuous use of such poor quality water is causing secondary salinization and ultimately reduces crop productivity [6, 51]. The groundwater quality results of Scenario I clearly indicated that an increase in the number of tube wells in the future could cause the problem of salinization; therefore, the groundwater regulation aimed at protecting the quality and quantity of groundwater resource must be implemented [6163].

The predicted results of the spatial variation of groundwater salinity up to year 2030, under Scenario II, are shown in Figure 14. It shows that 35.27% area of the aquifer will have good quality, 44.48% area will have marginal quality, and the remaining 20.23% area will have hazardous quality groundwater. The analysis of the groundwater quality results indicated that if the groundwater abstraction is increased by 35% and irrigation recharge is decreased by 35% in the upper part of the study, the overall increase in groundwater salinity will remain under the safe limit, while, in the lower part of the study area, if abstraction is decreased by 35% and irrigation recharge is increased by 35%, there will be an improvement in groundwater quality up to 500 mg/L (Figure 15). Under Scenario II, the area under good and hazardous quality increased by 15.08 and 10.35%, respectively, compared with current status.

The results of groundwater quality for Scenarios I and II are summarized in Table 5. The analysis of the model results of both scenarios up to year 2030 showed that 6.32% area has good quality water and 12.48% has less hazardous quality in Scenario II compared with Scenario I. In Scenario II, 6.13% more area is under marginal quality, possibly due to shifting of water quality from hazardous quality due to higher recharge from increased canal water supply in the lower part of study area. Hence the results of Scenario II revealed that the need of the hour is a need for a shift in surface irrigation water to the lower part of the study area to reduce groundwater salinity problems. This will also provide formers with a way forward for conjunctive surface and groundwater irrigation technology for sustainable water management for this region and to avoid secondary salinization [59]. Foster and van Steenbergen [64] concluded that the evolution to more planned conjunctive use of groundwater and surface-water resources offers great potential for increasing water-supply security in both irrigated agriculture and urban water supply across the developing world, especially on large alluvial plains which are often major centers of population and economic development. It was also reported that there are spontaneous conjunctive groundwater and surface-water use in Indian, Pakistani, Moroccan, and Argentinean irrigation-canal commands which have largely arisen due to inadequate surface-water supply to meet irrigation demand.

4. Conclusions

The MODFLOW and MT3D models were calibrated for the period of 2003–2007 and validated for the years 2008–2012 using the measured groundwater level and quality data, respectively, for regional groundwater management in Punjab province of Pakistan. The main conclusions drawn from this study were as follows:(i)The surface soil textures are largely fine and moderately medium with good permeability properties. In general, the aquifers of the study area are unconfined.(ii)Almost half of the aquifer of the study area has marginal to hazardous quality of groundwater. A 30 km wide strip along the Chenab River had groundwater of good quality. It was also observed that, at a distance of 40 km away from the river, a 20 km wide strip of hazardous quality was found.(iii)The volume of water withdrawn through tube wells increased by about 14% from the years 2003 to 2012 and caused the rapid lowering of water table.(iv)The predicted results of Scenario I revealed that the area good quality groundwater will reduce to 21.4%, while marginal and areas hazardous quality water will increase as 19.28 and 2%, respectively.(v)The most affected areas will be found in the middle and lower parts of the study area such as Faisalabad, Jhang, Gujra, and TT Singh.(vi)Under Scenario II, area of good quality groundwater will increase to 6.32% area and 12.48% area will possess less hazardous quality of groundwater.

These results indicated that the canal water supply should shift from good to poor quality groundwater region for proficient management of aquifer water. This can be achieved by regulating the surface-water supplies from irrigation headworks through irrigation department. It is also recommended that the detailed field survey is required to determine the exact advective and dispersive solute transport parameters for more accurate results of MT3D model.

Conflicts of Interest

The authors declared that there are no conflicts of interest regarding the publication of this paper.


The authors acknowledge the financial support from Higher Education Commission (HEC), Pakistan, under Ph.D. Indigenous Scholarship Program. They are also grateful to WAPDA, Pakistan Meteorology Department, and Department of Land and Reclamation, Faisalabad, for providing data to accomplish this research.