In recent decades, due to the population growth and low precipitation, the overexploitation of ground water resources has become an important issue. To ensure a sustainable scheme for these resources, understanding the behavior of the aquifers is a key step. This study takes a numerical modeling approach to investigate the behavior of an unconfined aquifer in an arid area located in the east of Iran. A novel hybrid model is proposed that couples the numerical modeling to a data assimilation model to remove the uncertainty in the hydrodynamic parameters of the aquifer including the hydraulic conductivity coefficients and specific yields. The uncertainty that exists in these parameters results in unreliability of the head values acquired from the models. Meshless local Petrov-Galerkin (MLPG) is used as the numerical model, and particle filter (PF) is our data assimilation model. These models are implemented in the MATLAB software. We have calibrated and validated our PF-MLPG model by the observation head data from the piezometers. The RMSE in head values for our model and other commonly used numerical models in the literature including the finite difference method and MPLG are calculated as 0.166, 1.197 and 0.757 m, respectively. This fact shows the necessity of using this method in each aquifer.

  • The PF-MLPG method is compatible with different kinds of aquifers and different conditions, including confined, unconfined, transient, steady, irregular geometry and regular geometry.

  • PF-MLPG simulates groundwater flow with satisfactory accuracy that helps decision makers to manage the aquifers properly.

  • This new calibration–simulation model removes high uncertainty of model and input parameters such as hydraulic conductivity coefficient and specific yield.

Aquifers have always been important water supply resources in the world, especially in semi-arid and arid regions with a low amount of rainfall and the scarcity of surface water resources. In recent decades, due to the population growth, the need for water increases dramatically and low precipitation, the aquifers have been seriously threatened by the overexploitation of groundwater. However, good planning and management as the result of good understanding of the aquifers' behavior can help resolving this issue and insuring their sustainable and optimal use for us and our next generations (Ghousehei 2010).

Aquifers are directly or indirectly fed by surface water and rainfall, and a good management of the supply–demand system can help their sustainable and optimal use. To this end, the aquifer behavior requires studies and discussions for which water engineers/managers need to have accurate information.

Undoubtedly, the best approach to understand the behavior of an aquifer system is to conduct a set of long-term field studies for each specific area. However, it is practically impossible due to limitations in research budget and time as it is time-consuming. Numerical models and simulations have shown to create satisfactory results in natural phenomena in a limited time. Especially, with the evident of super computers with high computational capabilities, using these numerical models has become very popular.

Aquifer models simulate the groundwater flow by solving the related governing equations using analytical methods (Akbarpour et al. 2012), but this is limited to very simple geometry cases and generally educational applications. When the geometry is irregular and complex, they are either non-existent or of little use; they are used only to evaluate and compare the numerical results to control and validate the method. To this end, they use numerical methods to solve the aquifers' governed equation. The finite elements, finite difference, boundary element and finite volume are among important numerical methods used for this purpose. Recently, meshless methods, although only three decades old, have also been added to the mentioned methods and have provided a wide and suitable platform for scientific, research and engineering activities.

Unlike mesh-dependent methods (e.g. finite element (FE)), meshless methods do not need to define a specific inter-node relation to interpret the problem's physical behavior; elements are replaced with a set of nodes and the problem's physical domain is, in fact, defined with enough number of points. In FE methods, points are related by elements, but in meshless methods the relation is through the number of common points in the common area of each sub-domain (Markopoulos et al. 2020). Researchers' great interest in these methods is due to the highly reduced mesh creation time in FE methods. Although meshless methods are not yet as widespread as finite element methods in engineering problems, they have been used here because of their efficiency and some other merits (Mohtashami et al. 2017a).

A number of authors, including those of this paper, have used meshless methods to model aquifers. In this regard, Swathi & Eldho (2017) used Petrov-Galerkin's meshless method that used the radial basis function (RBF) instead of the moving least squares (MLS) function, which had such limitations as changes in the aquifer boundary. They modeled two standard rectangular aquifers and a real one, found the groundwater level and compared their results with those of the meshless point collocation, finite difference and finite element methods; the main result of their research was the higher accuracy of Petrov-Galerkin's meshless method (Swathi & Eldho 2017). Mohtashami et al. (2017b) used the Petrov-Galerkin's meshless method to model aquifers in permanent and non-permanent states, calibrated and validated it and showed that the method was more accurate than the finite difference method. The important point in their work was using the MLS function that besides being used for the first time in this field, solved the problem of creating dummy nodes, which was the main drawback of the RBF (Mohtashami et al. 2017b). Pointing out the meshless models' superiority over mesh-dependent ones, Pathania et al. (2019) modeled three aquifers by the element-free Galerkin that was based on the MLS approximation function. Aquifers wherein the groundwater level was modeled were one 1D and two 2D aquifers; their real 2D aquifer was the Blue Lake in Canada (Pathania et al. 2019).

The uncertainty analysis of the aquifers' parameters is a serious issue beside the use of their simulation models. Uncertainty, a world's inherent characteristic widely seen in the man's daily activities usually stems from a lack of accuracy in the present or future situation due to limited data and information. Aquifer is a complex and open system that affects the hydrological and meteorological conditions, geological structure, topographic features, vegetation, human activities and so on. However, some features of this system are not directly measured; for instance, the hydraulic conductivity is found by analyzing the measured inlet and outlet values. Hence, studying uncertain parameters in the aquifer modeling is of great importance.

Among extensive studies on the uncertainty analysis, the share of those on groundwater models is not much (Abedini et al. 2017). Pohll et al. (2002) and Dettinger & Wilson (1981) studied uncertainty in aquifer simulations with a statistical view and proposed it to be used along with the Monte Carlo method as a useful tool in the analysis of uncertainty in future simulations.

Algorithms used in the uncertainty analysis of groundwater models include Monte Carlo, Markov chain and so on, but very limited studies have investigated uncertainties in groundwater modeling, most of which have been on the sensitivity analysis of model parameters and analysis of statistical indices. For instance, Chitsazan et al. (2008) used the statistical characteristics method for the uncertainty analysis of the mathematical model of Kazerun plain aquifer and observed a significant improvement in the model fitting criteria (Chitsazan et al. 2008). In the same year, Rasoulzadeh & Mousavi (2008) used the Levenburg algorithm (an inverse method) to study the uncertainty in estimating the aquifer model parameters of Arsanjan plain and concluded that although the simulated and observed groundwater levels conformed well in the calibration phase, the estimated parameters were not unique and had high uncertainty because some of them were highly correlated. GLUE1 is another method used to analyze uncertainties in aquifer models; Hamraz et al. (2015) used it along with the finite difference numerical model in the groundwater model system (GMS) environment to model the groundwater flow in an aquifer and evaluated the uncertainty of hydraulic conductivity and specific yield (Hamraz et al. 2015). Using the finite difference method in GMS and the GLUE method, Abedini et al. (2017) studied and evaluated the hydraulic conductivity and recharging parameters in the Bojnourd aquifer and showed that in the calibration process, recharging was more uncertain and less detectable than other parameters. Hydraulic conductivity too was close to its exact values; error was effectively reduced. Swathi & Eldho (2017) developed a model for estimating groundwater parameters. They applied a simulation–optimization model in a standard aquifer in order to identify the zonation structure and also groundwater parameters. Their achieved results showed satisfactory values (Swathi & Eldho 2018).

Using the Modflow model in GMS and such evaluation criteria as RMSE (root mean square error) and NRMSE (normalized root mean square error), fluctuations of groundwater level, discharge from aqueducts, high-discharge sites and status of the groundwater balance, Du et al. (2018) modeled and studied the reliability of the hydrodynamic parameters of an aquifer, China, and finally improved their modeling accuracy (Du et al. 2018). For a sustainable groundwater management in an aquifer in northwestern Bangladesh, Thouhidul Mostafa et al. (2019) predicted its groundwater level and showed the uncertainty effects of conceptual hydrogeological/climatic models in their predictions. Using a combination of GMS and Bayesian model, they found that 23% of the existing uncertainties depended on conceptual hydrogeological models and the rest were related to the flow and climatic models (Thouhidul Mostafa et al. 2019). Thomas et al. (2019) engaged some models to estimate parameters, e.g. transmissivity coefficient and longitudinal and transverse dispersivity coefficients. Their model included genetic algorithm, cat swarm optimization, particle swarm optimization and differential evolution. Results showed good performance of their model (Thomas et al. 2019). Using PMWIN (a finite difference-based numerical method) to study the input parameter uncertainties of their model, Nossent et al. (2020) simulated the groundwater level in an aquifer in Bangladesh and studied the uncertainties of hydraulic conductivity and specific yield parameters by combining the Differential Evolution Adaptive Metropolis (DREAM) method and the Bayesian model and using a new likelihood criterion. Results showed that their uncertainty model was highly efficient and enabled them to bring the modeling results closer to the observational data (Nossent et al. 2020). Eryiğit (2021) employed a modified clonal selection algorithm to estimate groundwater flow parameters. MODLFOW was engaged as simulation as well. The proposed model was applied in two hypothetical cases. Results presented the feasibility of the model in the estimation of groundwater parameters (Eryiğit 2021).

In this study, a new online calibration particle filter (PF) method has been coupled with a meshless groundwater flow model to find the acceptable hydrodynamic parameters of a real aquifer including hydraulic conductivity coefficient and specific yield. As the PF is a powerful method of estimating uncertain parameters, combining with a precise flow model is practical and predicts, highly accurately, how water moves in the aquifer. It is worth noting that it is the first time this method is used in a real field aquifer (Birjand aquifer) as it is challenging and since it is more accurate than old, traditional models, it examines the groundwater behavior with high accuracy and can also be used along with meshless aquifer models and be combined with the moving kriging (MK) approximation method to reduce the flow model-caused errors.

Case study (real aquifer)

Birjand unconfined aquifer is chosen as a real aquifer to investigate the performance of the proposed model in a challenging condition like irregular boundary, uncertainty of hydrodynamic parameters, etc. This aquifer is located in south Khorasan province, east of Iran. It is surrounded by the latitude of 32°34′ to 33°8′ North and longitude of 58°41′ to 59°44′ East. The area of this aquifer is 269 km2. Based on the Domarton climate index, this aquifer is counted as an arid region. The annual precipitation of this region is 170 mm (Mohtashami et al. 2017b). Figure 1 shows the location of the aquifer in Iran.

Figure 1

Geographical location of the Birjand aquifer in Iran. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.10.2166/hydro.2021.113.

Figure 1

Geographical location of the Birjand aquifer in Iran. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.10.2166/hydro.2021.113.

Close modal

There are 190 extraction wells and 10 observation wells (piezometers) in this aquifer which is determined with blue and red circle symbols in Figure 1(c), respectively.

Particle filter

A PF method known as a powerful estimation tool that computes the probability density function (PDF) of a random process and also estimates the state of the object in the future time precisely based on the states and observations of previous times. The PF makes some estimates for the state of the object to select the best one (Ruknudeen & Asokan 2017).

To this end, in the initial step, particles are scattered in the space which its dimension equals to the number of parameters that must be estimated. The particles are state vectors that are aggregated and use the probability distribution function. Each particle is assigned with a weight value. This weight value in the first step is equal for all particles around () (N is the number of the particles scattered in state space) (Zhang et al. 2010). In the next step, the weight value depending on the position of the particle is updated. Once the weight of each particle is determined, to prevent from degeneracy occurrence, a resampling method is carried out.

To be more specific, the particles with very low weight are emitted, and the particles with high weight are involved in the generation of new particles around their coordinates. This function reduces computation time and prevents the particles from spreading in vain (Ramgraber et al. 2019).

The particle filtering is an effective way to solve the estimation issues in nonlinear Gaussian noise systems. This method is a part of the Monte Carlo statistical methods. In order to describe the PF, consider the following nonlinear system:
(1)
(2)
where is the state vector, is the measurement vector, and f and g are nonlinear functions, respectively. Furthermore, it is assumed that and are processes called noise and measurement noise, respectively. The PF of the posterior PDF is defined as a set of weighted particles as follows (Figure 2):
(3)
where N, and are the number, locations and weights of particles. The PF approximates the posterior function of as follows (Manoli et al. 2015; Havangi 2018):
(4)
where is the Dirac Delta function and is the weight of and .
Figure 2

Descriptions of the posterior function (Havangi 2018).

Figure 2

Descriptions of the posterior function (Havangi 2018).

Close modal
Direct sampling from the PDF (known as the objective weighing function) is not possible; hence, the important sampling method is used. In the important sampling method, instead of directly sampling the target function, a proposed weighing function is proposed. Considering the proposed density function as , the particle weight is calculated as follows (Manoli et al. 2015):
(5)
By decomposing the proposed density function as follows, weights will be calculated recursively (Choe et al. 2015):
(6)
In this case, the particles of that are a sample from the proposed density distribution function of are obtained by the new state of that is sampled from . The weight of these particles is as follows (Choe et al. 2015; Li et al. 2015):
(7)
On the other hand, according to the Bayes law, the posterior PDF can be written as follows (Havangi 2018):
(8)
So, the particle weight is calculated as the following recursive relation (Arulampalam et al. 2002):
(9)

This method is known as the post probability of significant sequential important sampling (SIS) (Arulampalam et al. 2002; Havangi 2018). In order to avoid the corrosion phenomenon, with resampling, the SIS algorithm is introduced (Ruknudeen & Asokan 2017). This algorithm which is known as the sequential important resampling algorithm has three steps including sampling, calculating particle weight and resampling:

  1. Sampling: The particles are obtained by sampling from the proposed distribution.

  2. (10)
  3. Calculating particle weight with using Equation (9).

  4. Resampling: Higher-weight particles replace the lower-weight particles.

As we mentioned, particle degeneracy occurs in PF when the variance of weight values becomes large. This is the main drawback of PF in the computation procedure. In order to solve this problem, resampling method recommended. In this technique, the following concept introduced in the following equation (Havangi 2016):
(11)
In this equation, when the weight value () is lower than a threshold size (), resampling algorithm activates. The important issue is the value of (). Based on the references presented in this field, the best value for this parameter is (Manoli et al. 2015). The following equation shows the condition that resampling technique must be carried out (Havangi 2018):
(12)

A resampling technique makes the PF method compatible with the condition of problem and removes the particle degeneracy. In this paper, the particles considered the state space amount to 100, 80 and 60. The different numbers of particles considered in the PF-meshless local Petrov-Galerkin (MLPG) model and sensitivity analysis are carried out.

Meshless Local Petrov-Galerkin

Extensive research has been carried out in fluid mechanics related to the meshless methods (Liu & Gu 2005). These methods remove the drawbacks of the grid-based numerical methods, e.g. FEM and FDM, by just adding or removing a set of field nodes (Akbari et al. 2010). They overcome the limitations of the FEM numerical method in updating boundary conditions (Atluri & Zhu 1998). These methods have outperformed FEM in terms of accuracy (Porfiri 2006). Also, the defects of the finite difference method, which is restricted to the rectangular meshes, are eliminated in meshless numerical methods.

MLPG is a weak form of meshless methods presented in 1998 by Atluri & Zhu (1998) to solve the potential equation. This method is generally used in fluid mechanics like modeling the free surface, Laplace equation, etc., and includes two functions: the weight function (cubic spline) and MK which is used to construct a shape function. MLPG is based on the MLS approximation function. However, it suffers from satisfying Kronecker delta function. This property prevents MLPG to directly impose essential boundary conditions. In order to overcome this defect, MK function is applied. This function satisfies the Kronecker delta function and provides the model to enforce essential boundary conditions directly. A comprehensive explanation can be found in Mohtashami et al. (2017a, 2017b). It should be mentioned that MLPG has already been used by authors to simulate groundwater flow in field aquifer, and it is calibrated and verified by the observation data and FDM solutions. The higher accuracy of MLPG indicates its power and capability.

MK method

MK method formulation is similar to MLS2 approximation function which is used by MLPG in order to construct shape function (Khankham et al. 2015). This function satisfies the Kronecker delta property and considering the function is defined as groundwater head. Therefore, the formulation of shape function is presented in the following equation:
(13)
where is the shape function and L is the number of nodes scattered in the whole domain. In this equation, is computed with the following equation:
(14)
where A and B are expressed as follows:
(15)
(16)
In these equations, I is a unit matrix and is the basis function obtained from the monomials in the Khayyam–Pascal triangle. R and are defined in the following equation:
(17)
(18)
where is the correlation function between any pair of nodes located at and .

Groundwater flow equation in the unconfined aquifer

Groundwater equation in unconfined aquifer in is represented in the following equation (Wang & Anderson 1995):
(19)

In this equation, and are hydraulic conductivity coefficients, H is groundwater table, Q is the flow rate for wells, is the Dirac delta which for wells' place equals to 1 and q is the distributed flow like precipitation and evaporation. is also recognized as specific yield.

With substituting Equations (20) into (19):
(20)
Equation (21) will derive:
(21)
As the homogeneous and isotropic condition is governed ():
(22)
At last, a set of linear equations in the form of are obtained. The size of these equations equals to the number of nodes scattered in the aquifer. This algebraic set of equations is solved by the Gauss-Seidel iteration method. In equations (2325), each matrix determinates (Mohtashami et al. 2020):
(23)
(24)
(25)

In the previous equations, W, and are weight function, shape function and time interval, respectively.

PF, known as a data assimilation method, is a powerful tool for the estimation of uncertain parameters in dynamic systems. Coupling this method with meshless groundwater modeling creates a model named PF-MLPG and leads to accurate results with high correspondence to a real condition. In this method, firstly, particles with equal weight values scattered in the state space. The dimension of state space depends on the number of uncertain parameters. In this study, 2,350 uncertain parameters including 1,175 hydraulic conductivity coefficients and 1,175 specific yield, with the aim of finding the best groundwater head, are scattered in the state space. The case study was Birjand unconfined aquifer located in the east of Iran. With using PF-MLPG, précised values of uncertain parameters were achieved. Figure 3 shows the flowchart of this work.

Figure 3

Flowchart of the PF-MLPG model.

Figure 3

Flowchart of the PF-MLPG model.

Close modal

Figure 3 shows the flowchart of the PF-MLPG model. In the first step, the upper and lower bounds of uncertain parameters are determined. Then a set of particles generated in the state space. The number of particles distributed in the state space depends on the number of uncertain parameters. The dimensions of state space are equal to the number of uncertain parameters. In the first iteration which K = 0, particles scattered in the state space with random locations, however, in the next iterations, based on the value of maximum likelihood function, particles' location is updated. All particles have a weight value which is computed by Equation (30). In the first step, this value is equal for all particles (). In the next steps, based on the maximum likelihood function, new weight values are assigned to the particles. Note that particles with greater weight values are closer to their exact ones.

Maximum likelihood function is calculated in the MLPG-MK model. This function is presented in the following equation (Field et al. 2016):
(26)

In Equation (26), and are the weight values in steps k and k − 1. , and are the observation, simulation and standard deviation, respectively.

To prevent from degeneracy occurrence in the PF method, a resampling technique is employed. In this situation, the particles with degenerated weight values are corrected and modified. The stop condition that is defined for the PF-MLPG model is the difference between the current weight values and the previous weight values. If this difference is lower than a negligible value. The stop criterion is satisfied and the program is terminated.

In this study, the hydrodynamic parameters of a real unconfined aquifer were calculated accurately using the PF-MLPG method. Unlike the real aquifer, it makes the condition challenging. In the latter, as parameters are more uncertain, the proposed model needs to be tested under conditions where the aquifer is real. Hence, Birjand aquifer was used.

As it was mentioned previously, PF-MLPG has two separate models, such as data assimilation and simulation model. The data assimilation model plays as an ‘online calibration’ role, and it uses the simulation model in order to reduce the uncertainty of hydrodynamic parameters, e.g. specific yield and hydraulic conductivity parameters.

In the simulation model, the input data are: wells discharge rate, hydraulic conductivity coefficient, specific yield and precipitation rate during the simulation period (1 year with monthly time step). Based on the reports, 10% of precipitation in this aquifer is considered as the recharge value. Discharge variation represented by 190 pumping wells (including 139 agricultural, 31 drinking water and 20 industrial) located inside the aquifer boundary (Figure 1) and also 11 piezometric wells, in which monthly observed groundwater level data established in the model. The only MLPG model is validated in Birjand unconfined aquifer in two steady and unsteady states. For having pervasive information about the simulation process and data included, Mohtashami et al. (2017a) and Mohtashami et al. (2017b) studies can be helpful. Note that the achieved hydraulic head from the steady state is used as an initial condition for transient state.

In the data assimilation approach, uncertain parameters for hydraulic conductivity and specific yield were 1,175 each (totally 2,350), and their variations range was from the minimum to the maximum observed limits, that is [5–65] for the former (hydraulic conductivity) and [0.02–0.17] for the latter (specific yield). Since running the PF-MLPG MATLB code for these conditions and these parameters were heavy, use was made of a 64core CPU, 2 T Ram system that took about 17 days to run; particles considered here were 4,000, 5,000 and 6,000.

To see the precise values of two uncertain parameters more easily, two following figures depicted. As shown in Figure 4, the acquired hydraulic conductivity coefficient by PF-MLPG is plotted. The color of each node represents the k value. Changes of k are more noticeable in the central part of aquifer. The north and southwest parts are in the same values. The blue color shows this fact as well. This figure also shows the northeast parts which have higher k values than the rest of the parts.

Figure 4

Representation of the achieved hydraulic conductivity coefficient in the aquifer by PF-MLPG. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.10.2166/hydro.2021.113.

Figure 4

Representation of the achieved hydraulic conductivity coefficient in the aquifer by PF-MLPG. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.10.2166/hydro.2021.113.

Close modal
Figure 5

Representation of the achieved specific yield in the aquifer by PF-MLPG. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.10.2166/hydro.2021.113.

Figure 5

Representation of the achieved specific yield in the aquifer by PF-MLPG. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.10.2166/hydro.2021.113.

Close modal

As shown in Figure 5, the specific yield variation trend is normal and the estimated values vary between the upper and lower bounds set for the program. The maximum specific yield value is 0.1666 and is assigned for node 530 (yellow color), and the minimum is 0.0230 for node 255 in the zone (dark blue color).

Next, to investigate this issue better, the values obtained for each node were placed between its upper and lower bounds and its relative frequency was determined. Precise hydraulic conductivity coefficients and specific yields achieved from the PF-MLPG model are plotted in the following figures. The accurate values depicted in Figures 6 and 7 accompanied by their upper and lower bounds. As it can be seen most of the values of k are between 5 and 15 m/day, and a few of them are placed between 50 and 60 m/day. This figure also shows that the number of nodes becomes greater when the values of hydraulic conductivity coefficient are reduced.

Figure 6

Frequency of estimated k between upper and lower bounds.

Figure 6

Frequency of estimated k between upper and lower bounds.

Close modal
Figure 7

Frequency of estimated Sy between upper and lower bounds.

Figure 7

Frequency of estimated Sy between upper and lower bounds.

Close modal

Figure 6 shows the hydraulic conductivity variations between the upper and lower limits considered in the program, and the maximum and minimum values are 59.55 and 5.01 m/day, respectively. Figure 6(b) shows the relative frequency of hydraulic coefficients in the mentioned parts.

In Figure 7, the aquifer nodes have had their highest values between 0.06 and 0.04 and fewer points have had specific yields higher than 0.11. Around 320 nodes have specific values of between 0.04 and 0.06. However, the numbers of nodes which have specific yield values higher than 0.1 are really few.

To evaluate the obtained results, the MLPG-MK simulation code was engaged and the groundwater level was calculated in the entire aquifer. Figure 8 compares the groundwater level found from the PF-MLPG method with those of the MLPG and finite difference methods (the MLPG chart input data are based on those validated in Sadeghi Tabas et al.’s (2016) research and shows the observational, PF-MLPG, FDM and MLPG results in, respectively, black, green, blue and red.

Figure 8

Comparison of groundwater table computed with different methods with the observation data acquired from six different piezometers. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.10.2166/hydro.2021.113.

Figure 8

Comparison of groundwater table computed with different methods with the observation data acquired from six different piezometers. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.10.2166/hydro.2021.113.

Close modal

As shown, PF-MLPG results and observational data are close and have relatively high correlation with piezometer 1 as well as 8. The highest difference quite visible in the PF-MLPG model is for piezometer 4; next is MLPG with higher accuracy than FDM. It shows a significant fall in piezometers 2 and 6 in months 4, 5 and 6. Beside comparing different numerical methods, effects of the number of particles (in the PF) on the groundwater were also investigated in each piezometer (Figure 9). The right number of particles has always been important in PF-related studies; if insufficient, estimations could be inaccurate, and if in excess, computational costs would increase. Since no studies have addressed the relations between the number of uncertain parameters and the number of particles exactly, 4,000, 5,000 and 6,000 particles were used to find the acceptable values of hydraulic conductivity and specific yield coefficients.

Figure 9

Effect of particle size in the computed groundwater head.

Figure 9

Effect of particle size in the computed groundwater head.

Close modal
As shown, the 6,000-particle modeling results are much better than those with fewer particles; same is the case with 5,000 compared to 4,000 particles, concluding that an increase in the number of particles in the algorithm will lead to more accurate estimations and, hence, a better aquifer modeling. Since an increase in the number of particles highly increases the PF-MLPG running time, results of only three particle groups have been studied; studying more than 6,000 was not possible. However, as graph-evaluation of the numerical method accuracy is hard, it was assessed, for simplicity, based on several other evaluation criteria defined, respectively, in relations 26–29. Tables 1 and 2 show the comparison on these error values in the calibration stage (Sadeghi tabas et al. 2016).
(27)
(28)
(29)
(30)
Table 1

Comparison of error values in the calibration stage

Error criteriaPF-MLPG with 4,000 particlesPF-MLPG with 5,000 particlesPF-MLPG with 6,000 particlesMLPGFDM
SI error 1.79 × 10−6 1.54 × 10−6 1.047 × 10−6 4.77 × 10−6 8.59 × 10−6 
Mean error (m) −0.024 −0.008 −0.0241 −0.12 0.159 
Mean absolute Error (m) 0.254 0.218 0.123 0.573 1.434 
Root mean Square error (m) 0.284 0.244 0.166 0.757 1.197 
Error criteriaPF-MLPG with 4,000 particlesPF-MLPG with 5,000 particlesPF-MLPG with 6,000 particlesMLPGFDM
SI error 1.79 × 10−6 1.54 × 10−6 1.047 × 10−6 4.77 × 10−6 8.59 × 10−6 
Mean error (m) −0.024 −0.008 −0.0241 −0.12 0.159 
Mean absolute Error (m) 0.254 0.218 0.123 0.573 1.434 
Root mean Square error (m) 0.284 0.244 0.166 0.757 1.197 

A comparison of the results in both validation and calibration stages can conclude that PF-MLPG methods have, regardless of the number of particles in the algorithm, a higher accuracy than MLPG and FDM, revealing the necessity of using the calibration–simulation model. Hydrodynamic parameters' accuracy has reduced the root mean square error from 0.757 to 0.166 m in the calibration stage and 1.021 to 0.178 in the validation stage, which are a significant reduction.

Table 2

Comparison of error values in the validation stage

Error criteriaPF-MLPG with 6,000 particlesMLPGFDM
SI error 1.247 × 10−6 10.23 × 10−6 12.63 × 10−6 
Mean error (m) 0.0452 0.42 0.349 
Mean absolute error (m) 0.344 0.893 1.848 
Root mean square error (m) 0.178 1.021 1.534 
Error criteriaPF-MLPG with 6,000 particlesMLPGFDM
SI error 1.247 × 10−6 10.23 × 10−6 12.63 × 10−6 
Mean error (m) 0.0452 0.42 0.349 
Mean absolute error (m) 0.344 0.893 1.848 
Root mean square error (m) 0.178 1.021 1.534 

Figure 10 shows each piezometer values to investigate the obtained regression coefficient between the computed groundwater table and the observed groundwater table.

Figure 10

Comparison of the simulated groundwater head and the observation groundwater head from the piezometers.

Figure 10

Comparison of the simulated groundwater head and the observation groundwater head from the piezometers.

Close modal

As shown, piezometer 1 has the best solution due to the higher regression coefficient (0.95), and the lowest value is for piezometer 4 (0.24) shown clearly in Figure 10.

Figure 11 shows the simulated groundwater head by PF-MLPG with 6,000 particles in the first and last period of simulation. As it can be seen, the groundwater head is higher in the east region of the aquifer as the color bar represents. The red color regions indicate the groundwater table between 1,380 and 1,400 m, with traveling to the west the level of groundwater decreases and reduced to its minimum value around 1,263 m, determined by blue color.

Figure 11

(a) Simulated groundwater head in the first period of simulation and (b) the last period of simulation. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.10.2166/hydro.2021.113.

Figure 11

(a) Simulated groundwater head in the first period of simulation and (b) the last period of simulation. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.10.2166/hydro.2021.113.

Close modal

In this study, a novel method that couples the meshless numerical (MLPG) simulation to a data assimilation method named PF has been introduced to investigate the accurate behavior of an unconfined aquifer in Birjand, the east of Iran. We derived the acceptable values of uncertain hydrodynamic coefficients, e.g. specific yield and hydraulic conductivities for the aquifer. Our results show that most of the nodes have specific yield and hydraulic conductivity values less than 0.05 and 10 m/day, respectively. That means most of the aquifer covered with the low hydraulic conductivity soil. We have calibrated and validated our model by comparing the ground water head derived from the model with the ground water head from the piezometers. The RMSE for the PF-MLPG model is calculated as 0.166 m, while it is 1.197 and 0.757 m for FDM and MLPG methods, respectively. This clearly shows the good estimation and higher accuracy of the hybrid PF-MLPG in predicting the aquifer behavior. We also performed a sensitivity analysis for the number of particles. Three particle sets were considered for this including 4,000, 5,000 and 6,000 particles. This analysis depicts that Results represent that increment of the particle's number led to higher accuracy, in other words RMSE for 4,000 particles, 5,000 and 6,000 were 0.284, 0.244 and 0.166 m, respectively. Furthermore, comparing the acquired groundwater heads from our method with the common methods in the literature shows the high capability of our method in the investigation of aquifer behavior. This accurate model helps the water engineers to define many scenarios for having better management about the aquifer. Some packages, e.g. capture zone delineation and optimization algorithm, also can be added to this model and make it multi-applicable for all the aquifers.

The authors thank Mrs Elieh Mohtashami for her insightful comments during the preparation of this paper.

1

Generalized Likelihood Uncertainty Estimation.

2

Moving Least Squares.

All relevant data are included in the paper or its Supplementary Information.

Abedini
M.
,
Ziai
A. N.
,
Shafiei
M.
,
Ghahraman
B.
,
Ansari
H.
&
Meshkini
J.
2017
Uncertainty assessment of groundwater flow modeling by using generalized likelihood uncertainty estimation method (case study: Bojnourd Plain)
.
Iranian Journal of Irrigation and Drainage
10
(
6
),
755
769
.
Akbari
A. R.
,
Bagri
A.
,
Bordas
S. P. A.
&
Rabczuk
T.
2010
Analysis of thermoelastic waves in a two-dimensional functionally graded materials domain by the meshless local Petrov-Galerkin (MLPG) method
.
Computer Modeling in Engineering and Sciences
65
(
1
),
27
74
.
Akbarpour
A.
,
Azizi
M.
&
Shirazi
M.
2012
Groundwater Management of Mokhtaran Aquifer with Using Finite Difference Mathematical Finite Difference
.
First Confernece of Groundwater
,
Tehran
.
Arulampalam
S.
,
Maskell
S.
,
Gordon
N.
&
Clapp
T.
2002
Tutorialon particle filters for online nonlinear/nongaussian Bayesian tracking
.
IEEE Transaction Signal Process
50
(
2
),
174
188
.
Atluri
S. N.
&
Zhu
T.
1998
A new meshless local Petrov-Galerkin (MLPG) approach in computational mechanics
.
Computational Mechanics
22
(
2
),
117
127
.
Chitsazan
M.
,
Abedini
M. J.
&
Salek
M.
2008
The investigation and quantifying the uncertainty of groundwater models in Kazeroun aquifer with using statistic parameters
.
Journal of Irrigation Sciences and Engineering (JISE)
19
,
17
33
.
Choe
G.
,
Wang
T.
,
Liu
F.
,
Hyon
S.
&
Ha
J.
2015
Particle filter with spline resampling and global transition model
.
IET Computer Vision
19
(
2
),
184
197
.
Eryiğit
M.
2021
Estimation of parameters in groundwater modelling by modified clonalg
.
Journal of Hydroinformatics
23
(
2
),
298
306
.
Field
G.
,
Tavrisov
G.
,
Brown
C.
,
Harris
A.
&
Kreidl
O. P.
2016
Particle filters to estimate properties of confined aquifers
.
Water Resources Management
30
,
3175
3189
.
Ghousehei
M.
2010
Groundwater Balance Computaion of Damghan Aquifer
.
University of Kermanshah
,
Kermanshah
.
Hamraz
B. S.
,
Akbarpour
A.
,
Pourreza Bilondi
M.
&
Sadeghi Tabas
S.
2015
On the assessment of ground water parameter uncertainty over an arid aquifer
.
Arabian Journal of Geosciences
8
,
10759
10773
.
Havangi
R.
2016
Increasing consistency of particle filter using the classic method and particle swarm algorithm
.
Computational Intelligence in Electrical Engineering
7
(
2
),
77
88
.
Khankham
S.
,
Luadsong
A.
&
Aschariyaphotha
N.
2015
MLPG method based on moving kriging interpolation for solving convection–diffusion equations with integral condition
.
Journal of King Saud University – Science
27
(
4
),
292
301
.
Li
T.
,
Bolic
M.
&
Djuric
P. M.
2015
Resampling methods for particle filtering: classification, implementation and strategies
.
IEEE Signal Processing Magazine
32
(
3
),
70
86
.
Liu
G. R.
&
Gu
Y. T.
2005
An Introduction to Meshfree Methods and Their Programming
.
Springer
,
Singapore
.
Manoli
G.
,
Rossi
M.
,
Pasetto
D.
,
Deiana
R.
,
Ferraris
S.
,
Cassiani
G.
&
Putti
M.
2015
An iterative particle filter approach for coupled hydro-geophysical inversion of a controlled infiltration experiment
.
Journal of Computational Physics
283
,
37
51
.
Markopoulos
A. P.
,
Karkalos
N. E.
&
Papazoglou
E. L.
2020
Meshless methods for the simulation of machining and micro-machining: a review
.
Archives of Computational Methods in Engineering
27
(
3
),
831
853
.
Mohtashami
A.
,
Akbarpour
A.
&
Mollazadeh
M.
2017a
Modeling of groundwater flow in unconfined aquifer in steady state with meshless local Petrov-Galerkin
.
Modares Mechanical Engineering
17
(
2
),
393
403
.
Mohtashami
A.
,
Hashemi Monfared
S. A.
,
Azizyan
G.
&
Akbarpour
A.
2020
Determination of the optimal location of wells in aquifers with an accurate simulation-optimization model based on the meshless local Petrov-Galerkin
.
Arabian Journal of Geosciences
13
(
2
),
1
13
.
Nossent
J.
,
Thouhidul Mostafa
S. M.
,
Ghysels
G.
&
Huysmans
M.
2020
Integrated Bayesian Multi-model approach to quantify input, parameter and conceptual model structure uncertainty in groundwater modeling
.
Environmental Modelling and Software
126
,
104654
110467
.
Pathania
T.
,
Bottacin-Busolin
A.
,
Rastogi
A. K.
&
Eldho
T. I.
2019
Simulation of groundwater flow in an unconfined sloping aquifer using the element-free Galerkin method
.
Water Resources Management
33
,
2827
2845
.
Pholl
G.
,
Pohlmann
K.
,
Hassan
A.
,
Chapman
J.
&
Mihvec
T.
2002
Assessing groundwater model uncertainty for the central Nevada test area (No. A-2002-01). Desert Research Institute, University and Community College System of Nevada (USA)
.
Porfiri
M.
2006
Analysis by Meshless Local Petrov-Galerkin Method of Material Discontinuities, Pull-in Instability in MEMS, Vibrations of Cracked Beams, and Finite Deformations of Rubberlike Materials
.
Virginia Polytechnic Institute and State University
,
Virginia
.
Ramgraber
M.
,
Albert
C.
&
Schirmer
M.
2019
Data assimilation and online parameter optimization in groundwater modeling using nested particle filters
.
Water Resources Research
55
(
11
),
9724
9747
.
Rasoulzadeh
A.
&
Mousavi
S. A. A.
2008
Using Inverse WTF Uncertainty Method in Estimation of Groundwater Model Parameters
.
Tehran
.
Ruknudeen
F.
&
Asokan
S.
2017
Application particle filter to on-board life estimation of LED lights
.
IEEE Photonics Journal
9
(
3
),
2017
.
Sadeghi tabas
S.
,
Samadi
S. Z.
,
Akbarpour
A.
&
Pourreza Bilondi
M.
2016
Sustainable groundwater modeling using single-and multi-objective optimization algorithms
.
Journal of Hydroinformatics
18
(
5
),
1
18
.
Swathi
B.
&
Eldho
T. I.
2017
Groundwater flow simulation in unconfined aquifers using meshless local Petrov-Galerkin method
.
Engineering Analysis with Boundary Elements
48
,
43
52
.
Thomas
A.
,
Eldho
T. I.
,
Rastogi
A. K.
&
Majumder
P.
2019
A comparative study in aquifer parameter estimation using MFree point collocation method with evolutionary algorithms
.
Journal of Hydroinformatics
21
(
3
),
455
473
.
Thouhidul Mostafa
S. M.
,
Mudud Hasan
M.
,
Kumar Saha
A.
,
Parvin Rannu
R.
,
Uytven
E. V.
,
Willems
P.
&
Huysmans
M.
2019
Multi-model approach to quantify groundwater-level prediction uncertainty using an ensemble of global climate models and multiple abstraction scenarios
.
Hydrology and Earth System Sciences
23
,
2279
2303
.
Wang
H. F.
&
Anderson
M. P.
1995
Introduction to Groundwater Modeling: Finite Difference and Finite Element Methods
, 1st edn.
Academic Press
.
Zhang
G. Y.
,
Cheng
Y. M.
,
Yang
F.
,
Pan
Q.
&
Liang
Y.
2010
Design of an adaptive particle filter based on variance reduction technique
.
Acta Automatica Sinica
36
(
7
),
1020
1024
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).