Performance assessment of low-order versus high-order numerical schemes in the numerical simulation of aquifer flow

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.

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.

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.

Material properties and other groundwater flow parameters used by Samani et al. (2004) are: T= 8 × 10⁻⁵ m²/s, S= 8 × 10⁻³, h₀ = 100 m, Q= 6.28 × 10⁻⁴ m³/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.

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.

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

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.

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.

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.