Aquifer parameter and zonation structure estimation using meshless local Petrov–Galerkin method and particle swarm optimization

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.

Here, h(x, y, t) is hydraulic head (m); is transmissivity (m²/d); S is storage coefficient; x, y = horizontal space variables (m); Q_w is pumping/injection rate at well (−Q_w is pumping and +Q_w is injection) (m³/d/m²); t is time in days; Ω is flow region; Ω is the boundary region (⁠Ω₁Ω₂=Ω); is normal derivative; h_o(x, y) is initial head in the flow domain (m); h₁(x, y, t) is known head value of the boundary head (m); q(x, y, t) is known inflow rate (m³/d/m); δ is Dirac delta function, δ = 1 if x = x_i, y = y_i or δ = 0 if ⁠.

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

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.

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.

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.

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 m³/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 m³/d and P2 with 1,000 m³/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 m²/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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.