Numerical simulation of groundwater in an unconfined aquifer with a novel hybrid model (case study: Birjand Aquifer, Iran)

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. GLUE¹ 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.

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 km². 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.

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.

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).

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.

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.

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.

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.

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 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.

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.

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 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.

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.

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.

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

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.

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.

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