Application of Meshless local Petrov-Galerkin approach for steady state groundwater flow modeling

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.

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.

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 (m³/day/m²), 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.

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.

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

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.

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

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.

A simple rectangular unconfined aquifer is presented in Figure 2. The size of the aquifer is 15 × 15 m². 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.

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.

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.

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.

A standard unconfined aquifer with the area of 4,500 × 10,000 m² 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.

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.

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.

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.

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.

Birjand unconfined aquifer is located in South Khorasan province in Iran. The area of this aquifer is approximately 269 km². 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.

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.

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

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

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

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.

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.

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

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.