In recent decades, due to reduction in precipitation, groundwater resource management has become one of the most important issues considered to prevent loss of water. Many solutions are concerned with the investigation of groundwater flow behavior. In this regard, development of meshless methods for solving the groundwater flow system equations in both complex and simple aquifers' geometry make them useful tools for such investigations. The independency of these methods to meshing and remeshing, as well as its capability in both reducing the computation requirement and presenting accurate results, make them receive more attention than other numerical methods. In this study, meshless local Petrov–Galerkin (MLPG) is used to simulate groundwater flow in Birjand unconfined aquifer located in Iran in a transient state for 1 year with a monthly time step. Moving least squares and cubic spline are employed as approximation and weight functions respectively and the simulated head from MLPG is compared to the observation results and finite difference solutions. The results clearly reveal the capability and accuracy of MLPG in groundwater simulation as the acquired root mean square error is 0.757. Also, with using this method there is no need to change the geometry of aquifer in order to construct shape function.

INTRODUCTION

Understanding groundwater flow behavior is of high importance among hydrogeologists. They have tried to investigate the aquifer behavior in all aspects due to the low amount of precipitation in many places in the world. The way to recognize this action is to study the related governing equations. Therefore simulation of groundwater flow system, especially in arid zones is essential for the government and water scientists because of increasing rates of demand on groundwater reservoirs (Sadeghi Tabas et al. 2016).

Technically, the governed equations of groundwater flow can be solved by some numerical methods such as finite difference method (FDM) and finite element method (FEM). These methods are mesh based and in addition to their benefits they have some drawbacks related to the meshes. Many scientists have tried to resolve these shortages but FDM and FEM have their problems in the field of meshing. Using rectangular meshes makes the FDM amongst the most simple method for programming. However, it has its own limitations as meshes may not cover the whole domain precisely. Moreover, some problems concerned with changeable boundaries and constant update of the domain geometry is a time-consuming and complex task. In order to address such problems, a newly numerical method, called the meshless method, has been developed. The meshless method is a numerical technique to solve problems simply. This method is used to change partial differential equations into a system of equations for the entire domain without the use of a mesh (Mategaonkar & Eldho 2011). Recently, there has been great attention focused on the development of ‘meshless methods’ to reduce meshing problems (Li et al. 2003). The first idea of meshless methods was introduced by Gingold & Monaghan (1977). They used smoothed particle hydrodynamics for modeling specified phenomena. The main idea of these methods is to approximate a function in the entire domain by just using scattered nodes. Atluri & Zhu (1998) presented a new meshless method called meshless local Petrov–Galerkin (MLPG). They developed local symmetric weak form to solve Laplace and Poisson equations and found this method capable and accurate (Atluri & Zhu 1998).

There has been little research into groundwater flow modeling using meshless methods. Mategaonkar & Eldho (2011) used a meshless method based on point collocation called polynomial point collocation method for the groundwater flow simulation in unconfined aquifers (Mategaonkar & Eldho 2011). Their case study regions include two presumptive aquifers in one- and two-dimensional and one real aquifer. Their results in hypothetical aquifers were compared with the FEM and analytical solutions. Kovarik & Muzik (2013) used a combination of two numerical methods, radial basis function (RBF) and the local boundary integral element method, to solve the unsteady density-driven groundwater flow in a rectangular case study (Kovarik & Muzik 2013). Swathi & Eldho (2013) simulated groundwater flow in a two-dimensional (2D) presumptive rectangular confined aquifer with MLPG with Gaussian RBF. They found their used method effective (Swathi & Eldho 2014). Swathi & Eldho (2014) also simulated groundwater flow in an unconfined aquifer with MLPG with Gaussian RBF as weight and shape functions. They chose both weight function and approximation function from one space. Also, to keep the size of all support domains equal they considered some dummy nodes in the aquifer (Swathi & Eldho 2014).

Recently, numerous studies have been related to the simulation of groundwater flow in the Birjand unconfined aquifer. Sadeghi Tabas et al. (2016) presented research that linked multi-algorithm genetically adaptive search method with a groundwater model to define pumping rates within a well-distributed set of Pareto solutions (Sadeghi Tabas et al. 2016). Hamraz et al. (2015) assessed parameter uncertainty of groundwater simulation in the Birjand aquifer. They modeled Birjand aquifer in Matlab using MODFLOW. They evaluated parameter uncertainty using GLUE (generalized likelihood uncertainty estimation). Their results showed that the performance of the GLUE was satisfactory (Hamraz et al. 2015). Ghoochanian et al. (2013) developed a model that linked MODFLOW and WEAP to manage water resources in the Birjand aquifer.

In this paper, the MLPG method is employed to solve the time-dependent groundwater flow equation in Birjand unconfined aquifer located in the east of Iran. The shape and weight functions are chosen from two different spaces. The weight function is cubic spline, moreover, a moving least squares (MLS) approximation function has been chosen as interpolation function. Finally, the simulated head is compared with the observation results and FDM (MODFLOW) solutions for each time step.

METHODS

Geographical location of the case study

The Birjand unconfined aquifer is located in South Khorasan province in the east of Iran. The area of this aquifer is almost 265 km2. Mean annual precipitation in Birjand plain is 160 mm and generally this region is placed in an arid zone. The slope is shallow in the west part of the plain. The groundwater system is one-layer unconfined aquifer with changeable thickness of 5–225 m. (Sadeghi Tabas et al. 2016). Figure 1 shows the geographical location of the Birjand aquifer in Iran.
Figure 1

Geographical location of Birjand aquifer.

Figure 1

Geographical location of Birjand aquifer.

MLS method

MLS approximations were firstly introduced in the field of computational solid mechanics by means of the diffuse element method by Nayroles et al. (1992). Since then, this function has been utilized in the element-free Galerkin method by Belytschko et al. (1994). Now, it is an alternative approach for constructing meshless shape functions for approximation (Liu & Gu 2005). Due to two important characteristics, much attention has been attracted toward the MLS approximation: (1) the approximated field function is continuous and smooth in the whole problem domain when adequate nodes are used; and (2) it is able to produce an approximation with the desired order of consistency (Liu 2002).

The hydraulic head (UI) of groundwater in an aquifer is needed to compute. The MLS approximation of hydraulic head can be defined as (Park & Leap 2000): 
formula
1
is the number of the basis function and is a vector of unknown coefficient given by: 
formula
2
In Equation (1), the basis function is built using monomials from the Khayyam–Pascal triangle. In one-dimensional space, a complete polynomial basis of order l is given by: 
formula
3
and in 2D space: 
formula
4
In this study, the quadratic basis is used.
The functional J must be minimized with respect to: 
formula
5
where are the nodal unknowns associated with the neighbors nodes of point x and is a weight function of ith node whose value decreases as the distance between X and increases (Liu 2002).
The minimization condition for Equation (5) requires: 
formula
6
which leads to the linear equation system: 
formula
7
Assuming that the MLS moment matrix A is invertible, Equation (7) can then be solved for : 
formula
8
, and are computed: 
formula
9
 
formula
10
 
formula
11
Substituting Equation (9) back into Equation (1): 
formula
12
In the simpler form: 
formula
13
where is the matrix of MLS shape function.

Choice of the weight function

Choosing the weight function plays a significant role in the efficiency of the MLS approximation function (Liu & Gu 2005). The choice of is made. It has the following features: it is only strictly positive in a subdomain containing , but is zero outside this domain. is called the support of the function or the domain of influence of the node . reduces with the distance between and X. Several weight functions can be used. In this study, the cubic spline weight function is used: 
formula
14

In this equation and is the influence radius of node. For each node, must be selected in a way that the number of non-zero weight functions be more than each term in the polynomial.

MLPG method

This method is a truly meshless method, as it does not need a ‘mesh’, either for purposes of interpolation of the solution variables, or for the integration stage (Atluri & Zhu 2000). This method employs a local weak form to solve the equations. MLS approximation function is utilized for approximating. MLPG was first introduced by Atluri & Zhu (1998).

Discretization of groundwater flow equation in unconfined aquifer with MLPG method

The governing equation in groundwater flow modeling in unconfined aquifer is: 
formula
15
where , and indicate the hydraulic conductivity coefficients (L/T). H, and R show the pressure head (L), specific storage (1/L) and recharge or discharge rate per unit volume (1/T) (positive or negative) respectively. In the unconfined aquifer, the thickness of saturated layer varies by groundwater table. Some assumptions proposed by Dupuit (1863) are:
  • (a)

    the flow is horizontal;

  • (b)

    the hydraulic gradient is equal to slope of the free surface.

The equation is based on the Dupuit (1863) assumption and continuity equation as follows: 
formula
16
where is specific yield.
Since: 
formula
17
with substitution of Equation (17) in Equation (16): 
formula
18
Since Birjand unconfined aquifer is isotropic: 
formula
19
with using weighted residual method: 
formula
20
 
formula
21
By using integration of parts in Equation (21): 
formula
22
Since the Birjand aquifer has no normal flow over its boundary, the first and third term in the left side of Equation (22) vanish: 
formula
23
Substituting Equation (17) in Equation (23): 
formula
24
The approximate value for head is: 
formula
25
Substituting Equation (25) in Equation (24): 
formula
26
In this study, the FDM in time is used. Forward finite difference approximation is employed for computation of the time derivative: 
formula
27
Substituting Equation (27) in Equation (26): 
formula
28
 
formula
29
This equation is similar to the linear equation: 
formula
30
 
formula
31
 
formula
32

Equation (32) represents force body matrix. It indicates the rate of recharge or discharge of the aquifer. In other words, this matrix shows the interaction or extraction flow rate in two concentrated and distributed conditions.

The principles of groundwater modeling

Birjand aquifer was represented by using sets of nodes which were scattered uniformly in it. The word uniform means that the distance between two adjacent points in horizontal and vertical direction is 500 meters ). Input data, such as extraction wells, recharge and discharge rates, boundary conditions, hydraulic conductivity coefficient and specific yield of each node, were entered into the model. The algorithm of simulating with MLPG is presented in Figure 2.
Figure 2

Flowchart of groundwater flow modeling.

Figure 2

Flowchart of groundwater flow modeling.

The scattered nodal points in the domain are presented in Figure 3. Each blue point in the domain has four values: hydraulic head that should be computed, hydraulic conductivity coefficient, specific yield and the rate of recharge or discharge (please refer to the online version of this paper to see Figure 3 in color: http://dx.doi.org/10.2166/hydro.2017.024).
Figure 3

Scattering nodal points in simulated aquifer in Matlab software.

Figure 3

Scattering nodal points in simulated aquifer in Matlab software.

Boundary conditions

Generally, there are two classifications for boundary conditions in modeling of groundwater flow. One is nodes with a specified head (Dirichlet) boundary condition, and the other one is nodes with no flow or inactive (Neuman) boundary condition. Inactive nodes or no flow nodes are those for which no flow into or out of the nodes is permitted. Constant head nodes are those for which the head is specified in advance, and is held at this specified value through all time steps of the simulation. In the Birjand aquifer, there are 10 areas that have specified head boundary condition, nine inflow pathways and one outflow pathway. Figure 4 specifies these areas. The other boundary nodes have no flow boundary conditions, in other words they are defined as inactive nodes.
Figure 4

Presented inflow and outflow pathways in Birjand aquifer.

Figure 4

Presented inflow and outflow pathways in Birjand aquifer.

Wells

There are 190 extraction wells in the studied region. These pumping wells consist of 139 agricultural, 31 drinking water and 20 industrial wells. In this aquifer 10 observation wells were used to investigate the level of water table. Figures 5 and 6 show the location of extraction and observation wells with square and circle symbol respectively.
Figure 5

Extraction wells in Birjand aquifer.

Figure 5

Extraction wells in Birjand aquifer.

Figure 6

The location of piezometers (observation wells) in the aquifer.

Figure 6

The location of piezometers (observation wells) in the aquifer.

Recharge and discharge rate

Since the studied area is categorized as an arid region with low amounts of precipitation, this low precipitation is considered as the recharge value. This aquifer had 0.000727 m/day rain between 2011 and 2012 based on rain gauges employed in Birjand plain.

The volume of extracted water from extraction wells is used as the discharge rate in our proposed model.

Hydraulic conductivity coefficient and specific yield

In order to allocate hydraulic conductivity coefficient and specific yield in each node, the aquifer has been divided into polygons. This is done by ArcGIS software. For each polygon there is one value that represents the magnitude of the hydraulic conductivity coefficient and specific yield. All the nodes in each polygon have the same value as the polygon. The unit of hydraulic conductivity coefficient is m/day. Figures 7 and 8 indicate a divided aquifer for presenting hydraulic conductivity coefficient and specific yield parameters. These parameters are obtained from the results of pumping tests.
Figure 7

Hydraulic conductivity polygons in Birjand aquifer.

Figure 7

Hydraulic conductivity polygons in Birjand aquifer.

Figure 8

Representation of specific yield polygons in Birjand aquifer.

Figure 8

Representation of specific yield polygons in Birjand aquifer.

According to Figure 8, Birjand unconfined aquifer is divided into 48 polygons. Each polygon has a number as a legend. The value of the specific yield parameter for a polygon can be obtained based on its legend number from Table 1.

Table 1

Values of specific yield in polygons based on its number

Num Specific yield Num Specific yield Num Specific yield 
0.25 16 0.06 32 0.04 
0.02 17 0.05 33 0.07 
0.04 18 0.04 34 0.17 
0.03 19 0.064 35 0.04 
0.03 20 0.03 36 0.12 
0.054 21 0.06 37 0.04 
0.03 22 0.1 38 0.1 
0.054 23 0.03 39 0.05 
0.03 24 0.04 40 0.04 
0.065 25 0.1 41 0.08 
10 0.07 26 0.04 42 0.04 
11 0.04 27 0.03 43 0.1 
12 0.075 28 0.05 44 0.04 
13 0.045 29 0.1 45 0.1 
14 0.07 30 0.04 46 0.04 
15 0.043 31 0.04 47 0.035 
Num Specific yield Num Specific yield Num Specific yield 
0.25 16 0.06 32 0.04 
0.02 17 0.05 33 0.07 
0.04 18 0.04 34 0.17 
0.03 19 0.064 35 0.04 
0.03 20 0.03 36 0.12 
0.054 21 0.06 37 0.04 
0.03 22 0.1 38 0.1 
0.054 23 0.03 39 0.05 
0.03 24 0.04 40 0.04 
0.065 25 0.1 41 0.08 
10 0.07 26 0.04 42 0.04 
11 0.04 27 0.03 43 0.1 
12 0.075 28 0.05 44 0.04 
13 0.045 29 0.1 45 0.1 
14 0.07 30 0.04 46 0.04 
15 0.043 31 0.04 47 0.035 

MODFLOW model

The groundwater flow equation is solved by a finite difference scheme using a groundwater modeling system (GMS). MODFLOW is a three-dimensional finite-difference groundwater model and engages certain equations for computing the head of groundwater in different cells of the aquifer. The model was first released in 1984 by the United States Geological Survey (USGS). Now, it is used as a popular package for simulation of groundwater flow in both steady and unsteady states. Modeling groundwater flow with a MODFLOW package requires accurate data from the aquifer.

The first step to simulate groundwater flow in GMS is to define a conceptual model. This is done to simplify the analysis of the field data in real conditions (Anderson & Woessner 1991). All the meteorological, hydrological and geophysics studies, geological investigation, well logs, water level fluctuation, type of groundwater flow system, and hydrodynamic coefficients are employed in the conceptual model. The conceptual model consists of definition of boundary, type and structure of aquifer and sources and sinks, etc. The next step is to choose an algorithm for solving the related equation. In the fourth stage the conceptual model is converted to a numerical model. This stage includes griding domain, definition of time steps and implied boundary conditions. Finally, the model is calibrated. Calibration means the accommodation of simulation results with observation results.

Sadeghi Tabas et al. (2016) and Hamraz et al. (2015) modeled Birjand unconfined aquifer using the MODFLOW model. They developed their model based on available data, including well locations and measurements, geologic map, well logs, topography data, hydrography and recharge information.

Seven layers were used to construct the groundwater model. Aquifer boundary conditions, piezometers, pumping wells, surface recharge, drainage information, hydraulic conductivity and specific yield are applied to the MODFLOW model to create a groundwater numerical model (Sadeghi Tabas et al. 2016).

RESULTS AND DISCUSSION

The groundwater flow in the Birjand unconfined aquifer is simulated with the MLPG and finite difference (MODFLOW) methods. The water level is computed for every node in each time step. To investigate the accuracy of the results, observed head in piezometers (observation wells) was compared with MLPG and FDM solutions. Figure 9 shows the comparison of these numerical methods and observed results for four random piezometers.
Figure 9

Computed water level from MLPG, FDM in comparison with observed head.

Figure 9

Computed water level from MLPG, FDM in comparison with observed head.

Figure 10

The variability of hydraulic head in (a) the first stress period (b) the last stress period.

Figure 10

The variability of hydraulic head in (a) the first stress period (b) the last stress period.

Tables 2 and 3 show the simulated head with using MLPG and FDM in comparison of observation data in each time step for two piezometers (piezometer 1 and 3).

Table 2

Comparison in results of MLPG, FDM and observation in piezometer 1

Stress period (month) MLPG head (m) FDM head (m) Observed head (m) 
1,264.34 1,263.756 1,264.07 
1,264.33 1,263.341 1,264.08 
1,264.33 1,263.182 1,264.11 
1,264.27 1,263.203 1,264.14 
1,264.26 1,263.253 1,264.17 
1,264.26 1,263.31 1,264.29 
1,264.363 1,263.377 1,264.27 
1,264.362 1,263.444 1,264.12 
1,264.363 1,263.511 1,264.14 
10 1,264.366 1,263.576 1,264.22 
11 1,264.37 1,263.638 1,264.34 
12 1,264.37 1,263.698 1,264.39 
Stress period (month) MLPG head (m) FDM head (m) Observed head (m) 
1,264.34 1,263.756 1,264.07 
1,264.33 1,263.341 1,264.08 
1,264.33 1,263.182 1,264.11 
1,264.27 1,263.203 1,264.14 
1,264.26 1,263.253 1,264.17 
1,264.26 1,263.31 1,264.29 
1,264.363 1,263.377 1,264.27 
1,264.362 1,263.444 1,264.12 
1,264.363 1,263.511 1,264.14 
10 1,264.366 1,263.576 1,264.22 
11 1,264.37 1,263.638 1,264.34 
12 1,264.37 1,263.698 1,264.39 
Table 3

Comparison in results of MLPG, FDM and observation in piezometer 3

Stress period (month) MLPG head (m) FDM head (m) Observed head (m) 
1,307.145 1,308.44 1,307.29 
1,307.134 1,308.474 1,307.23 
1,307.134 1,308.479 1,307.18 
1,306.927 1,308.476 1,307.11 
1,306.928 1,308.473 1,307.07 
1,306.928 1,308.471 1,307.04 
1,307.168 1,308.466 1,306.97 
1,307.166 1,308.458 1,306.94 
1,307.166 1,308.448 1,306.89 
10 1,307.176 1,308.437 1,306.85 
11 1,307.175 1,308.427 1,306.8 
12 1,307.175 1,308.418 1,306.77 
Stress period (month) MLPG head (m) FDM head (m) Observed head (m) 
1,307.145 1,308.44 1,307.29 
1,307.134 1,308.474 1,307.23 
1,307.134 1,308.479 1,307.18 
1,306.927 1,308.476 1,307.11 
1,306.928 1,308.473 1,307.07 
1,306.928 1,308.471 1,307.04 
1,307.168 1,308.466 1,306.97 
1,307.166 1,308.458 1,306.94 
1,307.166 1,308.448 1,306.89 
10 1,307.176 1,308.437 1,306.85 
11 1,307.175 1,308.427 1,306.8 
12 1,307.175 1,308.418 1,306.77 

Figures 9 and 10 show the change of hydraulic head along the aquifer in the first and last stress period (month) of 2010–2011.

The groundwater level has a maximum value of approximately 1,390 m in the east of the aquifer, and it decreases gradually while traveling to the west. The south-west has a minimum value around 1,264 m.

Performance of the model

There are three criterion that were used to evaluate the error of the results. Mean error (ME), mean absolute error (MAE) and root mean square error (RMSE). Error estimation was computed using Equations (33)–(35) (Sadeghi Tabas et al. 2016): 
formula
33
 
formula
34
 
formula
35
where , indicates the observed and simulated heads respectively while ‘m’ is the number of monthly time steps and ‘n’ is the number of piezometers. These mentioned errors are computed in Table 4.
Table 4

ME, MAE and RMSE errors

 MLPG method (m) FDM method (MODFLOW) (m) 
ME –0.08 0.159 
MAE 0.573 1.434 
RMSE 0.757 1.197 
 MLPG method (m) FDM method (MODFLOW) (m) 
ME –0.08 0.159 
MAE 0.573 1.434 
RMSE 0.757 1.197 

Actually, the allowable error associated with the field and model data should be under +1.9 m. In this study, the objective error is RMSE, according to Table 4 the RMSE of using MLPG and FDM methods are 0.757 and 1.197. They have significant differences. So the MLPG method reveals results with more accuracy than the FDM method.

CONCLUSIONS

In recent decades, reduction in precipitation, population growth and civilization developments caused an increase in underground water withdrawal and endangered the groundwater flow system balance. In this regard, water resources management has become one of the main issues under investigation for better resources management. The best way for managing water resources is to use a numerical method. In this study, Birjand unconfined aquifer with complex geometry located in east of Iran is modeled using MLPG. The shape and weight functions are chosen from two different spaces and cubic spline is employed as the weight function. Moreover, the MLS approximation function is used as the approximation function. The simulated head by MLPG was compared with the results of observation data and FDM (MODFLOW) solution. For four piezometers as an example, comparison among the solutions of these methods was carried out. Also, the RMSE for MLPG solution was computed as 0.757. This value shows the high accuracy of the MLPG in simulating groundwater flow.

REFERENCES

REFERENCES
Anderson
M.
Woessner
W.
1991
Applied Groundwater Modeling Simulation of Flow and Advective Transport
,
1st edn
.
Academic Press
,
USA
.
Atluri
S. N.
Zhu
T. L.
1998
A new meshless method (MLPG) approach in computational mechanics
.
Comput. Mech.
22
(
2
),
117
127
.
Belytschko
T.
Lu
Y. Y.
Gu
L.
1994
Elements free Galerkin methods
.
Int. J. Numer. Meth. Eng.
30
(
2
),
229
256
.
Dupuit
J.
1863
Estudes Theoriques et Pratiques sur le Mouvement des Eaux
.
Dunod
,
Paris
.
Ghoochanian
E.
Etebari
B.
Akbarpour
A.
2013
Integrating groundwater management with WEAP and MODFLOW models (case study: Birjand Plain, east of Iran). Modflow and More: Translating Science into Practice – IGWMC
.
Gingold
R. A.
Monaghan
J. J.
1977
Smoothed particle hydrodynamics: theory and application to non-spherical stars
.
Mon. Not. R. Astron. Soc.
181
(
3
),
375
389
.
Hamraz
B. S.
Akbarpour
A.
Pourreza Bilondi
M.
Sadeghi Tabas
S.
2015
On the assessment of ground water parameter uncertainty over an arid aquifer
.
Arab. J. Geosci.
8
(
12
),
10759
10773
.
Kovarik
K.
Muzik
J.
2013
A meshless solution of two dimensional density-driven groundwater flow
.
Eng. Anal. Bound. Elem.
37
,
187
196
.
Li
J.
Chen
Y.
Pepper
D.
2003
Radial basis function method for 1-D and 2-D groundwater
.
Comput. Mech.
32
(
1
),
10
15
.
Liu
G.
2002
Mesh Free Methods: Moving Beyond the Finite Element Method
.
CRC Press
,
Boca Raton
.
Liu
G. R.
Gu
Y. T.
2005
An Introduction to Meshfree Methods and Their Programming
.
Springer
,
Singapore
.
Sadeghi Tabas
S.
Samadi
S. Z.
Akbarpour
A.
Pourreza Bilondi
M.
2016
Sustainable groundwater modeling using single-and multi-objective optimization algorithms
.
J. Hydroinform.
19
(
1
),
97
114
.
Swathi
B.
Eldho
T. I.
2013
Groundwater flow simulation in confined aquifers using meshless local Petrov-Galerkin (MLPG) method
.
ISH Journal of Hydraulic Engineering
19
(
3
),
335
348
.