Development of two-dimensional groundwater flow simulation model using meshless method based on MLS approximation function in unconfined aquifer in transient state

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.

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

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.

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

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.

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.

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.

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

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

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.

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.

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.