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

^{2}. 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.

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

*U*) of groundwater in an aquifer is needed to compute. The MLS approximation of hydraulic head can be defined as (Park & Leap 2000): is the number of the basis function and is a vector of unknown coefficient given by:

_{I}*J*must be minimized with respect to:where are the nodal unknowns associated with the neighbors nodes of point x and is a weight function of

*i*th node whose value decreases as the distance between

*X*and increases (Liu 2002).

*A*is invertible, Equation (7) can then be solved for :, and are computed:

### Choice of the weight function

*X*. Several weight functions can be used. In this study, the cubic spline weight function is used:

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

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

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

### Boundary conditions

### Wells

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

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.

Num . | Specific yield . | Num . | Specific yield . | Num . | Specific yield . |
---|---|---|---|---|---|

0 | 0.25 | 16 | 0.06 | 32 | 0.04 |

1 | 0.02 | 17 | 0.05 | 33 | 0.07 |

2 | 0.04 | 18 | 0.04 | 34 | 0.17 |

3 | 0.03 | 19 | 0.064 | 35 | 0.04 |

4 | 0.03 | 20 | 0.03 | 36 | 0.12 |

5 | 0.054 | 21 | 0.06 | 37 | 0.04 |

6 | 0.03 | 22 | 0.1 | 38 | 0.1 |

7 | 0.054 | 23 | 0.03 | 39 | 0.05 |

8 | 0.03 | 24 | 0.04 | 40 | 0.04 |

9 | 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 | 0.25 | 16 | 0.06 | 32 | 0.04 |

1 | 0.02 | 17 | 0.05 | 33 | 0.07 |

2 | 0.04 | 18 | 0.04 | 34 | 0.17 |

3 | 0.03 | 19 | 0.064 | 35 | 0.04 |

4 | 0.03 | 20 | 0.03 | 36 | 0.12 |

5 | 0.054 | 21 | 0.06 | 37 | 0.04 |

6 | 0.03 | 22 | 0.1 | 38 | 0.1 |

7 | 0.054 | 23 | 0.03 | 39 | 0.05 |

8 | 0.03 | 24 | 0.04 | 40 | 0.04 |

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

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

Stress period (month) . | MLPG head (m) . | FDM head (m) . | Observed head (m) . |
---|---|---|---|

1 | 1,264.34 | 1,263.756 | 1,264.07 |

2 | 1,264.33 | 1,263.341 | 1,264.08 |

3 | 1,264.33 | 1,263.182 | 1,264.11 |

4 | 1,264.27 | 1,263.203 | 1,264.14 |

5 | 1,264.26 | 1,263.253 | 1,264.17 |

6 | 1,264.26 | 1,263.31 | 1,264.29 |

7 | 1,264.363 | 1,263.377 | 1,264.27 |

8 | 1,264.362 | 1,263.444 | 1,264.12 |

9 | 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 | 1,264.34 | 1,263.756 | 1,264.07 |

2 | 1,264.33 | 1,263.341 | 1,264.08 |

3 | 1,264.33 | 1,263.182 | 1,264.11 |

4 | 1,264.27 | 1,263.203 | 1,264.14 |

5 | 1,264.26 | 1,263.253 | 1,264.17 |

6 | 1,264.26 | 1,263.31 | 1,264.29 |

7 | 1,264.363 | 1,263.377 | 1,264.27 |

8 | 1,264.362 | 1,263.444 | 1,264.12 |

9 | 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 | 1,307.145 | 1,308.44 | 1,307.29 |

2 | 1,307.134 | 1,308.474 | 1,307.23 |

3 | 1,307.134 | 1,308.479 | 1,307.18 |

4 | 1,306.927 | 1,308.476 | 1,307.11 |

5 | 1,306.928 | 1,308.473 | 1,307.07 |

6 | 1,306.928 | 1,308.471 | 1,307.04 |

7 | 1,307.168 | 1,308.466 | 1,306.97 |

8 | 1,307.166 | 1,308.458 | 1,306.94 |

9 | 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 | 1,307.145 | 1,308.44 | 1,307.29 |

2 | 1,307.134 | 1,308.474 | 1,307.23 |

3 | 1,307.134 | 1,308.479 | 1,307.18 |

4 | 1,306.927 | 1,308.476 | 1,307.11 |

5 | 1,306.928 | 1,308.473 | 1,307.07 |

6 | 1,306.928 | 1,308.471 | 1,307.04 |

7 | 1,307.168 | 1,308.466 | 1,306.97 |

8 | 1,307.166 | 1,308.458 | 1,306.94 |

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

*et al.*2016):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.

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