Inverse Parametrization of a Regional Groundwater Flow Model with the Aid of Modelling and GIS: Test and Application of Different Approaches

Abstract

The use of inverse methods allow efficient model calibration. This study employs PEST to calibrate a large catchment scale transient flow model. Results are demonstrated by comparing manually calibrated approaches with the automated approach. An advanced Tikhonov regularization algorithm was employed for carrying out the automated pilot point (PP) method. The results indicate that automated PP is more flexible and robust as compared to other approaches. Different statistical indicators show that this method yields reliable calibration as values of coefficient of determination (R²) range from 0.98 to 0.99, Nash Sutcliffe efficiency (ME) range from 0.964 to 0.976, and root mean square errors (RMSE) range from 1.68 m to 1.23 m, for manual and automated approaches, respectively. Validation results of automated PP show ME as 0.969 and RMSE as 1.31 m. The results of output sensitivity suggest that hydraulic conductivity is a more influential parameter. Considering the limitations of the current study, it is recommended to perform global sensitivity and linear uncertainty analysis for the better estimation of the modelling results.

1. Introduction

Groundwater supply is compulsory for the subsistence of agriculture in Pakistan. It is believed that about 40% of irrigation needs in Punjab, the main food-producing province in the country, are met from groundwater [1,2]. However, its present unregulated and uncontrolled application is repleted with serious threats and consequences to water resources [3]. Effects of ample use of groundwater can discern in form of declined water table in most of the canal commands in Punjab and Sindh provinces [1,4]. Rapidly falling groundwater levels are causing deterioration of its quality and lowering of crop yields against their potential levels [5], and hence need to be addressed by adequate management strategies.

Groundwater modelling could be a suitable tool to represent the natural groundwater flow in the environment and to predict the fate and movement of solutes and potential contaminants under natural or hypothetical scenarios [6]. Such models can be used to predict the effects of groundwater abstraction and irrigation developments with regard to the aquifer and to simulate various water management scenarios [7]. Nevertheless, groundwater modelling involves a variety of data requirements, whereas poor data availability may lead to fallacious results.

The process of model calibration is mandatory for acquiring reliable modelling results from prediction models. A good model calibration ensures that residuals between measured and computed data are minimized and parameter uncertainties are low. This can be verified by subsequent model validation (i.e., run the model with data that are not used for calibration). According to [8], a good calibration should only be acknowledged under a wide range of hydrologic conditions. Model calibration can be broadly classified into (1) traditional manual calibration [9] and (2) automated calibration using inverse techniques [10]. Historically and under consideration of the available computational power, model calibration was manually performed by adjusting model parameters with a more or less trial and error process. This approach has been widely used for several complex regional models, e.g., [11,12,13]. However, manual calibration is time-consuming, tedious and very subjective in nature [14]. The accuracy of results from this approach requires good experience of the modeller and thorough understanding of the system and is also characterized by the strategies employed to adjust the model parameters [13]. In recent times, manual calibration has been partly substituted by automated approaches, which are also recognized as nonlinear parameter estimation techniques [15]. The use of these techniques in groundwater modeling is now commonplace [16,17,18], because of their high speed to determine best fit parameters, their low subjectivity in the calibration procedure and the availability of generic software that can be easily linked to different hydrological models. It is very likely that the results from automated calibration techniques may be better than from manual techniques along with the simultaneous possibility of carrying out sensitivity and uncertainty analyses [19].

Despite the fact of advantages of automated calibration techniques, their application is still limited by a number of factors, which are well represented by [18] in the context of United Kingdom regional groundwater flow models. However, the concerns are common in many parts of the world, especially where field data sampling is not satisfactory. One notable factor is the limited computational capacity for regional models with distributed data inputs [19]. Moreover, such models generally possess high spatiotemporal variability which causes enormous nonlinearity and high correlation between different model input parameters. The result may be unrealistic model parameter distributions and also the calibration process descends into one local minimum without exploring other depressions [20]. This problem is pronounced in local search strategies as compared to global search strategies [8,20]. Nevertheless, computational costs are remarkably reduced by local strategies to explore objective function minima [8]. According to [13], it is very difficult to conclude about the effectiveness and efficiency of both techniques as their performance varies with the intended use of a model.

Uncertainty may not only emerge from imprecisely known model parameters, but also from inadequate simplification of the real system. In particular, calibration of any groundwater model is based on some method of spatial parameter characterization or zonation. This usually implies to subdivide the model domain into different hydrogeological units based on geological properties and other evidence. Then uniform hydraulic properties are assumed for each zone [21]. This simplification means that values of hydraulic properties at any location are weighted averages of their true values over a large area [18]. This use of discrete zones for hydraulic parameters may result in unnecessary structural uncertainties of the model system and, therefore, produce more unrealistic parameter values and also present geological heterogeneity in “unnatural” appearance [15]. An alternative approach can be the PP method [22], which interpolates hydraulic property values from a set of points distributed throughout the model domain [15]. The result is a smoother parameter distribution with reduced structural uncertainty as compared to conventional zone-based parameter assignment approaches. Nevertheless, the use of a larger number of pilot points could result in enhanced model run times. However, implementation of advanced regularization technique like Tikhonov regularization [23] can overcome this potential deficiency [24]. Further, an increase of hardware performance (e.g., multi-core processors, incorporation of graphics card processors) could facilitate the application of inverse methods for groundwater modelling [25].

Model calibration either from manual or automated techniques does not guarantee the full reliability of model predictions as groundwater flow and transport simulations are never closed systems and both independent and dependent variables of such systems are laden with inference and assumptions [26]. Parameter values obtained during calibration are only as realistic as the validity of the model assumptions and hence modelling is only a reflection of the real physical process [27,28]. Sensitivity analysis is a valuable tool to address uncertainties through identification of important model input parameters and to quantify the corresponding uncertainty in model output [29]. The results of sensitivity analysis are also important with regard to allocation of model sampling sites, a better understanding of the modelled system, improvement of the calibration process, and thus, model validation and reduction of uncertainties [30,31]. Sensitivity analysis without utilization of automated techniques (i.e., the manual variation of different input parameters one by one while keeping all other parameters unchanged) cannot truly reflect the real system behaviour. This approach exhibits model responses under non-calibrated conditions because changes in only one parameter at a time while keeping other parameters fixed can bring the model to non-calibrated state. This deficiency is well tackled by the use of automated calibration approaches which perform the sensitivity analysis of different model input parameters one by one while keeping the model in the calibrated state [18,32].

In the current study, the inverse parameter estimation tool PEST [32] was used for the automated calibration and sensitivity analysis of a regional groundwater flow modelling in the irrigated agricultural region of the Lower Chenab Canal (LCC), Rechna Doab, Punjab, Pakistan. The study aims to compare different model calibration techniques along with the use of some advanced regularization techniques for automated calibration methods. Detailed analyses of the sensitivity of model parameters are performed and areas of particular modelling importance are identified. Sensitivities of model results to different locations of observation points are also investigated and results are spatially presented for easy understanding. The arrangement of this manuscript is that the general information of the model region is described first; secondly, it presents the conceptual model and the theory of modelling approaches and, finally, the presentation of results and necessary discussion is done.

2. Description of the Study Area

LCC, Punjab, Pakistan has been chosen as the study region (Figure 1). The LCC irrigation system originates at the Khanki headworks which distribute water to its eastern and western parts through seven branch canals. This irrigation system was designed in 1892–1898 and its command area lies in Rechna Doab which comprises of the land mass between rivers Ravi and Chenab. The location of the area is between latitude 30°36′ and 32°09′ N and longitude 72°14′ and 77°44′ E. The present study mainly focuses on the eastern part of LCC. Two link canals namely Qadirabad-Balloki (QB) and Trimu-Sidhnai (TS) flow from north to south and fall into river Ravi. A major part of LCC east lies in the districts of Faisalabad and Toba Tek Singh. Administratively, the entire study area is split into 9 irrigation subdivisions; Chuharkana, Paccadala, Mohlan, Buchiana, Tandlianwala, Kanya, Tarkhani, Bhagat and Sultanpur. Irrigation subdivision is considered as the smallest management unit of the irrigation system in LCC. The structuring of these irrigation subdivisions ensures the equitable distribution of canal water among different consumers in the area.

LCC is mainly categorized as an agricultural area with a comprehensive irrigation canal network. Many different types of crops are grown throughout the year including rice, wheat, sugarcane, cotton, rabi fodder, kharif fodder etc. The whole cropping year can be subdivided into two seasons called kharif and rabi. The kharif season generally starts from May and ends in October, while the rabi season prevails from November to April. Rice and wheat are the two major crops during kharif and rabi seasons, respectively. The other crops cultivated during rabi season are rabi fodders (mainly barseem and oat), while cotton and kharif fodders (mainly sorghum, maize and millet) are grown in kharif season. Sugarcane is the annual crop which is cultivated in the months of September and February [33].

The climate of the area is arid to semi-arid with four different types of weather seasons including summer, winter, spring and autumn. The summers are hot and long lasting with temperatures fluctuating between 21 °C and 50 °C. During winters, daytime temperature ranges between 10 °C and 27 °C, whereas night temperature may drop to zero. The average annual precipitation in Rechna Doab varies from 290 mm in the south-west to 1046 mm in the north-east. Highest rainfalls occur during the monsoon period from July to September and account for about 60% of average annual rainfall [33]. There are three main weather stations in/near to the current study area operated by Pakistan Meteorological Department (PMD), which include Lahore (LHR), Faisalabad (FSD) and Toba Tek Singh (TTS).

3. Development of Conceptual Model

3.1. Processes

Figure 2 represents the conceptual model of the irrigation system in LCC. The model region is bound by two link canals on the eastern and western sides and along most of the southern side by river Ravi. The water flow from river Chenab to river Ravi is perennial as water rights of river Ravi are possessed by India under the Indus Water Treaty signed between the two countries in 1960. Therefore the river Ravi cannot cater the irrigation needs of the region alone. The northern boundary of the model domain acts like inflow/outflow at different sections, and the lateral flows were estimated by Darcy’s law using piezometric data [34]. Overall, the groundwater flow is parallel to this model boundary from north-east to south-west direction.

The irrigation network in LCC is quite large and comprised of main canals, branch canals, major distributaries and minor distributaries which are fed from Khanki headworks at the river Chenab. The transport of water through this irrigation network contributes a major recharge to groundwater at different stages of the irrigation system.

The elevation in the study region drops smoothly from north-east to south-west direction and hence causes regional groundwater flow movement along this drop. A large area is sown under different crops which are irrigated by canal water along with the support of groundwater pumping. As the cropping intensity in the region is quite high, and less rainfall is observed except during the monsoon season, the share of groundwater in irrigation needs is quite high. The major share to groundwater recharge takes place from monsoon rainfalls and also from field percolation losses. The climate of the region is mainly semi-arid with higher air temperatures especially during summers which cause an extensive loss of irrigation water to the atmosphere in form of evapotranspiration.

3.2. Hydraulic Properties and Geological Scheme

The study area is a part of an abandoned floodplain. The deeper part is formed by the underlying metamorphic and igneous rocks of Precambrian age. The area is underlain by highly stratified unconsolidated alluvial material composed of sands of various grades interbedded with discontinuous lenses of silt, clay and nodules of kanker, a calcium carbonate structure of secondary origin deposited by present and ancestral tributaries of the Indus River [4]. This forms the only one continuously interconnected aquifer and sediments at top mainly comprise of medium to fine sand, silt and clay. Overall, high percentage of silt and fine to very fine sand is dominant in the area with clay dominant in some depressions [35]. The origins of clay have not been identified specifically, but they are presumed to be the repeatedly reworked loess deposits of the hills at the north and northwest. Hydrogeological investigations in Rechna Doab were carried out during the 1957–1960 period wherein test holes were drilled throughout the Doab. The maximum thickness of alluvium is not accurately known although the logs of test wells show that thickness is over 150 m nearly everywhere (Figure 3).

Figure 3 illustrates the location of well-logs in LCC, while Figure 3 indicates lithological details of selected well-logs at different locations, according to which the thickness of the alluvium complex is relatively higher in the lower LCC parts compared to upper parts. It can also be seen that the aquifer is mainly composed of sand with deposits of clay, gravel and silt at different depths with no typical pattern in the arrangement of these materials. The aquifer material is highly porous and is capable of readily storing and transmitting water. The horizontal permeability is an order of magnitude larger than the vertical [36]. The total porosity of the water-bearing material ranges from 35% to 45% with an average specific yield of around 14%. The results of pumping tests and the lithological and mechanical analyses of test holes are presented by [37], according to which hydraulic conductivity varies from 24 m/day to 264 m/day and specific yield values vary from 1% to 33% in Rechna Doab.

3.3. Data Collection

3.3.1. Maps, Shape Files of Natural Features

The information about location and boundary of the model area is fundamental in the setup of the numerical model. The information regarding geometric features like boundaries of the model area, rivers, canals, pumping wells etc., is generally required in form of shape files. These shape files were obtained from Punjab Irrigation Department (PID), Pakistan.

3.3.2. Elevation and Well Log Data

For describing the geometry of the 3-D groundwater flow model, elevation data are required which were derived from a Digital Elevation Model (DEM) and from well log information. The Shuttle Radar Topography Mission (SRTM) digital elevation data of 90 m spatial resolution were retrieved cost-free from the website of the Consultative Group for International Agricultural Research (CGIAR) Consortium for Spatial Information (CGIAR-CSI) (http://srtm.csi.cgiar.org) and well log data were collected from Salinity Monitoring Organization (SMO), Pakistan.

3.3.3. Material Properties and Model Parameters

Piezometric water levels and river/canal gauges were attained from SMO and PID, respectively. The data of groundwater pumping and recharge were not readily available from any source and, hence, were estimated indirectly by utilizing different techniques described in detail by [34]. The major information about material properties in the saturated zone includes hydraulic conductivity, transmissivity, storativity, drain/fillable porosity etc. Soil texture details for different depth profiles were available from well-logs which were utilized to manipulate different hydraulic parameters against each material class [38].

4. Development of Mathematical Model

Figure 4 shows the detailed description of the adopted methodology for the current modelling study. The brief description of some important parts is given below.

4.1. Theory of Groundwater Flow

The groundwater modelling system FEFLOW v6.1 was used for the modelling procedure. FEFLOW has been successfully tested and applied in a number of benchmark studies around the world [39]. A 3-D finite element groundwater flow model is set up for the current study. The fundamental basis of FEFLOW is that it introduces the Darcy equation in the mass conservation equation to represent the mass conservation of any phase [39]:

$$\frac{\partial }{\partial \mathrm{t}}\left({\mathsf{\epsilon }}_{\mathsf{\alpha }}{\mathsf{\rho }}^{\mathsf{\alpha }}\right)+\frac{\partial }{\partial {\mathrm{X}}_{\mathrm{i}}}\left({\mathsf{\epsilon }}_{\mathsf{\alpha }}{\mathsf{\rho }}^{\mathsf{\alpha }}{\mathrm{V}}_{\mathrm{i}}^{\mathsf{\alpha }}\right)={\mathsf{\epsilon }}_{\mathsf{\alpha }}{\mathsf{\rho }}^{\mathsf{\alpha }}{\mathrm{Q}}_{\mathsf{\rho }}^{\mathsf{\alpha }}$$

(1)

where t is time [T], ${\mathsf{\epsilon }}_{\mathsf{\alpha }}$is volume fraction of $\mathsf{\alpha }$-phase [-],${\mathsf{\rho }}^{\mathsf{\alpha }}$ is density of $\mathsf{\alpha }$-phase [ML⁻³], X_i is Eulerian spatial coordinate vector [-], ${\mathrm{V}}_{\mathrm{i}}^{\mathsf{\alpha }}$ is velocity vector of $\mathsf{\alpha }$-phase [LT⁻¹], ${\mathrm{Q}}_{\mathsf{\rho }}^{\mathsf{\alpha }}$ is mass supply of $\mathsf{\alpha }$-phase [T⁻¹].

Further, the equation was modified for the mass conservation of water with uniform density and viscosity:

$$\frac{\partial }{\partial \mathrm{t}}\left({\mathsf{\epsilon }}_{\mathrm{f}}\right)+\frac{\partial }{\partial {\mathrm{X}}_{\mathrm{i}}}\left({\mathrm{q}}_{\mathrm{i}}^{\mathrm{f}}\right)={\mathrm{Q}}_{\mathsf{\rho }}^{\mathrm{f}}$$

(2)

The term ${\mathrm{q}}_{\mathrm{i}}^{\mathrm{f}}$ is known as the Darcy flux, which can be described for water with uniform density and viscosity as below:

$${\mathrm{q}}_{\mathrm{i}}^{\mathrm{f}}=-{\mathrm{K}}_{\mathrm{ij}}\left(\frac{\partial \mathrm{h}}{\partial {\mathrm{X}}_{\mathrm{j}}}+{\mathsf{\rho }}^{\mathrm{f}}{\mathrm{e}}_{\mathrm{j}}\right)$$

(3)

where ${\mathrm{K}}_{\mathrm{ij}}$ is tensor of hydraulic conductivity [LT⁻¹], ${\mathsf{\rho }}^{\mathrm{f}}$ is water density [ML⁻³], and ${\mathrm{e}}_{\mathrm{j}}$ is gravitational unit vector [-].

4.2. Setup of Numerical Model Using FEFLOW

4.2.1. Mesh Generation and Setting of Modelling Problem

FEFLOW involves a number of steps to set up a groundwater model which includes designing of super-element-mesh, generation of mesh and setting of modelling problem. In super-element-mesh design, conceptual definition of geometry features of lines and points are included which is followed by mesh generation. For the current problem, the triangular mesh algorithm [40] was selected owing to its fast speed and its capability to accommodate complex setups of polygons, lines and points. After refinement of triangular meshes, they were tested by Delaunay criteria for ensuring maximum stability of the model. Saturation case with a projection of 3-D phreatic aquifer was selected as the problem class. Each mesh of the model domain contains three dimensions with six nodes per element. The total number of elements per layer and nodes per slice were 24,137 and 12,674, respectively. The whole aquifer was subdivided into 4 layers and 5 slices. The upper-most layer is set as phreatic which means a fixed slice topping an unconfined layer. The other slices are set as dependent and the bottom slice is kept fixed. Different points in favor and against of multi-layer model were considered before its construction. The points in opposition include, the aquifer is very homogeneous and do not have a consistent clay layer. Moreover, only quantitative analysis was required without considering mass transport. Whereas the arguments in favor include as no prior information was available in literature on comprehensive groundwater system of the study region; secondly, the upper model regions show strata of relatively more heterogeneity; third sensitivity and parameter error in different boring depths was required. Both steady and transient models were prepared and groundwater heads of the steady-state model were utilized as initial conditions for the transient model. Constant and Adams-Bashforth time stepping schemes were used for steady and transient models, respectively. The Adams-Bashforth scheme is an automatic time step control scheme which gives easy convergence and stability to model runs. The PARDISO—Parallel Direct Solver [41] was chosen for solving the equation system.

4.2.2. Regionalization of Hydraulic Properties

Refer to Figure 3, there is no definite pattern in the arrangement of aquifer materials, therefore, the selection of a suitable regionalization method was crucial. The Akima interpolation algorithm [42] with linear interpolation type, with five neighbours and, zero over/undershooting criteria was used for interpolation of parameters for the whole modelling domain. Interpolation was performed on presumably logarithmized values as this can lead to better regionalization results in case of log-distributed parameters such as hydraulic conductivity and effective porosity. Figure S1 depicts the distributed hydraulic conductivity in different model layers along with information of maximum interior angles of mesh triangles and Delaunay criteria for the final model setup. Smaller angles of triangles are always considered better as they do not result in a violation of Delaunay criteria. A triangular finite element violates the Delaunay criterion if the circumcircle of the triangle includes a node not belonging to the finite element. Delaunay-compliant triangulations maximize the minimum angle of all the triangles in the mesh; elements with large angles that are potentially leading to instabilities are avoided. Elements violating the criterion get a value of 1 in this parameter; triangles being consistent with the criterion have a value of 0.

4.2.3. Setting up Different Model Boundary Conditions

Recharge is generally considered as a boundary condition for groundwater models but in case of FEFLOW, it is treated as a material property. As described earlier, recharge information was not available from any secondary source, therefore, it was estimated for a period from 2005 to 2012 at a spatial scale of 1 km × 1 km using water balance approaches and remotely sensed data. The detailed recharge estimation procedures and results of each water balance component can be read from [34]. Recharge is assigned as a time-variant material property to the transient FEFLOW problem.

Prescribed head boundary conditions were assigned along the QB and TS link canals based on their historical water level records [4]. The gauge data for river Ravi were available only at headworks of Balloki and Sidhnai, however, the distance between them is quite long. Therefore, water levels in the course of the river were carefully worked out with DEM before their assignment to the model. First kind or Dirichlet type, time-varying boundary condition was applied to the river instead of Cauchy type boundary condition. The Cauchy boundary condition requires hydraulic conductance of river bed material and its geometry which was inaccessible from any source. Nevertheless, the flow in the river is not dynamically variable and it is also flowing at lower elevation in comparison to the model domain which could minimize the effects of this boundary condition along the lower model boundary. The other boundary conditions applied include Neumann boundary conditions mainly to the northern boundaries of the model and also to parts of the southern boundary which were not in contact with the river Ravi. The inflow and outflow sections and flow to the model domain were identified using contours of piezometric water levels and by application of Darcy’s law. The detailed methodology of the procedure is accessible from [34]. The final boundary condition assigned to the model domain was groundwater pumping which was assigned as nodal source/sink type boundary condition and the information of groundwater pumping was taken from [34].

4.3. Model Calibration and Parameter Estimation

4.3.1. Functionality of PEST for Model Calibration and Parameters Sensitivities

PEST (which is an acronym for Parameter ESTimation) software [32,43] is a model-independent software for inverse parameter estimation and was used as a primary tool for optimization of the numerical model. The purpose of PEST is to assist in data interpretation, model calibration and predictive analysis [32]. PEST applies the Gauss-Marquart-Levenberg algorithm, which combines the advantages of the steepest descent and the inverse Hessian methods. This algorithm provides more efficient and faster convergence towards the objective function minimum [8].

The goodness of fit of model results to the observed data can be estimated and presented in form of residuals between measurements and simulated results. The sum of squares of these residuals is known as the objective function ($\mathsf{\Phi }$). Φ in groundwater models typically comprises of many different types of data, for instance, hydraulic heads, leakage rates or gauged flows [18,44]. The search algorithm used in PEST changes the model parameters until a minimum $\mathsf{\Phi }$ is achieved. During PEST execution, the user observes two steps per iteration which include the derivative calculation and the adjustment of parameter values aiming to reduce the objective function. The parameter estimation process requires a set of model runs to calculate the sensitivity of the simulated values to changes in each variable parameter [28]. The details of these sensitivities are recorded in form of a matrix, known as Jacobian matrix. PEST carries out a number of optimization iterations by which it attempts to vary the model parameter values to reduce the value of the objective function. The trialled changes made to any parameter are proportional to the respective sensitivity of the model results. Model runs are repeated with continuously updated parameters until no further improvement in model results can be obtained. The sensitivity matrix is recalculated at the start of each optimization iteration utilizing the former “best” parameter set [45]. This process continues until it reaches either the user specified targeted objective function, the maximum number of optimization iterations allowed, or execution fails.

4.3.2. Pilot Point Calibration Technique

The pilot point technique [22,46] defines parameters as a spatially variable distribution. In classical optimization approaches, the common assumption to each geologic formation is that it has spatially constant values, which is rarely true especially for regional models. To overcome this problem, the distribution of hydraulic properties within the model domain is described by a set of points at particular locations (the pilot points). A number of these pilot points are introduced to the model domain and PEST is asked to estimate the hydraulic properties of the aquifer at each such point by minimizing the objective function. These “point-hydraulic-properties” are spatially interpolated in all active cells/elements within the model domain using geostatistics (i.e., kriging).

Individual pilot points can be assigned to different zones within the model domain. Only those points assigned to a particular zone can be used in calculating hydraulic property values throughout that particular zone. Furthermore, the variogram can be different for each zone reflecting differences in the geology expected within each geological unit. For the current problem, each slice/layer of the model was treated as a single zone for performing the interpolation. Figure 5 shows the location of each pilot point used in the present inverse model calibration, which was set up according to the procedure defined by [15,32].

Conventional wisdom and the principle of parsimony dictate that the numbers of parameters involved in a parameter estimation process should be kept to a minimum. However, when using pilot points in conjunction with PEST’s “regularisation” mode, the opposite is often true. The use of pilot points in characterizing the spatial distribution of a hydraulic property was done through kriging. One of the benefits of kriging is that the factors by which hydraulic properties at the pilot points are multiplied to obtain the hydraulic property value are independent of the actual hydraulic property values at the pilot points. Hence a set of “kriging factors” pertaining to each of the cells/elements was calculated in advance of the actual interpolation process. The latter was undertaken again and again as the model was run repeatedly by PEST and it was therefore not necessary to repeat the calculation of the kriging factors within each model run. This resulted in a large saving of model run times to complete the overall parameter estimation process [10]. The added advantage of using kriging is its further use for the regularization processes based on the same variogram used for kriging [15].

4.3.3. Tikhonov Regularization

Application of regularization helps to reduce non-uniqueness in the parameters estimation. For the current study, Tikhonov regularization is implemented as the numbers of estimation parameters were quite large (360 pilot points in total). The Tikhonov regularization involves a number of “information” equations to define the initial values of the parameters; these equations can be defined by the user as prior information. The calibration process with Tikhonov regularization is formulated as constrained minimization process which minimizes the objective function to some user-specified target. If this condition is not met, PEST minimizes the objective function without considering conditions specified by the user and in the meantime, it adjusts weights to prior information. The use of prior information in PEST is quite sensitive especially for non-linear models (most models fall into this category) because specifying prior information usually results in parameter estimates that are close to values specified in the prior information [47].

4.4. Statistical Analysis

In order to evaluate the performance of the FEFLOW model, statistical measures were taken to quantify the differences in the measured and calculated groundwater heads. In the present study, the coefficient of determination (R²), root mean square error (RMSE), percent bias (PBIAS), and Nash Sutcliffe/model efficiency (ME) were used [34,48,49].

5. Results and Discussion

5.1. Selection of Calibration Parameters and Their Initial Values

The purpose of the reliable calibration process is to develop a model which reasonably represents water balance components including groundwater flow, recharge, and discharge, and reasonably matches observed groundwater levels [44]. The specifications of the optimization algorithm in PEST include model parameterization, the selection of parameters and defining their feasible range, assigning prior information to parameter groups and assigning weights to observation groups [19]. There is no rule of thumb about the decision on parameters for performing calibration, and finalization of calibration parameters is purely dependent on the modelling situation and available data [50,51].

It is noticed that if some calibration parameters are used with wrong limitations in inverse parametrization then it may still result in the best match to observations. This indicates a potential deficiency in model conceptualization [9]. Therefore, both the closeness between the observed and simulated conditions and the extent to which important aspects of simulation are incorporated in the model are important in the evaluation of model calibration. For the same reason, the number of calibration parameter types in the current modelling study is kept low. Some dynamically spatiotemporally variables like groundwater recharge and groundwater pumping along with lateral inflow/outflow were estimated through separate approaches with fair reliability [34], and they were not considered as calibration parameters. Hydraulic conductivity and effective porosity were only selected as primary calibration parameters with equal weights for all the points as according to [9], most models that are used to understand the past, understand the present, or to forecast the future are calibrated by matching observed groundwater heads. Both of the hydraulic conductivity and effective porosity was log transformed in accordance with the fact that most studies cited in the groundwater literature, which treat these parameters as regionalized variables, indicate that their distribution is better described by a log variogram [52].

Any initial model parameter values can be used for automated optimization operation, nevertheless, model fitting under such conditions represent untrue parameter values as compared to the real situation and thus it can affect both optimization and sensitivity analysis. Therefore, for current automated and manual calibration processes, site characterization information (local) was used as explained by [53]. In the first phase, the model was calibrated manually in steady state by altering values at the pilot points. The groundwater levels of October 2005 were used as initial condition for executing steady state calibration. The values of hydraulic conductivity are adjusted to bring calculated heads close to observed heads. Different statistical indicators show reasonable agreement between calculated and observed heads as R², ME, RMSE and correlation coefficient are 0.99 m, 0.98 m, 1.58 m and 0.98 m, respectively. For automated calibration, parameter values from the manually calibrated steady model (modeller’s pre-calibration expert judgment) were used as initial conditions in transient simulations as reported by [9]. Figure 6 (upper half) shows the distributed initial hydraulic conductivity and effective porosity utilized in the automated PP optimization at the head, middle and tail reaches of the 3-D groundwater model. For zone-based parameterization, separate zones were demarked for each model layer and distributions of initial hydraulic conductivity and effective porosity for this case at three model locations are represented in Figure 6 (lower half). The parameters bound set for all pilot points and zones ranged between 10⁻⁶ m · day⁻¹ to 300 m · day⁻¹, and 10⁻⁷ m · day⁻¹ to 30 m · day⁻¹, for horizontal and vertical conductivities, respectively, and for effective porosity it was ranged from 0.009 to 0.250.

The final values of variables for PEST operation can be found in Table S1, including lambda iteration, parameter variation, derivative method switch, termination criteria and optimization statistics. These values yielded smooth and successful convergence of the model setup.

5.2. Calibration and Validation Results for Transient Model under Pilot Points

The modelling results suggest that groundwater is fully connected in the region and it could have been modelled with 2D settings as well. Thus upper model layer could be well representative in explaining most of the current modeling results. Figure 7 demonstrates the development of the objective function with each iteration of the automated optimization. About 585 observations of groundwater heads were used for controlling the calibration process of the transient groundwater model. The objective function decreased from 30,932 m² to 1837 m² for PP optimization. In the initial phase, the decrease in objective function was greater but as the optimization process advanced, lower changes were observed. It is to be noted that PEST here estimated the best parameter values by minimizing the sum of squares of the differences between measured and simulated groundwater heads.

The performance of the groundwater flow model was evaluated by comparing the groundwater heads simulated by FEFLOW with measured piezometric heads calculated with reference to mean sea level (requirement of FEFLOW—3D model). The data of the piezometers from winter 2005 (post-monsoon) to winter 2008 (post-monsoon) were utilized for the model calibration. The simulated and measured heads were drawn on a scatter plot (Figure 8a), which indicates the majority of the points fall on/near to the main diagonal thus show successful calibration of the model. Moreover, Figure 8b represents the spatial comparison between simulated and measured groundwater heads for two different simulation times, one during kharif and other during rabi cropping season. At the majority of the locations, the observed and simulated groundwater heads are quite close except for some regions in the northern parts of the modelling domain. The possible reason for this relatively high difference could be due to the large number of tube wells in this region, which extensively pump groundwater, especially for rice cropping. Therefore, the possibility of occasional partial inflow through boundaries of the study area could not be well represented by the inflow/outflow boundary conditions. Similar types of modelling results were also reported by [4] in similar areas. Different error parameters including R², ME, PBIAS and RMSE indicate reliable calibration results. The value of R² for model calibration is 0.99. Similarly, ME, which is an indicator of the efficiency of calibration, is 0.976. This value indicates that deviation of simulated heads from measured heads is only 2.2%. The value of PBIAS is only 0.026 whereas the RMSE of the simulated heads from measured ones is 1.23 m.

Following successful calibration of transient groundwater flow model, it was validated for the piezometric data from summer, 2009 (pre-monsoon) to winter, 2011 (post-monsoon). The 1:1 chart for model validation is presented in Figure 8c. Same statistical measures as used for calibration were employed for validation of the simulated results. The value of ME is 0.969 which indicates that deviation of simulated heads from measured heads is only 3.0%. The values of PBIAS and RMSE are −0.205 and 1.31 m, respectively, pointing towards reasonable results from the initially calibrated transient model. It is to note that all of the comparison results presented herein were achieved by pilot point optimization.

5.3. Comparison of Calibration Results from Different Methods

Apart from the automated calibration of the groundwater model using PP, its comparison with conventional manual PP and zone-based models result is also performed. For zone-based manual calibration, the whole model domain was distributed into nine stratigraphic zones (Figure 6; model top view) based on well-logs information and textbook knowledge. Similar zones were demarked for each model layer which summed up to 36 zones in total for the current modelling problem.

The results of the manual calibration for the steady state model have been presented in the previous sections and the same model was extended to the transient state for both manual and automated optimizations using PP. Both manually calibrated, steady and transient, models showed relatively weak results in comparison to automated PP optimizations as the relatively higher error is observed for them. The statistical indicators, R², ME, PBIAS and RMSE, yielded values of 0.98, 0.964, 0.73 and 1.68 m, respectively, for manual calibration of the transient model with PP. The results of the manually calibrated zone-based model even produced weaker results as R², ME, PBIAS and RMSE were found to be 0.96, 0.934, 0.98 and 1.89 m, respectively. It is also to be noted that the manual calibration process of the transient models proved to be subjective and very cumbersome which resulted in the poor calibration of some model regions. For example, the simulations were run for about 58 and 45 times for PP and zone-based transient models before they were considered calibrated. Similarly, other researchers reported that they had to run the groundwater model for about 80 times to bring the difference between simulated and measured heads within 0.5% for manual calibration [52]. According to [51], it is common for manual calibration to make model runs from 20 to 50 times before acceptable results reach. It is also observed that model calibration in one particular model region could have brought major changes between measured and simulated heads in some other regions, which resulted in difficulties in model calibration. Nevertheless, manual calibration was proved helpful in advance of automated calibration to avoid a finding of local minima due to non-uniqueness of non-linear models by automated calibration, and also the number of iterations reduced remarkably from 13 to 7 to find minimum measurable objective function without failure of optimization process [19].

The calibration results of different methods also presented for a number of piezometers at different locations of the modelling domain to investigate the calibration performance spatially. Hydrographs of 6 different piezometers, two each at the head, middle and tail reaches of the modelling domain are presented in Figure 9. These piezometers are selected in a way to represent the whole model domain behaviour as locations of different piezometers are near to different boundaries of the modelling system like link canals, river, inflow and outflow boundaries. From the trends of these hydrographs, it is clear that automated PP calibration always yielded better results than manual calibration of both PP and zone-based model calibration. The differences of heads between manual calibrations and automated calibration were relatively higher in upper model regions as compared to lower model regions. This could be due to the possibility of occasional partial inflow/outflow from this region due to more pumping from this region as described earlier. The descent of residuals of piezometric heads was quicker in lower model regions as compared to upper regions, which is another indication of the previous statement. This also indicates that the model was well constructed especially for these regions. Similarly, the differences between hydrographs of observed and zone-based calibration were relatively larger at upper model locations. The possible major reason for this could be the use of generalized uniform hydraulic property values for zones which are not fully depicting possible local heterogeneous conditions. For instance, in the current problem the zones are mainly described based on the irrigation subdivisions which were redefined frequently for different model trials but, due to the larger heterogeneity of aquifer material especially in the upper model regions (refer to well-log information), no final combination with best results was achieved. The best calibration of the model could have been possible if the hydraulic parameter bounds, especially for hydraulic conductivity, were set even higher as the values of present upper bounds. Nevertheless, it would have been far away from the reality of the on ground strata.

The calibration at model locations near to the river (i.e., piezometer #236, Figure 9) was difficult by the manual zone-based approach but better calibration results were achieved by the PP techniques. It is reported by [34] that the current study area receives its major recharge during summer seasons due to more monsoon rainfalls. In winter seasons, both rainfall and canal flow decrease yielding low recharge in major parts of LCC. This phenomenon is well represented by the trend lines of measured heads which show an increase in height after summer seasons and falling trends after winter seasons. Both of the PP calibration methods (i.e., manual and automated) depicted this increasing/decreasing trends in groundwater heads but they were found missing in case of manual zone-based calibration approach. The difference between automated and manual PP calibration techniques is observed relatively less in the lower to middle reaches of the modelling domain.

5.4. Model Parameter Sensitivities and Parameter Error

A very useful feature of PEST optimization is automatic robust sensitivity analysis which was useful for user intervention to optimization as it helped to understand which model parameters were more influential and which were less influential and could be set to fixed values. The sensitivity outputs were easily available from PEST output files and were useful for guiding the calibration process through the exclusion of non-influential parameter and thus to avoid failing of PEST. The PEST sensitivity analysis also facilitated the optimization process by automating the tedious task of adjusting certain model inputs, running the model, reading the output of interest, recording their values, and then commencing the whole cycle again. According to [19], the results of PEST sensitivity analysis should be carefully interpreted, as the dimensionless scaled sensitivities depend on the parameter values and different initial parameter value sets lead to different model results [54].

PEST calculates a figure related to the sensitivity of each parameter with respect to all observations and lists the composite sensitivity to each parameter of all observation groups, as well as of each individual observation group. The composite parameter sensitivity of each observation group can be evaluated by calculating the magnitude of the respective column of the weighted Jacobian matrix with the summation confined to members of that particular observation group. The magnitude is then divided by the number of members of that observation group which have non-zero weights [32].

The sensitivities of all PP in different model layers were estimated for all calibration parameters (i.e., Hydraulic conductivity and effective porosity). However, Figure 10 shows only the results for selected pilot points to facilitate visibility and representation of the results. However, the complete results are submitted as a supplemental material in Table S2 and Figure S2. Overall, the sensitivity with respect to hydraulic conductivity is higher than for effective porosity. The sensitivity of model output to model parameters is found higher in most of the north-eastern parts in comparison to other model regions. It is even higher in 1st and 4th layer of the model domain. For PPs in south-west and south-east model regions (refer to Figure 5), the sensitivity with respect to hydraulic conductivity is higher for 3rd and 4th model layer. In case of effective porosity, there is no significant variation in sensitivity for different regions of the model domain, especially for 2nd and 4th model layer. Nevertheless, sensitivity is relatively higher for north-eastern model regions in the 1st model layer and for south-western and south-eastern regions of the 4th model layer. The spatial information of sensitivity of model output could be useful for a further decision about conducting data collection campaigns and to investigate impacts of abstraction in these areas [18]. It might also be an indication of problems with the model setup for such areas if parameter values have reached their upper bounds. Nevertheless, this is not a case for the current modelling setup.

It is also imperative to investigate the sensitivity results in accordance with parameter error results as the understanding of parameter influence is also critical to identify which parameters are not well defined by the observation data and user knowledge input to PEST [18]. Such analysis can help the modeller to identify regions where current field information is not sufficient and where future data assessments should be taken place in order to have better data. The parameter error is calculated as the difference between initial and optimized parameter values divided by the initial parameter value. For the current modelling problem, parameter error for hydraulic conductivity and effective porosity were estimated for all model layers which are provided as supplemental material in Table S2. Figure 11 shows the parameter error only for the selected layer. The evaluation for hydraulic conductivity showed that both sensitivity and parameter error is found higher in model layers 1 and 2 for pilot points in regions A, B and C (Figure 12 for parameter error of hydraulic conductivity in the first model layer). When this information is incorporated into model layer thickness then it is helpful in planning further field investigations (e.g., preferred pumping test site, drilled well logs etc.). Figure 11 is also marked with another region “D” where the parameter error is high but the sensitivity of current model output is low. This suggests that, as the model is not very sensitive to parameter values changes in such regions, therefore, there is very less potential for improvement in model results with further field investigations. Contrary to this, there are some regions in the 4th layer (upper model locations) where model sensitivity is quite high but in this case, the parameter error is low. This also guides researchers to expect less improvement in current modelling results with further field investigations at such locations. The other case is about moderate values of parameter error and model sensitivity for model layers 3 and 4. In such cases, improvement in the modelling results is possible but the decision about further field data collection is all about the purpose for which the model is constructed and of course on the availability of sufficient funds for fieldwork.

The higher sensitivity of parameters for upper model regions is attributed to the nature of groundwater flow through the model domain. As described earlier, the sensitivity with regard to effective porosity is generally not very high and less variable throughout the model domain which is also accompanied by lower parameter errors. From the current results of automated optimizations, it is obvious that equally good calibration of the model could be achieved from different parameter values due to model nonlinearity, model uncertainty and higher correlation between model parameters which is also reported by [44]. The problem is more pronounced for distributed models alike. However, this situation is really attempted to be countered by carefully applying available field information in describing of suitable initial parameter values. In many real modelling studies, it is quite usual that no evidence of initial parameter values is available. Consequently, it is a duty of the modeller to present multiple models which exhibit equally good results to predict scenarios related to decision making about aquifer management. The right conceptualization of modelling environment would also be valuable in this regard.

5.5. Sensitivities at Selected Observation Points

In the last section, we have demonstrated the results for parameter sensitivity analysis and demarked the regions where changes in the parameter values have strong/weak effects on the model optimization. We have also evaluated the deficiencies of available field information by utilizing sensitivity and parameter error to support future data assessments for further model improvements. However, in many cases, engineers need information on some local scales which demand the construction of local models with finer spatial discretization. For such models, generally, there is a need for more detailed field data which result in additional hydraulic tests. The decision about locations of such tests is crucial, both because of required time and money. The knowledge gained form sensitivity analysis performed from this particular study could help researchers to decide about important sites for conducting field assessments and to establish a successful model.

PEST estimates and keeps a record of the sensitivities for different observation points in the Jacobian matrix. An independent utility named “JROW2VEC” is utilized to extract this information from the Jacobian matrix for different individual observation points. According to [32], observation sensitivity is generally of less use. But if the general flow behaviour is not drastically variable in time and/or space then this information could be quite useful. In principle, sensitivity for all observation points could be deduced from the Jacobian matrix but for the simplicity of the results, it is presented only for few observation points (i.e., 22, 238, 242, 211 and 303), which are located at different regions of the modelling domain. Corresponding plots are quite informative due to their spatial representation of sensitivity results. Positive sensitivity values in such maps reflect the tendency of an increasing hydraulic head with a decrease in values of hydraulic parameters and vice versa. The sensitivity for each respective observation point is in reference to a particular model layer as can be seen in caption of Figure 12.

For instance, Figure 12a shows the sensitivity for observation point 22, which depicts that major regions of the model domain are not very sensitive to hydraulic head at this point but there are fragmented locations which show strong sensitivity either positive or negative at this observation point at all three different simulation times. Similarly, sensitivity distribution for observation point 242 can be observed from Figure 12c. This shows a very clear trend of sensitivity distribution as the majority of the model parts upstream of this location show positive sensitivity. It means if any hydraulic activity is to be performed at/near to this point; it must contain very good data about hydraulic properties in high sensitivity zones. The sensitivity distribution for observation point 211 (Figure 12d) has very clear information as two different behaviours of sensitivity can be observed both upstream and downstream of this location. Majority of the nearby downstream locations exhibit very high negative sensitivity and the majority of upstream location have weak to average positive sensitivity. This suggests that downstream areas need to be properly investigated for hydraulic properties if any hydraulic activity is to be planned at this location. The sensitivity distribution for observation point 238 does not show consistent results for all simulation times which means that flow trends at this location are variable from time to time and such locations can be identified by opposite arrow directions in Figure 12b. For such locations, the information of observation sensitivity is not very useful. The sensitivity distribution and demarcation of important location for each particular region for rest of the selected observation points can be seen from Figure 12. Moreover, the sensitivity in the south-east of model domain away from points 22, 238 and 242 could be associated to numerical artifact.

6. Conclusions and Outlook

Quality data are a prerequisite to corroborate model calibration and validation leading to reliable predictions. Frequently unavoidable simplifications may cause larger uncertainties in model outcomes. Model calibration methods range from manual to automated techniques. The application of automated techniques are gaining popularity and their advantages need to be addressed. The current study was conducted on a regional scale to compare both techniques. A comparison of manual and automated calibration techniques were performed. The following major conclusions can be drawn from this study:

• It is found that the automated pilot point calibration method is more flexible and robust in comparison to manual approaches using pilot point and zone-based parameterization of the model due to its lesser subjectivity on part of the modeller’s experience.
• Automated pilot point calibration results in a reliable model calibration and validation for a majority of model regions as different statistical indicators show reasonable values. For calibration of the transient case, the values of R², Nash Sutcliffe Efficiency, % BIAS and RMSE are 0.99, 0.976, 0.026 and 1.23 m, respectively, and for validation, the values are 0.987, 0.969, −0.205 and 1.31 m, respectively.
• Apart from the lower calibration efficiency, it is also observed that manual calibration is tedious and cumbersome due to more model runs to get reasonable results.
• The spatial comparison of model calibrations shows that the pilot point approach yields overall better results at different locations with some higher differences at upper locations as compared to zone-based model calibration.
• Parameter sensitivity analysis shows that overall hydraulic conductivity is more influential as compared to effective porosity. However, this sensitivity is quite variable for different model locations and model layers.
• Sensitivities and error parameter results also address limitations/deficiencies of current hydraulic field data and help to identify regions where further field investigations could be planned.
• Sensitivities of different observation points demark different regions of particular importance and therefore guide planners to perform field activities there in future.
• Present sensitivity analysis was performed by a local approach employed in PEST. For such methods, there is always a possibility that the entire parameter space might not be well represented which could be addressed in future by some global sensitivity analysis approach.
• Predictive analysis is another way to explore uncertainties of model results. For the current study, it was attempted; however, the use of such analysis is only limited to a well-posed problem which was not the case for the current model. If a problem is ill-posed, then it does not work because pertinent matrices become un-invertible. The only possibility then left is to explore non-linear uncertainty analysis options. PEST has provided utilities like PREDUNC and/or GENLINPRED for this purpose. Null space Monte Carlo and running model in “Pareto” mode could be alternative solutions. Hence, it is recommended to explore these different approaches in future studies.