Numerical methods have been widely used to simulate transient groundwater flow induced by pumping wells in geometrically and mathematically complex systems. However, flow and transport simulation using low-order numerical methods can be computationally expensive with a low rate of convergence in multi-scale problems where fine spatial discretization is required to ensure stability and desirable accuracy (for instance, close to a pumping well). Numerical approaches based on high-order test functions may better emulate the global behavior of parabolic and/or elliptic groundwater governing equations with and without the presence of pumping well(s). Here, we assess the appropriateness of high-order differential quadrature method (DQM) and radial basis function (RBF)-DQM approaches compared to low-order finite difference and finite element methods. This assessment is carried out using the exact analytical solution by Theis and observed head data as benchmarks. Numerical results show that high-order DQM and RBF-DQM are more efficient schemes compared to low-order numerical methods in the simulation of 1-D axisymmetric transient flow induced by a pumping well. Mesh-less RBF-DQM, with the ability to implement arbitrary (e.g., adaptive) node distribution, properly simulates 2-D transient flow induced by pumping wells in confined/unconfined aquifers with regular and irregular geometries, compared to the other high-order and low-order approaches presented in this paper.

## INTRODUCTION

Numerical approaches have been considered as efficient tools in the simulation of mathematically and geometrically complex groundwater flow problems (e.g., Ameli *et al.* 2015). However, for the modeling of multi-scale problems, widely used low-order and mesh-based numerical techniques including finite difference (FDM) and finite element (FEM) methods may suffer from large computational demands, stability issues, and low convergence rates (Samani *et al.* 2004; Rushton 2007; An *et al.* 2010). In well test analyses with high gradients in the vicinity of pumping wells, the low-order nature of numerical schemes means accurate results can only be achieved by adopting numerous nodes leading to high computational costs and low convergence rates. For example, Haitjema *et al.* (2010) used a MODFLOW model of 1,846,314 cells to properly obtain the drawdown–discharge relationship of a single horizontal well in a small homogenous aquifer with a regular geometry (a domain of 150 m × 480 m × 24 m). The requirements for using mesh-based, low-order numerical approaches for the design of a group of pumping wells are even more demanding: in the design process, the assessment of different scenarios including various layouts and length of wells is required (Moore *et al.* 2012).

Numerical methods with high-order test functions, on the other hand, are far more consistent with the global behavior of parabolic or elliptic groundwater governing equations. High-order test functions can minimize the truncation error significantly (without having to employ many nodes) and can therefore increase the rate of convergence. For example, Moridis & Kansa (1992), in simulation of the Theis problem, reported a comparable accuracy between the results obtained via high-order multi-quadratic functions with only seven nodes and FDM with 60 blocks. In addition, Golberg & Chen (1997) achieved the same level of accuracy using a FEM approach with 71,000 elements and 60 nodes built from the high-order radial basis function (RBF)-collocation method in numerical simulation of the 3-D Poisson equation. High-order spectral numerical methods were also applied in groundwater simulation (Fagherazzi *et al.* 2004; Giudici & Vassena 2008). These methods were reported to have better convergence and less CPU time as compared to the FEM. More recently, high-order analytical approaches have also been used to address naturally and mathematically complex 2-D and 3-D steady-state groundwater flow problems (Craig 2009; Ameli *et al.* 2013; Ameli & Craig 2014).

The differential quadrature method (DQM) as a high-order numerical technique has also been reported as an efficient alternative to low-order methods in various branches of science (Malekzadeh *et al.* 2007; Hashemi *et al.* 2008). Many researchers have indicated that DQM is a highly accurate scheme with minimum computational effort and high rate of convergence in simulation of boundary-value (Bert & Malik 1996; Malekzadeh *et al.* 2007; Hashemi *et al.* 2008) and initial-value problems (Fung 2001a, 2001b). However, it is known that DQM is only appropriate for simulation of the problems with regular computational domains (Khoshfetrat & Abedini 2012, 2013; Homayoon *et al.* 2013). In addition, DQM solution remains stable only upon employing certain types of node distribution, such as Chebyshev–Gauss–Lobatto (Fung 2001a, 2001b; Zong *et al.* 2005; Atkinson 2008). Atkinson (2008) and Fung (2001a) argued that even if the nodes are uniformly distributed, the DQM solution may not be numerically stable. Of course, this restriction on the type of node distribution together with the limitation on the incorporation of irregular geometry weakens the application of DQM in groundwater modeling.

Recent advances in DQM have relaxed the constraints on geometry by augmenting the standard polynomial-based basis function of DQM with mesh-less RBF; this coupled scheme – called RBF-DQM – was recently used to emulate the behavior of shallow surface water and Burger's equation (Khoshfetrat & Abedini 2012, 2013; Homayoon *et al.* 2013). Note that mesh-less in the context of discrete numerical methods implies that the computational domain is discretized using arbitrary-placed ‘unstructured’ nodes rather than pre-defined structured mesh (e.g., triangular or quadrilateral).

In this study, application of the RBF-DQM is extended to a real transient groundwater problem with a naturally complex domain, non-linear mathematical model and various types of boundary conditions. For this case study, the efficiency of RBF-DQM as a high-order mesh-less approach is compared to the low-order mesh-based FEM and observed data. Furthermore, the authors hypothesize that mesh-less RBF-DQM, in addition to relaxing standard DQM constraints on geometry, can also address the aforementioned node distribution restriction of standard DQM. To analyse this hypothesis, we consider 1-D and 2-D hypothetical well test problems in confined aquifers with regular domains. RBF-DQM and standard DQM efficiency, in terms of computational cost, accuracy, and rate of convergence are examined and compared to low-order FEM and FDM. The analytical Theis solution is employed as benchmark to assess these models and numerical experiments.

## METHODS

### Transient flow in a confined aquifer (Theis problem)

#### Cartesian coordinate system

*h*is the hydraulic head (L),

*T*is the transmissivity (L

^{2}/T),

*S*is the specific storage (dimensionless), and

*x*and

*y*are coordinates (L). Here, the efficiency of high-order and low-order numerical methods are assessed compared to the Theis analytical solution (Theis 1935), which is developed for a confined aquifer with the following characteristics: (1) homogeneous, (2) isotropic aquifer with (3) infinite horizontal extent and (4) constant thickness subject to: (a) a single fully-penetrated pumping well with (b) a constant pumping rate equal to

*Q*, and (c) negligible well diameter relative to the aquifer's horizontal dimensions. In light of the above assumptions, Equation (1) is subject to the following boundary and initial conditions in Cartesian coordinate system: where and are pumping well coordinates in plan view, and are grid spacing in

*x*and

*y*directions, and

*h*is initial head condition.

_{0}#### Polar coordinate system

### Transient flow in an unconfined aquifer

*R*is recharge rate,

*K*is isotropic hydraulic conductivity, is specific yield, and is equal to pumping flux for the points which represent the pumping wells and zero for the remaining points inside the modeled domain. The boundaries of the naturally complex unconfined aquifer used in this paper (test case (3) include: a river with constant head; a no-flow condition along the sides of the domain; and a ditch with a constant flux). The latter condition can be described by the following non-linear equation obtained based on the Dupuit–Forchheimer approximation. where is a constant flux along the ditch.

### Differential quadrature method (DQM) and RBF-DQM

*x*,

*y*,

*z*] is the spatial location where the derivative is expanded around, is the number of nodes in the ξ direction,

*m*is derivative order,

*r*is the basis function used to represent function

*f*, and is a weighting coefficient. Quan & Chang (1989) used a Lagrange interpolation polynomial as a basis function: where and obtained unknown weighting coefficients as follows: By using the aforementioned basis function, DQM is found to be most efficient only when the following Chebyshev–Gauss–Lobatto node distribution is employed: As noted earlier, this issue can be mitigated by replacing standard DQM basis function (Equation (6)) with a mesh-less RBF: where function

*f*is described as and function derivatives are In the above equation

*X*

*=*(

*x*), is the center of the RBF and

_{1},x_{2}… .x_{d}*c*and

*q*are a priori unknown shape parameters which depend upon node configuration and the number of nodes used to discretize the computational domain (Kansa 1990; Kansa & Hon 2000; Wertz

*et al.*2006; Fornberg & Zuev 2007). Shape parameters also affect the shape of basis functions (Wang & Liu 2002); this implies that the optimum shape parameters might be dependent upon the shape of the response surface. Note that here the function

*f*is estimated using all

*N*nodes (used to discretize the modeled domain) and not only the nodes in the derivative direction, as is common in standard DQM. In a manner similar to Quan & Chang (1989), weighting coefficients in Equation (8c) for RBF-DQM with arbitrary node distribution can be obtained.

### Evaluation of model performance

The validity and applicability of various numerical techniques are evaluated using conventional statistical measures. Note that there are no appropriate guidelines for choosing the best model if different error measures support different models. In the current study, differences between numerically simulated results and exact analytical (or observed) ones are summarized using measures including percent accuracy ratio (PAR(k)) and percent average accuracy ratio (PAAR) (Samani *et al.* 2004) as follows.

*k*th grid point of spatial domain, respectively,

*N*is the total number of discretization points in space, and

*J*

_{ref}is the reference head (Samani

*et al.*2004). Note that PAR(k) is the error at the

*k*th point of the spatial domain at a specified time step and PAAR is the average of all PAR(k) values at that time step. Estimates are considered to be unbiased if either PAR(k) or PAAR is close to zero. The Nash–Sutcliffe (NS) value is also obtained as where

**is the average of analytical or observed values of hydraulic head in the entire spatial domain.**

## RESULTS AND DISCUSSION

High-order DQM and RBF-DQM are here assessed and compared to low order FDM and FEM in the simulation of groundwater flow induced by pumping well(s) in confined and unconfined aquifers with regular and irregular geometries. For this purpose, three test cases have been devised. In test cases 1 and 2, the 1-D and 2-D Theis problem were considered. In test case 3, we assessed the efficiency of mesh-less RBF-DQM compared to SEEP/W (FEM-FDM) in the simulation of transient flow in an unconfined aquifer with a naturally complex geometry and boundary conditions, and a non-linear mathematical model.

### 1-D Theis problem

In test case 1, DQM and RBF-DQM were used to simulate a 1-D Theis problem with a negligible pumping well radius. This hypothetical example was originally solved by Samani *et al.* (2004) using two modified versions of MODFLOW (FDM). Samani *et al.* (2004) enhanced standard MODFLOW to simulate radial flow near pumping well(s) and assessed the efficiency of this modified version of MODFLOW, termed LSM, in the simulation of the Theis problem by comparing LSM results to another MODFLOW modification by Barrash & Dougherty (1997). In the current study, the efficiency of DQM with a Chebyshev–Gauss–Lobatto node distribution and mesh-less RBF-DQM with arbitrary node distribution were assessed and compared to FDM results presented by Samani *et al.* (2004). Figure 1(a) shows both RBF-DQM and standard DQM node distributions in spatial domain discretization. For adaptive node distribution used by the two modified versions of MODFLOW, readers are referred to Samani *et al.* (2004). It should be noted that in both RBF-DQM and standard DQM approaches used here the ‘temporal’ domain are discretized using standard DQM with a Chebyshev–Gauss–Lobatto node distribution (Equation (7)). To find RBF shape parameters, PEST (non-linear parameter optimization) software was used where the Theis analytical solution was considered the true solution.

Material properties and other groundwater flow parameters used by Samani *et al.* (2004) are: *T**=* 8 × 10^{−}^{5} m^{2}/s, *S**=* 8 × 10^{−}^{3}, *h*_{0} = 100 m, *Q**=* 6*.*28 × 10^{−}^{4} m^{3}/s, simulation period = 19,943 s, = 0*.*001 m, and *J*_{ref} = 12*.*5 m. RBF shape parameters for this example were found as *c**=* 0*.*0055 and *q**=* 1*.*57. Figure 1(b) shows the variation of PAR versus radial distance for both DQM and RBF-DQM at the end of the simulation period. The PAR of RBF-DQM in the simulation of head near the pumping well point is 0.05%, which is better than DQM (0.32%). It appears that augmenting conventional DQM with RBF removes the constraints on the application of adaptive node distribution within DQM, and improves the accuracy and applicability of the method for this multi-scale problem. The associated figures obtained for Barrash & Dougherty (1997) and LSM approaches reported by Samani *et al.* (2004) are 0.07% and 1%, respectively. With respect to PAAR criteria, Table 1 shows the superior performance of RBF-DQM compared to the LSM, DQM, and Barrash & Dougherty (1997) approaches in the whole computational domain and at the end of the simulation period. Note that for LSM and Barrash & Dougherty (1997) approaches, the number of spatial nodes (*N _{r}*) mentioned in Table 1 is equal to number of columns × number of rows × number of layers in their MODFLOW simulation (see Samani

*et al.*2004). Table 1 also indicates that DQM generally outperforms the LSM and Barrash and Dougherty approaches in the whole computational domain while using fewer nodes.

Method | N _{r} | N _{t} | PAAR (%) |
---|---|---|---|

Samani et al. (2004) (LSM) | 60 × 1 × 43 | 449 | 0.08 |

Barrash & Dougherty (1997) | 41 × 81 × 43 | 449 | 3.7 |

DQM | 30 | 6 | 2.68 |

DQM | 140 | 6 | 0.09 |

DQM | 160 | 6 | 0.07 |

DQM | 220 | 6 | 0.03 |

RBF-DQM | 140 | 6 | 0.02 |

Method | N _{r} | N _{t} | PAAR (%) |
---|---|---|---|

Samani et al. (2004) (LSM) | 60 × 1 × 43 | 449 | 0.08 |

Barrash & Dougherty (1997) | 41 × 81 × 43 | 449 | 3.7 |

DQM | 30 | 6 | 2.68 |

DQM | 140 | 6 | 0.09 |

DQM | 160 | 6 | 0.07 |

DQM | 220 | 6 | 0.03 |

RBF-DQM | 140 | 6 | 0.02 |

*N _{r}* and

*N*are the number of nodes in spatial and temporal domain discretization.

_{t}The optimized RBF's shape parameters were found to be highly dependent on node configuration and number of nodes as well as the shape of response function (Figure 2); the latter received scant attention in the relevant literature. Furthermore, the majority of research studies have considered *c* as the only shape parameter with a fixed *q* equal to 0.5 (Kansa 1990; Kansa & Hon 2000; Fornberg & Zuev 2007). Here, we consider the effect of response surface shape on optimized shape parameters *c* and *q* (Figure 2(a)) and also assess the role of *q* in the efficiency of the RBF-DQM approach (Figure 2(b)). As pumping well radius would affect the shape of response surface, we investigated the impact of well radius on optimized shape parameters. Upon increasing the well radius, the gradient near the well decreases significantly and the steep response surface in the vicinity of the pumping well changes to a milder one. The number of nodes (140 nodes) and node configuration (ascending distance) were assumed to be constant for all well radii to ensure an appropriate assessment of the effect of function response surface on shape parameters. Figure 2(a) depicts, upon increasing well radius from 0.001 m to 0.15 m, the optimized *c* and *q* values are leveled off after a sharp increase and decrease, respectively. This figure implies that as long as surface response varies very rapidly with well radius (when dealing with a very small range of radii), the optimized shape parameters change. However, once the response surface is stabilized for a larger pumping well radius, shape parameters do not show considerable changes. This behavior implies that although the optimized shape parameters are dependent upon response surface shape (here well radius), for a realistically large pumping well radius the optimized shape parameters are independent of pumping well radius. To assess the role of *q* parameter in RBF-DQM efficiency, we obtained the variation of PAAR with *q* for different values of the *c* parameter (Figure 2(b)). As can be seen, the widely used value of *q**=* 0*.*5 provides almost the lowest level of accuracy for various *c* values. Figure 2(b) also indicates that the sensitivity of PAAR to parameter *c* decreases significantly as *q* increases. It should be noted that the optimum range for *q* (0.98 to 1.03) reported by Wang & Liu (2002) in 2-D interpolation of the response surface using a RBF is not an optimal range for 1-D simulation of a Theis problem.

### 2-D Theis problem in Cartesian coordinates

The appropriateness of DQM and RBF-DQM is here assessed and compared to SEEP/W (FEM-FDM) and standard MODFLOW (FDM) in 2-D simulation of groundwater flow induced by a single pumping well (i.e., a 2-D representation of the Theis problem). The pumping well is located in the center of a 4,000 m × 4,000 m isotropic and homogeneous confined aquifer with a discharge of 4,000 m^{3}/day. Aquifer material properties are *T**=* 300 m^{2}/day and = 0.002. The pumping period considered for this example is 13.20 days with an assumed initial head of 10 m at the start of pumping. In a manner similar to Wang & Anderson (1995) and because of the symmetry of the flow around a single pumping well with a radius of 0.5 m, only one quarter of the domain with a pumping well at (*x _{p}* = 0,

*y*= 0) and two impermeable sides was considered (Figure 3). In addition, it was assumed that hydraulic head far away from the pumping well (i.e., outside the radius of influence) is equal to 10 m. DQM was implemented using the Chebyshev–Gauss–Lobatto node distribution in both directions while RBF-DQM was applied with an ascending distance node distribution from the pumping well (Figure 3(a)). In two separate simulations, triangular elements including 210 and 5,100 nodes with a size ratio of 6 were used within SEEP/W; Figure 3(b) shows the latter case with almost 9,000 triangular elements. The large size ratio used to develop the mesh of these models ensures proper incorporation of the pumping well where the elements in the vicinity of the pumping well (circle in Figure 3(b)) are considerably finer than those located farther away from the well. FDM blocks including 210 and 900 nodes were also implemented in MODFLOW software where both models were sufficiently refined in the vicinity of the pumping well.

_{p}The first four rows of Table 2 show that, using an identical number of nodes, RBF-DQM is more accurate than the remaining high-order and low-order approaches presented in this paper. Table 2 also depicts RBF-DQM with a coarse node distribution (210 nodes) works almost similar to MODFLOW (900 nodes) and better than SEEP/W with a very fine mesh (5,100 nodes). Additionally, Figure 4(a) suggests that RBF-DQM solution converges significantly faster than other numerical schemes. The optimized values for RBF shape parameters were obtained as 3 and 0.98 for *c* and *q*, respectively. An effort was also made to evaluate the feasibility of adapting the optimum range of *q* reported by Wang & Liu (2002) in 2-D interpolation of the response surface using a RBF. Figure 4(b) demonstrates the variation of PAAR with *q* for different values of *c*. PAAR for different values of *c* is within an acceptable range and the sensitivity of solution to the *c* parameter is small. This implies that the cited range for *q* (0.98 to1.03 excluding *q**=* 1) is optimum for the simulation of a 2-D Theis problem using the RBF-DQM numerical scheme. In recent years, RBF-DQM has been criticized for its dependency on shape parameters. This evaluation gives another advantage to RBF-DQM where sensitivity of the solution to *c* parameter is considerably reduced if the cited range of *q* is implemented.

Method | N _{r} | N _{t} | CPU time (sec) | PAAR (%) |
---|---|---|---|---|

MODFLOW | 210 | 13 | 1 | 0.52 |

SEEP/W | 210 | 13 | 7 | 1.86 |

DQM | 210 | 4 | 0.27 | 0.30 |

RBF-DQM | 210 | 4 | 0.62 | 0.16 |

SEEP/W | 5100 | 13 | 30 | 0.37 |

MODFLOW | 900 | 13 | 3 | 0.37 |

Method | N _{r} | N _{t} | CPU time (sec) | PAAR (%) |
---|---|---|---|---|

MODFLOW | 210 | 13 | 1 | 0.52 |

SEEP/W | 210 | 13 | 7 | 1.86 |

DQM | 210 | 4 | 0.27 | 0.30 |

RBF-DQM | 210 | 4 | 0.62 | 0.16 |

SEEP/W | 5100 | 13 | 30 | 0.37 |

MODFLOW | 900 | 13 | 3 | 0.37 |

Although based on the results in this example DQM generally outperforms SEEP/W and MODFLOW, the accuracy of 2-D DQM solution is less than the 1-D solution developed in the first example. This can be attributed to the fact that the restriction of DQM on the implementation of only specific node distributions is more problematic in 2-D simulation rather than 1-D one. This may explain why recent application of DQM in groundwater hydrology has been limited to 1-D well-behaved problems (Ghaheri & Meraji 2011; Kaya & Arisoy 2011). Augmenting standard DQM with RBFs removes the DQM constraint on employing only certain types of node distribution and improves its efficiency in the simulation of multi-scale 1-D, 2-D, and even 3-D groundwater problems.

### Real case study

In this example, the numerical behavior of RBF-DQM and SEEP/W (FEM-FDM) were assessed for the simulation of transient groundwater flow in a naturally complex heterogeneous unconfined aquifer (Figure 5). Observed head data collected at the monitoring points shown in Figure 5(a) (in alphabetical order) were used as a benchmark to calibrate and validate the developed models, and assess the efficiency of both approaches. The modelled domain encompasses a 1,500 m × 1,500 m sand and gravel aquifer in the presence of clay and silt lenses (Figure 5(b)). Figure 5 also depicts the boundaries of the unconfined aquifer. The northern boundary is formed by a 100 m wide river with constant head (heads are shown in Figure 5(a)) and irregular shape. This river is underlain by a mixture of sand and gravel of an approximately 30 m/day vertical hydraulic conductivity. The southern boundary is a gravel-bottomed ditch with a 100 m width and a constant flux of = 45,000 m^{3}/day. No-flow conditions represent bedrock along the sides (east and west) of the modeled domain. The entire domain is subject to a constant recharge of *R**=* 0.0001 m/day. Borehole tests at points A and M estimated a hydraulic conductivity of *K**=* 75 m/day at these two points. Specific yield for this sand and gravel aquifer was estimated as 0.10. Table 3 presents steady-state hydraulic heads collected along the entire aquifer at monitoring points shown in Figure 5(a) (column 7). Observed heads after 3 days of pumping at point A (fully penetrating pumping well) with a pumping rate of = 20,000 m^{3}/day are also presented in Table 3 (column 2).

Monitoring points | O-T | S-T RBF-DQM 100 × 100 | S-T SEEP/W 100 × 100 | S-T SEEP/W 50 × 50 | S-T SEEP/W 25 × 25 | O-SS | S-SS RBF-DQM 100 × 100 | S-SS SEEP/W 100 × 100 | S-SS SEEP/W 50 × 50 | S-SS SEEP/W 25 × 25 |
---|---|---|---|---|---|---|---|---|---|---|

S | 521.95 | 523.23 | 522.35 | 522.8 | 522.95 | 521.96 | 523.90 | 522.40 | 522.85 | 522.95 |

H | 519.55 | 520.55 | 519.40 | 520.00 | 520.05 | 519.70 | 521.24 | 519.85 | 520.05 | 520.10 |

O | 518.86 | 518.50 | 517.05 | 518.25 | 518.35 | 519.02 | 519.55 | 517.70 | 518.30 | 518.45 |

I | 519.25 | 518.60 | 517.50 | 518.70 | 519.00 | 519.28 | 518.64 | 517.55 | 518.80 | 519.10 |

D | 516.17 | 516.17 | 515.95 | 516.70 | 517.00 | 516.71 | 517.15 | 515.95 | 516.80 | 517.05 |

C | 515.66 | 515.71 | 515.65 | 516.50 | 517.00 | 516.03 | 516.07 | 515.75 | 516.55 | 517.05 |

Q | 518.18 | 517.63 | 515.65 | 516.60 | 517.00 | 518.32 | 517.65 | 515.70 | 516.65 | 517.05 |

M | 516.68 | 516.03 | 514.95 | 515.75 | 516.00 | 517.12 | 516.74 | 515.10 | 515.75 | 516.50 |

K | 512.21 | 512.37 | 513.65 | 514.40 | 515.00 | 513.88 | 513.44 | 513.70 | 514.40 | 514.65 |

J | 515.71 | 514.93 | 513.85 | 514.45 | 514.75 | 515.71 | 514.87 | 513.85 | 514.65 | 515.00 |

E | 513.04 | 512.29 | 512.75 | 513.45 | 514.00 | 513.17 | 512.72 | 512.85 | 513.50 | 514.15 |

A | 508.8 | 509.90 | 511.10 | 509.90 | 509.00 | 512.22 | 511.77 | 512.55 | 513.30 | 513.50 |

B | 511.29 | 510.93 | 512.30 | 513.00 | 513.15 | 511.95 | 511.38 | 512.40 | 513.20 | 513.30 |

N | 512.83 | 512.56 | 511.95 | 512.70 | 513.00 | 512.83 | 512.68 | 512.05 | 512.70 | 513.00 |

G | 507.99 | 509.38 | 510.30 | 511.00 | 511.00 | 508.19 | 509.57 | 510.35 | 511.10 | 511.15 |

F | 507.79 | 507.64 | 509.95 | 510.70 | 510.80 | 508.71 | 508.05 | 510.05 | 510.80 | 510.90 |

P | 509.11 | 509.80 | 510.05 | 510.40 | 510.50 | 509.12 | 509.81 | 510.05 | 510.40 | 510.50 |

Monitoring points | O-T | S-T RBF-DQM 100 × 100 | S-T SEEP/W 100 × 100 | S-T SEEP/W 50 × 50 | S-T SEEP/W 25 × 25 | O-SS | S-SS RBF-DQM 100 × 100 | S-SS SEEP/W 100 × 100 | S-SS SEEP/W 50 × 50 | S-SS SEEP/W 25 × 25 |
---|---|---|---|---|---|---|---|---|---|---|

S | 521.95 | 523.23 | 522.35 | 522.8 | 522.95 | 521.96 | 523.90 | 522.40 | 522.85 | 522.95 |

H | 519.55 | 520.55 | 519.40 | 520.00 | 520.05 | 519.70 | 521.24 | 519.85 | 520.05 | 520.10 |

O | 518.86 | 518.50 | 517.05 | 518.25 | 518.35 | 519.02 | 519.55 | 517.70 | 518.30 | 518.45 |

I | 519.25 | 518.60 | 517.50 | 518.70 | 519.00 | 519.28 | 518.64 | 517.55 | 518.80 | 519.10 |

D | 516.17 | 516.17 | 515.95 | 516.70 | 517.00 | 516.71 | 517.15 | 515.95 | 516.80 | 517.05 |

C | 515.66 | 515.71 | 515.65 | 516.50 | 517.00 | 516.03 | 516.07 | 515.75 | 516.55 | 517.05 |

Q | 518.18 | 517.63 | 515.65 | 516.60 | 517.00 | 518.32 | 517.65 | 515.70 | 516.65 | 517.05 |

M | 516.68 | 516.03 | 514.95 | 515.75 | 516.00 | 517.12 | 516.74 | 515.10 | 515.75 | 516.50 |

K | 512.21 | 512.37 | 513.65 | 514.40 | 515.00 | 513.88 | 513.44 | 513.70 | 514.40 | 514.65 |

J | 515.71 | 514.93 | 513.85 | 514.45 | 514.75 | 515.71 | 514.87 | 513.85 | 514.65 | 515.00 |

E | 513.04 | 512.29 | 512.75 | 513.45 | 514.00 | 513.17 | 512.72 | 512.85 | 513.50 | 514.15 |

A | 508.8 | 509.90 | 511.10 | 509.90 | 509.00 | 512.22 | 511.77 | 512.55 | 513.30 | 513.50 |

B | 511.29 | 510.93 | 512.30 | 513.00 | 513.15 | 511.95 | 511.38 | 512.40 | 513.20 | 513.30 |

N | 512.83 | 512.56 | 511.95 | 512.70 | 513.00 | 512.83 | 512.68 | 512.05 | 512.70 | 513.00 |

G | 507.99 | 509.38 | 510.30 | 511.00 | 511.00 | 508.19 | 509.57 | 510.35 | 511.10 | 511.15 |

F | 507.79 | 507.64 | 509.95 | 510.70 | 510.80 | 508.71 | 508.05 | 510.05 | 510.80 | 510.90 |

P | 509.11 | 509.80 | 510.05 | 510.40 | 510.50 | 509.12 | 509.81 | 510.05 | 510.40 | 510.50 |

(SS) refers to steady-state heads and (T) refers to transient heads after 3 days of pumping at point A.

Here we calibrated RBF-DQM and SEEP/W models to collected head data obtained after 3 days of pumping at point A. Observed steady-state head data were used to validate the developed models. The system of non-linear equations obtained based on Equation (4) was solved using a trust-region-reflective algorithm to simulate groundwater flow in an unconfined aquifer based on the Dupuit–Forchheimer approximation and boundary conditions noted above. RBF-DQM node distribution consists of a uniform node spacing of 100 m in both *x* and *y* directions; these nodes were located at the center of each cell shown in Figure 5(a). The triangular mesh used in the original SEEP/W model was also constructed with a uniform node spacing of 100 m in both *x* and *y* directions. The modeled domain was divided into seven hydrogeological zones (Figure 5(a)) where for each zone a homogenous isotropic hydraulic conductivity and a single specific yield parameter for the entire domain were obtained through calibration. In the RBF-DQM solution, shape parameters (i.e., *c* and *q*) in addition to the material parameters were also obtained in the calibration processes using the PEST software (Table 4). Again here the optimum shape parameter *q* is in the range of *q* reported by Wang & Liu (2002).

Parameters | RBF-DQM | SEEP/W |
---|---|---|

K_{1} | 213.07 | 195 |

K_{2} | 254.82 | 260 |

K_{3} | 38.45 | 35 |

K_{4} | 70.67 | 70 |

K_{5} | 30.00 | 35 |

K_{6} | 69.79 | 70 |

K_{7} | 40.16 | 40 |

S_{y} | 0.090 | 0.085 |

c | 0.005 | – |

q | 1.007 | – |

Parameters | RBF-DQM | SEEP/W |
---|---|---|

K_{1} | 213.07 | 195 |

K_{2} | 254.82 | 260 |

K_{3} | 38.45 | 35 |

K_{4} | 70.67 | 70 |

K_{5} | 30.00 | 35 |

K_{6} | 69.79 | 70 |

K_{7} | 40.16 | 40 |

S_{y} | 0.090 | 0.085 |

c | 0.005 | – |

q | 1.007 | – |

Table 3 presents the simulated heads after 3 days of pumping at A using RBF-DQM and SEEP/W models developed here (columns 3–6). In addition to the original node spacing (100 m × 100 m) used in the calibration phase within the SEEP/W model, simulated results using finer node spacing (i.e., 50 m × 50 m and 25 m × 25 m) were also obtained without further calibration to assess the effect of mesh refinement on SEEP/W model efficiency. Figure 6(a) compares the accuracy of RBF-DQM and SEEP/W in the calibration of simulated results to the observed head data after 3 days of pumping at point A. This figure shows that RBF-DQM using 100 m × 100 m node spacing outperforms SEEP/W in most monitoring points. Using RBF-DQM, PAAR and Nash–Sutcliffe values obtained based on all monitoring points are 0.11% and 0.97, respectively. For the SEEP/W model, PAAR values using 100 m × 100 m, 50 m × 50 m and 25 m × 25 m node spacing are 0.25%, 0.23%, and 0.22%, respectively; Nash–Sutcliffe values for the corresponding node spacing are 0.87, 0.88, and 0.88, respectively. Note that to obtain PAR and PAAR in this example, *J _{ref}* was considered as the average of observed head in the entire domain (514 m). The validation of the model was performed with respect to the observed steady-state head data. Simulated heads and PAR values at each monitoring points are shown in Table 3 (columns 8–11) and Figure 6(b), respectively. PAAR and Nash–Sutcliffe values, using RBF-DQM, are 0.14% and 0.95, respectively. Using SEEP/W model, PAAR values obtained based on 100 m × 100 m, 50 m × 50 m and 25 m × 25 m node spacing are 0.20%, 0.19%, and 0.19%, respectively; Nash–Sutcliffe values for the corresponding node spacing are 0.90, 0.91, and 0.91.

This example showed that RBF-DQM can properly simulate transient groundwater flow in the vicinity of a pumping well in an unconfined aquifer with a naturally complex geometry and non-linear governing equation and boundary conditions. RBF-DQM outperformed the SEEP/W (FEM-FDM) model using an identical node distribution in both calibration and validation phases. Optimum *q* value of the RBF-DQM approach was obtained within the range suggested by Wang & Liu (2002) similar to test case 2.

## CONCLUSIONS

Results of the current paper indicate that a mesh-less RBF-DQM approach outperforms SEEP/W (FEM-FDM) model in the simulation of transient groundwater flow near a vertical pumping well in an unconfined aquifer with naturally complex geometry and spatially varied material properties under non-linear Dupuit–Forchheimer condition. RBF-DQM shows also a superior performance compared to FEM and FDM approaches in the simulation of 1-D and 2-D transient groundwater flow near a vertical pumping well in confined aquifers. Based on high-order test functions used in the RBF-DQM algorithm, all nodes in the computational domain have direct impact on the estimation of derivatives and state variables at every node; this can aid in the rapid transfer of the effect of a pumping well and/or boundary conditions to the entire computational domain. Therefore, accurate results with rapid rate of convergence through a few nodes can be obtained using RBF-DQM compared to FEM and FDM methods. Similar to RBF-DQM, a high-order standard DQM approach can simulate groundwater flow with a high level of accuracy and rate of convergence; albeit the application of standard DQM is limited to 1-D well-behaved groundwater problems. However, the performance of DQM in temporal discretization is much better than traditional methods like FDM without any fear of instability and without limitations on the time step.

This paper also provides some insights into choosing RBF shape parameters used in the RBF-DQM algorithm. Results indicated that considering *q* as a variable shape parameter, instead of assigning a fixed value of 0.5, can increase the rate of convergence and level of accuracy for RBF-DQM, since the combination of *c* and *q* can construct the function response surface more accurately than only tuning *c*. Furthermore, in the presence of a large *q* (larger than standard *q* = 0.5)*,* the sensitivity of objective functions to *c* can be significantly reduced. This research also showed that the suitable range of *q* reported by Wang & Liu (2002) for 2-D point interpolation using RBF is an efficient range for solving the 2-D non-linear groundwater flow equation using RBF-DQM. Proposing a methodology to avoid optimization of shape parameters is underway and will be published in due course.

## ACKNOWLEDGEMENTS

Mark Ranjram is thanked for his comments that improved the readability of the manuscript. Matlab simulation codes, SEEP/W and MODFLOW data are available upon request from the corresponding author.