3-D numerical modelling of groundwater flow for scenario-based analysis and management

In recent years, extensive competition for groundwater use among different consumers has exploited major freshwater aquifers in Pakistan. There is an urgent need for appraisal of this precious resource followed by some mitigation strategies. This modelling study was conducted in the mixed cropping zone of the Punjab, Pakistan. Both remote sensing and secondary data were utilized to achieve objectives of this study. The data related to piezometric water levels, canal gauges, well logs, meteorological and lithological information were collected from Punjab Irrigation Department (PID), Water and Power Development Authority (WAPDA). Groundwater flow models for both steady and transient conditions were set-up using FEFLOW-3D. Water balance components and recharge were estimated using empirical relations and inverse modelling approaches. Both manual and automated approaches were utilized to calibrate the models. Moreover, sensitivity analysis was performed to see the response of model output against different model input parameters. Followed by calibration and validation, the model was run for different management scenarios, including lining of canal sections, minimization of field percolation, and change of groundwater abstraction. The study results show a drop in groundwater levels for almost all scenarios. The highest negative change was observed for the 4 scenario (i.e. 25% increase in groundwater pumping over a 10-year period), with a value of 3.73 m, by ignoring very wet summer and winter seasons. For normal weather conditions, the highest negative change was observed for the 4 scenario with a value of 2.91 m followed by 2.68 m for the 3 scenario (i.e. 50% reduction in canal seepage and field percolation over a 10-year period). For very wet summer and winter seasons, only one positive change was observed, for the 5 scenario (i.e. 25% decrease in groundwater pumping during 10 years period), with a value of 1.17 m. The changes for all other scenarios were negative. The mitigation strategy may include less groundwater pumping, by supporting cultivation of low delta crops and adjusting cropping patterns considering canal water supplies. It is further suggested to support current modelling results by incorporating more detailed information on cropping and by exploring the effect of climate change.


INTRODUCTION
About 70% of global groundwater resources are consumed by agriculture (Döll, 2009).The share of groundwater used is about 750-800 km3, which is about one-sixth of the total freshwater use (Shah, 2000).Although the share of groundwater used for agriculture is lower than that of surface water globally, the advantages of groundwater use, like on-demand availability, accessibility, reliability, and high productivity, outweighs those of surface water use (Qureshi et al., 2009).
Groundwater use is inevitable for agriculture in Pakistan.About 40 to 50% of irrigation water needs in the country are fulfilled by exploiting groundwater (Usman et al., 2015;Kazmi et al., 2012).The share of groundwater used for meeting domestic and industrial requirements is about 10%, and this value is growing due to urbanization, industrialization and improved living standards.Currently, in the Punjab Province, about 90% of daily domestic demands are met by groundwater.The consequences of this use include lowering of groundwater levels in many parts of the country, along with secondary secondary salinization and pumping costs (Qureshi et al., 2009).According to Qureshi (2010), there are about 21 canals commands (irrigated areas receiving water from a particular canal) in Punjab and Sindh Provinces where the groundwater table depth has fallen by up to 300 cm over a period of 10 years (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002).
The current study was conducted in the lower parts of the Lower Chenab Canal (LCC), Punjab, Pakistan.The region is considered important for the production of agricultural commodities but the sustainability of groundwater use is challenging due to water quantity and quality issues.Analysis of electrical conductivity (EC) shows that the majority of the area, including Bhagat and Sultanpur irrigation subdivisions, has saline groundwater in deeper aquifer depths, with a high concentration of salts.According to locally conducted groundwater sampling results, about 64.5% and 59.4% of the area of Bhagat irrigation subdivision encompasses unfit groundwater during the post-and pre-monsoon seasons.Similarly, about 90% of the total area of Sultanpur irrigation subdivision has groundwater quality that is marginal during both pre-and post-monsoon seasons.Due to this salinity problem, there is always competition for surface water use and groundwater abstraction is considered mandatory, which results in up-coning of saline water and consequently decreased agricultural productivity.
Considering the importance of groundwater use for agriculture, its proper distribution needs to be ensured, for meeting the future demands of food and fibre, which will increase in coming decades due to increases in population.This is only possible through proper understanding of the problem in order to devise workable strategies.Quantitative solution of the problem using 3-D numerical modelling is a potential means to conceptualize the groundwater resources and to make viable plans for its sustainability (Al-Zyoud et al., 2015).This particular study demonstrates the use of numerical modelling 147 for describing the problem and exploring the effects of different management scenarios on groundwater fluctuations.

STUDY AREA Location, climate and agriculture
Sultanpur irrigation subdivision was selected as a case study region in the LCC (Fig. 1).LCC off-takes from Khanki headwork on the river Chenab.The major crops cultivated include sugarcane, fodder and cotton in kharif and wheat during rabi cropping season.The two cropping seasons, rabi and kharif, usually prevail from November to April and from May to October, respectively.

Geological settings
The study area is a part of an abandoned floodplain, which is underlain by highly stratified unconsolidated alluvial material composed of sand of various grades interbedded with discontinuous lenses of silt, clay and nodules of kanker -a calcium carbonate of secondary origin deposited by present and ancestral tributaries of the Indus River.The bottom of the aquifer is formed by metamorphic and igneous rocks of Precambrian age (Sarwar and Eggers, 2006).The horizontal conductivity is generally an order greater than vertical (Usman et al., 2015).Khan (1978) has synthesized the results of pumping tests, lithological and mechanical analyses of the test boreholes, and reports that hydraulic conductivity varies from 24 to 264 m•day-1 and specific yield from 1 to 33% in different parts of LCC.

Hydrogeology/conceptual model
The conceptual model of the study area is shown in Fig. 2. The model area is bounded by a Trimu-Sidhnai link canal and the River Ravi on its western and southern sides.The northern and part of the eastern model boundaries are not connected with any water body and so interact with the modelling region through lateral groundwater inflow and outflow.The sections were identified using contour maps and lateral inflows/ outflows were estimated by Darcy's Law (Usman et al., 2015).Overall, the groundwater flow is parallel to this model boundary in a direction from north-east to south-west.The region is equipped with a wide irrigation network consisting of main canals, branch canals, major and minor distributaries.The irrigation water is supplied through the lower Chenab canal irrigation network.The transport of water through this irrigation network contributes a major recharge to groundwater at different stages of the irrigation system.The elevation of the region drops smoothly from north-east to south-west and causes regional groundwater flow along this elevation drop.As land use in the area is dominated by agriculture, groundwater pumping occurs along with canal water use.The major proportion of recharge takes place from monsoon rainfall and also from field percolation losses.The climate is regarded as semi-arid with high air temperatures which cause huge evapotranspiration losses.

Data acquisition and processing
Table 1 summarizes different data types utilized for the current study and its acquisition from different sources.

Setting of numerical model and scenario-based simulations
The triangular mesh algorithm (Shewchuk, 2005) was selected for finite element meshing owing to its fast speed and its capability to accommodate polygons, lines and points.The mesh quality was tested by Delaunay criteria for ensuring maximum model stability.Projection of 3-D unsaturated aquifer with saturation case was selected as the problem class.The aquifer was subdivided into 4 layers and 5 slices with the upper-most layer set as phreatic (i.e.unconfined).The thickness of different modelling layers was decided on the basis of the hydrogeological model and strata information from bore log data (Figure A1, Appendix).The total number of nodes per slice and elements per layer were 15 and 243, and 30 and 765, respectively.Both steady and transient models were prepared.The head results of the steady model were utilized as initial conditions for the transient model.Regionalization of both hydraulic conductivity and specific yield were performed using Akima interpolation algorithm (Akima, 1970) with linear interpolation type, with 5neighbours and zero over/ undershooting criteria.'Constant' (i.e.no time step) was used for the steady model, whereas for the transient model, the Adams-Bashforth time stepping scheme was used.The Adams-Bashforth scheme is an automatic time-step control scheme which gives good convergence and stability to model runs.PARDISO -Parallel Direct Solver (Schenk and Gärtner, 2004) was chosen for solving the equation system.

Model boundary conditions
In FEFLOW, recharge is considered as a material property instead of a boundary condition, which is estimated by using a water balance approach and modelling techniques, as can be seen in Usman et al. (2015).The recharge is assigned as a time-variant material property to the transient model.Prescribed head boundary conditions were assigned along the Trimu Sidhnai link canal and river Ravi based on their historical water-level records (Sarwar and Eggers, 2006), and by working out the water levels in the course of the channel with DEM.The Dirichlet type, time-varying boundary conditions were applied instead of Cauchy-type boundary conditions.The Cauchy boundary condition requires channel bed conductance and its geometry data, which were inaccessible from any source.With the eastern and northern model sides, Neumann-type boundary conditions were employed, for which the inflow and outflow sections and flow to the model domain were identified using piezometric water levels and by application of Darcy's law.The detailed methodology can be seen in Usman et al. (2015).Nodal source/sink boundary condition was applied for  groundwater pumping.The rates of pumping were estimated separately for all seasons employing a procedure described in Usman et al. (2015).

Model calibration, validation and sensitivity analysis
Figure 3 presents the flow diagram showing how calibration, validation and sensitivity analyses were conducted in the current study.The calibration of models was performed using both a manual approach and automated approach using PEST (Doherty, 2002).The pilot point technique (De Marsily et al., 1984) was used to define spatially variable parameters using kriging and each slice/layer of the model was treated as a single zone for performing the interpolation.The advanced Tikhonov regularization was implemented in PEST to encounter non-uniqueness of parameters.Tikhonov regularization is implemented as the number of estimation parameters was quite large (246 pilot points in total).The sensitivity analyses were performed manually for both steady and transient models.

Recharge components
Estimation of recharge at different segments of the irrigation system is necessary to set up different management scenarios.For this purpose, a water balance approach was utilized and input data for the water balance model were taken from remote sensing and secondary sources.Recharge was estimated individually from each irrigation system segment of LCC.The detailed methodology of recharge estimation for LCC irrigation subdivisions, including Sultanpur, can be seen in Usman et al. (2015).

Statistical analysis
Different statistical parameters were used to evaluate the model performance.These parameters include correlation coefficient (r), Nash Sutcliffe efficiency (NSE) and mean square error (MSE) as defined below: (4) (5) (6) where O is measured (observed) data; O m is the mean of measured data; S is simulated data and S m is the mean of the simulated data.

Calibration and validation of models
A steady-state model was developed and its calibration was performed in the first phase.The groundwater heads for April 2005 (pre-monsoon period) were used for manual followed by automated calibration operation using PEST.The hydraulic conductivity was taken as the adjustment parameter to match the model heads with observed heads.The values of hydraulic conductivity were retrieved from the bore-logs and used as the pilot points for PEST calibration.Before employing automated calibration using PEST, manual calibration was performed to bring simulated heads close to observed heads.This practice helps one to understand the model system better and also facilitates reducing the automated calibration convergence time considerably, thus providing a check against model failure.About 34 manual model runs were performed before the automated run.
Following the steady-state model calibration, a transient model was developed using groundwater heads of the steady model as initial conditions.For the transient case, both hydraulic conductivity and specific yield were considered as calibration parameters.The initial and parameter bound values for these hydrological parameters were taken from the literature against different depth of well-logs collected during drilling.Again, before running the model in automated mode, manual calibration was performed to facilitate the automation process and smooth convergence of the modelling system.Again about 15 manual model runs were performed.The automated calibration parameters for PEST operation were adjusted where needed and the model was run in 'estimation' mode.The performance of the calibration development was observed in the form of an objective function as residuals for each iteration of PEST run were noted.The model converged in 10 iterations and at some locations the residuals were higher, but dropped near to zero in 10 iterations.The model run time in this case was about 7 h on an Intel Core i5-2410 M CPU at 2.30 GHz with 8.0 GB installed memory and 64-bit operating system.For automated calibration of the transient model, data from October 2005 to 2009 were utilized.Figure 4 shows the development of the objective function in different iterations, whereas Fig. 5 (a,b,c) shows the locations of observation points and scatter plots for simulated and observed heads for both steady and transient models after the calibration process.
Different statistical indicators showed very good results for both the steady and transient models.In the case of the steady model, the values of NSE and MSE were 0.98 and 0.23 m, while for the transient case these values were 0.97 and 0.28 m with correlation coefficients of 0.98 and 0.97, respectively.After performing the model calibration, validation of the model was performed and the results can be seen in Fig. 5(c).The validation of model results for the period 2009-2012 showed that the model performed reasonably well, as NSE and MSE were 0.96 and 0.31 m, while the correlation coefficient was 0.97.In general, the model performed very well for the upper model regions as most of the points were placed very near to the 1:1 line, but relatively more variation was observed in the middle and lower model regions.The setting of initial hydrologic parameters, and setting of boundary conditions for the river, could be the reason for this variation.Nevertheless, overall the results were reasonable, to be adopted for further model simulations.

Sensitivity analysis
Variations in the hydrologic parameters and modelling simplifications are evaluated through sensitivity analysis of the modelling output.The sensitivity analysis was performed manually by uniformly multiplying hydraulic conductivity, specific yield and recharge by different multiplying factors.For each simulation, one parameter was changed at a time keeping other parameters unchanged.The multiplying factors used were 0.2, 0.4, 0.8, 1.2, 1.5, 2.0, 2.5, 3.0 throughout the interior nodes of the model and the model was rerun for each factor.
Figure 6 shows that higher sensitivity was observed for changes in hydraulic conductivity than recharge in the case of the steady model.The difference in head was higher for small multipliers up to 1.5, which became small for higher factor values.The sensitivity was less for smaller recharge factors, which became higher for greater factor values.
For the transient case (Fig. 7), sensitivity was observed to be relatively higher for hydraulic conductivity followed by specific yield and recharge.The sensitivity was found to be relatively higher for large multipliers for both recharge and specific yield.However, positive and negative changes in groundwater heads are observed for recharge and specific yield, respectively.The possible reason is availability of more water for the modelling domain from groundwater recharge and low storage due to higher values of specific yield.

Scenario analyses
To explore the future trends of groundwater behaviour under different management scenarios, simulations were run for a one-decade time period.Different weather conditions (i.e., precipitation) were considered for running the simulations, including normal weather (i.e., precipitation) conditions, very wet summer and winter conditions, and dry summer and winter conditions.values were compared to each season and year separately and three types of weather conditions identified, including normal weather conditions, weather conditions with higher rainfall and weather conditions with lower rainfall than average.Under each weather condition, five hypothetical water management scenarios were run, which included a 50% reduction in canal seepage in a 10-year period (1st), 50% reduction in field percolation in a 10-year period (2nd), 50% reduction in canal seepage and field percolation over 10 years (3rd), 25% increase in groundwater pumping over 10 years (4th), and 25% decrease in groundwater pumping over 10 years (5th).The baseline time considered for all the scenarios was 15 September, 2011 (i.e.Julian Day 40801).
The groundwater levels were never stable during the calibration and validation periods with a decrease and increase in trends during the first and second halves of the simulation period, respectively.Figure 8(a) shows that under normal weather conditions, groundwater levels decreased steadily in all of the scenarios, except for the 5th scenario, where groundwater level first decreased but then remained stable for the rest of the simulation period.In this scenario, the  In the second type of scenario run (Fig. 8b), ignoring very wet summer and winter weather conditions, groundwater drop was accelerated and a groundwater drop of about ~3.0 m was observed for all scenarios except the 5th scenario.However, under such weather conditions, the groundwater dropped below the level at the start of the model simulation run (i.e.Julian Day 38641).
For scenario runs considering very wet summer weather conditions (Fig. 8c), the groundwater remained almost stable within a range of about ~ 1.0 m for all scenarios except the 5th scenario, where groundwater developed positively during the simulation period.In all types of weather conditions, the drop in groundwater level was greatest in the 4th scenario, which indicate that if groundwater pumping should increase without introducing any water management intervention then groundwater resources could be at risk of diminishing.In the present study, all of these scenarios were run separately; however, results indicate that groundwater drop could be minimized if different scenarios were run in combination, which encourages the implementation of different on-farm and water management practices in the region.The increase in groundwater levels during wet years further strengthens the possibility of sustainable water management in this area.The possible water management solutions may include replacing of high delta crops like rice and sugarcane with less water consumptive crops.Moreover, water allowance needs to be revisited and cropping patterns may be adjusted to minimize the gap between demand and supply of canal water.

Figure 1
Figure 1 Location of the study area (in circle) in Pakistan, b) details of irrigation network

Figure 2
Figure 2Hydrogeology of the study area

Figure 3 Figure 4
Figure 3Flow diagram for current modelling study Figure 5(c) Validation of transient model

Figure 5
Figure 5 (a,b) Locations of the observation points and calibration results for steady and transient model

Figure 7
Figure 7Sensitivity analysis for transient model

Figure 6
Figure 6Sensitivity analysis for steady model