The complicated behavior of groundwater system in an arid aquifer is generally studied by solving the governing equations using either analytical or numerical methods. In this regard, analytical methods are just for some aquifers with regular boundaries. Numerical methods used for this aim are finite difference (FDM) and finite element methods (FEM) which are engaged for some simple aquifers. Using them in the complex cases with irregular boundaries has some shortcomings, dependent on meshes. In this study, meshless local Petrov-Galerkin (MLPG) method based on the moving kriging (MK) approximation function is used to simulate groundwater flow in steady state over three aquifers, two standard and a real field aquifer. Moving kriging function known as new function which reduces the uncertain parameter. For the first aquifer, a simple rectangular aquifer, MLPG-MK indicates good agreement with analytical solutions. In the second one, aquifer conditions get more complicated. However, MLPG-MK reveals results more accurate than FDM. RMSE for MLPG-MK and FDM is 0.066 and 0.322 m respectively. In the third aquifer, Birjand unconfined aquifer located in Iran is investigated. In this aquifer, there are 190 extraction wells. The geometry of the aquifer is irregular as well. With this challenging issues, MLPG-MK again shows satisfactory accuracy. As the RMSE for MLPG-MK and FDM are 0.483 m and 0.566 m. therefore, planning for this aquifer based on the MLPG-MK is closer to reality.

  • This Method (MLPG-MK) is applicable for all aquifers.

  • There is no limitation for this method.

     
  • MLPG

    Meshless local Petrov-Galerkin

  •  
  • MLS

    Moving Least Squares

  •  
  • MK

    Moving Kriging

  •  
  • FEM

    Finite Element Method

  •  
  • FDM

    Finite Difference Method

  •  
  • PPCM

    Polynomial Point Collocation Method

  •  
  • PCM

    Point Collocation Method

  •  
  • LRPIM

    Local Radial Point Interpolation Method

  •  
  • MQ-RBF

    Multi-Quartic-Radial Basis Function

  •  
  • ME

    Mean Error

  •  
  • MAE

    Mean Absolute Error

  •  
  • RMSE

    Root Mean Square Error

In recent years, due to population and industry growth, lack of precipitation caused by climate change and uncontrolled global and local pollution, finding clean, healthy and sufficient water resources is the big concern. In many parts of the world especially in Iran, having access to surface water resources is hard, therefore groundwater is counted as the main resource in this country.

As groundwater systems are complicated in nature, particularly in arid regions, it is important to recognize the groundwater temporal and spatial distribution when considering the pumping rate of extraction wells. Hence, simulation plays a vital role in planning and managing groundwater system (Swathi & Eldho 2014a, 2014b).

Groundwater model is a key element for presentation of the groundwater behavior in all type of aquifers confined, unconfined, standard or field aquifers (Anderson & Woessner 1992). This groundwater models solve governed equation of groundwater flow with using some numerical methods. Numerical methods e.g. FEM, FDM and boundary element method are the most common numerical methods which is used in this subject. Also, GMS (MODFLOW package) and FEFLOW are two popular software which are based on these numerical methods and engaged by some researchers for groundwater simulation (Harbaugh 2005; Zhao et al. 2005). Wang & Anderson (1982) used both FEM and FDM in order to simulate groundwater table for some aquifers with regular geometries. Their aquifers involved of confined and unconfined. They also investigated different type of boundary conditions on simulation results. Miller (2000) used MODFLOW for simulation of groundwater flow in a Karst aquifer in New York. The main purpose of his research was to provide the balance model of that aquifer (Miller 2000). Todd & Kenneth (2001) simulated groundwater flow in two transient and steady conditions. MODFLOW package was engaged for simulation purposes. They indicated that FDM was a suitable method for the modeling process. Rabinson & Reay (2002) determined the flow path of groundwater in a coastal aquifer in Virginia state by MODFLOW package. They planned a new strategy for management of that aquifer based on their simulation model. Zhao et al. (2005) investigated the vegetation changes in a dry aquifer based on the groundwater table in different depths. For this aim, they simulated groundwater flow by usage of FDM in FEFLOW software (Zhao et al. 2005). Trefry & Muffels (2007) stated that FEFLOW software is an extensive simulation software which helps water planners to have a better understanding about groundwater behavior. Zhou et al. (2010) used FDM for simulation of groundwater flow in order to compute drawdown. Based on that, the land subsidence was calculated and finally, endangered locations determined. Nguyen et al. (2013), Surinaidu et al. (2014), Karay & Hajnal (2015), Chen et al. (2016), Hashemi Monfared & Dehghani Darmian 2016; Hashemi Monfared et al. 2017; Darmian et al. 2018; Abdelhalim et al. (2019), Darmian et al. 2020; Al-Obaidi (2020a), Al-Obaidi (2020b), Al-Obaidi & Mishra (2020), Majidi Khalilabad et al. (2021) and Khodabandeh et al. (2021) used the same methods (FDM and FEM) for simulation of groundwater flow. As these are grid-based methods, despite their advantages, they have some dramatic shortcomings. The dependence of these methods on meshing over the entire domain is the main reason for the appearance of errors in their results (Grieble & Schweitzer 2000). In FEM, high time consumption and stability test, which forces the user to remesh the entire domain, are the shortcomings. In FDM, the incomplete coverage of mesh over the entire domain is the fundamental deficiency. These defects arise from the grid/mesh dependency of these methods in the solving process.

The idea of removing these problems is the beginning of meshless numerical methods (Liu 2002). Meshless methods are defined as a new numerical method employed to create a set of linear equations in the form of a matrix and solve these equations without considering a background or predefined mesh in the domain (Swathi & Eldho 2014a, 2014b). In these methods, the domain is comprised of a set of nodes, which are scattered uniformly or non-uniformly in the domain. The governed equation is enforced to each node. All the nodes are independent and have no connection among themselves. Recently, a myriad of meshfree methods have been employed for a wide range of engineering applications (Liu 2002). However, there are still some principal shortcomings that are related to the software packages employing meshless methods in the simulation process.

Few researches have used meshless models in groundwater simulation. Mategaonkar & Eldho (2011a) applied the polynomial point collocation method (PPCM) in two unconfined aquifers. These two were standard aquifers. The authors used the multiquadric radial basis function as an approximation function. Finally, they found their developed model satisfactory in comparison of FDM.

In another research, the point collocation method (PCM) was used for two hypothetical aquifers and one filed aquifer. Results in the hypothetical cases were verified by analytical solutions and they were acceptable. Simulation on the filed aquifer was also compared by boundary element method, and it was found that PCM was showing applicability in the simulation procedure (Mategaonkar & Eldho 2011b). Saeedpanah et al. (2011) developed a meshless model named the local radial point interpolation meshless (LRPIM) method with the aim of investigating leakage effects in a coastal confined aquifer. They employed two different functions for weight and approximation. Multi-quartic-radial basis function (MQ-RBF) and Heaviside were used for this purpose respectively. Kovarik & Muzik (2013) solved the density-driven groundwater flow equation using a combination of RBF interpolation and local boundary element (LBE) methods in a simple geometry region. Swathi & Eldho (2014b) used the MLPG method for simulation of groundwater flow in three unconfined aquifers. Exponential radial basis function was used for constructing of the shape function as well. The performance of MLPG-Exp-RBF was evaluated and it was stated that it could be engaged for all groundwater problems. In another study in the same year, Swathi and Eldho developed MLPG-Exp-RBF for confined aquifers. Three hypothetical cases were considered, a one-dimensional and two 2D aquifers with different boundary conditions. The same results as the previous study were achieved. In other words, MLPG-Exp-RBF showed good agreement results with the analytical ones. Abdi et al. (2021) simulated groundwater flow in a real aquifer with MLPG method; however, their approximation function was the moving least squares (MLS) method. They linked an optimization method to MLPG in order to find the best location of wells for excavation. Finally, they determined the best location for excavation of wells.

MLPG was firstly proposed by Atluri & Zhu (1998), for solving some problems, such as the potential equation and linear boundary problems. Application of the MLPG method indicated success of MLPG method in solving the advection-dispersion problems, Electro-Magnetics, Navier–Stokes flows, static and dynamic fracture mechanics problems, Plate bending problems, Water and shock wave problems and so on (Ching & Batra 2001; Cho & Atluri 2001; Long & Atluri 2002; Qian et al. 2004; Atluri 2005; Najafi et al. 2012; Mohtashami et al. 2019a; Al-Obaidi 2020c; Al-Obaidi 2021).

In the present study, an MLPG method based on the MK approximation function is developed for solving groundwater flow equations in an unconfined aquifer in steady state. The purpose of this research was to develop a new meshless based groundwater flow model in MALTAB platform. For testing its adequacy, a hypothetical problem and a real field problem is used and results of MLPG-MK is compared with finite different method and analytical solutions.

In this study, MLPG based on the MK approximation function is used for groundwater simulation in three unconfined aquifers. The independency of meshless methods to meshing stage, makes it more accurate and reliable for the simulation process. They remove errors related to meshes and make the model more compatible with irregular geometries. This method is applied in to two standard aquifers and one field aquifer. the properties of this cases will be explained in the next parts.

Meshless local Petrov-Galerkin (MLPG)

MLPG is known as a weak-form which enforces the governing equation in each node scattered in the domain and its boundary. This method was firstly introduced by Atluri & Zhu (1998) to solve some numerical equations. Cantilever beam, infinite plate with a circular hole, and a finite plate with an elliptical hole were the discussed examples. This is a truly meshless method, which does not depend on meshing in both integration and interpolation stages. MLPG is also allowed the users to define the test and trial functions differently, and it makes MLPG more compatible with the condition of the problem. In some problems, the changes in values of the filed variable is dramatic and, in some cases, it is gradually, MLPG with its capability make itself compatible with both conditions.

MLPG is based on the MLS approximation function. This function is used for constructing the shape function, and was firstly employed by Nayroles et al. (1992) for mathematical approximation. Despite its applicability, it suffers from satisfying Kronecker delta function. This matter prevents the essential boundary conditions to be enforced directly, unless some other techniques is used to overcome this drawback. Many studies used penalty function for enforcing essential boundary conditions, however this technique enters an uncertain parameter named penalty coefficient. In this study, to eliminate this limitation MK approximation function is engaged. MK and MLS have the same properties; however, MK satisfies the Kronecker delta function and provide the model to impose essential boundary conditions directly. The following section explains MK approximation function.

Groundwater flow equation

Groundwater flow equation in an unconfined aquifer in steady state for anisotropic condition is given in Equation (1):
formula
(1)

In which and are the hydraulic conductivity coefficients (m/day), H stands for the groundwater table (m), Q is denoted as discharge rate for extraction wells (m3/day/m2), is the Dirac delta function, which for the nodes that are defined by wells, equals 1 and, for the rest of the nodes becomes zero and q is the distributed discharge rate e.g. precipitation.

Regarding to Equation (2), that we have:
formula
(2)
Equation (1) can be written in the form of Equation (3):
formula
(3)
which is R defined as .
For isotropic condition, the hydraulic conductivity coefficient in both directions is equal, in other words .
formula
(4)
In order to establish the numerical procedure, weighted residual approach is employed. In this method the weighted integration is applied in the whole domain () (Mohtashami et al. 2021a):
formula
(5)
formula
(6)

In Equations (5) and (6), is denoted as a weighted value. It is computed based on the cubic Spline function, which is discussed in the next sections.

Two terms on the left side of Equation (6) are discretized with using integration by parts (Mohtashami et al. 2020):
formula
(7)
As the three investigated aquifers have Numann boundaries with zero value, the first and the third term of Equation (7), describing integration on boundary, equals to zero. Therefore Equation (8) is achieved:
formula
(8)
With substitution of Equation (2) in Equations (8) and (9) is derived:
formula
(9)
The approximated value for field variable of problem which in this study, it is groundwater table () is computed by a series term presented in Equation (10) (Zienkiewicz et al. 1977):
formula
(10)
which is shape function and is computed by MK, and, m is the number of gauss points. Using this equation into Equations (9) and (11) is presented as:
formula
(11)
Equation (11) is a set of linear equations and described as follows (Mohtashami et al. 2019b):
formula
(12)
where K is the stifness matrix, U is the field variable or groundwater table and F is the force body matrix, these matrices are explained in Equations (13)–(16) (Mohtashami et al. 2021b):
formula
(13)
formula
(14)
formula
(15)
formula
(16)

Noticed that is the number of all nodes scattered in the boundary and the domain.

Moving Kriging method (MK)

The Moving kriging method formulation is similar to MLS approximation function, which is used by MLPG in order to construct the shape function. This function satisfies the Kronecker delta property. Consider the function is defined as groundwater head. Therefore, the formulation of the shape function is presented in Equation (17) (Liu & Gu 2005):
formula
(17)
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 (Khankham et al. 2015):
formula
(18)
where A and B are expressed as (Khankham et al. 2015):
formula
(19)
formula
(20)
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 Equations (21) and (22) (Khankham et al. 2015):
formula
(21)
formula
(22)
where is the correlation function between any pair of nodes located at and .

Weight function

Choice of Weight function impacts on the performance of simulation. The Weight function in the simulation procedure must have at least three properties (Liu & Gu 2005). 1- It must be positive. 2- It must be zero outside the computational domain. 3- W(x) must be decreased monotonically from the interest point x.

In this study, based on the previous researches, the cubic Spline function is used. The formula of that is presented in Equation (23) (Liu & Gu 2005):
formula
(23)

Therefore, in a glance of view, the flowchart of simulation is determined in Figure 1.

Figure 1

Flowchart of MLPG-MK in the simulation process.

Figure 1

Flowchart of MLPG-MK in the simulation process.

Close modal

As can be seen, input data including discharge rate, hydraulic conductivity coefficient and boundary conditions are entered into the model. Based on Equations (14) and (16), stiffness and force body matrices are computed. Then they assemble in a matrix presented in Equation (13). After enforcement of boundary conditions, a set of linear equation system is solved by Gauss-seidel iteration method. At last, groundwater head for all the nodes is achieved.

In this section, the properties of three aquifers and their simulation procedure, including pre-processing, processing and post-processing are discussed. The first and the second aquifer are standard and the third is a real unconfined aquifer with more complexity. Firstly, the groundwater flow in standard aquifers would be simulated, and then the real aquifer would be modeled.

Case study 1

A simple rectangular unconfined aquifer is presented in Figure 2. The size of the aquifer is 15 × 15 m2. It has no extraction and observation wells. The boundary conditions for all sides of this aquifer are constant head boundary. The south and west edges have zero value; however the north and east sides are computed with and formulas respectively. The analytical solution for this problem is presented as . The steady condition is governed in this case study. Hydraulic conductivity coefficient is also equal to 1 m/day. Figure 2 shows the aquifer with its boundary conditions. The first step of the simulation is scattering nodes in both the boundary and domain of the problem. Twenty-five nodes are spread across the domain while two successive points have a uniform distance of 2.5 m in Figure 2, blue points are the nodes scattered across the boundary and the red points represent the nodes in the domain.

Figure 2

Representation of the first aquifer and uniform scattering of nodes.

Figure 2

Representation of the first aquifer and uniform scattering of nodes.

Close modal

Groundwater head is computed for the aquifer. Figure 3(a) shows head contours and their values on it. Moving to the north-east, the groundwater level increases significantly. Figure 3(b) is the 3-dimensional representation of the groundwater table.

Figure 3

Results of MLPG-MK method for the first aquifer.

Figure 3

Results of MLPG-MK method for the first aquifer.

Close modal

Head in the south and west edges remains zero due to the zero constant head, and it follows the and formulas for north and east boundaries, respectively.

In order to evaluate the performance of the MLPG-MK, the achieved results are compared with the analytical solution. Table 1 shows the comparison of the MLPG-MK results with the analytical solution.

Table 1

The comparison of MLPG solution and exact solution for the 2D rectangular aquifer

X(m)Y(m)MLPG-MK results (m)Analytical results (m)X(m)Y(m)MLPG-MK results (m)Analytical results (m)
2.5 166.67 166.66 7.5 500.01 500.00 
7.5 2.5 250.01 250.00 7.5 7.5 750.01 750.00 
10 2.5 333.34 333.33 10 7.5 1,000.01 1,000.00 
333.34 333.33 10 666.67 666.66 
7.5 500.01 500.00 7.5 10 1,000.01 1,000.00 
10 666.67 666.66 10 10 1,333.34 1,333.33 
X(m)Y(m)MLPG-MK results (m)Analytical results (m)X(m)Y(m)MLPG-MK results (m)Analytical results (m)
2.5 166.67 166.66 7.5 500.01 500.00 
7.5 2.5 250.01 250.00 7.5 7.5 750.01 750.00 
10 2.5 333.34 333.33 10 7.5 1,000.01 1,000.00 
333.34 333.33 10 666.67 666.66 
7.5 500.01 500.00 7.5 10 1,000.01 1,000.00 
10 666.67 666.66 10 10 1,333.34 1,333.33 

As can be seen, the MLPG-MK results are really close to the analytical solutions. The difference between MLPG-MK and the analytical solution is around 0.01 m, which is rare and not meaningful.

Case study 2

A standard unconfined aquifer with the area of 4,500 × 10,000 m2 is derived from Mckinney & Lin (1994) research (Figure 4). Boundary conditions for this aquifer in the south and north edges are a constant head boundary with 20 m value. East and west sides are no flow boundaries. The aquifer is investigated in steady state with isotropic and homogeneous conditions. Ten extraction wells are located in this aquifer, whose coordinates are specified in Table 2. The precipitation rate is 0.001 m/day. and the hydraulic conductivity value is 50 m/day as well.

Table 2

Extraction rate of wells in the second standard aquifer

WellLocation (m)Flow rate (m3/day)WellLocation (m)Flow rate (m3/day)
(1,250, 8,500) 7,000 (3,250, 5,000) 5,949.65 
(2,250, 8,500) 7,000 (1,250, 3,500) 6,729.79 
(3,250, 8,500) 7,000 (2,000, 3,500) 4,282.12 
(1,250, 5,000) 5,960.43 (2,500, 3,500) 4,230.97 
(2,250, 5,000) 4,530.70 (3,250, 3,500) 6,807.09 
WellLocation (m)Flow rate (m3/day)WellLocation (m)Flow rate (m3/day)
(1,250, 8,500) 7,000 (3,250, 5,000) 5,949.65 
(2,250, 8,500) 7,000 (1,250, 3,500) 6,729.79 
(3,250, 8,500) 7,000 (2,000, 3,500) 4,282.12 
(1,250, 5,000) 5,960.43 (2,500, 3,500) 4,230.97 
(2,250, 5,000) 4,530.70 (3,250, 3,500) 6,807.09 
Figure 4

Representation of the second aquifer with 10 extraction wells.

Figure 4

Representation of the second aquifer with 10 extraction wells.

Close modal

In this standard aquifer, firstly, nodes are distributed uniformly within 250 m distance in X and Y directions. Each node has information including the hydraulic conductivity coefficient, precipitation and flow rate. A discharge rate based on Table 1 is assigned to the extraction wells. Wells a, b and c have the highest amount of discharge among the other extraction wells. The value of precipitation which is given, allocated to all nodes, is scattered in the aquifer. The boundary conditions in this domain are constant head and constant flow. The left and right side of the aquifer is the impervious boundary condition known as constant flow with zero value. The north and south edges have constant head values. As is depicted, 116 nodes were placed on the boundaries with blue color, and 2,867 nodes were located in the aquifer with red color.

Groundwater head is computed for each node solving the set of linear equations with the Gauss-Seidel method. In Figure 5, the groundwater head near the northern and southern edges, which are constant head boundaries, remain constant and equal to 20 m. The red color also displays this fact. Drawdown near extraction wells is clearly determined. Extraction wells with the highest values of flow rate have the highest drawdown. Wells h and i with the least discharge rate have the lowest drawdown.

Figure 5

Groundwater head computation with MLPG model.

Figure 5

Groundwater head computation with MLPG model.

Close modal

In order to evaluate the accuracy of the MLPG model, the achieved results are compared with analytical and FDM solutions. Analytical result which are the exact solution are derived from McKinney & Lin (1994). FDM solutions also are presented in Mategaonkar & Eldho's (2011b) research. Table 3 presented this comparison.

Table 3

Comparison of MLPG and FDM solutions with analytical one

WellLocation (m)Flow rate (m3/day)Analytical solution (m) McKinney & Lin (1994) FDM (m) Mategaonkar & Eldho (2011b) MLPG-MK solution(m)
(1,250, 8,500) 7,000 12.18 11.99 12.10 
(2,250, 8,500) 7,000 11.4 11.34 11.42 
(3,250, 8,500) 7,000 12.18 11.62 12.08 
(1,250, 5,000) 5,960.43 0.12 0.89 0.21 
(2,250, 5,000) 4,530.70 0.33 0.47 0.37 
(3,250, 5,000) 5,949.65 0.03 0.08 0.05 
(1,250, 3,500) 6,729.79 1.22 1.03 1.11 
(2,000, 3,500) 4,282.12 0.34 0.41 0.38 
(2,500, 3,500) 4,230.97 0.68 0.51 0.72 
(3,250, 3,500) 6,807.09 0.17 0.23 0.21 
WellLocation (m)Flow rate (m3/day)Analytical solution (m) McKinney & Lin (1994) FDM (m) Mategaonkar & Eldho (2011b) MLPG-MK solution(m)
(1,250, 8,500) 7,000 12.18 11.99 12.10 
(2,250, 8,500) 7,000 11.4 11.34 11.42 
(3,250, 8,500) 7,000 12.18 11.62 12.08 
(1,250, 5,000) 5,960.43 0.12 0.89 0.21 
(2,250, 5,000) 4,530.70 0.33 0.47 0.37 
(3,250, 5,000) 5,949.65 0.03 0.08 0.05 
(1,250, 3,500) 6,729.79 1.22 1.03 1.11 
(2,000, 3,500) 4,282.12 0.34 0.41 0.38 
(2,500, 3,500) 4,230.97 0.68 0.51 0.72 
(3,250, 3,500) 6,807.09 0.17 0.23 0.21 

In Table 3, comparing FDM and MLPG-MK solutions with analytical responses, the high accuracy of the MLPG-MK model is clear. This fact indicates the capability of MLPG-MK in simulating groundwater flow in steady state aquifers. This matter is also shown in Figure 6.

Figure 6

Groundwater flow in the second aquifer.

Figure 6

Groundwater flow in the second aquifer.

Close modal
Three common error criteria are calculated. Mean error (ME), mean absolute error (MAE), and root mean square error (RMSE) are these criteria. These indices present error in the units (or squared units) and are computed using Equations (24)–(26) (Mohtashami et al. 2019b).
formula
(24)
formula
(25)
formula
(26)
where , and are the observed and simulated groundwater heads, respectively, and m is the number of piezometers in the aquifer. The calculated errors are presented in Table 4.
Table 4

Computed ME, MAE and RMSE errors

Error criteriaFDM (m)MLPG-MK (m)
ME 0.104 0.0043 
MAE 0.404 0.058 
RMSE 0.322 0.066 
Error criteriaFDM (m)MLPG-MK (m)
ME 0.104 0.0043 
MAE 0.404 0.058 
RMSE 0.322 0.066 

In these indices, RMSE is the main error criterion which shows the model performance more suitably. This criterion for MLPG-MK is lower than FDM, as it is 0.066 m and 0.322 m, respectively. MLPG-MK reduces the error around 80 percent. This reduction is meaningful and also shows the applicability of the MLPG-MK groundwater model.

As it is investigated in both standard aquifers, MLPG-MK has good and satisfactory results compared to the analytical solutions. Therefore this model must be tested in a real aquifer with an irregular boundary and challenging problems like changes in hydraulic conductivity coefficients and discharge rates. The next section is the enforcement of MLPG-MK in a field aquifer.

Real case study

Birjand unconfined aquifer is located in South Khorasan province in Iran. The area of this aquifer is approximately 269 km2. This aquifer is known as an arid region according to the Domarton index. Mean annual precipitation of this aquifer is 170 mm (Mohtashami et al. 2017). One hundred and ninety extraction and 11 observation wells are in this aquifer, which is shown in the following figure with blue and red symbols. There are ten regions, which have constant head boundaries in the aquifer and are determined with arrows in Figure 7.

Figure 7

(a) and (b) Geographical location of the Birjand aquifer, (c) wells, (d) hydraulic conductivity coefficients (Sadeghi-Tabas et al. 2017).

Figure 7

(a) and (b) Geographical location of the Birjand aquifer, (c) wells, (d) hydraulic conductivity coefficients (Sadeghi-Tabas et al. 2017).

Close modal

Hydraulic conductivity coefficients for this aquifer are derived from Sadeghi Tabas et al.’s (2017) research, which is calibrated by Cuckoo optimization algorithm and MODFLOW package (Figure 7(d)).

For simulation purposes, 1,175 nodes are scattered in the aquifer at a uniform distance. The difference between two adjacent points is 500 m. This value is chosen based on the Sadeghi Tabas et al. (2017) and Hamraz et al. (2016) researches, which provide us the possibility to compare our results with their works. All the input information, such as hydraulic conductivity coefficients, discharge rates, boundary conditions, and precipitation is assigned for each node. The hydraulic conductivity coefficients map is depicted in Figure 7(d). As was mentioned before, the boundary condition in this aquifer is only of the type of constant head boundary. Ten fronts are recognized as this type. These fronts including one outflow pathway and nine inflow pathways, which are determined in Figure 7(c) with black arrows. Based on the reports of the regional water company of South Khorasan, all precipitation is not infiltrated to the aquifer, and about 17 percent must be considered as recharge rate (Report 1132).

Finally, the groundwater table is computed for each node in the aquifer in a steady state. Notice that 1,175 nodes are placed on the boundary and the domain. To ascertain the performance of the MLPG-MK method, its results are compared with observed head and FDM solutions. Notice that FDM solutions are derived from Sadeghi-Tabas et al. (2017). Table 5 shows the comparison of these solutions.

Table 5

The comparison of MLPG-MK, FDM solutions, and observed data in simulation period (2012–2013)

WellX(m)Y(m)MLPG-MK head (m)FDM head (m) (Sadeghi Tabas et al. 2017)Observed head (m)
672,076.92 3,626,500 1,264.41 1,263.895 1,264.30 
673,616.684 3,629,000 1,291.85 1,290.572 1,291.55 
674,670.794 3,638,500 1,306.21 1,306.881 1,306.87 
675,659.263 3,634,500 1,296.93 1,296.745 1,296.43 
677,358.12 3,628,000 1,300.30 1,300.55 1,300.77 
681,191.541 3,638,000 1,310.13 1,309.874 1,309.98 
684,659.633 3,637,500 1,321.30 1,321.806 1,322.41 
693,716.317 3,641,500 1,342.05 1,341.247 1,342.24 
696,160.839 3,639,500 1,357.22 1,356.793 1,357.71 
701,775.426 3,639,000 1,362.87 1,362.955 1,362.9 
716,167.142 3,636,000 1,392.39 1,392.449 1,392.2 
WellX(m)Y(m)MLPG-MK head (m)FDM head (m) (Sadeghi Tabas et al. 2017)Observed head (m)
672,076.92 3,626,500 1,264.41 1,263.895 1,264.30 
673,616.684 3,629,000 1,291.85 1,290.572 1,291.55 
674,670.794 3,638,500 1,306.21 1,306.881 1,306.87 
675,659.263 3,634,500 1,296.93 1,296.745 1,296.43 
677,358.12 3,628,000 1,300.30 1,300.55 1,300.77 
681,191.541 3,638,000 1,310.13 1,309.874 1,309.98 
684,659.633 3,637,500 1,321.30 1,321.806 1,322.41 
693,716.317 3,641,500 1,342.05 1,341.247 1,342.24 
696,160.839 3,639,500 1,357.22 1,356.793 1,357.71 
701,775.426 3,639,000 1,362.87 1,362.955 1,362.9 
716,167.142 3,636,000 1,392.39 1,392.449 1,392.2 

The results obtained from the MLPG-MK method are very close, more than FDM, to the observed data. It means that the MLPG-MK solutions are more reliable and accurate than FDM.

In four piezometers (a, b, h and k), MLPG-MK model results are more accurate than FDM. These piezometers based on Figure 7(c) are placed in the south west, centre and east of the aquifer where they have extraction wells around themselves. (Figure 8).

Figure 8

Comparing MLPG-MK, FDM, and observation results in each piezometers.

Figure 8

Comparing MLPG-MK, FDM, and observation results in each piezometers.

Close modal

Finally, based on the groundwater table of nodes, the Birjand unconfined aquifer is interpolated. The results of interpolation are depicted in Figure 9. As can be seen, east of the aquifer has the highest groundwater table around 1,380 m, with blue color, and travelling to the west the groundwater table decreases gradually by around 1,280 m (determined by red color).

Figure 9

Groundwater level in the Birjand aquifer in a steady state.

Figure 9

Groundwater level in the Birjand aquifer in a steady state.

Close modal

Performance of MLPG-MK is again evaluated based on the previous error indices. Table 6 shows the error values.

Table 6

Computed ME, MAE and RMSE errors

Error criteriaMLPG-MK (m)FDM (m)
ME 0.234 0.321 
MAE 0.381 0.404 
RMSE 0.483 0.566 
Error criteriaMLPG-MK (m)FDM (m)
ME 0.234 0.321 
MAE 0.381 0.404 
RMSE 0.483 0.566 

According to Table 5 and Figure 8, RMSE for the MLPG-MK method is closer to zero than the FDM method; therefore, the MLPG method is more accurate and effective than FDM in simulation of groundwater in the Birjand aquifer.

As the purpose of the model is to predict groundwater behavior for the future, validation must be done.

In order to validate our model, we computed the groundwater table for the next year (2013–2014). Therefore, some of the wells extraction rate and precipitation rate change. Changes are applied into the model. Finally, the groundwater head is computed for the aquifer. the achieved results from the MLPG-MK model and observation data are in Table 7.

Table 7

MLPG-MK results in validation stage

WellX(m)Y(m)MLPG-MK head (m)Observed head (m)
672,076.92 3,626,500 1,264.45 1,264.29 
673,616.684 3,629,000 1,291.78 1,291.51 
674,670.794 3,638,500 1,305.45 1,306.77 
675,659.263 3,634,500 1,295.84 1,296.38 
677,358.12 3,628,000 1,299.21 1,300.69 
681,191.541 3,638,000 1,309.05 1,309.93 
684,659.633 3,637,500 1,320.18 1,322.32 
693,716.317 3,641,500 1,339.98 1,342.18 
696,160.839 3,639,500 1,354.12 1,357.68 
701,775.426 3,639,000 1,362.81 1,362.88 
716,167.142 3,636,000 1,391.29 1,392.16 
WellX(m)Y(m)MLPG-MK head (m)Observed head (m)
672,076.92 3,626,500 1,264.45 1,264.29 
673,616.684 3,629,000 1,291.78 1,291.51 
674,670.794 3,638,500 1,305.45 1,306.77 
675,659.263 3,634,500 1,295.84 1,296.38 
677,358.12 3,628,000 1,299.21 1,300.69 
681,191.541 3,638,000 1,309.05 1,309.93 
684,659.633 3,637,500 1,320.18 1,322.32 
693,716.317 3,641,500 1,339.98 1,342.18 
696,160.839 3,639,500 1,354.12 1,357.68 
701,775.426 3,639,000 1,362.81 1,362.88 
716,167.142 3,636,000 1,391.29 1,392.16 

Based on Table 7, MLPG-MK shows good agreement with the observation data of the next year. Table 7 also indicates the ascending trend of the groundwater table compared to the previous year. Regression coefficient for the simulated head in the verification and validation stage is computed and depicted in Figure 10.

Figure 10

Regression coefficient in (a) verification and (b) validation stages.

Figure 10

Regression coefficient in (a) verification and (b) validation stages.

Close modal

R2 for verification stage is 0.9998; however, it is 0.9992 for validation. The difference is ignorable and not meaningful; however, this occurred due to the changes in precipitation or discharge rates of some wells. Therefore, MLPG-MK presents a reliable, accurate and powerful method for groundwater simulation.

This model (MLPG-MK) is written in MATLAB v.2018 software and is run on a home use laptop with I7,4700/8Gb Ram properties. The time needed to get the results for the real aquifer is about 30 minutes, which is a bit high. This is the only concern about MLPG-MK. However, considering the reliability and accuracy of that, the consumed time would be ignorable.

Sensitivity analysis

In simulation procedure, in order to integrate the numerical integration explained in Equations (14) and (16), an integration domain named the quadrature domain is considered. The size of the quadrature domain, which is computed by Equation (27), plays an important role in the accuracy of the results. This size depends on an uncertain dimensionless parameter called . This coefficient, based on Liu & Gu (2005), shows better results when ranging from 1.5 to 2.5. To find the best value of this coefficient, we carried out a sensitivity analysis (Liu & Gu 2005).
formula
(27)
where, di is the size of the quadrature domain and is the distance between two successive points, which is 500 m for Birjand aquifer. In this research, different values of are tested to find its best value. The RMSE error is also computed. Table 8 determines the sensitivity analysis results for .
Table 8

The effect of quadrature domain size on accuracy

Well
1,264.45 1,264.32 1,264.41 1,265.07 1,265.91 
1,298.79 1,293.84 1,291.85 1,291.37 1,284.32 
1,307.68 1,306.54 1,306.21 1,306.13 1,305.41 
1,307.75 1,298.81 1,296.93 1,296.54 1,291.06 
1,304.30 1,303.71 1,300.30 1,298.97 1,282.30 
1,313.56 1,313.41 1,310.13 1,307.92 1,293.95 
1,285.58 1,324.21 1,321.30 1,320.79 1,313.67 
1,268.37 1,343.54 1,342.05 1,341.61 1,340.78 
1,337.71 1,357.51 1,357.22 1,356.67 1,354.98 
1,356.21 1,363.09 1,362.87 1,362.89 1,362.56 
1,391.96 1,392.41 1,392.39 1,392.24 1,392.32 
RMSE (m) 14.81 1.82 0.483 1.077 8.366 
Well
1,264.45 1,264.32 1,264.41 1,265.07 1,265.91 
1,298.79 1,293.84 1,291.85 1,291.37 1,284.32 
1,307.68 1,306.54 1,306.21 1,306.13 1,305.41 
1,307.75 1,298.81 1,296.93 1,296.54 1,291.06 
1,304.30 1,303.71 1,300.30 1,298.97 1,282.30 
1,313.56 1,313.41 1,310.13 1,307.92 1,293.95 
1,285.58 1,324.21 1,321.30 1,320.79 1,313.67 
1,268.37 1,343.54 1,342.05 1,341.61 1,340.78 
1,337.71 1,357.51 1,357.22 1,356.67 1,354.98 
1,356.21 1,363.09 1,362.87 1,362.89 1,362.56 
1,391.96 1,392.41 1,392.39 1,392.24 1,392.32 
RMSE (m) 14.81 1.82 0.483 1.077 8.366 

It is clear that for , the results are more accurate than the other values. This table proves this best choice of . The values in the range of 1.5–2.5 have fewer oscillations.

This research developed a meshless method in the 2D case in MATLAB software named MLPG-MK. MLPG, based on moving kriging, is used for the simulation of groundwater flow in three aquifers. The moving kriging function was employed as an approximation function to reduce the uncertain parameters of the old MLS approximation function. In addition to tackling the drawbacks of the meshes, this method can model domains with complex geometries. This method was applied in three aquifers, two standard aquifers and one real aquifer. In the first standard aquifer, which has a simple condition including no extraction well, no precipitation, and regular geometry, the groundwater table is computed. The achieved results were compared with analytical solutions and they show satisfactory accuracy. In the second aquifer, some conditions are more complicated. Ten extraction wells exist in the aquifer; boundary conditions are classified in two types. Precipitation is also non-zero in the aquifer. Simulation results were compared with analytical and FDM solutions. MLPG-MK is more accurate than FDM. In the real aquifer, Birjand unconfined aquifer, located in the east of Iran, has been investigated. In this aquifer, there are 190 extraction wells and 11 observation wells. the geometry is also irregular. In this challenging condition, the performance of the MLPG-MK is tested. Finally, the results were found to be satisfactory. The simulated groundwater table by MLPG-MK is compared with the results of FDM and observation wells. The acquired results show the high accuracy and capability of MLPG-MK more than grid-based methods as the RMSE for MLPG-MK and FDM are 0.483 m and 0.566 m respectively. This accurate model helps the decision makers to manage groundwater resources better because they can predict the behavior of aquifers accurately under many scenarios and can plan a schedule over the exploitation of groundwater.

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

Abdelhalim
A.
,
Sefelnasr
A.
&
Ismail
E.
2019
Numerical modeling technique for groundwater management in Samalut city, Minia Governature, Egypt
.
Arabian Journal of Geosciences
12
(
124
),
1
18
.
Abdi
M.
,
Ebrahimi
H.
&
Akbarpour
A.
2021
Optimal location of pumping wells by a mesh-free numerical method
.
Water Supply
2021, ws2021379. doi: https://doi.org/10.2166/ws.2021.379
.
Al-Obaidi
A. R.
2020b
Experimental comparative investigations to evaluate cavitation conditions within a centrifugal pump based on vibration and acoustic analyses techniques
.
Archives of Acoustics
45
(
3
),
541
556
.
Anderson
M.
,
Woessner
W.
&
Hunt
R.
2015
Applied Groundwater Modeling
, 2nd edn.
Academic Press, Cambridge, MA
.
Atluri
S. N.
2005
Methods of Computer Modeling in Engineering & the Sciences
, Vol.
1
.
Tech Science Press
,
Palmdale, FL
.
Atluri
S.
&
Zhu
T. A.
1998
A new MEshless method (MLPG) approach in computational mechanics
.
Computaional Mechanics
22
(
2
),
117
127
.
Chen
S.
,
Yang
W.
,
Huo
G.
&
Huang
Z.
2016
Groundwater simulation for efficient water resources management in Zhangye Oasis, Northwest China
.
Environmental Earth Science
75
(
647
),
1
13
.
Ching
H. K.
&
Batra
R. C.
2001
Determination of crack tip fields in linear elastostatics by the meshless local Petrov-Galerkin(MLPG) method
.
CMES- Computer Modeling in Engineering and Sciences
2
(
2
),
273
289
.
Cho
J. Y.
&
Atluri
S. N.
2001
Analysis of shear flexible beams, using the meshless local Petrov-Galerkin method, based on a locking-free formulation
.
Engineering Computations
18
(
1/2
),
215
240
.
Darmian
M. D.
,
Monfared
S. A. H.
,
Azizyan
G.
,
Snyder
S. A.
&
Giesy
J. P.
2018
Assessment of tools for protection of quality of water: uncontrollable discharges of pollutants
.
Ecotoxicology and Environmental Safety
161
,
190
197
.
Darmian
M. D.
,
Khodabandeh
F.
,
Azizyan
G.
,
Giesy
J. P.
&
Monfared
S. A. H.
2020
Analysis of assimilation capacity for conservation of water quality: controllable discharges of pollutants
.
Arabian Journal of Geosciences
13
(
17
),
1
13
.
Grieble
M.
&
Schweitzer
M. A.
2000
A particle-Partition of unity method for the solution of Eliptic, parabolic, and Hyperbolic PDES
.
SIAM Journal on Scientific Computing
22
(
3
),
853
890
.
Hamraz
B.
,
Akbarpour
A.
,
Bilondi
M. P.
&
Tabas
S. S.
2016
On the assessment of ground water parameter uncertainty over an arid aquifer
.
Arabian Journal of Geosciences
8
(
12
),
10759
10773
.
Harbaugh
A. W.
2005
MODFLOW-2005, the US Geological Survey Modular Ground-Water Model: the Ground-Water Flow Process
.
US Department of the Interior, US Geological Survey
,
Reston, VA
, pp.
6
A16
.
Hashemi Monfared
S. A.
&
Dehghani Darmian
M.
2016
Evaluation of appropriate advective transport function for one-dimensional pollutant simulation in rivers
.
International Journal of Environmental Research
10
(
1
),
77
84
.
Hashemi Monfared
S. A.
,
Darmian
M. D.
,
Snyder
S. A.
,
Azizyan
G.
,
Pirzadeh
B.
&
Moghaddam
M. A.
2017
Water quality planning in rivers: assimilative capacity and dilution flow
.
Bulletin of Environmental Contamination and Toxicology
99
(
5
),
531
541
.
Karay
G.
&
Hajnal
G.
2015
Modelling of groundwater flow in fractured Rocks
. In
7th Groundwater Symposium of the International Association for Hydro-Environment Engineering and Research (IAHR)
.
Procedia Environmental Sciences
.
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
.
Khodabandeh
F.
,
Darmian
M. D.
,
Moghaddam
M. A.
&
Monfared
S. A. H.
2021
Reservoir quality management with CE-QUAL-W2/ANN surrogate model and PSO algorithm (case study: Pishin Dam, Iran)
.
Arabian Journal of Geosciences
14
(
5
),
1
18
.
Kovarik
K.
&
Muzik
J.
2013
A meshless solution of two dimensional density-driven groundwater flow
.
Engineering Analysis with Boundary Elements
37
,
187
196
.
Liu
G.
2002
Mesh Free Methods: Moving Beyond the Finite Element Method
.
CRC Press
,
Boca Raton, FL
.
Liu
G. R.
&
Gu
Y. T.
2005
An Introduction to Meshfree Methods and Their Programming
.
Springer Science & Business Media, New York, NY
.
Long
S.
&
Atluri
S. N.
2002
A meshless local Petrov-Galerkin method for solving the bending problem of a thin plate
.
Computer Modeling in Engineering and Sciences
3
(
1
),
53
64
.
Majidi Khalilabad
N.
,
Mohtashami
A.
&
Akbarpour
A.
2021
Determination of well's capture zones using random walk algorithm and FeFlow simulation model
.
Iranian Journal of Irrigation & Drainage
14
(
6
),
1984
2002
.
Mategaonkar
M.
&
Eldho
T. I.
2011a
Meshless point collocation method for Ld and 2d groundwater flow simulation
.
Ish Journal Of Hydraulic Engineering
17
(
1
),
71
87
.
Mategaonkar
M.
&
Eldho
T. I.
2011b
Simulation of groundwater flow in unconfined aquifer using meshfree point collocation method
.
Engineering Analysis with Boundary Elements
35
(
4
),
700
707
.
McKinney
D. C.
&
Lin
M. D.
1994
Genetic algorithm solution of groundwater management models
.
Water Resources Research
30
(
6
),
1897
1906
.
Miller
T. S.
2000
Simulation of Ground-Water Flow in an Unconfined Sand and Gravel Aquifer at Marathon
.
USGS
,
New York
,
Cortland County
, pp.
1
29
.
Mohtashami
A.
,
Hashemi Monfared
S. A.
,
Azizyan
G.
&
Akbarpour
A.
2019a
Determination the capture zone of wells by using meshless local Petrov-Galerkin numerical model in confined aquifer in unsteady state (Case study: Birjand Aquifer)
.
Iranian Journal of Ecohydrology
6
(
1
),
239
255
.
Mohtashami
A.
,
Hashemi Monfared
S. A.
,
Azizyan
G.
&
Akbarpour
A.
2019b
Prediction of groundwater fluctuations using Meshless local Petrov-Galerkin numerical method in a field aquifer (Birjand Aquifer)
.
International Journal of Numerical Methods in Civil Engineering
3
(
4
),
33
41
.
Mohtashami
A.
,
Monfared
S. A. H.
,
Azizyan
G.
&
Akbarpour
A.
2021a
Usage of particle filter for exact estimation of constant head boundaries in unconfined aquifer
.
Amirkabir Journal of Civil Engineering
53
(
11
),
1
20
.
Mohtashami
A.
,
Hashemi Monfared
S. A.
,
Azizyan
G.
&
Akbarpour
A.
2021b
Estimation of parameters in groundwater modeling by particle filter linked to the meshless local Petrov-Galerkin numerical method
.
Journal of Hydraulic Structures
7
(
1
),
16
37
.
Najafi
M.
,
Arefmanesh
A.
&
Enjilela
V.
2012
Meshless local Petrov–Galerkin method-higher Reynolds numbers fluid flow applications
.
Engineering Analysis with Boundary Elements
36
(
11
),
1671
1685
.
Nayroles
B.
,
Touzot
G.
&
Villon
P.
1992
Generalizing the finite element method: diffuse approximation and diffuse elements
.
Computational Mechanics
10
,
307
318
.
Nguyen
T. T.
,
Tsujimura
M.
&
Naoaki
S.
2013
Groundwater flow modeling: considering water use in Tay, Island, Dong Thap Province, Southwest Vietnam
. In:
The 3th International Conference on Sustainable Future For Human Security SUSTAIN
.
Kyoto University
,
Japan
.
Sadeghi-Tabas
S.
,
Samadi
S. Z.
,
Akbarpour
A.
&
Pourreza-Bilondi
M.
2017
Sustainable groundwater modeling using single-and multi-objective optimization algorithms
.
Journal of Hydroinformatics
19
(
1
),
97
114
.
Saeedpanah
I.
,
Jabbari
E.
&
Shayanfar
M. A.
2011
Numerical simulation of ground water flow via a new approach to the local radial point interpolation meshless method
.
International Journal of Computational Fluid Dynamics
25
(
1
),
17
30
.
Swathi
B.
&
Eldho
T. I.
2014a
Groundwater flow simulation in confined aquifers using meshless Local Petrov-Galerkin
.
ISH Journal of Hydraulic Engineering
19
,
335
348
.
Swathi
B.
&
Eldho
T. I.
2014b
Groundwater flow simulation in unconfined aquifers using meshless local Petrov-Galerkin method
.
Engineering Analysis with Boundary Elements
48
,
43
52
.
Trefry
M. G.
&
Muffels
C.
2007
FEFLOW: A finite-element ground water flow and transport modeling tool
.
Groundwater
45
(
5
),
525
528
.
Wang
H. F.
&
Anderson
M. P.
1982
Introduction to Groundwater Modeling: Finite Difference and Finite Element Methods
, 1st edn.
Academic Press
.
Zhao
C.
,
Wang
Y.
,
Chen
X. I.
&
Li
B.
2005
Simulation of the effects of groundwater level on vegetation change by combining FEFLOW software
.
Ecological Modeling
187
(
2–3
),
341
351
.
Zhou
N.
,
Vermeer
P.
,
Lou
R.
,
Tang
Y.
&
Jiang
S.
2010
Numerical simulation of deep foundation pit dewatering and optimization of controling land subsidece
.
Engineering Geology
114
,
251
260
.
Zienkiewicz
O. C.
,
Taylor
R. L.
,
Nithiarasu
P.
&
Zhu
J. Z.
1977
The Finite Element Method
, Vol.
3
.
McGraw-hill
,
London
.
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/).