Abstract

The simulation-optimization (SO) modeling approach can be effectively used for aquifer parameter estimation. In this study, a numerical approach based on meshless local Petrov–Galerkin (MLPG) method is used for groundwater flow simulation and coupled with particle swarm optimization model for optimization. The study deals with the identification of the most suitable model structure for a hypothetical heterogeneous confined aquifer from a number of alternate models using zonation method of parameter estimation. A range of alternate models starting from homogeneous to more complex model structures are considered for the zonation. Inverse modeling of different model structures is carried out based on weighted least square performance criterion. The suitable models are selected and reliability analysis ascertained by computing three parameters of composite scaled sensitivity, coefficient of variation, and confidence interval, and the best model is selected. Sensitivity of estimated parameters is investigated by considering different sets of head data involving possible measurement errors. The solutions are compared with another inverse model using the MLPG and Levenberg–Marquardt algorithm. Based on the results, it is found that the proposed methodology can be utilized in the estimation of different unknown parameters in a regional groundwater system.

INTRODUCTION

The inverse problem of parameter estimation concerns the optimal determination of the system's hydro-geologic parameters based on the observations collected in the spatial and time domains. Hydraulic conductivity or transmissivity is often considered to be the most important calibration parameter in inverse modeling. Other parameters such as storativity, specific yield, contaminant transport parameters, stream-bed characteristics, or flows at model boundaries can also be estimated (Keidser & Rosbjerg 1991; Anderman et al. 1996; Rajanayaka & Kulasiri 2001; Straface 2009; Kourakos & Mantoglou 2012; Yeh et al. 2015; Onyari & Taigbenu 2017). Generally, inverse modeling is based on the ‘estimation’ of unknown parameters with a single best solution compatible with the available measurements. The inverse modeling approach proposed in the present study is based on spatial partitioning (zonation).

In the zonation approach, the entire aquifer domain is divided into a finite number of zones in which the parameter value is assumed to be constant. In the case of real field problems, the zonal structures may not be known and then the model structure is often determined by trial and error to minimize the error between observed and simulated head values (Ayvaz & Aral 2007). The simulation modeling is often combined with optimization algorithms to get deterministic solutions corresponding to different model structures. However, due to the ill-posed nature of the inverse problems, when gradient-based optimization algorithms are used, some difficulties such as convergence problems may occur. Therefore, meta-heuristic algorithms are usually preferred due to their ability in finding the global or near global optimum solutions without the necessity of working with gradients, as well as requiring information on an initial solution. During the last two decades, several meta-heuristic solution algorithms have been proposed by many researchers, such as genetic algorithm, tabu search, simulated annealing, swarm intelligent techniques, etc. (Maier et al. 2014). Generally, based on prior studies, real field data, computational tools' availability, and the user's discretion, the optimization algorithm is selected. The parameter structure identification procedures have been provided in earlier reviews and evaluations by Sun & Yeh (1985); Carrera et al. (2005); Hendricks Franssen et al. (2009); Kourakos & Mantoglou (2012), and in textbooks by Sun (1994) and Hill & Tiedeman (2007).

In the present study, we adopt an inverse modeling methodology which is based on zonation and is used to select an optimal zonation pattern based upon the weighted least squares performance criteria and covariance analysis. A coupled numerical model based on meshless local Petrov–Galerkin (MLPG) approach and particle swarm optimization (PSO) is developed to estimate the transmissivity in the various zones of an aquifer system. Here, 12 predesigned alternative models are considered and reliability analysis is performed to identify the best model structure for the assumed hypothetical confined aquifer. The sensitivity analysis of the estimated parameter is further investigated by considering different sets of head data involving measurement errors.

MESHLESS LOCAL PETROV–GALERKIN METHOD

The 2D groundwater flow governing equation for a heterogeneous, isotropic confined aquifer (Bear 1972) is considered as:  
formula
(1)
For transient state flow analysis, the initial condition used is h(x, y, 0) = ho(x, y) for Ω. The boundary conditions can be, Dirichlet boundary condition or Neumann boundary condition taken as h (x, y, t) = h1(x, y, t) for Ω1 and for Ω2, respectively.

Here, h(x, y, t) is hydraulic head (m); is transmissivity (m2/d); S is storage coefficient; x, y = horizontal space variables (m); Qw is pumping/injection rate at well (−Qw is pumping and +Qw is injection) (m3/d/m2); t is time in days; Ω is flow region; Ω is the boundary region (Ω1Ω2=Ω); is normal derivative; ho(x, y) is initial head in the flow domain (m); h1(x, y, t) is known head value of the boundary head (m); q(x, y, t) is known inflow rate (m3/d/m); δ is Dirac delta function, δ = 1 if x = xi, y = yi or δ = 0 if .

For groundwater flow simulation using meshless methods, we use a set of nodes within the problem domain as well as on the boundaries of the domain to represent (not discretize) the problem. As these nodes do not form a mesh they do not require prior information on the relationship between the nodes for the approximation of the unknown functions of the field variables (Swathi & Eldho 2014a). In the present study, a MLPG model based on MLS (moving least squares) is developed for inverse modeling in MATLAB. The shape functions are computed by applying the MLS approach as described in Swathi & Eldho (2017). To account for the time derivative, Crank–Nicolson time stepping scheme with weight factor, = 0.5, is used. The complete MLPG formulation for solving groundwater flow equation using a MLPG technique can be acquired from Swathi & Eldho (2014a; 2017). Based on the MLPG approach, Equation (1) is approximated and the resulting set of simultaneous equations are expressed in the matrix form as:  
formula
(2)
where  
formula
(3a)
 
formula
(3b)

The essential boundary conditions are enforced by the penalty parameter technique. Equation (2) is solved using the Gauss–Jordan method.

PARTICLE SWARM OPTIMIZATION

PSO is a swarm-based optimization which is based on the nature of flocking behavior of birds, fish, and other animals. This technique was first introduced by Kennedy & Eberhart (1995). The ‘particle’ in PSO is related to three main parameters: position, velocity, and fitness. Position represents the unknown variable of the problem, the rate of change of position (or change of variable) is given by velocity, and the fitness is a measure of how well the particle solves the objective function optimally (Kennedy & Eberhart 1995). The objective function is a measure of the quality of a solution, and in the present study, PSO is applied to get optimal solutions based on the minimization of the considered objective function.

The number of particles in a swarm is called population and ‘pbest’ is the coordinate in the search space associated with the best solution (fitness) the particle has achieved so far and ‘lbest’ is the best value obtained so far by any neighboring particle. When a particle takes all the population as its neighbors, the best value is a global best ‘gbest’. The PSO concept consists of, at each time step, changing the particle (velocity) towards its pbest and lbest locations. Acceleration is weighted by a random term, with separate random numbers being generated for acceleration toward pbest and lbest locations. The individual particle position is updated as follows (Kennedy & Eberhart 1995):  
formula
(4)
where velocity is calculated as:  
formula
(5)

where, is the particle position; is the particle velocity; is the best remembered individual particle position; is the best remembered swarm position; are the cognitive and social parameters; are the random numbers between 0 and 1 and w is the constriction factor.

MODEL DEVELOPMENT

Based on the MLPG-MLS methodology, a groundwater simulation model has been developed (Swathi & Eldho 2014a; 2014b; 2017) and verified. The MLPG-MLS model is coupled with the PSO-based optimization to get the MLPG-PSO model. The two-dimensional simulation-optimization (S-O) model is developed for steady and transient state conditions in MATLAB. The flowchart showing the methodology in the inverse modeling of aquifer parameters using the MLPG-PSO model is given in Figure 1.

Figure 1

Flowchart of the proposed inverse model.

Figure 1

Flowchart of the proposed inverse model.

MODEL APPLICATION

To investigate the performance of the proposed simulation-optimization model, a hypothetical confined aquifer of 36 sq km area (modified from Prasad & Rastogi (2001) and Snehalatha et al. (2006)) is considered. The aquifer domain (Figure 2) has a Dirichlet boundary condition of 100 m in the south, Neumann boundary condition in the west with an inflow rate of 0.2 m3/day/m and no flow boundaries in the east and north. On the northern part, the aquifer is recharged at the rate of 0.00015 m/d in AR1 and 0.00025 m/d in AR2, respectively, over the two strips through some aquitards. Two wells pumping at a rate of P1 with 600 m3/d and P2 with 1,000 m3/d are also considered within the flow domain. Eighteen observation wells are located within the flow region, locations of which are shown in Figure 3. The transmissivity values for the entire flow domain considered vary from 5 to 150 m2/day and storage coefficient is taken as 0.001. The considered flow problem is solved using MLPG-2D model with 169 nodes with Δx = Δy = 500 m (Figure 3) and radius of the sub-domain for the MLS approximation is 2Δ, where Δ is the distance between two consecutive nodes. The time step of 1 day is considered in the simulation. The hydraulic heads obtained from the MLPG model for the known transmissivity values for the flow problem are taken as observed heads in the present study.

Figure 2

Aquifer configuration with real transmissivity zonation (not to scale).

Figure 2

Aquifer configuration with real transmissivity zonation (not to scale).

Figure 3

Nodal arrangement for the problem with observation well location (shaded) and zonation (not to scale).

Figure 3

Nodal arrangement for the problem with observation well location (shaded) and zonation (not to scale).

The objective function (J) considered involves weighted least squares criterion, and the functional to be minimized (Swathi & Eldho 2014b) is:  
formula
(6)
subject to lower and upper bounds of the parameter: where is the MLPG head for the assumed parameter; is the MLPG head for the true parameter (these are replaced by observed head values in the field case studies); is transmissivity at block i; M is the parameter dimension; L is the number of observation wells; and are beginning and ending times of observations. Superscripts l and u are used to denote the lower and upper bounds of parameters and is the weighting coefficient = 1 for a parameter dimension. Also, the PSO model characteristics considered are: population size =100; maximum number of generations = 400; are 0.5 each;w is the constriction factor = 1.2.
For the present case study, 12 alternate models were constructed (Figure 4). The 12 models were decided based on trial and error such that the models vary from simplest homogeneous zone to complex heterogeneous zonal patterns. The aquifer parameters are estimated and the results are examined for their correctness and reliability. The error in parameter uncertainty is related to the variance-covariance matrix of the parameters. A well-estimated parameter is generally characterized by a small variance as compared to an insensitive parameter that is associated with large variance (Yeh & Yoon 1981). The diagonal elements of the variance-covariance matrix contain the individual parameter variances while the off-diagonal elements reflect the correlation between the parameters. The variance-covariance matrix of the estimated parameters is defined by Equation (7).  
formula
(7)
where is estimated parameters, is true parameters, E is mathematical expectation, and superscript T is a transpose vector. An approximation of the covariance matrix of the estimated parameters can be represented (Yeh & Yoon 1981) by the equation:  
formula
(8)
Figure 4

Zonation pattern for the hypothetical aquifer considered.

Figure 4

Zonation pattern for the hypothetical aquifer considered.

where is the least square error (LSE) (minimized value of objective function), L is the number of observations, M is parameter dimension, X is , and is the sensitivity matrix of h with respect to transmissivity T.

The objective function which is a measure of the modeling error in parameter estimation is computed. The trace of covariance matrix (TCM), a measure of uncertainty, is estimated as the sum of the diagonal components of covariance matrix. Based upon the LSE and TCM, the alternate models are compared (Beck & Arnold 1997; Kourakos & Mantoglou 2012). Parameter estimation is performed for each model for transient state at 100 days and is presented in Table 1.

Table 1

Estimated transmissivity for the 12 alternate models using transient state head data for t = 100 days and number of observation points: 18 (variance = 1): MLPG-PSO

Models
 
Zone 10 11 12 
150.2 150.2 159.2 143.6 149.4 149.9 150.2 148.3 172.4 173.2 180.1 111.2 
149.6 49.9 52.3 30.8 132.8 78.7 99.4 149.1 92.6 93.2  33.4 
50.4 15.2 14.7 4.9 29.9 49.8 33.4 73.2 30.1 27.0   
149.8 4.9   5.8 4.9 12.3 45.3     
50.5      4.2 26.5     
15.1       3.2     
50.2            
15.1            
4.9            
LSE 37 40 1,082 901 956 443 209 1,390 8,589 8,620 11,077 1,034 
TCM 730 22 535 208 3,586 582 721 8,994 4,032 4,002 189 823 
Runtime (min) 82 73 48 56 62 69 81 83 36 41 47 44 
Models
 
Zone 10 11 12 
150.2 150.2 159.2 143.6 149.4 149.9 150.2 148.3 172.4 173.2 180.1 111.2 
149.6 49.9 52.3 30.8 132.8 78.7 99.4 149.1 92.6 93.2  33.4 
50.4 15.2 14.7 4.9 29.9 49.8 33.4 73.2 30.1 27.0   
149.8 4.9   5.8 4.9 12.3 45.3     
50.5      4.2 26.5     
15.1       3.2     
50.2            
15.1            
4.9            
LSE 37 40 1,082 901 956 443 209 1,390 8,589 8,620 11,077 1,034 
TCM 730 22 535 208 3,586 582 721 8,994 4,032 4,002 189 823 
Runtime (min) 82 73 48 56 62 69 81 83 36 41 47 44 

Another inverse model with optimization based on the Levenberg–Marquardt algorithm (LMA) (Pujol 2007) with MLPG is developed in this study for the purpose of comparison. LMA is used to solve non-linear least squares problems and it interpolates between the Gauss–Newton algorithm and the method of gradient descent (Pujol 2007). Table 2 shows the parameter estimation results for inverse model based on MLPG-LMA. The simulations are carried out on a computer with Core i5 processor with speed of 3.1 Ghz and 8 Gb RAM and the computational runtime for different models are presented in Tables 1 and 2. From the results, as expected, the run time difference between the complex and simpler models is high. The time required to reach a best fit solution at all zones for complex zonal pattern models costs more computational time and vice versa.

Table 2

Estimated transmissivity for the 12 alternate models using transient state head data for t = 100 days and number of observation points: 18 (variance = 1): MLPG-LMA

Models
 
Zone 10 11 12 
149.3 149.1 152.8 140.2 149.8 148.9 149.9 148.2 175.7 173.2 178.3 109.9 
148.4 50.4 55.8 29.9 131.2 80.5 97.8 138.3 91.8 95.3  33.5 
51.3 15.4 12.3 4.8 27.3 40.4 32.8 75.3 26.5 26.1   
149.6 4.9   5.4 4.9 13.0 44.4     
51.1      3.9 26.2     
15.3       3.3     
50.2            
15.3            
4.9            
LSE 45 46 1,685 938 951 345 202 1,564 8,824 8,551 10,951 955 
TCM 1,087 37 544 251 4,529 665 624 9,511 4,127 4,458 126 841 
Runtime             
(min) 108 98 75 82 88 95 107 111 61 66 73 71 
Models
 
Zone 10 11 12 
149.3 149.1 152.8 140.2 149.8 148.9 149.9 148.2 175.7 173.2 178.3 109.9 
148.4 50.4 55.8 29.9 131.2 80.5 97.8 138.3 91.8 95.3  33.5 
51.3 15.4 12.3 4.8 27.3 40.4 32.8 75.3 26.5 26.1   
149.6 4.9   5.4 4.9 13.0 44.4     
51.1      3.9 26.2     
15.3       3.3     
50.2            
15.3            
4.9            
LSE 45 46 1,685 938 951 345 202 1,564 8,824 8,551 10,951 955 
TCM 1,087 37 544 251 4,529 665 624 9,511 4,127 4,458 126 841 
Runtime             
(min) 108 98 75 82 88 95 107 111 61 66 73 71 

Based on the results, models 1 and 2 are considered superior to all other models due to many less objective function values. In general, where both LSE and TCM values are lower, the model is more preferable. Models 1 and 2 represent the zonation pattern correctly although the latter consists of fewer zones. These two models are considered for reliability analysis to select the best model. However, one more model, say model 4, which does not represent the true zonation pattern, is also considered for further analysis to assess the inverse model performance.

RELIABILITY OF ESTIMATED PARAMETERS

The reliability of estimated parameters is ascertained by computing three parameters, composite scaled sensitivity (CSS), coefficient of variation (CV) and the confidence interval (CI), which are subsequently explained. The above parameters are estimated for models 1, 2, and 4, which are presented in Tables 3 and 4.

Table 3

Estimated transmissivities and their reliability for different models using MLPG-PSO

  Model 1
 
Model 2
 
Model 4
 
Zone True T (m2/d) Estimated T (m2/d) CV CI 95% CSS True T (m2/d) Estimated T (m2/d) CV CI 95% CSS Estimated T (m2/d) CV CI 95% CSS 
150.0 150.2 0.049 ±7.2 15.23 150.0 150.2 0.025 ±7.1 28.33 143.6 0.532 ±74.7 4.03 
150.0 149.6 0.262 ±21.4 8.65 50.0 49.9 0.020 ±5.1 28.08 30.8 0.486 ±63.8 12.64 
50.0 50.4 0.098 ±22.8 9.32 15.0 15.2 0.032 ±7.3 18.72 4.9 0.019 ±3.5 26.45 
150.0 149.8 0.034 ±5.9 21.82 5.0 4.9 0.019 ±3.5 26.77     
50.0 50.5 0.262 ±31.4 7.72          
15.0 15.1 0.012 ±3.5 16.64          
50.0 50.2 0.034 ±7.6 25.48          
15.0 15.1 0.018 ±3.8 26.08          
5.0 4.9 0.019 ±3.5 26.36          
Average % error: 3.076 2.171 43.187 
  Model 1
 
Model 2
 
Model 4
 
Zone True T (m2/d) Estimated T (m2/d) CV CI 95% CSS True T (m2/d) Estimated T (m2/d) CV CI 95% CSS Estimated T (m2/d) CV CI 95% CSS 
150.0 150.2 0.049 ±7.2 15.23 150.0 150.2 0.025 ±7.1 28.33 143.6 0.532 ±74.7 4.03 
150.0 149.6 0.262 ±21.4 8.65 50.0 49.9 0.020 ±5.1 28.08 30.8 0.486 ±63.8 12.64 
50.0 50.4 0.098 ±22.8 9.32 15.0 15.2 0.032 ±7.3 18.72 4.9 0.019 ±3.5 26.45 
150.0 149.8 0.034 ±5.9 21.82 5.0 4.9 0.019 ±3.5 26.77     
50.0 50.5 0.262 ±31.4 7.72          
15.0 15.1 0.012 ±3.5 16.64          
50.0 50.2 0.034 ±7.6 25.48          
15.0 15.1 0.018 ±3.8 26.08          
5.0 4.9 0.019 ±3.5 26.36          
Average % error: 3.076 2.171 43.187 
Table 4

Estimated transmissivities and their reliability for different models using MLPG-LMA

  Model 1
 
Model 2
 
Model 4
 
Zone True T (m2/d) Estimated T (m2/d) CV CI 95% CSS True T (m2/d) Estimated T (m2/d) CV CI 95% CSS Estimated T (m2/d) CV CI 95% CSS 
150.0 149.3 0.132 ±9.7 12.89 150.0 149.1 0.146 ±18.4 13.35 140.2 0.501 ±70.2 4.52 
150.0 148.4 0.308 ±36.1 8.87 50.0 50.4 0.042 ±11.3 16.78 29.9 0.473 ±58.6 11.37 
50.0 51.3 0.272 ±34.7 6.63 15.0 15.4 0.021 ±8.9 21.26 4.8 0.028 ±5.1 18.34 
150.0 149.6 0.084 ±14.7 15.71 5.0 4.9 0.019 ±3.4 26.95     
50.0 51.1 0.240 ±29.4 7.83          
15.0 15.3 0.056 ±9.1 12.54          
50.0 50.2 0.018 ±7.3 16.32          
15.0 15.3 0.047 ±8.8 17.45          
5.0 4.9 0.018 ±3.5 26.63          
Average % error: 4.834 3.321 44.842 
  Model 1
 
Model 2
 
Model 4
 
Zone True T (m2/d) Estimated T (m2/d) CV CI 95% CSS True T (m2/d) Estimated T (m2/d) CV CI 95% CSS Estimated T (m2/d) CV CI 95% CSS 
150.0 149.3 0.132 ±9.7 12.89 150.0 149.1 0.146 ±18.4 13.35 140.2 0.501 ±70.2 4.52 
150.0 148.4 0.308 ±36.1 8.87 50.0 50.4 0.042 ±11.3 16.78 29.9 0.473 ±58.6 11.37 
50.0 51.3 0.272 ±34.7 6.63 15.0 15.4 0.021 ±8.9 21.26 4.8 0.028 ±5.1 18.34 
150.0 149.6 0.084 ±14.7 15.71 5.0 4.9 0.019 ±3.4 26.95     
50.0 51.1 0.240 ±29.4 7.83          
15.0 15.3 0.056 ±9.1 12.54          
50.0 50.2 0.018 ±7.3 16.32          
15.0 15.3 0.047 ±8.8 17.45          
5.0 4.9 0.018 ±3.5 26.63          
Average % error: 4.834 3.321 44.842 

Coefficient of variation

The coefficient of variation (CV) is equal to the ratio of the standard deviation of estimated parameter and parameter estimate, which is dimensionless. A relatively smaller value of CV corresponds to better estimate of the zonal parameter.  
formula
(9)

From Tables 3 and 4, as expected, model 4 gives higher values and thus this model can be discarded. The CV values for model 1 (Table 3) in zones 2 and 5 are higher. This suggest that the reliability of the estimated parameters in these zones is less. From Table 3, it is also found that the CV values for model 2 are far less compared to model 1 which suggest its better reliability. Also, PSO algorithm is able to estimate parameters better than LMA algorithm.

Confidence intervals

The confidence with which the parameters are estimated can be determined by computing the CI based on approximating the first order methods. The individual CIs presented in this study are calculated as the final estimated parameter value ± the product of the standard deviation of that estimate and Student's t-statistic. The Student's t-statistic is dependent upon the degrees of freedom of the problem (total number of observations minus the number of estimated parameters) and the desired level of confidence. Finally, the CI is a function of the quantity of observed data, parameterization and the uncertainty of observed data to estimated parameters. In the present analysis, the 95% CI values are used.

The results from Tables 3 and 4 show that model 4 is clearly not a good fit to represent the aquifer domain. Then, when model 1 (Table 3) is considered, the CI values for the transmissivity estimates of zones 2, 3, and 5 are higher indicating that the parameter uncertainty is greater. The results of model 2 show that the CI is very much less. This suggests that higher CI for some zones in model 1 may be due to insufficient observed head data whereas re-parameterization (model 2) has improved the reliability of estimated parameters. The higher CI (model 1) shows that the uncertainty of head observations with the estimated parameters in some of the zones is very low, hence such type of zones can be combined and thereby the number of zones can be reduced. The above results show that model 2 represents the aquifer zonation pattern more accurately. Model 1 represents over-parameterization, which can be avoided by zonation-based studies. The CI results from Tables 3 and 4 prove that the meta-heuristic PSO algorithm is able to estimate parameters better than the LMA algorithm.

Composite scaled sensitivity

CSS is a statistic which summarizes all the sensitivities for one parameter and therefore indicates the cumulative amount of information that the measurements contain for the estimation of that parameter. CSS (Hill et al. 1998) for parameter j, is computed as:  
formula
(10)
where is the estimated transmissivity, is the weighting coefficient, is the sensitivity coefficient, and L is the number of observations. Parameters with relatively large CSS values are likely to be better estimated. Between models 1 and 2, the CSS values in the case of model 2 are higher compared to model 1. It is observed that the CSS values of the estimated transmissivity parameter in zones 2, 3, and 5 are less compared to CSS values for the remaining zones (Tables 3 and 4). This shows that sensitivity in the estimated transmissivity parameter in zones 2, 3, and 5 is more in the case of model 1. This sensitivity is reflected in terms of higher percentage error in estimated transmissivity of that zone in model 1. Similarly, it is found that the CSS values of estimated transmissivity of zones 7, 8, and 9 are higher compared to CSS values for the remaining zones, which shows that the reliability of the estimated transmissivity parameter in these zones is higher. Thus, based on CSS results, model 2 is better representing the problem domain. For model 4, the CSS values are lower indicating it as a poorer zonation model and has to be restructured for improvement. Based on the above reliability analysis, model 2 is considered as the best model and is taken for further sensitivity studies.

Sensitivity analysis

Further, for sensitivity analysis for model 2, different data sets of observed head are obtained by adding certain noise to the MLPG solutions for both steady and transient cases. Data set 1 represents the error-free head data, whereas data sets 2 to 4 are obtained by adding normally distributed random errors with zero mean and standard deviations of 0.01, 0.1, and 1, respectively, to the data set 1. This analysis is carried out for both steady and transient state to establish the influence of the quality of head measurement in estimation of parameters using the MLPG-PSO model. The data sets with noise are considered to be equivalent to the observation heads with measurement errors, which are common in field cases.

The observed head data at 100 days of pumping are obtained from the MLPG model at the 18 observation points. From Tables 5 and 6, the results from the model show that the parameter values are estimated satisfactorily for all four data sets within ±5 units of true estimates which is acceptable in groundwater systems. This analysis demonstrates the applicability of the present methodology for real field studies.

Table 5

Sensitivity analysis for estimation of transmissivity in various zones for steady state observed head data sets and measurement errors for 18 observation points

Zone Estimated parameter (T)
 
True parameter (m2/d) 
Data set 1 (Error free) Data set 2 N (0, 0.01) Data set 3 N (0, 0.1) Data set 4 N (0, 1) 
150.2 150.42 149.88 148.70 150 
49.9 50.12 51.09 53.82 50 
15.2 15.32 15.68 16.01 15 
4.9 5.03 5.24 6.89 
Zone Estimated parameter (T)
 
True parameter (m2/d) 
Data set 1 (Error free) Data set 2 N (0, 0.01) Data set 3 N (0, 0.1) Data set 4 N (0, 1) 
150.2 150.42 149.88 148.70 150 
49.9 50.12 51.09 53.82 50 
15.2 15.32 15.68 16.01 15 
4.9 5.03 5.24 6.89 
Table 6

Sensitivity analysis for estimation of transmissivity in various zones for transient state observed head data sets and measurement errors for 18 observation points at t = 100 days

Zone Estimated parameter (T)
 
True parameter (m2/d) 
Data set 1 (Error free) Data set 2 N (0, 0.01) Data set 3 N (0, 0.1) Data set 4 N (0, 1) 
150.2 150.40 151.12 152.84 150 
49.9 50.15 51.24 52.99 50 
15.2 15.34 16.08 16.86 15 
4.9 5.08 5.13 6.48 
Zone Estimated parameter (T)
 
True parameter (m2/d) 
Data set 1 (Error free) Data set 2 N (0, 0.01) Data set 3 N (0, 0.1) Data set 4 N (0, 1) 
150.2 150.40 151.12 152.84 150 
49.9 50.15 51.24 52.99 50 
15.2 15.34 16.08 16.86 15 
4.9 5.08 5.13 6.48 

DISCUSSION

The MLS-based MLPG model has been proven to be an effective alternative approach in simulating groundwater flow. The MLPG is a truly meshless formulation which significantly reduces the manual efforts of meshing and re-meshing leading to lesser computational efforts, which was a major disadvantage posed by grid-based methods such as finite difference (FDM) and finite element (FEM) methods. Other advantages of the MLPG model include: no global background cell for integration, local form formulation leading to simple numerical implementation, no need of global compatibility of trail (shape) function, and better rate of convergence and accuracy. Further, for inverse modeling and parameter estimation, the PSO-based optimization model is considered. PSO-based optimization model has a number of advantages such as achieving global optimal solution in few iterations, simplicity, lower number of parameters, ease of implementation in comparison with other evolutionary algorithms, and less computational efforts. For the PSO model developed, a sensitivity analysis of parameters is carried out before its application and thus faster convergence is achieved.

In groundwater hydrology, generally the hydraulic conductivity or transmissivity is considered zonal wise as continuous variation rarely occurs and it is difficult to estimate. In this study, for the problem considered, the most suitable model structure is identified from a number of alternate models using the zonation method of parameter estimation. A range of 12 alternate models were identified based on trial and error, starting from homogeneous to more complex model structures, and considered for the zonation. The inverse modeling of different model structures is carried out based on the LSE and TCM for each model. Based on this, the most suitable models are selected and reliability analysis ascertained by computing the three parameters of CSS, CV, and CI, with the best model being selected.

However, in real-life scenarios where prior knowledge about an aquifer is sparse, the general approach to obtain the zonation pattern is to couple the zonation algorithm with a clustering approach, such as using hierarchical or fuzzy K-means, where the optimum parameter structure associated with the groundwater inverse problem is determined. This step, wherein an optimization technique which can be used to determine optimum zone structure where initially the number of zones selected are one and are increased over iterations until the optimum parameter structure is identified, is not attempted here, as this was the first attempt made to solve inverse modeling using the MLPG-PSO model. This can be attempted in future studies.

CONCLUSIONS

In this study, the application of a simulation optimization model based on MLPG method and PSO for aquifer parameter identification and zonation structure estimation is presented. For the considered case study, a number of possible alternate model structures using zonation method are identified through a trial and error approach. In the inverse modeling procedure, weighted least square performance criterion is considered for the estimation of transmissivity for the confined aquifer considered. Reliability analysis is performed to identify the best model structure. Models 1 and 2 are identified as better model structures out of the 12 considered models. However, based on reliability of estimated parameters, it can be concluded that over-parameterization occurs in the case of model 1 which leads to more uncertainty in the estimated parameters. Hence, model 2 is considered to be the best model for the studied aquifer system.

From the present study it can be concluded that estimation of parameters needs to be supported by reliability analysis for the appropriate inverse modeling of an aquifer system. To demonstrate the effectiveness of the MLPG-PSO model, the results are compared with another optimization technique, the LMA. Comparison of model results shows that the MLPG-PSO results are better, proving the superiority of meta-heuristic algorithms. Further, the transmissivity values are computed for both steady and transient states using error free data and different sets of noisy data. This indicates that the model can be applied to regional aquifers where the data available may be sparse and with errors. The present study demonstrates the effectiveness of the proposed methodology for aquifer parameter and zonation structure estimation in groundwater systems.

ACKNOWLEDGEMENTS

The authors are grateful to the anonymous reviewers and the editor for their constructive comments and their suggestions which helped to improve the final presentation of this paper.

REFERENCES

REFERENCES
Anderman
,
E. R.
,
Hill
,
M. C.
&
Poeter
,
E. P.
1996
Two-dimensional advective transport in ground-water flow parameter estimation
.
Ground Water
34
(
6
),
1001
1009
.
Bear
,
J.
1972
Hydraulics of Groundwater
.
McGraw Hill Publishing
,
New York
.
Beck
,
J. V.
&
Arnold
,
K. J.
1997
Parameter Estimation in Engineering and Science
.
Wiley
,
New York
.
Carrera
,
J.
,
Alcolea
,
A.
,
Medina
,
A.
,
Hidalgo
,
J.
&
Slooten
,
L. J.
2005
Inverse problem in hydrology
.
Hydrogeol. J.
13
,
206
222
.
Hendricks Franssen
,
H. J.
,
Alcolea
,
A.
,
Riva
,
M.
,
Bakr
,
M.
,
van der Wiel
,
N.
,
Stauffer
,
F.
&
Guadagnini
,
A.
2009
A comparison of seven methods for the inverse modeling of groundwater flow: application to the characterization of well catchments
.
Adv. Water Resour.
32
,
851
872
.
Hill
,
M. C.
&
Tiedeman
,
C. R.
2007
Groundwater Model Calibration
.
Wiley
,
New Jersey
.
Hill
,
M. C.
,
Cooley
,
R. L.
&
Pollock
,
D. W.
1998
A controlled experiment in groundwater flow model calibration
.
Ground Water
36
(
3
),
520
535
.
Kennedy
,
J.
&
Eberhart
,
R.
1995
Particle swarm optimization
. In:
Proceedings of the IEEE International Conference on Neural Networks
.
IEEE Service Center
,
Piscataway, NJ
, pp.
1942
1948
.
Kourakos
,
G.
&
Mantoglou
,
A.
2012
Inverse groundwater modeling with emphasis on model parameterization
.
Water Resour. Res.
48
(
5
),
W05540
.
Maier
,
H. R.
,
Kapelan
,
Z.
,
Kasprzyk
,
J.
,
Kollat
,
J.
,
Matott
,
L. S.
,
Cunha
,
M. C.
,
Dandy
,
G. C.
,
Gibbs
,
M. S.
,
Keedwell
,
E.
,
Marchi
,
A.
,
Ostfeld
,
A.
,
Savic
,
D.
,
Solomatine
,
D. P.
,
Vrugt
,
J. A.
,
Zecchin
,
A. C.
,
Minsker
,
B. S.
,
Barbour
,
E. J.
,
Kuczera
,
G.
,
Pasha
,
F.
,
Castelletti
,
A.
,
Giuliani
,
M.
&
Reed
,
P. M.
2014
Evolutionary algorithms and other metaheuristics in water resources: current status, research challenges and future directions
.
Environ. Model. Softw.
62
,
271
299
.
Rajanayaka
,
C.
&
Kulasiri
,
D.
2001
Investigation of a parameter estimation method for contaminant transport in aquifers
.
J. Hydroinformatics
3
,
203
213
.
Sun
,
N.-Z.
1994
Inverse Problems in Groundwater Modeling
.
Kluwer Academic
,
Dordrecht
.
Sun
,
N.-Z.
&
Yeh
,
W. W.-G.
1985
Identification of parameter structure in groundwater inverse problem
.
Water Resour. Res.
21
(
6
),
869
883
.
Swathi
,
B.
&
Eldho
,
T. I.
2014b
Inverse modeling of groundwater system using coupled PSO-MLPG techniques
. In:
Proceedings of the 11th International Conference on Hydroinformatics
,
CUNY Academic Works
,
New York
.
Yeh
,
W. W.-G.
&
Yoon
,
Y. S.
1981
Parameter identification with optimum dimension in parameterization
.
Water Resour. Res.
17
(
3
),
664
672
.
Yeh
,
T. J.
,
Mao
,
D.
,
Zha
,
Y.
,
Wen
,
J.
,
Wan
,
L.
,
Hsu
,
K.
&
Lee
,
C.
2015
Uniqueness, scale, and resolution issues in groundwater model parameter identification
.
Water Sci. Eng.
8
,
175
194
.