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

The 2-D transient equation which describes groundwater movement in an isotropic and homogeneous confined aquifer is as follows: 
formula
1
where h is the hydraulic head (L), T is the transmissivity (L2/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: 
formula
2a
 
formula
2b
 
formula
2c
where and are pumping well coordinates in plan view, and are grid spacing in x and y directions, and h0 is initial head condition.

Polar coordinate system

Equation (1) can be represented in 1-D radial coordinate system by ignoring the directional variation of hydraulic head (due to the symmetric radial flow through a single well) as follows: 
formula
3a
Similarly, boundary conditions along the well screen and far away from the well screen are as follows: 
formula
3b
 
formula
3c
 
formula
3d

Transient flow in an unconfined aquifer

The transient non-linear partial differential equation which describes groundwater movement in a homogenous and isotropic unconfined aquifer with a priori unknown water table location (variable aquifer thickness), can be obtained by introducing Darcy's equation into the continuity equation and employing the Dupuit–Forchheimer approximation as follows: 
formula
4a
where 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. 
formula
4b
where is a constant flux along the ditch.

Differential quadrature method (DQM) and RBF-DQM

The basic idea behind DQM is the approximation of function derivatives by using weighted sums of function values: 
formula
5a
where
Here, ξ = [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: 
formula
6a
where  
formula
and obtained unknown weighting coefficients as follows: 
formula
6b
 
formula
6c
By using the aforementioned basis function, DQM is found to be most efficient only when the following Chebyshev–Gauss–Lobatto node distribution is employed: 
formula
7
As noted earlier, this issue can be mitigated by replacing standard DQM basis function (Equation (6)) with a mesh-less RBF: 
formula
8a
where function f is described as 
formula
8b
and function derivatives are 
formula
8c
In the above equation X= (x1 ,x2… .xd), is the center of the RBF and 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.

For given t 
formula
9
 
formula
10
where and are the simulated and Theis analytical (or observed) values of hydraulic head in the kth grid point of spatial domain, respectively, N is the total number of discretization points in space, and Jref is the reference head (Samani et al. 2004). Note that PAR(k) is the error at the kth 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 
formula
11
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.

Figure 1

Test case 1. (a) DQM node distribution based on Chebyshev–Gauss–Lobatto type (Equation (7)) using Nr = 220 nodes and ascending distance node distribution of RBF-DQM using Nr = 140 nodes, black circle shows the pumping well with a radius of 0.001 m located at r = 0. (b) Percent accuracy ratio versus radial distance using both methods at the end of the simulation period, values in parentheses depict the number of nodes in the spatial and temporal domain discretization, respectively.

Figure 1

Test case 1. (a) DQM node distribution based on Chebyshev–Gauss–Lobatto type (Equation (7)) using Nr = 220 nodes and ascending distance node distribution of RBF-DQM using Nr = 140 nodes, black circle shows the pumping well with a radius of 0.001 m located at r = 0. (b) Percent accuracy ratio versus radial distance using both methods at the end of the simulation period, values in parentheses depict the number of nodes in the spatial and temporal domain discretization, respectively.

Material properties and other groundwater flow parameters used by Samani et al. (2004) are: T= 8 × 105 m2/s, S= 8 × 103, h0 = 100 m, Q= 6.28 × 104 m3/s, simulation period = 19,943 s, = 0.001 m, and Jref = 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 (Nr) 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.

Table 1

PAAR performance statistics over the entire spatial domain for test case 1 (1-D simulation of Theis problem) using different approaches at the end of the simulation period

Method Nr Nt PAAR (%) 
Samani et al. (2004) (LSM) 60 × 1 × 43 449 0.08 
Barrash & Dougherty (1997)  41 × 81 × 43 449 3.7 
DQM 30 2.68 
DQM 140 0.09 
DQM 160 0.07 
DQM 220 0.03 
RBF-DQM 140 0.02 
Method Nr Nt PAAR (%) 
Samani et al. (2004) (LSM) 60 × 1 × 43 449 0.08 
Barrash & Dougherty (1997)  41 × 81 × 43 449 3.7 
DQM 30 2.68 
DQM 140 0.09 
DQM 160 0.07 
DQM 220 0.03 
RBF-DQM 140 0.02 

Nr and Nt are the number of nodes in spatial and temporal domain discretization.

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.

Figure 2

The relationship between RBF shape parameters, shape of response surface, and the accuracy of RBF-DQM. (a) Effect of pumping well radius on optimized RBF shape parameters (c and q). (b) Variation of PAAR with q at the end of simulation period for various c parameter values.

Figure 2

The relationship between RBF shape parameters, shape of response surface, and the accuracy of RBF-DQM. (a) Effect of pumping well radius on optimized RBF shape parameters (c and q). (b) Variation of PAAR with q at the end of simulation period for various c parameter values.

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 m3/day. Aquifer material properties are T= 300 m2/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 (xp = 0, yp = 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.

Figure 3

Mesh and nodes topology used to develop models for test case 2. (a) DQM Chebyshev–Gauss–Lobatto node distribution (gridlines) and RBF-DQM ascending distance node distribution and (b) SEEP/W mesh topology with 9,000 triangular elements (5,100 nodes). Bold black lines show the impermeable sides and circles show the location of the pumping well.

Figure 3

Mesh and nodes topology used to develop models for test case 2. (a) DQM Chebyshev–Gauss–Lobatto node distribution (gridlines) and RBF-DQM ascending distance node distribution and (b) SEEP/W mesh topology with 9,000 triangular elements (5,100 nodes). Bold black lines show the impermeable sides and circles show the location of the pumping well.

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.

Table 2

Performance indices of test case 2 for various approaches in 2-D simulation of Theis problem at the end of the simulation period

Method Nr Nt CPU time (sec) PAAR (%) 
MODFLOW 210 13 0.52 
SEEP/W 210 13 1.86 
DQM 210 0.27 0.30 
RBF-DQM 210 0.62 0.16 
SEEP/W 5100 13 30 0.37 
MODFLOW 900 13 0.37 
Method Nr Nt CPU time (sec) PAAR (%) 
MODFLOW 210 13 0.52 
SEEP/W 210 13 1.86 
DQM 210 0.27 0.30 
RBF-DQM 210 0.62 0.16 
SEEP/W 5100 13 30 0.37 
MODFLOW 900 13 0.37 
Figure 4

Performance of the numerical methods in the 2-D simulation of the Theis problem. (a) Convergence behavior of various numerical schemes where N refers to the number of nodes. (b) Variation of PAAR with q at the end of simulation period for various values of c parameter in the RBF-DQM solution.

Figure 4

Performance of the numerical methods in the 2-D simulation of the Theis problem. (a) Convergence behavior of various numerical schemes where N refers to the number of nodes. (b) Variation of PAAR with q at the end of simulation period for various values of c parameter in the RBF-DQM solution.

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 m3/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 m3/day are also presented in Table 3 (column 2).

Table 3

Observed (O) and simulated (S) heads using RBF-DQM and SEEP/W (with different mesh sizes) at monitoring points shown in Figure 5(a)

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 
521.95 523.23 522.35 522.8 522.95 521.96 523.90 522.40 522.85 522.95 
519.55 520.55 519.40 520.00 520.05 519.70 521.24 519.85 520.05 520.10 
518.86 518.50 517.05 518.25 518.35 519.02 519.55 517.70 518.30 518.45 
519.25 518.60 517.50 518.70 519.00 519.28 518.64 517.55 518.80 519.10 
516.17 516.17 515.95 516.70 517.00 516.71 517.15 515.95 516.80 517.05 
515.66 515.71 515.65 516.50 517.00 516.03 516.07 515.75 516.55 517.05 
518.18 517.63 515.65 516.60 517.00 518.32 517.65 515.70 516.65 517.05 
516.68 516.03 514.95 515.75 516.00 517.12 516.74 515.10 515.75 516.50 
512.21 512.37 513.65 514.40 515.00 513.88 513.44 513.70 514.40 514.65 
515.71 514.93 513.85 514.45 514.75 515.71 514.87 513.85 514.65 515.00 
513.04 512.29 512.75 513.45 514.00 513.17 512.72 512.85 513.50 514.15 
508.8 509.90 511.10 509.90 509.00 512.22 511.77 512.55 513.30 513.50 
511.29 510.93 512.30 513.00 513.15 511.95 511.38 512.40 513.20 513.30 
512.83 512.56 511.95 512.70 513.00 512.83 512.68 512.05 512.70 513.00 
507.99 509.38 510.30 511.00 511.00 508.19 509.57 510.35 511.10 511.15 
507.79 507.64 509.95 510.70 510.80 508.71 508.05 510.05 510.80 510.90 
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 
521.95 523.23 522.35 522.8 522.95 521.96 523.90 522.40 522.85 522.95 
519.55 520.55 519.40 520.00 520.05 519.70 521.24 519.85 520.05 520.10 
518.86 518.50 517.05 518.25 518.35 519.02 519.55 517.70 518.30 518.45 
519.25 518.60 517.50 518.70 519.00 519.28 518.64 517.55 518.80 519.10 
516.17 516.17 515.95 516.70 517.00 516.71 517.15 515.95 516.80 517.05 
515.66 515.71 515.65 516.50 517.00 516.03 516.07 515.75 516.55 517.05 
518.18 517.63 515.65 516.60 517.00 518.32 517.65 515.70 516.65 517.05 
516.68 516.03 514.95 515.75 516.00 517.12 516.74 515.10 515.75 516.50 
512.21 512.37 513.65 514.40 515.00 513.88 513.44 513.70 514.40 514.65 
515.71 514.93 513.85 514.45 514.75 515.71 514.87 513.85 514.65 515.00 
513.04 512.29 512.75 513.45 514.00 513.17 512.72 512.85 513.50 514.15 
508.8 509.90 511.10 509.90 509.00 512.22 511.77 512.55 513.30 513.50 
511.29 510.93 512.30 513.00 513.15 511.95 511.38 512.40 513.20 513.30 
512.83 512.56 511.95 512.70 513.00 512.83 512.68 512.05 512.70 513.00 
507.99 509.38 510.30 511.00 511.00 508.19 509.57 510.35 511.10 511.15 
507.79 507.64 509.95 510.70 510.80 508.71 508.05 510.05 510.80 510.90 
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.

Figure 5

Layout of the unconfined aquifer used in test case 3. (a) Plan view along with seven hydrogeological zones used in the model development phase and the monitoring points shown in alphabetical order. Constant heads in the river are also shown at the north of the domain. (b) Cross section of the unconfined aquifer at x= 750 m. Figure was modified from Anderson & Woessner (1992).

Figure 5

Layout of the unconfined aquifer used in test case 3. (a) Plan view along with seven hydrogeological zones used in the model development phase and the monitoring points shown in alphabetical order. Constant heads in the river are also shown at the north of the domain. (b) Cross section of the unconfined aquifer at x= 750 m. Figure was modified from Anderson & Woessner (1992).

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

Table 4

Calibrated material and shape parameters for RBF-DQM and SEEP/W models (test case 3)

Parameters RBF-DQM SEEP/W 
K1 213.07 195 
K2 254.82 260 
K3 38.45 35 
K4 70.67 70 
K5 30.00 35 
K6 69.79 70 
K7 40.16 40 
Sy 0.090 0.085 
c 0.005 – 
q 1.007 – 
Parameters RBF-DQM SEEP/W 
K1 213.07 195 
K2 254.82 260 
K3 38.45 35 
K4 70.67 70 
K5 30.00 35 
K6 69.79 70 
K7 40.16 40 
Sy 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, Jref 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.

Figure 6

Comparison between the accuracy of RBF-DQM and SEEP/W (with different mesh sizes) for the simulation of groundwater flow in an unconfined aquifer in test case 3. (a) Calibration and (b) validation phases.

Figure 6

Comparison between the accuracy of RBF-DQM and SEEP/W (with different mesh sizes) for the simulation of groundwater flow in an unconfined aquifer in test case 3. (a) Calibration and (b) validation phases.

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.

REFERENCES

REFERENCES
An
H.
Ichikawa
Y.
Tachikawa
Y.
Shiiba
M.
2010
Three-dimensional finite difference saturated-unsaturated flow modeling with nonorthogonal grids using a coordinate transformation method
.
Water Resour. Res.
46
(
11
), W11521. doi:10.1029/2009wr009024.
Anderson
M. P.
Woessner
W. M.
1992
Applied Groundwater Modeling
.
Academic Press, San Diego, 381 pp.
Atkinson
K. E.
2008
An Introduction to Numerical Analysis
.
John Wiley & Sons
,
New York
.
Craig
J.
2009
Analytic elements for flow in harmonically heterogeneous aquifers
.
Water Resour. Res.
45
(
6
), W06422.
Fagherazzi
S.
Furbish
D. J.
Rasetarinera
P.
Hussaini
M. Y.
2004
Application of the discontinuous spectral Galerkin method to groundwater flow
.
Adv. Water Resour.
27
(
2
),
129
140
.
Giudici
M.
Vassena
C.
2008
Spectral analysis of the balance equation of ground water hydrology
.
Transport Porous Media
72
(
2
),
171
178
.
Golberg
M. A.
Chen
C.-S.
1997
Discrete Projection Methods for Integral Equations
.
Computational Mechanics Publications
,
Southampton
.
Haitjema
H.
Kuzin
S.
Kelson
V.
Abrams
D.
2010
Modeling flow into horizontal wells in a Dupuit-Forchheimer model
.
Ground Water
48
(
6
),
878
883
.
Homayoon
L.
Abedini
M. J.
Hashemi
S. M. R.
2013
RBF-DQ solution for shallow water equations
.
J. Waterway Port Coast. Ocean Eng.
139
(
1
),
45
60
.
Kaya
B.
Arisoy
Y.
2011
Differential quadrature solution for one-dimensional aquifer flow
.
Math. Comput. Applic.
16
(
2
),
524
.
Malekzadeh
P.
Rahideh
H.
Setoodeh
A.
2007
Optimization of non-symmetric convective–radiative annular fins by differential quadrature method
.
Energy Conversion Manage.
48
(
5
),
1671
1677
.
Moore
R.
Kelson
V.
Wittman
J.
Rash
V.
2012
A modeling framework for the design of collector wells
.
Ground Water
50
(
3
),
355
366
.
Moridis
G.
Kansa
E.
1992
The method of Laplace Transform MultiQuadrics (LTMQ) for the solution of the groundwater flow equation
.
NASA STI/Recon Technical Report No. 94, 14053
.
Theis
C. V.
1935
The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using ground water storage. US Department of the Interior, Geological Survey, Water Resources Division, Ground Water Branch
.
Wang
H. F.
Anderson
M. P.
1995
Introduction to Groundwater Modeling: Finite Difference and Finite Element Methods
.
Academic Press
,
San Diego, CA
.